Index: trunk/psModules/src/imcombine/Makefile.am =================================================================== --- trunk/psModules/src/imcombine/Makefile.am (revision 15745) +++ trunk/psModules/src/imcombine/Makefile.am (revision 15756) @@ -9,8 +9,11 @@ pmStackReject.c \ pmSubtraction.c \ + pmSubtractionEquation.c \ pmSubtractionKernels.c \ + pmSubtractionMask.c \ + pmSubtractionMatch.c \ + pmSubtractionParams.c \ pmSubtractionStamps.c \ - pmSubtractionMatch.c \ - pmSubtractionParams.c + pmPSFEnvelope.c pkginclude_HEADERS = \ @@ -20,8 +23,11 @@ pmStackReject.h \ pmSubtraction.h \ + pmSubtractionEquation.h \ pmSubtractionKernels.h \ + pmSubtractionMask.h \ + pmSubtractionMatch.h \ + pmSubtractionParams.h \ pmSubtractionStamps.h \ - pmSubtractionMatch.h \ - pmSubtractionParams.h + pmPSFEnvelope.h CLEANFILES = *~ Index: trunk/psModules/src/imcombine/pmStackReject.c =================================================================== --- trunk/psModules/src/imcombine/pmStackReject.c (revision 15745) +++ trunk/psModules/src/imcombine/pmStackReject.c (revision 15756) @@ -57,7 +57,5 @@ for (int i = 0; i < numRegions; i++) { psRegion *region = regions->data[i]; // Region of interest - psVector *solution = solutions->data[i]; // Solution of interest - if (!pmSubtractionConvolve(convRO, inRO, NULL, NULL, 0, region, solution, kernels, - PM_SUBTRACTION_MODE_1, true)) { + if (!pmSubtractionConvolve(convRO, NULL, inRO, NULL, NULL, 0, region, kernels, true)) { psError(PS_ERR_UNKNOWN, false, "Unable to convolve mask image in region %d.", i); psFree(convRO); @@ -72,5 +70,5 @@ float xNorm = (region->x0 + 0.5 * (region->x1 - region->x0) - numCols/2.0) / (float)numCols; float yNorm = (region->y0 + 0.5 * (region->y1 - region->y0) - numRows/2.0) / (float)numRows; - psImage *kernel = pmSubtractionKernelImage(solution, kernels, xNorm, yNorm); + psImage *kernel = pmSubtractionKernelImage(kernels, xNorm, yNorm, false); if (!kernel) { psError(PS_ERR_UNKNOWN, false, "Unable to generate kernel image."); Index: trunk/psModules/src/imcombine/pmSubtraction.c =================================================================== --- trunk/psModules/src/imcombine/pmSubtraction.c (revision 15745) +++ trunk/psModules/src/imcombine/pmSubtraction.c (revision 15756) @@ -4,6 +4,6 @@ * @author GLG, MHPCC * - * @version $Revision: 1.73 $ $Name: not supported by cvs2svn $ - * @date $Date: 2007-11-05 23:43:38 $ + * @version $Revision: 1.74 $ $Name: not supported by cvs2svn $ + * @date $Date: 2007-12-07 01:57:15 $ * * Copyright 2004-2007 Institute for Astronomy, University of Hawaii @@ -23,4 +23,5 @@ #include "pmFPA.h" #include "pmSubtractionStamps.h" +#include "pmSubtractionEquation.h" #include "pmSubtraction.h" @@ -34,44 +35,4 @@ // Private (file-static) functions ////////////////////////////////////////////////////////////////////////////////////////////////////////////// - -// Return the number of polynomial terms there are for a given order -static inline long polyTerms(int order // Order of polynomial - ) -{ - return (order + 1) * (order + 2) / 2; -} - -// Given a stamp coordinates (x,y), generate a matrix where the elements (i,j) are x^i * y^j -static psImage *spatialPolyValues(int spatialOrder, // Maximum spatial polynomial order - float x, float y // Normalised position of interest, [-1,1] - ) -{ - assert(spatialOrder >= 0); - assert(x >= -1 && x <= 1); - assert(y >= -1 && y <= 1); - - psImage *polyValues = psImageAlloc(spatialOrder + 1, spatialOrder + 1, PS_TYPE_F64); // Matrix to return - polyValues->data.F64[0][0] = 1.0; - - double value = 1.0; - for (int i = 1; i <= spatialOrder; i++) { - value *= x; - polyValues->data.F64[0][i] = value; - } - - value = 1.0; - for (int j = 1; j <= spatialOrder; j++) { - value *= y; - polyValues->data.F64[j][0] = value; - } - - for (int j = 1; j <= spatialOrder; j++) { - for (int i = 1; i <= spatialOrder - j; i++) { - polyValues->data.F64[j][i] = polyValues->data.F64[j][0] * polyValues->data.F64[0][i]; - } - } - - return polyValues; -} // Generate the kernel to apply to the variance from the normal kernel @@ -108,24 +69,14 @@ // Generate an image of the solved kernel -static inline psKernel *solvedKernel(psKernel *kernel, // Kernel, to return - const psVector *solution, // Solution to the least-squares problem - const pmSubtractionKernels *kernels, // Kernel basis functions - const psImage *polyValues // Spatial polynomial values - ) -{ +static psKernel *solvedKernel(psKernel *kernel, // Kernel, to return + const pmSubtractionKernels *kernels, // Kernel basis functions + const psImage *polyValues, // Spatial polynomial values + bool wantDual // Want the dual (second) kernel? + ) +{ + assert(kernels); + assert(polyValues); + int numKernels = kernels->num; // Number of kernel basis functions - assert(kernels->u->n == numKernels); - assert(kernels->v->n == numKernels); - int spatialOrder = kernels->spatialOrder; // Maximum spatial polynomial order - assert(solution->n == numKernels * polyTerms(spatialOrder) + polyTerms(kernels->bgOrder)); - - // Ensure the subIndex for POIS kernels is what is expected - assert((kernels->type != PM_SUBTRACTION_KERNEL_POIS && - kernels->type != PM_SUBTRACTION_KERNEL_SPAM && - kernels->type != PM_SUBTRACTION_KERNEL_FRIES && - kernels->type != PM_SUBTRACTION_KERNEL_GUNK && - kernels->type != PM_SUBTRACTION_KERNEL_RINGS) || - (kernels->u->data.S32[kernels->subIndex] == 0 && kernels->v->data.S32[kernels->subIndex] == 0)); - int size = kernels->size; // Kernel half-size if (!kernel) { @@ -135,10 +86,5 @@ for (int i = 0; i < numKernels; i++) { - double value = 0.0; // Value of this component, from the sum of multiple polynomial terms - for (int yOrder = 0, index = i; yOrder <= spatialOrder; yOrder++) { - for (int xOrder = 0; xOrder <= spatialOrder - yOrder; xOrder++, index += numKernels) { - value += polyValues->data.F64[yOrder][xOrder] * solution->data.F64[index]; - } - } + double value = p_pmSubtractionSolutionCoeff(kernels, polyValues, i, wantDual); // Polynomial value switch (kernels->type) { @@ -147,8 +93,5 @@ int v = kernels->v->data.S32[i]; // Offset in y kernel->kernel[v][u] += value; - if (kernels->spatialOrder > 0 && i != kernels->subIndex) { - // The (0,0) element is subtracted from most kernels to preserve photometric scaling - kernel->kernel[0][0] -= value; - } + kernel->kernel[0][0] -= value; break; } @@ -167,8 +110,5 @@ for (int u = uStart; u <= uStop; u++) { kernel->kernel[v][u] += norm * value; - if (kernels->spatialOrder > 0 && i != kernels->subIndex) { - // The (0,0) element is subtracted from most kernels to preserve photometric scaling - kernel->kernel[0][0] -= value; - } + kernel->kernel[0][0] -= value; } } @@ -188,12 +128,8 @@ } else { // Using delta function - bool subtract = (kernels->spatialOrder > 0 && i != kernels->subIndex); // Subtract (0,0)? int u = kernels->u->data.S32[i]; // Offset in x int v = kernels->v->data.S32[i]; // Offset in y kernel->kernel[v][u] += value; - if (subtract) { - // The (0,0) element is subtracted from most kernels to preserve photometric scaling - kernel->kernel[0][0] -= value; - } + kernel->kernel[0][0] -= value; } break; @@ -211,8 +147,4 @@ } case PM_SUBTRACTION_KERNEL_RINGS: { - if (i == kernels->subIndex) { - kernel->kernel[0][0] += value; - break; - } psArray *preCalc = kernels->preCalc->data[i]; // Precalculated data psVector *uCoords = preCalc->data[0]; // u coordinates @@ -232,4 +164,10 @@ } } + + if (!wantDual) { + // Put in the normalisation component + kernel->kernel[0][0] += p_pmSubtractionSolutionNorm(kernels); + } + return kernel; @@ -328,19 +266,39 @@ } -// Mark a pixel as blank in the image, mask and weight -static inline void markBlank(psImage *image, // Image to mark as blank - psImage *mask, // Mask to mark as blank (or NULL) - psImage *weight, // Weight map to mark as blank (or NULL) - int x, int y, // Coordinates to mark blank - psMaskType blank // Blank mask value +// Convolve a region of an image +static inline void convolveRegion(psImage *convImage, // Convolved image (output) + psImage *convWeight, // Convolved weight map (output), or NULL + psKernel **kernelImage, // Convolution kernel for the image + psKernel **kernelWeight, // Convolution kernel for the weight map, or NULL + const psImage *image, // Image to convolve + const psImage *weight, // Weight map to convolve, or NULL + const psImage *subMask, // Subtraction mask + psMaskType maskVal, // Mask value to apply in convolution + const pmSubtractionKernels *kernels, // Kernels + psImage *polyValues, // Polynomial values + float background, // Background value to apply + psRegion region, // Region to convolve + bool useFFT // Use FFT to convolve? ) { - image->data.F32[y][x] = NAN; - if (mask) { - mask->data.PS_TYPE_MASK_DATA[y][x] |= blank; - } + *kernelImage = solvedKernel(*kernelImage, kernels, polyValues, false); if (weight) { - weight->data.F32[y][x] = NAN; - } + *kernelWeight = varianceKernel(*kernelWeight, *kernelImage); + } + + if (useFFT) { + // Use Fast Fourier Transform to do the convolution + // This provides a big speed-up for large kernels + convolveFFT(convImage, image, subMask, maskVal, *kernelImage, region, background, kernels->size); + if (weight) { + convolveFFT(convWeight, weight, subMask, maskVal, *kernelWeight, region, 0.0, kernels->size); + } + } else { + convolveDirect(convImage, image, *kernelImage, region, background, kernels->size); + if (weight) { + convolveDirect(convWeight, weight, *kernelWeight, region, 0.0, kernels->size); + } + } + return; } @@ -350,5 +308,4 @@ ////////////////////////////////////////////////////////////////////////////////////////////////////////////// -// Generate the convolution given a precalculated kernel psKernel *p_pmSubtractionConvolveStampPrecalc(const psKernel *image, const psKernel *kernel) { @@ -363,4 +320,33 @@ } +psImage *p_pmSubtractionPolynomial(psImage *output, int spatialOrder, float x, float y) +{ + assert(spatialOrder >= 0); + assert(x >= -1 && x <= 1); + assert(y >= -1 && y <= 1); + + output = psImageRecycle(output, spatialOrder + 1, spatialOrder + 1, PS_TYPE_F64); + output->data.F64[0][0] = 1.0; + + double value = 1.0; + for (int i = 1; i <= spatialOrder; i++) { + value *= x; + output->data.F64[0][i] = value; + } + + value = 1.0; + for (int j = 1; j <= spatialOrder; j++) { + value *= y; + output->data.F64[j][0] = value; + } + + for (int j = 1; j <= spatialOrder; j++) { + for (int i = 1; i <= spatialOrder - j; i++) { + output->data.F64[j][i] = output->data.F64[j][0] * output->data.F64[0][i]; + } + } + + return output; +} ////////////////////////////////////////////////////////////////////////////////////////////////////////////// @@ -368,166 +354,4 @@ ////////////////////////////////////////////////////////////////////////////////////////////////////////////// -psImage *pmSubtractionMask(const psImage *mask1, const psImage *mask2, psMaskType maskVal, - int size, int footprint, float badFrac, bool useFFT) -{ - PS_ASSERT_IMAGE_NON_NULL(mask1, NULL); - PS_ASSERT_IMAGE_TYPE(mask1, PS_TYPE_MASK, NULL); - if (mask2) { - PS_ASSERT_IMAGE_NON_NULL(mask2, NULL); - PS_ASSERT_IMAGE_TYPE(mask2, PS_TYPE_MASK, NULL); - PS_ASSERT_IMAGES_SIZE_EQUAL(mask2, mask1, NULL); - } - PS_ASSERT_INT_NONNEGATIVE(size, NULL); - PS_ASSERT_INT_NONNEGATIVE(footprint, NULL); - if (isfinite(badFrac)) { - PS_ASSERT_FLOAT_LARGER_THAN(badFrac, 0.0, NULL); - PS_ASSERT_FLOAT_LESS_THAN_OR_EQUAL(badFrac, 1.0, NULL); - } - - int numCols = mask1->numCols, numRows = mask1->numRows; // Size of the images - - // Dereference inputs for convenience - psMaskType **data1 = mask1->data.PS_TYPE_MASK_DATA; - psMaskType **data2 = NULL; - if (mask2) { - data2 = mask2->data.PS_TYPE_MASK_DATA; - } - - // First, a pass through to determine the fraction of bad pixels - if (isfinite(badFrac) && badFrac != 1.0) { - int numBad = 0; // Number of bad pixels - for (int y = 0; y < numRows; y++) { - for (int x = 0; x < numCols; x++) { - if (data1[y][x] & maskVal) { - numBad++; - continue; - } - if (data2 && data2[y][x] & maskVal) { - numBad++; - } - } - } - if (numBad > badFrac * numCols * numRows) { - psError(PS_ERR_BAD_PARAMETER_VALUE, true, - "Fraction of bad pixels (%d/%d=%f) exceeds limit (%f)\n", - numBad, numCols * numRows, (float)numBad/(float)(numCols * numRows), badFrac); - return NULL; - } - } - - // Worried about the masks for bad pixels and bad stamps colliding, so make our own mask - psImage *mask = psImageAlloc(numCols, numRows, PS_TYPE_MASK); // The global mask - psImageInit(mask, 0); - psMaskType **maskData = mask->data.PS_TYPE_MASK_DATA; // Dereference for convenience - - // Block out a border around the edge of the image - - // Bottom stripe - for (int y = 0; y < PS_MIN(size + footprint, numRows); y++) { - for (int x = 0; x < numCols; x++) { - maskData[y][x] |= PM_SUBTRACTION_MASK_BORDER; - } - } - // Either side - for (int y = PS_MIN(size + footprint, numRows); y < numRows - size - footprint; y++) { - for (int x = 0; x < PS_MIN(size + footprint, numCols); x++) { - maskData[y][x] |= PM_SUBTRACTION_MASK_BORDER; - } - for (int x = PS_MAX(numCols - size - footprint, 0); x < numCols; x++) { - maskData[y][x] |= PM_SUBTRACTION_MASK_BORDER; - } - } - // Top stripe - for (int y = PS_MAX(numRows - size - footprint, 0); y < numRows; y++) { - for (int x = 0; x < numCols; x++) { - maskData[y][x] |= PM_SUBTRACTION_MASK_BORDER; - } - } - - // XXX Could do something smarter here --- we will get images that are predominantly masked (where the - // skycell isn't overlapped by a large fraction by the observation), so that convolving around every bad - // pixel is wasting time. As a first cut, I've put in a check on the fraction of bad pixels, but we could - // imagine looking for the edge of big regions and convolving just at the edge. As a second cut, allow - // use of FFT convolution. - - for (int y = 0; y < numRows; y++) { - for (int x = 0; x < numCols; x++) { - if (data1[y][x] & maskVal) { - maskData[y][x] |= PM_SUBTRACTION_MASK_BAD_1; - } - if (data2 && data2[y][x] & maskVal) { - maskData[y][x] |= PM_SUBTRACTION_MASK_BAD_2; - } - } - } - - // Block out the entire stamp footprint around bad input pixels. - - // We want to block out with the CONVOLVE mask anything that would be bad if we convolved with a bad - // reference pixel (within 'size'). Then we want to block out with the FOOTPRINT mask everything within a - // footprint's distance of those (within 'footprint'). - - if (useFFT) { - if (!psImageConvolveMaskFFT(mask, mask, PM_SUBTRACTION_MASK_BAD_1 | PM_SUBTRACTION_MASK_BAD_2, - PM_SUBTRACTION_MASK_CONVOLVE_1 | PM_SUBTRACTION_MASK_FOOTPRINT_1, - -size, size, -size, size, 0.5)) { - psError(PS_ERR_UNKNOWN, false, "Unable to convolve bad pixels from mask 1."); - psFree(mask); - return NULL; - } - if (!psImageConvolveMaskFFT(mask, mask, PM_SUBTRACTION_MASK_BAD_2, - PM_SUBTRACTION_MASK_CONVOLVE_2 | PM_SUBTRACTION_MASK_FOOTPRINT_2, - -size, size, -size, size, 0.5)) { - psError(PS_ERR_UNKNOWN, false, "Unable to convolve bad pixels from mask 2."); - psFree(mask); - return NULL; - } - if (!psImageConvolveMaskFFT(mask, mask, PM_SUBTRACTION_MASK_CONVOLVE_1, - PM_SUBTRACTION_MASK_FOOTPRINT_1, - -footprint, footprint, -footprint, footprint, 0.5)) { - psError(PS_ERR_UNKNOWN, false, "Unable to reconvolve bad pixels from mask 1."); - psFree(mask); - return NULL; - } - if (!psImageConvolveMaskFFT(mask, mask, PM_SUBTRACTION_MASK_CONVOLVE_2, - PM_SUBTRACTION_MASK_FOOTPRINT_2, - -footprint, footprint, -footprint, footprint, 0.5)) { - psError(PS_ERR_UNKNOWN, false, "Unable to reconvolve bad pixels from mask 2."); - psFree(mask); - return NULL; - } - } else { - if (!psImageConvolveMaskDirect(mask, mask, PM_SUBTRACTION_MASK_BAD_1, - PM_SUBTRACTION_MASK_CONVOLVE_1 | PM_SUBTRACTION_MASK_FOOTPRINT_1, - -size, size, -size, size)) { - psError(PS_ERR_UNKNOWN, false, "Unable to convolve bad pixels from mask 1."); - psFree(mask); - return NULL; - } - if (!psImageConvolveMaskDirect(mask, mask, PM_SUBTRACTION_MASK_BAD_2, - PM_SUBTRACTION_MASK_CONVOLVE_2 | PM_SUBTRACTION_MASK_FOOTPRINT_2, - -size, size, -size, size)) { - psError(PS_ERR_UNKNOWN, false, "Unable to convolve bad pixels from mask 2."); - psFree(mask); - return NULL; - } - if (!psImageConvolveMaskDirect(mask, mask, PM_SUBTRACTION_MASK_CONVOLVE_1, - PM_SUBTRACTION_MASK_FOOTPRINT_1, - -footprint, footprint, -footprint, footprint)) { - psError(PS_ERR_UNKNOWN, false, "Unable to reconvolve bad pixels from mask 1."); - psFree(mask); - return NULL; - } - if (!psImageConvolveMaskDirect(mask, mask, PM_SUBTRACTION_MASK_CONVOLVE_2, - PM_SUBTRACTION_MASK_FOOTPRINT_2, - -footprint, footprint, -footprint, footprint)) { - psError(PS_ERR_UNKNOWN, false, "Unable to reconvolve bad pixels from mask 2."); - psFree(mask); - return NULL; - } - } - - return mask; -} // Convolve a stamp by a single kernel basis function @@ -538,4 +362,11 @@ ) { +#if 0 + if (index == kernels->num) { + // Normalisation component + return convolveOffset(image, 0, 0, footprint); + } +#endif + switch (kernels->type) { case PM_SUBTRACTION_KERNEL_POIS: { @@ -543,7 +374,5 @@ int v = kernels->v->data.S32[index]; // Offset in y psKernel *convolved = convolveOffset(image, u, v, footprint); // Convolved image - if (kernels->spatialOrder > 0 && index != kernels->subIndex) { - convolveSub(convolved, image, footprint); - } + convolveSub(convolved, image, footprint); return convolved; } @@ -569,7 +398,5 @@ } } - if (kernels->spatialOrder > 0 && index != kernels->subIndex) { - convolveSub(convolved, image, footprint); - } + convolveSub(convolved, image, footprint); return convolved; } @@ -583,7 +410,5 @@ int v = kernels->v->data.S32[index]; // Offset in y psKernel *convolved = convolveOffset(image, u, v, footprint); // Convolved image - if (kernels->spatialOrder > 0 && index != kernels->subIndex) { - convolveSub(convolved, image, footprint); - } + convolveSub(convolved, image, footprint); return convolved; } @@ -593,8 +418,4 @@ } case PM_SUBTRACTION_KERNEL_RINGS: { - if (index == kernels->subIndex) { - // Simply copying over the image - return convolveOffset(image, 0, 0, footprint); - } psKernel *convolved = psKernelAlloc(-footprint, footprint, -footprint, footprint); // Convolved image @@ -650,10 +471,8 @@ -bool pmSubtractionConvolveStamp(pmSubtractionStamp *stamp, const pmSubtractionKernels *kernels, int footprint, - pmSubtractionMode mode) +bool pmSubtractionConvolveStamp(pmSubtractionStamp *stamp, const pmSubtractionKernels *kernels, int footprint) { PS_ASSERT_PTR_NON_NULL(stamp, false); - PS_ASSERT_PTR_NON_NULL(kernels, false); - PS_ASSERT_VECTORS_SIZE_EQUAL(kernels->u, kernels->v, false); + PM_ASSERT_SUBTRACTION_KERNELS_NON_NULL(kernels, false); PS_ASSERT_INT_NONNEGATIVE(footprint, false); @@ -663,5 +482,5 @@ } - switch (mode) { + switch (kernels->mode) { case PM_SUBTRACTION_MODE_TARGET: case PM_SUBTRACTION_MODE_1: @@ -677,5 +496,5 @@ break; default: - psAbort("Unsupported subtraction mode: %x", mode); + psAbort("Unsupported subtraction mode: %x", kernels->mode); } @@ -684,382 +503,4 @@ -bool pmSubtractionCalculateEquation(pmSubtractionStampList *stamps, const pmSubtractionKernels *kernels, - int footprint, pmSubtractionMode mode) -{ - PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(stamps, false); - PS_ASSERT_PTR_NON_NULL(kernels, false); - PS_ASSERT_VECTORS_SIZE_EQUAL(kernels->u, kernels->v, false); - PS_ASSERT_INT_NONNEGATIVE(footprint, false); - - int spatialOrder = kernels->spatialOrder; // Maximum order of spatial variation - int numKernels = kernels->num; // Number of kernel basis functions - int numSpatial = polyTerms(spatialOrder); // Number of spatial variations - int numBackground = polyTerms(kernels->bgOrder); // Number of background terms - - // Total number of parameters to solve for: coefficient of each kernel basis function, multipled by the - // number of coefficients for the spatial polynomial, and a constant background offset. - int numParams = numKernels * numSpatial + numBackground; - int bgIndex = numParams - numBackground; // Index in matrix for the background - - // We iterate over each stamp, allocate the matrix and vectors if - // necessary, and then calculate those matrix/vectors. - 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; - } - - psImage *stampMatrix = stamp->matrix; // Least-squares matrix for this stamp - psVector *stampVector = stamp->vector; // Least-squares vector for this stamp - - if (!stampMatrix) { - stamp->matrix = stampMatrix = psImageAlloc(numParams, numParams, PS_TYPE_F64); - } - assert(stampMatrix->type.type == PS_TYPE_F64); - assert(stampMatrix->numCols == numParams && stampMatrix->numRows == numParams); - psImageInit(stampMatrix, 0.0); - - if (!stampVector) { - stamp->vector = stampVector = psVectorAlloc(numParams, PS_TYPE_F64); - } - assert(stampVector->type.type == PS_TYPE_F64); - assert(stampVector->n == numParams); - psVectorInit(stampVector, 0.0); - - // Generate convolutions - if (!pmSubtractionConvolveStamp(stamp, kernels, footprint, mode)) { - psError(PS_ERR_UNKNOWN, false, "Unable to convolve stamp %d.", i); - return NULL; - } - - psKernel *weight = stamp->weight; // Weight map postage stamp - psKernel *target; // Target postage stamp - psArray *convolutions; // Convolutions for each kernel - - switch (mode) { - case PM_SUBTRACTION_MODE_TARGET: - case PM_SUBTRACTION_MODE_1: - target = stamp->image2; - convolutions = stamp->convolutions1; - break; - case PM_SUBTRACTION_MODE_2: - target = stamp->image1; - convolutions = stamp->convolutions2; - break; - default: - psAbort("Unsupported subtraction mode: %x", mode); - } - - psImage *polyValues = spatialPolyValues(spatialOrder, stamp->xNorm, stamp->yNorm); // Polynomial terms - -#ifdef TESTING - for (int j = 0; j < numKernels; j++) { - psKernel *conv = convolutions->data[j]; // Convolution of interest - psString filename = NULL; // Filename for output - psStringAppend(&filename, "conv_%03d_%03d.fits", i, j); - psFits *fits = psFitsOpen(filename, "w"); // FITS file pointer - psFree(filename); - psFitsWriteImage(fits, NULL, conv->image, 0, NULL); - psFitsClose(fits); - } -#endif - - // Generate least-squares vector and matrix - bool bad = false; // Are there bad values? - psString badString = NULL; // Details on bad values - for (int j = 0; j < numKernels && !bad; j++) { - psKernel *jConv = convolutions->data[j]; // Convolution for j-th element - - // Generate matrix - for (int k = j; k < numKernels && !bad; k++) { - psKernel *kConv = convolutions->data[k]; // Convolution for k-th element - double sumCC = 0.0; // Sum of the convolution products - for (int y = - footprint; y <= footprint; y++) { - for (int x = - footprint; x <= footprint; x++) { - sumCC += jConv->kernel[y][x] * kConv->kernel[y][x] / weight->kernel[y][x]; - } - } - if (!isfinite(sumCC)) { - psStringAppend(&badString, "\nBad sumCC at %d, %d", j, k); - bad = true; - continue; - } - // Generate the pseudo-convolutions from the spatial polynomial terms - for (int jyOrder = 0, jIndex = j; jyOrder <= spatialOrder; jyOrder++) { - for (int jxOrder = 0; jxOrder <= spatialOrder - jyOrder; - jxOrder++, jIndex += numKernels) { - for (int kyOrder = 0, kIndex = k; kyOrder <= spatialOrder; kyOrder++) { - for (int kxOrder = 0; kxOrder <= spatialOrder - kyOrder; - kxOrder++, kIndex += numKernels) { - double convPoly = sumCC * polyValues->data.F64[jyOrder][jxOrder] * - polyValues->data.F64[kyOrder][kxOrder]; // Temporary value - stampMatrix->data.F64[jIndex][kIndex] = convPoly; - stampMatrix->data.F64[kIndex][jIndex] = convPoly; - } - } - } - } - } - - // Vector and background term for matrix - double sumC = 0.0; // Sum of the convolution - double sumTC = 0.0; // Sum of the convolution/target product - for (int y = - footprint; y <= footprint; y++) { - for (int x = - footprint; x <= footprint; x++) { - double convProduct = jConv->kernel[y][x] / weight->kernel[y][x]; // Convolution / noise^2 - sumC += convProduct; - sumTC += target->kernel[y][x] * convProduct; - } - } - if (!isfinite(sumC)) { - psStringAppend(&badString, "\nBad sumC at %d", j); - bad = true; - continue; - } - if (!isfinite(sumTC)) { - psStringAppend(&badString, "\nBad sumIC detected at %d", j); - bad = true; - continue; - } - - for (int jyOrder = 0, jIndex = j; jyOrder <= spatialOrder; jyOrder++) { - for (int jxOrder = 0; jxOrder <= spatialOrder - jyOrder; - jxOrder++, jIndex += numKernels) { - stampVector->data.F64[jIndex] = sumTC * polyValues->data.F64[jyOrder][jxOrder]; - - double convPoly = sumC * polyValues->data.F64[jyOrder][jxOrder]; // Value - stampMatrix->data.F64[jIndex][bgIndex] = convPoly; - stampMatrix->data.F64[bgIndex][jIndex] = convPoly; - } - } - - } - psFree(polyValues); - - // Background only terms - double sum1 = 0.0; // Sum of the weighting - double sumT = 0.0; // Sum of the target - for (int y = - footprint; y <= footprint; y++) { - for (int x = - footprint; x <= footprint; x++) { - double invNoise2 = 1.0 / weight->kernel[y][x]; // Inverse noise, squared - sum1 += invNoise2; - sumT += target->kernel[y][x] * invNoise2; - } - } - if (!isfinite(sum1)) { - psStringAppend(&badString, "\nBad sum1 detected"); - bad = true; - } else if (!isfinite(sumT)) { - psStringAppend(&badString, "\nBad sumI detected"); - bad = true; - } - - stampMatrix->data.F64[bgIndex][bgIndex] = sum1; - stampVector->data.F64[bgIndex] = sumT; - - if (bad) { - stamp->status = PM_SUBTRACTION_STAMP_REJECTED; - psWarning("Rejecting stamp %d (%d,%d) because of bad equation:%s", - i, (int)(stamp->x + 0.5), (int)(stamp->y + 0.5), badString); - psFree(badString); - } else { - stamp->status = PM_SUBTRACTION_STAMP_USED; - } - - if (psTraceGetLevel("psModules.imcombine.equation") >= 10) { - psString matrixName = NULL; - psStringAppend(&matrixName, "matrix%d.fits", i); - psFits *matrixFile = psFitsOpen(matrixName, "w"); - psFree(matrixName); - psFitsWriteImage(matrixFile, NULL, stampMatrix, 0, NULL); - psFitsClose(matrixFile); - } - - } - - return true; -} - - -psVector *pmSubtractionSolveEquation(psVector *solution, const pmSubtractionStampList *stamps) -{ - PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(stamps, NULL); - - // Check inputs, while summing the stamp matrices and vectors - long numParams = -1; // Number of parameters - for (int i = 0; i < stamps->num; i++) { - pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest - PS_ASSERT_PTR_NON_NULL(stamp, NULL); - if (stamp->status != PM_SUBTRACTION_STAMP_USED) { - continue; - } - if (numParams == -1) { - numParams = stamp->vector->n; - } else { - PS_ASSERT_VECTOR_SIZE(stamp->vector, numParams, NULL); - } - PS_ASSERT_VECTOR_TYPE(stamp->vector, PS_TYPE_F64, NULL); - PS_ASSERT_IMAGE_TYPE(stamp->matrix, PS_TYPE_F64, NULL); - PS_ASSERT_IMAGE_SIZE(stamp->matrix, numParams, numParams, NULL); - } - if (numParams == -1) { - psError(PS_ERR_BAD_PARAMETER_VALUE, true, "No suitable stamps found."); - return NULL; - } - - if (solution) { - PS_ASSERT_VECTOR_TYPE(solution, PS_TYPE_F64, NULL); - PS_ASSERT_VECTOR_SIZE(solution, numParams, NULL); - } - - // Accumulate the least-squares matricies and vectors - psImage *sumMatrix = psImageAlloc(numParams, numParams, PS_TYPE_F64); // Combined matrix for all stamps - psVector *sumVector = psVectorAlloc(numParams, PS_TYPE_F64); // Combined vector for all stamps - psVectorInit(sumVector, 0.0); - psImageInit(sumMatrix, 0.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) { - (void)psBinaryOp(sumMatrix, sumMatrix, "+", stamp->matrix); - (void)psBinaryOp(sumVector, sumVector, "+", stamp->vector); - } - } - - psVector *permutation = NULL; // Permutation vector, required for LU decomposition - psImage *luMatrix = psMatrixLUD(NULL, &permutation, sumMatrix); - psFree(sumMatrix); - if (!luMatrix) { - psError(PS_ERR_UNKNOWN, true, "LU Decomposition of least-squares matrix failed.\n"); - psFree(sumVector); - psFree(luMatrix); - psFree(permutation); - return NULL; - } - - solution = psMatrixLUSolve(solution, luMatrix, sumVector, permutation); - psFree(sumVector); - psFree(luMatrix); - psFree(permutation); - if (!solution) { - psError(PS_ERR_UNKNOWN, true, "Failed to solve the least-squares system.\n"); - return NULL; - } - - if (psTraceGetLevel("psModules.imcombine") >= 7) { - for (int i = 0; i < solution->n; i++) { - psTrace("psModules.imcombine", 7, "Solution %d: %f\n", i, solution->data.F64[i]); - } - } - - return solution; -} - -psVector *pmSubtractionCalculateDeviations(pmSubtractionStampList *stamps, const psVector *solution, - int footprint, const pmSubtractionKernels *kernels, - pmSubtractionMode mode) -{ - PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(stamps, NULL); - PS_ASSERT_VECTOR_NON_NULL(solution, NULL); - PS_ASSERT_VECTOR_TYPE(solution, PS_TYPE_F64, NULL); - PS_ASSERT_PTR_NON_NULL(kernels, NULL); - PS_ASSERT_VECTOR_SIZE(solution, kernels->num * polyTerms(kernels->spatialOrder) + - polyTerms(kernels->bgOrder), NULL); - PS_ASSERT_VECTOR_SIZE(kernels->u, kernels->num, NULL); - PS_ASSERT_VECTORS_SIZE_EQUAL(kernels->u, kernels->v, NULL); - - psVector *deviations = psVectorAlloc(stamps->num, PS_TYPE_F32); // Mean deviation for stamps - double devNorm = 1.0 / PS_SQR(2 * footprint + 1); // Normalisation for deviations - int numKernels = kernels->num; // Number of kernels - int spatialOrder = kernels->spatialOrder; // Order of kernel spatial variations - float background = solution->data.F64[solution->n-1]; // The difference in background - - psVector *coeff = psVectorAlloc(numKernels, PS_TYPE_F64); // Coefficients for the kernel basis functions - - 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) { - deviations->data.F32[i] = NAN; - continue; - } - - // Calculate coefficients of the kernel basis functions - psImage *polyValues = spatialPolyValues(kernels->spatialOrder, - stamp->xNorm, stamp->yNorm); // Polynomial terms - for (int j = 0; j < numKernels; j++) { - double polynomial = 0.0; // Value of the polynomial - for (int yOrder = 0, index = j; yOrder <= spatialOrder; yOrder++) { - for (int xOrder = 0; xOrder <= spatialOrder - yOrder; xOrder++, index += numKernels) { - polynomial += solution->data.F64[index] * polyValues->data.F64[yOrder][xOrder]; - } - } - coeff->data.F64[j] = polynomial; - } - psFree(polyValues); - - // Calculate residuals - psKernel *weight = stamp->weight; // Weight postage stamp - psKernel *target; // Target postage stamp - psArray *convolutions; // Convolution postage stamps for each kernel basis function - switch (mode) { - case PM_SUBTRACTION_MODE_TARGET: - case PM_SUBTRACTION_MODE_1: - target = stamp->image2; - convolutions = stamp->convolutions1; - break; - case PM_SUBTRACTION_MODE_2: - target = stamp->image1; - convolutions = stamp->convolutions2; - break; - default: - psAbort("Unsupported subtraction mode: %x", mode); - } - -#ifdef TESTING - psKernel *residual = psKernelAlloc(-footprint, footprint, -footprint, footprint); -#endif - float deviation = 0.0; // Deviation for this stamp - for (int y = - footprint; y <= footprint; y++) { - for (int x = - footprint; x <= footprint; x++) { - float conv = background; // The value of the convolution - for (int j = 0; j < numKernels; j++) { - psKernel *convolution = convolutions->data[j]; // Convolution - conv += convolution->kernel[y][x] * coeff->data.F64[j]; - } - float diff = target->kernel[y][x] - conv; - deviation += PS_SQR(diff) / weight->kernel[y][x]; -#ifdef TESTING - residual->kernel[y][x] = diff; -#endif - } - } - - deviations->data.F32[i] = sqrtf(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]); - 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)); - 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); - } - psFree(residual); -#endif - - } - psFree(coeff); - - return deviations; -} @@ -1148,22 +589,16 @@ } -psImage *pmSubtractionKernelImage(const psVector *solution, const pmSubtractionKernels *kernels, - float x, float y - ) -{ - PS_ASSERT_VECTOR_NON_NULL(solution, NULL); - PS_ASSERT_PTR_NON_NULL(kernels, NULL); - PS_ASSERT_VECTOR_SIZE(solution, kernels->num * polyTerms(kernels->spatialOrder) + - polyTerms(kernels->bgOrder), NULL); +psImage *pmSubtractionKernelImage(const pmSubtractionKernels *kernels, float x, float y, bool wantDual) +{ + PM_ASSERT_SUBTRACTION_KERNELS_NON_NULL(kernels, NULL); + PM_ASSERT_SUBTRACTION_KERNELS_SOLUTION(kernels, NULL); PS_ASSERT_FLOAT_WITHIN_RANGE(x, -1.0, 1.0, NULL); PS_ASSERT_FLOAT_WITHIN_RANGE(y, -1.0, 1.0, NULL); - PS_ASSERT_VECTOR_SIZE(kernels->u, kernels->num, false); - PS_ASSERT_VECTORS_SIZE_EQUAL(kernels->u, kernels->v, NULL); // Precalulate polynomial values - psImage *polyValues = spatialPolyValues(kernels->spatialOrder, x, y); + psImage *polyValues = p_pmSubtractionPolynomial(NULL, kernels->spatialOrder, x, y); // The appropriate kernel - psKernel *kernel = solvedKernel(NULL, solution, kernels, polyValues); + psKernel *kernel = solvedKernel(NULL, kernels, polyValues, wantDual); psFree(polyValues); @@ -1175,16 +610,11 @@ } -psArray *pmSubtractionKernelSolutions(const psVector *solution, const pmSubtractionKernels *kernels, - float x, float y - ) -{ - PS_ASSERT_VECTOR_NON_NULL(solution, NULL); - PS_ASSERT_PTR_NON_NULL(kernels, NULL); - PS_ASSERT_VECTOR_SIZE(solution, kernels->num * polyTerms(kernels->spatialOrder) + - polyTerms(kernels->bgOrder), NULL); +#if 0 +psArray *pmSubtractionKernelSolutions(const pmSubtractionKernels *kernels, float x, float y, bool wantDual) +{ + PM_ASSERT_SUBTRACTION_KERNELS_NON_NULL(kernels, NULL); + PM_ASSERT_SUBTRACTION_KERNELS_SOLUTION(kernels, NULL); PS_ASSERT_FLOAT_WITHIN_RANGE(x, -1.0, 1.0, NULL); PS_ASSERT_FLOAT_WITHIN_RANGE(y, -1.0, 1.0, NULL); - PS_ASSERT_VECTOR_SIZE(kernels->u, kernels->num, false); - PS_ASSERT_VECTORS_SIZE_EQUAL(kernels->u, kernels->v, NULL); psArray *images = psArrayAlloc(solution->n - 1); // Images of each kernel to return @@ -1194,5 +624,5 @@ for (int i = 0; i < solution->n - 1; i++) { fakeSolution->data.F64[i] = solution->data.F64[i]; - images->data[i] = pmSubtractionKernelImage(fakeSolution, kernels, x, y); + images->data[i] = pmSubtractionKernelImage(kernels, x, y, wantDual); fakeSolution->data.F64[i] = 0.0; } @@ -1202,11 +632,13 @@ return images; } - -bool pmSubtractionConvolve(pmReadout *out, const pmReadout *ro1, const pmReadout *ro2, +#endif + + + +bool pmSubtractionConvolve(pmReadout *out1, pmReadout *out2, const pmReadout *ro1, const pmReadout *ro2, const psImage *subMask, psMaskType blank, const psRegion *region, - const psVector *solution, const pmSubtractionKernels *kernels, - pmSubtractionMode mode, bool useFFT) -{ - PS_ASSERT_PTR_NON_NULL(out, false); + const pmSubtractionKernels *kernels, bool useFFT) +{ + PS_ASSERT_PTR_NON_NULL(out1, false); PS_ASSERT_PTR_NON_NULL(ro1, false); PS_ASSERT_IMAGE_NON_NULL(ro1->image, false); @@ -1217,5 +649,8 @@ PS_ASSERT_IMAGES_SIZE_EQUAL(ro1->weight, ro1->image, false); } - if (mode != PM_SUBTRACTION_MODE_1) { + PM_ASSERT_SUBTRACTION_KERNELS_NON_NULL(kernels, false); + PM_ASSERT_SUBTRACTION_KERNELS_SOLUTION(kernels, false); + if (kernels->mode != PM_SUBTRACTION_MODE_1) { + PS_ASSERT_PTR_NON_NULL(out2, false); PS_ASSERT_PTR_NON_NULL(ro2, false); PS_ASSERT_IMAGE_NON_NULL(ro2->image, false); @@ -1232,27 +667,20 @@ PS_ASSERT_IMAGES_SIZE_EQUAL(subMask, ro1->image, false); } - if (out->image) { - PS_ASSERT_IMAGE_NON_NULL(out->image, false); - PS_ASSERT_IMAGES_SIZE_EQUAL(out->image, ro1->image, false); - PS_ASSERT_IMAGE_TYPE(out->image, PS_TYPE_F32, false); - } - if (out->mask) { - PS_ASSERT_IMAGE_NON_NULL(out->mask, false); - PS_ASSERT_IMAGES_SIZE_EQUAL(out->mask, ro1->image, false); - PS_ASSERT_IMAGE_TYPE(out->mask, PS_TYPE_MASK, false); + if (out1->image) { + PS_ASSERT_IMAGE_NON_NULL(out1->image, false); + PS_ASSERT_IMAGES_SIZE_EQUAL(out1->image, ro1->image, false); + PS_ASSERT_IMAGE_TYPE(out1->image, PS_TYPE_F32, false); + } + if (out1->mask) { + PS_ASSERT_IMAGE_NON_NULL(out1->mask, false); + PS_ASSERT_IMAGES_SIZE_EQUAL(out1->mask, ro1->image, false); + PS_ASSERT_IMAGE_TYPE(out1->mask, PS_TYPE_MASK, false); PS_ASSERT_IMAGE_NON_NULL(subMask, false); } - if (out->weight) { - PS_ASSERT_IMAGE_NON_NULL(out->weight, false); - PS_ASSERT_IMAGES_SIZE_EQUAL(out->weight, ro1->image, false); - PS_ASSERT_IMAGE_TYPE(out->weight, PS_TYPE_F32, false); - } - PS_ASSERT_VECTOR_NON_NULL(solution, false); - PS_ASSERT_VECTOR_TYPE(solution, PS_TYPE_F64, false); - PS_ASSERT_PTR_NON_NULL(kernels, false); - PS_ASSERT_VECTOR_SIZE(solution, kernels->num * polyTerms(kernels->spatialOrder) + - polyTerms(kernels->bgOrder), -1); - PS_ASSERT_VECTOR_SIZE(kernels->u, kernels->num, false); - PS_ASSERT_VECTORS_SIZE_EQUAL(kernels->u, kernels->v, false); + if (out1->weight) { + PS_ASSERT_IMAGE_NON_NULL(out1->weight, false); + PS_ASSERT_IMAGES_SIZE_EQUAL(out1->weight, ro1->image, false); + PS_ASSERT_IMAGE_TYPE(out1->weight, PS_TYPE_F32, false); + } if (region) { if (psRegionIsNaN(*region)) { @@ -1273,7 +701,8 @@ psImage *image, *weight; // Image and weight map to convolve - switch (mode) { + switch (kernels->mode) { case PM_SUBTRACTION_MODE_TARGET: case PM_SUBTRACTION_MODE_1: + case PM_SUBTRACTION_MODE_DUAL: image = ro1->image; weight = ro1->weight; @@ -1281,42 +710,56 @@ case PM_SUBTRACTION_MODE_2: image = ro2->image; - weight = ro2->image; + weight = ro2->weight; break; default: - psAbort("Unsupported subtraction mode: %x", mode); - } - - int numBackground = polyTerms(kernels->bgOrder); // Number of background terms - float background = solution->data.F64[solution->n - numBackground]; // The difference in background + psAbort("Unsupported subtraction mode: %x", kernels->mode); + } + int numCols = ro1->image->numCols, numRows = ro1->image->numRows; // Image dimensions - psImage *convImage = out->image; // Convolved image - if (!convImage) { - out->image = psImageAlloc(numCols, numRows, PS_TYPE_F32); - convImage = out->image; + psImage *convImage1 = out1->image; // Convolved image + if (!convImage1) { + convImage1 = out1->image = psImageAlloc(numCols, numRows, PS_TYPE_F32); } psImage *convMask = NULL; // Convolved mask image if (subMask) { - if (!out->mask) { - out->mask = psImageAlloc(numCols, numRows, PS_TYPE_MASK); - } - convMask = out->mask; + if (!out1->mask) { + out1->mask = psImageAlloc(numCols, numRows, PS_TYPE_MASK); + } + convMask = out1->mask; psImageInit(convMask, 0); } - psImage *convWeight = NULL; // Convolved weight (variance) image + psImage *convWeight1 = NULL; // Convolved weight (variance) image if (weight) { - if (!out->weight) { - out->weight = psImageAlloc(numCols, numRows, PS_TYPE_F32); - } - convWeight = out->weight; - psImageInit(convWeight, 0.0); + if (!out1->weight) { + out1->weight = psImageAlloc(numCols, numRows, PS_TYPE_F32); + } + convWeight1 = out1->weight; + psImageInit(convWeight1, 0.0); + } + + psImage *convImage2, *convWeight2; // Convolved products for dual mode + if (kernels->mode == PM_SUBTRACTION_MODE_DUAL) { + convImage2 = out2->image; + if (!convImage2) { + convImage2 = out2->image = psImageAlloc(numCols, numRows, PS_TYPE_F32); + } + convWeight2 = NULL; // Convolved weight (variance) image + if (ro2->weight) { + if (!out2->weight) { + out2->weight = psImageAlloc(numCols, numRows, PS_TYPE_F32); + } + convWeight2 = out2->weight; + psImageInit(convWeight2, 0.0); + } + if (subMask) { + // Copying mask --- they're the same! + psFree(out2->mask); + out2->mask = psMemIncrRefCounter(convMask); + } } int size = kernels->size; // Half-size of kernel int fullSize = 2 * size + 1; // Full size of kernel - psImage *polyValues = NULL; // Pre-calculated polynomial values - - psKernel *kernelImage = NULL; // Kernel for the image - psKernel *kernelWeight = NULL; // Kernel for the weight map // Get region for convolution: [xMin:xMax,yMin:yMax] @@ -1330,6 +773,7 @@ } - psMaskType maskSource, maskTarget; // Mask values for source and target - switch (mode) { + psMaskType maskSource; // Mask these pixels when convolving + psMaskType maskTarget; // Mark these pixels as bad when propagating the subtractionMask + switch (kernels->mode) { case PM_SUBTRACTION_MODE_TARGET: case PM_SUBTRACTION_MODE_1: @@ -1341,40 +785,36 @@ maskTarget = PM_SUBTRACTION_MASK_BAD_1 | PM_SUBTRACTION_MASK_CONVOLVE_2; break; + case PM_SUBTRACTION_MODE_DUAL: + maskSource = PM_SUBTRACTION_MASK_BAD_1 | PM_SUBTRACTION_MASK_BAD_2; + maskTarget = PM_SUBTRACTION_MASK_CONVOLVE_1 | PM_SUBTRACTION_MASK_CONVOLVE_2; + break; default: - psAbort("Unsupported subtraction mode: %x", mode); - } + psAbort("Unsupported subtraction mode: %x", kernels->mode); + } + + psImage *polyValues = NULL; // Pre-calculated polynomial values + psKernel *kernelImage = NULL; // Kernel for the images + psKernel *kernelWeight = NULL; // Kernel for the weight maps for (int j = yMin; j < yMax; j += fullSize) { int ySubMax = PS_MIN(j + fullSize, yMax); // Range for subregion of interest + float yNorm = 2.0 * (float)(j + size + 1 - numRows/2.0) / (float)numRows; // Normalised coordinate for (int i = xMin; i < xMax; i += fullSize) { int xSubMax = PS_MIN(i + fullSize, xMax); // Range for subregion of interest + float xNorm = 2.0 * (float)(i + size + 1 - numCols/2.0) / (float)numCols; // Normalised coordinate // Only generate polynomial values every kernel footprint, since we have already assumed // (with the stamps) that it does not vary rapidly on this scale. - psFree(polyValues); - polyValues = spatialPolyValues(kernels->spatialOrder, - 2.0 * (float)(i + size + 1 - numCols/2.0) / (float)numCols, - 2.0 * (float)(j + size + 1 - numRows/2.0) / (float)numRows); - - kernelImage = solvedKernel(kernelImage, solution, kernels, polyValues); - if (weight) { - kernelWeight = varianceKernel(kernelWeight, kernelImage); - } + polyValues = p_pmSubtractionPolynomial(polyValues, kernels->spatialOrder, xNorm, yNorm); + float background = p_pmSubtractionSolutionBackground(kernels, polyValues); // Background term psRegion subRegion = psRegionSet(i, xSubMax, j, ySubMax); // Sub-region to convolve - if (useFFT) { - // Use Fast Fourier Transform to do the convolution - // This provides a big speed-up for large kernels - convolveFFT(convImage, image, subMask, maskSource, kernelImage, subRegion, - background, size); - if (weight) { - convolveFFT(convWeight, weight, subMask, maskSource, kernelWeight, subRegion, - 0.0, size); - } - } else { - convolveDirect(convImage, image, kernelImage, subRegion, background, size); - if (weight) { - convolveDirect(convWeight, weight, kernelWeight, subRegion, 0.0, size); - } + + convolveRegion(convImage1, convWeight1, &kernelImage, &kernelWeight, image, weight, + subMask, maskSource, kernels, polyValues, background, subRegion, useFFT); + + if (kernels->mode == PM_SUBTRACTION_MODE_DUAL) { + convolveRegion(convImage2, convWeight2, &kernelImage, &kernelWeight, ro2->image, ro2->weight, + subMask, maskSource, kernels, polyValues, background, subRegion, useFFT); } @@ -1385,7 +825,7 @@ if (subMask->data.PS_TYPE_MASK_DATA[y][x] & maskTarget) { convMask->data.PS_TYPE_MASK_DATA[y][x] |= blank; - convImage->data.F32[y][x] = NAN; + convImage1->data.F32[y][x] = NAN; if (weight) { - convWeight->data.F32[y][x] = NAN; + convWeight1->data.F32[y][x] = NAN; } } @@ -1402,44 +842,2 @@ return true; } - -bool pmSubtractionBorder(psImage *image, psImage *weight, psImage *mask, - int size, psMaskType blank) -{ - PS_ASSERT_IMAGE_NON_NULL(image, false); - PS_ASSERT_IMAGE_TYPE(image, PS_TYPE_F32, false); - if (mask) { - PS_ASSERT_IMAGE_NON_NULL(mask, false); - PS_ASSERT_IMAGES_SIZE_EQUAL(mask, image, false); - PS_ASSERT_IMAGE_TYPE(mask, PS_TYPE_MASK, false); - } - if (weight) { - PS_ASSERT_IMAGE_NON_NULL(weight, false); - PS_ASSERT_IMAGES_SIZE_EQUAL(weight, image, false); - PS_ASSERT_IMAGE_TYPE(weight, PS_TYPE_F32, false); - } - - int numCols = image->numCols, numRows = image->numRows; // Image dimensions - - for (int y = size; y < numRows - size; y++) { - for (int x = 0; x < size; x++) { - markBlank(image, mask, weight, x, y, blank); - } - for (int x = numCols - size; x < numCols; x++) { - markBlank(image, mask, weight, x, y, blank); - } - } - for (int y = 0; y < size; y++) { - for (int x = 0; x < numCols; x++) { - markBlank(image, mask, weight, x, y, blank); - } - } - for (int y = numRows - size; y < numRows; y++) { - for (int x = 0; x < numCols; x++) { - markBlank(image, mask, weight, x, y, blank); - } - } - - return true; -} - - Index: trunk/psModules/src/imcombine/pmSubtraction.h =================================================================== --- trunk/psModules/src/imcombine/pmSubtraction.h (revision 15745) +++ trunk/psModules/src/imcombine/pmSubtraction.h (revision 15756) @@ -6,6 +6,6 @@ * @author GLG, MHPCC * - * @version $Revision: 1.21 $ $Name: not supported by cvs2svn $ - * @date $Date: 2007-11-10 01:09:20 $ + * @version $Revision: 1.22 $ $Name: not supported by cvs2svn $ + * @date $Date: 2007-12-07 01:57:15 $ * Copyright 2004-207 Institute for Astronomy, University of Hawaii */ @@ -37,39 +37,28 @@ } pmSubtractionMasks; -/// Generate a mask for use in the subtraction process -psImage *pmSubtractionMask(const psImage *refMask, ///< Mask for the reference image (will be convolved) - const psImage *inMask, ///< Mask for the input image, or NULL - psMaskType maskVal, ///< Value to mask out - int size, ///< Half-size of the kernel (pmSubtractionKernels.size) - int footprint, ///< Half-size of the kernel footprint - float badFrac, ///< Maximum fraction of bad input pixels to accept - bool useFFT ///< Use FFT to do convolution? - ); + +/// Number of terms in a polynomial +#define PM_SUBTRACTION_POLYTERMS(ORDER) (((ORDER) + 1) * ((ORDER) + 2) / 2) + +/// Set the indices for the normalisation and background terms +#define PM_SUBTRACTION_INDICES(NORM,BG,KERNELS) { \ + int numSpatial = PM_SUBTRACTION_POLYTERMS((KERNELS)->spatialOrder); /* Number of spatial terms */ \ + NORM = (KERNELS)->num * numSpatial; \ + BG = NORM + 1; \ +} + +/// Return the index for the start of the normalisation terms +#define PM_SUBTRACTION_INDEX_NORM(KERNELS) \ + ((KERNELS)->num * PM_SUBTRACTION_POLYTERMS((KERNELS)->spatialOrder)) + +/// Return the index for the start of the background terms +#define PM_SUBTRACTION_INDEX_BG(KERNELS) \ + (((KERNELS)->num * PM_SUBTRACTION_POLYTERMS((KERNELS)->spatialOrder)) + 1) + /// Convolve the reference stamp with the kernel components bool pmSubtractionConvolveStamp(pmSubtractionStamp *stamp, ///< Stamp to convolve const pmSubtractionKernels *kernels, ///< Kernel parameters - int footprint, ///< Half-size of region over which to calculate equation - pmSubtractionMode mode ///< Mode for subtraction (which to convolve) - ); - -/// Calculate the least-squares equation to match the image quality -bool pmSubtractionCalculateEquation(pmSubtractionStampList *stamps, ///< Stamps - const pmSubtractionKernels *kernels, ///< Kernel parameters - int footprint, ///< Half-size of region over which to calculate equation - pmSubtractionMode mode ///< Mode for subtraction (which to convolve) - ); - -/// Solve the least-squares equation to match the image quality -psVector *pmSubtractionSolveEquation(psVector *solution, ///< Solution vector, or NULL - const pmSubtractionStampList *stamps ///< Stamps - ); - -/// Calculate deviations -psVector *pmSubtractionCalculateDeviations(pmSubtractionStampList *stamps, ///< Stamps - const psVector *solution, ///< Solution vector - int footprint, ///< Half-size of stamp - const pmSubtractionKernels *kernels, ///< Kernel parameters - pmSubtractionMode mode ///< Mode for subtraction + int footprint ///< Half-size of region over which to calculate equation ); @@ -83,7 +72,7 @@ /// Generate an image of the convolution kernel -psImage *pmSubtractionKernelImage(const psVector *solution, ///< Solution vector - const pmSubtractionKernels *kernels, ///< Kernel parameters - float x, float y ///< Normalised position [-1,1] for which to generate image +psImage *pmSubtractionKernelImage(const pmSubtractionKernels *kernels, ///< Kernel parameters + float x, float y,///< Normalised position [-1,1] for which to generate image + bool wantDual ///< Calculate for the dual kernel? ); @@ -95,5 +84,6 @@ /// Convolve image in preparation for subtraction -bool pmSubtractionConvolve(pmReadout *out, ///< Output image +bool pmSubtractionConvolve(pmReadout *out1, ///< Output image 1 + pmReadout *out2, ///< Output image 2 (DUAL mode only) const pmReadout *ro1, // Input image 1 const pmReadout *ro2, // Input image 2 @@ -101,16 +91,6 @@ psMaskType blank, ///< Mask value for blank regions const psRegion *region, ///< Region to convolve (or NULL) - const psVector *solution, ///< The solution vector const pmSubtractionKernels *kernels, ///< Kernel parameters - pmSubtractionMode mode, ///< Mode for subtraction bool useFFT ///< Use Fast Fourier Transform for the convolution? - ); - -/// Mark the non-convolved part of the image as blank -bool pmSubtractionBorder(psImage *image,///< Image - psImage *weight, ///< Weight map (or NULL) - psImage *mask, ///< Mask (or NULL) - int size, ///< Kernel half-size - psMaskType blank ///< Mask value for blank regions ); @@ -122,4 +102,10 @@ ); +/// Given (normalised) coordinates (x,y), generate a matrix where the elements (i,j) are x^i * y^j +psImage *p_pmSubtractionPolynomial(psImage *output, ///< Output matrix, or NULL + int spatialOrder, ///< Maximum spatial polynomial order + float x, float y ///< Normalised position of interest, [-1,1] + ); + /// @} #endif Index: trunk/psModules/src/imcombine/pmSubtractionEquation.c =================================================================== --- trunk/psModules/src/imcombine/pmSubtractionEquation.c (revision 15756) +++ trunk/psModules/src/imcombine/pmSubtractionEquation.c (revision 15756) @@ -0,0 +1,936 @@ +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include +#include + +#include "pmSubtraction.h" +#include "pmSubtractionKernels.h" +#include "pmSubtractionStamps.h" + +#include "pmSubtractionEquation.h" + + +#define TESTING + +////////////////////////////////////////////////////////////////////////////////////////////////////////////// +// Private (file-static) functions +////////////////////////////////////////////////////////////////////////////////////////////////////////////// + +// Calculate the sum over a stamp product +static inline double calculateSumProduct(const psKernel *image1, // First image in multiplication + const psKernel *image2, // Second image in multiplication + const psKernel *weight, // Weight image + int footprint // (Half-)Size of stamp + ) +{ + double sum = 0.0; // Sum of the image products + for (int y = - footprint; y <= footprint; y++) { + for (int x = - footprint; x <= footprint; x++) { + sum += image1->kernel[y][x] * image2->kernel[y][x] / weight->kernel[y][x]; + } + } + return sum; +} + +// Calculate a single element of the least-squares matrix, with the polynomial expansions in one direction +static inline bool calculateMatrixElement1(psImage *matrix, // Matrix to calculate + int i, int j, // Coordinates of element + const psKernel *image1, // First image in multiplication + const psKernel *image2, // Second image in multiplication + const psKernel *weight, // Weight image + const psImage *polyValues, // Spatial polynomial values + int numKernels, // Number of kernel basis functions + int footprint, // (Half-)Size of stamp + int spatialOrder, // Maximum order of spatial variation + bool symmetric // Is the matrix symmetric? + ) +{ + double sum = calculateSumProduct(image1, image2, weight, footprint); // Sum of the image products + if (!isfinite(sum)) { + return false; + } + + // Generate the pseudo-convolutions from the spatial polynomial terms + for (int iyOrder = 0, iIndex = i; iyOrder <= spatialOrder; iyOrder++) { + for (int ixOrder = 0; ixOrder <= spatialOrder - iyOrder; ixOrder++, iIndex += numKernels) { + double convPoly = sum * polyValues->data.F64[iyOrder][ixOrder]; + // XXX Not sure I've got this the right way around + matrix->data.F64[iIndex][j] = convPoly; + if (symmetric) { + matrix->data.F64[j][iIndex] = convPoly; + } + } + } + return true; +} + +// Calculate a single element of the least-squares matrix, with the polynomial expansions in both directions +static inline bool calculateMatrixElement2(psImage *matrix, // Matrix to calculate + int i, int j, // Coordinates of element + const psKernel *image1, // First image in multiplication + const psKernel *image2, // Second image in multiplication + const psKernel *weight, // Weight image + const psImage *polyValues, // Spatial polynomial values + int numKernels, // Number of kernel basis functions + int footprint, // (Half-)Size of stamp + int spatialOrder, // Maximum order of spatial variation + bool symmetric // Is the matrix symmetric? + ) +{ + double sum = calculateSumProduct(image1, image2, weight, footprint); // Sum of the image products + if (!isfinite(sum)) { + return false; + } + + // Generate the pseudo-convolutions from the spatial polynomial terms + for (int iyOrder = 0, iIndex = i; iyOrder <= spatialOrder; iyOrder++) { + for (int ixOrder = 0; ixOrder <= spatialOrder - iyOrder; ixOrder++, iIndex += numKernels) { + double iPoly = polyValues->data.F64[iyOrder][ixOrder]; // Value of polynomial + for (int jyOrder = 0, jIndex = j; jyOrder <= spatialOrder; jyOrder++) { + for (int jxOrder = 0; jxOrder <= spatialOrder - jyOrder; jxOrder++, jIndex += numKernels) { + double convPoly = sum * iPoly * polyValues->data.F64[jyOrder][jxOrder]; + matrix->data.F64[iIndex][jIndex] = convPoly; + if (symmetric) { + matrix->data.F64[jIndex][iIndex] = convPoly; + } + } + } + } + } + return true; +} + +// Calculate the square part of the matrix derived from multiplying convolutions +static bool calculateMatrixSquare(psImage *matrix, // Matrix to calculate + const psArray *convolutions1, // Convolutions for element 1 + const psArray *convolutions2, // Convolutions for element 2 + const psKernel *weight, // Weight image + const psImage *polyValues, // Polynomial values + int numKernels, // Number of kernel basis functions + int spatialOrder, // Order of spatial variation + int footprint // Half-size of stamp + ) +{ + bool symmetric = (convolutions1 == convolutions2 ? true : false); // Is matrix symmetric? + + for (int i = 0; i < numKernels; i++) { + psKernel *iConv = convolutions1->data[i]; // Convolution for i-th element + + for (int j = (symmetric ? i : 0); j < numKernels; j++) { + psKernel *jConv = convolutions2->data[j]; // Convolution for j-th element + + if (!calculateMatrixElement2(matrix, i, j, iConv, jConv, weight, polyValues, numKernels, + footprint, spatialOrder, symmetric)) { + psTrace("psModules.imcombine", 2, "Bad sumCC at %d, %d", i, j); + return false; + } + } + } + + return true; +} + +// Calculate least-squares matrix and vector +static bool calculateMatrix(psImage *matrix, // Matrix to calculate + const pmSubtractionKernels *kernels, // Kernel components + const psArray *convolutions, // Convolutions of source with kernels + const psKernel *input, // Input stamp, or NULL + const psKernel *target, // Target stamp, or NULL + const psKernel *weight, // Weight stamp + const psImage *polyValues, // Spatial polynomial values + int footprint, // (Half-)Size of stamp + bool normAndBG // Calculate normalisation and background terms? + ) +{ + int numKernels = kernels->num; // Number of kernel components + int spatialOrder = kernels->spatialOrder; // Maximum order of spatial variation + int numSpatial = PM_SUBTRACTION_POLYTERMS(spatialOrder); // Number of spatial variation terms + int bgOrder = kernels->bgOrder; // Maximum order of background fit + int numBackground = normAndBG ? PM_SUBTRACTION_POLYTERMS(bgOrder) : 0; // Number of background terms + int numTerms = numKernels * numSpatial + (normAndBG ? 1 + numBackground : 0); // Total number of terms + assert(matrix); + assert(matrix->numCols == matrix->numRows); + assert(matrix->numCols == numTerms); + assert(convolutions && convolutions->n == numKernels); + assert(target); + assert(polyValues); + assert(!normAndBG || input); // If we want the normalisation and BG, then we need the input image + + // Square part of the matrix (convolution-convolution products) + if (!calculateMatrixSquare(matrix, convolutions, convolutions, weight, polyValues, numKernels, + spatialOrder, footprint)) { + return false; + } + + // XXX To support higher-order background model than simply constant, the below code needs to be updated. + if (normAndBG) { + // Normalisation and background terms + int normIndex, bgIndex; // Indices in matrix for normalisation and background + PM_SUBTRACTION_INDICES(normIndex, bgIndex, kernels); + + for (int i = 0; i < numKernels; i++) { + psKernel *conv = convolutions->data[i]; // Convolution for i-th element + // Normalisation + if (!calculateMatrixElement1(matrix, i, normIndex, conv, input, weight, polyValues, numKernels, + footprint, spatialOrder, true)) { + psTrace("psModules.imcombine", 2, "Bad sumIC at %d", i); + return false; + } + + // Background + double sumC = 0.0; // Sum of the convolution + for (int y = - footprint; y <= footprint; y++) { + for (int x = - footprint; x <= footprint; x++) { + sumC += conv->kernel[y][x] / weight->kernel[y][x]; + } + } + if (!isfinite(sumC)) { + psTrace("psModules.imcombine", 2, "Bad sumC at %d", i); + return false; + } + for (int yOrder = 0, index = i; yOrder <= spatialOrder; yOrder++) { + for (int xOrder = 0; xOrder <= spatialOrder - yOrder; xOrder++, index += numKernels) { + matrix->data.F64[index][bgIndex] = matrix->data.F64[bgIndex][index] = + sumC * polyValues->data.F64[yOrder][xOrder]; + } + } + } + + // Normalisation only and background only terms in the matrix + double sum1 = 0.0; // Sum of the weighting + double sumI = 0.0; // Sum of the input + double sumII = 0.0; // Sum of the input squared + for (int y = - footprint; y <= footprint; y++) { + for (int x = - footprint; x <= footprint; x++) { + float invNoise2 = 1.0 / weight->kernel[y][x]; + float value = input->kernel[y][x] * invNoise2; + sumI += value; + sumII += input->kernel[y][x] * value; + sum1 += invNoise2; + } + } + if (!isfinite(sumII)) { + psTrace("psModules.imcombine", 2, "Bad sumII detected"); + return false; + } + if (!isfinite(sumI)) { + psTrace("psModules.imcombine", 2, "Bad sumI detected"); + return false; + } + if (!isfinite(sum1)) { + psTrace("psModules.imcombine", 2, "Bad sum1 detected"); + return false; + } + matrix->data.F64[normIndex][normIndex] = sumII; + matrix->data.F64[bgIndex][bgIndex] = sum1; + matrix->data.F64[normIndex][bgIndex] = sumI; + matrix->data.F64[bgIndex][normIndex] = sumI; + } + + return true; +} + + +// Calculate least-squares matrix and vector +static bool calculateVector(psVector *vector, // Vector to calculate, or NULL + const pmSubtractionKernels *kernels, // Kernel components + const psArray *convolutions, // Convolutions of source with kernels + const psKernel *input, // Input stamp, or NULL + const psKernel *target, // Target stamp, or NULL + const psKernel *weight, // Weight stamp + const psImage *polyValues, // Spatial polynomial values + int footprint // (Half-)Size of stamp + ) +{ + int numKernels = kernels->num; // Number of kernel components + int spatialOrder = kernels->spatialOrder; // Maximum order of spatial variation + int numSpatial = PM_SUBTRACTION_POLYTERMS(spatialOrder); // Number of spatial variation terms + int bgOrder = kernels->bgOrder; // Maximum order of background fit + int numBackground = PM_SUBTRACTION_POLYTERMS(bgOrder); // Number of background terms + int numTerms = numKernels * numSpatial + 1 + numBackground; // Total number of terms + assert(vector && vector->n == numTerms); + assert(convolutions && convolutions->n == numKernels); + assert(input); + assert(target); + assert(polyValues); + + // Convolution terms + for (int i = 0; i < numKernels; i++) { + psKernel *conv = convolutions->data[i]; // Convolution for i-th element + double sumTC = 0.0; // Sum of the target and convolution + for (int y = - footprint; y <= footprint; y++) { + for (int x = - footprint; x <= footprint; x++) { + sumTC += target->kernel[y][x] * conv->kernel[y][x] / weight->kernel[y][x]; + } + } + if (!isfinite(sumTC)) { + psTrace("psModules.imcombine", 2, "Bad sumTC at %d", i); + return false; + } + for (int yOrder = 0, index = i; yOrder <= spatialOrder; yOrder++) { + for (int xOrder = 0; xOrder <= spatialOrder - yOrder; xOrder++, index += numKernels) { + vector->data.F64[index] = sumTC * polyValues->data.F64[yOrder][xOrder]; + } + } + } + + // Normalisation and background terms + double sumIT = 0.0; // Sum of the target and input + double sumT = 0.0; // Sum of the target + for (int y = - footprint; y <= footprint; y++) { + for (int x = - footprint; x <= footprint; x++) { + float value = target->kernel[y][x] / weight->kernel[y][x]; + sumIT += value * input->kernel[y][x]; + sumT += value; + } + } + if (!isfinite(sumIT)) { + psTrace("psModules.imcombine", 2, "Bad sumIT detected"); + return false; + } + if (!isfinite(sumT)) { + psTrace("psModules.imcombine", 2, "Bad sumI detected"); + return false; + } + + int normIndex, bgIndex; // Indices in vector for normalisation and background terms + PM_SUBTRACTION_INDICES(normIndex, bgIndex, kernels); + vector->data.F64[normIndex] = sumIT; + vector->data.F64[bgIndex] = sumT; + + return true; +} + + + +// Calculate the cross-matrix, composed of convolutions of each image +// Note that the cross-matrix is NOT square +static bool calculateMatrixCross(psImage *matrix, // Matrix to calculate + const pmSubtractionKernels *kernels, // Kernel components + const psArray *convolutions1, // Convolutions of image 1 + const psArray *convolutions2, // Convolutions of image 2 + const psKernel *image1, // Image 1 stamp + const psKernel *weight, // Weight stamp + const psImage *polyValues, // Spatial polynomial values + int footprint // (Half-)Size of stamp + ) +{ + assert(matrix); + int numKernels = kernels->num; // Number of kernel components + int spatialOrder = kernels->spatialOrder; // Maximum order of spatial variation + int numSpatial = PM_SUBTRACTION_POLYTERMS(spatialOrder); // Number of spatial polynomial terms + int numBackground = PM_SUBTRACTION_POLYTERMS(kernels->bgOrder); // Number of background terms + int numCols = numKernels * numSpatial + 1 + numBackground; // Number of columns + int numRows = numKernels * numSpatial; // Number of rows + assert(matrix->numCols == numCols && matrix->numRows == numRows); + assert(convolutions1 && convolutions1->n == numKernels); + assert(convolutions2 && convolutions2->n == numKernels); + + int normIndex, bgIndex; // Indices in matrix for normalisation and background terms + PM_SUBTRACTION_INDICES(normIndex, bgIndex, kernels); + + if (!calculateMatrixSquare(matrix, convolutions1, convolutions2, weight, polyValues, numKernels, + spatialOrder, footprint)) { + return false; + } + + for (int i = 0; i < numKernels; i++) { + // Normalisation + psKernel *conv = convolutions2->data[i]; // Convolution + if (!calculateMatrixElement1(matrix, i, normIndex, conv, image1, weight, polyValues, numKernels, + footprint, spatialOrder, false)) { + psTrace("psModules.imcombine", 2, "Bad sumIC at %d", i); + return false; + } + + // Background + double sumC = 0.0; // Sum of the weighting + for (int y = - footprint; y <= footprint; y++) { + for (int x = - footprint; x <= footprint; x++) { + sumC += conv->kernel[y][x] / weight->kernel[y][x]; + } + } + if (!isfinite(sumC)) { + psTrace("psModules.imcombine", 2, "Bad sumC detected at %d", i); + return false; + } + + matrix->data.F64[i][bgIndex] = sumC; + } + + return true; +} + +////////////////////////////////////////////////////////////////////////////////////////////////////////////// +// Semi-public functions +////////////////////////////////////////////////////////////////////////////////////////////////////////////// + +// Calculate the value of a polynomial, specified by coefficients and polynomial values +inline double p_pmSubtractionCalculatePolynomial(const psVector *coeff, // Coefficients + const psImage *polyValues, // Polynomial values + int order, // Order of polynomials + int index, // Index at which to begin + int step // Step between subsequent indices + ) +{ + double sum = 0.0; // Value of the polynomial sum + for (int yOrder = 0; yOrder <= order; yOrder++) { + for (int xOrder = 0; xOrder <= order - yOrder; xOrder++, index += step) { + sum += coeff->data.F64[index] * polyValues->data.F64[yOrder][xOrder]; + } + } + return sum; +} + +inline double p_pmSubtractionSolutionCoeff(const pmSubtractionKernels *kernels, const psImage *polyValues, + int index, bool wantDual) +{ +#if 0 + // This is probably in a tight loop, so don't check inputs + PM_ASSERT_SUBTRACTION_KERNELS_NON_NULL(kernels, NAN); + PM_ASSERT_SUBTRACTION_KERNELS_SOLUTION(kernels, NAN); + PS_ASSERT_IMAGE_NON_NULL(polyValues, NAN); + PS_ASSERT_INT_POSITIVE(index, NAN); +#endif + + psVector *solution = wantDual ? kernels->solution2 : kernels->solution1; // Solution vector + return p_pmSubtractionCalculatePolynomial(solution, polyValues, kernels->spatialOrder, index, + kernels->num); +} + +inline double p_pmSubtractionSolutionNorm(const pmSubtractionKernels *kernels) +{ +#if 0 + // This is probably in a tight loop, so don't check inputs + PM_ASSERT_SUBTRACTION_KERNELS_NON_NULL(kernels, NAN); + PM_ASSERT_SUBTRACTION_KERNELS_SOLUTION(kernels, NAN); + PS_ASSERT_IMAGE_NON_NULL(polyValues, NAN); +#endif + + int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation + return kernels->solution1->data.F64[normIndex]; +} + +inline double p_pmSubtractionSolutionBackground(const pmSubtractionKernels *kernels, + const psImage *polyValues) +{ +#if 0 + // This is probably in a tight loop, so don't check inputs + PM_ASSERT_SUBTRACTION_KERNELS_NON_NULL(kernels, NAN); + PM_ASSERT_SUBTRACTION_KERNELS_SOLUTION(kernels, NAN); + PS_ASSERT_IMAGE_NON_NULL(polyValues, NAN); +#endif + + int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index for background + return p_pmSubtractionCalculatePolynomial(kernels->solution1, polyValues, kernels->bgOrder, bgIndex, 1); +} + +////////////////////////////////////////////////////////////////////////////////////////////////////////////// +// Public functions +////////////////////////////////////////////////////////////////////////////////////////////////////////////// + +bool pmSubtractionCalculateEquation(pmSubtractionStampList *stamps, const pmSubtractionKernels *kernels) +{ + PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(stamps, false); + PM_ASSERT_SUBTRACTION_KERNELS_NON_NULL(kernels, false); + + int footprint = stamps->footprint; // Half-size of stamps + int spatialOrder = kernels->spatialOrder; // Maximum order of spatial variation + int numKernels = kernels->num; // Number of kernel basis functions + int numSpatial = PM_SUBTRACTION_POLYTERMS(spatialOrder); // Number of spatial variations + int numBackground = PM_SUBTRACTION_POLYTERMS(kernels->bgOrder); // Number of background terms + + // Total number of parameters to solve for: coefficient of each kernel basis function, multipled by the + // number of coefficients for the spatial polynomial, normalisation and a constant background offset. + int numParams = numKernels * numSpatial + 1 + numBackground; + + psImage *polyValues = NULL; // Polynomial terms + + // We iterate over each stamp, allocate the matrix and vectors if + // necessary, and then calculate those matrix/vectors. + 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; + } + + // Generate convolutions + if (!pmSubtractionConvolveStamp(stamp, kernels, footprint)) { + psError(PS_ERR_UNKNOWN, false, "Unable to convolve stamp %d.", i); + psFree(polyValues); + return NULL; + } + + polyValues = p_pmSubtractionPolynomial(polyValues, spatialOrder, stamp->xNorm, stamp->yNorm); + + stamp->vector = psVectorAlloc(numParams, PS_TYPE_F64); +#ifdef TESTING + psVectorInit(stamp->vector, NAN); +#endif + + bool status; // Status of least-squares matrix/vector calculation + switch (kernels->mode) { + case PM_SUBTRACTION_MODE_TARGET: + case PM_SUBTRACTION_MODE_1: + stamp->matrix1 = psImageAlloc(numParams, numParams, PS_TYPE_F64); +#ifdef TESTING + psImageInit(stamp->matrix1, NAN); +#endif + status = calculateMatrix(stamp->matrix1, kernels, stamp->convolutions1, stamp->image1, + stamp->image2, stamp->weight, polyValues, footprint, true); + status &= calculateVector(stamp->vector, kernels, stamp->convolutions1, stamp->image1, + stamp->image2, stamp->weight, polyValues, footprint); + break; + case PM_SUBTRACTION_MODE_2: + stamp->matrix1 = psImageAlloc(numParams, numParams, PS_TYPE_F64); +#ifdef TESTING + psImageInit(stamp->matrix1, NAN); +#endif + status = calculateMatrix(stamp->matrix1, kernels, stamp->convolutions2, stamp->image2, + stamp->image1, stamp->weight, polyValues, footprint, true); + status &= calculateVector(stamp->vector, kernels, stamp->convolutions2, stamp->image2, + stamp->image1, stamp->weight, polyValues, footprint); + break; + case PM_SUBTRACTION_MODE_DUAL: + stamp->matrix1 = psImageAlloc(numParams, numParams, PS_TYPE_F64); + stamp->matrix2 = psImageAlloc(numKernels * numSpatial, numKernels * numSpatial, PS_TYPE_F64); + stamp->matrixX = psImageAlloc(numParams, numKernels * numSpatial, PS_TYPE_F64); +#ifdef TESTING + psImageInit(stamp->matrix1, NAN); + psImageInit(stamp->matrix2, NAN); + psImageInit(stamp->matrixX, NAN); +#endif + status = calculateMatrix(stamp->matrix1, kernels, stamp->convolutions1, stamp->image1, + stamp->image2, stamp->weight, polyValues, footprint, true); + status &= calculateMatrix(stamp->matrix2, kernels, stamp->convolutions2, stamp->image2, + stamp->image1, stamp->weight, polyValues, footprint, false); + // XXX Not sure I'm passing the correct stamp->image to calculateMatrixCross: + status &= calculateMatrixCross(stamp->matrixX, kernels, stamp->convolutions1, + stamp->convolutions2, stamp->image2, stamp->weight, polyValues, + footprint); + status &= calculateVector(stamp->vector, kernels, stamp->convolutions1, stamp->image1, + stamp->image2, stamp->weight, polyValues, footprint); + break; + default: + psAbort("Unsupported subtraction mode: %x", kernels->mode); + } + + if (!status) { + stamp->status = PM_SUBTRACTION_STAMP_REJECTED; + psWarning("Rejecting stamp %d (%d,%d) because of bad equation", + i, (int)(stamp->x + 0.5), (int)(stamp->y + 0.5)); + } else { + stamp->status = PM_SUBTRACTION_STAMP_USED; + } + +#ifdef TESTING + if (psTraceGetLevel("psModules.imcombine.equation") >= 10) { + psString matrixName = NULL; + psStringAppend(&matrixName, "matrix1_%d.fits", i); + psFits *matrixFile = psFitsOpen(matrixName, "w"); + psFree(matrixName); + psFitsWriteImage(matrixFile, NULL, stamp->matrix1, 0, NULL); + psFitsClose(matrixFile); + + matrixName = NULL; + psStringAppend(&matrixName, "vector_%d.fits", i); + psImage *dummy = psImageAlloc(stamp->vector->n, 1, PS_TYPE_F64); + memcpy(dummy->data.F64[0], stamp->vector->data.F64, + PSELEMTYPE_SIZEOF(PS_TYPE_F64) * stamp->vector->n); + matrixFile = psFitsOpen(matrixName, "w"); + psFree(matrixName); + psFitsWriteImage(matrixFile, NULL, dummy, 0, NULL); + psFree(dummy); + psFitsClose(matrixFile); + + + if (kernels->mode == PM_SUBTRACTION_MODE_DUAL) { + matrixName = NULL; + psStringAppend(&matrixName, "matrix2_%d.fits", i); + matrixFile = psFitsOpen(matrixName, "w"); + psFree(matrixName); + psFitsWriteImage(matrixFile, NULL, stamp->matrix2, 0, NULL); + psFitsClose(matrixFile); + + matrixName = NULL; + psStringAppend(&matrixName, "matrixX_%d.fits", i); + matrixFile = psFitsOpen(matrixName, "w"); + psFree(matrixName); + psFitsWriteImage(matrixFile, NULL, stamp->matrixX, 0, NULL); + psFitsClose(matrixFile); + } + } +#endif + + } + psFree(polyValues); + + return true; +} + +bool pmSubtractionSolveEquation(pmSubtractionKernels *kernels, const pmSubtractionStampList *stamps, + float tolerance) +{ + PM_ASSERT_SUBTRACTION_KERNELS_NON_NULL(kernels, false); + PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(stamps, false); + if (kernels->mode == PM_SUBTRACTION_MODE_DUAL) { + PS_ASSERT_FLOAT_LARGER_THAN(tolerance, 0.0, false); + } + + // Check inputs + long numParams = -1; // Number of parameters + long numParams2 = 0; // Number of parameters for part solution (DUAL mode) + for (int i = 0; i < stamps->num; i++) { + pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest + PS_ASSERT_PTR_NON_NULL(stamp, false); + if (stamp->status != PM_SUBTRACTION_STAMP_USED) { + continue; + } + + if (kernels->mode == PM_SUBTRACTION_MODE_DUAL) { + PS_ASSERT_IMAGE_NON_NULL(stamp->matrix1, false); + PS_ASSERT_IMAGE_NON_NULL(stamp->matrix2, false); + PS_ASSERT_IMAGE_NON_NULL(stamp->matrixX, false); + if (numParams == -1) { + numParams = stamp->matrix1->numCols; + numParams2 = stamp->matrix2->numCols; + } + PS_ASSERT_IMAGE_SIZE(stamp->matrix1, numParams, numParams, false); + PS_ASSERT_IMAGE_SIZE(stamp->matrix2, numParams2, numParams2, false); + // XXX Not sure I've got this the right way around + PS_ASSERT_IMAGE_SIZE(stamp->matrixX, numParams2, numParams, false); + PS_ASSERT_IMAGE_TYPE(stamp->matrix1, PS_TYPE_F64, false); + PS_ASSERT_IMAGE_TYPE(stamp->matrix2, PS_TYPE_F64, false); + PS_ASSERT_IMAGE_TYPE(stamp->matrixX, PS_TYPE_F64, false); + PS_ASSERT_VECTOR_NON_NULL(stamp->vector, false); + PS_ASSERT_VECTOR_SIZE(stamp->vector, numParams, false); + PS_ASSERT_VECTOR_TYPE(stamp->vector, PS_TYPE_F64, false); + } else { + PS_ASSERT_VECTOR_NON_NULL(stamp->vector, false); + if (numParams == -1) { + numParams = stamp->vector->n; + } + PS_ASSERT_VECTOR_SIZE(stamp->vector, numParams, false); + PS_ASSERT_VECTOR_TYPE(stamp->vector, PS_TYPE_F64, false); + PS_ASSERT_IMAGE_NON_NULL(stamp->matrix1, false); + PS_ASSERT_IMAGE_SIZE(stamp->matrix1, numParams, numParams, false); + PS_ASSERT_IMAGE_TYPE(stamp->matrix1, PS_TYPE_F64, false); + } + } + if (numParams == -1) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, "No suitable stamps found."); + return NULL; + } + + if (kernels->mode != PM_SUBTRACTION_MODE_DUAL) { + // Accumulate the least-squares matricies and vectors + psImage *sumMatrix = psImageAlloc(numParams, numParams, PS_TYPE_F64); // Combined matrix + psVector *sumVector = psVectorAlloc(numParams, PS_TYPE_F64); // Combined vector + psVectorInit(sumVector, 0.0); + psImageInit(sumMatrix, 0.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) { + (void)psBinaryOp(sumMatrix, sumMatrix, "+", stamp->matrix1); + (void)psBinaryOp(sumVector, sumVector, "+", stamp->vector); + } + } + + psVector *permutation = NULL; // Permutation vector, required for LU decomposition + psImage *luMatrix = psMatrixLUD(NULL, &permutation, sumMatrix); + psFree(sumMatrix); + if (!luMatrix) { + psError(PS_ERR_UNKNOWN, true, "LU Decomposition of least-squares matrix failed.\n"); + psFree(sumVector); + psFree(luMatrix); + psFree(permutation); + return NULL; + } + kernels->solution1 = psMatrixLUSolve(kernels->solution1, luMatrix, sumVector, permutation); + psFree(sumVector); + psFree(luMatrix); + psFree(permutation); + if (!kernels->solution1) { + psError(PS_ERR_UNKNOWN, true, "Failed to solve the least-squares system.\n"); + return NULL; + } + } else { + // Dual convolution solution + + // Accumulation of stamp matrices/vectors + psImage *sumMatrix1 = psImageAlloc(numParams, numParams, PS_TYPE_F64); + psImage *sumMatrix2 = psImageAlloc(numParams2, numParams2, PS_TYPE_F64); + psImage *sumMatrixX = psImageAlloc(numParams, numParams2, PS_TYPE_F64); + psVector *sumVector = psVectorAlloc(numParams, PS_TYPE_F64); + psImageInit(sumMatrix1, 0.0); + psImageInit(sumMatrix2, 0.0); + psImageInit(sumMatrixX, 0.0); + psVectorInit(sumVector, 0.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) { + (void)psBinaryOp(sumMatrix1, sumMatrix1, "+", stamp->matrix1); + (void)psBinaryOp(sumMatrix2, sumMatrix2, "+", stamp->matrix2); + (void)psBinaryOp(sumMatrixX, sumMatrixX, "+", stamp->matrixX); + (void)psBinaryOp(sumVector, sumVector, "+", stamp->vector); + } + } + + psVector *permutation1, *permutation2; // Permutation vectors for LUD + psImage *lu1 = psMatrixLUD(NULL, &permutation1, sumMatrix1); // LUD of matrix 1 + psImage *lu2 = psMatrixLUD(NULL, &permutation2, sumMatrix2); // LUD of matrix 2 + psFree(sumMatrix1); + psFree(sumMatrix2); + if (!lu1 || !lu2) { + psError(PS_ERR_UNKNOWN, false, "Unable to generate LU decomposed matrices."); + psFree(permutation1); + psFree(permutation2); + psFree(lu1); + psFree(lu2); + psFree(sumMatrixX); + psFree(sumVector); + return NULL; + } + + psVector *lastSol2 = psVectorAlloc(numParams2, PS_TYPE_F64); // Value of sol2 in last iteration + psVectorInit(lastSol2, 0.0); + + double diff = NAN; // Difference between iterations + psVector *vector = NULL; // RHS for equations + do { + if (!vector) { + // Start with traditional Alard-Lupton solution + vector = psVectorCopy(NULL, sumVector, PS_TYPE_F64); + } else { + // Generate vector for RHS of equation 1 + vector->n = numParams; + psVector *sol2 = kernels->solution2; + assert(sol2); + for (int i = 0; i < numParams; i++) { + double value = 0.0; // Value of vector element + for (int j = 0; j < numParams2; j++) { + // XXX Not sure I've got this the right way around + value += sol2->data.F64[j] * sumMatrixX->data.F64[j][i]; + } + vector->data.F64[i] = value; + } + } + + kernels->solution1 = psMatrixLUSolve(kernels->solution1, lu1, vector, permutation1); + if (!kernels->solution1) { + psError(PS_ERR_UNKNOWN, false, "Unable to solve initial matrix equation."); + psFree(permutation1); + psFree(permutation2); + psFree(lu1); + psFree(lu2); + psFree(sumMatrixX); + psFree(sumVector); + psFree(vector); + psFree(lastSol2); + return NULL; + } + + // Generate vector for RHS of equation 2 + vector->n = numParams2; + psVector *sol1 = kernels->solution1; // Solution of first equation + for (int i = 0; i < numParams2; i++) { + double value = 0.0; // Value of vector element + for (int j = 0; j < numParams; j++) { + // XXX Not sure I've got this the right way around + value += sol1->data.F64[j] * sumMatrixX->data.F64[i][j]; + } + vector->data.F64[i] = value; + } + + kernels->solution2 = psMatrixLUSolve(kernels->solution2, lu2, vector, permutation2); + if (!kernels->solution2) { + psError(PS_ERR_UNKNOWN, false, "Unable to solve matrix equation 2."); + psFree(permutation1); + psFree(permutation2); + psFree(lu1); + psFree(lu2); + psFree(sumMatrixX); + psFree(sumVector); + psFree(vector); + psFree(lastSol2); + return NULL; + } + + // Get difference produced by iteration + diff = 0.0; + psVector *sol2 = kernels->solution2; // Solution of second equation + for (int i = 0; i < numParams2; i++) { + diff += fabs(sol2->data.F64[i] - lastSol2->data.F64[i]) / sol2->data.F64[i]; + } + lastSol2 = psVectorCopy(lastSol2, sol2, PS_TYPE_F64); + } while (diff > tolerance); + + psFree(permutation1); + psFree(permutation2); + psFree(lu1); + psFree(lu2); + psFree(sumMatrixX); + psFree(sumVector); + psFree(vector); + psFree(lastSol2); + } + + 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]); + } + 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]); + } + } + } + + return true; +} + +psVector *pmSubtractionCalculateDeviations(pmSubtractionStampList *stamps, + const 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); + + psVector *deviations = psVectorAlloc(stamps->num, PS_TYPE_F32); // Mean deviation for stamps + int footprint = stamps->footprint; // Half-size of stamps + long numPixels = PS_SQR(2 * footprint + 1); // Number of pixels in footprint + double devNorm = 1.0 / (double)numPixels; // Normalisation for deviations + int numKernels = kernels->num; // Number of kernels + + psVector *coeff1 = psVectorAlloc(numKernels, PS_TYPE_F64); // Coefficients for first image + psVector *coeff2 = kernels->mode == PM_SUBTRACTION_MODE_DUAL ? + psVectorAlloc(numKernels, PS_TYPE_F64) : NULL; // Coefficients for second image + psImage *polyValues = NULL; // Polynomial values + psKernel *residual = psKernelAlloc(-footprint, footprint, -footprint, footprint); // Residual image + + 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) { + deviations->data.F32[i] = NAN; + continue; + } + + // Calculate coefficients of the kernel basis functions + polyValues = p_pmSubtractionPolynomial(polyValues, kernels->spatialOrder, stamp->xNorm, stamp->yNorm); + for (int j = 0; j < numKernels; j++) { + coeff1->data.F64[j] = p_pmSubtractionSolutionCoeff(kernels, polyValues, j, false); + } + double norm = p_pmSubtractionSolutionNorm(kernels); // Normalisation + double background = p_pmSubtractionSolutionBackground(kernels, polyValues);// Difference in background + + // Calculate residuals + psKernel *weight = stamp->weight; // Weight postage stamp + 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_TARGET: + 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 = coeff1->data.F64[j]; // 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++) { + residual->kernel[y][x] += target->kernel[y][x] - background - source->kernel[y][x] * norm; + } + } + } 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 + + for (int j = 0; j < numKernels; j++) { + coeff2->data.F64[j] = p_pmSubtractionSolutionCoeff(kernels, polyValues, j, true); + } + +#ifdef TESTING + psImageInit(residual->image, 0.0); +#endif + + 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 coefficient1 = coeff1->data.F64[j]; // Coefficient for convolution 1 + double coefficient2 = coeff2->data.F64[j]; // Coefficient for convolution 2 + + for (int y = - footprint; y <= footprint; y++) { + for (int x = - footprint; x <= footprint; x++) { + residual->kernel[y][x] += conv1->kernel[y][x] * coefficient1 - + conv2->kernel[y][x] * coefficient2; + } + } + } + for (int y = - footprint; y <= footprint; y++) { + for (int x = - footprint; x <= footprint; x++) { + residual->kernel[y][x] += image1->kernel[y][x] * norm - background; + } + } + } + + double deviation = 0.0; // Sum of differences + for (int y = - footprint; y <= footprint; y++) { + for (int x = - footprint; x <= footprint; x++) { + deviation += PS_SQR(residual->kernel[y][x]) / weight->kernel[y][x]; + } + } + deviations->data.F32[i] = sqrtf(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]); + 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)); + 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); + } +#endif + + } + psFree(residual); + psFree(polyValues); + psFree(coeff1); + + return deviations; +} Index: trunk/psModules/src/imcombine/pmSubtractionEquation.h =================================================================== --- trunk/psModules/src/imcombine/pmSubtractionEquation.h (revision 15756) +++ trunk/psModules/src/imcombine/pmSubtractionEquation.h (revision 15756) @@ -0,0 +1,47 @@ +#ifndef PM_SUBTRACTION_EQUATION_H +#define PM_SUBTRACTION_EQUATION_H + +#include "pmSubtractionStamps.h" +#include "pmSubtractionKernels.h" + +/// Calculate the least-squares equation to match the image quality +bool pmSubtractionCalculateEquation(pmSubtractionStampList *stamps, ///< Stamps + const pmSubtractionKernels *kernels ///< Kernel parameters + ); + +/// Solve the least-squares equation to match the image quality +bool pmSubtractionSolveEquation(pmSubtractionKernels *kernels, ///< Kernel parameters + const pmSubtractionStampList *stamps, ///< Stamps + float tolerance ///< Tolerance for iteration (dual mode only) + ); + +/// Calculate deviations +psVector *pmSubtractionCalculateDeviations(pmSubtractionStampList *stamps, ///< Stamps + const pmSubtractionKernels *kernels ///< Kernel parameters + ); + +/// Calculate the value of a polynomial, specified by coefficients and polynomial values +inline double p_pmSubtractionCalculatePolynomial(const psVector *coeff, ///< Coefficients + const psImage *polyValues, ///< Polynomial values + int order, ///< Order of polynomials + int index, ///< Index at which to begin + int step ///< Step between subsequent indices + ); + +/// Return the specified coefficient in the solution +inline double p_pmSubtractionSolutionCoeff(const pmSubtractionKernels *kernels, ///< Kernel parameters + const psImage *polyValues, ///< Polynomial values + int index, ///< Coefficient index to calculate + bool wantDual ///< Calculate the coefficient for the dual solution? + ); + +/// Return the normalisation in the solution +inline double p_pmSubtractionSolutionNorm(const pmSubtractionKernels *kernels ///< Kernel parameters + ); + +/// Return the background (difference) in the solution +inline double p_pmSubtractionSolutionBackground(const pmSubtractionKernels *kernels, ///< Kernel parameters + const psImage *polyValues ///< Polynomial values + ); + +#endif Index: trunk/psModules/src/imcombine/pmSubtractionKernels.c =================================================================== --- trunk/psModules/src/imcombine/pmSubtractionKernels.c (revision 15745) +++ trunk/psModules/src/imcombine/pmSubtractionKernels.c (revision 15756) @@ -4,7 +4,9 @@ #include +#include #include #include +#include "pmSubtraction.h" #include "pmSubtractionKernels.h" @@ -22,7 +24,9 @@ psFree(kernels->vStop); psFree(kernels->preCalc); -} - -// Raise a number to an integer power + psFree(kernels->solution1); + psFree(kernels->solution2); +} + +// Raise an integer to an integer power static inline long power(int value, // Value int exp // Exponent @@ -45,16 +49,9 @@ bool p_pmSubtractionKernelsAddGrid(pmSubtractionKernels *kernels, int start, int size) { - PS_ASSERT_PTR_NON_NULL(kernels, false); - PS_ASSERT_VECTOR_NON_NULL(kernels->u, false); - PS_ASSERT_VECTOR_NON_NULL(kernels->v, false); - PS_ASSERT_VECTOR_NON_NULL(kernels->widths, false); - PS_ASSERT_VECTOR_TYPE(kernels->u, PS_TYPE_S32, false); - PS_ASSERT_VECTOR_TYPE(kernels->v, PS_TYPE_S32, false); - PS_ASSERT_VECTOR_TYPE(kernels->widths, PS_TYPE_F32, false); - PS_ASSERT_ARRAY_NON_NULL(kernels->preCalc, false); + PM_ASSERT_SUBTRACTION_KERNELS_NON_NULL(kernels, false); PS_ASSERT_INT_NONNEGATIVE(start, false); PS_ASSERT_INT_NONNEGATIVE(size, false); - int numNew = PS_SQR(2 * size + 1); // Number of new kernel parameters to add + int numNew = PS_SQR(2 * size + 1) - 1; // Number of new kernel parameters to add // Ensure the sizes match @@ -68,4 +65,9 @@ for (int v = - size, index = start; v <= size; v++) { for (int u = - size; u <= size; u++, index++) { + if (v == 0 && u == 0) { + // Skip normalisation component: added explicitly + index--; + continue; + } kernels->widths->data.F32[index] = NAN; kernels->u->data.S32[index] = u; @@ -81,5 +83,6 @@ pmSubtractionKernels *p_pmSubtractionKernelsRawISIS(int size, int spatialOrder, - const psVector *fwhms, const psVector *orders) + const psVector *fwhms, const psVector *orders, + pmSubtractionMode mode) { PS_ASSERT_VECTOR_NON_NULL(fwhms, NULL); @@ -97,10 +100,10 @@ for (int i = 0; i < numGaussians; i++) { int gaussOrder = orders->data.S32[i]; // Polynomial order to apply to Gaussian - psStringAppend(¶ms, "(%.2f,%d)", fwhms->data.F32[i], orders->data.S32[i]); + 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); // The kernels + size, spatialOrder, mode); // The kernels psStringAppend(&kernels->description, "ISIS(%d,%s,%d)", size, params, spatialOrder); @@ -149,6 +152,4 @@ } - kernels->subIndex = -1; - return kernels; } @@ -159,5 +160,5 @@ pmSubtractionKernels *pmSubtractionKernelsAlloc(int numBasisFunctions, pmSubtractionKernelsType type, - int size, int spatialOrder) + int size, int spatialOrder, pmSubtractionMode mode) { pmSubtractionKernels *kernels = psAlloc(sizeof(pmSubtractionKernels)); // Kernels, to return @@ -173,22 +174,24 @@ kernels->uStop = NULL; kernels->vStop = NULL; - kernels->subIndex = 0; kernels->size = size; kernels->inner = 0; kernels->spatialOrder = spatialOrder; kernels->bgOrder = 0; - - return kernels; -} - -pmSubtractionKernels *pmSubtractionKernelsPOIS(int size, int spatialOrder) + kernels->mode = mode; + kernels->solution1 = NULL; + kernels->solution2 = NULL; + + return kernels; +} + +pmSubtractionKernels *pmSubtractionKernelsPOIS(int size, int spatialOrder, pmSubtractionMode mode) { PS_ASSERT_INT_POSITIVE(size, NULL); PS_ASSERT_INT_NONNEGATIVE(spatialOrder, NULL); - int num = PS_SQR(2 * size + 1); // Number of basis functions + int num = PS_SQR(2 * size + 1) - 1; // Number of basis functions pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_POIS, - size, spatialOrder); // The kernels + size, spatialOrder, mode); // The kernels psStringAppend(&kernels->description, "POIS(%d,%d)", size, spatialOrder); psLogMsg("psModules.imcombine", PS_LOG_INFO, "POIS kernel: %d,%d --> %d elements", @@ -199,8 +202,4 @@ } - kernels->subIndex = num / 2; - assert(kernels->u->data.S32[kernels->subIndex] == 0 && - kernels->v->data.S32[kernels->subIndex] == 0); - return kernels; } @@ -208,23 +207,22 @@ pmSubtractionKernels *pmSubtractionKernelsISIS(int size, int spatialOrder, - const psVector *fwhms, const psVector *orders) + const psVector *fwhms, const psVector *orders, + pmSubtractionMode mode) { pmSubtractionKernels *kernels = p_pmSubtractionKernelsRawISIS(size, spatialOrder, - fwhms, orders); // Kernels + fwhms, orders, mode); // Kernels if (!kernels) { return NULL; } - // Subtract a reference kernel from all the others to maintain flux scaling + // Subtract normalisation from all the others to maintain flux scaling if (spatialOrder > 0) { - int subIndex = 0; // Index of kernel to subtract from others - psKernel *subtract = kernels->preCalc->data[subIndex]; // Kernel to subtract from others int numGaussians = fwhms->n; // Number of Gaussians for (int i = 0, index = 0; i < numGaussians; i++) { for (int uOrder = 0; uOrder <= orders->data.S32[i]; uOrder++) { for (int vOrder = 0; vOrder <= orders->data.S32[i] - uOrder; vOrder++, index++) { - if (uOrder % 2 == 0 && vOrder % 2 == 0 && index != subIndex) { + if (uOrder % 2 == 0 && vOrder % 2 == 0) { psKernel *kernel = kernels->preCalc->data[index]; // Kernel to subtract - psBinaryOp(kernel->image, kernel->image, "-", subtract->image); + kernel->kernel[0][0] -= 1.0; } } @@ -242,16 +240,19 @@ psFitsWriteImage(kernelFile, NULL, kernel->image, 0, NULL); psFitsClose(kernelFile); - } - } - - return kernels; -} - -/// Generate SPAM kernels -pmSubtractionKernels *pmSubtractionKernelsSPAM(int size, ///< Half-size of the kernel - int spatialOrder, ///< Order of spatial variations - int inner, ///< Inner radius to preserve unbinned - int binning ///< Kernel binning factor - ) + double sum = 0.0; + for (int y = 0; y < kernel->image->numRows; y++) { + for (int x = 0; x < kernel->image->numCols; x++) { + sum += kernel->image->data.F32[y][x]; + } + } + psTrace("psModules.imcombine.kernel", 10, "Kernel %d sum: %lf\n", i, sum); + } + } + + return kernels; +} + +pmSubtractionKernels *pmSubtractionKernelsSPAM(int size, int spatialOrder, int inner, int binning, + pmSubtractionMode mode) { PS_ASSERT_INT_POSITIVE(size, NULL); @@ -269,10 +270,10 @@ psTrace("psModules.imcombine", 3, "Inner: %d Outer: %d\n", numInner, numOuter); - int num = PS_SQR(2 * numTotal + 1); // Number of basis functions + int num = PS_SQR(2 * numTotal + 1) - 1; // Number of basis functions psTrace("psModules.imcombine", 3, "Number of basis functions: %d\n", num); pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_SPAM, - size, spatialOrder); // The kernels + size, spatialOrder, mode); // The kernels psStringAppend(&kernels->description, "SPAM(%d,%d,%d,%d)", size, inner, binning, spatialOrder); @@ -280,6 +281,6 @@ size, inner, binning, spatialOrder, num); - kernels->uStop = psVectorAlloc(num, PS_TYPE_F32); - kernels->vStop = psVectorAlloc(num, PS_TYPE_F32); + kernels->uStop = psVectorAlloc(num, PS_TYPE_S32); + kernels->vStop = psVectorAlloc(num, PS_TYPE_S32); psVector *locations = psVectorAlloc(2 * numTotal + 1, PS_TYPE_S32); // Locations for each kernel element @@ -314,4 +315,9 @@ for (int j = - numTotal; j <= numTotal; j++, index++) { + if (i == 0 && j == 0) { + // Skip normalisation component: added explicitly + index--; + continue; + } int v = locations->data.S32[numTotal + j]; // Location of pixel int vStop = v + widths->data.S32[numTotal + j]; // Width of pixel @@ -327,10 +333,4 @@ } - kernels->subIndex = num / 2; - assert(kernels->u->data.S32[kernels->subIndex] == 0 && - kernels->v->data.S32[kernels->subIndex] == 0 && - kernels->uStop->data.S32[kernels->subIndex] == 0 && - kernels->vStop->data.S32[kernels->subIndex] == 0); - psFree(locations); psFree(widths); @@ -340,9 +340,5 @@ -/// Generate FRIES kernels -pmSubtractionKernels *pmSubtractionKernelsFRIES(int size, ///< Half-size of the kernel - int spatialOrder, ///< Order of spatial variations - int inner ///< Inner radius to preserve unbinned - ) +pmSubtractionKernels *pmSubtractionKernelsFRIES(int size, int spatialOrder, int inner, pmSubtractionMode mode) { PS_ASSERT_INT_POSITIVE(size, NULL); @@ -366,10 +362,10 @@ psTrace("psModules.imcombine", 3, "Inner: %d Outer: %d\n", numInner, numOuter); - int num = PS_SQR(2 * numTotal + 1); // Number of basis functions + int num = PS_SQR(2 * numTotal + 1) - 1; // Number of basis functions psTrace("psModules.imcombine", 3, "Number of basis functions: %d\n", num); pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_FRIES, - size, spatialOrder); // The kernels + size, spatialOrder, mode); // The kernels psStringAppend(&kernels->description, "FRIES(%d,%d,%d)", size, inner, spatialOrder); @@ -377,6 +373,6 @@ size, inner, spatialOrder, num); - kernels->uStop = psVectorAlloc(num, PS_TYPE_F32); - kernels->vStop = psVectorAlloc(num, PS_TYPE_F32); + kernels->uStop = psVectorAlloc(num, PS_TYPE_S32); + kernels->vStop = psVectorAlloc(num, PS_TYPE_S32); psVector *start = psVectorAlloc(2 * numTotal + 1, PS_TYPE_S32); @@ -409,4 +405,9 @@ int uStop = stop->data.S32[numTotal + i]; // Width of pixel for (int j = - numTotal; j <= numTotal; j++, index++) { + if (i == 0 && j == 0) { + // Skip normalisation component: added explicitly + index--; + continue; + } int v = start->data.S32[numTotal + j]; // Location of pixel int vStop = stop->data.S32[numTotal + j]; // Width of pixel @@ -422,10 +423,4 @@ } - kernels->subIndex = num / 2; - assert(kernels->u->data.S32[kernels->subIndex] == 0 && - kernels->v->data.S32[kernels->subIndex] == 0 && - kernels->uStop->data.S32[kernels->subIndex] == 0 && - kernels->vStop->data.S32[kernels->subIndex] == 0); - psFree(start); psFree(stop); @@ -436,5 +431,5 @@ // Grid United with Normal Kernel pmSubtractionKernels *pmSubtractionKernelsGUNK(int size, int spatialOrder, const psVector *fwhms, - const psVector *orders, int inner) + const psVector *orders, int inner, pmSubtractionMode mode) { PS_ASSERT_INT_POSITIVE(size, NULL); @@ -449,5 +444,5 @@ pmSubtractionKernels *kernels = p_pmSubtractionKernelsRawISIS(size, spatialOrder, - fwhms, orders); // Kernels + fwhms, orders, mode); // Kernels psStringPrepend(&kernels->description, "GUNK="); psStringAppend(&kernels->description, "+POIS(%d,%d)", inner, spatialOrder); @@ -469,5 +464,4 @@ int numISIS = kernels->num; // Number of ISIS kernels - int numInner = PS_SQR(2 * inner + 1); // Number of inner kernel elements if (!p_pmSubtractionKernelsAddGrid(kernels, numISIS, inner)) { @@ -475,13 +469,10 @@ } - kernels->subIndex = numInner/2 + numISIS; - assert(kernels->u->data.S32[kernels->subIndex] == 0 && - kernels->v->data.S32[kernels->subIndex] == 0); - return kernels; } // RINGS --- just what it says -pmSubtractionKernels *pmSubtractionKernelsRINGS(int size, int spatialOrder, int inner, int ringsOrder) +pmSubtractionKernels *pmSubtractionKernelsRINGS(int size, int spatialOrder, int inner, int ringsOrder, + pmSubtractionMode mode) { PS_ASSERT_INT_POSITIVE(size, NULL); @@ -505,14 +496,14 @@ } - int numInner = inner; // Number of pixels in the inner region + int numInner = inner - 1; // Number of pixels in the inner region int numOuter = fibNum; // Number of summed pixels in the outer region int numRings = numOuter + numInner; // Number of rings (not including the central pixel) - int numPoly = (ringsOrder + 1) * (ringsOrder + 2) / 2; // Number of polynomial variants of each ring - - int num = numRings * numPoly + 1; // Total number of basis functions + int numPoly = PM_SUBTRACTION_POLYTERMS(ringsOrder); // Number of polynomial variants of each ring + + int num = numRings * numPoly; // Total number of basis functions pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_RINGS, - size, spatialOrder); // The kernels + size, spatialOrder, mode); // The kernels psStringAppend(&kernels->description, "RINGS(%d,%d,%d,%d)", size, inner, ringsOrder, spatialOrder); @@ -522,13 +513,9 @@ // Set the Gaussian kernel parameters int fibIndex = 1, fibIndexMinus1 = 0; // Fibonnacci parameters - int radiusLast = 0; // Last radius - for (int i = 0, index = 0; i < numRings + 1; i++) { - + int radiusLast = 1; // Last radius + for (int i = 1, index = 0; i < numRings + 1; i++) { float lower2; // Lower limit of radius^2 float upper2; // Upper limit of radius^2 - if (i == 0) { - lower2 = 0; - upper2 = PS_SQR(0.5); - } else if (i <= numInner) { + if (i <= inner) { // A ring every pixel width float radius = i; @@ -572,5 +559,4 @@ for (int v = -size; v <= size; v++) { int v2 = PS_SQR(v); // Square of v -// float vPoly = power(v, vOrder); // Value of v^vOrder float vPoly = powf(v/(float)size, vOrder); // Value of v^vOrder @@ -579,5 +565,4 @@ int distance2 = u2 + v2; // Distance from the centre if (distance2 > lower2 && distance2 < upper2) { -// float uPoly = power(u, uOrder); // Value of u^uOrder float uPoly = powf(u/(float)size, uOrder); // Value of u^uOrder @@ -627,6 +612,4 @@ } - kernels->subIndex = 0; - return kernels; } @@ -634,19 +617,19 @@ pmSubtractionKernels *pmSubtractionKernelsGenerate(pmSubtractionKernelsType type, int size, int spatialOrder, const psVector *fwhms, const psVector *orders, int inner, - int binning, int ringsOrder) + int binning, int ringsOrder, pmSubtractionMode mode) { switch (type) { case PM_SUBTRACTION_KERNEL_POIS: - return pmSubtractionKernelsPOIS(size, spatialOrder); + return pmSubtractionKernelsPOIS(size, spatialOrder, mode); case PM_SUBTRACTION_KERNEL_ISIS: - return pmSubtractionKernelsISIS(size, spatialOrder, fwhms, orders); + return pmSubtractionKernelsISIS(size, spatialOrder, fwhms, orders, mode); case PM_SUBTRACTION_KERNEL_SPAM: - return pmSubtractionKernelsSPAM(size, spatialOrder, inner, binning); + return pmSubtractionKernelsSPAM(size, spatialOrder, inner, binning, mode); case PM_SUBTRACTION_KERNEL_FRIES: - return pmSubtractionKernelsFRIES(size, spatialOrder, inner); + return pmSubtractionKernelsFRIES(size, spatialOrder, inner, mode); case PM_SUBTRACTION_KERNEL_GUNK: - return pmSubtractionKernelsGUNK(size, spatialOrder, fwhms, orders, inner); + return pmSubtractionKernelsGUNK(size, spatialOrder, fwhms, orders, inner, mode); case PM_SUBTRACTION_KERNEL_RINGS: - return pmSubtractionKernelsRINGS(size, spatialOrder, inner, ringsOrder); + return pmSubtractionKernelsRINGS(size, spatialOrder, inner, ringsOrder, mode); default: psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Unknown kernel type: %x", type); Index: trunk/psModules/src/imcombine/pmSubtractionKernels.h =================================================================== --- trunk/psModules/src/imcombine/pmSubtractionKernels.h (revision 15745) +++ trunk/psModules/src/imcombine/pmSubtractionKernels.h (revision 15756) @@ -31,7 +31,6 @@ 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) or ring radius (RINGS) + psVector *widths; ///< Gaussian FWHMs (ISIS) psVector *uStop, *vStop; ///< Width of kernel element (SPAM,FRIES only) - int subIndex; ///< Index of kernel to be subtracted (to maintain flux conservation) psArray *preCalc; ///< Array of images containing pre-calculated kernel (for ISIS) int size; ///< The half-size of the kernel @@ -39,5 +38,58 @@ 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 } pmSubtractionKernels; + +// Assertion to check pmSubtractionKernels +#define PM_ASSERT_SUBTRACTION_KERNELS_NON_NULL(KERNELS, RETURNVALUE) { \ + PS_ASSERT_PTR_NON_NULL(KERNELS, RETURNVALUE); \ + PS_ASSERT_STRING_NON_EMPTY((KERNELS)->description, RETURNVALUE); \ + PS_ASSERT_INT_POSITIVE((KERNELS)->num, RETURNVALUE); \ + PS_ASSERT_VECTOR_NON_NULL((KERNELS)->u, RETURNVALUE); \ + PS_ASSERT_VECTOR_NON_NULL((KERNELS)->v, RETURNVALUE); \ + PS_ASSERT_VECTOR_TYPE((KERNELS)->u, PS_TYPE_S32, RETURNVALUE); \ + PS_ASSERT_VECTOR_TYPE((KERNELS)->v, PS_TYPE_S32, RETURNVALUE); \ + PS_ASSERT_VECTOR_SIZE((KERNELS)->u, (KERNELS)->num, RETURNVALUE); \ + PS_ASSERT_VECTOR_SIZE((KERNELS)->v, (KERNELS)->num, RETURNVALUE); \ + if ((KERNELS)->type == PM_SUBTRACTION_KERNEL_ISIS) { \ + PS_ASSERT_VECTOR_NON_NULL((KERNELS)->widths, RETURNVALUE); \ + PS_ASSERT_VECTOR_TYPE((KERNELS)->widths, PS_TYPE_F32, RETURNVALUE); \ + PS_ASSERT_VECTOR_SIZE((KERNELS)->widths, (KERNELS)->num, RETURNVALUE); \ + } \ + if ((KERNELS)->uStop || (KERNELS)->vStop) { \ + PS_ASSERT_VECTOR_NON_NULL((KERNELS)->uStop, RETURNVALUE); \ + PS_ASSERT_VECTOR_NON_NULL((KERNELS)->vStop, RETURNVALUE); \ + PS_ASSERT_VECTOR_TYPE((KERNELS)->uStop, PS_TYPE_S32, RETURNVALUE); \ + PS_ASSERT_VECTOR_TYPE((KERNELS)->vStop, PS_TYPE_S32, RETURNVALUE); \ + PS_ASSERT_VECTOR_SIZE((KERNELS)->uStop, (KERNELS)->num, RETURNVALUE); \ + PS_ASSERT_VECTOR_SIZE((KERNELS)->vStop, (KERNELS)->num, RETURNVALUE); \ + } \ + if ((KERNELS)->preCalc) { \ + PS_ASSERT_ARRAY_NON_NULL((KERNELS)->preCalc, RETURNVALUE); \ + PS_ASSERT_ARRAY_SIZE((KERNELS)->preCalc, (KERNELS)->num, RETURNVALUE); \ + } \ + PS_ASSERT_INT_NONNEGATIVE((KERNELS)->size, RETURNVALUE); \ + PS_ASSERT_INT_NONNEGATIVE((KERNELS)->inner, RETURNVALUE); \ + PS_ASSERT_INT_NONNEGATIVE((KERNELS)->spatialOrder, RETURNVALUE); \ + PS_ASSERT_INT_NONNEGATIVE((KERNELS)->bgOrder, RETURNVALUE); \ +} + +// Assertion to check that the solution is attached +#define PM_ASSERT_SUBTRACTION_KERNELS_SOLUTION(KERNELS, RETURNVALUE) { \ + PS_ASSERT_VECTOR_NON_NULL((KERNELS)->solution1, RETURNVALUE); \ + PS_ASSERT_VECTOR_TYPE((KERNELS)->solution1, PS_TYPE_F64, RETURNVALUE); \ + PS_ASSERT_VECTOR_SIZE((KERNELS)->solution1, \ + (KERNELS)->num * PM_SUBTRACTION_POLYTERMS((KERNELS)->spatialOrder) + 1 + \ + PM_SUBTRACTION_POLYTERMS((KERNELS)->bgOrder), \ + RETURNVALUE); \ + if (kernels->mode == PM_SUBTRACTION_MODE_DUAL) { \ + PS_ASSERT_VECTOR_NON_NULL(kernels->solution2, RETURNVALUE); \ + PS_ASSERT_VECTOR_TYPE((KERNELS)->solution2, PS_TYPE_F64, RETURNVALUE); \ + PS_ASSERT_VECTOR_SIZE((KERNELS)->solution2, \ + (KERNELS)->num * PM_SUBTRACTION_POLYTERMS((KERNELS)->spatialOrder), \ + RETURNVALUE); \ + } \ +} /// Generate a delta-function grid for subtraction kernels (like the POIS kernel) @@ -54,10 +106,12 @@ pmSubtractionKernelsType type, ///< Kernel type int size, ///< Half-size of kernel - int spatialOrder ///< Order of spatial variations + int spatialOrder, ///< Order of spatial variations + pmSubtractionMode mode ///< Mode for subtraction ); /// Generate POIS kernels pmSubtractionKernels *pmSubtractionKernelsPOIS(int size, ///< Half-size of the kernel (in both dims) - int spatialOrder ///< Order of spatial variations + int spatialOrder, ///< Order of spatial variations + pmSubtractionMode mode ///< Mode for subtraction ); @@ -66,5 +120,6 @@ int spatialOrder, ///< Order of spatial variations const psVector *fwhms, ///< Gaussian FWHMs - const psVector *orders ///< Polynomial order of gaussians + const psVector *orders, ///< Polynomial order of gaussians + pmSubtractionMode mode ///< Mode for subtraction ); @@ -73,5 +128,6 @@ int spatialOrder, ///< Order of spatial variations const psVector *fwhms, ///< Gaussian FWHMs - const psVector *orders ///< Polynomial order of gaussians + const psVector *orders, ///< Polynomial order of gaussians + pmSubtractionMode mode ///< Mode for subtraction ); @@ -80,5 +136,6 @@ int spatialOrder, ///< Order of spatial variations int inner, ///< Inner radius to preserve unbinned - int binning ///< Kernel binning factor + int binning, ///< Kernel binning factor + pmSubtractionMode mode ///< Mode for subtraction ); @@ -86,5 +143,6 @@ pmSubtractionKernels *pmSubtractionKernelsFRIES(int size, ///< Half-size of the kernel int spatialOrder, ///< Order of spatial variations - int inner ///< Inner radius to preserve unbinned + int inner, ///< Inner radius to preserve unbinned + pmSubtractionMode mode ///< Mode for subtraction ); @@ -94,5 +152,6 @@ const psVector *fwhms, ///< Gaussian FWHMs const psVector *orders, ///< Polynomial order of gaussians - int inner ///< Inner radius containing grid of delta functions + int inner, ///< Inner radius containing grid of delta functions + pmSubtractionMode mode ///< Mode for subtraction ); @@ -101,5 +160,6 @@ int spatialOrder, ///< Order of spatial variations int inner, ///< Inner radius to preserve unbinned - int ringsOrder ///< Polynomial order + int ringsOrder, ///< Polynomial order + pmSubtractionMode mode ///< Mode for subtraction ); @@ -113,5 +173,6 @@ int inner, ///< Inner radius to preserve unbinned int binning, ///< Kernel binning factor - int ringsOrder ///< Polynomial order for RINGS + int ringsOrder, ///< Polynomial order for RINGS + pmSubtractionMode mode ///< Mode for subtraction ); Index: trunk/psModules/src/imcombine/pmSubtractionMask.c =================================================================== --- trunk/psModules/src/imcombine/pmSubtractionMask.c (revision 15756) +++ trunk/psModules/src/imcombine/pmSubtractionMask.c (revision 15756) @@ -0,0 +1,240 @@ +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include + +#include "pmSubtraction.h" + +#include "pmSubtractionMask.h" + +////////////////////////////////////////////////////////////////////////////////////////////////////////////// +// Private (file-static) functions +////////////////////////////////////////////////////////////////////////////////////////////////////////////// + +// Mark a pixel as blank in the image, mask and weight +static inline void markBlank(psImage *image, // Image to mark as blank + psImage *mask, // Mask to mark as blank (or NULL) + psImage *weight, // Weight map to mark as blank (or NULL) + int x, int y, // Coordinates to mark blank + psMaskType blank // Blank mask value + ) +{ + image->data.F32[y][x] = NAN; + if (mask) { + mask->data.PS_TYPE_MASK_DATA[y][x] |= blank; + } + if (weight) { + weight->data.F32[y][x] = NAN; + } + return; +} + +////////////////////////////////////////////////////////////////////////////////////////////////////////////// +// Public functions +////////////////////////////////////////////////////////////////////////////////////////////////////////////// + +psImage *pmSubtractionMask(const psImage *mask1, const psImage *mask2, psMaskType maskVal, + int size, int footprint, float badFrac, bool useFFT) +{ + PS_ASSERT_IMAGE_NON_NULL(mask1, NULL); + PS_ASSERT_IMAGE_TYPE(mask1, PS_TYPE_MASK, NULL); + if (mask2) { + PS_ASSERT_IMAGE_NON_NULL(mask2, NULL); + PS_ASSERT_IMAGE_TYPE(mask2, PS_TYPE_MASK, NULL); + PS_ASSERT_IMAGES_SIZE_EQUAL(mask2, mask1, NULL); + } + PS_ASSERT_INT_NONNEGATIVE(size, NULL); + PS_ASSERT_INT_NONNEGATIVE(footprint, NULL); + if (isfinite(badFrac)) { + PS_ASSERT_FLOAT_LARGER_THAN(badFrac, 0.0, NULL); + PS_ASSERT_FLOAT_LESS_THAN_OR_EQUAL(badFrac, 1.0, NULL); + } + + int numCols = mask1->numCols, numRows = mask1->numRows; // Size of the images + + // Dereference inputs for convenience + psMaskType **data1 = mask1->data.PS_TYPE_MASK_DATA; + psMaskType **data2 = NULL; + if (mask2) { + data2 = mask2->data.PS_TYPE_MASK_DATA; + } + + // First, a pass through to determine the fraction of bad pixels + if (isfinite(badFrac) && badFrac != 1.0) { + int numBad = 0; // Number of bad pixels + for (int y = 0; y < numRows; y++) { + for (int x = 0; x < numCols; x++) { + if (data1[y][x] & maskVal) { + numBad++; + continue; + } + if (data2 && data2[y][x] & maskVal) { + numBad++; + } + } + } + if (numBad > badFrac * numCols * numRows) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, + "Fraction of bad pixels (%d/%d=%f) exceeds limit (%f)\n", + numBad, numCols * numRows, (float)numBad/(float)(numCols * numRows), badFrac); + return NULL; + } + } + + // Worried about the masks for bad pixels and bad stamps colliding, so make our own mask + psImage *mask = psImageAlloc(numCols, numRows, PS_TYPE_MASK); // The global mask + psImageInit(mask, 0); + psMaskType **maskData = mask->data.PS_TYPE_MASK_DATA; // Dereference for convenience + + // Block out a border around the edge of the image + + // Bottom stripe + for (int y = 0; y < PS_MIN(size + footprint, numRows); y++) { + for (int x = 0; x < numCols; x++) { + maskData[y][x] |= PM_SUBTRACTION_MASK_BORDER; + } + } + // Either side + for (int y = PS_MIN(size + footprint, numRows); y < numRows - size - footprint; y++) { + for (int x = 0; x < PS_MIN(size + footprint, numCols); x++) { + maskData[y][x] |= PM_SUBTRACTION_MASK_BORDER; + } + for (int x = PS_MAX(numCols - size - footprint, 0); x < numCols; x++) { + maskData[y][x] |= PM_SUBTRACTION_MASK_BORDER; + } + } + // Top stripe + for (int y = PS_MAX(numRows - size - footprint, 0); y < numRows; y++) { + for (int x = 0; x < numCols; x++) { + maskData[y][x] |= PM_SUBTRACTION_MASK_BORDER; + } + } + + // XXX Could do something smarter here --- we will get images that are predominantly masked (where the + // skycell isn't overlapped by a large fraction by the observation), so that convolving around every bad + // pixel is wasting time. As a first cut, I've put in a check on the fraction of bad pixels, but we could + // imagine looking for the edge of big regions and convolving just at the edge. As a second cut, allow + // use of FFT convolution. + + for (int y = 0; y < numRows; y++) { + for (int x = 0; x < numCols; x++) { + if (data1[y][x] & maskVal) { + maskData[y][x] |= PM_SUBTRACTION_MASK_BAD_1; + } + if (data2 && data2[y][x] & maskVal) { + maskData[y][x] |= PM_SUBTRACTION_MASK_BAD_2; + } + } + } + + // Block out the entire stamp footprint around bad input pixels. + + // We want to block out with the CONVOLVE mask anything that would be bad if we convolved with a bad + // reference pixel (within 'size'). Then we want to block out with the FOOTPRINT mask everything within a + // footprint's distance of those (within 'footprint'). + + if (useFFT) { + if (!psImageConvolveMaskFFT(mask, mask, PM_SUBTRACTION_MASK_BAD_1 | PM_SUBTRACTION_MASK_BAD_2, + PM_SUBTRACTION_MASK_CONVOLVE_1 | PM_SUBTRACTION_MASK_FOOTPRINT_1, + -size, size, -size, size, 0.5)) { + psError(PS_ERR_UNKNOWN, false, "Unable to convolve bad pixels from mask 1."); + psFree(mask); + return NULL; + } + if (!psImageConvolveMaskFFT(mask, mask, PM_SUBTRACTION_MASK_BAD_2, + PM_SUBTRACTION_MASK_CONVOLVE_2 | PM_SUBTRACTION_MASK_FOOTPRINT_2, + -size, size, -size, size, 0.5)) { + psError(PS_ERR_UNKNOWN, false, "Unable to convolve bad pixels from mask 2."); + psFree(mask); + return NULL; + } + if (!psImageConvolveMaskFFT(mask, mask, PM_SUBTRACTION_MASK_CONVOLVE_1, + PM_SUBTRACTION_MASK_FOOTPRINT_1, + -footprint, footprint, -footprint, footprint, 0.5)) { + psError(PS_ERR_UNKNOWN, false, "Unable to reconvolve bad pixels from mask 1."); + psFree(mask); + return NULL; + } + if (!psImageConvolveMaskFFT(mask, mask, PM_SUBTRACTION_MASK_CONVOLVE_2, + PM_SUBTRACTION_MASK_FOOTPRINT_2, + -footprint, footprint, -footprint, footprint, 0.5)) { + psError(PS_ERR_UNKNOWN, false, "Unable to reconvolve bad pixels from mask 2."); + psFree(mask); + return NULL; + } + } else { + if (!psImageConvolveMaskDirect(mask, mask, PM_SUBTRACTION_MASK_BAD_1, + PM_SUBTRACTION_MASK_CONVOLVE_1 | PM_SUBTRACTION_MASK_FOOTPRINT_1, + -size, size, -size, size)) { + psError(PS_ERR_UNKNOWN, false, "Unable to convolve bad pixels from mask 1."); + psFree(mask); + return NULL; + } + if (!psImageConvolveMaskDirect(mask, mask, PM_SUBTRACTION_MASK_BAD_2, + PM_SUBTRACTION_MASK_CONVOLVE_2 | PM_SUBTRACTION_MASK_FOOTPRINT_2, + -size, size, -size, size)) { + psError(PS_ERR_UNKNOWN, false, "Unable to convolve bad pixels from mask 2."); + psFree(mask); + return NULL; + } + if (!psImageConvolveMaskDirect(mask, mask, PM_SUBTRACTION_MASK_CONVOLVE_1, + PM_SUBTRACTION_MASK_FOOTPRINT_1, + -footprint, footprint, -footprint, footprint)) { + psError(PS_ERR_UNKNOWN, false, "Unable to reconvolve bad pixels from mask 1."); + psFree(mask); + return NULL; + } + if (!psImageConvolveMaskDirect(mask, mask, PM_SUBTRACTION_MASK_CONVOLVE_2, + PM_SUBTRACTION_MASK_FOOTPRINT_2, + -footprint, footprint, -footprint, footprint)) { + psError(PS_ERR_UNKNOWN, false, "Unable to reconvolve bad pixels from mask 2."); + psFree(mask); + return NULL; + } + } + + return mask; +} + + +bool pmSubtractionBorder(psImage *image, psImage *weight, psImage *mask, + int size, psMaskType blank) +{ + PS_ASSERT_IMAGE_NON_NULL(image, false); + PS_ASSERT_IMAGE_TYPE(image, PS_TYPE_F32, false); + if (mask) { + PS_ASSERT_IMAGE_NON_NULL(mask, false); + PS_ASSERT_IMAGES_SIZE_EQUAL(mask, image, false); + PS_ASSERT_IMAGE_TYPE(mask, PS_TYPE_MASK, false); + } + if (weight) { + PS_ASSERT_IMAGE_NON_NULL(weight, false); + PS_ASSERT_IMAGES_SIZE_EQUAL(weight, image, false); + PS_ASSERT_IMAGE_TYPE(weight, PS_TYPE_F32, false); + } + + int numCols = image->numCols, numRows = image->numRows; // Image dimensions + + for (int y = size; y < numRows - size; y++) { + for (int x = 0; x < size; x++) { + markBlank(image, mask, weight, x, y, blank); + } + for (int x = numCols - size; x < numCols; x++) { + markBlank(image, mask, weight, x, y, blank); + } + } + for (int y = 0; y < size; y++) { + for (int x = 0; x < numCols; x++) { + markBlank(image, mask, weight, x, y, blank); + } + } + for (int y = numRows - size; y < numRows; y++) { + for (int x = 0; x < numCols; x++) { + markBlank(image, mask, weight, x, y, blank); + } + } + + return true; +} Index: trunk/psModules/src/imcombine/pmSubtractionMask.h =================================================================== --- trunk/psModules/src/imcombine/pmSubtractionMask.h (revision 15756) +++ trunk/psModules/src/imcombine/pmSubtractionMask.h (revision 15756) @@ -0,0 +1,25 @@ +#ifndef PM_SUBTRACTION_MASK_H +#define PM_SUBTRACTION_MASK_H + +#include + +/// Generate a mask for use in the subtraction process +psImage *pmSubtractionMask(const psImage *refMask, ///< Mask for the reference image (will be convolved) + const psImage *inMask, ///< Mask for the input image, or NULL + psMaskType maskVal, ///< Value to mask out + int size, ///< Half-size of the kernel (pmSubtractionKernels.size) + int footprint, ///< Half-size of the kernel footprint + float badFrac, ///< Maximum fraction of bad input pixels to accept + bool useFFT ///< Use FFT to do convolution? + ); + +/// Mark the non-convolved part of the image as blank +bool pmSubtractionBorder(psImage *image,///< Image + psImage *weight, ///< Weight map (or NULL) + psImage *mask, ///< Mask (or NULL) + int size, ///< Kernel half-size + psMaskType blank ///< Mask value for blank regions + ); + + +#endif Index: trunk/psModules/src/imcombine/pmSubtractionMatch.c =================================================================== --- trunk/psModules/src/imcombine/pmSubtractionMatch.c (revision 15745) +++ trunk/psModules/src/imcombine/pmSubtractionMatch.c (revision 15756) @@ -13,7 +13,11 @@ #include "pmSubtractionKernels.h" #include "pmSubtractionStamps.h" +#include "pmSubtractionEquation.h" #include "pmSubtraction.h" +#include "pmSubtractionMask.h" #include "pmSubtractionMatch.h" +#define TOL 1.0e-3 // Tolerance for subtraction solution +#define KERNEL_MOSAIC 2 // Half-number of kernel instances in the mosaic image static bool useFFT = true; // Do convolutions using FFT @@ -63,5 +67,6 @@ { psTrace("psModules.imcombine", 3, "Finding stamps...\n"); - *stamps = pmSubtractionStampsFind(*stamps, ro1->image, subMask, region, threshold, stampSpacing, mode); + *stamps = pmSubtractionStampsFind(*stamps, ro1->image, subMask, region, threshold, footprint, + stampSpacing, mode); if (!*stamps) { psError(PS_ERR_UNKNOWN, false, "Unable to find stamps."); @@ -72,5 +77,5 @@ psTrace("psModules.imcombine", 3, "Extracting stamps...\n"); - if (!pmSubtractionStampsExtract(*stamps, ro1->image, ro2 ? ro2->image : NULL, weight, footprint, size)) { + if (!pmSubtractionStampsExtract(*stamps, ro1->image, ro2 ? ro2->image : NULL, weight, size)) { psError(PS_ERR_UNKNOWN, false, "Unable to extract stamps."); return false; @@ -86,5 +91,5 @@ -bool pmSubtractionMatch(pmReadout *convolved, const pmReadout *ro1, const pmReadout *ro2, +bool pmSubtractionMatch(pmReadout *conv1, pmReadout *conv2, const pmReadout *ro1, const pmReadout *ro2, int footprint, float regionSize, float stampSpacing, float threshold, const psArray *sources, const char *stampsName, @@ -95,5 +100,8 @@ psMaskType maskBlank, float badFrac, pmSubtractionMode mode) { - PS_ASSERT_PTR_NON_NULL(convolved, false); + PS_ASSERT_PTR_NON_NULL(conv1, false); + if (mode == PM_SUBTRACTION_MODE_DUAL) { + PS_ASSERT_PTR_NON_NULL(conv2, false); + } PS_ASSERT_PTR_NON_NULL(ro1, false); PS_ASSERT_IMAGE_NON_NULL(ro1->image, false); @@ -174,15 +182,15 @@ // Reset the output readout, just in case - if (convolved->image) { - psFree(convolved->image); - convolved->image = NULL; - } - if (convolved->mask) { - psFree(convolved->mask); - convolved->mask = NULL; - } - if (convolved->weight) { - psFree(convolved->weight); - convolved->weight = NULL; + if (conv1->image) { + psFree(conv1->image); + conv1->image = NULL; + } + if (conv1->mask) { + psFree(conv1->mask); + conv1->mask = NULL; + } + if (conv1->weight) { + psFree(conv1->weight); + conv1->weight = NULL; } @@ -250,7 +258,9 @@ if (sources) { - stamps = pmSubtractionStampsSetFromSources(sources, subMask, region, stampSpacing, mode); + stamps = pmSubtractionStampsSetFromSources(sources, subMask, region, footprint, + stampSpacing, mode); } else if (stampsName && strlen(stampsName) > 0) { - stamps = pmSubtractionStampsSetFromFile(stampsName, subMask, region, stampSpacing, mode); + stamps = pmSubtractionStampsSetFromFile(stampsName, subMask, region, footprint, + stampSpacing, mode); } @@ -295,7 +305,7 @@ // Not an ISIS/GUNK kernel, or the optimum kernel search failed kernels = pmSubtractionKernelsGenerate(type, size, spatialOrder, isisWidths, isisOrders, - inner, binning, ringsOrder); - } - psMetadataAddPtr(convolved->analysis, PS_LIST_TAIL, "SUBTRACTION.KERNEL", + inner, binning, ringsOrder, mode); + } + psMetadataAddPtr(conv1->analysis, PS_LIST_TAIL, "SUBTRACTION.KERNEL", PS_DATA_UNKNOWN | PS_META_DUPLICATE_OK, "Subtraction kernels", kernels); @@ -312,5 +322,5 @@ psTrace("psModules.imcombine", 3, "Calculating equation...\n"); - if (!pmSubtractionCalculateEquation(stamps, kernels, footprint, mode)) { + if (!pmSubtractionCalculateEquation(stamps, kernels)) { psError(PS_ERR_UNKNOWN, false, "Unable to calculate least-squares equation."); goto MATCH_ERROR; @@ -320,6 +330,6 @@ psTrace("psModules.imcombine", 3, "Solving equation...\n"); - solution = pmSubtractionSolveEquation(solution, stamps); - if (!solution) { + + if (!pmSubtractionSolveEquation(kernels, stamps, TOL)) { psError(PS_ERR_UNKNOWN, false, "Unable to calculate least-squares equation."); goto MATCH_ERROR; @@ -328,6 +338,5 @@ memCheck(" solve equation"); - psVector *deviations = pmSubtractionCalculateDeviations(stamps, solution, footprint, kernels, - mode); // Deviations for each stamp + psVector *deviations = pmSubtractionCalculateDeviations(stamps, kernels); // Stamp deviations if (!deviations) { psError(PS_ERR_UNKNOWN, false, "Unable to calculate deviations."); @@ -351,11 +360,9 @@ if (numRejected > 0) { psTrace("psModules.imcombine", 3, "Solving equation...\n"); - solution = pmSubtractionSolveEquation(solution, stamps); - if (!solution) { + if (!pmSubtractionSolveEquation(kernels, stamps, TOL)) { psError(PS_ERR_UNKNOWN, false, "Unable to calculate least-squares equation."); goto MATCH_ERROR; } - psVector *deviations = pmSubtractionCalculateDeviations(stamps, solution, footprint, kernels, - mode); // Deviations for each stamp + psVector *deviations = pmSubtractionCalculateDeviations(stamps, kernels); // Stamp deviations if (!deviations) { psError(PS_ERR_UNKNOWN, false, "Unable to calculate deviations."); @@ -376,8 +383,9 @@ psImage *convKernels = psImageAlloc(5 * fullSize - 1, 5 * fullSize - 1, PS_TYPE_F32); psImageInit(convKernels, NAN); - for (int j = -2; j <= 2; j++) { - for (int i = -2; i <= 2; i++) { - psImage *kernel = pmSubtractionKernelImage(solution, kernels, (float)i / 2.0, - (float)j / 2.0); // Image of the kernel + for (int j = -KERNEL_MOSAIC; j <= KERNEL_MOSAIC; j++) { + for (int i = -KERNEL_MOSAIC; i <= KERNEL_MOSAIC; i++) { + psImage *kernel = pmSubtractionKernelImage(kernels, (float)i / (float)KERNEL_MOSAIC, + (float)j / (float)KERNEL_MOSAIC, + false); // Image of the kernel if (!kernel) { psError(PS_ERR_UNKNOWN, false, "Unable to generate kernel image."); @@ -386,6 +394,6 @@ } - if (psImageOverlaySection(convKernels, kernel, (i + 2) * fullSize, - (j + 2) * fullSize, "=") == 0) { + if (psImageOverlaySection(convKernels, kernel, (i + KERNEL_MOSAIC) * fullSize, + (j + KERNEL_MOSAIC) * fullSize, "=") == 0) { psError(PS_ERR_UNKNOWN, false, "Unable to overlay kernel image."); psFree(kernel); @@ -399,5 +407,5 @@ psString comment = NULL; // Comment for metadata psStringAppend(&comment, "Subtraction kernel for region %s", regionString); - psMetadataAddImage(convolved->analysis, PS_LIST_TAIL, "SUBTRACTION.KERNEL.IMAGE", + psMetadataAddImage(conv1->analysis, PS_LIST_TAIL, "SUBTRACTION.KERNEL.IMAGE", PS_META_DUPLICATE_OK, comment, convKernels); psFree(comment); @@ -427,6 +435,5 @@ psTrace("psModules.imcombine", 2, "Convolving...\n"); - if (!pmSubtractionConvolve(convolved, ro1, ro2, subMask, maskBlank, region, - solution, kernels, mode, useFFT)) { + if (!pmSubtractionConvolve(conv1, conv2, ro1, ro2, subMask, maskBlank, region, kernels, useFFT)) { psError(PS_ERR_UNKNOWN, false, "Unable to convolve image."); goto MATCH_ERROR; @@ -439,13 +446,13 @@ psString comment = NULL; // Comment for metadata psStringAppend(&comment, "Subtraction solution for region %s", regionString); - psMetadataAddVector(convolved->analysis, PS_LIST_TAIL, "SUBTRACTION.SOLUTION", + psMetadataAddVector(conv1->analysis, PS_LIST_TAIL, "SUBTRACTION.SOLUTION", PS_META_DUPLICATE_OK, comment, solution); psFree(comment); if (region) { - psMetadataAddPtr(convolved->analysis, PS_LIST_TAIL, "SUBTRACTION.REGION", + psMetadataAddPtr(conv1->analysis, PS_LIST_TAIL, "SUBTRACTION.REGION", PS_META_DUPLICATE_OK | PS_DATA_REGION, comment, region); } else { region = psRegionAlloc(0, numCols, 0, numRows); - psMetadataAddPtr(convolved->analysis, PS_LIST_TAIL, "SUBTRACTION.REGION", + psMetadataAddPtr(conv1->analysis, PS_LIST_TAIL, "SUBTRACTION.REGION", PS_META_DUPLICATE_OK | PS_DATA_REGION, comment, region); psFree(region); @@ -458,8 +465,15 @@ // There is data in the readout now - convolved->data_exists = true; - if (convolved->parent) { - convolved->parent->data_exists = true; - convolved->parent->parent->data_exists = true; + conv1->data_exists = true; + if (conv1->parent) { + conv1->parent->data_exists = true; + conv1->parent->parent->data_exists = true; + } + if (mode == PM_SUBTRACTION_MODE_DUAL) { + conv2->data_exists = true; + if (conv2->parent) { + conv2->parent->data_exists = true; + conv2->parent->parent->data_exists = true; + } } } @@ -474,5 +488,5 @@ weight = NULL; - if (!pmSubtractionBorder(convolved->image, convolved->weight, convolved->mask, size, maskBlank)) { + if (!pmSubtractionBorder(conv1->image, conv1->weight, conv1->mask, size, maskBlank)) { psError(PS_ERR_UNKNOWN, false, "Unable to set border of convolved image."); goto MATCH_ERROR; Index: trunk/psModules/src/imcombine/pmSubtractionMatch.h =================================================================== --- trunk/psModules/src/imcombine/pmSubtractionMatch.h (revision 15745) +++ trunk/psModules/src/imcombine/pmSubtractionMatch.h (revision 15756) @@ -10,5 +10,6 @@ /// Match two images -bool pmSubtractionMatch(pmReadout *convolved, ///< Output convolved data +bool pmSubtractionMatch(pmReadout *conv1, ///< Output convolved data for image 1 + pmReadout *conv2, ///< Output convolved data for image 2 const pmReadout *ro1, ///< Image 1 const pmReadout *ro2, ///< Image 2 Index: trunk/psModules/src/imcombine/pmSubtractionParams.c =================================================================== --- trunk/psModules/src/imcombine/pmSubtractionParams.c (revision 15745) +++ trunk/psModules/src/imcombine/pmSubtractionParams.c (revision 15756) @@ -232,5 +232,5 @@ psVectorInit(orders, maxOrder); pmSubtractionKernels *kernels = p_pmSubtractionKernelsRawISIS(size, spatialOrder, - fwhms, orders); // Kernels + fwhms, orders, mode); // Kernels psFree(orders); psFree(kernels->description); @@ -279,5 +279,5 @@ continue; } - if (!pmSubtractionConvolveStamp(stamp, kernels, footprint, mode)) { + if (!pmSubtractionConvolveStamp(stamp, kernels, footprint)) { psError(PS_ERR_UNKNOWN, false, "Unable to convolve stamp %d.", i); psFree(targets); @@ -505,8 +505,4 @@ kernels->type = PM_SUBTRACTION_KERNEL_GUNK; - - kernels->subIndex = PS_SQR(2 * inner + 1) / 2 + numGaussians; - assert(kernels->u->data.S32[kernels->subIndex] == 0 && - kernels->v->data.S32[kernels->subIndex] == 0); } Index: trunk/psModules/src/imcombine/pmSubtractionStamps.c =================================================================== --- trunk/psModules/src/imcombine/pmSubtractionStamps.c (revision 15745) +++ trunk/psModules/src/imcombine/pmSubtractionStamps.c (revision 15756) @@ -44,6 +44,4 @@ ) { - psFree(stamp->matrix); - psFree(stamp->vector); psFree(stamp->image1); psFree(stamp->image2); @@ -51,4 +49,10 @@ psFree(stamp->convolutions1); psFree(stamp->convolutions2); + + psFree(stamp->matrix1); + psFree(stamp->matrix2); + psFree(stamp->matrixX); + psFree(stamp->vector); + } @@ -103,5 +107,5 @@ pmSubtractionStampList *pmSubtractionStampListAlloc(int numCols, int numRows, const psRegion *region, - float spacing) + int footprint, float spacing) { pmSubtractionStampList *list = psAlloc(sizeof(pmSubtractionStampList)); // Stamp list to return @@ -124,14 +128,16 @@ for (int y = 0, index = 0; y < yStamps; y++) { - int yStart = yMin + y * ySize / (yStamps + 1); // Subregion starts here - int yStop = yMin + (y + 1) * ySize / (yStamps + 1) - 1; // Subregion stops here + int yStart = yMin + y * ((float)ySize / (float)(yStamps + 1)); // Subregion starts here + int yStop = yMin + (y + 1) * ((float)ySize / (float)(yStamps + 1)) - 1; // Subregion stops here assert(yStart >= yMin && yStop < yMax); for (int x = 0; x < xStamps; x++, index++) { - int xStart = xMin + x * xSize / (xStamps + 1); // Subregion starts here - int xStop = xMin + (x + 1) * xSize / (xStamps + 1) - 1; // Subregion stops here + int xStart = xMin + x * ((float)xSize / (float)(xStamps + 1)); // Subregion starts here + int xStop = xMin + (x + 1) * ((float)xSize / (float)(xStamps + 1)) - 1; // Subregion stops here assert(xStart >= xMin && xStop < xMax); list->stamps->data[index] = pmSubtractionStampAlloc(); + psTrace("psModules.imcombine", 6, "Stamp region %d: [%d:%d,%d:%d]\n", + index, xStart, xStop, yStart, yStop); list->regions->data[index] = psRegionAlloc(xStart, xStop, yStart, yStop); } @@ -141,4 +147,5 @@ list->y = NULL; list->flux = NULL; + list->footprint = footprint; return list; @@ -155,6 +162,4 @@ stamp->xNorm = NAN; stamp->yNorm = NAN; - stamp->matrix = NULL; - stamp->vector = NULL; stamp->status = PM_SUBTRACTION_STAMP_INIT; @@ -165,4 +170,9 @@ stamp->convolutions2 = NULL; + stamp->matrix1 = NULL; + stamp->matrix2 = NULL; + stamp->matrixX = NULL; + stamp->vector = NULL; + return stamp; } @@ -171,5 +181,6 @@ pmSubtractionStampList *pmSubtractionStampsFind(pmSubtractionStampList *stamps, const psImage *image, const psImage *subMask, const psRegion *region, - float threshold, float spacing, pmSubtractionMode mode) + float threshold, int footprint, float spacing, + pmSubtractionMode mode) { PS_ASSERT_IMAGE_NON_NULL(image, NULL); @@ -180,4 +191,5 @@ PS_ASSERT_IMAGE_TYPE(subMask, PS_TYPE_MASK, NULL); } + PS_ASSERT_INT_NONNEGATIVE(footprint, NULL); PS_ASSERT_FLOAT_LARGER_THAN(spacing, 0.0, NULL); if (region) { @@ -201,5 +213,5 @@ if (!stamps) { - stamps = pmSubtractionStampListAlloc(numCols, numRows, region, spacing); + stamps = pmSubtractionStampListAlloc(numCols, numRows, region, footprint, spacing); } @@ -244,6 +256,6 @@ fluxStamp = threshold; psRegion *subRegion = stamps->regions->data[i]; // Sub-region of interest - for (int y = subRegion->y0; y <= subRegion->y1 ; y++) { - for (int x = subRegion->x0; x <= subRegion->x1 ; x++) { + for (int y = subRegion->y0; y <= subRegion->y1; y++) { + for (int x = subRegion->x0; x <= subRegion->x1; x++) { if (checkStampMask(x, y, subMask, mode) && image->data.F32[y][x] > fluxStamp) { fluxStamp = image->data.F32[y][x]; @@ -289,5 +301,6 @@ pmSubtractionStampList *pmSubtractionStampsSet(const psVector *x, const psVector *y, const psVector *flux, const psImage *image, const psImage *subMask, - const psRegion *region, float spacing, pmSubtractionMode mode) + const psRegion *region, int footprint, float spacing, + pmSubtractionMode mode) { @@ -316,5 +329,5 @@ int numStars = x->n; // Number of stars pmSubtractionStampList *stamps = pmSubtractionStampListAlloc(subMask->numCols, subMask->numRows, - region, spacing); // Stamp list + region, footprint, spacing); // Stamp list int numStamps = stamps->num; // Number of stamps @@ -401,5 +414,5 @@ bool pmSubtractionStampsExtract(pmSubtractionStampList *stamps, psImage *image1, psImage *image2, - psImage *weight, int footprint, int kernelSize) + psImage *weight, int kernelSize) { PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(stamps, false); @@ -414,9 +427,8 @@ PS_ASSERT_IMAGES_SIZE_EQUAL(weight, image1, false); PS_ASSERT_IMAGE_TYPE(weight, PS_TYPE_F32, false); - PS_ASSERT_INT_NONNEGATIVE(footprint, false); PS_ASSERT_INT_NONNEGATIVE(kernelSize, false); int numCols = image1->numCols, numRows = image1->numRows; // Size of images - int size = kernelSize + footprint; // Size of postage stamps + int size = kernelSize + stamps->footprint; // Size of postage stamps for (int i = 0; i < stamps->num; i++) { @@ -467,12 +479,11 @@ #if 0 -bool pmSubtractionStampsGenerate(pmSubtractionStampList *stamps, float fwhm, int footprint, int kernelSize) +bool pmSubtractionStampsGenerate(pmSubtractionStampList *stamps, float fwhm, int kernelSize) { PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(stamps, false); PS_ASSERT_FLOAT_LARGER_THAN(fwhm, 0.0, false); - PS_ASSERT_INT_NONNEGATIVE(footprint, false); PS_ASSERT_INT_NONNEGATIVE(kernelSize, false); - int size = kernelSize + footprint; // Size of postage stamps + int size = kernelSize + stamps->footprint; // Size of postage stamps int num = stamps->num; // Number of stamps float sigma = fwhm / (2.0 * sqrtf(2.0 * log(2.0))); // Gaussian sigma @@ -514,6 +525,6 @@ pmSubtractionStampList *pmSubtractionStampsSetFromSources(const psArray *sources, const psImage *subMask, - const psRegion *region, float spacing, - pmSubtractionMode mode) + const psRegion *region, int footprint, + float spacing, pmSubtractionMode mode) { PS_ASSERT_ARRAY_NON_NULL(sources, NULL); @@ -539,5 +550,5 @@ pmSubtractionStampList *stamps = pmSubtractionStampsSet(x, y, flux, NULL, subMask, region, - spacing, mode); // Stamps to return + footprint, spacing, mode); // Stamps to return psFree(x); psFree(y); @@ -553,5 +564,5 @@ pmSubtractionStampList *pmSubtractionStampsSetFromFile(const char *filename, const psImage *subMask, - const psRegion *region, float spacing, + const psRegion *region, int footprint, float spacing, pmSubtractionMode mode) { @@ -572,6 +583,6 @@ psBinaryOp(y, y, "-", psScalarAlloc(1.0, PS_TYPE_F32)); - pmSubtractionStampList *stamps = pmSubtractionStampsSet(x, y, flux, NULL, subMask, region, spacing, - mode); + pmSubtractionStampList *stamps = pmSubtractionStampsSet(x, y, flux, NULL, subMask, region, footprint, + spacing, mode); psFree(data); Index: trunk/psModules/src/imcombine/pmSubtractionStamps.h =================================================================== --- trunk/psModules/src/imcombine/pmSubtractionStamps.h (revision 15745) +++ trunk/psModules/src/imcombine/pmSubtractionStamps.h (revision 15756) @@ -23,4 +23,5 @@ psArray *x, *y; ///< Coordinates for possible stamps (or NULL) psArray *flux; ///< Fluxes for possible stamps (or NULL) + int footprint; ///< Half-size of stamps } pmSubtractionStampList; @@ -29,4 +30,5 @@ int numRows, // Number of rows in image const psRegion *region, // Region for stamps, or NULL + int footprint, // Half-size of stamps float spacing // Rough average spacing between stamps ); @@ -40,4 +42,5 @@ PS_ASSERT_ARRAY_SIZE((LIST)->stamps, (LIST)->num, RETURNVALUE); \ PS_ASSERT_ARRAY_SIZE((LIST)->regions, (LIST)->num, RETURNVALUE); \ + PS_ASSERT_INT_NONNEGATIVE((LIST)->footprint, RETURNVALUE); \ if ((LIST)->x || (LIST)->y || (LIST)->flux) { \ PS_ASSERT_ARRAY_NON_NULL((LIST)->x, RETURNVALUE); \ @@ -57,9 +60,10 @@ psKernel *image1; ///< Reference image postage stamp psKernel *image2; ///< Input image postage stamp - psKernel *weight; ///< Weight image postage stamp - psArray *convolutions1; ///< Convolutions of image 1 for each kernel component - psArray *convolutions2; ///< Convolutions of image 2 for each kernel component - psImage *matrix; ///< Associated matrix - psVector *vector; ///< Assoicated vector + psKernel *weight; ///< Weight 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; ///< Matrices for each image, or NULL + psImage *matrixX; ///< Cross-matrix (for mode DUAL), or NULL + psVector *vector; ///< Associated vector (when mode not DUAL), or NULL pmSubtractionStampStatus status; ///< Status of stamp } pmSubtractionStamp; @@ -74,4 +78,5 @@ const psRegion *region, ///< Region to search, or NULL float threshold, ///< Threshold for stamps in the image + int footprint, ///< Half-size for stamps float spacing, ///< Rough spacing for stamps pmSubtractionMode mode ///< Mode for subtraction @@ -85,4 +90,5 @@ const psImage *mask, ///< Mask, or NULL const psRegion *region, ///< Region to search, or NULL + int footprint, ///< Half-size for stamps float spacing, ///< Rough spacing for stamps pmSubtractionMode mode ///< Mode for subtraction @@ -94,4 +100,5 @@ const psImage *subMask, ///< Mask, or NULL const psRegion *region, ///< Region to search, or NULL + int footprint, ///< Half-size for stamps float spacing, ///< Rough spacing for stamps pmSubtractionMode mode ///< Mode for subtraction @@ -103,4 +110,5 @@ const psImage *subMask, ///< Mask, or NULL const psRegion *region, ///< Region to search, or NULL + int footprint, ///< Half-size for stamps float spacing, ///< Rough spacing for stamps pmSubtractionMode mode ///< Mode for subtraction @@ -109,8 +117,7 @@ /// Extract stamps from the images bool pmSubtractionStampsExtract(pmSubtractionStampList *stamps, ///< Stamps - psImage *reference, ///< Reference image - psImage *input, ///< Input image (or NULL) + psImage *image1, ///< Reference image + psImage *image2, ///< Input image (or NULL) psImage *weight, ///< Weight (variance) map - int footprint, ///< Stamp footprint size int kernelSize ///< Kernel half-size );