Index: trunk/psModules/src/imcombine/pmSubtraction.c =================================================================== --- trunk/psModules/src/imcombine/pmSubtraction.c (revision 14477) +++ trunk/psModules/src/imcombine/pmSubtraction.c (revision 14480) @@ -4,6 +4,6 @@ * @author GLG, MHPCC * - * @version $Revision: 1.35 $ $Name: not supported by cvs2svn $ - * @date $Date: 2007-08-09 21:20:35 $ + * @version $Revision: 1.36 $ $Name: not supported by cvs2svn $ + * @date $Date: 2007-08-14 02:14:30 $ * * Copyright 2004-2007 Institute for Astronomy, University of Hawaii @@ -217,5 +217,5 @@ int index, // Kernel basis function index int x, int y, // Pixel around which to convolve - const psImage *image // Image to convolve + const psKernel *image // Image to convolve (a kernel for convenience) ) { @@ -225,8 +225,8 @@ int u = kernels->u->data.S32[index]; // Offset in x int v = kernels->v->data.S32[index]; // Offset in y - float value = image->data.F32[y + v][x + u]; // Value of convolution + float value = image->kernel[y + v][x + u]; // Value of convolution if (kernels->spatialOrder > 0 && index != kernels->subIndex) { // The (0,0) element is subtracted from most kernels to preserve photometric scaling - value -= image->data.F32[y][x]; + value -= image->kernel[y][x]; } return value; @@ -242,5 +242,5 @@ for (int v = vStart; v <= vStop; v++) { for (int u = uStart; u <= uStop; u++) { - sum += image->data.F32[y + v][x + u]; + sum += image->kernel[y + v][x + u]; } } @@ -248,5 +248,5 @@ if (kernels->spatialOrder > 0 && index != kernels->subIndex) { // The (0,0) element is subtracted from most kernels to preserve photometric scaling - sum -= image->data.F32[y][x]; + sum -= image->kernel[y][x]; } return sum; @@ -260,5 +260,5 @@ for (int v = -size; v <= size; v++) { for (int u = -size; u <= size; u++) { - sum += kernel->kernel[v][u] * image->data.F32[y + v][x + u]; + sum += kernel->kernel[v][u] * image->kernel[y + v][x + u]; } } @@ -268,8 +268,8 @@ int u = kernels->u->data.S32[index]; // Offset in x int v = kernels->v->data.S32[index]; // Offset in y - float value = image->data.F32[y + v][x + u]; // Value of convolution + float value = image->kernel[y + v][x + u]; // Value of convolution // The (0,0) delta function is subtracted from most kernels to preserve photometric scaling if (kernels->spatialOrder > 0 && index != kernels->subIndex) { - value -= image->data.F32[y][x]; + value -= image->kernel[y][x]; } return value; @@ -281,5 +281,5 @@ for (int v = -size; v <= size; v++) { for (int u = -size; u <= size; u++) { - sum += kernel->kernel[v][u] * image->data.F32[y + v][x + u]; + sum += kernel->kernel[v][u] * image->kernel[y + v][x + u]; // Photometric scaling is already built in to the precalculated kernel } @@ -289,5 +289,5 @@ case PM_SUBTRACTION_KERNEL_RINGS: { if (index == kernels->subIndex) { - return image->data.F32[y][x]; + return image->kernel[y][x]; } psArray *preCalc = kernels->preCalc->data[index]; // Precalculated data @@ -299,8 +299,8 @@ for (int j = 0; j < num; j++) { int u = uCoords->data.S32[j], v = vCoords->data.S32[j]; // Kernel coordinates - sum += image->data.F32[y + v][x + u] * poly->data.F32[j]; + sum += image->kernel[y + v][x + u] * poly->data.F32[j]; } // The (0,0) kernel is subtracted from other kernels to preserve photometric scaling - return sum - num * image->data.F32[y][x]; + return sum - num * image->kernel[y][x]; } default: @@ -436,25 +436,10 @@ -bool pmSubtractionCalculateEquation(psArray *stamps, const psImage *reference, const psImage *input, - const psImage *weight, const pmSubtractionKernels *kernels, int footprint) +bool pmSubtractionCalculateEquation(psArray *stamps, const pmSubtractionKernels *kernels, int footprint) { PS_ASSERT_ARRAY_NON_NULL(stamps, false); - PS_ASSERT_IMAGE_NON_NULL(reference, false); - PS_ASSERT_IMAGE_TYPE(reference, PS_TYPE_F32, false); - PS_ASSERT_IMAGE_NON_NULL(input, false); - PS_ASSERT_IMAGE_TYPE(input, PS_TYPE_F32, false); - PS_ASSERT_IMAGES_SIZE_EQUAL(reference, input, false); PS_ASSERT_PTR_NON_NULL(kernels, false); PS_ASSERT_VECTORS_SIZE_EQUAL(kernels->u, kernels->v, false); - if (weight) { - PS_ASSERT_IMAGE_NON_NULL(weight, false); - PS_ASSERT_IMAGES_SIZE_EQUAL(weight, input, false); - PS_ASSERT_IMAGE_TYPE(weight, PS_TYPE_F32, false); - } PS_ASSERT_INT_NONNEGATIVE(footprint, false); - - // Image size - int numCols = reference->numCols; - int numRows = reference->numRows; int spatialOrder = kernels->spatialOrder; // Maximum order of spatial variation @@ -491,12 +476,14 @@ psVectorInit(stampVector, 0.0); - float xNorm = 2.0 * (float)(stamp->x - numCols/2.0) / (float)numCols; // Normalised x coord - float yNorm = 2.0 * (float)(stamp->y - numRows/2.0) / (float)numRows; // Normalised y coord - psImage *polyValues = spatialPolyValues(spatialOrder, xNorm, yNorm); // Spatial polynomial terms - - for (int y = stamp->y - footprint; y <= stamp->y + footprint; y++) { - for (int x = stamp->x - footprint; x <= stamp->x + footprint; x++) { - float invNoise2 = 1.0 / (weight ? weight->data.F32[y][x] : - input->data.F32[y][x]); // Inverse square noise + // Spatial polynomial terms + psImage *polyValues = spatialPolyValues(spatialOrder, stamp->xNorm, stamp->yNorm); + + psKernel *reference = stamp->reference; // Reference postage stamp + psKernel *input = stamp->input; // Input postage stamp + psKernel *weight = stamp->weight; // Weight map postage stamp + + for (int y = - footprint; y <= footprint; y++) { + for (int x = - footprint; x <= footprint; x++) { + float invNoise2 = 1.0 / weight->kernel[y][x]; // Inverse square noise // Generate the convolutions @@ -519,5 +506,5 @@ convolutions->data.F64[j] * invNoise2; } - stampVector->data.F64[i] += input->data.F32[y][x] * convolutions->data.F64[i] * + stampVector->data.F64[i] += input->kernel[y][x] * convolutions->data.F64[i] * invNoise2; @@ -527,5 +514,5 @@ // Background only terms stampMatrix->data.F64[bgIndex][bgIndex] += invNoise2; - stampVector->data.F64[bgIndex] += input->data.F32[y][x] * invNoise2; + stampVector->data.F64[bgIndex] += input->kernel[y][x] * invNoise2; } } Index: trunk/psModules/src/imcombine/pmSubtraction.h =================================================================== --- trunk/psModules/src/imcombine/pmSubtraction.h (revision 14477) +++ trunk/psModules/src/imcombine/pmSubtraction.h (revision 14480) @@ -6,6 +6,6 @@ * @author GLG, MHPCC * - * @version $Revision: 1.6 $ $Name: not supported by cvs2svn $ - * @date $Date: 2007-08-09 20:25:52 $ + * @version $Revision: 1.7 $ $Name: not supported by cvs2svn $ + * @date $Date: 2007-08-14 02:14:30 $ * Copyright 2004-207 Institute for Astronomy, University of Hawaii */ @@ -40,7 +40,4 @@ /// Calculate the least-squares equation to match the image quality bool pmSubtractionCalculateEquation(psArray *stamps, ///< The stamps for which to calculate the equation, - const psImage *reference, ///< Reference image - const psImage *input, ///< Input image - const psImage *weight, ///< Weight image, or NULL const pmSubtractionKernels *kernels, ///< Kernel parameters int footprint ///< Half-size of region over which to calculate equation Index: trunk/psModules/src/imcombine/pmSubtractionStamps.c =================================================================== --- trunk/psModules/src/imcombine/pmSubtractionStamps.c (revision 14477) +++ trunk/psModules/src/imcombine/pmSubtractionStamps.c (revision 14480) @@ -33,7 +33,13 @@ stamp->x = 0; stamp->y = 0; + stamp->xNorm = 0.0; + stamp->yNorm = 0.0; stamp->matrix = NULL; stamp->vector = NULL; stamp->status = status; + + stamp->reference = NULL; + stamp->input = NULL; + stamp->weight = NULL; return stamp; @@ -137,4 +143,13 @@ stamp->x = xBest; stamp->y = yBest; + stamp->xNorm = 2.0 * (float)(xBest - numCols/2.0) / (float)numCols; + stamp->yNorm = 2.0 * (float)(yBest - numRows/2.0) / (float)numRows; + + // Reset the postage stamps since we're making a new stamp + psFree(stamp->reference); + psFree(stamp->input); + psFree(stamp->weight); + stamp->reference = stamp->input = stamp->weight = NULL; + if (fluxBest > threshold) { stamp->status = PM_SUBTRACTION_STAMP_CALCULATE; @@ -154,2 +169,61 @@ return stamps; } + + +bool pmSubtractionExtractStamps(psArray *stamps, psImage *reference, psImage *input, psImage *weight, + int footprint, const pmSubtractionKernels *kernels) +{ + PS_ASSERT_ARRAY_NON_NULL(stamps, false); + PS_ASSERT_IMAGE_NON_NULL(reference, false); + PS_ASSERT_IMAGE_TYPE(reference, PS_TYPE_F32, false); + if (input) { + PS_ASSERT_IMAGE_NON_NULL(input, false); + PS_ASSERT_IMAGES_SIZE_EQUAL(input, reference, false); + PS_ASSERT_IMAGE_TYPE(input, PS_TYPE_F32, false); + } + if (weight) { + PS_ASSERT_IMAGE_NON_NULL(weight, false); + PS_ASSERT_IMAGES_SIZE_EQUAL(weight, reference, false); + PS_ASSERT_IMAGE_TYPE(weight, PS_TYPE_F32, false); + } + + int numCols = reference->numCols, numRows = reference->numRows; // Size of images + int size = kernels->size + footprint; // Size of postage stamps + + if (!weight) { + // Use the input as a rough approximation to the variance map, and HOPE that it's positive. + weight = input; + } + + for (int i = 0; i < stamps->n; i++) { + pmSubtractionStamp *stamp = stamps->data[i]; // Stamp of interest + + + int x = stamp->x, y = stamp->y; // Stamp coordinates + if (x < size || x > numCols - size || y < size || y > numRows - size) { + psError(PS_ERR_UNKNOWN, false, "Stamp %d (%d,%d) is within the image border.\n", i, x, y); + return false; + } + + psRegion region = psRegionSet(x - size, x + size + 1, y - size, y + size + 1); // Region of interest + + psImage *refSub = psImageSubset(reference, region); // Subimage with stamp + stamp->reference = psKernelAllocFromImage(refSub, size, size); + psFree(refSub); // Drop reference + + if (input) { + psImage *inSub = psImageSubset(input, region); // Subimage with stamp + stamp->input = psKernelAllocFromImage(inSub, size, size); + psFree(inSub); // Drop reference + } + + psImage *wtSub = psImageSubset(weight, region); // Subimage with stamp + stamp->weight = psKernelAllocFromImage(wtSub, size, size); + psFree(wtSub); // Drop reference + + } + + return true; +} + + Index: trunk/psModules/src/imcombine/pmSubtractionStamps.h =================================================================== --- trunk/psModules/src/imcombine/pmSubtractionStamps.h (revision 14477) +++ trunk/psModules/src/imcombine/pmSubtractionStamps.h (revision 14480) @@ -4,4 +4,5 @@ #include +#include "pmSubtractionKernels.h" /// Status of stamp @@ -17,4 +18,8 @@ typedef struct { int x, y; ///< Position + float xNorm, yNorm; ///< Normalised position + psKernel *reference; ///< Reference image postage stamp + psKernel *input; ///< Input image postage stamp + psKernel *weight; ///< Weight image postage stamp psImage *matrix; ///< Associated matrix psVector *vector; ///< Assoicated vector @@ -31,4 +36,12 @@ ); +/// Extract stamps from the images +bool pmSubtractionExtractStamps(psArray *stamps, ///< Stamps + psImage *reference, ///< Reference image + psImage *input, ///< Input image (or NULL) + psImage *weight, ///< Weight (variance) map + int footprint, ///< Stamp footprint size + const pmSubtractionKernels *kernels ///< Kernel basis functions (for size) + ); #endif