Index: trunk/psModules/src/imcombine/Makefile.am =================================================================== --- trunk/psModules/src/imcombine/Makefile.am (revision 30621) +++ trunk/psModules/src/imcombine/Makefile.am (revision 30622) @@ -30,4 +30,5 @@ pmStackReject.h \ pmSubtraction.h \ + pmSubtractionTypes.h \ pmSubtractionAnalysis.h \ pmSubtractionEquation.h \ Index: trunk/psModules/src/imcombine/pmReadoutCombine.c =================================================================== --- trunk/psModules/src/imcombine/pmReadoutCombine.c (revision 30621) +++ trunk/psModules/src/imcombine/pmReadoutCombine.c (revision 30622) @@ -33,5 +33,5 @@ params->blank = 0; params->nKeep = 0; - params->fracHigh = 0.0; + params->fracLow = 0.0; params->fracHigh = 0.0; params->iter = 1; Index: trunk/psModules/src/imcombine/pmStackReject.c =================================================================== --- trunk/psModules/src/imcombine/pmStackReject.c (revision 30621) +++ trunk/psModules/src/imcombine/pmStackReject.c (revision 30622) @@ -7,4 +7,6 @@ #include +#include "pmFPA.h" +#include "pmSubtractionTypes.h" #include "pmSubtraction.h" #include "pmSubtractionThreads.h" Index: trunk/psModules/src/imcombine/pmSubtraction.c =================================================================== --- trunk/psModules/src/imcombine/pmSubtraction.c (revision 30621) +++ trunk/psModules/src/imcombine/pmSubtraction.c (revision 30622) @@ -20,8 +20,11 @@ #include "pmHDU.h" // Required for pmFPA.h #include "pmFPA.h" +#include "pmSubtractionTypes.h" #include "pmSubtractionStamps.h" #include "pmSubtractionEquation.h" #include "pmSubtractionVisual.h" #include "pmSubtractionThreads.h" + +bool psFitsWriteImageSimple (char *filename, psImage *image, psMetadata *header); #include "pmSubtraction.h" @@ -773,5 +776,4 @@ if (convolutions) { - // Already done return convolutions; } @@ -787,4 +789,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 +831,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 +943,35 @@ 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 model: %e + %e flux + %e flux^2\n", model->coeff[0], model->coeff[1], model->coeff[2]); + psLogMsg ("pmPSFtry", 4, "chisq vs flux resid: %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 +985,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 +998,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 +1045,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); } @@ -1373,4 +1439,7 @@ psFree(kernelErr2); + static int nOut1 = 0; + static int nOut2 = 0; + // Calculate covariances // This can be fairly involved, so we only do it for a small number of instances @@ -1386,4 +1455,23 @@ psKernelTruncate(kernel, covarFrac); covars->data[i] = psImageCovarianceCalculate(kernel, ro1->covariance); + if (0) { + char name[128]; + snprintf (name, 128, "covar.sample1.%03d.fits", nOut1); + psKernel *cov = covars->data[i]; + psFitsWriteImageSimple (name, cov->image, NULL); + + snprintf (name, 128, "incovar.sample1.%03d.fits", nOut1); + psFitsWriteImageSimple (name, ro1->covariance->image, NULL); + + snprintf (name, 128, "conv.sample1.%03d.fits", nOut1); + psFitsWriteImageSimple (name, kernel->image, NULL); + + fprintf (stderr, "incov: %d,%d; kern: %d,%d, outcov: %d,%d\n", + ro1->covariance->image->numCols, ro1->covariance->image->numRows, + kernel->image->numCols, kernel->image->numRows, + cov->image->numCols, cov->image->numRows); + + nOut1 ++; + } psFree(kernel); } @@ -1407,4 +1495,23 @@ psKernelTruncate(kernel, covarFrac); covars->data[i] = psImageCovarianceCalculate(kernel, ro2->covariance); + if (0) { + char name[128]; + snprintf (name, 128, "covar.sample2.%03d.fits", nOut2); + psKernel *cov = covars->data[i]; + psFitsWriteImageSimple (name, cov->image, NULL); + + snprintf (name, 128, "incovar.sample2.%03d.fits", nOut2); + psFitsWriteImageSimple (name, ro2->covariance->image, NULL); + + snprintf (name, 128, "conv.sample2.%03d.fits", nOut2); + psFitsWriteImageSimple (name, kernel->image, NULL); + + fprintf (stderr, "incov: %d,%d; kern: %d,%d, outcov: %d,%d\n", + ro2->covariance->image->numCols, ro2->covariance->image->numRows, + kernel->image->numCols, kernel->image->numRows, + cov->image->numCols, cov->image->numRows); + + nOut2 ++; + } psFree(kernel); } @@ -1422,8 +1529,12 @@ psImageCovarianceSetThreads(oldThreads); - // Copy anything that wasn't convolved + // Copy anything that wasn't convolved (they may have been allocated though, so free them) switch (kernels->mode) { case PM_SUBTRACTION_MODE_1: if (out2) { + psFree(out2->image); + psFree(out2->variance); + psFree(out2->mask); + psFree(out2->covariance); out2->image = psMemIncrRefCounter(ro2->image); out2->variance = psMemIncrRefCounter(ro2->variance); @@ -1434,4 +1545,8 @@ case PM_SUBTRACTION_MODE_2: if (out1) { + psFree(out1->image); + psFree(out1->variance); + psFree(out1->mask); + psFree(out1->covariance); out1->image = psMemIncrRefCounter(ro1->image); out1->variance = psMemIncrRefCounter(ro1->variance); @@ -1479,2 +1594,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: trunk/psModules/src/imcombine/pmSubtraction.h =================================================================== --- trunk/psModules/src/imcombine/pmSubtraction.h (revision 30621) +++ trunk/psModules/src/imcombine/pmSubtraction.h (revision 30622) @@ -14,10 +14,9 @@ #define PM_SUBTRACTION_H -#include - -#include -#include -#include -#include +// #include +// #include +// #include +// #include +// #include // if we use the original ppSub implementation, we subtract a central delta-function for all @@ -30,18 +29,4 @@ /// @addtogroup imcombine Image Combinations /// @{ - -/// Mask values for the subtraction mask -typedef enum { - PM_SUBTRACTION_MASK_CLEAR = 0x00, // No masking - PM_SUBTRACTION_MASK_BAD_1 = 0x01, // Image 1 is bad - PM_SUBTRACTION_MASK_BAD_2 = 0x02, // Image 2 is bad - PM_SUBTRACTION_MASK_CONVOLVE_1 = 0x04, // If image 1 is convolved, would be poor or bad - PM_SUBTRACTION_MASK_CONVOLVE_2 = 0x08, // If image 2 is convolved, would be poor or bad - PM_SUBTRACTION_MASK_CONVOLVE_BAD_1 = 0x10, // If image 1 is convolved, would be bad - PM_SUBTRACTION_MASK_CONVOLVE_BAD_2 = 0x20, // If image 2 is convolved, would be bad - PM_SUBTRACTION_MASK_BORDER = 0x40, // Image border - PM_SUBTRACTION_MASK_REJ = 0x80, // Previously tried as a stamp, and rejected -} pmSubtractionMasks; - /// Number of terms in a polynomial @@ -70,8 +55,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 +156,6 @@ bool pmSubtractionSetFWHMs(float fwhm1, float fwhm2); +pmSubtractionQuality *pmSubtractionQualityAlloc(); + /// @} #endif Index: trunk/psModules/src/imcombine/pmSubtractionAnalysis.c =================================================================== --- trunk/psModules/src/imcombine/pmSubtractionAnalysis.c (revision 30621) +++ trunk/psModules/src/imcombine/pmSubtractionAnalysis.c (revision 30622) @@ -7,4 +7,5 @@ #include "pmFPA.h" +#include "pmSubtractionTypes.h" #include "pmSubtraction.h" #include "pmSubtractionKernels.h" Index: trunk/psModules/src/imcombine/pmSubtractionDeconvolve.c =================================================================== --- trunk/psModules/src/imcombine/pmSubtractionDeconvolve.c (revision 30621) +++ trunk/psModules/src/imcombine/pmSubtractionDeconvolve.c (revision 30622) @@ -10,4 +10,5 @@ #include "pmFPA.h" +#include "pmSubtractionTypes.h" #include "pmSubtractionKernels.h" #include "pmSubtractionDeconvolve.h" Index: trunk/psModules/src/imcombine/pmSubtractionEquation.c =================================================================== --- trunk/psModules/src/imcombine/pmSubtractionEquation.c (revision 30621) +++ trunk/psModules/src/imcombine/pmSubtractionEquation.c (revision 30622) @@ -8,4 +8,7 @@ #include "pmErrorCodes.h" +#include "pmVisual.h" +#include "pmFPA.h" +#include "pmSubtractionTypes.h" #include "pmSubtraction.h" #include "pmSubtractionKernels.h" @@ -16,14 +19,16 @@ #include "pmSubtractionVisual.h" -//#define TESTING // TESTING output for debugging; may not work with threads! - -//#define USE_WEIGHT // Include weight (1/variance) in equation? - -// XXX TEST: -//# define USE_WINDOW // window to avoid neighbor contamination +//# define TESTING // TESTING output for debugging; may not work with threads! +# define USE_WEIGHT // Include weight (1/variance) in equation? +# define USE_WINDOW // window to avoid neighbor contamination + +/* I believe we want to apply the WEIGHT to the chisq portions of the calculation (but not the WINDOW), + * and the WINDOW to the moments portiosn of the calculations (but not the WEIGHT) + * + */ # define PENALTY false # define MOMENTS (!PENALTY) -# define MOMENTS_PENALTY_SCALE 2e-5 // XXX this value is not completely arbitrary, but I don't understand why it needs to be this value... +# define MOMENTS_PENALTY_SCALE 20 // up-weight the moments somewhat ////////////////////////////////////////////////////////////////////////////////////////////////////////////// @@ -107,5 +112,6 @@ cc *= weight->kernel[y][x]; } - if (window) { + // XXX NOTE: do NOT apply the window to the chisq portions of the calculation + if (false && window) { cc *= window->kernel[y][x]; } @@ -138,5 +144,6 @@ rc *= wtVal; } - if (window) { + // XXX NOTE: do NOT apply the window to the chisq portions of the calculation + if (false && window) { float winVal = window->kernel[y][x]; ic *= winVal; @@ -173,4 +180,5 @@ } +# define RENORM_BY_FLUX 0 // Calculate the least-squares matrix and vector for dual convolution @@ -268,32 +276,34 @@ for (int y = - footprint; y <= footprint; y++) { for (int x = - footprint; x <= footprint; x++) { + + // XXX NOTE: clipping low S/N pixels does not seem to work very well + if (false && weight) { + float i1 = image1->kernel[y][x]; + float i2 = image2->kernel[y][x]; + float sn = (i1 + i2) / sqrt (weight->kernel[y][x]); + if (sn < 0.5) continue; + } + double aa = iConv1->kernel[y][x] * jConv1->kernel[y][x] * PS_SQR(normValue); double bb = iConv2->kernel[y][x] * jConv2->kernel[y][x]; double ab = iConv1->kernel[y][x] * jConv2->kernel[y][x] * normValue; - if (weight) { - float wtVal = weight->kernel[y][x]; - aa *= wtVal; - bb *= wtVal; - ab *= wtVal; - } - if (window) { - float wtVal = window->kernel[y][x]; - aa *= wtVal; - bb *= wtVal; - ab *= wtVal; - } - sumAA += aa; - sumBB += bb; - sumAB += ab; + + float wtVal = (weight) ? weight->kernel[y][x] : 1.0; + sumAA += wtVal*aa; + sumBB += wtVal*bb; + sumAB += wtVal*ab; if (MOMENTS) { - MxxAA += x*x*aa; - MyyAA += y*y*aa; - MxxBB += x*x*bb; - MyyBB += y*y*bb; + float winVal = (window) ? window->kernel[y][x] : 1.0; + MxxAA += winVal*x*x*aa; + MyyAA += winVal*y*y*aa; + MxxBB += winVal*x*x*bb; + MyyBB += winVal*y*y*bb; } } } + // XXX does normSquare1,2 mess up the relative scaling? + // XXX no: normSquare1,2 is the sum of the flux^2 for the source if (MOMENTS) { MxxAA /= stamp->normSquare1 * PS_SQR(normValue); @@ -305,7 +315,9 @@ // XXX this makes the Chisq portion independent of the normalization and star flux // but may be mis-scaling between stars of different fluxes +# if (RENORM_BY_FLUX) sumAA /= PS_SQR(stamp->normI2); sumAB /= PS_SQR(stamp->normI2); sumBB /= PS_SQR(stamp->normI2); +# endif // fprintf (stderr, "i,j : %d %d : M(xx,yy)(AA,BB) : %f %f %f %f\n", i, j, MxxAA, MyyAA, MxxBB, MyyBB); @@ -328,5 +340,4 @@ matrix->data.F64[iIndex][jIndex] += kernels->penalty * MxxAA * MOMENTS_PENALTY_SCALE; matrix->data.F64[iIndex][jIndex] += kernels->penalty * MyyAA * MOMENTS_PENALTY_SCALE; - matrix->data.F64[jIndex][iIndex] += kernels->penalty * MxxAA * MOMENTS_PENALTY_SCALE; matrix->data.F64[jIndex][iIndex] += kernels->penalty * MyyAA * MOMENTS_PENALTY_SCALE; @@ -334,5 +345,4 @@ matrix->data.F64[iIndex + numParams][jIndex + numParams] += kernels->penalty * MxxBB * MOMENTS_PENALTY_SCALE; matrix->data.F64[iIndex + numParams][jIndex + numParams] += kernels->penalty * MyyBB * MOMENTS_PENALTY_SCALE; - matrix->data.F64[jIndex + numParams][iIndex + numParams] += kernels->penalty * MxxBB * MOMENTS_PENALTY_SCALE; matrix->data.F64[jIndex + numParams][iIndex + numParams] += kernels->penalty * MyyBB * MOMENTS_PENALTY_SCALE; @@ -354,5 +364,6 @@ ab *= weight->kernel[y][x]; } - if (window) { + // XXX NOTE: do NOT apply the window to the chisq portions of the calculation + if (false && window) { ab *= window->kernel[y][x]; } @@ -363,5 +374,7 @@ // XXX this makes the Chisq portion independent of the normalization and star flux // but may be mis-scaling between stars of different fluxes +# if (RENORM_BY_FLUX) sumAB /= PS_SQR(stamp->normI2); +# endif // Spatial variation of kernel coeffs @@ -391,4 +404,10 @@ float i2 = image2->kernel[y][x]; + // XXX NOTE: clipping low S/N pixels does not seem to work very well + if (false && weight) { + float sn = (i1 + i2) / sqrt (weight->kernel[y][x]); + if (sn < 0.5) continue; + } + double ai2 = a * i2 * normValue; double bi2 = b * i2; @@ -396,28 +415,16 @@ double bi1 = b * i1 * normValue; - if (weight) { - float wtVal = weight->kernel[y][x]; - ai2 *= wtVal; - bi2 *= wtVal; - ai1 *= wtVal; - bi1 *= wtVal; - } - if (window) { - float wtVal = window->kernel[y][x]; - ai2 *= wtVal; - bi2 *= wtVal; - ai1 *= wtVal; - bi1 *= wtVal; - } - sumAI2 += ai2; - sumBI2 += bi2; - sumAI1 += ai1; - sumBI1 += bi1; + float wtVal = (weight) ? weight->kernel[y][x] : 1.0; + sumAI2 += wtVal*ai2; + sumBI2 += wtVal*bi2; + sumAI1 += wtVal*ai1; + sumBI1 += wtVal*bi1; if (MOMENTS) { - MxxAI1 += x*x*ai1; - MyyAI1 += y*y*ai1; - MxxBI2 += x*x*bi2; - MyyBI2 += y*y*bi2; + float winVal = (window) ? window->kernel[y][x] : 1.0; + MxxAI1 += winVal*x*x*ai1; + MyyAI1 += winVal*y*y*ai1; + MxxBI2 += winVal*x*x*bi2; + MyyBI2 += winVal*y*y*bi2; } } @@ -435,8 +442,10 @@ // XXX this makes the Chisq portion independent of the normalization and star flux // but may be mis-scaling between stars of different fluxes +# if (RENORM_BY_FLUX) sumAI1 /= PS_SQR(stamp->normI2); sumBI1 /= PS_SQR(stamp->normI2); sumAI2 /= PS_SQR(stamp->normI2); sumBI2 /= PS_SQR(stamp->normI2); +# endif // Spatial variation @@ -788,5 +797,5 @@ stamp->normSquare2 = normSquare2; - psLogMsg ("psModules.imcombine", PS_LOG_DETAIL, "normValue: %f %f %f (%f %f)\n", normI1, normI2, stamp->norm, normSquare1, normSquare2); + // psLogMsg ("psModules.imcombine", PS_LOG_DETAIL, "normValue: %f %f %f (%f %f)\n", normI1, normI2, stamp->norm, normSquare1, normSquare2); return true; @@ -800,11 +809,9 @@ pmSubtractionKernels *kernels = job->args->data[1]; // Kernels int index = PS_SCALAR_VALUE(job->args->data[2], S32); // Stamp index - pmSubtractionEquationCalculationMode mode = PS_SCALAR_VALUE(job->args->data[3], S32); // calculation model - - return pmSubtractionCalculateEquationStamp(stamps, kernels, index, mode); -} - -bool pmSubtractionCalculateEquationStamp(pmSubtractionStampList *stamps, pmSubtractionKernels *kernels, - int index, const pmSubtractionEquationCalculationMode mode) + + return pmSubtractionCalculateEquationStamp(stamps, kernels, index); +} + +bool pmSubtractionCalculateEquationStamp(pmSubtractionStampList *stamps, pmSubtractionKernels *kernels, int index) { PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(stamps, false); @@ -833,40 +840,14 @@ psAssert(stamp->status == PM_SUBTRACTION_STAMP_CALCULATE, "We only operate on stamps with this state."); - // Generate convolutions: these are generated once and saved - if (!pmSubtractionConvolveStamp(stamp, kernels, footprint)) { - psError(psErrorCodeLast(), false, "Unable to convolve stamp %d.", index); - return NULL; - } - -#ifdef TESTING - for (int j = 0; j < numKernels; 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 - - // XXX visualize the set of convolved stamps - - psImage *polyValues = p_pmSubtractionPolynomial(NULL, spatialOrder, - stamp->xNorm, stamp->yNorm); // Polynomial terms - - bool new = stamp->vector ? false : true; // Is this a new run? + psImage *polyValues = p_pmSubtractionPolynomial(NULL, spatialOrder, stamp->xNorm, stamp->yNorm); // Polynomial terms + + // Is this a new run? Have we allocated the correct sized vector/matrix? + bool new = stamp->vector ? false : true; + if (!new && (stamp->vector->n != numParams)) { + psFree (stamp->vector); + psFree (stamp->matrix); + new = true; + } + if (new) { stamp->matrix = psImageAlloc(numParams, numParams, PS_TYPE_F64); @@ -941,6 +922,5 @@ } -bool pmSubtractionCalculateEquation(pmSubtractionStampList *stamps, pmSubtractionKernels *kernels, - const pmSubtractionEquationCalculationMode mode) +bool pmSubtractionCalculateEquation(pmSubtractionStampList *stamps, pmSubtractionKernels *kernels) { PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(stamps, false); @@ -953,4 +933,5 @@ for (int i = 0; i < stamps->num; i++) { pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest + if (stamp->status != PM_SUBTRACTION_STAMP_CALCULATE) { continue; @@ -969,10 +950,9 @@ psArrayAdd(job->args, 1, (pmSubtractionKernels*)kernels); // Casting away const to put on array PS_ARRAY_ADD_SCALAR(job->args, i, PS_TYPE_S32); - PS_ARRAY_ADD_SCALAR(job->args, mode, PS_TYPE_S32); if (!psThreadJobAddPending(job)) { return false; } } else { - pmSubtractionCalculateEquationStamp(stamps, kernels, i, mode); + pmSubtractionCalculateEquationStamp(stamps, kernels, i); } } @@ -983,11 +963,9 @@ } - pmSubtractionVisualPlotLeastSquares(stamps); pmSubtractionVisualShowKernels((pmSubtractionKernels *)kernels); pmSubtractionVisualShowBasis(stamps); - - psLogMsg("psModules.imcombine", PS_LOG_INFO, "Calculate equation: %f sec", - psTimerClear("pmSubtractionCalculateEquation")); - + pmSubtractionVisualPlotLeastSquares(stamps); + + psLogMsg("psModules.imcombine", PS_LOG_INFO, "Calculate equation: %f sec", psTimerClear("pmSubtractionCalculateEquation")); return true; @@ -998,27 +976,10 @@ bool psFitsWriteImageSimple (char *filename, psImage *image, psMetadata *header); -psImage *p_pmSubSolve_wUt (psVector *w, psImage *U); -psImage *p_pmSubSolve_VwUt (psImage *V, psImage *wUt); - -bool p_pmSubSolve_SetWeights (psVector *wApply, psVector *w, psVector *wMask); - -bool p_pmSubSolve_UtB (psVector **UtB, psImage *U, psVector *B); -bool p_pmSubSolve_wUtB (psVector **wUtB, psVector *w, psVector *UtB); -bool p_pmSubSolve_VwUtB (psVector **VwUtB, psImage *V, psVector *wUtB); - -bool p_pmSubSolve_Ax (psVector **B, psImage *A, psVector *x); -bool p_pmSubSolve_VdV (double *value, psVector *x, psVector *y); -bool p_pmSubSolve_y2 (double *y2, pmSubtractionKernels *kernels, const pmSubtractionStampList *stamps); - -psImage *p_pmSubSolve_Xvar (psImage *V, psVector *w); - -double p_pmSubSolve_ChiSquare (pmSubtractionKernels *kernels, const pmSubtractionStampList *stamps); - -bool pmSubtractionSolveEquation(pmSubtractionKernels *kernels, - const pmSubtractionStampList *stamps, - const pmSubtractionEquationCalculationMode mode) +bool pmSubtractionSolveEquation(pmSubtractionKernels *kernels, const pmSubtractionStampList *stamps) { PM_ASSERT_SUBTRACTION_KERNELS_NON_NULL(kernels, false); PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(stamps, false); + + psTimerStart("pmSubtractionSolveEquation"); // Check inputs @@ -1042,4 +1003,7 @@ } + if (stamp->vector->n != numParams) { + fprintf (stderr, "mismatch length\n"); + } PS_ASSERT_VECTOR_NON_NULL(stamp->vector, false); PS_ASSERT_VECTOR_SIZE(stamp->vector, (long)numParams, false); @@ -1080,4 +1044,6 @@ } + pmSubtractionVisualPlotLeastSquaresResid(stamps, sumMatrix, numStamps); + #if 0 psImage *save = psImageCopy(NULL, sumMatrix, PS_TYPE_F32); @@ -1087,4 +1053,15 @@ #endif + // XXX TEST : print the matrix & vector + if (0) { + for (int iy = 0; iy < sumMatrix->numRows; iy++) { + for (int ix = 0; ix < sumMatrix->numCols; ix++) { + fprintf (stderr, "%e ", sumMatrix->data.F64[iy][ix]); + } + fprintf (stderr, " : %e\n", sumVector->data.F64[iy]); + } + } + + psImage *invMatrix = NULL; psVector *solution = NULL; // Solution to equation! solution = psVectorAlloc(numParams, PS_TYPE_F64); @@ -1094,30 +1071,38 @@ // solution = psMatrixSolveSVD(solution, sumMatrix, sumVector, NAN); // SINGLE solution - if (1) { - solution = psMatrixSolveSVD(solution, sumMatrix, sumVector, NAN); - } else { - psVector *PERM = NULL; - psImage *LU = psMatrixLUDecomposition(NULL, &PERM, sumMatrix); - solution = psMatrixLUSolution(solution, LU, sumVector, PERM); - psFree (LU); - psFree (PERM); - } +# if (1) + solution = psMatrixSolveSVD(solution, sumMatrix, sumVector, 1e-10); + invMatrix = psMatrixInvert(NULL, sumMatrix, NULL); +# endif +# if (0) + psMatrixLUSolve(sumMatrixLU, sumVector); + solution = psMemIncrRefCounter(sumVector); + invMatrix = psMemIncrRefCounter(sumMatrix); +# endif +# if (0) + psMatrixGJSolve(sumMatrix, sumVector); + invMatrix = psMemIncrRefCounter(sumMatrix); + solution = psMemIncrRefCounter(sumVector); +# endif # if (0) for (int i = 0; i < solution->n; i++) { - psLogMsg("psModules.imcombine", PS_LOG_DETAIL, "Single solution %d: %lf\n", i, solution->data.F64[i]); + psLogMsg("psModules.imcombine", PS_LOG_DETAIL, "Single solution %d: %lf +/- %lf\n", i, solution->data.F64[i], sqrt(fabs(invMatrix->data.F64[i][i]))); } # endif - if (!kernels->solution1) { - kernels->solution1 = psVectorAlloc(sumVector->n + 2, PS_TYPE_F64); // 1 for norm, 1 for bg - psVectorInit(kernels->solution1, 0.0); - } + // ensure we have a solution vector of the right size + kernels->solution1 = psVectorRecycle(kernels->solution1, sumVector->n + 2, PS_TYPE_F64); // 1 for norm, 1 for bg + kernels->solution1err = psVectorRecycle(kernels->solution1err, sumVector->n + 2, PS_TYPE_F64); // 1 for norm, 1 for bg + psVectorInit(kernels->solution1, 0.0); + psVectorInit(kernels->solution1err, 0.0); int numKernels = kernels->num; int spatialOrder = kernels->spatialOrder; // Order of spatial variation int numPoly = PM_SUBTRACTION_POLYTERMS(spatialOrder); // Number of polynomial terms + for (int i = 0; i < numKernels * numPoly; i++) { kernels->solution1->data.F64[i] = solution->data.F64[i]; + kernels->solution1err->data.F64[i] = sqrt(invMatrix->data.F64[i][i]); } @@ -1131,4 +1116,5 @@ psFree(sumVector); psFree(sumMatrix); + psFree(invMatrix); } else { @@ -1160,4 +1146,6 @@ } + pmSubtractionVisualPlotLeastSquaresResid(stamps, sumMatrix, numStamps); + #if 0 psImage *save = psImageCopy(NULL, sumMatrix, PS_TYPE_F32); @@ -1171,4 +1159,15 @@ } + // XXX TEST : print the matrix & vector + if (0) { + for (int iy = 0; iy < sumMatrix->numRows; iy++) { + for (int ix = 0; ix < sumMatrix->numCols; ix++) { + fprintf (stderr, "%e ", sumMatrix->data.F64[iy][ix]); + } + fprintf (stderr, " : %e\n", sumVector->data.F64[iy]); + } + } + + psImage *invMatrix = NULL; psVector *solution = NULL; // Solution to equation! solution = psVectorAlloc(numParams, PS_TYPE_F64); @@ -1176,20 +1175,36 @@ // DUAL solution - solution = psMatrixSolveSVD(solution, sumMatrix, sumVector, NAN); +# if (1) + solution = psMatrixSolveSVD(solution, sumMatrix, sumVector, 1e-10); + invMatrix = psMatrixInvert(NULL, sumMatrix, NULL); +# endif +# if (0) + psMatrixLUSolve(sumMatrix, sumVector); + solution = psMemIncrRefCounter(sumVector); + invMatrix = psMemIncrRefCounter(sumMatrix); +# endif #if (0) for (int i = 0; i < solution->n; i++) { - fprintf(stderr, "Dual solution %d: %lf\n", i, solution->data.F64[i]); + fprintf(stderr, "Dual solution %d: %lf +/- %lf\n", i, solution->data.F64[i], sqrt(invMatrix->data.F64[i][i])); } #endif - if (!kernels->solution1) { - kernels->solution1 = psVectorAlloc(numSolution1 + 2, PS_TYPE_F64); - psVectorInit (kernels->solution1, 0.0); - } - if (!kernels->solution2) { - kernels->solution2 = psVectorAlloc(numSolution2, PS_TYPE_F64); - psVectorInit (kernels->solution2, 0.0); - } + // XXX TEST: manually set the coeffs to a desired solution + // solution->data.F64[0] = +1.826; + // solution->data.F64[1] = -0.115; + // solution->data.F64[2] = 0.0; + // solution->data.F64[3] = 0.0; + + // ensure we have solution vectors of the right size + kernels->solution1 = psVectorRecycle(kernels->solution1, numSolution1 + 2, PS_TYPE_F64); // 1 for norm, 1 for bg + kernels->solution1err = psVectorRecycle(kernels->solution1err, numSolution1 + 2, PS_TYPE_F64); // 1 for norm, 1 for bg + kernels->solution2 = psVectorRecycle(kernels->solution2, numSolution2, PS_TYPE_F64); // 1 for norm, 1 for bg + kernels->solution2err = psVectorRecycle(kernels->solution2err, numSolution2, PS_TYPE_F64); // 1 for norm, 1 for bg + + psVectorInit(kernels->solution1, 0.0); + psVectorInit(kernels->solution1err, 0.0); + psVectorInit(kernels->solution2, 0.0); + psVectorInit(kernels->solution2err, 0.0); // for DUAL convolution analysis, we apply the normalization to I1 as follows: @@ -1205,4 +1220,8 @@ kernels->solution1->data.F64[i] = solution->data.F64[i] * stamps->normValue; kernels->solution2->data.F64[i] = solution->data.F64[i + numSolution1]; + + kernels->solution1err->data.F64[i] = sqrt(invMatrix->data.F64[i][i]) * stamps->normValue; + int i2 = i + numSolution1; + kernels->solution2err->data.F64[i] = sqrt(invMatrix->data.F64[i2][i2]); } @@ -1213,7 +1232,8 @@ kernels->solution1->data.F64[bgIndex] = 0.0; + psFree(solution); + psFree(sumVector); psFree(sumMatrix); - psFree(sumVector); - psFree(solution); + psFree(invMatrix); } @@ -1224,12 +1244,14 @@ if (psTraceGetLevel("psModules.imcombine") >= 7) { for (int i = 0; i < kernels->solution1->n; i++) { - psTrace("psModules.imcombine", 7, "Solution 1 %d: %f\n", i, kernels->solution1->data.F64[i]); + psTrace("psModules.imcombine", 7, "Solution 1 %d: %f +/- %f\n", i, kernels->solution1->data.F64[i], kernels->solution1err->data.F64[i]); } if (kernels->mode == PM_SUBTRACTION_MODE_DUAL) { for (int i = 0; i < kernels->solution2->n; i++) { - psTrace("psModules.imcombine", 7, "Solution 2 %d: %f\n", i, kernels->solution2->data.F64[i]); + psTrace("psModules.imcombine", 7, "Solution 2 %d: %f +/- %f\n", i, kernels->solution2->data.F64[i], kernels->solution2err->data.F64[i]); } } } + + psLogMsg("psModules.imcombine", PS_LOG_INFO, "Solve equation: %f sec", psTimerClear("pmSubtractionSolveEquation")); // pmSubtractionVisualPlotLeastSquares((pmSubtractionStampList *) stamps); //casting away const @@ -1283,10 +1305,417 @@ } -psVector *pmSubtractionCalculateDeviations(pmSubtractionStampList *stamps, - pmSubtractionKernels *kernels) +// given the convolved image(s) and the residual image, calculate the second moment(s) and the chisq +bool pmSubtractionChisqStats(psVector *fluxesVector, psVector *chisqDVector, psVector *chisqRVector, psVector *momentVector, psVector *stampMask, psKernel *convolved1, psKernel *convolved2, psKernel *difference, psKernel *residual, psKernel *weight, psKernel *window) { + +# ifndef USE_WEIGHT + psAssert(weight == NULL, "impossible!"); +# endif +# ifndef USE_WINDOW + psAssert(window == NULL, "impossible!"); +# endif + + int npix = 0; + float chisqR = 0; + float chisqD = 0; + + // get the chisq + for (int y = residual->yMin; y <= residual->yMax; y++) { + for (int x = residual->xMin; x <= residual->xMax; x++) { + float valueR = PS_SQR(residual->kernel[y][x]); + if (weight) { + valueR *= weight->kernel[y][x]; + } + // XXX NOTE: do NOT apply the window to the chisq portions of the calculation (that would bias the chisq) + chisqR += valueR; + + float valueD = PS_SQR(difference->kernel[y][x]); + if (weight) { + valueD *= weight->kernel[y][x]; + } + chisqD += valueD; + npix ++; + } + } + psVectorAppend(chisqRVector, chisqR / npix); + psVectorAppend(chisqDVector, chisqD / npix); + + float value1 = 0; + float value2 = 0; + float flux2 = 0; + float fluxX = 0; + float fluxY = 0; + float fluxX2 = 0; + float fluxY2 = 0; + + float fluxC1 = 0; + float fluxC2 = 0; + + float moment = 0; + + // get the moments from convolved1 + if (convolved1) { + for (int y = residual->yMin; y <= residual->yMax; y++) { + for (int x = residual->xMin; x <= residual->xMax; x++) { + value1 = convolved1->kernel[y][x]; + value2 = PS_SQR(value1); + + if (window) { + value1 *= window->kernel[y][x]; + value2 *= window->kernel[y][x]; + } + + fluxC1 += value1; + flux2 += value2; + fluxX += x * value2; + fluxY += y * value2; + // fluxX2 += PS_SQR(x) * value2; + // fluxY2 += PS_SQR(y) * value2; + fluxX2 += PS_SQR(x) * value1; + fluxY2 += PS_SQR(y) * value1; + } + } + // float Mx = fluxX / flux2; + // float My = fluxY / flux2; + // float Mxx = fluxX2 / flux2; + // float Myy = fluxY2 / flux2; + float Mxx = fluxX2 / fluxC1; + float Myy = fluxY2 / fluxC1; + + // fprintf (stderr, "conv1, flux2: %f, Mx: %f, My: %f, Mxx: %f, Myy: %f, chisq: %f, npix: %d\n", flux2, Mx, My, Mxx, Myy, chisq, npix); + moment += Mxx + Myy; + } + + // get the moments from convolved1 + if (convolved2) { + for (int y = residual->yMin; y <= residual->yMax; y++) { + for (int x = residual->xMin; x <= residual->xMax; x++) { + value1 = convolved2->kernel[y][x]; + value2 = PS_SQR(value1); + + // XXX NOTE: do NOT apply the weight to the moments calculation + if (false && weight) { + value2 *= weight->kernel[y][x]; + } + if (window) { + value1 *= window->kernel[y][x]; + value2 *= window->kernel[y][x]; + } + + fluxC2 += value1; + flux2 += value2; + fluxX += x * value2; + fluxY += y * value2; + // fluxX2 += PS_SQR(x) * value2; + // fluxY2 += PS_SQR(y) * value2; + fluxX2 += PS_SQR(x) * value1; + fluxY2 += PS_SQR(y) * value1; + } + } + // float Mx = fluxX / flux2; + // float My = fluxY / flux2; + // float Mxx = fluxX2 / flux2; + // float Myy = fluxY2 / flux2; + float Mxx = fluxX2 / fluxC2; + float Myy = fluxY2 / fluxC2; + + // fprintf (stderr, "conv2, flux2: %f, Mx: %f, My: %f, Mxx: %f, Myy: %f, chisq: %f, npix: %d\n", flux2, Mx, My, Mxx, Myy, chisq, npix); + moment += Mxx + Myy; + } + + float flux = fluxC1 + fluxC2; + + if (convolved1 && convolved2) { + moment *= 0.5; + flux *= 0.5; + } + psVectorAppend(momentVector, moment); + psVectorAppend(fluxesVector, flux); + + // check that the last appended values are ok: + int Nelem = fluxesVector->n - 1; + bool valid = true; + valid &= isfinite(chisqRVector->data.F32[Nelem]); + valid &= isfinite(fluxesVector->data.F32[Nelem]); + valid &= isfinite(momentVector->data.F32[Nelem]); + if (valid) { + psVectorAppend(stampMask, 0); + } else { + psVectorAppend(stampMask, 0x02); + } + return true; +} + +bool pmSubtractionCalculateChisqAndMoments(pmSubtractionQuality **bestMatch, + pmSubtractionStampList *stamps, + pmSubtractionKernels *kernels) { PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(stamps, NULL); PM_ASSERT_SUBTRACTION_KERNELS_NON_NULL(kernels, NULL); PM_ASSERT_SUBTRACTION_KERNELS_SOLUTION(kernels, NULL); + + psTimerStart("pmSubtractionCalculateChisqAndMoments"); + + // XXX need to save these somewhere + psVector *fluxes = psVectorAllocEmpty(stamps->num, PS_TYPE_F32); + psVector *chisqD = psVectorAllocEmpty(stamps->num, PS_TYPE_F32); + psVector *chisqR = psVectorAllocEmpty(stamps->num, PS_TYPE_F32); + psVector *moments = psVectorAllocEmpty(stamps->num, PS_TYPE_F32); + psVector *stampMask = psVectorAllocEmpty(stamps->num, PS_TYPE_VECTOR_MASK); + + int footprint = stamps->footprint; // Half-size of stamps + int numKernels = kernels->num; // Number of kernels + + psImage *polyValues = NULL; // Polynomial values + + // storage for the image (convolved2 is not used in SINGLE mode) + psKernel *residual = psKernelAlloc(-footprint, footprint, -footprint, footprint); // Residual image + psKernel *difference = psKernelAlloc(-footprint, footprint, -footprint, footprint); // Residual image + psKernel *convolved1 = psKernelAlloc(-footprint, footprint, -footprint, footprint); // Residual image + psKernel *convolved2 = psKernelAlloc(-footprint, footprint, -footprint, footprint); // Residual image + + int nGood = 0; + for (int i = 0; i < stamps->num; i++) { + pmSubtractionStamp *stamp = stamps->stamps->data[i]; // The stamp of interest + if (stamp->status != PM_SUBTRACTION_STAMP_USED) { + // mark this stamp as unused (note that we have to append NANs to the other vectors to keep the lengths in sync) + psVectorAppend(moments, NAN); + psVectorAppend(fluxes, NAN); + psVectorAppend(chisqD, NAN); + psVectorAppend(chisqR, NAN); + psVectorAppend(stampMask, 0x01); + continue; + } + nGood ++; + + // Calculate coefficients of the kernel basis functions + polyValues = p_pmSubtractionPolynomial(polyValues, kernels->spatialOrder, stamp->xNorm, stamp->yNorm); + double norm = p_pmSubtractionSolutionNorm(kernels); // Normalisation + double background = p_pmSubtractionSolutionBackground(kernels, polyValues);// Difference in background + + // Calculate residuals + psImageInit(residual->image, 0.0); + psImageInit(difference->image, 0.0); + + psKernel *weight = NULL; + psKernel *window = NULL; + +#ifdef USE_WEIGHT + weight = stamp->weight; +#endif +#ifdef USE_WINDOW + window = stamps->window; +#endif + + if (kernels->mode != PM_SUBTRACTION_MODE_DUAL) { + + // the single-direction psf match code attempts to find the kernel such that: + // source * kernel = target. we need to assign 'source' and 'target' correctly + // depending on which of image1 or image2 we asked to be convolved. + + psKernel *target; // Target postage stamp (convolve source to match the target) + psKernel *source; // Source postage stamp (convolve source to match the target) + psArray *convolutions; // Convolution postage stamps for each kernel basis function + + // init the accumulation image + psImageInit(convolved1->image, 0.0); + + switch (kernels->mode) { + case PM_SUBTRACTION_MODE_1: + target = stamp->image2; + source = stamp->image1; + convolutions = stamp->convolutions1; + break; + case PM_SUBTRACTION_MODE_2: + target = stamp->image1; + source = stamp->image2; + convolutions = stamp->convolutions2; + break; + default: + psAbort("Unsupported subtraction mode: %x", kernels->mode); + } + + // generate the convolved source image (sum over kernels) + for (int j = 0; j < numKernels; j++) { + psKernel *convolution = convolutions->data[j]; // Convolution + double coefficient = p_pmSubtractionSolutionCoeff(kernels, polyValues, j, false); // Coefficient + for (int y = - footprint; y <= footprint; y++) { + for (int x = - footprint; x <= footprint; x++) { + convolved1->kernel[y][x] += convolution->kernel[y][x] * coefficient; + } + } + } + + // Generate the difference, residual, and convolved source images. Note the we + // accumulate the convolution of (A-B), so we need to replace it to generate the + // images of the convolved source image. + for (int y = - footprint; y <= footprint; y++) { + for (int x = - footprint; x <= footprint; x++) { + difference->kernel[y][x] = target->kernel[y][x] - source->kernel[y][x] * norm - background; + residual->kernel[y][x] = difference->kernel[y][x] - convolved1->kernel[y][x]; + convolved1->kernel[y][x] += source->kernel[y][x] * norm; + } + } + + // XXX if we want to have a weight and window, we'll need to pass through to here + pmSubtractionChisqStats(fluxes, chisqD, chisqR, moments, stampMask, convolved1, NULL, difference, residual, weight, window); + + } else { + + // Dual convolution + psArray *convolutions1 = stamp->convolutions1; // Convolutions of the first image + psArray *convolutions2 = stamp->convolutions2; // Convolutions of the second image + psKernel *image1 = stamp->image1; // The first image + psKernel *image2 = stamp->image2; // The second image + + // init the accumulation images + psImageInit(convolved1->image, 0.0); + psImageInit(convolved2->image, 0.0); + + for (int j = 0; j < numKernels; j++) { + psKernel *conv1 = convolutions1->data[j]; // Convolution of first image + psKernel *conv2 = convolutions2->data[j]; // Convolution of second image + double coeff1 = p_pmSubtractionSolutionCoeff(kernels, polyValues, j, false); // Coefficient 1 + double coeff2 = p_pmSubtractionSolutionCoeff(kernels, polyValues, j, true); // Coefficient 2 + + for (int y = - footprint; y <= footprint; y++) { + for (int x = - footprint; x <= footprint; x++) { + // NOTE sign for coeff2 + convolved1->kernel[y][x] += +conv1->kernel[y][x] * coeff1; + convolved2->kernel[y][x] += -conv2->kernel[y][x] * coeff2; + } + } + } + + // Generate the difference, residual, and convolved source images. Note the we + // accumulate the convolutions of (A-B), so we need to replace (A or B) to generate + // the images of the convolved source images. + for (int y = - footprint; y <= footprint; y++) { + for (int x = - footprint; x <= footprint; x++) { + difference->kernel[y][x] = image2->kernel[y][x] - image1->kernel[y][x] * norm - background; + residual->kernel[y][x] = difference->kernel[y][x] + convolved2->kernel[y][x] - convolved1->kernel[y][x]; + convolved1->kernel[y][x] += image1->kernel[y][x] * norm; + convolved2->kernel[y][x] += image2->kernel[y][x]; + } + } + + if (0) { + psFitsWriteImageSimple("conv1.fits", convolved1->image, NULL); + psFitsWriteImageSimple("conv2.fits", convolved2->image, NULL); + psFitsWriteImageSimple("resid.fits", residual->image, NULL); + pmVisualAskUser(NULL); + } + + pmSubtractionChisqStats(fluxes, chisqD, chisqR, moments, stampMask, convolved1, convolved2, difference, residual, weight, window); + } + } + + // find the mean chisq and mean moment + psStats *stats = psStatsAlloc(PS_STAT_SAMPLE_MEAN); + psVectorStats (stats, chisqD, NULL, stampMask, 0xff); + float chisqDValue = stats->sampleMean; + + psStatsInit(stats); + psVectorStats (stats, chisqR, NULL, stampMask, 0xff); + float chisqRValue = stats->sampleMean; + + psStatsInit(stats); + psVectorStats (stats, moments, NULL, stampMask, 0xff); + float momentValue = stats->sampleMean; + + double sumKernel1 = 0.0, sumKernel2 = 0.0; // Sum of the kernel + + // calculate the variance contribution from this smoothing kernel + psKernel *modelKernel = pmSubtractionKernel(kernels, 0.0, 0.0, false); + for (int y = modelKernel->yMin; y <= modelKernel->yMax; y++) { + for (int x = modelKernel->xMin; x <= modelKernel->xMax; x++) { + if (!isfinite(modelKernel->kernel[y][x])) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Non-finite covariance matrix element in kernel at %d,%d", x, y); + return NULL; + } + sumKernel1 += PS_SQR(modelKernel->kernel[y][x]); + } + } + psFree (modelKernel); + + if (kernels->mode == PM_SUBTRACTION_MODE_DUAL) { + psKernel *modelKernel = pmSubtractionKernel(kernels, 0.0, 0.0, true); + for (int y = modelKernel->yMin; y <= modelKernel->yMax; y++) { + for (int x = modelKernel->xMin; x <= modelKernel->xMax; x++) { + if (!isfinite(modelKernel->kernel[y][x])) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Non-finite covariance matrix element in kernel at %d,%d", x, y); + return NULL; + } + sumKernel2 += PS_SQR(modelKernel->kernel[y][x]); + } + } + psFree (modelKernel); + } else { + sumKernel2 = 1.0; + } + + // if we modify the chisq value by the (sumKernel1 + sumKernel2), we account for the + // smoothing coming from larger kernels adding additional spatial fit terms should be + // penalized by increasing the score somewhat. the 0.01 value is not well-chosen. + float orderFactor = 0.01 * kernels->spatialOrder; + float score = 2.0 * chisqRValue / (sumKernel1 + sumKernel2) + orderFactor; + psLogMsg("psModules.imcombine", PS_LOG_INFO, "chisq: %6.3f, chisqD: %6.3f, moment: %6.3f, sumKernel_1: %6.3f, sumKernel_2, score: %6.3f: %6.3f\n", chisqRValue, chisqDValue, momentValue, sumKernel1, sumKernel2, score); + + // save this result if it is the first or the best (skip if bestMatch is NULL) + if (bestMatch) { + pmSubtractionQuality *match = *bestMatch; + bool keep = false; + if (match == NULL) { + *bestMatch = match = pmSubtractionQualityAlloc(); + keep = true; + } else { + if (score < match->score) { + psFree(match->fluxes); + psFree(match->chisq); + psFree(match->moments); + psFree(match->stampMask); + keep = true; + } + } + if (keep) { + psLogMsg("psModules.imcombine", PS_LOG_INFO, "keeping order: %d, mode: %d, score: %f\n", kernels->spatialOrder, kernels->mode, score); + match->score = score; + match->spatialOrder = kernels->spatialOrder; + match->mode = kernels->mode; + match->nGood = nGood; + match->fluxes = psMemIncrRefCounter(fluxes); + match->chisq = psMemIncrRefCounter(chisqR); + match->moments = psMemIncrRefCounter(moments); + match->stampMask = psMemIncrRefCounter(stampMask); + } + } + + pmSubtractionVisualPlotChisqAndMoments(fluxes, chisqR, moments); + + psFree(stats); + psFree(chisqR); + psFree(chisqD); + psFree(fluxes); + psFree(moments); + psFree(stampMask); + + psFree(residual); + psFree(difference); + psFree(convolved1); + psFree(convolved2); + psFree(polyValues); + + psLogMsg("psModules.imcombine", PS_LOG_INFO, "Calculate Chisq and Moments: %f sec", psTimerClear("pmSubtractionCalculateChisqAndMoments")); + + return true; +} + +// XXX for now, let's not use this, and let's instead just use values from pmSubtractionCalculateChisqAndMoments +psVector *pmSubtractionCalculateDeviations(pmSubtractionStampList *stamps, pmSubtractionKernels *kernels) +{ + PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(stamps, NULL); + PM_ASSERT_SUBTRACTION_KERNELS_NON_NULL(kernels, NULL); + PM_ASSERT_SUBTRACTION_KERNELS_SOLUTION(kernels, NULL); + + psTimerStart("pmSubtractionCalculateDeviations"); psVector *deviations = psVectorAlloc(stamps->num, PS_TYPE_F32); // Mean deviation for stamps @@ -1298,7 +1727,4 @@ psImage *polyValues = NULL; // Polynomial values psKernel *residual = psKernelAlloc(-footprint, footprint, -footprint, footprint); // Residual image - - // set up holding images for the visualization - pmSubtractionVisualShowFitInit (stamps); psVector *fResSigma = psVectorAllocEmpty(stamps->num, PS_TYPE_F32); @@ -1401,7 +1827,4 @@ } - // XXX visualize the target, source, convolution and residual - pmSubtractionVisualShowFitAddStamp (target, source, residual, background, norm, i); - for (int y = - footprint; y <= footprint; y++) { for (int x = - footprint; x <= footprint; x++) { @@ -1438,7 +1861,4 @@ } - // XXX visualize the target, source, convolution and residual - pmSubtractionVisualShowFitAddStamp (image2, image1, residual, background, norm, i); - for (int y = - footprint; y <= footprint; y++) { for (int x = - footprint; x <= footprint; x++) { @@ -1455,4 +1875,5 @@ } + double flux = 0.0; double deviation = 0.0; // Sum of differences for (int y = - footprint; y <= footprint; y++) { @@ -1460,62 +1881,22 @@ double dev = PS_SQR(residual->kernel[y][x]) * weight->kernel[y][x]; deviation += dev; -#ifdef TESTING - residual->kernel[y][x] = dev; -#endif + flux += stamp->image1->kernel[y][x] + stamp->image2->kernel[y][x]; } } deviations->data.F32[i] = devNorm * deviation; - psTrace("psModules.imcombine", 5, "Deviation for stamp %d (%d,%d): %f\n", - i, (int)(stamp->x - 0.5), (int)(stamp->y - 0.5), deviations->data.F32[i]); + psTrace("psModules.imcombine", 5, "Deviation and Flux for stamp %d (%d,%d): %f %f\n", + i, (int)(stamp->x - 0.5), (int)(stamp->y - 0.5), deviations->data.F32[i], flux); psStringAppend(&log, "Stamp %d (%d,%d): %f\n", i, (int)(stamp->x - 0.5), (int)(stamp->y - 0.5), deviations->data.F32[i]); if (!isfinite(deviations->data.F32[i])) { stamp->status = PM_SUBTRACTION_STAMP_REJECTED; - psTrace("psModules.imcombine", 5, - "Rejecting stamp %d (%d,%d) because of non-finite deviation\n", - i, (int)(stamp->x - 0.5), (int)(stamp->y - 0.5)); + psTrace("psModules.imcombine", 5, "Rejecting stamp %d (%d,%d) because of non-finite deviation\n", i, (int)(stamp->x - 0.5), (int)(stamp->y - 0.5)); continue; } - -#ifdef TESTING - { - psString filename = NULL; - psStringAppend(&filename, "resid_%03d.fits", i); - psFits *fits = psFitsOpen(filename, "w"); - psFree(filename); - psFitsWriteImage(fits, NULL, residual->image, 0, NULL); - psFitsClose(fits); - } - if (stamp->image1) { - psString filename = NULL; - psStringAppend(&filename, "stamp_image1_%03d.fits", i); - psFits *fits = psFitsOpen(filename, "w"); - psFree(filename); - psFitsWriteImage(fits, NULL, stamp->image1->image, 0, NULL); - psFitsClose(fits); - } - if (stamp->image2) { - psString filename = NULL; - psStringAppend(&filename, "stamp_image2_%03d.fits", i); - psFits *fits = psFitsOpen(filename, "w"); - psFree(filename); - psFitsWriteImage(fits, NULL, stamp->image2->image, 0, NULL); - psFitsClose(fits); - } - if (stamp->weight) { - psString filename = NULL; - psStringAppend(&filename, "stamp_weight_%03d.fits", i); - psFits *fits = psFitsOpen(filename, "w"); - psFree(filename); - psFitsWriteImage(fits, NULL, stamp->weight->image, 0, NULL); - psFitsClose(fits); - } -#endif - } psFree(keepStamps); - psLogMsg("psModules.imcombine", PS_LOG_DETAIL, "%s", log); + psLogMsg("psModules.imcombine", PS_LOG_MINUTIA, "%s", log); psFree(log); @@ -1527,7 +1908,4 @@ psLogMsg("psModules.imcombine", PS_LOG_INFO, "normalization: %f, background: %f", norm, background); - pmSubtractionVisualShowFit(norm); - pmSubtractionVisualPlotFit(kernels); - psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV); psVectorStats (stats, fResSigma, NULL, NULL, 0); @@ -1560,303 +1938,7 @@ psFree(polyValues); + psLogMsg("psModules.imcombine", PS_LOG_INFO, "Calculate Deviations: %f sec", psTimerClear("pmSubtractionCalculateDeviations")); + return deviations; -} - -// we are supplied U, not Ut; w represents a diagonal matrix (also, we apply 1/w instead of w) -psImage *p_pmSubSolve_wUt (psVector *w, psImage *U) { - - psAssert (w->n == U->numCols, "w and U dimensions do not match"); - - // wUt has dimensions transposed relative to Ut. - psImage *wUt = psImageAlloc (U->numRows, U->numCols, PS_TYPE_F64); - psImageInit (wUt, 0.0); - - for (int i = 0; i < wUt->numCols; i++) { - for (int j = 0; j < wUt->numRows; j++) { - if (!isfinite(w->data.F64[j])) continue; - if (w->data.F64[j] == 0.0) continue; - wUt->data.F64[j][i] = U->data.F64[i][j] / w->data.F64[j]; - } - } - return wUt; -} - -// XXX this is just standard matrix multiplication: use psMatrixMultiply? -psImage *p_pmSubSolve_VwUt (psImage *V, psImage *wUt) { - - psAssert (V->numCols == wUt->numRows, "matrix dimensions do not match"); - - psImage *Ainv = psImageAlloc (wUt->numCols, V->numRows, PS_TYPE_F64); - - for (int i = 0; i < Ainv->numCols; i++) { - for (int j = 0; j < Ainv->numRows; j++) { - double sum = 0.0; - for (int k = 0; k < V->numCols; k++) { - sum += V->data.F64[j][k] * wUt->data.F64[k][i]; - } - Ainv->data.F64[j][i] = sum; - } - } - return Ainv; -} - -// we are supplied U, not Ut -bool p_pmSubSolve_UtB (psVector **UtB, psImage *U, psVector *B) { - - psAssert (U->numRows == B->n, "U and B dimensions do not match"); - - UtB[0] = psVectorRecycle (UtB[0], U->numCols, PS_TYPE_F64); - - for (int i = 0; i < U->numCols; i++) { - double sum = 0.0; - for (int j = 0; j < U->numRows; j++) { - sum += B->data.F64[j] * U->data.F64[j][i]; - } - UtB[0]->data.F64[i] = sum; - } - return true; -} - -// w is diagonal -bool p_pmSubSolve_wUtB (psVector **wUtB, psVector *w, psVector *UtB) { - - psAssert (w->n == UtB->n, "w and UtB dimensions do not match"); - - // wUt has dimensions transposed relative to Ut. - wUtB[0] = psVectorRecycle (wUtB[0], w->n, PS_TYPE_F64); - psVectorInit (wUtB[0], 0.0); - - for (int i = 0; i < w->n; i++) { - if (!isfinite(w->data.F64[i])) continue; - if (w->data.F64[i] == 0.0) continue; - wUtB[0]->data.F64[i] = UtB->data.F64[i] / w->data.F64[i]; - } - return true; -} - -// this is basically matrix * vector -bool p_pmSubSolve_VwUtB (psVector **VwUtB, psImage *V, psVector *wUtB) { - - psAssert (V->numCols == wUtB->n, "V and wUtB dimensions do not match"); - - VwUtB[0] = psVectorRecycle (*VwUtB, V->numRows, PS_TYPE_F64); - - for (int j = 0; j < V->numRows; j++) { - double sum = 0.0; - for (int i = 0; i < V->numCols; i++) { - sum += V->data.F64[j][i] * wUtB->data.F64[i]; - } - VwUtB[0]->data.F64[j] = sum; - } - return true; -} - -// this is basically matrix * vector -bool p_pmSubSolve_Ax (psVector **B, psImage *A, psVector *x) { - - psAssert (A->numCols == x->n, "A and x dimensions do not match"); - - B[0] = psVectorRecycle (*B, A->numRows, PS_TYPE_F64); - - for (int j = 0; j < A->numRows; j++) { - double sum = 0.0; - for (int i = 0; i < A->numCols; i++) { - sum += A->data.F64[j][i] * x->data.F64[i]; - } - B[0]->data.F64[j] = sum; - } - return true; -} - -// this is basically Vector * vector -bool p_pmSubSolve_VdV (double *value, psVector *x, psVector *y) { - - psAssert (x->n == y->n, "x and y dimensions do not match"); - - double sum = 0.0; - for (int i = 0; i < x->n; i++) { - sum += x->data.F64[i] * y->data.F64[i]; - } - *value = sum; - return true; -} - -bool p_pmSubSolve_y2 (double *y2, pmSubtractionKernels *kernels, const pmSubtractionStampList *stamps) { - - int footprint = stamps->footprint; // Half-size of stamps - - double sum = 0.0; - for (int i = 0; i < stamps->num; i++) { - - pmSubtractionStamp *stamp = stamps->stamps->data[i]; - if (stamp->status != PM_SUBTRACTION_STAMP_USED) continue; - - psKernel *weight = NULL; - psKernel *window = NULL; - psKernel *input = NULL; - -#ifdef USE_WEIGHT - weight = stamp->weight; -#endif -#ifdef USE_WINDOW - window = stamps->window; -#endif - - switch (kernels->mode) { - // MODE_1 : convolve image 1 to match image 2 (and vice versa) - case PM_SUBTRACTION_MODE_1: - input = stamp->image2; - break; - case PM_SUBTRACTION_MODE_2: - input = stamp->image1; - break; - default: - psAbort ("programming error"); - } - - for (int y = - footprint; y <= footprint; y++) { - for (int x = - footprint; x <= footprint; x++) { - double in = input->kernel[y][x]; - double value = in*in; - if (weight) { - float wtVal = weight->kernel[y][x]; - value *= wtVal; - } - if (window) { - float winVal = window->kernel[y][x]; - value *= winVal; - } - sum += value; - } - } - } - *y2 = sum; - return true; -} - -double p_pmSubSolve_ChiSquare (pmSubtractionKernels *kernels, const pmSubtractionStampList *stamps) { - - int footprint = stamps->footprint; // Half-size of stamps - int numKernels = kernels->num; // Number of kernels - - double sum = 0.0; - - psKernel *residual = psKernelAlloc(-footprint, footprint, -footprint, footprint); // Residual image - psImageInit(residual->image, 0.0); - - psImage *polyValues = NULL; // Polynomial values - - for (int i = 0; i < stamps->num; i++) { - - pmSubtractionStamp *stamp = stamps->stamps->data[i]; - if (stamp->status != PM_SUBTRACTION_STAMP_USED) continue; - - psKernel *weight = NULL; - psKernel *window = NULL; - psKernel *target = NULL; - psKernel *source = NULL; - - psArray *convolutions = NULL; - -#ifdef USE_WEIGHT - weight = stamp->weight; -#endif -#ifdef USE_WINDOW - window = stamps->window; -#endif - - switch (kernels->mode) { - // MODE_1 : convolve image 1 to match image 2 (and vice versa) - case PM_SUBTRACTION_MODE_1: - target = stamp->image2; - source = stamp->image1; - convolutions = stamp->convolutions1; - break; - case PM_SUBTRACTION_MODE_2: - target = stamp->image1; - source = stamp->image2; - convolutions = stamp->convolutions2; - break; - default: - psAbort ("programming error"); - } - - // Calculate coefficients of the kernel basis functions - polyValues = p_pmSubtractionPolynomial(polyValues, kernels->spatialOrder, stamp->xNorm, stamp->yNorm); - double norm = p_pmSubtractionSolutionNorm(kernels); // Normalisation - double background = p_pmSubtractionSolutionBackground(kernels, polyValues);// Difference in background - - psImageInit(residual->image, 0.0); - for (int j = 0; j < numKernels; j++) { - psKernel *convolution = convolutions->data[j]; // Convolution - double coefficient = p_pmSubtractionSolutionCoeff(kernels, polyValues, j, false); // Coefficient - for (int y = - footprint; y <= footprint; y++) { - for (int x = - footprint; x <= footprint; x++) { - residual->kernel[y][x] -= convolution->kernel[y][x] * coefficient; - } - } - } - - for (int y = - footprint; y <= footprint; y++) { - for (int x = - footprint; x <= footprint; x++) { - double resid = target->kernel[y][x] - background - source->kernel[y][x] * norm + residual->kernel[y][x]; - double value = PS_SQR(resid); - if (weight) { - float wtVal = weight->kernel[y][x]; - value *= wtVal; - } - if (window) { - float winVal = window->kernel[y][x]; - value *= winVal; - } - sum += value; - } - } - } - psFree (polyValues); - psFree (residual); - - return sum; -} - -bool p_pmSubSolve_SetWeights (psVector *wApply, psVector *w, psVector *wMask) { - - for (int i = 0; i < w->n; i++) { - wApply->data.F64[i] = wMask->data.U8[i] ? 0.0 : w->data.F64[i]; - } - return true; -} - -// we are supplied V and w; w represents a diagonal matrix (also, we apply 1/w instead of w) -psImage *p_pmSubSolve_Xvar (psImage *V, psVector *w) { - - psAssert (w->n == V->numCols, "w and U dimensions do not match"); - - psImage *Vn = psImageAlloc (V->numCols, V->numRows, PS_TYPE_F64); - psImageInit (Vn, 0.0); - - // generate Vn = V * w^{-1} - for (int j = 0; j < Vn->numRows; j++) { - for (int i = 0; i < Vn->numCols; i++) { - if (!isfinite(w->data.F64[i])) continue; - if (w->data.F64[i] == 0.0) continue; - Vn->data.F64[j][i] = V->data.F64[j][i] / w->data.F64[i]; - } - } - - psImage *Xvar = psImageAlloc (V->numCols, V->numRows, PS_TYPE_F64); - psImageInit (Xvar, 0.0); - - // generate Xvar = Vn * Vn^T - for (int j = 0; j < Vn->numRows; j++) { - for (int i = 0; i < Vn->numCols; i++) { - double sum = 0.0; - for (int k = 0; k < Vn->numCols; k++) { - sum += Vn->data.F64[k][i]*Vn->data.F64[k][j]; - } - Xvar->data.F64[j][i] = sum; - } - } - return Xvar; } @@ -1864,5 +1946,4 @@ // of the elements of an image A(x,y) = A->data.F64[y][x] = A_x,y, a matrix // multiplication is: A_k,j * B_i,k = C_i,j - bool psFitsWriteImageSimple (char *filename, psImage *image, psMetadata *header) { Index: trunk/psModules/src/imcombine/pmSubtractionEquation.h =================================================================== --- trunk/psModules/src/imcombine/pmSubtractionEquation.h (revision 30621) +++ trunk/psModules/src/imcombine/pmSubtractionEquation.h (revision 30622) @@ -4,12 +4,5 @@ #include "pmSubtractionStamps.h" #include "pmSubtractionKernels.h" - -typedef enum { - PM_SUBTRACTION_EQUATION_NONE = 0x00, - PM_SUBTRACTION_EQUATION_NORM = 0x01, - PM_SUBTRACTION_EQUATION_BG = 0x02, - PM_SUBTRACTION_EQUATION_KERNELS = 0x04, - PM_SUBTRACTION_EQUATION_ALL = 0x07, // value should be NORM | BG | KERNELS -} pmSubtractionEquationCalculationMode; +#include "pmSubtraction.h" /// Execute a thread job to calculate the least-squares equation for a stamp @@ -20,18 +13,15 @@ bool pmSubtractionCalculateEquationStamp(pmSubtractionStampList *stamps, ///< Stamps pmSubtractionKernels *kernels, ///< Kernel parameters - int index, ///< Index of stamp - const pmSubtractionEquationCalculationMode mode + int index ///< Index of stamp ); /// Calculate the least-squares equation to match the image quality bool pmSubtractionCalculateEquation(pmSubtractionStampList *stamps, ///< Stamps - pmSubtractionKernels *kernels, ///< Kernel parameters - const pmSubtractionEquationCalculationMode mode + pmSubtractionKernels *kernels ///< Kernel parameters ); /// Solve the least-squares equation to match the image quality bool pmSubtractionSolveEquation(pmSubtractionKernels *kernels, ///< Kernel parameters - const pmSubtractionStampList *stamps, ///< Stamps - const pmSubtractionEquationCalculationMode mode + const pmSubtractionStampList *stamps ///< Stamps ); @@ -92,3 +82,6 @@ bool pmSubtractionCalculateMomentsKernel(double *Mxx, double *Myy, psKernel *image, int footprint, int window); +bool pmSubtractionChisqStats(psVector *fluxesVector, psVector *chisqDVector, psVector *chisqRVector, psVector *momentVector, psVector *stampMask, psKernel *convolved1, psKernel *convolved2, psKernel *difference, psKernel *residual, psKernel *weight, psKernel *window); + +bool pmSubtractionCalculateChisqAndMoments(pmSubtractionQuality **bestMatch, pmSubtractionStampList *stamps, pmSubtractionKernels *kernels); #endif Index: trunk/psModules/src/imcombine/pmSubtractionEquation.v0.c =================================================================== --- trunk/psModules/src/imcombine/pmSubtractionEquation.v0.c (revision 30621) +++ trunk/psModules/src/imcombine/pmSubtractionEquation.v0.c (revision 30622) @@ -9,4 +9,5 @@ #include "pmErrorCodes.h" #include "pmSubtraction.h" +#include "pmSubtractionTypes.h" #include "pmSubtractionKernels.h" #include "pmSubtractionStamps.h" Index: trunk/psModules/src/imcombine/pmSubtractionHermitian.c =================================================================== --- trunk/psModules/src/imcombine/pmSubtractionHermitian.c (revision 30621) +++ trunk/psModules/src/imcombine/pmSubtractionHermitian.c (revision 30622) @@ -8,4 +8,5 @@ #include +#include "pmSubtractionTypes.h" #include "pmSubtractionHermitian.h" Index: trunk/psModules/src/imcombine/pmSubtractionIO.c =================================================================== --- trunk/psModules/src/imcombine/pmSubtractionIO.c (revision 30621) +++ trunk/psModules/src/imcombine/pmSubtractionIO.c (revision 30622) @@ -12,4 +12,5 @@ #include "pmConceptsRead.h" +#include "pmSubtractionTypes.h" #include "pmSubtraction.h" #include "pmSubtractionKernels.h" Index: trunk/psModules/src/imcombine/pmSubtractionKernels.c =================================================================== --- trunk/psModules/src/imcombine/pmSubtractionKernels.c (revision 30621) +++ trunk/psModules/src/imcombine/pmSubtractionKernels.c (revision 30622) @@ -8,4 +8,6 @@ #include +#include "pmFPA.h" +#include "pmSubtractionTypes.h" #include "pmSubtraction.h" #include "pmSubtractionKernels.h" @@ -20,4 +22,6 @@ { psFree(kernels->description); + psFree(kernels->fwhms); + psFree(kernels->orders); psFree(kernels->u); psFree(kernels->v); @@ -30,4 +34,6 @@ psFree(kernels->solution1); psFree(kernels->solution2); + psFree(kernels->solution1err); + psFree(kernels->solution2err); psFree(kernels->sampleStamps); } @@ -419,18 +425,12 @@ int num = 0; // Number of basis functions - psString params = NULL; // List of parameters for (int i = 0; i < numGaussians; i++) { int gaussOrder = orders->data.S32[i]; // Polynomial order to apply to Gaussian - psStringAppend(¶ms, "(%.1f,%d)", fwhms->data.F32[i], orders->data.S32[i]); num += (gaussOrder + 1) * (gaussOrder + 2) / 2; } - pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_ISIS, size, - spatialOrder, penalty, bounds, mode); // Kernels - psStringAppend(&kernels->description, "ISIS(%d,%s,%d,%.2e)", size, params, spatialOrder, penalty); - - psLogMsg("psModules.imcombine", PS_LOG_INFO, "ISIS kernel: %s,%d --> %d elements", - params, spatialOrder, num); - psFree(params); + pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_ISIS, size, fwhms, orders, spatialOrder, penalty, bounds, mode); // Kernels + pmSubtractionKernelsMakeDescription(kernels); + psLogMsg("psModules.imcombine", PS_LOG_INFO, "kernel: %s --> %d elements", kernels->description, num); # if (!CENTRAL_DELTA && !ZERO_KERNEL_ZERO_FLUX) @@ -503,18 +503,13 @@ int num = 0; // Number of basis functions - psString params = NULL; // List of parameters for (int i = 0; i < numGaussians; i++) { int gaussOrder = orders->data.S32[i]; // Polynomial order to apply to Gaussian - psStringAppend(¶ms, "(%.1f,%d)", fwhms->data.F32[i], orders->data.S32[i]); num += (gaussOrder + 1) * (gaussOrder + 2) / 2; num += (11 - gaussOrder - 1); // include all higher order radial terms } - pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_ISIS_RADIAL, size, - spatialOrder, penalty, bounds, mode); // Kernels - psStringAppend(&kernels->description, "ISIS_RADIAL(%d,%s,%d,%.2e)", size, params, spatialOrder, penalty); - - psLogMsg("psModules.imcombine", PS_LOG_INFO, "ISIS_RADIAL kernel: %s,%d --> %d elements", params, spatialOrder, num); - psFree(params); + pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_ISIS_RADIAL, size, fwhms, orders, spatialOrder, penalty, bounds, mode); // Kernels + pmSubtractionKernelsMakeDescription(kernels); + psLogMsg("psModules.imcombine", PS_LOG_INFO, "kernel: %s --> %d elements", kernels->description, num); // Set the kernel parameters @@ -569,18 +564,12 @@ int num = 0; // Number of basis functions - psString params = NULL; // List of parameters for (int i = 0; i < numGaussians; i++) { int gaussOrder = orders->data.S32[i]; // Polynomial order to apply to Gaussian - psStringAppend(¶ms, "(%.1f,%d)", fwhms->data.F32[i], orders->data.S32[i]); num += (gaussOrder + 1) * (gaussOrder + 2) / 2; } - pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_HERM, size, - spatialOrder, penalty, bounds, mode); // Kernels - psStringAppend(&kernels->description, "HERM(%d,%s,%d,%.2e)", size, params, spatialOrder, penalty); - - psLogMsg("psModules.imcombine", PS_LOG_INFO, "HERM kernel: %s,%d --> %d elements", - params, spatialOrder, num); - psFree(params); + pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_HERM, size, fwhms, orders, spatialOrder, penalty, bounds, mode); // Kernels + pmSubtractionKernelsMakeDescription(kernels); + psLogMsg("psModules.imcombine", PS_LOG_INFO, "kernel: %s --> %d elements", kernels->description, num); // Set the kernel parameters @@ -627,17 +616,12 @@ int num = 0; // Number of basis functions - psString params = NULL; // List of parameters for (int i = 0; i < numGaussians; i++) { int gaussOrder = orders->data.S32[i]; // Polynomial order to apply to Gaussian - psStringAppend(¶ms, "(%.1f,%d)", fwhms->data.F32[i], orders->data.S32[i]); num += PS_SQR(gaussOrder + 1); } - pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_DECONV_HERM, size, - spatialOrder, penalty, bounds, mode); // Kernels - psStringAppend(&kernels->description, "DECONV_HERM(%d,%s,%d,%.2e)", size, params, spatialOrder, penalty); - - psLogMsg("psModules.imcombine", PS_LOG_INFO, "DECONVOLVED HERM kernel: %s,%d --> %d elements", params, spatialOrder, num); - psFree(params); + pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_DECONV_HERM, size, fwhms, orders, spatialOrder, penalty, bounds, mode); // Kernels + pmSubtractionKernelsMakeDescription(kernels); + psLogMsg("psModules.imcombine", PS_LOG_INFO, "kernel: %s --> %d elements", kernels->description, num); // XXXXX hard-wired reference sigma for now of 1.7 pix (== 4.0 pix fwhm == 1.0 arcsec in simtest) @@ -713,5 +697,5 @@ pmSubtractionKernels *pmSubtractionKernelsAlloc(int numBasisFunctions, pmSubtractionKernelsType type, - int size, int spatialOrder, float penalty, psRegion bounds, + int size, psVector *fwhms, psVector *orders, int spatialOrder, float penalty, psRegion bounds, pmSubtractionMode mode) { @@ -722,4 +706,6 @@ kernels->description = NULL; kernels->num = numBasisFunctions; + kernels->fwhms = psMemIncrRefCounter(fwhms); + kernels->orders = psMemIncrRefCounter(orders); kernels->u = psVectorAlloc(numBasisFunctions, PS_TYPE_S32); kernels->v = psVectorAlloc(numBasisFunctions, PS_TYPE_S32); @@ -740,4 +726,6 @@ kernels->size = size; kernels->inner = 0; + kernels->binning = 0; + kernels->ringsOrder = 0; kernels->spatialOrder = spatialOrder; kernels->bgOrder = 0; @@ -745,4 +733,6 @@ kernels->solution1 = NULL; kernels->solution2 = NULL; + kernels->solution1err = NULL; + kernels->solution2err = NULL; kernels->mean = NAN; kernels->rms = NAN; @@ -823,9 +813,7 @@ int num = PS_SQR(2 * size + 1) - 1; // Number of basis functions - pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(0, PM_SUBTRACTION_KERNEL_POIS, size, - spatialOrder, penalty, bounds, mode); // Kernels - psStringAppend(&kernels->description, "POIS(%d,%d,%.2e)", size, spatialOrder, penalty); - psLogMsg("psModules.imcombine", PS_LOG_INFO, "POIS kernel: %d,%d --> %d elements", - size, spatialOrder, num); + pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(0, PM_SUBTRACTION_KERNEL_POIS, size, NULL, NULL, spatialOrder, penalty, bounds, mode); // Kernels + pmSubtractionKernelsMakeDescription(kernels); + psLogMsg("psModules.imcombine", PS_LOG_INFO, "kernel: %s --> %d elements", kernels->description, num); if (!p_pmSubtractionKernelsAddGrid(kernels, 0, size)) { @@ -871,12 +859,9 @@ psTrace("psModules.imcombine", 3, "Number of basis functions: %d\n", num); - pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_SPAM, size, - spatialOrder, penalty, bounds, mode); // Kernels + pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_SPAM, size, NULL, NULL, spatialOrder, penalty, bounds, mode); // Kernels kernels->inner = inner; - psStringAppend(&kernels->description, "SPAM(%d,%d,%d,%d,%.2e)", size, inner, binning, spatialOrder, - penalty); - - psLogMsg("psModules.imcombine", PS_LOG_INFO, "SPAM kernel: %d,%d,%d,%d --> %d elements", - size, inner, binning, spatialOrder, num); + kernels->binning = binning; + pmSubtractionKernelsMakeDescription(kernels); + psLogMsg("psModules.imcombine", PS_LOG_INFO, "kernel: %s --> %d elements", kernels->description, num); kernels->uStop = psVectorAlloc(num, PS_TYPE_S32); @@ -970,11 +955,8 @@ psTrace("psModules.imcombine", 3, "Number of basis functions: %d\n", num); - pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_FRIES, size, - spatialOrder, penalty, bounds, mode); // Kernels + pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_FRIES, size, NULL, NULL, spatialOrder, penalty, bounds, mode); // Kernels kernels->inner = inner; - psStringAppend(&kernels->description, "FRIES(%d,%d,%d,%.2e)", size, inner, spatialOrder, penalty); - - psLogMsg("psModules.imcombine", PS_LOG_INFO, "FRIES kernel: %d,%d,%d --> %d elements", - size, inner, spatialOrder, num); + pmSubtractionKernelsMakeDescription(kernels); + psLogMsg("psModules.imcombine", PS_LOG_INFO, "kernel: %s --> %d elements", kernels->description, num); kernels->uStop = psVectorAlloc(num, PS_TYPE_S32); @@ -1053,9 +1035,9 @@ PS_ASSERT_INT_LESS_THAN(inner, size, NULL); - pmSubtractionKernels *kernels = p_pmSubtractionKernelsRawISIS(size, spatialOrder, fwhms, orders, - penalty, bounds, mode); // Kernels + pmSubtractionKernels *kernels = p_pmSubtractionKernelsRawISIS(size, spatialOrder, fwhms, orders, penalty, bounds, mode); // Kernels kernels->type = PM_SUBTRACTION_KERNEL_GUNK; - psStringPrepend(&kernels->description, "GUNK="); - psStringAppend(&kernels->description, "+POIS(%d,%d)", inner, spatialOrder); + kernels->inner = inner; + pmSubtractionKernelsMakeDescription(kernels); + psLogMsg("psModules.imcombine", PS_LOG_INFO, "kernel: %s --> %d elements", kernels->description, (int) kernels->num); int numISIS = kernels->num; // Number of ISIS kernels @@ -1100,12 +1082,9 @@ int num = numRings * numPoly; // Total number of basis functions - pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_RINGS, size, - spatialOrder, penalty, bounds, mode); // Kernels + pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_RINGS, size, NULL, NULL, spatialOrder, penalty, bounds, mode); // Kernels kernels->inner = inner; - psStringAppend(&kernels->description, "RINGS(%d,%d,%d,%d,%.2e)", size, inner, ringsOrder, spatialOrder, - penalty); - - psLogMsg("psModules.imcombine", PS_LOG_INFO, "RINGS kernel: %d,%d,%d,%d --> %d elements", - size, inner, ringsOrder, spatialOrder, num); + kernels->ringsOrder = ringsOrder; + pmSubtractionKernelsMakeDescription(kernels); + psLogMsg("psModules.imcombine", PS_LOG_INFO, "kernel: %s --> %d elements", kernels->description, num); // Set the Gaussian kernel parameters @@ -1405,4 +1384,63 @@ } +bool pmSubtractionKernelsMakeDescription(pmSubtractionKernels *kernels) { + + // free if it exists + psFree (kernels->description); + + // generate the description parameter string + psString params = NULL; + if (kernels->fwhms) { + for (int i = 0; i < kernels->fwhms->n; i++) { + psStringAppend(¶ms, "(%.1f,%d)", kernels->fwhms->data.F32[i], kernels->orders->data.S32[i]); + } + } + + switch (kernels->type) { + case PM_SUBTRACTION_KERNEL_ISIS: + psStringAppend (&kernels->description, "ISIS(%d,%s,%d,%.2e)", kernels->size, params, kernels->spatialOrder, kernels->penalty); + break; + + case PM_SUBTRACTION_KERNEL_ISIS_RADIAL: + psStringAppend(&kernels->description, "ISIS_RADIAL(%d,%s,%d,%.2e)", kernels->size, params, kernels->spatialOrder, kernels->penalty); + break; + + case PM_SUBTRACTION_KERNEL_HERM: + psStringAppend(&kernels->description, "HERM(%d,%s,%d,%.2e)", kernels->size, params, kernels->spatialOrder, kernels->penalty); + break; + + case PM_SUBTRACTION_KERNEL_DECONV_HERM: + psStringAppend(&kernels->description, "DECONV_HERM(%d,%s,%d,%.2e)", kernels->size, params, kernels->spatialOrder, kernels->penalty); + break; + + case PM_SUBTRACTION_KERNEL_POIS: + psStringAppend(&kernels->description, "POIS(%d,%d,%.2e)", kernels->size, kernels->spatialOrder, kernels->penalty); + break; + + case PM_SUBTRACTION_KERNEL_SPAM: + psStringAppend(&kernels->description, "SPAM(%d,%d,%d,%d,%.2e)", kernels->size, kernels->inner, kernels->binning, kernels->spatialOrder, kernels->penalty); + break; + + case PM_SUBTRACTION_KERNEL_FRIES: + psStringAppend(&kernels->description, "FRIES(%d,%d,%d,%.2e)", kernels->size, kernels->inner, kernels->spatialOrder, kernels->penalty); + break; + + // Grid United with Normal Kernel [description: GUNK=ISIS(...)+POIS(...)] + case PM_SUBTRACTION_KERNEL_GUNK: + psStringAppend(&kernels->description, "GUNK=ISIS(%d,%s,%d,%.2e)", kernels->size, params, kernels->spatialOrder, kernels->penalty); + psStringAppend(&kernels->description, "+POIS(%d,%d)", kernels->inner, kernels->spatialOrder); + break; + + case PM_SUBTRACTION_KERNEL_RINGS: + psStringAppend(&kernels->description, "RINGS(%d,%d,%d,%d,%.2e)", kernels->size, kernels->inner, kernels->ringsOrder, kernels->spatialOrder, kernels->penalty); + break; + + default: + psAbort("unknown kernel"); + } + psFree (params); + return true; +} + pmSubtractionKernels *pmSubtractionKernelsCopy(const pmSubtractionKernels *in) { @@ -1435,4 +1473,6 @@ out->solution1 = in->solution1 ? psVectorCopy(NULL, in->solution1, PS_TYPE_F64) : NULL; out->solution2 = in->solution2 ? psVectorCopy(NULL, in->solution2, PS_TYPE_F64) : NULL; + out->solution1err = in->solution1err ? psVectorCopy(NULL, in->solution1err, PS_TYPE_F64) : NULL; + out->solution2err = in->solution2err ? psVectorCopy(NULL, in->solution2err, PS_TYPE_F64) : NULL; out->sampleStamps = psMemIncrRefCounter(in->sampleStamps); Index: trunk/psModules/src/imcombine/pmSubtractionKernels.h =================================================================== --- trunk/psModules/src/imcombine/pmSubtractionKernels.h (revision 30621) +++ trunk/psModules/src/imcombine/pmSubtractionKernels.h (revision 30622) @@ -2,72 +2,6 @@ #define PM_SUBTRACTION_KERNELS_H -#include -#include - -/// Type of subtraction kernel -typedef enum { - PM_SUBTRACTION_KERNEL_NONE, ///< Nothing --- an error - PM_SUBTRACTION_KERNEL_POIS, ///< Pan-STARRS Optimal Image Subtraction --- delta functions - PM_SUBTRACTION_KERNEL_ISIS, ///< Traditional kernel --- gaussians modified by polynomials - PM_SUBTRACTION_KERNEL_ISIS_RADIAL, ///< ISIS + higher-order radial Hermitians - PM_SUBTRACTION_KERNEL_HERM, ///< Hermitian polynomial kernels - PM_SUBTRACTION_KERNEL_DECONV_HERM, ///< Deconvolved Hermitian polynomial kernels - PM_SUBTRACTION_KERNEL_SPAM, ///< Summed Pixels for Advanced Matching --- summed delta functions - PM_SUBTRACTION_KERNEL_FRIES, ///< Fibonacci Radius Increases Excellence of Subtraction - PM_SUBTRACTION_KERNEL_GUNK, ///< Grid United with Normal Kernel --- POIS and ISIS hybrid - PM_SUBTRACTION_KERNEL_RINGS, ///< Rings Instead of the Normal Gaussian Subtraction -} pmSubtractionKernelsType; - -/// Modes --- specifies which image to convolve -typedef enum { - PM_SUBTRACTION_MODE_ERR, // Error in the mode - PM_SUBTRACTION_MODE_1, // Convolve image 1 - PM_SUBTRACTION_MODE_2, // Convolve image 2 - PM_SUBTRACTION_MODE_UNSURE, // Not sure yet which image to convolve so try to satisfy both - PM_SUBTRACTION_MODE_DUAL, // Dual convolution -} pmSubtractionMode; - -/// Kernels specification -typedef struct { - pmSubtractionKernelsType type; ///< Type of kernels --- allowing the use of multiple kernels - psString description; ///< Description of the kernel parameters - int xMin, xMax, yMin, yMax; ///< Bounds of image (for normalisation) - long num; ///< Number of kernel components (not including the spatial ones) - psVector *u, *v; ///< Offset (for POIS) or polynomial order (for ISIS, HERM or DECONV_HERM) - psVector *widths; ///< Gaussian FWHMs (ISIS, HERM or DECONV_HERM) - psVector *uStop, *vStop; ///< Width of kernel element (SPAM,FRIES only) - psArray *preCalc; ///< Array of images containing pre-calculated kernel (for ISIS, HERM or DECONV_HERM) - float penalty; ///< Penalty for wideness - psVector *penalties1; ///< Penalty for each kernel component - psVector *penalties2; ///< Penalty for each kernel component - bool havePenalties; ///< flag to test if we have already calculated the penalties or not. - int size; ///< The half-size of the kernel - int inner; ///< The size of an inner region - int spatialOrder; ///< The spatial order of the kernels - int bgOrder; ///< The order for the background fitting - pmSubtractionMode mode; ///< Mode for subtraction - psVector *solution1, *solution2; ///< Solution for the PSF matching - // Quality information - float mean, rms; ///< Mean and RMS of chi^2 from stamps - int numStamps; ///< Number of good stamps - float fResSigmaMean; ///< mean fractional stdev of residuals - float fResSigmaStdev; ///< stdev of fractional stdev of residuals - float fResOuterMean; ///< mean fractional positive swing in residuals - float fResOuterStdev; ///< stdev of fractional positive swing in residuals - float fResTotalMean; ///< mean fractional negative swing in residuals - float fResTotalStdev; ///< stdev of fractional negative swing in residuals - psArray *sampleStamps; ///< array of brightest set of stamps for output visualizations -} pmSubtractionKernels; - -// pmSubtractionKernels->preCalc is an array of pmSubtractionKernelPreCalc structures -typedef struct { - psVector *uCoords; // used by RINGS - psVector *vCoords; // used by RINGS - psVector *poly; // used by RINGS - - psVector *xKernel; // used by ISIS, HERM, DECONV_HERM - psVector *yKernel; // used by ISIS, HERM, DECONV_HERM - psKernel *kernel; // used by ISIS, HERM, DECONV_HERM -} pmSubtractionKernelPreCalc; +// #include +// #include // Assertion to check pmSubtractionKernels @@ -162,4 +96,6 @@ pmSubtractionKernelsType type, ///< Kernel type int size, ///< Half-size of kernel + psVector *fwhms, ///< requested kernel basis function + psVector *orders, int spatialOrder, ///< Order of spatial variations float penalty, ///< Penalty for wideness @@ -303,4 +239,7 @@ ); +bool pmSubtractionKernelsMakeDescription(pmSubtractionKernels *kernels); + + /// Copy kernels /// Index: trunk/psModules/src/imcombine/pmSubtractionMask.c =================================================================== --- trunk/psModules/src/imcombine/pmSubtractionMask.c (revision 30621) +++ trunk/psModules/src/imcombine/pmSubtractionMask.c (revision 30622) @@ -7,4 +7,6 @@ #include "pmErrorCodes.h" +#include "pmFPA.h" +#include "pmSubtractionTypes.h" #include "pmSubtraction.h" #include "pmSubtractionKernels.h" Index: trunk/psModules/src/imcombine/pmSubtractionMatch.c =================================================================== --- trunk/psModules/src/imcombine/pmSubtractionMatch.c (revision 30621) +++ trunk/psModules/src/imcombine/pmSubtractionMatch.c (revision 30622) @@ -11,9 +11,10 @@ #include "pmFPA.h" #include "pmHDUUtils.h" +#include "pmSubtractionTypes.h" +#include "pmSubtraction.h" #include "pmSubtractionParams.h" #include "pmSubtractionKernels.h" #include "pmSubtractionStamps.h" #include "pmSubtractionEquation.h" -#include "pmSubtraction.h" #include "pmSubtractionAnalysis.h" #include "pmSubtractionMask.h" @@ -27,6 +28,4 @@ static bool useFFT = true; // Do convolutions using FFT - -# define SUBMODE PM_SUBTRACTION_EQUATION_ALL //#define TESTING @@ -59,53 +58,4 @@ } - -static bool subtractionGetStamps(pmSubtractionStampList **stamps, // Stamps to read - const pmReadout *ro1, // Readout 1 - const pmReadout *ro2, // Readout 2 - const psImage *subMask, // Mask for subtraction, or NULL - psImage *variance, // Variance map - const psRegion *region, // Region of interest - float thresh1, // Threshold for stamp finding on readout 1 - float thresh2, // Threshold for stamp finding on readout 2 - float stampSpacing, // Spacing between stamps - float normFrac, // Fraction of flux in window for normalisation window - float sysError, // Relative systematic error in images - float skyError, // Relative systematic error in images - int size, // Kernel half-size - int footprint, // Convolution footprint for stamps - pmSubtractionMode mode // Mode for subtraction - ) -{ - PS_ASSERT_PTR_NON_NULL(stamps, false); - PM_ASSERT_READOUT_NON_NULL(ro1, false); - PM_ASSERT_READOUT_NON_NULL(ro2, false); - PS_ASSERT_IMAGE_NON_NULL(subMask, false); - PS_ASSERT_IMAGE_NON_NULL(variance, false); - PS_ASSERT_PTR_NON_NULL(region, false); - - psTrace("psModules.imcombine", 3, "Finding stamps...\n"); - - psImage *image1 = ro1 ? ro1->image : NULL, *image2 = ro2 ? ro2->image : NULL; // Images of interest - - *stamps = pmSubtractionStampsFind(*stamps, image1, image2, subMask, region, thresh1, thresh2, - size, footprint, stampSpacing, normFrac, sysError, skyError, mode); - if (!*stamps) { - psError(psErrorCodeLast(), false, "Unable to find stamps."); - return false; - } - - memCheck(" find stamps"); - - psTrace("psModules.imcombine", 3, "Extracting stamps...\n"); - if (!pmSubtractionStampsExtract(*stamps, ro1->image, ro2->image, variance, size, *region)) { - psError(psErrorCodeLast(), false, "Unable to extract stamps."); - return false; - } - - memCheck(" extract stamps"); - pmSubtractionVisualPlotStamps(*stamps, (pmReadout *) ro1); - return true; -} - // Check input arguments static bool subtractionMatchCheck(pmReadout *conv1, pmReadout *conv2, // Convolved images @@ -123,5 +73,5 @@ float badFrac, // Maximum fraction of bad input pixels to accept pmSubtractionMode subMode // Mode of subtraction - ) + ) { if (subMode != PM_SUBTRACTION_MODE_2) { @@ -473,4 +423,55 @@ } +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; + } else { + pmSubtractionKernelsMakeDescription(kernels); + psLogMsg("psModules.imcombine", PS_LOG_INFO, "final kernel: %s", kernels->description); + } + return true; +} bool pmSubtractionMatch(pmReadout *conv1, pmReadout *conv2, const pmReadout *ro1, const pmReadout *ro2, @@ -478,5 +479,5 @@ const psArray *sources, const char *stampsName, pmSubtractionKernelsType type, int size, int spatialOrder, - const psVector *isisWidths, const psVector *isisOrders, + psVector *isisWidths, const psVector *isisOrders, int inner, int ringsOrder, int binning, float penalty, bool optimum, const psVector *optFWHMs, int optOrder, float optThreshold, @@ -559,4 +560,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"); @@ -628,5 +655,5 @@ regionString = psRegionToString(*region); psLogMsg("psModules.imcombine", PS_LOG_DETAIL, "Iso-kernel region: %s out of %d,%d\n", - regionString, numCols, numRows); + regionString, numCols, numRows); if (stampsName && strlen(stampsName) > 0) { @@ -640,167 +667,182 @@ } - // We get the stamps here; we will also attempt to get stamps at the first iteration, but it - // doesn't matter. - if (!subtractionGetStamps(&stamps, ro1, ro2, subMask, variance, region, stampThresh1, stampThresh2, - stampSpacing, normFrac, sysError, skyError, size, footprint, subMode)) { - goto MATCH_ERROR; - } - - - // generate the window function from the set of stamps - if (!pmSubtractionStampsGetWindow(stamps, size)) { - psError(psErrorCodeLast(), false, "Unable to get stamp window."); - goto MATCH_ERROR; - } - - // Define kernel basis functions - if (optimum && (type == PM_SUBTRACTION_KERNEL_ISIS || type == PM_SUBTRACTION_KERNEL_GUNK)) { - kernels = pmSubtractionKernelsOptimumISIS(type, size, inner, spatialOrder, - optFWHMs, optOrder, stamps, footprint, - optThreshold, penalty, bounds, subMode); - if (!kernels) { - psErrorClear(); - psWarning("Unable to derive optimum ISIS kernel --- switching to default."); - } - } - if (kernels == NULL) { - // Not an ISIS/GUNK kernel, or the optimum kernel search failed - kernels = pmSubtractionKernelsGenerate(type, size, spatialOrder, isisWidths, isisOrders, - inner, binning, ringsOrder, penalty, bounds, subMode); - } - - memCheck("kernels"); - - if (subMode == PM_SUBTRACTION_MODE_UNSURE) { - pmSubtractionMode newMode = pmSubtractionBestMode(&stamps, &kernels, subMask, rej); - switch (newMode) { - case PM_SUBTRACTION_MODE_1: - psLogMsg("psModules.imcombine", PS_LOG_INFO, "Convolving image 1 to match image 2."); - break; - case PM_SUBTRACTION_MODE_2: - psLogMsg("psModules.imcombine", PS_LOG_INFO, "Convolving image 2 to match image 1."); - break; - default: - psError(psErrorCodeLast(), false, "Unable to determine subtraction order."); - goto MATCH_ERROR; - } - subMode = newMode; - } - - int numRejected = -1; // Number of rejected stamps in each iteration - for (int k = 0; k < iter && numRejected != 0; k++) { - psLogMsg("psModules.imcombine", PS_LOG_INFO, "Iteration %d.", k); - - if (!subtractionGetStamps(&stamps, ro1, ro2, subMask, variance, region, - stampThresh1, stampThresh2, stampSpacing, normFrac, - sysError, skyError, size, footprint, subMode)) { - goto MATCH_ERROR; - } - - // generate the window function from the set of stamps - if (!pmSubtractionStampsGetWindow(stamps, size)) { - psError(psErrorCodeLast(), false, "Unable to get stamps window."); - goto MATCH_ERROR; - } + bool tryAgain = true; + while (tryAgain) { + // We get the stamps here; we will also attempt to get stamps at the first iteration, but it + // doesn't matter. + if (!pmSubtractionStampsSelect(&stamps, ro1, ro2, subMask, variance, region, stampThresh1, stampThresh2, + stampSpacing, normFrac, sysError, skyError, size, footprint, subMode)) { + goto MATCH_ERROR; + } + + // generate the window function from the set of stamps + if (!pmSubtractionStampsGetWindow(&tryAgain, stamps, size)) { + // if we failed, it might be due to the desired normWindow being larger than the current footprint. + // in this case, just adjust the footprint and try again. + if (tryAgain) { + // keep the border constant + int border = footprint - size; + size = PS_MAX(stamps->normWindow1, stamps->normWindow2) + 2; + footprint = size + border; + + // we need to reconstruct everything, so just free the stamps here and retry + psFree(stamps); + } else { + // unrecoverable error + psError(psErrorCodeLast(), false, "Unable to get stamp window."); + goto MATCH_ERROR; + } + } + } + + // check on the kernel scaling -- if the kron-based radial moments are very different, adjust to match them + { + // float fwhm1; + // float fwhm2; + // pmSubtractionGetFWHMs(&fwhm1, &fwhm2); + // psAssert(isfinite(fwhm1), "fwhm 1 not set"); + // psAssert(isfinite(fwhm2), "fwhm 2 not set"); + + // XXX this is BAD: depends on the relationship below: + // stamps->normWindow1 = 2.75*R1; + // stamps->normWindow2 = 2.75*R2; + float radMoment1 = stamps->normWindow1 / 2.75; + float radMoment2 = stamps->normWindow2 / 2.75; + pmSubtractionParamsScale(NULL, NULL, isisWidths, radMoment1, radMoment2); + + // float maxFWHM = PS_MAX(fwhm1, fwhm2); + // float maxRadial = PS_MAX(radMoment1, radMoment2); + + // if (fabs(2.0*(maxFWHM - maxRadial)/(maxFWHM + maxRadial)) > 0.25) { + // if (1) { + // + // float scale = maxRadial / maxFWHM; + // psLogMsg ("psModules.imcombine", PS_LOG_INFO, "Kron and FWHM scales are quite different, re-scale by %f to use Kron", scale); + // + // for (int i = 0; i < isisWidths->n; i++) { + // isisWidths->data.F32[i] *= scale; + // } + // } + } + + // Define kernel basis functions + if (optimum && (type == PM_SUBTRACTION_KERNEL_ISIS || type == PM_SUBTRACTION_KERNEL_GUNK)) { + kernels = pmSubtractionKernelsOptimumISIS(type, size, inner, spatialOrder, + optFWHMs, optOrder, stamps, footprint, + optThreshold, penalty, bounds, subMode); + if (!kernels) { + psErrorClear(); + psWarning("Unable to derive optimum ISIS kernel --- switching to default."); + } + } + if (kernels == NULL) { + // Not an ISIS/GUNK kernel, or the optimum kernel search failed + kernels = pmSubtractionKernelsGenerate(type, size, spatialOrder, isisWidths, isisOrders, + inner, binning, ringsOrder, penalty, bounds, subMode); + } + + memCheck("kernels"); + + if (subMode == PM_SUBTRACTION_MODE_UNSURE) { + pmSubtractionMode newMode = pmSubtractionBestMode(&stamps, &kernels, subMask, rej); + switch (newMode) { + case PM_SUBTRACTION_MODE_1: + psLogMsg("psModules.imcombine", PS_LOG_INFO, "Convolving image 1 to match image 2."); + break; + case PM_SUBTRACTION_MODE_2: + psLogMsg("psModules.imcombine", PS_LOG_INFO, "Convolving image 2 to match image 1."); + break; + default: + psError(psErrorCodeLast(), false, "Unable to determine subtraction order."); + goto MATCH_ERROR; + } + subMode = newMode; + } + + int numRejected = -1; // Number of rejected stamps in each iteration + for (int k = 0; (k < iter) && (numRejected != 0); k++) { + psLogMsg("psModules.imcombine", PS_LOG_INFO, "Iteration %d.", k); + + bool tryAgain = true; + while (tryAgain) { + // We get the stamps here; we will also attempt to get stamps at the first iteration, but it + // doesn't matter. + if (!pmSubtractionStampsSelect(&stamps, ro1, ro2, subMask, variance, region, stampThresh1, stampThresh2, + stampSpacing, normFrac, sysError, skyError, size, footprint, subMode)) { + goto MATCH_ERROR; + } + + // generate the window function from the set of stamps + if (!pmSubtractionStampsGetWindow(&tryAgain, stamps, size)) { + // if we failed, it might be due to the desired normWindow being larger than the current footprint. + // in this case, just adjust the footprint and try again. + if (tryAgain) { + footprint = PS_MAX(stamps->normWindow1, stamps->normWindow2) + 2; + + // we need to reconstruct everything, so just free the stamps here and retry + psFree(stamps); + } else { + // unrecoverable error + psError(psErrorCodeLast(), false, "Unable to get stamp window."); + goto MATCH_ERROR; + } + } + } // step 0 : calculate the normalizations, pass along to the next steps via stamps->normValue - psTrace("psModules.imcombine", 3, "Calculating normalization...\n"); - if (!pmSubtractionCalculateNormalization(stamps, kernels->mode)) { - psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation."); - goto MATCH_ERROR; - } - - // 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."); - 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"); - - 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); - } - psFree(stamps); - stamps = NULL; - - memCheck("solution"); - - if (!pmSubtractionAnalysis(analysis, header, kernels, region, numCols, numRows)) { - psError(psErrorCodeLast(), false, "Unable to generate QA data"); - goto MATCH_ERROR; - } - memCheck("diag outputs"); - - psTrace("psModules.imcombine", 2, "Convolving...\n"); - if (!pmSubtractionConvolve(conv1, conv2, ro1, ro2, subMask, stride, maskBad, maskPoor, poorFrac, - kernelError, covarFrac, region, kernels, true, useFFT)) { - psError(psErrorCodeLast(), false, "Unable to convolve image."); - goto MATCH_ERROR; - } - - psFree(kernels); - kernels = NULL; - } + psTrace("psModules.imcombine", 3, "Calculating normalization...\n"); + if (!pmSubtractionCalculateNormalization(stamps, kernels->mode)) { + psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation."); + goto MATCH_ERROR; + } + + // 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; + } + memCheck(" reject stamps"); + } + + // apply the best fit so we are ready to roll + psLogMsg("psModules.imcombine", PS_LOG_INFO, "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); + psFree(bestMatch); + memCheck("solution"); + + if (!pmSubtractionAnalysis(analysis, header, kernels, region, numCols, numRows)) { + psError(psErrorCodeLast(), false, "Unable to generate QA data"); + goto MATCH_ERROR; + } + memCheck("diag outputs"); + + psTrace("psModules.imcombine", 2, "Convolving...\n"); + if (!pmSubtractionConvolve(conv1, conv2, ro1, ro2, subMask, stride, maskBad, maskPoor, poorFrac, + kernelError, covarFrac, region, kernels, true, useFFT)) { + psError(psErrorCodeLast(), false, "Unable to convolve image."); + goto MATCH_ERROR; + } + + psFree(kernels); + kernels = NULL; + } } psFree(rng); @@ -816,10 +858,10 @@ if (conv1 && !pmSubtractionBorder(conv1->image, conv1->variance, conv1->mask, size, maskBad)) { - psError(psErrorCodeLast(), false, "Unable to set border of convolved image."); - goto MATCH_ERROR; + psError(psErrorCodeLast(), false, "Unable to set border of convolved image."); + goto MATCH_ERROR; } if (conv2 && !pmSubtractionBorder(conv2->image, conv2->variance, conv2->mask, size, maskBad)) { - psError(psErrorCodeLast(), false, "Unable to set border of convolved image."); - goto MATCH_ERROR; + psError(psErrorCodeLast(), false, "Unable to set border of convolved image."); + goto MATCH_ERROR; } @@ -832,15 +874,15 @@ #ifdef TESTING { - if (subMode == PM_SUBTRACTION_MODE_1 || subMode == PM_SUBTRACTION_MODE_DUAL) { - psFits *fits = psFitsOpen("convolved1.fits", "w"); - psFitsWriteImage(fits, NULL, conv1->image, 0, NULL); - psFitsClose(fits); - } - - if (subMode == PM_SUBTRACTION_MODE_2 || subMode == PM_SUBTRACTION_MODE_DUAL) { - psFits *fits = psFitsOpen("convolved2.fits", "w"); - psFitsWriteImage(fits, NULL, conv2->image, 0, NULL); - psFitsClose(fits); - } + if (subMode == PM_SUBTRACTION_MODE_1 || subMode == PM_SUBTRACTION_MODE_DUAL) { + psFits *fits = psFitsOpen("convolved1.fits", "w"); + psFitsWriteImage(fits, NULL, conv1->image, 0, NULL); + psFitsClose(fits); + } + + if (subMode == PM_SUBTRACTION_MODE_2 || subMode == PM_SUBTRACTION_MODE_DUAL) { + psFits *fits = psFitsOpen("convolved2.fits", "w"); + psFitsWriteImage(fits, NULL, conv2->image, 0, NULL); + psFitsClose(fits); + } } #endif @@ -858,4 +900,5 @@ psFree(variance); psFree(rng); + psFree(bestMatch); return false; } @@ -866,8 +909,8 @@ // increment). static int subtractionOrderWidth(const psKernel *kernel, // Image - float bg, // Background in image - int size, // Maximum size - const psArray *models, // Buffer of models - const psVector *modelSums // Buffer of model sums + float bg, // Background in image + int size, // Maximum size + const psArray *models, // Buffer of models + const psVector *modelSums // Buffer of model sums ) { @@ -882,20 +925,20 @@ psVector *chi2 = psVectorAlloc(size, PS_TYPE_F32); // chi^2 as a function of radius for (int sigma = 0; sigma < size; sigma++) { - double sumFG = 0.0; // Sum for calculating the normalisation of the Gaussian - psKernel *model = models->data[sigma]; // Model of interest - for (int y = yMin; y <= yMax; y++) { - for (int x = xMin; x <= xMax; x++) { - sumFG += model->kernel[y][x] * (kernel->kernel[y][x] - bg); - } - } - float norm = sumFG * modelSums->data.F64[sigma]; // Normalisation for Gaussian - double sumDev2 = 0.0; // Sum of square deviations - for (int y = yMin; y <= yMax; y++) { - for (int x = xMin; x <= xMax; x++) { - float dev = kernel->kernel[y][x] - bg - norm * model->kernel[y][x]; // Deviation - sumDev2 += PS_SQR(dev); - } - } - chi2->data.F32[sigma] = sumDev2; + double sumFG = 0.0; // Sum for calculating the normalisation of the Gaussian + psKernel *model = models->data[sigma]; // Model of interest + for (int y = yMin; y <= yMax; y++) { + for (int x = xMin; x <= xMax; x++) { + sumFG += model->kernel[y][x] * (kernel->kernel[y][x] - bg); + } + } + float norm = sumFG * modelSums->data.F64[sigma]; // Normalisation for Gaussian + double sumDev2 = 0.0; // Sum of square deviations + for (int y = yMin; y <= yMax; y++) { + for (int x = xMin; x <= xMax; x++) { + float dev = kernel->kernel[y][x] - bg - norm * model->kernel[y][x]; // Deviation + sumDev2 += PS_SQR(dev); + } + } + chi2->data.F32[sigma] = sumDev2; } @@ -904,8 +947,8 @@ float bestChi2 = INFINITY; // Best chi^2 for (int i = 0; i < size; i++) { - if (chi2->data.F32[i] < bestChi2) { - bestChi2 = chi2->data.F32[i]; - bestIndex = i; - } + if (chi2->data.F32[i] < bestChi2) { + bestChi2 = chi2->data.F32[i]; + bestIndex = i; + } } psFree(chi2); @@ -916,6 +959,6 @@ bool pmSubtractionOrderStamp(psVector *ratios, psVector *mask, const pmSubtractionStampList *stamps, - const psArray *models, const psVector *modelSums, - int index, float bg1, float bg2) + const psArray *models, const psVector *modelSums, + int index, float bg1, float bg2) { PS_ASSERT_VECTOR_NON_NULL(ratios, false); @@ -931,5 +974,5 @@ pmSubtractionStamp *stamp = stamps->stamps->data[index]; // Stamp of interest psAssert(stamp->status == PM_SUBTRACTION_STAMP_CALCULATE || stamp->status == PM_SUBTRACTION_STAMP_USED, - "We checked this earlier."); + "We checked this earlier."); // Widths of stars @@ -938,11 +981,11 @@ if (width1 == 0 || width2 == 0) { - ratios->data.F32[index] = NAN; - mask->data.PS_TYPE_VECTOR_MASK_DATA[index] = 0xff; + ratios->data.F32[index] = NAN; + mask->data.PS_TYPE_VECTOR_MASK_DATA[index] = 0xff; } else { - ratios->data.F32[index] = (float)width1 / (float)width2; - mask->data.PS_TYPE_VECTOR_MASK_DATA[index] = 0; - psTrace("psModules.imcombine", 3, "Stamp %d (%.1f,%.1f) widths: %d, %d --> %f\n", - index, stamp->x, stamp->y, width1, width2, ratios->data.F32[index]); + ratios->data.F32[index] = (float)width1 / (float)width2; + mask->data.PS_TYPE_VECTOR_MASK_DATA[index] = 0; + psTrace("psModules.imcombine", 3, "Stamp %d (%.1f,%.1f) widths: %d, %d --> %f\n", + index, stamp->x, stamp->y, width1, width2, ratios->data.F32[index]); } @@ -978,60 +1021,60 @@ psVector *modelSums = psVectorAlloc(size, PS_TYPE_F64); // Gaussian model sums for (int sigma = 0; sigma < size; sigma++) { - psKernel *model = psKernelAlloc(-size, size, -size, size); // Gaussian model - float invSigma2 = 1.0 / (float)PS_SQR(1 + sigma); // Inverse sigma squared - double sumGG = 0.0; // Sum of square of Gaussian - for (int y = -size; y <= size; y++) { - int y2 = PS_SQR(y); // y squared - for (int x = -size; x <= size; x++) { - float rad2 = PS_SQR(x) + y2; // Radius squared - float value = expf(-rad2 * invSigma2); // Model value - model->kernel[y][x] = value; - sumGG += PS_SQR(value); - } - } - models->data[sigma] = model; - modelSums->data.F64[sigma] = 1.0 / sumGG; + psKernel *model = psKernelAlloc(-size, size, -size, size); // Gaussian model + float invSigma2 = 1.0 / (float)PS_SQR(1 + sigma); // Inverse sigma squared + double sumGG = 0.0; // Sum of square of Gaussian + for (int y = -size; y <= size; y++) { + int y2 = PS_SQR(y); // y squared + for (int x = -size; x <= size; x++) { + float rad2 = PS_SQR(x) + y2; // Radius squared + float value = expf(-rad2 * invSigma2); // Model value + model->kernel[y][x] = value; + sumGG += PS_SQR(value); + } + } + models->data[sigma] = model; + modelSums->data.F64[sigma] = 1.0 / sumGG; } // Fit models to stamps for (int i = 0; i < stamps->num; i++) { - pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest - if (stamp->status != PM_SUBTRACTION_STAMP_CALCULATE && stamp->status != PM_SUBTRACTION_STAMP_USED) { - mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = 0xff; - continue; - } - - if (pmSubtractionThreaded()) { - psThreadJob *job = psThreadJobAlloc("PSMODULES_SUBTRACTION_ORDER"); - psArrayAdd(job->args, 1, ratios); - psArrayAdd(job->args, 1, mask); - psArrayAdd(job->args, 1, stamps); - psArrayAdd(job->args, 1, models); - psArrayAdd(job->args, 1, modelSums); - PS_ARRAY_ADD_SCALAR(job->args, i, PS_TYPE_S32); - PS_ARRAY_ADD_SCALAR(job->args, bg1, PS_TYPE_F32); - PS_ARRAY_ADD_SCALAR(job->args, bg2, PS_TYPE_F32); - if (!psThreadJobAddPending(job)) { - return false; - } - } else { - if (!pmSubtractionOrderStamp(ratios, mask, stamps, models, modelSums, i, bg1, bg2)) { - psError(psErrorCodeLast(), false, "Unable to measure PSF width for stamp %d", i); - psFree(models); - psFree(modelSums); - psFree(ratios); - psFree(mask); - return false; - } - } + pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest + if (stamp->status != PM_SUBTRACTION_STAMP_CALCULATE && stamp->status != PM_SUBTRACTION_STAMP_USED) { + mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = 0xff; + continue; + } + + if (pmSubtractionThreaded()) { + psThreadJob *job = psThreadJobAlloc("PSMODULES_SUBTRACTION_ORDER"); + psArrayAdd(job->args, 1, ratios); + psArrayAdd(job->args, 1, mask); + psArrayAdd(job->args, 1, stamps); + psArrayAdd(job->args, 1, models); + psArrayAdd(job->args, 1, modelSums); + PS_ARRAY_ADD_SCALAR(job->args, i, PS_TYPE_S32); + PS_ARRAY_ADD_SCALAR(job->args, bg1, PS_TYPE_F32); + PS_ARRAY_ADD_SCALAR(job->args, bg2, PS_TYPE_F32); + if (!psThreadJobAddPending(job)) { + return false; + } + } else { + if (!pmSubtractionOrderStamp(ratios, mask, stamps, models, modelSums, i, bg1, bg2)) { + psError(psErrorCodeLast(), false, "Unable to measure PSF width for stamp %d", i); + psFree(models); + psFree(modelSums); + psFree(ratios); + psFree(mask); + return false; + } + } } if (!psThreadPoolWait(true)) { - psError(psErrorCodeLast(), false, "Error waiting for threads."); - psFree(models); - psFree(modelSums); - psFree(ratios); - psFree(mask); - return false; + psError(psErrorCodeLast(), false, "Error waiting for threads."); + psFree(models); + psFree(modelSums); + psFree(ratios); + psFree(mask); + return false; } @@ -1041,9 +1084,9 @@ psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN); if (!psVectorStats(stats, ratios, NULL, mask, 0xff)) { - psError(psErrorCodeLast(), false, "Unable to calculate statistics for moments ratio."); - psFree(mask); - psFree(ratios); - psFree(stats); - return PM_SUBTRACTION_MODE_ERR; + psError(psErrorCodeLast(), false, "Unable to calculate statistics for moments ratio."); + psFree(mask); + psFree(ratios); + psFree(stats); + return PM_SUBTRACTION_MODE_ERR; } psFree(ratios); @@ -1052,6 +1095,6 @@ // XXX raise an error here or not? if (isnan(stats->robustMedian)) { - psFree(stats); - return PM_SUBTRACTION_MODE_ERR; + psFree(stats); + return PM_SUBTRACTION_MODE_ERR; } @@ -1066,23 +1109,31 @@ // Test a subtraction mode by performing a single iteration static bool subtractionModeTest(pmSubtractionStampList *stamps, // Stamps to use to find best mode - pmSubtractionKernels *kernels, // Kernel description - const char *description, // Description for trace - psImage *subMask, // Subtraction mask - float rej // Rejection threshold - ) + pmSubtractionKernels *kernels, // Kernel description + const char *description, // Description for trace + psImage *subMask, // Subtraction mask + float rej // Rejection threshold + ) { assert(stamps); 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)) { - psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation."); - return false; + if (!pmSubtractionCalculateEquation(stamps, kernels)) { + psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation."); + return false; } psTrace("psModules.imcombine", 3, "Solving %s normalization equation...\n", description); - if (!pmSubtractionSolveEquation(kernels, stamps, SUBMODE)) { - psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation."); - return false; + if (!pmSubtractionSolveEquation(kernels, stamps)) { + psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation."); + return false; } @@ -1090,47 +1141,48 @@ psVector *deviations = pmSubtractionCalculateDeviations(stamps, kernels); // Stamp deviations if (!deviations) { - psError(psErrorCodeLast(), false, "Unable to calculate deviations."); - return false; - } - + psError(psErrorCodeLast(), false, "Unable to calculate deviations."); + return false; + } + + // 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); if (numRejected < 0) { - psError(psErrorCodeLast(), false, "Unable to reject stamps."); - psFree(deviations); - return false; + psError(psErrorCodeLast(), false, "Unable to reject stamps."); + psFree(deviations); + return false; } psFree(deviations); if (numRejected > 0) { - // 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)) { - psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation."); - return false; - } - - psTrace("psModules.imcombine", 3, "Resolving %s equation...\n", description); - if (!pmSubtractionSolveEquation(kernels, stamps, PM_SUBTRACTION_EQUATION_ALL)) { - psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation."); - return false; - } - psTrace("psModules.imcombine", 3, "Recalculate %s deviations...\n", description); - - psVector *deviations = pmSubtractionCalculateDeviations(stamps, kernels); // Stamp deviations - if (!deviations) { - psError(psErrorCodeLast(), false, "Unable to calculate deviations."); - return false; - } - psTrace("psModules.imcombine", 3, "Measuring %s quality...\n", description); - long numRejected = pmSubtractionRejectStamps(kernels, stamps, deviations, subMask, NAN); - if (numRejected < 0) { - psError(psErrorCodeLast(), false, "Unable to reject stamps."); - psFree(deviations); - return false; - } - psFree(deviations); - } - + // Allow re-fit with reduced stamps set + psTrace("psModules.imcombine", 3, "Calculating %s normalization equation...\n", description); + if (!pmSubtractionCalculateEquation(stamps, kernels)) { + psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation."); + return false; + } + + psTrace("psModules.imcombine", 3, "Resolving %s equation...\n", description); + if (!pmSubtractionSolveEquation(kernels, stamps)) { + psError(psErrorCodeLast(), false, "Unable to calculate least-squares equation."); + return false; + } + psTrace("psModules.imcombine", 3, "Recalculate %s deviations...\n", description); + + psVector *deviations = pmSubtractionCalculateDeviations(stamps, kernels); // Stamp deviations + if (!deviations) { + psError(psErrorCodeLast(), false, "Unable to calculate deviations."); + return false; + } + psTrace("psModules.imcombine", 3, "Measuring %s quality...\n", description); + long numRejected = pmSubtractionRejectStamps(kernels, stamps, deviations, subMask, NAN); + if (numRejected < 0) { + psError(psErrorCodeLast(), false, "Unable to reject stamps."); + psFree(deviations); + return false; + } + psFree(deviations); + } +# endif return true; } @@ -1138,5 +1190,5 @@ pmSubtractionMode pmSubtractionBestMode(pmSubtractionStampList **stamps, pmSubtractionKernels **kernels, - const psImage *subMask, float rej) + const psImage *subMask, float rej) { PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(*stamps, PM_SUBTRACTION_MODE_ERR); @@ -1150,9 +1202,9 @@ if (!subtractionModeTest(stamps1, kernels1, "convolve 1", subMask1, rej)) { - psError(psErrorCodeLast(), false, "Unable to test subtraction with convolution of image 1"); - psFree(stamps1); - psFree(kernels1); - psFree(subMask1); - return PM_SUBTRACTION_MODE_ERR; + psError(psErrorCodeLast(), false, "Unable to test subtraction with convolution of image 1"); + psFree(stamps1); + psFree(kernels1); + psFree(subMask1); + return PM_SUBTRACTION_MODE_ERR; } psFree(subMask1); @@ -1165,11 +1217,11 @@ if (!subtractionModeTest(stamps2, kernels2, "convolve 2", subMask2, rej)) { - psError(psErrorCodeLast(), false, "Unable to test subtraction with convolution of image 2"); - psFree(stamps2); - psFree(kernels2); - psFree(subMask2); - psFree(stamps1); - psFree(kernels1); - return PM_SUBTRACTION_MODE_ERR; + psError(psErrorCodeLast(), false, "Unable to test subtraction with convolution of image 2"); + psFree(stamps2); + psFree(kernels2); + psFree(subMask2); + psFree(stamps1); + psFree(kernels1); + return PM_SUBTRACTION_MODE_ERR; } psFree(subMask2); @@ -1179,14 +1231,14 @@ pmSubtractionKernels *bestKernels = NULL; // Best choice for kernels psLogMsg("psModules.imcombine", PS_LOG_INFO, - "Image 1: %f +/- %f from %d stamps\nImage 2: %f +/- %f from %d stamps\n", - kernels1->mean, kernels1->rms, kernels1->numStamps, - kernels2->mean, kernels2->rms, kernels2->numStamps); + "Image 1: %f +/- %f from %d stamps\nImage 2: %f +/- %f from %d stamps\n", + kernels1->mean, kernels1->rms, kernels1->numStamps, + kernels2->mean, kernels2->rms, kernels2->numStamps); if (kernels1->mean < kernels2->mean) { - bestStamps = stamps1; - bestKernels = kernels1; + bestStamps = stamps1; + bestKernels = kernels1; } else { - bestStamps = stamps2; - bestKernels = kernels2; + bestStamps = stamps2; + bestKernels = kernels2; } @@ -1204,42 +1256,62 @@ } - -bool pmSubtractionParamsScale(int *kernelSize, int *stampSize, psVector *widths, - float scaleRef, float scaleMin, float scaleMax) +static float scaleRefOption = NAN; +static float scaleMinOption = NAN; +static float scaleMaxOption = NAN; +static bool scaleOption = false; + +bool pmSubtractionParamScaleOptions(bool scale, float scaleRef, float scaleMin, float scaleMax) { + + if (scale) { + PS_ASSERT_FLOAT_LARGER_THAN(scaleRef, 0.0, false); + PS_ASSERT_FLOAT_LARGER_THAN(scaleMin, 0.0, false); + PS_ASSERT_FLOAT_LARGER_THAN(scaleMax, 0.0, false); + PS_ASSERT_FLOAT_LARGER_THAN(scaleMax, scaleMin, false); + } + + scaleRefOption = scaleRef; + scaleMinOption = scaleMin; + scaleMaxOption = scaleMax; + scaleOption = scale; + + return true; +} + +bool pmSubtractionParamsScale(int *kernelSize, int *stampSize, psVector *widths, float fwhm1, float fwhm2) { - PS_ASSERT_PTR_NON_NULL(kernelSize, false); - PS_ASSERT_PTR_NON_NULL(stampSize, false); + // PS_ASSERT_PTR_NON_NULL(kernelSize, false); + // PS_ASSERT_PTR_NON_NULL(stampSize, false); PS_ASSERT_VECTOR_NON_NULL(widths, false); PS_ASSERT_VECTOR_TYPE(widths, PS_TYPE_F32, false); - PS_ASSERT_FLOAT_LARGER_THAN(scaleRef, 0.0, false); - PS_ASSERT_FLOAT_LARGER_THAN(scaleMin, 0.0, false); - PS_ASSERT_FLOAT_LARGER_THAN(scaleMax, 0.0, false); - PS_ASSERT_FLOAT_LARGER_THAN(scaleMax, scaleMin, false); - - float fwhm1; - float fwhm2; - - pmSubtractionGetFWHMs(&fwhm1, &fwhm2); - psAssert(isfinite(fwhm1), "fwhm 1 not set"); - psAssert(isfinite(fwhm2), "fwhm 2 not set"); + + if (!scaleOption) return true; + + // pmSubtractionGetFWHMs(&fwhm1, &fwhm2); + // psAssert(isfinite(fwhm1), "fwhm 1 not set"); + // psAssert(isfinite(fwhm2), "fwhm 2 not set"); // float diff = sqrtf(PS_SQR(PS_MAX(fwhm1, fwhm2)) - PS_SQR(PS_MIN(fwhm1, fwhm2))); // Difference - float scale = PS_MAX(fwhm1, fwhm2) / scaleRef; // Scaling factor - - if (isfinite(scaleMin) && scale < scaleMin) { - scale = scaleMin; - } - if (isfinite(scaleMax) && scale > scaleMax) { - scale = scaleMax; + float scale = PS_MAX(fwhm1, fwhm2) / scaleRefOption; // Scaling factor + + if (isfinite(scaleMinOption) && scale < scaleMinOption) { + scale = scaleMinOption; + } + if (isfinite(scaleMaxOption) && scale > scaleMaxOption) { + scale = scaleMaxOption; } for (int i = 0; i < widths->n; i++) { - widths->data.F32[i] *= scale; - } - *kernelSize = *kernelSize * scale + 0.5; - *stampSize = *stampSize * scale + 0.5; - - psLogMsg("psModules.imcombine", PS_LOG_INFO, - "Scaling kernel parameters by %f: %d %d", scale, *kernelSize, *stampSize); + widths->data.F32[i] *= scale; + } + if (kernelSize) { + *kernelSize = *kernelSize * scale + 0.5; + } + if (stampSize) { + *stampSize = *stampSize * scale + 0.5; + } + + psLogMsg("psModules.imcombine", PS_LOG_INFO, "Scaling kernel parameters by %f", scale); + if (kernelSize) psLogMsg("psModules.imcombine", PS_LOG_INFO, " modified kernel size %d", *kernelSize); + if (stampSize) psLogMsg("psModules.imcombine", PS_LOG_INFO, " modified stamp size %d", *stampSize); return true; Index: trunk/psModules/src/imcombine/pmSubtractionMatch.h =================================================================== --- trunk/psModules/src/imcombine/pmSubtractionMatch.h (revision 30621) +++ trunk/psModules/src/imcombine/pmSubtractionMatch.h (revision 30622) @@ -8,4 +8,5 @@ #include #include +#include /// Match two images @@ -26,5 +27,5 @@ int size, ///< Kernel half-size int order, ///< Spatial polynomial order - const psVector *widths, ///< ISIS Gaussian widths + psVector *widths, ///< ISIS Gaussian widths const psVector *orders, ///< ISIS Polynomial orders int inner, ///< Inner radius for various kernel types @@ -100,11 +101,24 @@ /// Scale subtraction parameters according to the FWHMs of the inputs -bool pmSubtractionParamsScale( - int *kernelSize, ///< Half-size of the kernel - int *stampSize, ///< Half-size of the stamp (footprint) - psVector *widths, ///< ISIS widths - float scaleRef, ///< Reference width for scaling - float scaleMin, ///< Minimum scaling ratio, or NAN - float scaleMax ///< Maximum scaling ratio, or NAN +// bool pmSubtractionParamsScale( +// int *kernelSize, ///< Half-size of the kernel +// int *stampSize, ///< Half-size of the stamp (footprint) +// psVector *widths, ///< ISIS widths +// float scaleRef, ///< Reference width for scaling +// float scaleMin, ///< Minimum scaling ratio, or NAN +// float scaleMax ///< Maximum scaling ratio, or NAN +// ); + +bool pmSubtractionParamsScale(int *kernelSize, int *stampSize, psVector *widths, float fwhm1, float fwhm2); + +bool pmSubtractionParamScaleOptions(bool scale, float scaleRef, float scaleMin, float scaleMax); + +bool pmSubtractionMatchAttempt( + pmSubtractionQuality **bestMatch, + pmSubtractionKernels *kernels, + pmSubtractionStampList *stamps, + pmSubtractionMode mode, + int spatialOrder, + bool final ); Index: trunk/psModules/src/imcombine/pmSubtractionParams.c =================================================================== --- trunk/psModules/src/imcombine/pmSubtractionParams.c (revision 30621) +++ trunk/psModules/src/imcombine/pmSubtractionParams.c (revision 30622) @@ -9,6 +9,8 @@ #include "pmErrorCodes.h" +#include "pmFPA.h" +#include "pmSubtractionTypes.h" +#include "pmSubtraction.h" #include "pmSubtractionStamps.h" -#include "pmSubtraction.h" #include "pmSubtractionParams.h" @@ -16,4 +18,5 @@ #if 0 +// XXX this was moved to pmSubtraction.c in r15443 -- delete // Convolve the reference stamp by the kernel static psKernel *convolveStamp(const pmSubtractionStamp *stamp, // Stamp to be convolved Index: trunk/psModules/src/imcombine/pmSubtractionStamps.c =================================================================== --- trunk/psModules/src/imcombine/pmSubtractionStamps.c (revision 30621) +++ trunk/psModules/src/imcombine/pmSubtractionStamps.c (revision 30622) @@ -27,6 +27,8 @@ #include "pmSource.h" +#include "pmSubtractionTypes.h" #include "pmSubtraction.h" #include "pmSubtractionStamps.h" +#include "pmSubtractionVisual.h" #define STAMP_LIST_BUFFER 20 // Number of stamps to add to list at a time @@ -122,4 +124,5 @@ if ((image1 && image1->data.F32[y][x] < thresh1) || (image2 && image2->data.F32[y][x] < thresh2)) { + // fprintf (stderr, "%f,%f : thresh\n", xRaw, yRaw); continue; } @@ -366,10 +369,63 @@ } - -pmSubtractionStampList *pmSubtractionStampsFind(pmSubtractionStampList *stamps, const psImage *image1, - const psImage *image2, const psImage *subMask, - const psRegion *region, float thresh1, float thresh2, - int size, int footprint, float spacing, float normFrac, - float sysErr, float skyErr, pmSubtractionMode mode) +bool pmSubtractionStampsSelect(pmSubtractionStampList **stamps, // Stamps to read + const pmReadout *ro1, // Readout 1 + const pmReadout *ro2, // Readout 2 + const psImage *subMask, // Mask for subtraction, or NULL + psImage *variance, // Variance map + const psRegion *region, // Region of interest + float thresh1, // Threshold for stamp finding on readout 1 + float thresh2, // Threshold for stamp finding on readout 2 + float stampSpacing, // Spacing between stamps + float normFrac, // Fraction of flux in window for normalisation window + float sysError, // Relative systematic error in images + float skyError, // Relative systematic error in images + int size, // Kernel half-size + int footprint, // Convolution footprint for stamps + pmSubtractionMode mode // Mode for subtraction + ) +{ + PS_ASSERT_PTR_NON_NULL(stamps, false); + PM_ASSERT_READOUT_NON_NULL(ro1, false); + PM_ASSERT_READOUT_NON_NULL(ro2, false); + PS_ASSERT_IMAGE_NON_NULL(subMask, false); + PS_ASSERT_IMAGE_NON_NULL(variance, false); + PS_ASSERT_PTR_NON_NULL(region, false); + + psTrace("psModules.imcombine", 3, "Finding stamps...\n"); + + psImage *image1 = ro1 ? ro1->image : NULL, *image2 = ro2 ? ro2->image : NULL; // Images of interest + + *stamps = pmSubtractionStampsFind(*stamps, image1, image2, subMask, region, thresh1, thresh2, + size, footprint, stampSpacing, normFrac, sysError, skyError, mode); + if (!*stamps) { + psError(psErrorCodeLast(), false, "Unable to find stamps."); + return false; + } + + psTrace("psModules.imcombine", 3, "Extracting stamps...\n"); + if (!pmSubtractionStampsExtract(*stamps, ro1->image, ro2->image, variance, size, *region)) { + psError(psErrorCodeLast(), false, "Unable to extract stamps."); + return false; + } + + pmSubtractionVisualPlotStamps(*stamps, (pmReadout *) ro1); + return true; +} + +pmSubtractionStampList *pmSubtractionStampsFind(pmSubtractionStampList *stamps, + const psImage *image1, + const psImage *image2, + const psImage *subMask, + const psRegion *region, + float thresh1, + float thresh2, + int size, + int footprint, + float spacing, + float normFrac, + float sysErr, + float skyErr, + pmSubtractionMode mode) { if (!image1 && !image2) { @@ -429,4 +485,19 @@ stamps = pmSubtractionStampListAlloc(numCols, numRows, region, footprint, spacing, normFrac, sysErr, skyErr); + } + + // XXX TEST : dump all stars in the stamps here + if (0) { + FILE *f = fopen ("stamp.dat", "w"); + for (int i = 0; i < stamps->num; i++) { + psVector *xList = stamps->x->data[i]; + psVector *yList = stamps->y->data[i]; // Coordinate lists + psVector *fluxList = stamps->flux->data[i]; // List of stamp fluxes + + for (int j = 0; j < xList->n; j++) { + fprintf (f, "%d %d %f %f %f\n", i, j, xList->data.F32[j], yList->data.F32[j], fluxList->data.F32[j]); + } + } + fclose (f); } @@ -636,4 +707,6 @@ } + int nTotal = 0; + // Sort the list by flux, with the brightest last for (int i = 0; i < numStamps; i++) { @@ -662,17 +735,35 @@ stamps->y->data[i] = ySorted; stamps->flux->data[i] = fluxSorted; - } - + nTotal += num; + } + // fprintf (stderr, "nTotal %d\n", nTotal); + return stamps; } - -bool pmSubtractionStampsGetWindow(pmSubtractionStampList *stamps, int kernelSize) +// we are essentially using aperture photometry to determine the photometric match between the +// images. we need to choose an appropriate-sized aperture for this analysis. If it is too +// large, the measurement will be noisy (and possibly biased) due to the sky noise. If it is +// too small, or inconsistent, the measurement will be biased. We use Kron-mag like aperture +// scaled by the first radial moment. +bool pmSubtractionStampsGetWindow(bool *tryAgain, pmSubtractionStampList *stamps, int kernelSize) { PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(stamps, false); PS_ASSERT_INT_NONNEGATIVE(kernelSize, false); + // if we succeed, or fail with an unrecoverable error, do not try again + if (tryAgain) { + *tryAgain = false; + } + int size = stamps->footprint; // Size of postage stamps + // window for moments calculations downstream + psFree (stamps->window); + stamps->window = psKernelAlloc(-size, size, -size, size); + psImageInit(stamps->window->image, 0.0); + + // window1 and window2 are mean stamp images used here to measure the + // first radial moment, and thus the normalization window psFree (stamps->window1); stamps->window1 = psKernelAlloc(-size, size, -size, size); @@ -683,22 +774,16 @@ psImageInit(stamps->window2->image, 0.0); - // Generate a weighting window based on the fwhms (20% larger than the largest) - { - float fwhm1, fwhm2; - - // XXX this is annoyingly hack-ish - pmSubtractionGetFWHMs(&fwhm1, &fwhm2); - - float sigma = 1.5 * PS_MAX(fwhm1, fwhm2) / 2.35; - - psFree (stamps->window); - stamps->window = psKernelAlloc(-size, size, -size, size); - psImageInit(stamps->window->image, 0.0); - - for (int y = -size; y <= size; y++) { - for (int x = -size; x <= size; x++) { - stamps->window->kernel[y][x] = exp(-0.5*(x*x + y*y)/(sigma*sigma)); - } - } + // Generate an initial weighting window based on the fwhms (50% larger than the largest) + float fwhm1, fwhm2; + + // XXX this is annoyingly hack-ish + pmSubtractionGetFWHMs(&fwhm1, &fwhm2); + + float sigma = 1.5 * PS_MAX(fwhm1, fwhm2) / 2.35; + + for (int y = -size; y <= size; y++) { + for (int x = -size; x <= size; x++) { + stamps->window->kernel[y][x] = exp(-0.5*(x*x + y*y)/(sigma*sigma)); + } } @@ -790,6 +875,5 @@ psTrace("psModules.imcombine", 3, "Window total (2): %f, threshold: %f\n", sum2, (1.0 - stamps->normFrac) * sum2); -# if (1) - // this block attempts to calculate the radius based on the first radial moment + // attempt to calculate the normalization window based on the first radial moment double Sr1 = 0.0; double Sr2 = 0.0; @@ -809,71 +893,56 @@ float R2 = Sr2 / Sf2; - stamps->normWindow1 = 2.0*R1; - stamps->normWindow2 = 2.0*R2; - psLogMsg ("psModules.imcombine", PS_LOG_DETAIL, "Kron Radii: %f for 1, %f for 2\n", stamps->normWindow1, stamps->normWindow2); - -# else - // XXX : this block attempts to calculate the radius by looking at the curve of growth (or something vaguely equivalent). - // It did not do very well (though a true curve-of-growth analysis might be better...) - bool done1 = false; - bool done2 = false; - double prior1 = 0.0; - double prior2 = 0.0; - double delta1o = 1.0; - double delta2o = 1.0; - for (int radius = 1; radius <= size && !(done1 && done2); radius++) { - double within1 = 0.0; - double within2 = 0.0; - for (int y = -radius; y <= radius; y++) { - for (int x = -radius; x <= radius; x++) { - if (PS_SQR(x) + PS_SQR(y) <= PS_SQR(radius)) { - within1 += stamps->window1->kernel[y][x]; - } - if (PS_SQR(x) + PS_SQR(y) <= PS_SQR(radius)) { - within2 += stamps->window2->kernel[y][x]; - } - } - } - double delta1 = (within1 - prior1) / within1; - if (!done1 && (fabs(delta1) < stamps->normFrac)) { - // interpolate to the radius at which delta2 is normFrac: - stamps->normWindow1 = radius - (stamps->normFrac - delta1) / (delta1o - delta1); - psLogMsg ("psModules.imcombine", PS_LOG_DETAIL, "choosing %f (%d) for 1 (%f : %f)\n", stamps->normWindow1, radius, within1, delta1); - done1 = true; - } - double delta2 = (within2 - prior2) / within2; - if (!done2 && (fabs(delta2) < stamps->normFrac)) { - // interpolate to the radius at which delta2 is normFrac: - stamps->normWindow2 = radius - (stamps->normFrac - delta2) / (delta2o - delta2); - psLogMsg ("psModules.imcombine", PS_LOG_DETAIL, "choosing %f (%d) for 2 (%f : %f)\n", stamps->normWindow2, radius, within2, delta2); - done2 = true; - } - psTrace("psModules.imcombine", 5, "Radius %d: %f (%f) and %f (%f)\n", radius, within1, delta1, within2, delta2); - - prior1 = within1; - prior2 = within2; - delta1o = delta1; - delta2o = delta2; - - // if (!done1 && (within1 > (1.0 - stamps->normFrac) * sum1)) { - // stamps->normWindow1 = radius; - // done1 = true; - // } - // if (!done2 && (within2 > (1.0 - stamps->normFrac) * sum2)) { - // stamps->normWindow2 = radius; - // done2 = true; - // } - - } -# endif - - psTrace("psModules.imcombine", 3, "Normalisation window radii set to %f and %f\n", stamps->normWindow1, stamps->normWindow2); - if (stamps->normWindow1 == 0 || stamps->normWindow1 >= size) { + // Compare the Kron Radii (R1 & R2) to above to the FWHMs : if they are too discrepant, we will need to rescale + psLogMsg ("psModules.imcombine", PS_LOG_DETAIL, "Kron Radii vs FWHMs 1: fwhm: %f, kron %f\n", fwhm1, R1); + psLogMsg ("psModules.imcombine", PS_LOG_DETAIL, "Kron Radii vs FWHMs 2: fwhm: %f, kron %f\n", fwhm2, R2); + + // XXX CAREFUL : in pmSubtractionMatch.c:703, we rely on this factor of 2.75.. + stamps->normWindow1 = 2.75*R1; + stamps->normWindow2 = 2.75*R2; + psLogMsg ("psModules.imcombine", PS_LOG_DETAIL, "Windows from Kron Radii: %f for 1, %f for 2\n", stamps->normWindow1, stamps->normWindow2); + + // if the calculated normWindows are too large, we will fall off the stamps. In this case, we need to try again. + if ((stamps->normWindow1 > size) || (stamps->normWindow2 > size)) { + if (tryAgain) { + *tryAgain = true; + } + psFree (stats); + psFree (flux1); + psFree (flux2); + psFree (norm1); + psFree (norm2); + return false; + } + + // this is an unrecoverable error : something really bogus in the data + if (stamps->normWindow1 == 0) { psError(PM_ERR_STAMPS, true, "Unable to determine normalisation window size (1)."); + psFree (stats); + psFree (flux1); + psFree (flux2); + psFree (norm1); + psFree (norm2); return false; } - if (stamps->normWindow2 == 0 || stamps->normWindow2 >= size) { + if (stamps->normWindow2 == 0) { psError(PM_ERR_STAMPS, true, "Unable to determine normalisation window size (2)."); + psFree (stats); + psFree (flux1); + psFree (flux2); + psFree (norm1); + psFree (norm2); return false; + } + + // Generate a weighting window based on the kron radii + float radius = 2.0 * PS_MAX(R1, R2); + psImageInit(stamps->window->image, 0.0); + + // we use a top-hat window for the moments analysis + for (int y = -size; y <= size; y++) { + for (int x = -size; x <= size; x++) { + if (hypot(x,y) > radius) continue; + stamps->window->kernel[y][x] = 1.0; + } } @@ -890,16 +959,4 @@ } } - -#if 0 - { - psFits *fits = NULL; - fits = psFitsOpen ("window1.norm.fits", "w"); - psFitsWriteImage (fits, NULL, stamps->window1->image, 0, NULL); - psFitsClose (fits); - fits = psFitsOpen ("window2.norm.fits", "w"); - psFitsWriteImage (fits, NULL, stamps->window2->image, 0, NULL); - psFitsClose (fits); - } -#endif psFree (stats); @@ -1128,4 +1185,9 @@ continue; } + + // XXX this is somewhat arbitrary... + if (source->errMag > 0.05) continue; + if (fabs(source->psfMag - source->apMag) > 0.5) continue; + if (source->modelPSF) { x->data.F32[numOut] = source->modelPSF->params->data.F32[PM_PAR_XPOS]; @@ -1135,5 +1197,4 @@ y->data.F32[numOut] = source->peak->yf; } - // fprintf (stderr, "stamp: %5.1f %5.1f\n", x->data.F32[numOut], y->data.F32[numOut]); numOut++; } Index: trunk/psModules/src/imcombine/pmSubtractionStamps.h =================================================================== --- trunk/psModules/src/imcombine/pmSubtractionStamps.h (revision 30621) +++ trunk/psModules/src/imcombine/pmSubtractionStamps.h (revision 30622) @@ -1,38 +1,4 @@ #ifndef PM_SUBTRACTION_STAMPS_H #define PM_SUBTRACTION_STAMPS_H - -#include - -#include "pmSubtractionKernels.h" - -/// Status of stamp -typedef enum { - PM_SUBTRACTION_STAMP_INIT, ///< Initial state - PM_SUBTRACTION_STAMP_FOUND, ///< Found a suitable source for this stamp - PM_SUBTRACTION_STAMP_CALCULATE, ///< Calculate matrix and vector values for this stamp - PM_SUBTRACTION_STAMP_USED, ///< Use this stamp - PM_SUBTRACTION_STAMP_REJECTED, ///< This stamp has been rejected - PM_SUBTRACTION_STAMP_NONE ///< No stamp in this region -} pmSubtractionStampStatus; - -/// A list of stamps -typedef struct { - long num; ///< Number of stamps - psArray *stamps; ///< The stamps - psArray *regions; ///< Regions for each stamp - psArray *x, *y; ///< Coordinates for possible stamps (or NULL) - psArray *flux; ///< Fluxes for possible stamps (or NULL) - int footprint; ///< Half-size of stamps - float normFrac; ///< Fraction of flux in window for normalisation window - float normValue; ///< calculated normalization - float normValue2; ///< calculated normalization - psKernel *window1; ///< window function generated from ensemble of stamps (input 1) - psKernel *window2; ///< window function generated from ensemble of stamps (input 2) - psKernel *window; ///< weighting window function (sigma = 1.1 * MAX(fwhm)) - float normWindow1; ///< Size of window for measuring normalisation - float normWindow2; ///< Size of window for measuring normalisation - float sysErr; ///< Systematic error - float skyErr; ///< increase effective readnoise -} pmSubtractionStampList; /// Allocate a list of stamps @@ -74,35 +40,25 @@ ); - -/// A stamp for image subtraction -typedef struct { - float x, y; ///< Position - float flux; ///< Flux - float xNorm, yNorm; ///< Normalised position - psKernel *image1; ///< Reference image postage stamp - psKernel *image2; ///< Input image postage stamp - psKernel *weight; ///< Weight image (1/variance) postage stamp, or NULL - psArray *convolutions1; ///< Convolutions of image 1 for each kernel component, or NULL - psArray *convolutions2; ///< Convolutions of image 2 for each kernel component, or NULL - psImage *matrix; ///< Least-squares matrix, or NULL - psVector *vector; ///< Least-squares vector, or NULL - double norm; ///< Normalisation difference - double normI1; ///< Sum(flux) for image 1 - double normI2; ///< Sum(flux) for image 2 - double normSquare1; ///< Sum(flux^2) for image 1 (used for penalty) - double normSquare2; ///< Sum(flux^2) for image 2 (used for penalty) - pmSubtractionStampStatus status; ///< Status of stamp - psVector *MxxI1; ///< second moments of convolution images - psVector *MyyI1; ///< second moments of convolution images - psVector *MxxI2; ///< second moments of convolution images - psVector *MyyI2; ///< second moments of convolution images - double MxxI1raw; - double MyyI1raw; - double MxxI2raw; - double MyyI2raw; -} pmSubtractionStamp; - /// Allocate a stamp pmSubtractionStamp *pmSubtractionStampAlloc(void); + +// find and extract the stamps +bool pmSubtractionStampsSelect(pmSubtractionStampList **stamps, // Stamps to read + const pmReadout *ro1, // Readout 1 + const pmReadout *ro2, // Readout 2 + const psImage *subMask, // Mask for subtraction, or NULL + psImage *variance, // Variance map + const psRegion *region, // Region of interest + float thresh1, // Threshold for stamp finding on readout 1 + float thresh2, // Threshold for stamp finding on readout 2 + float stampSpacing, // Spacing between stamps + float normFrac, // Fraction of flux in window for normalisation window + float sysError, // Relative systematic error in images + float skyError, // Relative systematic error in images + int size, // Kernel half-size + int footprint, // Convolution footprint for stamps + pmSubtractionMode mode // Mode for subtraction + ); + /// Find stamps on an image @@ -172,4 +128,5 @@ /// Calculate the window and normalisation window from the stamps bool pmSubtractionStampsGetWindow( + bool *tryAgain, ///< re-try with new stamp size? pmSubtractionStampList *stamps, ///< List of stamps int kernelSize ///< Half-size of kernel Index: trunk/psModules/src/imcombine/pmSubtractionThreads.c =================================================================== --- trunk/psModules/src/imcombine/pmSubtractionThreads.c (revision 30621) +++ trunk/psModules/src/imcombine/pmSubtractionThreads.c (revision 30622) @@ -6,4 +6,5 @@ #include +#include "pmSubtractionTypes.h" #include "pmSubtractionMatch.h" #include "pmSubtractionEquation.h" @@ -33,6 +34,13 @@ { - psThreadTask *task = psThreadTaskAlloc("PSMODULES_SUBTRACTION_CALCULATE_EQUATION", 4); + psThreadTask *task = psThreadTaskAlloc("PSMODULES_SUBTRACTION_CALCULATE_EQUATION", 3); task->function = &pmSubtractionCalculateEquationThread; + psThreadTaskAdd(task); + psFree(task); + } + + { + psThreadTask *task = psThreadTaskAlloc("PSMODULES_SUBTRACTION_CONVOLVE_STAMP", 3); + task->function = &pmSubtractionConvolveStampThread; psThreadTaskAdd(task); psFree(task); Index: trunk/psModules/src/imcombine/pmSubtractionTypes.h =================================================================== --- trunk/psModules/src/imcombine/pmSubtractionTypes.h (revision 30622) +++ trunk/psModules/src/imcombine/pmSubtractionTypes.h (revision 30622) @@ -0,0 +1,172 @@ +/* @file pmSubtraction.h + * + * PSF-matched image subtraction, based on the Alard & Lupton (1998) and Alard (2000) methods. + * + * @author Paul Price, IfA + * @author GLG, MHPCC + * + * @version $Revision: 1.36 $ $Name: not supported by cvs2svn $ + * @date $Date: 2009-02-06 02:31:25 $ + * Copyright 2004-207 Institute for Astronomy, University of Hawaii + */ + +#ifndef PM_SUBTRACTION_TYPES_H +#define PM_SUBTRACTION_TYPES_H + +/// @addtogroup imcombine Image Combinations +/// @{ + +/// Mask values for the subtraction mask +typedef enum { + PM_SUBTRACTION_MASK_CLEAR = 0x00, // No masking + PM_SUBTRACTION_MASK_BAD_1 = 0x01, // Image 1 is bad + PM_SUBTRACTION_MASK_BAD_2 = 0x02, // Image 2 is bad + PM_SUBTRACTION_MASK_CONVOLVE_1 = 0x04, // If image 1 is convolved, would be poor or bad + PM_SUBTRACTION_MASK_CONVOLVE_2 = 0x08, // If image 2 is convolved, would be poor or bad + PM_SUBTRACTION_MASK_CONVOLVE_BAD_1 = 0x10, // If image 1 is convolved, would be bad + PM_SUBTRACTION_MASK_CONVOLVE_BAD_2 = 0x20, // If image 2 is convolved, would be bad + PM_SUBTRACTION_MASK_BORDER = 0x40, // Image border + PM_SUBTRACTION_MASK_REJ = 0x80, // Previously tried as a stamp, and rejected +} pmSubtractionMasks; + +/// Type of subtraction kernel +typedef enum { + PM_SUBTRACTION_KERNEL_NONE, ///< Nothing --- an error + PM_SUBTRACTION_KERNEL_POIS, ///< Pan-STARRS Optimal Image Subtraction --- delta functions + PM_SUBTRACTION_KERNEL_ISIS, ///< Traditional kernel --- gaussians modified by polynomials + PM_SUBTRACTION_KERNEL_ISIS_RADIAL, ///< ISIS + higher-order radial Hermitians + PM_SUBTRACTION_KERNEL_HERM, ///< Hermitian polynomial kernels + PM_SUBTRACTION_KERNEL_DECONV_HERM, ///< Deconvolved Hermitian polynomial kernels + PM_SUBTRACTION_KERNEL_SPAM, ///< Summed Pixels for Advanced Matching --- summed delta functions + PM_SUBTRACTION_KERNEL_FRIES, ///< Fibonacci Radius Increases Excellence of Subtraction + PM_SUBTRACTION_KERNEL_GUNK, ///< Grid United with Normal Kernel --- POIS and ISIS hybrid + PM_SUBTRACTION_KERNEL_RINGS, ///< Rings Instead of the Normal Gaussian Subtraction +} pmSubtractionKernelsType; + +/// Modes --- specifies which image to convolve +typedef enum { + PM_SUBTRACTION_MODE_ERR, // Error in the mode + PM_SUBTRACTION_MODE_1, // Convolve image 1 + PM_SUBTRACTION_MODE_2, // Convolve image 2 + PM_SUBTRACTION_MODE_UNSURE, // Not sure yet which image to convolve so try to satisfy both + PM_SUBTRACTION_MODE_DUAL, // Dual convolution +} pmSubtractionMode; + +/// Status of stamp +typedef enum { + PM_SUBTRACTION_STAMP_INIT, ///< Initial state + PM_SUBTRACTION_STAMP_FOUND, ///< Found a suitable source for this stamp + PM_SUBTRACTION_STAMP_CALCULATE, ///< Calculate matrix and vector values for this stamp + PM_SUBTRACTION_STAMP_USED, ///< Use this stamp + PM_SUBTRACTION_STAMP_REJECTED, ///< This stamp has been rejected + PM_SUBTRACTION_STAMP_NONE ///< No stamp in this region +} pmSubtractionStampStatus; + +typedef struct { + double score; + pmSubtractionMode mode; + int spatialOrder; + int nGood; + psVector *fluxes; + psVector *chisq; + psVector *moments; + psVector *stampMask; +} pmSubtractionQuality; + +/// Kernels specification +typedef struct { + pmSubtractionKernelsType type; ///< Type of kernels --- allowing the use of multiple kernels + psString description; ///< Description of the kernel parameters + int xMin, xMax, yMin, yMax; ///< Bounds of image (for normalisation) + long num; ///< Number of kernel components (not including the spatial ones) + psVector *fwhms; ///< requested fwhms of the kernel Gaussians (ISIS, HERM or DECONV_HERM) + psVector *orders; ///< polynomial orders for each Gaussian (ISIS, HERM or DECONV_HERM) + psVector *u, *v; ///< Offset (for POIS) or polynomial order (for ISIS, HERM or DECONV_HERM) + psVector *widths; ///< measured Gaussian FWHMs of Gauss*poly (ISIS, HERM or DECONV_HERM) + psVector *uStop, *vStop; ///< Width of kernel element (SPAM,FRIES only) + psArray *preCalc; ///< Array of images containing pre-calculated kernel (for ISIS, HERM or DECONV_HERM) + float penalty; ///< Penalty for wideness + psVector *penalties1; ///< Penalty for each kernel component + psVector *penalties2; ///< Penalty for each kernel component + bool havePenalties; ///< flag to test if we have already calculated the penalties or not. + int size; ///< The half-size of the kernel + int inner; ///< The size of an inner region + int binning; ///< Binning used for the SPAM kernels + int ringsOrder; ///< + int spatialOrder; ///< The spatial order of the kernels + int bgOrder; ///< The order for the background fitting + pmSubtractionMode mode; ///< Mode for subtraction + psVector *solution1, *solution2; ///< Solution for the PSF matching + psVector *solution1err, *solution2err; ///< error in solution for the PSF matching + // Quality information + float mean, rms; ///< Mean and RMS of chi^2 from stamps + int numStamps; ///< Number of good stamps + float fResSigmaMean; ///< mean fractional stdev of residuals + float fResSigmaStdev; ///< stdev of fractional stdev of residuals + float fResOuterMean; ///< mean fractional positive swing in residuals + float fResOuterStdev; ///< stdev of fractional positive swing in residuals + float fResTotalMean; ///< mean fractional negative swing in residuals + float fResTotalStdev; ///< stdev of fractional negative swing in residuals + psArray *sampleStamps; ///< array of brightest set of stamps for output visualizations +} pmSubtractionKernels; + +// pmSubtractionKernels->preCalc is an array of pmSubtractionKernelPreCalc structures +typedef struct { + psVector *uCoords; // used by RINGS + psVector *vCoords; // used by RINGS + psVector *poly; // used by RINGS + + psVector *xKernel; // used by ISIS, HERM, DECONV_HERM + psVector *yKernel; // used by ISIS, HERM, DECONV_HERM + psKernel *kernel; // used by ISIS, HERM, DECONV_HERM +} pmSubtractionKernelPreCalc; + +/// A list of stamps +typedef struct { + long num; ///< Number of stamps + psArray *stamps; ///< The stamps + psArray *regions; ///< Regions for each stamp + psArray *x, *y; ///< Coordinates for possible stamps (or NULL) + psArray *flux; ///< Fluxes for possible stamps (or NULL) + int footprint; ///< Half-size of stamps + float normFrac; ///< Fraction of flux in window for normalisation window + float normValue; ///< calculated normalization + float normValue2; ///< calculated normalization + psKernel *window1; ///< window function generated from ensemble of stamps (input 1) + psKernel *window2; ///< window function generated from ensemble of stamps (input 2) + psKernel *window; ///< weighting window function (sigma = 1.1 * MAX(fwhm)) + float normWindow1; ///< Size of window for measuring normalisation + float normWindow2; ///< Size of window for measuring normalisation + float sysErr; ///< Systematic error + float skyErr; ///< increase effective readnoise +} pmSubtractionStampList; + +/// A stamp for image subtraction +typedef struct { + float x, y; ///< Position + float flux; ///< Flux + float xNorm, yNorm; ///< Normalised position + psKernel *image1; ///< Reference image postage stamp + psKernel *image2; ///< Input image postage stamp + psKernel *weight; ///< Weight image (1/variance) postage stamp, or NULL + psArray *convolutions1; ///< Convolutions of image 1 for each kernel component, or NULL + psArray *convolutions2; ///< Convolutions of image 2 for each kernel component, or NULL + psImage *matrix; ///< Least-squares matrix, or NULL + psVector *vector; ///< Least-squares vector, or NULL + double norm; ///< Normalisation difference + double normI1; ///< Sum(flux) for image 1 + double normI2; ///< Sum(flux) for image 2 + double normSquare1; ///< Sum(flux^2) for image 1 (used for penalty) + double normSquare2; ///< Sum(flux^2) for image 2 (used for penalty) + pmSubtractionStampStatus status; ///< Status of stamp + psVector *MxxI1; ///< second moments of convolution images + psVector *MyyI1; ///< second moments of convolution images + psVector *MxxI2; ///< second moments of convolution images + psVector *MyyI2; ///< second moments of convolution images + double MxxI1raw; + double MyyI1raw; + double MxxI2raw; + double MyyI2raw; +} pmSubtractionStamp; + +#endif Index: trunk/psModules/src/imcombine/pmSubtractionVisual.c =================================================================== --- trunk/psModules/src/imcombine/pmSubtractionVisual.c (revision 30621) +++ trunk/psModules/src/imcombine/pmSubtractionVisual.c (revision 30622) @@ -16,8 +16,11 @@ #include "pmKapaPlots.h" +#include "pmFPA.h" +#include "pmSubtractionTypes.h" #include "pmSubtraction.h" #include "pmSubtractionStamps.h" #include "pmSubtractionEquation.h" #include "pmSubtractionKernels.h" +#include "pmSubtractionVisual.h" #include "pmVisual.h" @@ -210,4 +213,18 @@ psImageOverlaySection(canvas, im, x0, y0, "="); stampNum++; + + // renormalize the section + float maxValue = 0; + for (int iy = 0; iy < im->numRows; iy++) { + for (int ix = 0; ix < im->numCols; ix++) { + maxValue = PS_MAX(maxValue, im->data.F64[iy][ix]); + } + } + if (maxValue == 0.0) continue; + for (int iy = 0; iy < im->numRows; iy++) { + for (int ix = 0; ix < im->numCols; ix++) { + canvas->data.F64[y0 + iy][x0 + ix] /= maxValue; + } + } } @@ -231,4 +248,109 @@ } +/** Plot the least-squares matrix of each stamp */ +bool pmSubtractionVisualPlotLeastSquaresResid (const pmSubtractionStampList *stamps, psImage *matrixIn, int nUsed) { + + if (!pmVisualTestLevel("ppsub.chisq", 1)) return true; + + if (!plotLeastSquares) return true; + + if (!pmVisualInitWindow (&kapa1, "PPSub:Images")) { + return false; + } + + psImage *matrixNorm = psImageCopy(NULL, matrixIn, PS_TYPE_F64); + { + // renormalize the matrix + float maxValue = 0; + for (int iy = 0; iy < matrixNorm->numRows; iy++) { + for (int ix = 0; ix < matrixNorm->numCols; ix++) { + maxValue = PS_MAX(maxValue, matrixNorm->data.F64[iy][ix]); + } + } + for (int iy = 0; iy < matrixNorm->numRows; iy++) { + for (int ix = 0; ix < matrixNorm->numCols; ix++) { + matrixNorm->data.F64[iy][ix] /= maxValue; + } + } + } + + // Find the stamp size + int imageMax = -1; + int numStamps = 0; + psElemType type = PS_TYPE_F64; + + for(int i = 0; i < stamps->num; i++) { + pmSubtractionStamp *stamp = stamps->stamps->data[i]; + if (stamp == NULL) continue; + + psImage *im = stamp->matrix; + if (im == NULL) continue; + + imageMax = PS_MAX(imageMax, im->numCols); + imageMax = PS_MAX(imageMax, im->numRows); + numStamps++; + type = im->type.type; + } + if (imageMax == -1) return false; + + int border = 15; + imageMax += border; + int tileRowCount = (int) ceil(sqrt(numStamps)); + int canvasX = tileRowCount * (imageMax); + int canvasY = tileRowCount * (imageMax); + psImage *canvas = psImageAlloc (canvasX, canvasY, type); + psImageInit (canvas, NAN); + + // overlay the images + int stampNum = 0; + int stampListNum = 0; + while (stampNum < numStamps) { + int x0 = (imageMax) * (stampNum % tileRowCount); + int y0 = (imageMax) * (stampNum / tileRowCount); + + pmSubtractionStamp *stamp = stamps->stamps->data[stampListNum++]; + if (stamp == NULL) continue; + + psImage *im = stamp->matrix; + if (im == NULL) continue; + + stampNum++; + + if (stamp->status != PM_SUBTRACTION_STAMP_USED) continue; + + // renormalize the section + float maxValue = 0; + for (int iy = 0; iy < im->numRows; iy++) { + for (int ix = 0; ix < im->numCols; ix++) { + maxValue = PS_MAX(maxValue, im->data.F64[iy][ix]); + } + } + if (maxValue == 0.0) continue; + for (int iy = 0; iy < im->numRows; iy++) { + for (int ix = 0; ix < im->numCols; ix++) { + canvas->data.F64[y0 + iy][x0 + ix] = (im->data.F64[iy][ix] / maxValue - matrixNorm->data.F64[iy][ix]) / matrixNorm->data.F64[iy][ix]; + } + } + } + + psImage *canvas32 = pmVisualImageToFloat(canvas); + pmVisualRangeImage(kapa2, canvas32, "Least_Squares", 0, -100.0, 100.0); + + if (0) { + static int count = 0; + char filename[64]; + sprintf (filename, "chisq.%02d.fits", count); + count ++; + psFits *fits = psFitsOpen (filename, "w"); + psFitsWriteImage (fits, NULL, canvas32, 0, NULL); + psFitsClose (fits); + } + + pmVisualAskUser(&plotLeastSquares); + psFree(canvas); + psFree(canvas32); + return true; +} + bool pmSubtractionVisualShowSubtraction(psImage *image, psImage *ref, psImage *sub) { @@ -304,5 +426,6 @@ } - // XXX clear the overlay(s) (red at least!) + // clear the overlay (red at least!) + KiiEraseOverlay (kapa2, "red"); // get the kernel sizes @@ -337,4 +460,5 @@ int nKernels = 0; + // paste in the kernel images, scaled by sum2 if (maxStamp->convolutions1) { // output image is a grid of NXsub by NYsub sub-images @@ -359,14 +483,16 @@ int yPix = ySub * (2*footprint + 1 + 3) + footprint; - double sum = 0.0; double sum2 = 0.0; for (int y = -footprint; y <= footprint; y++) { for (int x = -footprint; x <= footprint; x++) { - output->data.F32[y + yPix][x + xPix] = kernel->kernel[y][x]; - sum += kernel->kernel[y][x]; sum2 += PS_SQR(kernel->kernel[y][x]); } } - // fprintf (stderr, "kernel %d, sum %f, sum2: %e\n", i, sum, sum2); + float scale = sqrt(sum2) / PS_SQR(2*footprint + 1); + for (int y = -footprint; y <= footprint; y++) { + for (int x = -footprint; x <= footprint; x++) { + output->data.F32[y + yPix][x + xPix] = kernel->kernel[y][x] / scale; + } + } } pmVisualScaleImage(kapa2, output, "Image", 0, true); @@ -401,14 +527,16 @@ int yPix = ySub * (2*footprint + 1 + 3) + footprint; - double sum = 0.0; double sum2 = 0.0; for (int y = -footprint; y <= footprint; y++) { for (int x = -footprint; x <= footprint; x++) { - output->data.F32[y + yPix][x + xPix] = kernel->kernel[y][x]; - sum += kernel->kernel[y][x]; sum2 += PS_SQR(kernel->kernel[y][x]); } } - // fprintf (stderr, "kernel %d, sum %f, sum2: %e\n", i, sum, sum2); + float scale = sqrt(sum2) / PS_SQR(2*footprint + 1); + for (int y = -footprint; y <= footprint; y++) { + for (int x = -footprint; x <= footprint; x++) { + output->data.F32[y + yPix][x + xPix] = kernel->kernel[y][x] / scale; + } + } } pmVisualScaleImage(kapa2, output, "Image", 1, true); @@ -468,122 +596,4 @@ KiiLoadOverlay (kapa2, overlay, Noverlay, "red"); FREE (overlay); - return true; -} - -static int footprint = 0; -static int NX = 0; -static int NY = 0; -static psImage *sourceImage = NULL; -static psImage *targetImage = NULL; -static psImage *residualImage = NULL; -static psImage *fresidualImage = NULL; -static psImage *differenceImage = NULL; -static psImage *convolutionImage = NULL; - -bool pmSubtractionVisualShowFitInit(pmSubtractionStampList *stamps) { - - if (!pmVisualTestLevel("ppsub.fit", 1)) return true; - - // generate 4 storage images large enough to hold the N stamps: - - footprint = stamps->footprint; - - float NXf = sqrt(stamps->num); - NX = (int) NXf == NXf ? NXf : NXf + 1.0; - - float NYf = stamps->num / NX; - NY = (int) NYf == NY ? NYf : NYf + 1.0; - - int NXpix = (2*footprint + 1) * NX; - NXpix += (NX > 1) ? 3 * NX : 0; - - int NYpix = (2*footprint + 1) * NY; - NYpix += (NY > 1) ? 3 * NY : 0; - - sourceImage = psImageAlloc (NXpix, NYpix, PS_TYPE_F32); - targetImage = psImageAlloc (NXpix, NYpix, PS_TYPE_F32); - residualImage = psImageAlloc (NXpix, NYpix, PS_TYPE_F32); - fresidualImage = psImageAlloc (NXpix, NYpix, PS_TYPE_F32); - differenceImage = psImageAlloc (NXpix, NYpix, PS_TYPE_F32); - convolutionImage = psImageAlloc (NXpix, NYpix, PS_TYPE_F32); - - psImageInit (sourceImage, 0.0); - psImageInit (targetImage, 0.0); - psImageInit (residualImage, 0.0); - psImageInit (fresidualImage, 0.0); - psImageInit (differenceImage, 0.0); - psImageInit (convolutionImage, 0.0); - - return true; -} - -bool pmSubtractionVisualShowFitAddStamp(psKernel *target, psKernel *source, psKernel *convolution, double background, double norm, int index) { - - if (!pmVisualTestLevel("ppsub.stamp", 1)) return true; - - double sum; - - int NXoff = index % NX; - int NYoff = index / NX; - - int NXpix = NXoff * (2*footprint + 1 + 3) + footprint; - int NYpix = NYoff * (2*footprint + 1 + 3) + footprint; - - // insert the (target) kernel into the (target) image: - sum = 0.0; - for (int y = -footprint; y <= footprint; y++) { - for (int x = -footprint; x <= footprint; x++) { - targetImage->data.F32[y + NYpix][x + NXpix] = target->kernel[y][x]; - sum += targetImage->data.F32[y + NYpix][x + NXpix]; - } - } - targetImage->data.F32[footprint + 1 + NYpix][NXpix] = sum; - - // insert the (source) kernel into the (source) image: - sum = 0.0; - for (int y = -footprint; y <= footprint; y++) { - for (int x = -footprint; x <= footprint; x++) { - sourceImage->data.F32[y + NYpix][x + NXpix] = source->kernel[y][x]; - sum += sourceImage->data.F32[y + NYpix][x + NXpix]; - } - } - sourceImage->data.F32[footprint + 1 + NYpix][NXpix] = sum; - - // insert the (convolution) kernel into the (convolution) image: - sum = 0.0; - for (int y = -footprint; y <= footprint; y++) { - for (int x = -footprint; x <= footprint; x++) { - convolutionImage->data.F32[y + NYpix][x + NXpix] = convolution->kernel[y][x]; - sum += convolutionImage->data.F32[y + NYpix][x + NXpix]; - } - } - convolutionImage->data.F32[footprint + 1 + NYpix][NXpix] = sum; - - // insert the (difference) kernel into the (difference) image: - sum = 0.0; - for (int y = -footprint; y <= footprint; y++) { - for (int x = -footprint; x <= footprint; x++) { - differenceImage->data.F32[y + NYpix][x + NXpix] = target->kernel[y][x] - background - source->kernel[y][x] * norm; - sum += differenceImage->data.F32[y + NYpix][x + NXpix]; - } - } - differenceImage->data.F32[footprint + 1 + NYpix][NXpix] = sum; - - // insert the (residual) kernel into the (residual) image: - sum = 0.0; - for (int y = -footprint; y <= footprint; y++) { - for (int x = -footprint; x <= footprint; x++) { - residualImage->data.F32[y + NYpix][x + NXpix] = target->kernel[y][x] - background - source->kernel[y][x] * norm - convolution->kernel[y][x]; - sum += residualImage->data.F32[y + NYpix][x + NXpix]; - } - } - residualImage->data.F32[footprint + 1 + NYpix][NXpix] = sum; - - // insert the (fresidual) kernel into the (fresidual) image: - for (int y = -footprint; y <= footprint; y++) { - for (int x = -footprint; x <= footprint; x++) { - fresidualImage->data.F32[y + NYpix][x + NXpix] = residualImage->data.F32[y + NYpix][x + NXpix] / sqrt(PS_MAX(target->kernel[y][x], 100.0)); - } - } return true; } @@ -611,15 +621,15 @@ } -bool pmSubtractionVisualShowFit(double norm) { - - // for testing, dump the residual image and exit - if (0) { - psMetadata *header = psMetadataAlloc(); - psMetadataAddF32 (header, PS_LIST_TAIL, "NORM", 0, "Normalization", norm); - psFits *fits = psFitsOpen("resid.fits", "w"); - psFitsWriteImage(fits, header, residualImage, 0, NULL); - psFitsClose(fits); - // exit (0); - } +static int footprint = 0; +static int NX = 0; +static int NY = 0; +static psImage *sourceImage = NULL; +static psImage *targetImage = NULL; +static psImage *residualImage = NULL; +static psImage *fresidualImage = NULL; +static psImage *differenceImage = NULL; +static psImage *convolutionImage = NULL; + +bool pmSubtractionVisualShowFit(pmSubtractionStampList *stamps, pmSubtractionKernels *kernels) { if (!pmVisualTestLevel("ppsub.fit", 1)) return true; @@ -627,4 +637,187 @@ if (!pmVisualInitWindow(&kapa1, "ppSub:Images")) return false; if (!pmVisualInitWindow(&kapa2, "ppSub:Misc")) return false; + + // set up holding images for the visualization + pmSubtractionVisualShowFitInit (stamps); + + int numKernels = kernels->num; // Number of kernels + + psImage *polyValues = NULL; // Polynomial values + psKernel *residual = psKernelAlloc(-stamps->footprint, stamps->footprint, -stamps->footprint, stamps->footprint); // Residual image + + double norm = p_pmSubtractionSolutionNorm(kernels); // Normalisation + + for (int i = 0; i < stamps->num; i++) { + pmSubtractionStamp *stamp = stamps->stamps->data[i]; // The stamp of interest + if (stamp->status != PM_SUBTRACTION_STAMP_USED) { continue; } + + // Calculate coefficients of the kernel basis functions + polyValues = p_pmSubtractionPolynomial(polyValues, kernels->spatialOrder, stamp->xNorm, stamp->yNorm); + double background = p_pmSubtractionSolutionBackground(kernels, polyValues); // Difference in background + + psImageInit(residual->image, 0.0); + + if (kernels->mode != PM_SUBTRACTION_MODE_DUAL) { + psKernel *target; // Target postage stamp + psKernel *source; // Source postage stamp + psArray *convolutions; // Convolution postage stamps for each kernel basis function + switch (kernels->mode) { + case PM_SUBTRACTION_MODE_1: + target = stamp->image2; + source = stamp->image1; + convolutions = stamp->convolutions1; + break; + case PM_SUBTRACTION_MODE_2: + target = stamp->image1; + source = stamp->image2; + convolutions = stamp->convolutions2; + break; + default: + psAbort("Unsupported subtraction mode: %x", kernels->mode); + } + + for (int j = 0; j < numKernels; j++) { + psKernel *convolution = convolutions->data[j]; // Convolution + double coefficient = p_pmSubtractionSolutionCoeff(kernels, polyValues, j, false); // Coefficient + for (int y = - footprint; y <= footprint; y++) { + for (int x = - footprint; x <= footprint; x++) { + residual->kernel[y][x] += convolution->kernel[y][x] * coefficient; + } + } + } + // visualize the target, source, convolution and residual + pmSubtractionVisualShowFitAddStamp (target, source, residual, background, norm, i); + } else { + // Dual convolution + psArray *convolutions1 = stamp->convolutions1; // Convolutions of the first image + psArray *convolutions2 = stamp->convolutions2; // Convolutions of the second image + psKernel *image1 = stamp->image1; // The first image + psKernel *image2 = stamp->image2; // The second image + + for (int j = 0; j < numKernels; j++) { + psKernel *conv1 = convolutions1->data[j]; // Convolution of first image + psKernel *conv2 = convolutions2->data[j]; // Convolution of second image + double coeff1 = p_pmSubtractionSolutionCoeff(kernels, polyValues, j, false); // Coefficient 1 + double coeff2 = p_pmSubtractionSolutionCoeff(kernels, polyValues, j, true); // Coefficient 2 + + for (int y = - footprint; y <= footprint; y++) { + for (int x = - footprint; x <= footprint; x++) { + residual->kernel[y][x] += conv2->kernel[y][x] * coeff2 + conv1->kernel[y][x] * coeff1; + } + } + } + // visualize the target, source, convolution and residual + pmSubtractionVisualShowFitAddStamp (image2, image1, residual, background, norm, i); + } + } + pmSubtractionVisualShowFitImage(norm); + + return true; +} + +// generate 4 storage images large enough to hold the stamps: +bool pmSubtractionVisualShowFitInit(pmSubtractionStampList *stamps) { + + footprint = stamps->footprint; + + float NXf = sqrt(stamps->num); + NX = (int) NXf == NXf ? NXf : NXf + 1.0; + + float NYf = stamps->num / NX; + NY = (int) NYf == NY ? NYf : NYf + 1.0; + + int NXpix = (2*footprint + 1) * NX; + NXpix += (NX > 1) ? 3 * NX : 0; + + int NYpix = (2*footprint + 1) * NY; + NYpix += (NY > 1) ? 3 * NY : 0; + + sourceImage = psImageAlloc (NXpix, NYpix, PS_TYPE_F32); + targetImage = psImageAlloc (NXpix, NYpix, PS_TYPE_F32); + residualImage = psImageAlloc (NXpix, NYpix, PS_TYPE_F32); + fresidualImage = psImageAlloc (NXpix, NYpix, PS_TYPE_F32); + differenceImage = psImageAlloc (NXpix, NYpix, PS_TYPE_F32); + convolutionImage = psImageAlloc (NXpix, NYpix, PS_TYPE_F32); + + psImageInit (sourceImage, 0.0); + psImageInit (targetImage, 0.0); + psImageInit (residualImage, 0.0); + psImageInit (fresidualImage, 0.0); + psImageInit (differenceImage, 0.0); + psImageInit (convolutionImage, 0.0); + + return true; +} + +bool pmSubtractionVisualShowFitAddStamp(psKernel *target, psKernel *source, psKernel *convolution, double background, double norm, int index) { + + double sum; + + int NXoff = index % NX; + int NYoff = index / NX; + + int NXpix = NXoff * (2*footprint + 1 + 3) + footprint; + int NYpix = NYoff * (2*footprint + 1 + 3) + footprint; + + // insert the (target) kernel into the (target) image: + sum = 0.0; + for (int y = -footprint; y <= footprint; y++) { + for (int x = -footprint; x <= footprint; x++) { + targetImage->data.F32[y + NYpix][x + NXpix] = target->kernel[y][x]; + sum += targetImage->data.F32[y + NYpix][x + NXpix]; + } + } + targetImage->data.F32[footprint + 1 + NYpix][NXpix] = sum; + + // insert the (source) kernel into the (source) image: + sum = 0.0; + for (int y = -footprint; y <= footprint; y++) { + for (int x = -footprint; x <= footprint; x++) { + sourceImage->data.F32[y + NYpix][x + NXpix] = source->kernel[y][x]; + sum += sourceImage->data.F32[y + NYpix][x + NXpix]; + } + } + sourceImage->data.F32[footprint + 1 + NYpix][NXpix] = sum; + + // insert the (convolution) kernel into the (convolution) image: + sum = 0.0; + for (int y = -footprint; y <= footprint; y++) { + for (int x = -footprint; x <= footprint; x++) { + convolutionImage->data.F32[y + NYpix][x + NXpix] = convolution->kernel[y][x]; + sum += convolutionImage->data.F32[y + NYpix][x + NXpix]; + } + } + convolutionImage->data.F32[footprint + 1 + NYpix][NXpix] = sum; + + // insert the (difference) kernel into the (difference) image: + sum = 0.0; + for (int y = -footprint; y <= footprint; y++) { + for (int x = -footprint; x <= footprint; x++) { + differenceImage->data.F32[y + NYpix][x + NXpix] = target->kernel[y][x] - background - source->kernel[y][x] * norm; + sum += differenceImage->data.F32[y + NYpix][x + NXpix]; + } + } + differenceImage->data.F32[footprint + 1 + NYpix][NXpix] = sum; + + // insert the (residual) kernel into the (residual) image: + sum = 0.0; + for (int y = -footprint; y <= footprint; y++) { + for (int x = -footprint; x <= footprint; x++) { + residualImage->data.F32[y + NYpix][x + NXpix] = target->kernel[y][x] - background - source->kernel[y][x] * norm - convolution->kernel[y][x]; + sum += residualImage->data.F32[y + NYpix][x + NXpix]; + } + } + residualImage->data.F32[footprint + 1 + NYpix][NXpix] = sum; + + // insert the (fresidual) kernel into the (fresidual) image: + for (int y = -footprint; y <= footprint; y++) { + for (int x = -footprint; x <= footprint; x++) { + fresidualImage->data.F32[y + NYpix][x + NXpix] = residualImage->data.F32[y + NYpix][x + NXpix] / sqrt(PS_MAX(target->kernel[y][x], 100.0)); + } + } + return true; +} + +bool pmSubtractionVisualShowFitImage(double norm) { KiiEraseOverlay (kapa1, "red"); @@ -673,4 +866,5 @@ psVector *x = psVectorAllocEmpty (kernels->num, PS_TYPE_F32); psVector *y = psVectorAllocEmpty (kernels->num, PS_TYPE_F32); + psVector *dy = psVectorAllocEmpty (kernels->num, PS_TYPE_F32); graphdata.xmin = -1.0; @@ -685,8 +879,9 @@ x->data.F32[i] = i; y->data.F32[i] = p_pmSubtractionSolutionCoeff(kernels, polyValues, i, false); + dy->data.F32[i] = kernels->solution1err->data.F64[i]; graphdata.ymin = PS_MIN(graphdata.ymin, y->data.F32[i]); graphdata.ymax = PS_MAX(graphdata.ymax, y->data.F32[i]); } - x->n = y->n = kernels->num; + x->n = y->n = dy->n = kernels->num; float range; @@ -709,12 +904,128 @@ graphdata.size = 0.5; graphdata.style = 2; + graphdata.etype |= 0x01; KapaPrepPlot (kapa3, x->n, &graphdata); KapaPlotVector (kapa3, x->n, x->data.F32, "x"); KapaPlotVector (kapa3, x->n, y->data.F32, "y"); + KapaPlotVector (kapa3, x->n, dy->data.F32, "dym"); + KapaPlotVector (kapa3, x->n, dy->data.F32, "dyp"); psFree (x); psFree (y); + psFree (dy); psFree (polyValues); + + pmVisualAskUser(NULL); + return true; +} + +// plot log(flux) vs log(chisq), log(flux) vs log(moments), log(chisq) vs log(moments) +bool pmSubtractionVisualPlotChisqAndMoments(psVector *fluxes, psVector *chisq, psVector *moments) { + + Graphdata graphdata; + KapaSection section; + + if (!pmVisualTestLevel("ppsub.fit", 1)) return true; + + if (!pmVisualInitWindow(&kapa3, "ppSub:plots")) return false; + + KapaClearSections (kapa3); + KapaInitGraph (&graphdata); + KiiResize(kapa3, 1500, 500); + + psVector *lchi = psVectorAlloc (fluxes->n, PS_TYPE_F32); + psVector *lflx = psVectorAlloc (fluxes->n, PS_TYPE_F32); + psVector *lMxx = psVectorAlloc (fluxes->n, PS_TYPE_F32); + + // construct the plot vectors + for (int i = 0; i < fluxes->n; i++) { + lchi->data.F32[i] = log10(chisq->data.F32[i]); + lflx->data.F32[i] = log10(fluxes->data.F32[i]); + lMxx->data.F32[i] = log10(moments->data.F32[i]); + } + + section.bg = KapaColorByName ("none"); // XXX probably should be 'none' + + // section 1: lflux vs lchi + section.dx = 0.33; + section.dy = 1.00; + section.x = 0.00; + section.y = 0.00; + section.name = psStringCopy ("flux.v.chi"); + KapaSetSection (kapa3, §ion); + psFree (section.name); + + graphdata.color = KapaColorByName ("black"); + pmVisualScaleGraphdata(&graphdata, lflx, lchi, false); + KapaSetLimits (kapa3, &graphdata); + + KapaSetFont (kapa3, "helvetica", 14); + KapaBox (kapa3, &graphdata); + KapaSendLabel (kapa3, "log(flux)", KAPA_LABEL_XM); + KapaSendLabel (kapa3, "log(chisq)", KAPA_LABEL_YM); + + graphdata.color = KapaColorByName ("black"); + graphdata.ptype = 2; + graphdata.size = 0.5; + graphdata.style = 2; + + KapaPrepPlot (kapa3, lflx->n, &graphdata); + KapaPlotVector (kapa3, lflx->n, lflx->data.F32, "x"); + KapaPlotVector (kapa3, lflx->n, lchi->data.F32, "y"); + + // section 2: lflux vs lMxx + section.dx = 0.33; + section.dy = 1.00; + section.x = 0.33; + section.y = 0.00; + section.name = psStringCopy ("flux.v.mom"); + KapaSetSection (kapa3, §ion); + psFree (section.name); + + graphdata.color = KapaColorByName ("black"); + pmVisualScaleGraphdata(&graphdata, lflx, lMxx, false); + KapaSetLimits (kapa3, &graphdata); + + KapaSetFont (kapa3, "helvetica", 14); + KapaBox (kapa3, &graphdata); + KapaSendLabel (kapa3, "log(flux)", KAPA_LABEL_XM); + KapaSendLabel (kapa3, "log(moments)", KAPA_LABEL_YM); + + graphdata.color = KapaColorByName ("black"); + graphdata.ptype = 2; + graphdata.size = 0.5; + graphdata.style = 2; + + KapaPrepPlot (kapa3, lflx->n, &graphdata); + KapaPlotVector (kapa3, lflx->n, lflx->data.F32, "x"); + KapaPlotVector (kapa3, lflx->n, lMxx->data.F32, "y"); + + // section 1: lflux vs lchi + section.dx = 0.33; + section.dy = 1.00; + section.x = 0.66; + section.y = 0.00; + section.name = psStringCopy ("chi.v.mom"); + KapaSetSection (kapa3, §ion); + psFree (section.name); + + graphdata.color = KapaColorByName ("black"); + pmVisualScaleGraphdata(&graphdata, lchi, lMxx, false); + KapaSetLimits (kapa3, &graphdata); + + KapaSetFont (kapa3, "helvetica", 14); + KapaBox (kapa3, &graphdata); + KapaSendLabel (kapa3, "log(chisq)", KAPA_LABEL_XM); + KapaSendLabel (kapa3, "log(moments)", KAPA_LABEL_YM); + + graphdata.color = KapaColorByName ("black"); + graphdata.ptype = 2; + graphdata.size = 0.5; + graphdata.style = 2; + + KapaPrepPlot (kapa3, lflx->n, &graphdata); + KapaPlotVector (kapa3, lflx->n, lchi->data.F32, "x"); + KapaPlotVector (kapa3, lflx->n, lMxx->data.F32, "y"); pmVisualAskUser(NULL); Index: trunk/psModules/src/imcombine/pmSubtractionVisual.h =================================================================== --- trunk/psModules/src/imcombine/pmSubtractionVisual.h (revision 30621) +++ trunk/psModules/src/imcombine/pmSubtractionVisual.h (revision 30622) @@ -6,11 +6,16 @@ bool pmSubtractionVisualPlotStamps(pmSubtractionStampList *stamps, pmReadout *ro); bool pmSubtractionVisualPlotLeastSquares(pmSubtractionStampList *stamps); +bool pmSubtractionVisualPlotLeastSquaresResid (const pmSubtractionStampList *stamps, psImage *matrixIn, int nUsed); bool pmSubtractionVisualShowSubtraction(psImage *image, psImage *ref, psImage *sub); + +bool pmSubtractionVisualShowFit(pmSubtractionStampList *stamps, pmSubtractionKernels *kernels); bool pmSubtractionVisualShowFitInit(pmSubtractionStampList *stamps); bool pmSubtractionVisualShowFitAddStamp(psKernel *target, psKernel *source, psKernel *convolution, double background, double norm, int index); -bool pmSubtractionVisualShowFit(double norm); +bool pmSubtractionVisualShowFitImage(double norm); + bool pmSubtractionVisualPlotFit(const pmSubtractionKernels *kernels); bool pmSubtractionVisualShowKernels(pmSubtractionKernels *kernels); bool pmSubtractionVisualShowBasis(pmSubtractionStampList *stamps); +bool pmSubtractionVisualPlotChisqAndMoments(psVector *fluxes, psVector *chisq, psVector *moments); #endif