Index: trunk/psModules/src/imcombine/pmSubtractionParams.c =================================================================== --- trunk/psModules/src/imcombine/pmSubtractionParams.c (revision 15247) +++ trunk/psModules/src/imcombine/pmSubtractionParams.c (revision 15443) @@ -50,4 +50,21 @@ #endif +/// Select the appropriate convolution, given the kernel basis function and subtraction mode +static inline psKernel *selectConvolution(const pmSubtractionStamp *stamp, // Stamp + int kernelIndex, // Index for kernel component + pmSubtractionMode mode // Mode of subtraction + ) +{ + switch (mode) { + case PM_SUBTRACTION_MODE_1: + return stamp->convolutions1->data[kernelIndex]; + case PM_SUBTRACTION_MODE_2: + return stamp->convolutions2->data[kernelIndex]; + default: + psAbort("Unsupported subtraction mode: %x", mode); + } + return NULL; // Unreached +} + // Accumulate cross-term sums for a stamp static void accumulateCross(double *sumI, // Sum of I(x)/sigma(x)^2 @@ -55,14 +72,15 @@ double *sumIC, // Sum of I(x)conv(x)/sigma(x)^2 const pmSubtractionStamp *stamp, // Stamp with weight - const psKernel *input, // Input image, I(x) + const psKernel *target, // Target stamp int kernelIndex, // Index for kernel component - int footprint // Size of region of interest + int footprint, // Size of region of interest + pmSubtractionMode mode // Mode of subtraction ) { psKernel *weight = stamp->weight; // Weight, sigma(x)^2 - psKernel *convolution = stamp->convolutions->data[kernelIndex]; // Convolution of interest + psKernel *convolution = selectConvolution(stamp, kernelIndex, mode); // Convolution of interest for (int y = -footprint; y <= footprint; y++) { - psF32 *in = &input->kernel[y][-footprint]; // Dereference input + psF32 *in = &target->kernel[y][-footprint]; // Dereference input psF32 *wt = &weight->kernel[y][-footprint]; // Dereference weight psF32 *conv = &convolution->kernel[y][-footprint]; // Dereference convolution @@ -82,9 +100,10 @@ const pmSubtractionStamp *stamp, // Stamp with input and weight int kernelIndex, // Index for kernel component - int footprint // Size of region of interest + int footprint, // Size of region of interest + pmSubtractionMode mode // Mode of subtraction ) { psKernel *weight = stamp->weight; // Weight, sigma(x)^2 - psKernel *convolution = stamp->convolutions->data[kernelIndex]; // Convolution of interest + psKernel *convolution = selectConvolution(stamp, kernelIndex, mode); // Convolution of interest for (int y = -footprint; y <= footprint; y++) { @@ -100,18 +119,19 @@ } -static double accumulateChi2(psKernel *input, // Input stamp +static double accumulateChi2(const psKernel *target, // Target stamp pmSubtractionStamp *stamp, // Stamp with weight int kernelIndex, // Index for kernel component double coeff, // Coefficient of convolution double bg, // Background term - int footprint // Size of region of interest + int footprint, // Size of region of interest + pmSubtractionMode mode // Mode of subtraction ) { double chi2 = 0.0; psKernel *weight = stamp->weight; // Weight, sigma(x)^2 - psKernel *convolution = stamp->convolutions->data[kernelIndex]; // Convolution of interest + psKernel *convolution = selectConvolution(stamp, kernelIndex, mode); // Convolution of interest for (int y = -footprint; y <= footprint; y++) { - psF32 *in = &input->kernel[y][-footprint]; // Dereference input + psF32 *in = &target->kernel[y][-footprint]; // Dereference input psF32 *wt = &weight->kernel[y][-footprint]; // Dereference weight psF32 *conv = &convolution->kernel[y][-footprint]; // Dereference convolution @@ -125,17 +145,28 @@ // Return the initial value of chi^2 -static double initialChi2(psKernel *input, // Input stamp +static double initialChi2(const psKernel *target, // Target stamp const pmSubtractionStamp *stamp, // Stamp with weight - int footprint // Size of convolution + int footprint, // Size of convolution + pmSubtractionMode mode // Mode of subtraction ) { psKernel *weight = stamp->weight; // Weight map - psKernel *reference = stamp->reference; // Reference stamp + psKernel *source; // Source stamp + switch (mode) { + case PM_SUBTRACTION_MODE_1: + source = stamp->image1; + break; + case PM_SUBTRACTION_MODE_2: + source = stamp->image2; + break; + default: + psAbort("Unsupported subtraction mode: %x", mode); + } double chi2 = 0.0; // Chi^2 for (int y = -footprint; y <= footprint; y++) { - psF32 *in = &input->kernel[y][-footprint]; // Dereference input + psF32 *in = &target->kernel[y][-footprint]; // Dereference input psF32 *wt = &weight->kernel[y][-footprint]; // Dereference weight - psF32 *ref = &reference->kernel[y][-footprint]; // Derference reference + psF32 *ref = &source->kernel[y][-footprint]; // Derference reference for (int x = -footprint; x <= footprint; x++, in++, wt++, ref++) { float diff = *in - *ref; // Temporary value @@ -148,16 +179,16 @@ // Subtract a convolution from the input -static void subtractConvolution(psKernel *input, // Input stamp +static void subtractConvolution(psKernel *target, // Target stamp const pmSubtractionStamp *stamp, // Stamp with weight int kernelIndex, // Index for kernel component float coeff, // Coefficient of subtraction float bg, // Background term - int footprint // Size of region of interest - ) -{ - psKernel *convolution = stamp->convolutions->data[kernelIndex]; // Convolution of interest - + int footprint, // Size of region of interest + pmSubtractionMode mode // Mode of subtraction + ) +{ + psKernel *convolution = selectConvolution(stamp, kernelIndex, mode); // Convolution of interest for (int y = -footprint; y <= footprint; y++) { - psF32 *in = &input->kernel[y][-footprint]; // Dereference input + psF32 *in = &target->kernel[y][-footprint]; // Dereference input psF32 *conv = &convolution->kernel[y][-footprint]; // Dereference convolution for (int x = -footprint; x <= footprint; x++, in++, conv++) { @@ -173,5 +204,5 @@ int spatialOrder, const psVector *fwhms, int maxOrder, const pmSubtractionStampList *stamps, int footprint, - float tolerance) + float tolerance, pmSubtractionMode mode) { if (type != PM_SUBTRACTION_KERNEL_ISIS && type != PM_SUBTRACTION_KERNEL_GUNK) { @@ -210,5 +241,5 @@ // Need to save the stamp inputs --- we're changing the values! int numStamps = stamps->num; // Number of stamps - psArray *inputs = psArrayAlloc(numStamps); // Deep copies of the inputs + psArray *targets = psArrayAlloc(numStamps); // Deep copies of the targets psVector *badStamps = psVectorAlloc(numStamps, PS_TYPE_U8); // Mark the bad stamps psVectorInit(badStamps, 0); @@ -219,7 +250,17 @@ continue; } - psKernel *input = stamp->input; // Input image of interest - psImage *copy = psImageCopy(NULL, input->image, PS_TYPE_F32); // Copy of the image - inputs->data[i] = psKernelAllocFromImage(copy, size + footprint, size + footprint); + psKernel *target; // Target image of interest + switch (mode) { + case PM_SUBTRACTION_MODE_1: + target = stamp->image2; + break; + case PM_SUBTRACTION_MODE_2: + target = stamp->image1; + break; + default: + psAbort("Unsupported subtraction mode: %x", mode); + } + psImage *copy = psImageCopy(NULL, target->image, PS_TYPE_F32); // Copy of the image + targets->data[i] = psKernelAllocFromImage(copy, size + footprint, size + footprint); psFree(copy); // Drop reference } @@ -238,7 +279,7 @@ continue; } - if (!pmSubtractionConvolveStamp(stamp, kernels, footprint)) { + if (!pmSubtractionConvolveStamp(stamp, kernels, footprint, mode)) { psError(PS_ERR_UNKNOWN, false, "Unable to convolve stamp %d.", i); - psFree(inputs); + psFree(targets); psFree(kernels); psFree(badStamps); @@ -258,5 +299,5 @@ "Sum of 1/sigma^2 is non-finite for stamp %d (%d,%d)\n", i, (int)stamp->x, (int)stamp->y); - psFree(inputs); + psFree(targets); psFree(kernels); psFree(badStamps); @@ -265,8 +306,8 @@ for (int j = 0; j < numKernels; j++) { - accumulateConvolutions(&sumC->data.F64[j], &sumCC->data.F64[j], stamp, j, footprint); - } - - lastChi2 += initialChi2(inputs->data[i], stamp, footprint); + accumulateConvolutions(&sumC->data.F64[j], &sumCC->data.F64[j], stamp, j, footprint, mode); + } + + lastChi2 += initialChi2(targets->data[i], stamp, footprint, mode); numPixels += PS_SQR(2 * footprint + 1); } @@ -297,5 +338,5 @@ } pmSubtractionStamp *stamp = stamps->stamps->data[j]; // Stamp of interest - accumulateCross(&sumI, &sumII, &sumIC, stamp, inputs->data[j], i, footprint); + accumulateCross(&sumI, &sumII, &sumIC, stamp, targets->data[j], i, footprint, mode); } @@ -310,5 +351,5 @@ } pmSubtractionStamp *stamp = stamps->stamps->data[j]; // Stamp of interest - chi2 += accumulateChi2(inputs->data[j], stamp, i, coeff, bg, footprint); + chi2 += accumulateChi2(targets->data[j], stamp, i, coeff, bg, footprint, mode); } @@ -328,5 +369,5 @@ if (bestIndex == -1) { psError(PS_ERR_UNKNOWN, false, "Unable to find best kernel component in round %d.", iter); - psFree(inputs); + psFree(targets); psFree(sumC); psFree(sumCC); @@ -345,5 +386,5 @@ } pmSubtractionStamp *stamp = stamps->stamps->data[j]; // Stamp of interest - subtractConvolution(inputs->data[j], stamp, bestIndex, bestCoeff, bestBG, footprint); + subtractConvolution(targets->data[j], stamp, bestIndex, bestCoeff, bestBG, footprint, mode); } @@ -361,5 +402,5 @@ lastChi2 = bestChi2; } - psFree(inputs); + psFree(targets); psFree(sumC); psFree(sumCC); @@ -401,5 +442,5 @@ pmSubtractionStamp *stamp = stamps->stamps->data[j]; // Stamp of interest psArray *convolutions = convNew->data[j]; // Convolutions for this stamp - convolutions->data[rank] = psMemIncrRefCounter(stamp->convolutions->data[i]); + convolutions->data[rank] = psMemIncrRefCounter(selectConvolution(stamp, i, mode)); } } @@ -420,6 +461,18 @@ } pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest - psFree(stamp->convolutions); - stamp->convolutions = convNew->data[i]; + psFree(stamp->convolutions1); + psFree(stamp->convolutions2); + switch (mode) { + case PM_SUBTRACTION_MODE_1: + stamp->convolutions1 = convNew->data[i]; + stamp->convolutions2 = NULL; + break; + case PM_SUBTRACTION_MODE_2: + stamp->convolutions1 = NULL; + stamp->convolutions2 = convNew->data[i]; + break; + default: + psAbort("Unsupported subtraction mode: %x", mode); + } }