Index: /branches/eam_branches/ipp-20101205/psModules/src/imcombine/pmSubtraction.c =================================================================== --- /branches/eam_branches/ipp-20101205/psModules/src/imcombine/pmSubtraction.c (revision 30287) +++ /branches/eam_branches/ipp-20101205/psModules/src/imcombine/pmSubtraction.c (revision 30288) @@ -773,5 +773,4 @@ if (convolutions) { - // Already done return convolutions; } @@ -787,4 +786,15 @@ } + +bool pmSubtractionConvolveStampThread(psThreadJob *job) +{ + PS_ASSERT_THREAD_JOB_NON_NULL(job, false); + + pmSubtractionStamp *stamp = job->args->data[0]; // List of stamps + pmSubtractionKernels *kernels = job->args->data[1]; // Kernels + int footprint = PS_SCALAR_VALUE(job->args->data[2], S32); // Stamp index + + return pmSubtractionConvolveStamp(stamp, kernels, footprint); +} bool pmSubtractionConvolveStamp (pmSubtractionStamp *stamp, pmSubtractionKernels *kernels, int footprint) @@ -818,21 +828,111 @@ } +#ifdef TESTING + for (int j = 0; j < kernels->num; j++) { + if (stamp->convolutions1) { + psString convName = NULL; + psStringAppend(&convName, "conv1_%03d_%03d.fits", index, j); + psFits *fits = psFitsOpen(convName, "w"); + psFree(convName); + psKernel *conv = stamp->convolutions1->data[j]; + psFitsWriteImage(fits, NULL, conv->image, 0, NULL); + psFitsClose(fits); + } + + if (stamp->convolutions2) { + psString convName = NULL; + psStringAppend(&convName, "conv2_%03d_%03d.fits", index, j); + psFits *fits = psFitsOpen(convName, "w"); + psFree(convName); + psKernel *conv = stamp->convolutions2->data[j]; + psFitsWriteImage(fits, NULL, conv->image, 0, NULL); + psFitsClose(fits); + } + } +#endif + return true; } +bool pmSubtractionConvolveStamps(pmSubtractionStampList *stamps, pmSubtractionKernels *kernels) +{ + PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(stamps, false); + PM_ASSERT_SUBTRACTION_KERNELS_NON_NULL(kernels, false); + + psTimerStart("pmSubtractionConvolveStamps"); + + int footprint = stamps->footprint; // Half-size of stamps + + // We iterate over each stamp and generate the convolution if needed. We do NOT need the + // convolution if (a) it has already been calculated or (b) the stamp is not available for + // use (available = USED or CALCULATE) + + for (int i = 0; i < stamps->num; i++) { + pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest + + bool keep = false; + keep |= (stamp->status == PM_SUBTRACTION_STAMP_USED); + keep |= (stamp->status == PM_SUBTRACTION_STAMP_CALCULATE); + if (!keep) continue; + + bool haveConvolutions = false; + if (kernels->mode == PM_SUBTRACTION_MODE_1) { + haveConvolutions = (stamp->convolutions1 != NULL); + } + if (kernels->mode == PM_SUBTRACTION_MODE_2) { + haveConvolutions = (stamp->convolutions2 != NULL); + } + if (kernels->mode == PM_SUBTRACTION_MODE_DUAL) { + haveConvolutions = (stamp->convolutions1 != NULL) && (stamp->convolutions2 != NULL); + } + if (haveConvolutions) { + continue; + } + + if (pmSubtractionThreaded()) { + psThreadJob *job = psThreadJobAlloc("PSMODULES_SUBTRACTION_CONVOLVE_STAMP"); + psArrayAdd(job->args, 1, stamp); + psArrayAdd(job->args, 1, kernels); + PS_ARRAY_ADD_SCALAR(job->args, footprint, PS_TYPE_S32); + if (!psThreadJobAddPending(job)) { + return false; + } + } else { + pmSubtractionConvolveStamp(stamp, kernels, footprint); + } + } + if (!psThreadPoolWait(true)) { + psError(psErrorCodeLast(), false, "Error waiting for threads."); + return false; + } + psLogMsg("psModules.imcombine", PS_LOG_INFO, "Convolve stamps: %f sec", psTimerClear("pmSubtractionConvolveStamps")); + return true; +} int pmSubtractionRejectStamps(pmSubtractionKernels *kernels, pmSubtractionStampList *stamps, - const psVector *deviations, psImage *subMask, float sigmaRej) + pmSubtractionQuality *match, psImage *subMask, float sigmaRej) { PM_ASSERT_SUBTRACTION_KERNELS_NON_NULL(kernels, false); PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(stamps, -1); - PS_ASSERT_VECTOR_NON_NULL(deviations, -1); - PS_ASSERT_VECTOR_TYPE(deviations, PS_TYPE_F32, -1); PS_ASSERT_IMAGE_NON_EMPTY(subMask, -1); PS_ASSERT_IMAGE_TYPE(subMask, PS_TYPE_IMAGE_MASK, -1); - // I used to measure the rms deviation about zero, and use that as the sigma against which to clip, but - // the distribution is actually something like a chi^2 or Student's t, both of which become Gaussian-like - // with large N. Therefore, let's just treat this as a Gaussian distribution. + // Comment from PAP (r18287): I used to measure the rms deviation about zero, and use that as the + // sigma against which to clip, but the distribution is actually something like a chi^2 or + // Student's t, both of which become Gaussian-like with large N. Therefore, let's just + // treat this as a Gaussian distribution. + + // Comment from EAM (r29777): The residual distribution is only chisq-like if the model is + // a good fit to the data. In the (likely) case that there is a systematic difference + // between the model and the data, the squared-residual distribution grows quadratically + // with increasing flux: the systematic residual flux is a constant factor times the source + // flux; the squared-residual is then of the form (k0 + k1*flux)^2, where k0 comes from the + // Gaussian distributed residual and k1*flux is the systematic residual error. + + // By rejecting sources with the largest squared-residuals, the rejection biases against + // the brighter sources; in severe cases, this pushes the measurement to the weakest + // sources with the most noise. To account for this, let's fit a 2nd order polynomial to + // the distribution of flux vs squared-residual, subtract that fit, and reject sources + // which are significantly deviant from that distribution. kernels->mean = NAN; @@ -840,67 +940,34 @@ kernels->numStamps = -1; - int numStamps = 0; // Number of used stamps - psVector *mask = psVectorAlloc(stamps->num, PS_TYPE_VECTOR_MASK); // Mask, for statistics - psVectorInit(mask, 0); - for (int i = 0; i < stamps->num; i++) { - pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest - if (stamp->status != PM_SUBTRACTION_STAMP_USED) { - mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = 0xff; - continue; - } - numStamps++; - } - psTrace("psModules.imcombine", 1, "Number of good stamps: %d\n", numStamps); - - if (numStamps == 0) { - psError(PM_ERR_STAMPS, true, "No good stamps found."); - psFree(mask); - return -1; - } - - psStats *stats = psStatsAlloc(PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV | - PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_QUARTILE); // Statistics for deviatns - if (!psVectorStats(stats, deviations, NULL, mask, 0xff)) { - psError(PM_ERR_DATA, false, "Unable to measure statistics for deviations."); + psTrace("psModules.imcombine", 1, "Number of good stamps: %d\n", match->nGood); + + // the chisq & flux vectors are calculated by pmSubtractionCalculateChisqAndMoments + + // use 3hi/3lo sigma clipping on the chisq fit + psStats *stats = psStatsAlloc(PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV); + stats->clipSigma = 5.0; + stats->clipIter = 2; + psPolynomial1D *model = psPolynomial1DAlloc (PS_POLYNOMIAL_ORD, 2); + + bool result = psVectorClipFitPolynomial1D(model, stats, match->stampMask, 0xff, match->chisq, NULL, match->fluxes); + if (!result) { + psError(PM_ERR_DATA, false, "Unable to measure statistics for deviations."); + psFree(model); psFree(stats); - psFree(mask); - return -1; - } - psFree(mask); - - // XXX raise an error? + return -1; + } if (isnan(stats->sampleMean)) { + psError(PM_ERR_DATA, false, "Unable to measure statistics for deviations."); + psFree(model); psFree(stats); return -1; } - double mean, rms; // Mean and RMS of deviations - if (numStamps < MIN_SAMPLE_STATS) { - mean = stats->sampleMean; - rms = stats->sampleStdev; - } else { - mean = stats->sampleMedian; - rms = 0.74 * (stats->sampleUQ - stats->sampleLQ); - } - psFree(stats); - - psTrace("psModules.imcombine", 1, "Mean: %f\n", mean); - psTrace("psModules.imcombine", 1, "RMS deviation: %f\n", rms); - - kernels->mean = mean; - kernels->rms = rms; - kernels->numStamps = numStamps; - - psLogMsg("psModules.imcombine", PS_LOG_INFO, "Mean deviation from %d stamps: %lf +/- %lf", - numStamps, mean, rms); - - if (!isfinite(sigmaRej) || sigmaRej <= 0.0) { - // User just wanted to calculate and record the deviation for posterity - return 0; - } - - float limit = sigmaRej * rms; // Limit on maximum deviation - psTrace("psModules.imcombine", 1, "Deviation limit: %f\n", limit); - + kernels->mean = stats->sampleMean; + kernels->rms = stats->sampleStdev; + kernels->numStamps = stats->clippedNvalues; + + psLogMsg ("pmPSFtry", 4, "chisq vs flux fit: %f +/- %f\n", stats->sampleMean, stats->sampleStdev); + psLogMsg("psModules.imcombine", PS_LOG_INFO, "Mean deviation from %d stamps: %lf +/- %lf", kernels->numStamps, kernels->mean, kernels->rms); psString ds9name = NULL; // Filename for ds9 region file @@ -914,10 +981,11 @@ int numRejected = 0; // Number of stamps rejected int numGood = 0; // Number of good stamps - double newMean = 0.0; // New mean psString log = NULL; // Log message - psStringAppend(&log, "Rejecting stamps, mean = %f, threshold = %f\n", mean, limit); + + // save DS9 region files for the stamps and mark for rejection and replacement for (int i = 0; i < stamps->num; i++) { pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest - if (stamp->status == PM_SUBTRACTION_STAMP_USED) { + if (stamp->status != PM_SUBTRACTION_STAMP_USED) { continue; } + if (match->stampMask->data.PS_TYPE_VECTOR_MASK_DATA[i]) { // Should we reject stars with low deviation? Well, if this is really a Gaussian-like // distribution and they're low, then we have the right to ask why. Isn't it suspicious that @@ -926,45 +994,40 @@ // subtract well, in which case very few (if any) stars will be legitimately rejected for being // low. - if (fabsf(deviations->data.F32[i] - mean) > limit) { - // Mask out the stamp in the image so you it's not found again - psTrace("psModules.imcombine", 3, "Rejecting stamp %d (%d,%d)\n", i, - (int)(stamp->x - 0.5), (int)(stamp->y - 0.5)); - psStringAppend(&log, "Stamp %d (%d,%d): %f\n", i, - (int)(stamp->x - 0.5), (int)(stamp->y - 0.5), - fabsf(deviations->data.F32[i] - mean)); - numRejected++; - for (int y = stamp->y - footprint; y <= stamp->y + footprint; y++) { - for (int x = stamp->x - footprint; x <= stamp->x + footprint; x++) { - subMask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] |= PM_SUBTRACTION_MASK_REJ; - } - } - pmSubtractionStampPrint(ds9, stamp->x, stamp->y, stamps->footprint, "red"); - - // Set stamp for replacement - stamp->x = 0; - stamp->y = 0; - stamp->xNorm = NAN; - stamp->yNorm = NAN; - stamp->status = PM_SUBTRACTION_STAMP_REJECTED; - // Recalculate convolutions - psFree(stamp->convolutions1); - psFree(stamp->convolutions2); - stamp->convolutions1 = stamp->convolutions2 = NULL; - psFree(stamp->image1); - psFree(stamp->image2); - psFree(stamp->weight); - stamp->image1 = stamp->image2 = stamp->weight = NULL; - psFree(stamp->matrix); - stamp->matrix = NULL; - psFree(stamp->vector); - stamp->vector = NULL; - } else { - numGood++; - newMean += deviations->data.F32[i]; - pmSubtractionStampPrint(ds9, stamp->x, stamp->y, stamps->footprint, "green"); - } - } - } - newMean /= numGood; + psTrace("psModules.imcombine", 3, "Rejecting stamp %d (%d,%d)\n", i, + (int)(stamp->x - 0.5), (int)(stamp->y - 0.5)); + psStringAppend(&log, "Stamp %d (%d,%d): %f : %f : %f\n", + i, (int)(stamp->x - 0.5), (int)(stamp->y - 0.5), + match->chisq->data.F32[i], match->fluxes->data.F32[i], match->chisq->data.F32[i] - psPolynomial1DEval(model, match->fluxes->data.F32[i])); + numRejected++; + for (int y = stamp->y - footprint; y <= stamp->y + footprint; y++) { + for (int x = stamp->x - footprint; x <= stamp->x + footprint; x++) { + subMask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] |= PM_SUBTRACTION_MASK_REJ; + } + } + pmSubtractionStampPrint(ds9, stamp->x, stamp->y, stamps->footprint, "red"); + + // Set stamp for replacement + stamp->x = 0; + stamp->y = 0; + stamp->xNorm = NAN; + stamp->yNorm = NAN; + stamp->status = PM_SUBTRACTION_STAMP_REJECTED; + // Recalculate convolutions + psFree(stamp->convolutions1); + psFree(stamp->convolutions2); + stamp->convolutions1 = stamp->convolutions2 = NULL; + psFree(stamp->image1); + psFree(stamp->image2); + psFree(stamp->weight); + stamp->image1 = stamp->image2 = stamp->weight = NULL; + psFree(stamp->matrix); + stamp->matrix = NULL; + psFree(stamp->vector); + stamp->vector = NULL; + } else { + numGood++; + pmSubtractionStampPrint(ds9, stamp->x, stamp->y, stamps->footprint, "green"); + } + } if (numRejected == 0) { @@ -978,12 +1041,11 @@ } + psFree(model); + psFree(stats); + if (numRejected > 0) { - psLogMsg("psModules.imcombine", PS_LOG_INFO, - "%d good stamps; %d rejected.\nMean deviation: %lf --> %lf\n", - numGood, numRejected, mean, newMean); + psLogMsg("psModules.imcombine", PS_LOG_INFO, "%d good stamps; %d rejected.\n", numGood, numRejected); } else { - psLogMsg("psModules.imcombine", PS_LOG_INFO, - "%d good stamps; 0 rejected.\nMean deviation: %lf\n", - numGood, mean); + psLogMsg("psModules.imcombine", PS_LOG_INFO, "%d good stamps; 0 rejected.\n", numGood); } @@ -1479,2 +1541,28 @@ return true; } + +static void pmSubtractionQualityFree(pmSubtractionQuality *quality) { + + psFree (quality->fluxes); + psFree (quality->chisq); + psFree (quality->moments); + psFree (quality->stampMask); +} + +pmSubtractionQuality *pmSubtractionQualityAlloc() { + + pmSubtractionQuality *quality = psAlloc(sizeof(pmSubtractionQuality)); // Stamp list to return + psMemSetDeallocator(quality, (psFreeFunc)pmSubtractionQualityFree); + + quality->fluxes = NULL; + quality->chisq = NULL; + quality->moments = NULL; + quality->stampMask = NULL; + + quality->score = NAN; + quality->mode = PM_SUBTRACTION_MODE_ERR; + quality->spatialOrder = -1; + quality->nGood = 0; + + return quality; +} Index: /branches/eam_branches/ipp-20101205/psModules/src/imcombine/pmSubtraction.h =================================================================== --- /branches/eam_branches/ipp-20101205/psModules/src/imcombine/pmSubtraction.h (revision 30287) +++ /branches/eam_branches/ipp-20101205/psModules/src/imcombine/pmSubtraction.h (revision 30288) @@ -44,4 +44,14 @@ } pmSubtractionMasks; +typedef struct { + double score; + pmSubtractionMode mode; + int spatialOrder; + int nGood; + psVector *fluxes; + psVector *chisq; + psVector *moments; + psVector *stampMask; +} pmSubtractionQuality; /// Number of terms in a polynomial @@ -70,8 +80,12 @@ ); +bool pmSubtractionConvolveStamps(pmSubtractionStampList *stamps, pmSubtractionKernels *kernels); + +bool pmSubtractionConvolveStampThread(psThreadJob *job); + /// Reject stamps int pmSubtractionRejectStamps(pmSubtractionKernels *kernels, ///< Kernel parameters to update pmSubtractionStampList *stamps, ///< Stamps - const psVector *deviations, ///< Deviations for each stamp + pmSubtractionQuality *match, ///< data on the subtraction quality psImage *subMask, ///< Subtraction mask float sigmaRej ///< Number of RMS deviations above zero at which to reject @@ -167,4 +181,6 @@ bool pmSubtractionSetFWHMs(float fwhm1, float fwhm2); +pmSubtractionQuality *pmSubtractionQualityAlloc(); + /// @} #endif Index: /branches/eam_branches/ipp-20101205/psModules/src/imcombine/pmSubtractionMatch.c =================================================================== --- /branches/eam_branches/ipp-20101205/psModules/src/imcombine/pmSubtractionMatch.c (revision 30287) +++ /branches/eam_branches/ipp-20101205/psModules/src/imcombine/pmSubtractionMatch.c (revision 30288) @@ -27,6 +27,4 @@ static bool useFFT = true; // Do convolutions using FFT - -# define SUBMODE PM_SUBTRACTION_EQUATION_ALL //#define TESTING @@ -473,4 +471,52 @@ } +bool pmSubtractionMatchAttempt(pmSubtractionQuality **bestMatch, pmSubtractionKernels *kernels, pmSubtractionStampList *stamps, pmSubtractionMode mode, int spatialOrder, bool final) { + + pmSubtractionMode nativeMode = kernels->mode; + pmSubtractionMode nativeOrder = kernels->spatialOrder; + + kernels->mode = mode; + kernels->spatialOrder = spatialOrder; + + // we always need to recalculate the matrix equation elements... + pmSubtractionStampsResetStatus(stamps); + + psTrace("psModules.imcombine", 3, "Convolving stamps as needed...\n"); + if (!pmSubtractionConvolveStamps(stamps, kernels)) { + psError(psErrorCodeLast(), false, "Unable to convolve stamps."); + return false; + } + + // step 1: generate the elements of the matrix equation Ax = B + psTrace("psModules.imcombine", 3, "Calculating kernel equations...\n"); + if (!pmSubtractionCalculateEquation(stamps, kernels)) { + psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation."); + return false; + } + + // step 2: solve the matrix equation Ax = B + psTrace("psModules.imcombine", 3, "Solving kernel equations...\n"); + if (!pmSubtractionSolveEquation(kernels, stamps)) { + psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation."); + return false; + } + memCheck(" solve equation"); + + // calculate the score for this model fit attempt + // XXX store the chisq, flux and moments for stamp rejection + pmSubtractionCalculateChisqAndMoments(bestMatch, stamps, kernels); // Stamp deviations + + // display the input and model stamps + pmSubtractionVisualShowFit(stamps, kernels); + pmSubtractionVisualPlotFit(kernels); + pmSubtractionVisualPlotConvKernels(kernels); + + // reset the kernel if desired (on final pass, do not reset) + if (!final) { + kernels->mode = nativeMode; + kernels->spatialOrder = nativeOrder; + } + return true; +} bool pmSubtractionMatch(pmReadout *conv1, pmReadout *conv2, const pmReadout *ro1, const pmReadout *ro2, @@ -559,4 +605,30 @@ psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS); // Random number generator + pmSubtractionQuality *bestMatch = NULL; + + int N_TEST_MODES; + int N_TEST_ORDER = spatialOrder; + + pmSubtractionMode TestModes[3]; + switch (subMode) { + case PM_SUBTRACTION_MODE_1: + N_TEST_MODES = 1; + TestModes[0] = PM_SUBTRACTION_MODE_1; + break; + case PM_SUBTRACTION_MODE_2: + N_TEST_MODES = 1; + TestModes[0] = PM_SUBTRACTION_MODE_2; + break; + case PM_SUBTRACTION_MODE_DUAL: + N_TEST_MODES = 3; + TestModes[0] = PM_SUBTRACTION_MODE_1; + TestModes[1] = PM_SUBTRACTION_MODE_2; + TestModes[2] = PM_SUBTRACTION_MODE_DUAL; + break; + default: + psError(psErrorCodeLast(), false, "For now, only modes 1, 2, and DUAL are supported."); + goto MATCH_ERROR; + } + memCheck("start"); @@ -689,5 +761,5 @@ int numRejected = -1; // Number of rejected stamps in each iteration - for (int k = 0; k < iter && numRejected != 0; k++) { + for (int k = 0; (k < iter) && (numRejected != 0); k++) { psLogMsg("psModules.imcombine", PS_LOG_INFO, "Iteration %d.", k); @@ -711,80 +783,34 @@ } - // step 0b : calculate the moments, pass along to the next steps via stamps->normValue - // XXX this step is now skipped -- delete - psTrace("psModules.imcombine", 3, "Calculating normalization...\n"); - if (!pmSubtractionCalculateMoments(kernels, stamps)) { - psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation."); + // on each iteration, we start from scratch + psFree(bestMatch); + + // choose the spatial order and subtraction direction (1, 2, dual) + // XXX need to make these respect recipe somewhat + for (int order = 0; order <= N_TEST_ORDER; order++) { + for (int j = 0; j < N_TEST_MODES; j++) { + if (!pmSubtractionMatchAttempt(&bestMatch, kernels, stamps, TestModes[j], order, false)) { + goto MATCH_ERROR; + } + } + } + + // reject the deviant stamps based on the stats of the best match + psTrace("psModules.imcombine", 3, "Rejecting stamps...\n"); + numRejected = pmSubtractionRejectStamps(kernels, stamps, bestMatch, subMask, rej); + if (numRejected < 0) { + psError(psErrorCodeLast(), false, "Unable to reject stamps."); goto MATCH_ERROR; } - - // step 1: generate the elements of the matrix equation Ax = B - psTrace("psModules.imcombine", 3, "Calculating kernel equations...\n"); - if (!pmSubtractionCalculateEquation(stamps, kernels, SUBMODE)) { - psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation."); - goto MATCH_ERROR; - } - - // step 2: solve the matrix equation Ax = B - psTrace("psModules.imcombine", 3, "Solving kernel equations...\n"); - if (!pmSubtractionSolveEquation(kernels, stamps, SUBMODE)) { - psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation."); - goto MATCH_ERROR; - } - memCheck(" solve equation"); - - pmSubtractionVisualPlotConvKernels(kernels); - - psVector *deviations = pmSubtractionCalculateDeviations(stamps, kernels); // Stamp deviations - if (!deviations) { - psError(psErrorCodeLast(), false, "Unable to calculate deviations."); - goto MATCH_ERROR; - } - memCheck(" calculate deviations"); - - pmSubtractionCalculateChisqAndMoments(stamps, kernels); // Stamp deviations - - psTrace("psModules.imcombine", 3, "Rejecting stamps...\n"); - numRejected = pmSubtractionRejectStamps(kernels, stamps, deviations, subMask, rej); - if (numRejected < 0) { - psError(psErrorCodeLast(), false, "Unable to reject stamps."); - psFree(deviations); - goto MATCH_ERROR; - } - psFree(deviations); - memCheck(" reject stamps"); } - // if we hit the max number of iterations and we have rejected stamps, re-solve - if (numRejected > 0) { - - // step 1: generate the elements of the matrix equation Ax = B - psTrace("psModules.imcombine", 3, "Calculating equation for normalization...\n"); - if (!pmSubtractionCalculateEquation(stamps, kernels, SUBMODE)) { - psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation."); - goto MATCH_ERROR; - } - - // step 2: solve the matrix equation Ax = B - psTrace("psModules.imcombine", 3, "Solving equation for kernels...\n"); - if (!pmSubtractionSolveEquation(kernels, stamps, SUBMODE)) { - psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation."); - goto MATCH_ERROR; - } - - pmSubtractionVisualPlotConvKernels(kernels); - - psVector *deviations = pmSubtractionCalculateDeviations(stamps, kernels); // Stamp deviations - if (!deviations) { - psError(psErrorCodeLast(), false, "Unable to calculate deviations."); - goto MATCH_ERROR; - } - pmSubtractionRejectStamps(kernels, stamps, deviations, subMask, NAN); - psFree(deviations); - } + // apply the best fit so we are ready to roll + fprintf (stderr, "applying order: %d, mode: %d\n", bestMatch->spatialOrder, bestMatch->mode); + if (!pmSubtractionMatchAttempt(NULL, kernels, stamps, bestMatch->mode, bestMatch->spatialOrder, true)) { + goto MATCH_ERROR; + } psFree(stamps); - stamps = NULL; - + psFree(bestMatch); memCheck("solution"); @@ -860,4 +886,5 @@ psFree(variance); psFree(rng); + psFree(bestMatch); return false; } @@ -1077,6 +1104,14 @@ assert(kernels); + psAbort("this function is not working"); +# if (0) + psTrace("psModules.imcombine", 3, "Convolving stamps as needed...\n"); + if (!pmSubtractionConvolveStamps(stamps, kernels)) { + psError(psErrorCodeLast(), false, "Unable to convolve stamps."); + return false; + } + psTrace("psModules.imcombine", 3, "Calculating %s normalization equation...\n", description); - if (!pmSubtractionCalculateEquation(stamps, kernels, SUBMODE)) { + if (!pmSubtractionCalculateEquation(stamps, kernels)) { psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation."); return false; @@ -1084,5 +1119,5 @@ psTrace("psModules.imcombine", 3, "Solving %s normalization equation...\n", description); - if (!pmSubtractionSolveEquation(kernels, stamps, SUBMODE)) { + if (!pmSubtractionSolveEquation(kernels, stamps)) { psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation."); return false; @@ -1096,4 +1131,5 @@ } + // XXX this needs to be made consistent with the modified 'reject stamps' function psTrace("psModules.imcombine", 3, "Rejecting %s stamps...\n", description); long numRejected = pmSubtractionRejectStamps(kernels, stamps, deviations, subMask, rej); @@ -1108,5 +1144,5 @@ // Allow re-fit with reduced stamps set psTrace("psModules.imcombine", 3, "Calculating %s normalization equation...\n", description); - if (!pmSubtractionCalculateEquation(stamps, kernels, PM_SUBTRACTION_EQUATION_ALL)) { + if (!pmSubtractionCalculateEquation(stamps, kernels)) { psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation."); return false; @@ -1114,5 +1150,5 @@ psTrace("psModules.imcombine", 3, "Resolving %s equation...\n", description); - if (!pmSubtractionSolveEquation(kernels, stamps, PM_SUBTRACTION_EQUATION_ALL)) { + if (!pmSubtractionSolveEquation(kernels, stamps)) { psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation."); return false; @@ -1134,5 +1170,5 @@ psFree(deviations); } - +# endif return true; } Index: /branches/eam_branches/ipp-20101205/psModules/src/imcombine/pmSubtractionMatch.h =================================================================== --- /branches/eam_branches/ipp-20101205/psModules/src/imcombine/pmSubtractionMatch.h (revision 30287) +++ /branches/eam_branches/ipp-20101205/psModules/src/imcombine/pmSubtractionMatch.h (revision 30288) @@ -8,4 +8,5 @@ #include #include +#include /// Match two images @@ -109,3 +110,12 @@ ); +bool pmSubtractionMatchAttempt( + pmSubtractionQuality **bestMatch, + pmSubtractionKernels *kernels, + pmSubtractionStampList *stamps, + pmSubtractionMode mode, + int spatialOrder, + bool final + ); + #endif