Index: trunk/psModules/src/imcombine/pmSubtraction.c =================================================================== --- trunk/psModules/src/imcombine/pmSubtraction.c (revision 14456) +++ 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; } }