Index: trunk/psModules/src/imcombine/pmSubtractionKernels.c =================================================================== --- trunk/psModules/src/imcombine/pmSubtractionKernels.c (revision 25060) +++ trunk/psModules/src/imcombine/pmSubtractionKernels.c (revision 25101) @@ -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 = 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: trunk/psModules/src/imcombine/pmSubtractionKernels.h =================================================================== --- trunk/psModules/src/imcombine/pmSubtractionKernels.h (revision 25060) +++ trunk/psModules/src/imcombine/pmSubtractionKernels.h (revision 25101) @@ -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: trunk/psModules/src/imcombine/pmSubtractionMatch.c =================================================================== --- trunk/psModules/src/imcombine/pmSubtractionMatch.c (revision 25060) +++ trunk/psModules/src/imcombine/pmSubtractionMatch.c (revision 25101) @@ -869,2 +869,133 @@ return mode; } + + +#if 0 +/// 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 +} 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 *variance; ///< Variance image 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 *matrix1, *matrix2; ///< Least-squares matrices for each image, or NULL + psImage *matrixX; ///< Cross-matrix (for mode DUAL), or NULL + psVector *vector1, *vector2; ///< Least-squares vectors for each image, or NULL + pmSubtractionStampStatus status; ///< Status of stamp +} pmSubtractionStamp; + +/// Kernels specification +typedef struct { + pmSubtractionKernelsType type; ///< Type of kernels --- allowing the use of multiple kernels + psString description; ///< Description of the kernel parameters + long num; ///< Number of kernel components (not including the spatial ones) + psVector *u, *v; ///< Offset (for POIS) or polynomial order (for ISIS) + psVector *widths; ///< Gaussian FWHMs (ISIS) + psVector *uStop, *vStop; ///< Width of kernel element (SPAM,FRIES only) + psArray *preCalc; ///< Array of images containing pre-calculated kernel (for ISIS) + float penalty; ///< Penalty for wideness + psVector *penalties; ///< Penalty for each kernel component + 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 + int numCols, numRows; ///< Size of image (for normalisation), or zero to use image provided + 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 +} pmSubtractionKernels; + +// 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 + ) +{ + 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."); + kernels->mode = oldMode; + 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, footprint); + if (numRejected < 0) { + psError(PS_ERR_UNKNOWN, false, "Unable to reject stamps."); + psFree(deviations); + return false; + } + psFree(deviations); + + if (numRejected > 0) { + psTrace("psModules.imcombine", 3, "Solving equation...\n"); + if (!pmSubtractionSolveEquation(kernels, stamps)) { + psError(PS_ERR_UNKNOWN, false, "Unable to calculate least-squares equation."); + return false; + } + psVector *deviations = pmSubtractionCalculateDeviations(stamps, kernels); // Stamp deviations + if (!deviations) { + psError(PS_ERR_UNKNOWN, false, "Unable to calculate deviations."); + return false; + } + psFree(deviations); + } + + return true; +} + + +pmSubtractionMode pmSubtractionBestMode(pmSubtractionStampList *stamps, pmSubtractionKernels *kernels) +{ + 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 so we can try each way + pmSubtractionStampList *stamps1 = pmSubtractionStampsListCopy(stamps); + pmSubtractionStampList *stamps2 = pmSubtractionStampsListCopy(stamps); + + pmSubtractionKernels *kernels1 = pmSubtractionKernelsCopy(kernels); + pmSubtractionKernels *kernels2 = pmSubtractionKernelsCopy(kernels); + + kernels1->mode = PM_SUBTRACTION_MODE_1; + kernels2->mode = PM_SUBTRACTION_MODE_2; + + + subtractionModeTest(stamps1, kernels1, "forward"); + subtractionModeTest(stamps2, kernels2, "backward"); + + // XXX Compare kernels1->mean, kernels2->mean +} + +#endif Index: trunk/psModules/src/imcombine/pmSubtractionStamps.c =================================================================== --- trunk/psModules/src/imcombine/pmSubtractionStamps.c (revision 25060) +++ trunk/psModules/src/imcombine/pmSubtractionStamps.c (revision 25101) @@ -225,4 +225,64 @@ } +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 = psKernelCopy(inStamp->image1); + outStamp->image2 = psKernelCopy(inStamp->image2); + outStamp->variance = psKernelCopy(inStamp->variance); + + if (inStamp->convolutions1) { + int size = inStamp->convolutions1->n; // Size of array + outStamp->convolutions1 = psArrayAlloc(size); + for (int j = 0; j < size; j++) { + outStamp->convolutions1->data[j] = psKernelCopy(inStamp->convolutions1->data[j]); + } + } + if (inStamp->convolutions2) { + int size = inStamp->convolutions2->n; // Size of array + outStamp->convolutions2 = psArrayAlloc(size); + for (int j = 0; j < size; j++) { + outStamp->convolutions2->data[j] = psKernelCopy(inStamp->convolutions2->data[j]); + } + } + + 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: trunk/psModules/src/imcombine/pmSubtractionStamps.h =================================================================== --- trunk/psModules/src/imcombine/pmSubtractionStamps.h (revision 25060) +++ trunk/psModules/src/imcombine/pmSubtractionStamps.h (revision 25101) @@ -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