IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 30288


Ignore:
Timestamp:
Jan 17, 2011, 5:07:46 PM (15 years ago)
Author:
eugene
Message:

test chisq improvements from increasing the spatial order and DUAL vs SINGLE1,2; reject stamps on the basis of the chisq vs flux2 model fit (allowing for systematic errors in the fit); track the best fit with the new pmSubtractionQuality structure; split out pmSubtractionConvolveStamps from pmSubtractionCalculateEquation; replace pmSubtractionCalculateDeviations with pmSubtractionCalculateChisqAndMoments

Location:
branches/eam_branches/ipp-20101205/psModules/src/imcombine
Files:
4 edited

Legend:

Unmodified
Added
Removed

  • branches/eam_branches/ipp-20101205/psModules/src/imcombine/pmSubtraction.c

    r29777 r30288  
    773773
    774774    if (convolutions) {
    775         // Already done
    776775        return convolutions;
    777776    }
    … …  
    787786}
    788787
     788
     789bool pmSubtractionConvolveStampThread(psThreadJob *job)
     790{
     791    PS_ASSERT_THREAD_JOB_NON_NULL(job, false);
     792
     793    pmSubtractionStamp *stamp = job->args->data[0]; // List of stamps
     794    pmSubtractionKernels *kernels = job->args->data[1]; // Kernels
     795    int footprint = PS_SCALAR_VALUE(job->args->data[2], S32); // Stamp index
     796
     797    return pmSubtractionConvolveStamp(stamp, kernels, footprint);
     798}
    789799
    790800bool pmSubtractionConvolveStamp (pmSubtractionStamp *stamp, pmSubtractionKernels *kernels, int footprint)
    … …  
    818828    }
    819829
     830#ifdef TESTING
     831    for (int j = 0; j < kernels->num; j++) {
     832        if (stamp->convolutions1) {
     833            psString convName = NULL;
     834            psStringAppend(&convName, "conv1_%03d_%03d.fits", index, j);
     835            psFits *fits = psFitsOpen(convName, "w");
     836            psFree(convName);
     837            psKernel *conv = stamp->convolutions1->data[j];
     838            psFitsWriteImage(fits, NULL, conv->image, 0, NULL);
     839            psFitsClose(fits);
     840        }
     841
     842        if (stamp->convolutions2) {
     843            psString convName = NULL;
     844            psStringAppend(&convName, "conv2_%03d_%03d.fits", index, j);
     845            psFits *fits = psFitsOpen(convName, "w");
     846            psFree(convName);
     847            psKernel *conv = stamp->convolutions2->data[j];
     848            psFitsWriteImage(fits, NULL, conv->image, 0, NULL);
     849            psFitsClose(fits);
     850        }
     851    }
     852#endif
     853
    820854    return true;
    821855}
    822856
     857bool pmSubtractionConvolveStamps(pmSubtractionStampList *stamps, pmSubtractionKernels *kernels)
     858{
     859    PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(stamps, false);
     860    PM_ASSERT_SUBTRACTION_KERNELS_NON_NULL(kernels, false);
     861
     862    psTimerStart("pmSubtractionConvolveStamps");
     863
     864    int footprint = stamps->footprint;  // Half-size of stamps
     865
     866    // We iterate over each stamp and generate the convolution if needed.  We do NOT need the
     867    // convolution if (a) it has already been calculated or (b) the stamp is not available for
     868    // use (available = USED or CALCULATE)
     869   
     870    for (int i = 0; i < stamps->num; i++) {
     871        pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest
     872
     873        bool keep = false;
     874        keep |= (stamp->status == PM_SUBTRACTION_STAMP_USED);
     875        keep |= (stamp->status == PM_SUBTRACTION_STAMP_CALCULATE);
     876        if (!keep) continue;
     877
     878        bool haveConvolutions = false;
     879        if (kernels->mode == PM_SUBTRACTION_MODE_1) {
     880            haveConvolutions = (stamp->convolutions1 != NULL);
     881        }
     882        if (kernels->mode == PM_SUBTRACTION_MODE_2) {
     883            haveConvolutions = (stamp->convolutions2 != NULL);
     884        }
     885        if (kernels->mode == PM_SUBTRACTION_MODE_DUAL) {
     886            haveConvolutions = (stamp->convolutions1 != NULL) && (stamp->convolutions2 != NULL);
     887        }
     888        if (haveConvolutions) {
     889            continue;
     890        }
     891
     892        if (pmSubtractionThreaded()) {
     893            psThreadJob *job = psThreadJobAlloc("PSMODULES_SUBTRACTION_CONVOLVE_STAMP");
     894            psArrayAdd(job->args, 1, stamp);
     895            psArrayAdd(job->args, 1, kernels);
     896            PS_ARRAY_ADD_SCALAR(job->args, footprint, PS_TYPE_S32);
     897            if (!psThreadJobAddPending(job)) {
     898                return false;
     899            }
     900        } else {
     901            pmSubtractionConvolveStamp(stamp, kernels, footprint);
     902        }
     903    }
     904    if (!psThreadPoolWait(true)) {
     905        psError(psErrorCodeLast(), false, "Error waiting for threads.");
     906        return false;
     907    }
     908    psLogMsg("psModules.imcombine", PS_LOG_INFO, "Convolve stamps: %f sec", psTimerClear("pmSubtractionConvolveStamps"));
     909    return true;
     910}
    823911
    824912int pmSubtractionRejectStamps(pmSubtractionKernels *kernels, pmSubtractionStampList *stamps,
    825                               const psVector *deviations, psImage *subMask, float sigmaRej)
     913                              pmSubtractionQuality *match, psImage *subMask, float sigmaRej)
    826914{
    827915    PM_ASSERT_SUBTRACTION_KERNELS_NON_NULL(kernels, false);
    828916    PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(stamps, -1);
    829     PS_ASSERT_VECTOR_NON_NULL(deviations, -1);
    830     PS_ASSERT_VECTOR_TYPE(deviations, PS_TYPE_F32, -1);
    831917    PS_ASSERT_IMAGE_NON_EMPTY(subMask, -1);
    832918    PS_ASSERT_IMAGE_TYPE(subMask, PS_TYPE_IMAGE_MASK, -1);
    833919
    834     // I used to measure the rms deviation about zero, and use that as the sigma against which to clip, but
    835     // the distribution is actually something like a chi^2 or Student's t, both of which become Gaussian-like
    836     // with large N.  Therefore, let's just treat this as a Gaussian distribution.
     920    // Comment from PAP (r18287): I used to measure the rms deviation about zero, and use that as the
     921    // sigma against which to clip, but the distribution is actually something like a chi^2 or
     922    // Student's t, both of which become Gaussian-like with large N.  Therefore, let's just
     923    // treat this as a Gaussian distribution.
     924
     925    // Comment from EAM (r29777): The residual distribution is only chisq-like if the model is
     926    // a good fit to the data.  In the (likely) case that there is a systematic difference
     927    // between the model and the data, the squared-residual distribution grows quadratically
     928    // with increasing flux: the systematic residual flux is a constant factor times the source
     929    // flux; the squared-residual is then of the form (k0 + k1*flux)^2, where k0 comes from the
     930    // Gaussian distributed residual and k1*flux is the systematic residual error.
     931
     932    // By rejecting sources with the largest squared-residuals, the rejection biases against
     933    // the brighter sources; in severe cases, this pushes the measurement to the weakest
     934    // sources with the most noise.  To account for this, let's fit a 2nd order polynomial to
     935    // the distribution of flux vs squared-residual, subtract that fit, and reject sources
     936    // which are significantly deviant from that distribution.
    837937
    838938    kernels->mean = NAN;
    … …  
    840940    kernels->numStamps = -1;
    841941
    842     int numStamps = 0;                  // Number of used stamps
    843     psVector *mask = psVectorAlloc(stamps->num, PS_TYPE_VECTOR_MASK); // Mask, for statistics
    844     psVectorInit(mask, 0);
    845     for (int i = 0; i < stamps->num; i++) {
    846         pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest
    847         if (stamp->status != PM_SUBTRACTION_STAMP_USED) {
    848             mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = 0xff;
    849             continue;
    850         }
    851         numStamps++;
    852     }
    853     psTrace("psModules.imcombine", 1, "Number of good stamps: %d\n", numStamps);
    854 
    855     if (numStamps == 0) {
    856         psError(PM_ERR_STAMPS, true, "No good stamps found.");
    857         psFree(mask);
    858         return -1;
    859     }
    860 
    861     psStats *stats = psStatsAlloc(PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV |
    862                                   PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_QUARTILE); // Statistics for deviatns
    863     if (!psVectorStats(stats, deviations, NULL, mask, 0xff)) {
    864         psError(PM_ERR_DATA, false, "Unable to measure statistics for deviations.");
     942    psTrace("psModules.imcombine", 1, "Number of good stamps: %d\n", match->nGood);
     943
     944    // the chisq & flux vectors are calculated by pmSubtractionCalculateChisqAndMoments
     945
     946    // use 3hi/3lo sigma clipping on the chisq fit
     947    psStats *stats = psStatsAlloc(PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV);
     948    stats->clipSigma = 5.0;
     949    stats->clipIter = 2;
     950    psPolynomial1D *model = psPolynomial1DAlloc (PS_POLYNOMIAL_ORD, 2);
     951
     952    bool result = psVectorClipFitPolynomial1D(model, stats, match->stampMask, 0xff, match->chisq, NULL, match->fluxes);
     953    if (!result) {
     954        psError(PM_ERR_DATA, false, "Unable to measure statistics for deviations.");
     955        psFree(model);
    865956        psFree(stats);
    866         psFree(mask);
    867         return -1;
    868     }
    869     psFree(mask);
    870 
    871     // XXX raise an error?
     957        return -1;
     958    }
    872959    if (isnan(stats->sampleMean)) {
     960        psError(PM_ERR_DATA, false, "Unable to measure statistics for deviations.");
     961        psFree(model);
    873962        psFree(stats);
    874963        return -1;
    875964    }
    876965
    877     double mean, rms;                 // Mean and RMS of deviations
    878     if (numStamps < MIN_SAMPLE_STATS) {
    879         mean = stats->sampleMean;
    880         rms = stats->sampleStdev;
    881     } else {
    882         mean = stats->sampleMedian;
    883         rms = 0.74 * (stats->sampleUQ - stats->sampleLQ);
    884     }
    885     psFree(stats);
    886 
    887     psTrace("psModules.imcombine", 1, "Mean: %f\n", mean);
    888     psTrace("psModules.imcombine", 1, "RMS deviation: %f\n", rms);
    889 
    890     kernels->mean = mean;
    891     kernels->rms = rms;
    892     kernels->numStamps = numStamps;
    893 
    894     psLogMsg("psModules.imcombine", PS_LOG_INFO, "Mean deviation from %d stamps: %lf +/- %lf",
    895              numStamps, mean, rms);
    896 
    897     if (!isfinite(sigmaRej) || sigmaRej <= 0.0) {
    898         // User just wanted to calculate and record the deviation for posterity
    899         return 0;
    900     }
    901 
    902     float limit = sigmaRej * rms; // Limit on maximum deviation
    903     psTrace("psModules.imcombine", 1, "Deviation limit: %f\n", limit);
    904 
     966    kernels->mean = stats->sampleMean;
     967    kernels->rms = stats->sampleStdev;
     968    kernels->numStamps = stats->clippedNvalues;
     969
     970    psLogMsg ("pmPSFtry", 4, "chisq vs flux fit: %f +/- %f\n", stats->sampleMean, stats->sampleStdev);
     971    psLogMsg("psModules.imcombine", PS_LOG_INFO, "Mean deviation from %d stamps: %lf +/- %lf",  kernels->numStamps, kernels->mean, kernels->rms);
    905972
    906973    psString ds9name = NULL;            // Filename for ds9 region file
    … …  
    914981    int numRejected = 0;                // Number of stamps rejected
    915982    int numGood = 0;                    // Number of good stamps
    916     double newMean = 0.0;               // New mean
    917983    psString log = NULL;                // Log message
    918     psStringAppend(&log, "Rejecting stamps, mean = %f, threshold = %f\n", mean, limit);
     984
     985    // save DS9 region files for the stamps and mark for rejection and replacement
    919986    for (int i = 0; i < stamps->num; i++) {
    920987        pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest
    921         if (stamp->status == PM_SUBTRACTION_STAMP_USED) {
     988        if (stamp->status  != PM_SUBTRACTION_STAMP_USED) { continue; }
     989        if (match->stampMask->data.PS_TYPE_VECTOR_MASK_DATA[i]) {
    922990            // Should we reject stars with low deviation?  Well, if this is really a Gaussian-like
    923991            // distribution and they're low, then we have the right to ask why.  Isn't it suspicious that
    … …  
    926994            // subtract well, in which case very few (if any) stars will be legitimately rejected for being
    927995            // low.
    928             if (fabsf(deviations->data.F32[i] - mean) > limit) {
    929                 // Mask out the stamp in the image so you it's not found again
    930                 psTrace("psModules.imcombine", 3, "Rejecting stamp %d (%d,%d)\n", i,
    931                         (int)(stamp->x - 0.5), (int)(stamp->y - 0.5));
    932                 psStringAppend(&log, "Stamp %d (%d,%d): %f\n", i,
    933                                (int)(stamp->x - 0.5), (int)(stamp->y - 0.5),
    934                                fabsf(deviations->data.F32[i] - mean));
    935                 numRejected++;
    936                 for (int y = stamp->y - footprint; y <= stamp->y + footprint; y++) {
    937                     for (int x = stamp->x - footprint; x <= stamp->x + footprint; x++) {
    938                         subMask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] |= PM_SUBTRACTION_MASK_REJ;
    939                     }
    940                 }
    941                 pmSubtractionStampPrint(ds9, stamp->x, stamp->y, stamps->footprint, "red");
    942 
    943                 // Set stamp for replacement
    944                 stamp->x = 0;
    945                 stamp->y = 0;
    946                 stamp->xNorm = NAN;
    947                 stamp->yNorm = NAN;
    948                 stamp->status = PM_SUBTRACTION_STAMP_REJECTED;
    949                 // Recalculate convolutions
    950                 psFree(stamp->convolutions1);
    951                 psFree(stamp->convolutions2);
    952                 stamp->convolutions1 = stamp->convolutions2 = NULL;
    953                 psFree(stamp->image1);
    954                 psFree(stamp->image2);
    955                 psFree(stamp->weight);
    956                 stamp->image1 = stamp->image2 = stamp->weight = NULL;
    957                 psFree(stamp->matrix);
    958                 stamp->matrix = NULL;
    959                 psFree(stamp->vector);
    960                 stamp->vector = NULL;
    961             } else {
    962                 numGood++;
    963                 newMean += deviations->data.F32[i];
    964                 pmSubtractionStampPrint(ds9, stamp->x, stamp->y, stamps->footprint, "green");
    965             }
    966         }
    967     }
    968     newMean /= numGood;
     996            psTrace("psModules.imcombine", 3, "Rejecting stamp %d (%d,%d)\n", i,
     997                    (int)(stamp->x - 0.5), (int)(stamp->y - 0.5));
     998            psStringAppend(&log, "Stamp %d (%d,%d): %f : %f : %f\n",
     999                           i, (int)(stamp->x - 0.5), (int)(stamp->y - 0.5),
     1000                           match->chisq->data.F32[i], match->fluxes->data.F32[i], match->chisq->data.F32[i] - psPolynomial1DEval(model, match->fluxes->data.F32[i]));
     1001            numRejected++;
     1002            for (int y = stamp->y - footprint; y <= stamp->y + footprint; y++) {
     1003                for (int x = stamp->x - footprint; x <= stamp->x + footprint; x++) {
     1004                    subMask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] |= PM_SUBTRACTION_MASK_REJ;
     1005                }
     1006            }
     1007            pmSubtractionStampPrint(ds9, stamp->x, stamp->y, stamps->footprint, "red");
     1008
     1009            // Set stamp for replacement
     1010            stamp->x = 0;
     1011            stamp->y = 0;
     1012            stamp->xNorm = NAN;
     1013            stamp->yNorm = NAN;
     1014            stamp->status = PM_SUBTRACTION_STAMP_REJECTED;
     1015            // Recalculate convolutions
     1016            psFree(stamp->convolutions1);
     1017            psFree(stamp->convolutions2);
     1018            stamp->convolutions1 = stamp->convolutions2 = NULL;
     1019            psFree(stamp->image1);
     1020            psFree(stamp->image2);
     1021            psFree(stamp->weight);
     1022            stamp->image1 = stamp->image2 = stamp->weight = NULL;
     1023            psFree(stamp->matrix);
     1024            stamp->matrix = NULL;
     1025            psFree(stamp->vector);
     1026            stamp->vector = NULL;
     1027        } else {
     1028            numGood++;
     1029            pmSubtractionStampPrint(ds9, stamp->x, stamp->y, stamps->footprint, "green");
     1030        }
     1031    }
    9691032
    9701033    if (numRejected == 0) {
    … …  
    9781041    }
    9791042
     1043    psFree(model);
     1044    psFree(stats);
     1045
    9801046    if (numRejected > 0) {
    981         psLogMsg("psModules.imcombine", PS_LOG_INFO,
    982                  "%d good stamps; %d rejected.\nMean deviation: %lf --> %lf\n",
    983                  numGood, numRejected, mean, newMean);
     1047        psLogMsg("psModules.imcombine", PS_LOG_INFO, "%d good stamps; %d rejected.\n", numGood, numRejected);
    9841048    } else {
    985         psLogMsg("psModules.imcombine", PS_LOG_INFO,
    986                  "%d good stamps; 0 rejected.\nMean deviation: %lf\n",
    987                  numGood, mean);
     1049        psLogMsg("psModules.imcombine", PS_LOG_INFO, "%d good stamps; 0 rejected.\n", numGood);
    9881050    }
    9891051
    … …  
    14791541  return true;
    14801542}
     1543
     1544static void pmSubtractionQualityFree(pmSubtractionQuality *quality) {
     1545
     1546    psFree (quality->fluxes);
     1547    psFree (quality->chisq);
     1548    psFree (quality->moments);
     1549    psFree (quality->stampMask);
     1550}   
     1551
     1552pmSubtractionQuality *pmSubtractionQualityAlloc() {
     1553
     1554    pmSubtractionQuality *quality = psAlloc(sizeof(pmSubtractionQuality)); // Stamp list to return
     1555    psMemSetDeallocator(quality, (psFreeFunc)pmSubtractionQualityFree);
     1556
     1557    quality->fluxes = NULL;
     1558    quality->chisq = NULL;
     1559    quality->moments = NULL;
     1560    quality->stampMask = NULL;
     1561
     1562    quality->score = NAN;
     1563    quality->mode = PM_SUBTRACTION_MODE_ERR;
     1564    quality->spatialOrder = -1;
     1565    quality->nGood = 0;
     1566   
     1567    return quality;
     1568}

  • branches/eam_branches/ipp-20101205/psModules/src/imcombine/pmSubtraction.h

    r29601 r30288  
    4444} pmSubtractionMasks;
    4545
     46typedef struct {
     47    double score;
     48    pmSubtractionMode mode;
     49    int spatialOrder;
     50    int nGood;
     51    psVector *fluxes;
     52    psVector *chisq;
     53    psVector *moments;
     54    psVector *stampMask;
     55} pmSubtractionQuality;
    4656
    4757/// Number of terms in a polynomial
    … …  
    7080    );
    7181
     82bool pmSubtractionConvolveStamps(pmSubtractionStampList *stamps, pmSubtractionKernels *kernels);
     83
     84bool pmSubtractionConvolveStampThread(psThreadJob *job);
     85
    7286/// Reject stamps
    7387int pmSubtractionRejectStamps(pmSubtractionKernels *kernels, ///< Kernel parameters to update
    7488                              pmSubtractionStampList *stamps, ///< Stamps
    75                               const psVector *deviations, ///< Deviations for each stamp
     89                              pmSubtractionQuality *match, ///< data on the subtraction quality
    7690                              psImage *subMask, ///< Subtraction mask
    7791                              float sigmaRej ///< Number of RMS deviations above zero at which to reject
    … …  
    167181bool pmSubtractionSetFWHMs(float fwhm1, float fwhm2);
    168182
     183pmSubtractionQuality *pmSubtractionQualityAlloc();
     184
    169185/// @}
    170186#endif

  • branches/eam_branches/ipp-20101205/psModules/src/imcombine/pmSubtractionMatch.c

    r30266 r30288  
    2727
    2828static bool useFFT = true;              // Do convolutions using FFT
    29 
    30 # define SUBMODE PM_SUBTRACTION_EQUATION_ALL
    3129
    3230//#define TESTING
    … …  
    473471}
    474472
     473bool pmSubtractionMatchAttempt(pmSubtractionQuality **bestMatch, pmSubtractionKernels *kernels, pmSubtractionStampList *stamps, pmSubtractionMode mode, int spatialOrder, bool final) {
     474
     475    pmSubtractionMode nativeMode = kernels->mode;
     476    pmSubtractionMode nativeOrder = kernels->spatialOrder;
     477
     478    kernels->mode = mode;
     479    kernels->spatialOrder = spatialOrder;
     480
     481    // we always need to recalculate the matrix equation elements...
     482    pmSubtractionStampsResetStatus(stamps);
     483
     484    psTrace("psModules.imcombine", 3, "Convolving stamps as needed...\n");
     485    if (!pmSubtractionConvolveStamps(stamps, kernels)) {
     486        psError(psErrorCodeLast(), false, "Unable to convolve stamps.");
     487        return false;
     488    }
     489
     490    // step 1: generate the elements of the matrix equation Ax = B
     491    psTrace("psModules.imcombine", 3, "Calculating kernel equations...\n");
     492    if (!pmSubtractionCalculateEquation(stamps, kernels)) {
     493        psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation.");
     494        return false;
     495    }
     496                   
     497    // step 2: solve the matrix equation Ax = B
     498    psTrace("psModules.imcombine", 3, "Solving kernel equations...\n");
     499    if (!pmSubtractionSolveEquation(kernels, stamps)) {
     500        psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation.");
     501        return false;
     502    }
     503    memCheck("  solve equation");
     504
     505    // calculate the score for this model fit attempt
     506    // XXX store the chisq, flux and moments for stamp rejection
     507    pmSubtractionCalculateChisqAndMoments(bestMatch, stamps, kernels); // Stamp deviations
     508
     509    // display the input and model stamps
     510    pmSubtractionVisualShowFit(stamps, kernels);
     511    pmSubtractionVisualPlotFit(kernels);
     512    pmSubtractionVisualPlotConvKernels(kernels);
     513
     514    // reset the kernel if desired (on final pass, do not reset)
     515    if (!final) {
     516        kernels->mode = nativeMode;
     517        kernels->spatialOrder = nativeOrder;
     518    }
     519    return true;
     520}
    475521
    476522bool pmSubtractionMatch(pmReadout *conv1, pmReadout *conv2, const pmReadout *ro1, const pmReadout *ro2,
    … …  
    559605    psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS); // Random number generator
    560606
     607    pmSubtractionQuality *bestMatch = NULL;
     608
     609    int N_TEST_MODES;
     610    int N_TEST_ORDER = spatialOrder;
     611
     612    pmSubtractionMode TestModes[3];
     613    switch (subMode) {
     614      case PM_SUBTRACTION_MODE_1:
     615        N_TEST_MODES = 1;
     616        TestModes[0] = PM_SUBTRACTION_MODE_1;
     617        break;
     618      case PM_SUBTRACTION_MODE_2:
     619        N_TEST_MODES = 1;
     620        TestModes[0] = PM_SUBTRACTION_MODE_2;
     621        break;
     622      case PM_SUBTRACTION_MODE_DUAL:
     623        N_TEST_MODES = 3;
     624        TestModes[0] = PM_SUBTRACTION_MODE_1;
     625        TestModes[1] = PM_SUBTRACTION_MODE_2;
     626        TestModes[2] = PM_SUBTRACTION_MODE_DUAL;
     627        break;
     628      default:
     629        psError(psErrorCodeLast(), false, "For now, only modes 1, 2, and DUAL are supported.");
     630        goto MATCH_ERROR;
     631    }
     632   
    561633    memCheck("start");
    562634
    … …  
    689761
    690762            int numRejected = -1;       // Number of rejected stamps in each iteration
    691             for (int k = 0; k < iter && numRejected != 0; k++) {
     763            for (int k = 0; (k < iter) && (numRejected != 0); k++) {
    692764                psLogMsg("psModules.imcombine", PS_LOG_INFO, "Iteration %d.", k);
    693765
    … …  
    711783                }
    712784
    713                 // step 0b : calculate the moments, pass along to the next steps via stamps->normValue
    714                 // XXX this step is now skipped -- delete
    715                 psTrace("psModules.imcombine", 3, "Calculating normalization...\n");
    716                 if (!pmSubtractionCalculateMoments(kernels, stamps)) {
    717                     psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation.");
     785                // on each iteration, we start from scratch
     786                psFree(bestMatch);
     787
     788                // choose the spatial order and subtraction direction (1, 2, dual)
     789                // XXX need to make these respect recipe somewhat
     790                for (int order = 0; order <= N_TEST_ORDER; order++) {
     791                    for (int j = 0; j < N_TEST_MODES; j++) {
     792                        if (!pmSubtractionMatchAttempt(&bestMatch, kernels, stamps, TestModes[j], order, false)) {
     793                            goto MATCH_ERROR;
     794                        }
     795                    }
     796                }
     797               
     798                // reject the deviant stamps based on the stats of the best match
     799                psTrace("psModules.imcombine", 3, "Rejecting stamps...\n");
     800                numRejected = pmSubtractionRejectStamps(kernels, stamps, bestMatch, subMask, rej);
     801                if (numRejected < 0) {
     802                    psError(psErrorCodeLast(), false, "Unable to reject stamps.");
    718803                    goto MATCH_ERROR;
    719804                }
    720 
    721                 // step 1: generate the elements of the matrix equation Ax = B
    722                 psTrace("psModules.imcombine", 3, "Calculating kernel equations...\n");
    723                 if (!pmSubtractionCalculateEquation(stamps, kernels, SUBMODE)) {
    724                     psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation.");
    725                     goto MATCH_ERROR;
    726                 }
    727 
    728                 // step 2: solve the matrix equation Ax = B
    729                 psTrace("psModules.imcombine", 3, "Solving kernel equations...\n");
    730                 if (!pmSubtractionSolveEquation(kernels, stamps, SUBMODE)) {
    731                     psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation.");
    732                     goto MATCH_ERROR;
    733                 }
    734                 memCheck("  solve equation");
    735 
    736                 pmSubtractionVisualPlotConvKernels(kernels);
    737 
    738                 psVector *deviations = pmSubtractionCalculateDeviations(stamps, kernels); // Stamp deviations
    739                 if (!deviations) {
    740                     psError(psErrorCodeLast(), false, "Unable to calculate deviations.");
    741                     goto MATCH_ERROR;
    742                 }
    743                 memCheck("   calculate deviations");
    744 
    745                 pmSubtractionCalculateChisqAndMoments(stamps, kernels); // Stamp deviations
    746 
    747                 psTrace("psModules.imcombine", 3, "Rejecting stamps...\n");
    748                 numRejected = pmSubtractionRejectStamps(kernels, stamps, deviations, subMask, rej);
    749                 if (numRejected < 0) {
    750                     psError(psErrorCodeLast(), false, "Unable to reject stamps.");
    751                     psFree(deviations);
    752                     goto MATCH_ERROR;
    753                 }
    754                 psFree(deviations);
    755 
    756805                memCheck("  reject stamps");
    757806            }
    758807
    759             // if we hit the max number of iterations and we have rejected stamps, re-solve
    760             if (numRejected > 0) {
    761 
    762                 // step 1: generate the elements of the matrix equation Ax = B
    763                 psTrace("psModules.imcombine", 3, "Calculating equation for normalization...\n");
    764                 if (!pmSubtractionCalculateEquation(stamps, kernels, SUBMODE)) {
    765                     psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation.");
    766                     goto MATCH_ERROR;
    767                 }
    768 
    769                 // step 2: solve the matrix equation Ax = B
    770                 psTrace("psModules.imcombine", 3, "Solving equation for kernels...\n");
    771                 if (!pmSubtractionSolveEquation(kernels, stamps, SUBMODE)) {
    772                     psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation.");
    773                     goto MATCH_ERROR;
    774                 }
    775 
    776                 pmSubtractionVisualPlotConvKernels(kernels);
    777 
    778                 psVector *deviations = pmSubtractionCalculateDeviations(stamps, kernels); // Stamp deviations
    779                 if (!deviations) {
    780                     psError(psErrorCodeLast(), false, "Unable to calculate deviations.");
    781                     goto MATCH_ERROR;
    782                 }
    783                 pmSubtractionRejectStamps(kernels, stamps, deviations, subMask, NAN);
    784                 psFree(deviations);
    785             }
     808            // apply the best fit so we are ready to roll
     809            fprintf (stderr, "applying order: %d, mode: %d\n", bestMatch->spatialOrder, bestMatch->mode);
     810            if (!pmSubtractionMatchAttempt(NULL, kernels, stamps, bestMatch->mode, bestMatch->spatialOrder, true)) {
     811                goto MATCH_ERROR;
     812            }
    786813            psFree(stamps);
    787             stamps = NULL;
    788 
     814            psFree(bestMatch);
    789815            memCheck("solution");
    790816
    … …  
    860886    psFree(variance);
    861887    psFree(rng);
     888    psFree(bestMatch);
    862889    return false;
    863890}
    … …  
    10771104    assert(kernels);
    10781105
     1106    psAbort("this function is not working");
     1107# if (0)
     1108    psTrace("psModules.imcombine", 3, "Convolving stamps as needed...\n");
     1109    if (!pmSubtractionConvolveStamps(stamps, kernels)) {
     1110        psError(psErrorCodeLast(), false, "Unable to convolve stamps.");
     1111        return false;
     1112    }
     1113
    10791114    psTrace("psModules.imcombine", 3, "Calculating %s normalization equation...\n", description);
    1080     if (!pmSubtractionCalculateEquation(stamps, kernels, SUBMODE)) {
     1115    if (!pmSubtractionCalculateEquation(stamps, kernels)) {
    10811116        psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation.");
    10821117        return false;
    … …  
    10841119
    10851120    psTrace("psModules.imcombine", 3, "Solving %s normalization equation...\n", description);
    1086     if (!pmSubtractionSolveEquation(kernels, stamps, SUBMODE)) {
     1121    if (!pmSubtractionSolveEquation(kernels, stamps)) {
    10871122        psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation.");
    10881123        return false;
    … …  
    10961131    }
    10971132
     1133    // XXX this needs to be made consistent with the modified 'reject stamps' function
    10981134    psTrace("psModules.imcombine", 3, "Rejecting %s stamps...\n", description);
    10991135    long numRejected = pmSubtractionRejectStamps(kernels, stamps, deviations, subMask, rej);
    … …  
    11081144        // Allow re-fit with reduced stamps set
    11091145        psTrace("psModules.imcombine", 3, "Calculating %s normalization equation...\n", description);
    1110         if (!pmSubtractionCalculateEquation(stamps, kernels, PM_SUBTRACTION_EQUATION_ALL)) {
     1146        if (!pmSubtractionCalculateEquation(stamps, kernels)) {
    11111147            psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation.");
    11121148            return false;
    … …  
    11141150
    11151151        psTrace("psModules.imcombine", 3, "Resolving %s equation...\n", description);
    1116         if (!pmSubtractionSolveEquation(kernels, stamps, PM_SUBTRACTION_EQUATION_ALL)) {
     1152        if (!pmSubtractionSolveEquation(kernels, stamps)) {
    11171153            psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation.");
    11181154            return false;
    … …  
    11341170        psFree(deviations);
    11351171    }
    1136 
     1172# endif
    11371173    return true;
    11381174}

  • branches/eam_branches/ipp-20101205/psModules/src/imcombine/pmSubtractionMatch.h

    r29543 r30288  
    88#include <pmSubtractionKernels.h>
    99#include <pmSubtractionStamps.h>
     10#include <pmSubtraction.h>
    1011
    1112/// Match two images
    … …  
    109110    );
    110111
     112bool pmSubtractionMatchAttempt(
     113    pmSubtractionQuality **bestMatch,
     114    pmSubtractionKernels *kernels,
     115    pmSubtractionStampList *stamps,
     116    pmSubtractionMode mode,
     117    int spatialOrder,
     118    bool final
     119    );
     120
    111121#endif
Note: See TracChangeset for help on using the changeset viewer.