Index: trunk/psModules/src/imcombine/pmSubtractionMatch.c =================================================================== --- trunk/psModules/src/imcombine/pmSubtractionMatch.c (revision 15329) +++ trunk/psModules/src/imcombine/pmSubtractionMatch.c (revision 15443) @@ -50,6 +50,6 @@ static bool getStamps(pmSubtractionStampList **stamps, // Stamps to read - const pmReadout *reference, // Reference readout - const pmReadout *input, // Input readout, or NULL to generate fake stamps + const pmReadout *ro1, // Readout 1 + const pmReadout *ro2, // Readout 2 const psImage *subMask, // Mask for subtraction, or NULL psImage *weight, // Weight map @@ -57,11 +57,11 @@ float threshold, // Threshold for stamp finding float stampSpacing, // Spacing between stamps - float targetWidth,// Target width for fake stamps int size, // Kernel half-size - int footprint // Convolution footprint for stamps + int footprint, // Convolution footprint for stamps + pmSubtractionMode mode // Mode for subtraction ) { psTrace("psModules.imcombine", 3, "Finding stamps...\n"); - *stamps = pmSubtractionStampsFind(*stamps, reference->image, subMask, region, threshold, stampSpacing); + *stamps = pmSubtractionStampsFind(*stamps, ro1->image, subMask, region, threshold, stampSpacing, mode); if (!*stamps) { psError(PS_ERR_UNKNOWN, false, "Unable to find stamps."); @@ -71,14 +71,6 @@ memCheck(" find stamps"); - if (!input && !pmSubtractionStampsGenerate(*stamps, targetWidth, footprint, size)) { - psError(PS_ERR_UNKNOWN, false, "Unable to generate target stamps."); - return false; - } - - memCheck(" generate stamps"); - psTrace("psModules.imcombine", 3, "Extracting stamps...\n"); - if (!pmSubtractionStampsExtract(*stamps, reference->image, input ? input->image : NULL, - weight, footprint, size)) { + if (!pmSubtractionStampsExtract(*stamps, ro1->image, ro2 ? ro2->image : NULL, weight, footprint, size)) { psError(PS_ERR_UNKNOWN, false, "Unable to extract stamps."); return false; @@ -94,40 +86,40 @@ -bool pmSubtractionMatch(pmReadout *convolved, const pmReadout *reference, const pmReadout *input, +bool pmSubtractionMatch(pmReadout *convolved, const pmReadout *ro1, const pmReadout *ro2, int footprint, float regionSize, float stampSpacing, float threshold, - const psArray *sources, const char *stampsName, float targetWidth, + const psArray *sources, const char *stampsName, pmSubtractionKernelsType type, int size, int spatialOrder, const psVector *isisWidths, const psVector *isisOrders, int inner, int ringsOrder, int binning, bool optimum, const psVector *optFWHMs, int optOrder, float optThreshold, int iter, float rej, psMaskType maskBad, - psMaskType maskBlank, float badFrac) + psMaskType maskBlank, float badFrac, pmSubtractionMode mode) { PS_ASSERT_PTR_NON_NULL(convolved, false); - PS_ASSERT_PTR_NON_NULL(reference, false); - PS_ASSERT_IMAGE_NON_NULL(reference->image, false); - PS_ASSERT_IMAGE_TYPE(reference->image, PS_TYPE_F32, false); - if (reference->mask) { - PS_ASSERT_IMAGE_NON_NULL(reference->mask, false); - PS_ASSERT_IMAGE_TYPE(reference->mask, PS_TYPE_MASK, false); - PS_ASSERT_IMAGES_SIZE_EQUAL(reference->mask, reference->image, false); - } - if (reference->weight) { - PS_ASSERT_IMAGE_NON_NULL(reference->weight, false); - PS_ASSERT_IMAGE_TYPE(reference->weight, PS_TYPE_F32, false); - PS_ASSERT_IMAGES_SIZE_EQUAL(reference->weight, reference->image, false); - } - if (input) { - PS_ASSERT_IMAGE_NON_NULL(input->image, false); - PS_ASSERT_IMAGE_TYPE(input->image, PS_TYPE_F32, false); - PS_ASSERT_IMAGES_SIZE_EQUAL(input->image, reference->image, false); - if (input->mask) { - PS_ASSERT_IMAGE_NON_NULL(input->mask, false); - PS_ASSERT_IMAGE_TYPE(input->mask, PS_TYPE_MASK, false); - PS_ASSERT_IMAGES_SIZE_EQUAL(input->mask, reference->image, false); - } - if (input->weight) { - PS_ASSERT_IMAGE_NON_NULL(input->weight, false); - PS_ASSERT_IMAGE_TYPE(input->weight, PS_TYPE_F32, false); - PS_ASSERT_IMAGES_SIZE_EQUAL(input->weight, reference->image, false); + PS_ASSERT_PTR_NON_NULL(ro1, false); + PS_ASSERT_IMAGE_NON_NULL(ro1->image, false); + PS_ASSERT_IMAGE_TYPE(ro1->image, PS_TYPE_F32, false); + if (ro1->mask) { + PS_ASSERT_IMAGE_NON_NULL(ro1->mask, false); + PS_ASSERT_IMAGE_TYPE(ro1->mask, PS_TYPE_MASK, false); + PS_ASSERT_IMAGES_SIZE_EQUAL(ro1->mask, ro1->image, false); + } + if (ro1->weight) { + PS_ASSERT_IMAGE_NON_NULL(ro1->weight, false); + PS_ASSERT_IMAGE_TYPE(ro1->weight, PS_TYPE_F32, false); + PS_ASSERT_IMAGES_SIZE_EQUAL(ro1->weight, ro1->image, false); + } + if (ro2) { + PS_ASSERT_IMAGE_NON_NULL(ro2->image, false); + PS_ASSERT_IMAGE_TYPE(ro2->image, PS_TYPE_F32, false); + PS_ASSERT_IMAGES_SIZE_EQUAL(ro2->image, ro1->image, false); + if (ro2->mask) { + PS_ASSERT_IMAGE_NON_NULL(ro2->mask, false); + PS_ASSERT_IMAGE_TYPE(ro2->mask, PS_TYPE_MASK, false); + PS_ASSERT_IMAGES_SIZE_EQUAL(ro2->mask, ro1->image, false); + } + if (ro2->weight) { + PS_ASSERT_IMAGE_NON_NULL(ro2->weight, false); + PS_ASSERT_IMAGE_TYPE(ro2->weight, PS_TYPE_F32, false); + PS_ASSERT_IMAGES_SIZE_EQUAL(ro2->weight, ro1->image, false); } } else if (!stampsName && !sources) { @@ -144,5 +136,4 @@ } // stampsName may be anything - // targetWidth can be just about anything (except maybe negative, but it can be NAN) // We'll check kernel type when we allocate the kernels PS_ASSERT_INT_POSITIVE(size, false); @@ -196,20 +187,18 @@ } - psImage *inImage = NULL, *inMask = NULL; // Input image, mask, weight - if (input) { - inImage = input->image; - inMask = input->mask; - } - // Where does our weight map come from? psImage *weight = NULL; // Weight image to use - if (input && input->weight) { - weight = input->weight; - } else if (reference->weight) { - weight = reference->weight; - } else if (input) { - weight = input->image; + if (ro1->weight && ro2 && ro2->weight) { + weight = (psImage*)psBinaryOp(NULL, ro1->weight, "+", ro2->weight); + } else if (ro1->weight) { + weight = psMemIncrRefCounter(ro1->weight); + } else if (ro2) { + if (ro2->weight) { + weight = psMemIncrRefCounter(ro2->weight); + } else { + weight = (psImage*)psBinaryOp(NULL, ro1->image, "+", ro2->image); + } } else { - weight = reference->image; + weight = psMemIncrRefCounter(ro1->image); } @@ -221,12 +210,13 @@ psVector *solution = NULL; // Solution to match PSF pmSubtractionKernels *kernels = NULL; // Kernel basis functions - - int numCols = reference->image->numCols, numRows = reference->image->numRows; // Image dimensions + int numCols = ro1->image->numCols, numRows = ro1->image->numRows; // Image dimensions memCheck("start"); - subMask = pmSubtractionMask(reference->mask, inMask, maskBad, size, footprint, badFrac, useFFT); + subMask = pmSubtractionMask(ro1->mask, ro2 ? ro2->mask : NULL, maskBad, size, footprint, + badFrac, useFFT); if (!subMask) { psError(PS_ERR_UNKNOWN, false, "Unable to generate subtraction mask."); + psFree(weight); return false; } @@ -260,39 +250,41 @@ if (sources) { - stamps = pmSubtractionStampsSetFromSources(sources, subMask, region, stampSpacing, - input ? 0 : 2 * footprint); + stamps = pmSubtractionStampsSetFromSources(sources, subMask, region, stampSpacing, mode); } else if (stampsName && strlen(stampsName) > 0) { - // Read stamps from file - psTrace("psModules.imcombine", 3, "Reading stamps from %s...\n", stampsName); - const char *stampFormat = input ? "%f %f" : "%f %f %f"; // Format for reading stamp file - psArray *stampsData = stampsData = psVectorsReadFromFile(stampsName, stampFormat); - psVector *xStamp = NULL, *yStamp = NULL, *fluxStamp = NULL; // Stamp positions and fluxes - if (!stampsData) { - psError(PS_ERR_IO, false, "Unable to read stamps file %s", stampsName); - goto ERROR; - } - xStamp = stampsData->data[0]; - yStamp = stampsData->data[1]; - if (!input) { - fluxStamp = stampsData->data[2]; - } - - // Correct for IRAF/FITS (unit-offset) positions to C (zero-offset) positions - psBinaryOp(xStamp, xStamp, "-", psScalarAlloc(1.0, PS_TYPE_F32)); - psBinaryOp(yStamp, yStamp, "-", psScalarAlloc(1.0, PS_TYPE_F32)); - - stamps = pmSubtractionStampsSet(xStamp, yStamp, fluxStamp, reference->image, subMask, - region, stampSpacing, input ? 0 : footprint + size); - psFree(stampsData); + stamps = pmSubtractionStampsSetFromFile(stampsName, subMask, region, stampSpacing, mode); + } + + // We get the stamps here; we will also attempt to get stamps at the first iteration, but it + // doesn't matter. + if (!getStamps(&stamps, ro1, ro2, subMask, weight, NULL, threshold, stampSpacing, + size, footprint, mode)) { + goto MATCH_ERROR; + } + + if (mode == PM_SUBTRACTION_MODE_UNSURE || mode == PM_SUBTRACTION_MODE_TARGET) { + pmSubtractionMode newMode = pmSubtractionOrder(stamps, footprint); // Subtraction mode + 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: + if (mode == PM_SUBTRACTION_MODE_TARGET) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, + "Input PSF is larger than target PSF --- can't match image."); + goto MATCH_ERROR; + } + psLogMsg("psModules.imcombine", PS_LOG_INFO, "Convolving image 2 to match image 1."); + break; + default: + psError(PS_ERR_UNKNOWN, false, "Unable to determine subtraction order."); + goto MATCH_ERROR; + } + mode = newMode; } // Define kernel basis functions if (optimum && (type == PM_SUBTRACTION_KERNEL_ISIS || type == PM_SUBTRACTION_KERNEL_GUNK)) { - if (!getStamps(&stamps, reference, input, subMask, weight, NULL, - threshold, stampSpacing, targetWidth, size, footprint)) { - goto ERROR; - } kernels = pmSubtractionKernelsOptimumISIS(type, size, inner, spatialOrder, optFWHMs, optOrder, - stamps, footprint, optThreshold); + stamps, footprint, optThreshold, mode); if (!kernels) { psErrorClear(); @@ -314,13 +306,13 @@ psLogMsg("psModules.imcombine", PS_LOG_INFO, "Iteration %d.", k); - if (!getStamps(&stamps, reference, input, subMask, weight, region, - threshold, stampSpacing, targetWidth, size, footprint)) { - goto ERROR; + if (!getStamps(&stamps, ro1, ro2, subMask, weight, region, threshold, stampSpacing, + size, footprint, mode)) { + goto MATCH_ERROR; } psTrace("psModules.imcombine", 3, "Calculating equation...\n"); - if (!pmSubtractionCalculateEquation(stamps, kernels, footprint)) { + if (!pmSubtractionCalculateEquation(stamps, kernels, footprint, mode)) { psError(PS_ERR_UNKNOWN, false, "Unable to calculate least-squares equation."); - goto ERROR; + goto MATCH_ERROR; } @@ -331,15 +323,27 @@ if (!solution) { psError(PS_ERR_UNKNOWN, false, "Unable to calculate least-squares equation."); - goto ERROR; + goto MATCH_ERROR; } memCheck(" solve equation"); + psVector *deviations = pmSubtractionCalculateDeviations(stamps, solution, footprint, kernels, + mode); // Deviations for each stamp + if (!deviations) { + psError(PS_ERR_UNKNOWN, false, "Unable to calculate deviations."); + goto MATCH_ERROR; + } + + memCheck(" calculate deviations"); + psTrace("psModules.imcombine", 3, "Rejecting stamps...\n"); - numRejected = pmSubtractionRejectStamps(stamps, subMask, solution, footprint, rej, kernels); + numRejected = pmSubtractionRejectStamps(stamps, deviations, subMask, rej, footprint); if (numRejected < 0) { psError(PS_ERR_UNKNOWN, false, "Unable to reject stamps."); - goto ERROR; - } + psFree(deviations); + goto MATCH_ERROR; + } + psFree(deviations); + memCheck(" reject stamps"); } @@ -350,7 +354,14 @@ if (!solution) { psError(PS_ERR_UNKNOWN, false, "Unable to calculate least-squares equation."); - goto ERROR; - } - (void)pmSubtractionRejectStamps(stamps, subMask, solution, footprint, NAN, kernels); + goto MATCH_ERROR; + } + psVector *deviations = pmSubtractionCalculateDeviations(stamps, solution, footprint, kernels, + mode); // Deviations for each stamp + if (!deviations) { + psError(PS_ERR_UNKNOWN, false, "Unable to calculate deviations."); + goto MATCH_ERROR; + } + (void)pmSubtractionRejectStamps(stamps, deviations, subMask, footprint, NAN); + psFree(deviations); } psFree(stamps); @@ -372,5 +383,5 @@ psError(PS_ERR_UNKNOWN, false, "Unable to generate kernel image."); psFree(convKernels); - goto ERROR; + goto MATCH_ERROR; } @@ -380,5 +391,5 @@ psFree(kernel); psFree(convKernels); - goto ERROR; + goto MATCH_ERROR; } psFree(kernel); @@ -416,9 +427,8 @@ psTrace("psModules.imcombine", 2, "Convolving...\n"); - if (!pmSubtractionConvolve(&convolved->image, &convolved->weight, &convolved->mask, - reference->image, reference->weight, subMask, maskBlank, region, - solution, kernels, useFFT)) { - psError(PS_ERR_UNKNOWN, false, "Unable to convolve reference image."); - goto ERROR; + if (!pmSubtractionConvolve(convolved, ro1, ro2, subMask, maskBlank, region, + solution, kernels, mode, useFFT)) { + psError(PS_ERR_UNKNOWN, false, "Unable to convolve image."); + goto MATCH_ERROR; } psFree(kernels); @@ -461,8 +471,10 @@ psFree(subMask); subMask = NULL; + psFree(weight); + weight = NULL; if (!pmSubtractionBorder(convolved->image, convolved->weight, convolved->mask, size, maskBlank)) { psError(PS_ERR_UNKNOWN, false, "Unable to set border of convolved image."); - goto ERROR; + goto MATCH_ERROR; } @@ -472,5 +484,5 @@ return true; -ERROR: +MATCH_ERROR: psFree(region); psFree(regionString); @@ -479,4 +491,131 @@ psFree(stamps); psFree(solution); + psFree(weight); return false; } + +// Calculate the second order moments for an image +static float subtractionOrderMoment(const psKernel *kernel, // Image for which to measure moments + int radius // Maximum radius + ) +{ + assert(kernel && kernel->kernel); + + int xMin = PS_MAX(kernel->xMin, -radius), xMax = PS_MIN(kernel->xMax, radius); // Bounds in x + int yMin = PS_MAX(kernel->yMin, -radius), yMax = PS_MIN(kernel->yMax, radius); // Bounds in y + + float xCentroid = 0.0, yCentroid = 0.0; // Centroid (first moment) + float sum = 0.0; // Sum (zero-th moment) + for (int y = yMin; y <= yMax; y++) { + for (int x = xMin; x <= xMax; x++) { + xCentroid += kernel->kernel[y][x] * x; + yCentroid += kernel->kernel[y][x] * y; + sum += kernel->kernel[y][x]; + } + } + xCentroid /= sum; + yCentroid /= sum; + + float eta20 = 0.0, eta02 = 0.0; // Second moments + for (int y = yMin; y <= yMax; y++) { + float yDiff = y - yCentroid; + for (int x = xMin; x <= xMax; x++) { + float xDiff = x - xCentroid; + eta20 += PS_SQR(xDiff) * kernel->kernel[y][x]; + eta02 += PS_SQR(yDiff) * kernel->kernel[y][x]; + } + } + + // Normalise to calculate the scale-invariant + float sum2 = PS_SQR(sum); + eta20 /= sum2; + eta02 /= sum2; + // eta11 /= sum2; + + return eta20 + eta02; +} + +#if 0 +// Calculate the deviations for a particular subtraction order +static psVector *subtractionOrderDeviation(float *sumKernel, // Sum of the kernel + pmSubtractionStampList *stamps, // Stamps to convolve + const pmSubtractionKernels *kernels, // Kernel basis functions + int footprint, // Stamp footprint + pmSubtractionMode mode // Mode of subtraction + ) +{ + assert(stamps); + assert(footprint >= 0); + assert(mode == PM_SUBTRACTION_MODE_1 || mode == PM_SUBTRACTION_MODE_2); + + if (!pmSubtractionCalculateEquation(stamps, kernels, footprint, mode)) { + psError(PS_ERR_UNKNOWN, false, "Unable to calculate least-squares equation."); + return NULL; + } + + psVector *solution = pmSubtractionSolveEquation(NULL, stamps); + if (!solution) { + psError(PS_ERR_UNKNOWN, false, "Unable to calculate least-squares equation."); + return NULL; + } + + if (sumKernel) { + float sum = 0.0; // Sum of the kernel + psImage *image = pmSubtractionKernelImage(solution, kernels, 0.0, 0.0); // Image of kernel + for (int y = 0; y < image->numRows; y++) { + for (int x = 0; x < image->numCols; x++) { + sum += image->data.F32[y][x]; + } + } + psFree(image); + *sumKernel = sum; + } + + psVector *deviations = pmSubtractionCalculateDeviations(stamps, solution, footprint, kernels, mode); + psFree(solution); + if (!deviations) { + psError(PS_ERR_UNKNOWN, false, "Unable to calculate deviations."); + return NULL; + } + + return deviations; +} +#endif + +pmSubtractionMode pmSubtractionOrder(pmSubtractionStampList *stamps, int radius) +{ + PS_ASSERT_INT_POSITIVE(radius, PM_SUBTRACTION_MODE_ERR); + + psVector *mask = psVectorAlloc(stamps->num, PS_TYPE_MASK); // Mask for stamps + psVector *moments = psVectorAlloc(stamps->num, PS_TYPE_F32); // Moments + 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_MASK_DATA[i] = 0xff; + continue; + } + mask->data.PS_TYPE_MASK_DATA[i] = 0; + moments->data.F32[i] = subtractionOrderMoment(stamp->image1, radius) / + subtractionOrderMoment(stamp->image2, radius); + } + + psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN); + if (!psVectorStats(stats, moments, NULL, mask, 0xff)) { + psError(PS_ERR_UNKNOWN, false, "Unable to calculate statistics for moments ratio."); + psFree(mask); + psFree(moments); + psFree(stats); + return PM_SUBTRACTION_MODE_ERR; + } + psFree(moments); + psFree(mask); + + psLogMsg("psModules.imcombine", PS_LOG_INFO, "Median ratio of second moments: %lf", stats->robustMedian); + pmSubtractionMode mode = (stats->robustMedian <= 1.0 ? PM_SUBTRACTION_MODE_1 : PM_SUBTRACTION_MODE_2); + psFree(stats); + + return mode; +} + + +