Index: /trunk/psModules/src/imcombine/pmSubtraction.c =================================================================== --- /trunk/psModules/src/imcombine/pmSubtraction.c (revision 14318) +++ /trunk/psModules/src/imcombine/pmSubtraction.c (revision 14319) @@ -4,6 +4,6 @@ * @author GLG, MHPCC * - * @version $Revision: 1.28 $ $Name: not supported by cvs2svn $ - * @date $Date: 2007-07-19 21:26:12 $ + * @version $Revision: 1.29 $ $Name: not supported by cvs2svn $ + * @date $Date: 2007-07-19 21:42:12 $ * * Copyright 2004-2007 Institute for Astronomy, University of Hawaii @@ -274,108 +274,4 @@ #endif -#ifndef USE_FUNCTIONS_INSTEAD_OF_MACROS -// Generate the convolved pixel value -// Meant to replace the following function declaration: -// static inline double convolvePixel(const pmSubtractionKernels *kernels, // Kernel basis functions -// int index, // Kernel basis function index -// int x, int y, // Pixel around which to convolve -// const psImage *image, // Image to convolve -// const psImage *polyValues, // Spatial polynomial values -// double (*weightFunc)(double value) // Function for weighting -// ); -// because even though weightFunc is "static inline", it still appears in the assembly. -// -// TARGET: The value to 'return' (double) -// KERNELS: Kernel basis functions (pmSubtractionKernels*) -// INDEX: Kernel basis function index (int) -// X, Y: Pixel around which to convolve (int) -// IMAGE: Image to convolve (psImage*) -// POLYVALUES: Spatial polynomial values (psImage*) -// FUNC: Function for weighting (double (*weightFunc)(double value)) -#define CONVOLVE_PIXEL(TARGET, KERNELS, INDEX, X, Y, IMAGE, POLYVALUES, FUNC) { \ - int xOrder = (KERNELS)->xOrder->data.S32[(INDEX)]; /* Polynomial order in x */ \ - int yOrder = (KERNELS)->yOrder->data.S32[(INDEX)]; /* Polynomial order in y */ \ - double polyValue = (POLYVALUES)->data.F64[yOrder][xOrder]; /* Value of spatial polynomial */ \ - \ - switch ((KERNELS)->type) { \ - case PM_SUBTRACTION_KERNEL_POIS: { \ - /* Convolution with a delta function is just the value specified by the offset */ \ - int u = (KERNELS)->u->data.S32[(INDEX)]; /* Offset in x */ \ - int v = (KERNELS)->v->data.S32[(INDEX)]; /* Offset in y */ \ - (TARGET) = FUNC(polyValue) * (IMAGE)->data.F32[(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 */ \ - (TARGET) += FUNC(-1.0) * (IMAGE)->data.F32[(Y)][(X)]; \ - } \ - break; \ - } \ - /* Method for SPAM and FRIES is the same */ \ - case PM_SUBTRACTION_KERNEL_SPAM: \ - case PM_SUBTRACTION_KERNEL_FRIES: { \ - int uStart = (KERNELS)->u->data.S32[(INDEX)]; \ - int uStop = (KERNELS)->uStop->data.S32[(INDEX)]; \ - int vStart = (KERNELS)->v->data.S32[(INDEX)]; \ - int vStop = (KERNELS)->vStop->data.S32[(INDEX)]; \ - double sum = 0.0; \ - for (int v = vStart; v <= vStop; v++) { \ - for (int u = uStart; u <= uStop; u++) { \ - sum += FUNC(polyValue) * (IMAGE)->data.F32[(Y) + v][(X) + u]; \ - } \ - } \ - if ((KERNELS)->spatialOrder > 0 && (INDEX) != (KERNELS)->subIndex) { \ - /* The (0,0) element is subtracted from most kernels to preserve photometric scaling */ \ - sum += FUNC(-1.0) * (IMAGE)->data.F32[(Y)][(X)]; \ - } \ - (TARGET) = sum; \ - break; \ - } \ - case PM_SUBTRACTION_KERNEL_GUNK: { \ - double value; /* The value to return */ \ - if ((INDEX) < (KERNELS)->inner) { \ - /* Using pre-calculated function */ \ - psKernel *kernel = (KERNELS)->preCalc->data[(INDEX)]; /* The convolution kernel */ \ - int size = (KERNELS)->size; /* Kernel half-size */ \ - double sum = 0.0; /* Accumulated sum from convolution */ \ - for (int v = -size; v <= size; v++) { \ - for (int u = -size; u <= size; u++) { \ - sum += FUNC(kernel->kernel[v][u]) * (IMAGE)->data.F32[(Y) + v][(X) + u]; \ - } \ - } \ - value = FUNC(polyValue) * sum; \ - } else { \ - /* Using delta function */ \ - int u = (KERNELS)->u->data.S32[(INDEX)]; /* Offset in x */ \ - int v = (KERNELS)->v->data.S32[(INDEX)]; /* Offset in y */ \ - value = FUNC(polyValue) * (IMAGE)->data.F32[(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 += FUNC(-1.0) * (IMAGE)->data.F32[(Y)][(X)]; \ - } \ - (TARGET) = value; \ - } \ - case PM_SUBTRACTION_KERNEL_ISIS: { \ - psKernel *kernel = (KERNELS)->preCalc->data[(INDEX)]; /* The convolution kernel */ \ - psKernel *subKernel = (KERNELS)->preCalc->data[(KERNELS)->subIndex]; /* Kernel to subtract */ \ - int size = (KERNELS)->size; /* Kernel half-size */ \ - double sum = 0.0; /* Accumulated sum from convolution */ \ - double sub = 0.0; /* Accumulated sum to subtract */ \ - for (int v = -size, yLocation = (Y) - size; v <= size; v++, yLocation++) { \ - for (int u = -size, xLocation = (X) - size; u <= size; u++, xLocation++) { \ - sum += FUNC(kernel->kernel[v][u]) * (IMAGE)->data.F32[yLocation][xLocation]; \ - /* The (0,0) kernel is subtracted from other kernels to preserve photometric scaling */ \ - if ((KERNELS)->spatialOrder > 0 && (INDEX) != (KERNELS)->subIndex) { \ - sub += FUNC(subKernel->kernel[v][u]) * (IMAGE)->data.F32[yLocation][xLocation]; \ - } \ - } \ - } \ - TARGET = FUNC(polyValue) * sum + FUNC(-1.0) * sub; \ - break; \ - } \ - default: \ - psAbort("Should never get here."); \ - } \ -} -#else // Generate the convolved pixel value static inline double convolvePixel(const pmSubtractionKernels *kernels, // Kernel basis functions @@ -383,7 +279,6 @@ int x, int y, // Pixel around which to convolve const psImage *image, // Image to convolve - const psImage *polyValues, // Spatial polynomial values - double (*weightFunc)(double value) // Function for weighting - ) + const psImage *polyValues // Spatial polynomial values + ) { int xOrder = kernels->xOrder->data.S32[index]; // Polynomial order in x @@ -396,12 +291,12 @@ int u = kernels->u->data.S32[index]; // Offset in x int v = kernels->v->data.S32[index]; // Offset in y - double value = weightFunc(polyValue) * image->data.F32[y + v][x + u]; // Value of convolution + double value = polyValue * image->data.F32[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 += weightFunc(-1.0) * image->data.F32[y][x]; + value -= image->data.F32[y][x]; } return value; } - /* Method for SPAM and FRIES is the same */ + // Method for SPAM and FRIES is the same case PM_SUBTRACTION_KERNEL_SPAM: case PM_SUBTRACTION_KERNEL_FRIES: { @@ -413,10 +308,10 @@ for (int v = vStart; v <= vStop; v++) { for (int u = uStart; u <= uStop; u++) { - sum += weightFunc(polyValue) * image->data.F32[y + v][x + u]; + sum += polyValue * image->data.F32[y + v][x + u]; } } if (kernels->spatialOrder > 0 && index != kernels->subIndex) { - /* The (0,0) element is subtracted from most kernels to preserve photometric scaling */ - sum += weightFunc(-1.0) * image->data.F32[y][x]; + // The (0,0) element is subtracted from most kernels to preserve photometric scaling + sum -= image->data.F32[y][x]; } return sum; @@ -431,17 +326,17 @@ for (int v = -size; v <= size; v++) { for (int u = -size; u <= size; u++) { - sum += weightFunc(kernel->kernel[v][u]) * image->data.F32[y + v][x + u]; + sum += kernel->kernel[v][u] * image->data.F32[y + v][x + u]; } } - value = weightFunc(polyValue) * sum; + value = polyValue * sum; } else { // Using delta function int u = kernels->u->data.S32[index]; // Offset in x int v = kernels->v->data.S32[index]; // Offset in y - value = weightFunc(polyValue) * image->data.F32[y + v][x + u]; // Value of convolution + value = polyValue * image->data.F32[y + v][x + u]; // Value of convolution } - /* The (0,0) delta function is subtracted from most kernels to preserve photometric scaling */ + // The (0,0) delta function is subtracted from most kernels to preserve photometric scaling if (kernels->spatialOrder > 0 && index != kernels->subIndex) { - value += weightFunc(-1.0) * image->data.F32[y][x]; + value -= image->data.F32[y][x]; } return value; @@ -455,12 +350,12 @@ for (int v = -size; v <= size; v++) { for (int u = -size; u <= size; u++) { - sum += weightFunc(kernel->kernel[v][u]) * image->data.F32[y + v][x + u]; + sum += kernel->kernel[v][u] * image->data.F32[y + v][x + u]; // The (0,0) kernel is subtracted from other kernels to preserve photometric scaling if (kernels->spatialOrder > 0 && index != kernels->subIndex) { - sub += weightFunc(subKernel->kernel[v][u]) * image->data.F32[y + v][x + u]; + sub += subKernel->kernel[v][u] * image->data.F32[y + v][x + u]; } } } - return weightFunc(polyValue) * sum + weightFunc(-1.0) * sub; + return polyValue * sum - sub; } default: @@ -469,7 +364,6 @@ return NAN; } -#endif - -// Weighting function for use with convolvePixel: no weighting applied, suitable for combining image pixels + +// Weighting function for use with solvedKernel: no weighting applied, suitable for combining image pixels static inline double imageWeighting(double value) { @@ -477,5 +371,5 @@ } -// Weighting function for use with convolvePixel: weighting suitable for combining variances +// Weighting function for use with solvedKernel: weighting suitable for combining variances static inline double varianceWeighting(double value) { @@ -674,15 +568,9 @@ 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 + reference->data.F32[y][x]); // Inverse square noise // Generate the convolutions for (int i = 0; i < numKernels; i++) { -#ifdef USE_FUNCTIONS_INSTEAD_OF_MACROS - convolutions->data.F64[i] = convolvePixel(kernels, i, x, y, reference, polyValues, - imageWeighting); -#else - CONVOLVE_PIXEL(convolutions->data.F64[i], kernels, i, x, y, reference, polyValues, - imageWeighting); -#endif + convolutions->data.F64[i] = convolvePixel(kernels, i, x, y, reference, polyValues); }