Index: /branches/czw_branch/cleanup/psModules/src/imcombine/pmSubtraction.c =================================================================== --- /branches/czw_branch/cleanup/psModules/src/imcombine/pmSubtraction.c (revision 25147) +++ /branches/czw_branch/cleanup/psModules/src/imcombine/pmSubtraction.c (revision 25148) @@ -733,5 +733,5 @@ int pmSubtractionRejectStamps(pmSubtractionKernels *kernels, pmSubtractionStampList *stamps, - const psVector *deviations, psImage *subMask, float sigmaRej, int footprint) + const psVector *deviations, psImage *subMask, float sigmaRej) { PM_ASSERT_SUBTRACTION_KERNELS_NON_NULL(kernels, false); @@ -821,4 +821,5 @@ ds9num++; + int footprint = stamps->footprint; // Half-size of stamp region of interest int numRejected = 0; // Number of stamps rejected int numGood = 0; // Number of good stamps Index: /branches/czw_branch/cleanup/psModules/src/imcombine/pmSubtraction.h =================================================================== --- /branches/czw_branch/cleanup/psModules/src/imcombine/pmSubtraction.h (revision 25147) +++ /branches/czw_branch/cleanup/psModules/src/imcombine/pmSubtraction.h (revision 25148) @@ -68,6 +68,5 @@ const psVector *deviations, ///< Deviations for each stamp psImage *subMask, ///< Subtraction mask - float sigmaRej, ///< Number of RMS deviations above zero at which to reject - int footprint ///< Half-size of stamp + float sigmaRej ///< Number of RMS deviations above zero at which to reject ); Index: /branches/czw_branch/cleanup/psModules/src/imcombine/pmSubtractionAnalysis.c =================================================================== --- /branches/czw_branch/cleanup/psModules/src/imcombine/pmSubtractionAnalysis.c (revision 25147) +++ /branches/czw_branch/cleanup/psModules/src/imcombine/pmSubtractionAnalysis.c (revision 25148) @@ -17,10 +17,10 @@ -bool pmSubtractionAnalysis(psMetadata *analysis, pmSubtractionKernels *kernels, psRegion *region, +bool pmSubtractionAnalysis(psMetadata *analysis, psMetadata *header, + pmSubtractionKernels *kernels, psRegion *region, int numCols, int numRows) { - if (analysis) { - PS_ASSERT_METADATA_NON_NULL(analysis, false); - } + PS_ASSERT_METADATA_NON_NULL(analysis, false); + PS_ASSERT_METADATA_NON_NULL(header, false); PM_ASSERT_SUBTRACTION_KERNELS_NON_NULL(kernels, false); PM_ASSERT_SUBTRACTION_KERNELS_SOLUTION(kernels, false); @@ -38,5 +38,11 @@ PS_DATA_REGION | PS_META_DUPLICATE_OK, "Region over which subtraction was performed", subRegion); + + psString string = psRegionToString(*subRegion); psFree(subRegion); + + psMetadataAddStr(header, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_REGION, PS_META_DUPLICATE_OK, + "Region over which subtraction was performed", string); + psFree(string); } @@ -45,5 +51,7 @@ PS_DATA_UNKNOWN | PS_META_DUPLICATE_OK, "Subtraction kernels", kernels); psMetadataAddS32(analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_MODE, - PS_META_DUPLICATE_OK, "Subtraction kernels", kernels->mode); + PS_META_DUPLICATE_OK, "Subtraction mode", kernels->mode); + psMetadataAddS32(header, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_MODE, + PS_META_DUPLICATE_OK, "Subtraction mode", kernels->mode); // Realisations of kernel @@ -163,10 +171,22 @@ psFree(image); - psMetadataItem *item = psMetadataLookup(analysis, PM_SUBTRACTION_ANALYSIS_DECONV_MAX); // Previous - if (item) { - item->data.F32 = PS_MAX(item->data.F32, max); - } else { - psMetadataAddF32(analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_DECONV_MAX, 0, - "Maximum deconvolution fraction", max); + { + psMetadataItem *item = psMetadataLookup(analysis, PM_SUBTRACTION_ANALYSIS_DECONV_MAX); // Previous + if (item) { + max = item->data.F32 = PS_MAX(item->data.F32, max); + } else { + psMetadataAddF32(analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_DECONV_MAX, 0, + "Maximum deconvolution fraction", max); + } + } + + { + psMetadataItem *item = psMetadataLookup(header, PM_SUBTRACTION_ANALYSIS_DECONV_MAX); // Previous + if (item) { + item->data.F32 = max; + } else { + psMetadataAddF32(header, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_DECONV_MAX, 0, + "Maximum deconvolution fraction", max); + } } } @@ -216,4 +236,17 @@ psMetadataAddF32(analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_MYY, PS_META_DUPLICATE_OK, "Moment in yy", m02); + + psMetadataAddF32(header, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_NORM, + PS_META_DUPLICATE_OK, "Normalisation", m00); + psMetadataAddF32(header, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_MX, + PS_META_DUPLICATE_OK, "Moment in x", m10); + psMetadataAddF32(header, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_MY, + PS_META_DUPLICATE_OK, "Moment in y", m01); + psMetadataAddF32(header, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_MXX, + PS_META_DUPLICATE_OK, "Moment in xx", m20); + psMetadataAddF32(header, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_MXY, + PS_META_DUPLICATE_OK, "Moment in xy", m11); + psMetadataAddF32(header, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_MYY, + PS_META_DUPLICATE_OK, "Moment in yy", m02); } @@ -225,4 +258,6 @@ psMetadataAddF32(analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_BGDIFF, + PS_META_DUPLICATE_OK, "Background difference", bg); + psMetadataAddF32(header, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_BGDIFF, PS_META_DUPLICATE_OK, "Background difference", bg); psFree(polyValues); @@ -237,4 +272,11 @@ psMetadataAddF32(analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_DEV_RMS, 0, "RMS stamp deviation", kernels->rms); + + psMetadataAddS32(header, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_STAMPS, 0, "Number of stamps", + kernels->numStamps); + psMetadataAddF32(header, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_DEV_MEAN, 0, "Mean stamp deviation", + kernels->mean); + psMetadataAddF32(header, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_DEV_RMS, 0, "RMS stamp deviation", + kernels->rms); } Index: /branches/czw_branch/cleanup/psModules/src/imcombine/pmSubtractionKernels.c =================================================================== --- /branches/czw_branch/cleanup/psModules/src/imcombine/pmSubtractionKernels.c (revision 25147) +++ /branches/czw_branch/cleanup/psModules/src/imcombine/pmSubtractionKernels.c (revision 25148) @@ -765,2 +765,33 @@ return PM_SUBTRACTION_KERNEL_NONE; } + +pmSubtractionKernels *pmSubtractionKernelsCopy(const pmSubtractionKernels *in) +{ + PM_ASSERT_SUBTRACTION_KERNELS_NON_NULL(in, NULL); + + pmSubtractionKernels *out = psAlloc(sizeof(pmSubtractionKernels)); // Kernels, to return + psMemSetDeallocator(out, (psFreeFunc)subtractionKernelsFree); + + out->type = in->type; + out->description = psMemIncrRefCounter(in->description); + out->num = in->num; + out->u = psMemIncrRefCounter(in->u); + out->v = psMemIncrRefCounter(in->v); + out->widths = psMemIncrRefCounter(in->widths); + out->preCalc = psMemIncrRefCounter(in->preCalc); + out->penalty = in->penalty; + out->penalties = psMemIncrRefCounter(in->penalties); + out->uStop = psMemIncrRefCounter(in->uStop); + out->vStop = psMemIncrRefCounter(in->vStop); + out->size = in->size; + out->inner = in->inner; + out->spatialOrder = in->spatialOrder; + out->bgOrder = in->bgOrder; + out->mode = in->mode; + out->numCols = in->numCols; + out->numRows = in->numRows; + out->solution1 = in->solution1 ? psVectorCopy(NULL, in->solution1, PS_TYPE_F64) : NULL; + out->solution2 = in->solution2 ? psVectorCopy(NULL, in->solution2, PS_TYPE_F64) : NULL; + + return out; +} Index: /branches/czw_branch/cleanup/psModules/src/imcombine/pmSubtractionKernels.h =================================================================== --- /branches/czw_branch/cleanup/psModules/src/imcombine/pmSubtractionKernels.h (revision 25147) +++ /branches/czw_branch/cleanup/psModules/src/imcombine/pmSubtractionKernels.h (revision 25148) @@ -203,4 +203,11 @@ ); +/// Copy kernels +/// +/// A deep copy is performed on the solution only; the other components are merely pointers. +pmSubtractionKernels *pmSubtractionKernelsCopy( + const pmSubtractionKernels *in // Kernels to copy + ); + #endif Index: /branches/czw_branch/cleanup/psModules/src/imcombine/pmSubtractionMatch.c =================================================================== --- /branches/czw_branch/cleanup/psModules/src/imcombine/pmSubtractionMatch.c (revision 25147) +++ /branches/czw_branch/cleanup/psModules/src/imcombine/pmSubtractionMatch.c (revision 25148) @@ -447,5 +447,24 @@ } + + // 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, 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, subMode); + } + + memCheck("kernels"); + if (subMode == PM_SUBTRACTION_MODE_UNSURE) { +#if 0 // Get backgrounds psStats *bgStats = psStatsAlloc(BG_STAT); // Statistics for background @@ -469,4 +488,6 @@ pmSubtractionMode newMode = pmSubtractionOrder(stamps, bg1, bg2); // Subtraction mode to use +#endif + pmSubtractionMode newMode = pmSubtractionBestMode(&stamps, &kernels, subMask, rej); switch (newMode) { case PM_SUBTRACTION_MODE_1: @@ -483,21 +504,4 @@ } - // 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, 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, subMode); - } - - memCheck("kernels"); - int numRejected = -1; // Number of rejected stamps in each iteration for (int k = 0; k < iter && numRejected != 0; k++) { @@ -535,5 +539,5 @@ psTrace("psModules.imcombine", 3, "Rejecting stamps...\n"); - numRejected = pmSubtractionRejectStamps(kernels, stamps, deviations, subMask, rej, footprint); + numRejected = pmSubtractionRejectStamps(kernels, stamps, deviations, subMask, rej); if (numRejected < 0) { psError(PS_ERR_UNKNOWN, false, "Unable to reject stamps."); @@ -557,5 +561,5 @@ goto MATCH_ERROR; } - pmSubtractionRejectStamps(kernels, stamps, deviations, subMask, NAN, footprint); + pmSubtractionRejectStamps(kernels, stamps, deviations, subMask, NAN); psFree(deviations); } @@ -842,2 +846,137 @@ return mode; } + + +// 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 + ) +{ + assert(stamps); + assert(kernels); + + psTrace("psModules.imcombine", 3, "Calculating %s equation...\n", description); + if (!pmSubtractionCalculateEquation(stamps, kernels)) { + psError(PS_ERR_UNKNOWN, false, "Unable to calculate least-squares equation."); + return false; + } + + psTrace("psModules.imcombine", 3, "Solving %s equation...\n", description); + if (!pmSubtractionSolveEquation(kernels, stamps)) { + psError(PS_ERR_UNKNOWN, false, "Unable to calculate least-squares equation."); + return false; + } + + psTrace("psModules.imcombine", 3, "Calculate %s deviations...\n", description); + psVector *deviations = pmSubtractionCalculateDeviations(stamps, kernels); // Stamp deviations + if (!deviations) { + psError(PS_ERR_UNKNOWN, false, "Unable to calculate deviations."); + return false; + } + + psTrace("psModules.imcombine", 3, "Rejecting %s stamps...\n", description); + long numRejected = pmSubtractionRejectStamps(kernels, stamps, deviations, subMask, rej); + if (numRejected < 0) { + psError(PS_ERR_UNKNOWN, 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, "Resolving %s equation...\n", description); + if (!pmSubtractionSolveEquation(kernels, stamps)) { + psError(PS_ERR_UNKNOWN, 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(PS_ERR_UNKNOWN, 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(PS_ERR_UNKNOWN, false, "Unable to reject stamps."); + psFree(deviations); + return false; + } + psFree(deviations); + } + + return true; +} + + +pmSubtractionMode pmSubtractionBestMode(pmSubtractionStampList **stamps, pmSubtractionKernels **kernels, + const psImage *subMask, float rej) +{ + PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(*stamps, PM_SUBTRACTION_MODE_ERR); + PM_ASSERT_SUBTRACTION_KERNELS_NON_NULL(*kernels, PM_SUBTRACTION_MODE_ERR); + + // Copies of the inputs + pmSubtractionStampList *stamps1 = pmSubtractionStampListCopy(*stamps); + pmSubtractionKernels *kernels1 = pmSubtractionKernelsCopy(*kernels); + psImage *subMask1 = psImageCopy(NULL, subMask, subMask->type.type); + kernels1->mode = PM_SUBTRACTION_MODE_1; + + if (!subtractionModeTest(stamps1, kernels1, "forward", subMask1, rej)) { + psError(PS_ERR_UNKNOWN, false, "Unable to test forward subtraction"); + psFree(stamps1); + psFree(kernels1); + psFree(subMask1); + return PM_SUBTRACTION_MODE_ERR; + } + psFree(subMask1); + + // Copies of the inputs + pmSubtractionStampList *stamps2 = pmSubtractionStampListCopy(*stamps); + pmSubtractionKernels *kernels2 = pmSubtractionKernelsCopy(*kernels); + psImage *subMask2 = psImageCopy(NULL, subMask, subMask->type.type); + kernels2->mode = PM_SUBTRACTION_MODE_2; + + if (!subtractionModeTest(stamps2, kernels2, "backward", subMask2, rej)) { + psError(PS_ERR_UNKNOWN, false, "Unable to test backward subtraction"); + psFree(stamps2); + psFree(kernels2); + psFree(subMask2); + psFree(stamps1); + psFree(kernels1); + return PM_SUBTRACTION_MODE_ERR; + } + psFree(subMask2); + + + pmSubtractionStampList *bestStamps = NULL; // Best choice for stamps + pmSubtractionKernels *bestKernels = NULL; // Best choice for kernels + psLogMsg("psModules.imcombine", PS_LOG_INFO, + "Forward: %f +/- %f from %d stamps\nBackward: %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; + } else { + bestStamps = stamps2; + bestKernels = kernels2; + } + + psFree(*stamps); + psFree(*kernels); + *stamps = psMemIncrRefCounter(bestStamps); + *kernels = psMemIncrRefCounter(bestKernels); + + psFree(stamps1); + psFree(stamps2); + psFree(kernels1); + psFree(kernels2); + + return bestKernels->mode; +} + Index: /branches/czw_branch/cleanup/psModules/src/imcombine/pmSubtractionMatch.h =================================================================== --- /branches/czw_branch/cleanup/psModules/src/imcombine/pmSubtractionMatch.h (revision 25147) +++ /branches/czw_branch/cleanup/psModules/src/imcombine/pmSubtractionMatch.h (revision 25148) @@ -83,4 +83,13 @@ ); +/// Determine best subtraction mode to use +/// +/// Subtractions are attempted each way, and the mode with the lower residual is taken to be the best +pmSubtractionMode pmSubtractionBestMode( + pmSubtractionStampList **stamps, ///< Stamps to use for solution + pmSubtractionKernels **kernels, ///< Kernels to use for solution + const psImage *subMask, ///< Subtraction mask + float rej ///< Rejection threshold for stamps + ); #endif Index: /branches/czw_branch/cleanup/psModules/src/imcombine/pmSubtractionStamps.c =================================================================== --- /branches/czw_branch/cleanup/psModules/src/imcombine/pmSubtractionStamps.c (revision 25147) +++ /branches/czw_branch/cleanup/psModules/src/imcombine/pmSubtractionStamps.c (revision 25148) @@ -225,4 +225,70 @@ } +pmSubtractionStampList *pmSubtractionStampListCopy(const pmSubtractionStampList *in) +{ + PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(in, NULL); + + pmSubtractionStampList *out = psAlloc(sizeof(pmSubtractionStampList)); // Copied stamp list to return + psMemSetDeallocator(out, (psFreeFunc)subtractionStampListFree); + + int num = out->num = in->num; // Number of stamps + out->stamps = psArrayAlloc(num); + out->regions = psArrayAlloc(num); + out->x = NULL; + out->y = NULL; + out->flux = NULL; + out->footprint = in->footprint; + + for (int i = 0; i < num; i++) { + psRegion *inRegion = in->regions->data[i]; // Input region + out->regions->data[i] = psRegionAlloc(inRegion->x0, inRegion->x1, inRegion->y0, inRegion->y1); + + pmSubtractionStamp *inStamp = in->stamps->data[i]; // Input stamp + pmSubtractionStamp *outStamp = psAlloc(sizeof(pmSubtractionStamp)); + psMemSetDeallocator(outStamp, (psFreeFunc)subtractionStampFree); + outStamp->x = inStamp->x; + outStamp->y = inStamp->y; + outStamp->flux = inStamp->flux; + outStamp->xNorm = inStamp->xNorm; + outStamp->yNorm = inStamp->yNorm; + outStamp->status = inStamp->status; + + outStamp->image1 = inStamp->image1 ? psKernelCopy(inStamp->image1) : NULL; + outStamp->image2 = inStamp->image2 ? psKernelCopy(inStamp->image2) : NULL; + outStamp->variance = inStamp->variance ? psKernelCopy(inStamp->variance) : NULL; + + if (inStamp->convolutions1) { + int size = inStamp->convolutions1->n; // Size of array + outStamp->convolutions1 = psArrayAlloc(size); + for (int j = 0; j < size; j++) { + psKernel *conv = inStamp->convolutions1->data[j]; // Convolution + outStamp->convolutions1->data[j] = conv ? psKernelCopy(conv) : NULL; + } + } else { + outStamp->convolutions1 = NULL; + } + if (inStamp->convolutions2) { + int size = inStamp->convolutions2->n; // Size of array + outStamp->convolutions2 = psArrayAlloc(size); + for (int j = 0; j < size; j++) { + psKernel *conv = inStamp->convolutions2->data[j]; // Convolution + outStamp->convolutions2->data[j] = conv ? psKernelCopy(conv) : NULL; + } + } else { + outStamp->convolutions2 = NULL; + } + + outStamp->matrix1 = inStamp->matrix1 ? psImageCopy(NULL, inStamp->matrix1, PS_TYPE_F64) : NULL; + outStamp->matrix2 = inStamp->matrix2 ? psImageCopy(NULL, inStamp->matrix2, PS_TYPE_F64) : NULL; + outStamp->matrixX = inStamp->matrixX ? psImageCopy(NULL, inStamp->matrixX, PS_TYPE_F64) : NULL; + outStamp->vector1 = inStamp->vector1 ? psVectorCopy(NULL, inStamp->vector1, PS_TYPE_F64) : NULL; + outStamp->vector2 = inStamp->vector2 ? psVectorCopy(NULL, inStamp->vector2, PS_TYPE_F64) : NULL; + + out->stamps->data[i] = outStamp; + } + + return out; +} + pmSubtractionStamp *pmSubtractionStampAlloc(void) { Index: /branches/czw_branch/cleanup/psModules/src/imcombine/pmSubtractionStamps.h =================================================================== --- /branches/czw_branch/cleanup/psModules/src/imcombine/pmSubtractionStamps.h (revision 25147) +++ /branches/czw_branch/cleanup/psModules/src/imcombine/pmSubtractionStamps.h (revision 25148) @@ -52,4 +52,12 @@ } \ } + +/// Copy a list of stamps +/// +/// A deep copy is performed of the stamp list and the component stamps +pmSubtractionStampList *pmSubtractionStampListCopy( + const pmSubtractionStampList *in // Stamp list to copy + ); + /// A stamp for image subtraction