Index: trunk/psModules/src/imcombine/pmSubtraction.c =================================================================== --- trunk/psModules/src/imcombine/pmSubtraction.c (revision 14455) +++ trunk/psModules/src/imcombine/pmSubtraction.c (revision 14456) @@ -4,6 +4,6 @@ * @author GLG, MHPCC * - * @version $Revision: 1.34 $ $Name: not supported by cvs2svn $ - * @date $Date: 2007-08-09 20:25:52 $ + * @version $Revision: 1.35 $ $Name: not supported by cvs2svn $ + * @date $Date: 2007-08-09 21:20:35 $ * * Copyright 2004-2007 Institute for Astronomy, University of Hawaii @@ -25,6 +25,4 @@ #include "pmSubtraction.h" -#define USE_FUNCTIONS_INSTEAD_OF_MACROS // Can I pass the address of a static inline function? - ////////////////////////////////////////////////////////////////////////////////////////////////////////////// // Private (file-static) functions @@ -71,135 +69,13 @@ } -#ifndef USE_FUNCTIONS_INSTEAD_OF_MACROS -// Generate an image of the solved kernel -// Meant to replace the following function declaration: -// 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 -// double (*weightFunc)(double value) // Function for weighting -// ); -// because even though the weightFunc is "static inline", it still appears in the assembly. -// TARGET: Kernel, to return (psKernel*) -// SOLUTION: Solution to the least-squares problem (psVector*) -// KERNELS: Kernel basis functions (pmSubtractionKernels*) -// POLYVALUES: Spatial polynomial values (psImage*) -// FUNC: Function for weighting (double (*weightFunc)(double value)) -#define SOLVED_KERNEL(TARGET, SOLUTION, KERNELS, POLYVALUES, FUNC) { \ - int numKernels = (SOLUTION)->n - 1; /* Number of kernel basis functions */ \ - assert((KERNELS)->u->n == numKernels); \ - assert((KERNELS)->v->n == numKernels); \ - assert((KERNELS)->xOrder->n == numKernels); \ - assert((KERNELS)->yOrder->n == numKernels); \ - \ - /* 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)->u->data.S32[(KERNELS)->subIndex] == 0 && \ - (KERNELS)->v->data.S32[(KERNELS)->subIndex] == 0 && \ - (KERNELS)->xOrder->data.S32[(KERNELS)->subIndex] == 0 && \ - (KERNELS)->yOrder->data.S32[(KERNELS)->subIndex] == 0)); \ - \ - int size = (KERNELS)->size; /* Kernel half-size */ \ - if (!(TARGET)) { \ - (TARGET) = psKernelAlloc(-size, size, -size, size); \ - } \ - psImageInit((TARGET)->image, 0.0); \ - \ - for (int i = 0; i < numKernels; i++) { \ - int xOrder = (KERNELS)->xOrder->data.S32[i]; /* Polynomial order in x */ \ - int yOrder = (KERNELS)->yOrder->data.S32[i]; /* Polynomial order in y */ \ - float value = FUNC((POLYVALUES)->data.F64[yOrder][xOrder] * (SOLUTION)->data.F64[i]); \ - float subValue = FUNC(-(SOLUTION)->data.F64[i]); /* Value to subtract (add because of minus) */ \ - switch ((KERNELS)->type) { \ - case PM_SUBTRACTION_KERNEL_POIS: { \ - int u = (KERNELS)->u->data.S32[i]; /* Offset in x */ \ - int v = (KERNELS)->v->data.S32[i]; /* Offset in y */ \ - (TARGET)->kernel[v][u] += value; \ - if ((KERNELS)->spatialOrder > 0 && i != (KERNELS)->subIndex) { \ - /* The (0,0) element is subtracted from most kernels to preserve photometric scaling */ \ - (TARGET)->kernel[0][0] += subValue; \ - } \ - break; \ - } \ - /* SPAM and FRIES use the same method */ \ - case PM_SUBTRACTION_KERNEL_SPAM: \ - case PM_SUBTRACTION_KERNEL_FRIES: { \ - int uStart = (KERNELS)->u->data.S32[i]; \ - int uStop = (KERNELS)->uStop->data.S32[i]; \ - int vStart = (KERNELS)->v->data.S32[i]; \ - int vStop = (KERNELS)->vStop->data.S32[i]; \ - /* Normalising sum of kernel component to unity */ \ - value /= FUNC((uStop - uStart) * (vStop - vStart)); \ - for (int v = vStart; v <= vStop; v++) { \ - for (int u = uStart; u <= uStop; u++) { \ - (TARGET)->kernel[v][u] += value; \ - } \ - } \ - if ((KERNELS)->spatialOrder > 0 && i != (KERNELS)->subIndex) { \ - /* The (0,0) element is subtracted from most kernels to preserve photometric scaling */ \ - (TARGET)->kernel[0][0] += subValue; \ - } \ - break; \ - } \ - case PM_SUBTRACTION_KERNEL_GUNK: { \ - if (i < (KERNELS)->inner) { \ - /* Using pre-calculated function */ \ - psKernel *preCalc = (KERNELS)->preCalc->data[i]; /* Precalculated values */ \ - for (int v = -size; v <= size; v++) { \ - for (int u = -size; u <= size; u++) { \ - (TARGET)->kernel[v][u] += value * FUNC(preCalc->kernel[v][u]); \ - } \ - } \ - } else { \ - /* Using delta function */ \ - int u = (KERNELS)->u->data.S32[i]; /* Offset in x */ \ - int v = (KERNELS)->v->data.S32[i]; /* Offset in y */ \ - (TARGET)->kernel[v][u] += value; \ - } \ - /* The (0,0) kernel is subtracted from other kernels to preserve photometric scaling */ \ - if ((KERNELS)->spatialOrder > 0 && i != (KERNELS)->subIndex) { \ - (TARGET)->kernel[0][0] += subValue; \ - } \ - break; \ - } \ - case PM_SUBTRACTION_KERNEL_ISIS: { \ - psKernel *preCalc = (KERNELS)->preCalc->data[i]; /* Precalculated values */ \ - for (int v = -size; v <= size; v++) { \ - for (int u = -size; u <= size; u++) { \ - (TARGET)->kernel[v][u] += value * FUNC(preCalc->kernel[v][u]); \ - /* Photometric scaling is already built in to the precalculated kernel */ \ - } \ - } \ - break; \ - } \ - case PM_SUBTRACTION_KERNEL_RINGS: { \ - if (i == (KERNELS)->subIndex) { \ - (TARGET)->kernel[0][0] += value; \ - break; \ - } \ - psArray *preCalc = (KERNELS)->preCalc->data[i]; /* Precalculated data */ \ - psVector *uCoords = preCalc->data[0]; /* u coordinates */ \ - psVector *vCoords = preCalc->data[1]; /* v coordinates */ \ - psVector *poly = preCalc->data[2]; /* Polynomial values */ \ - int num = uCoords->n; /* Number of pixels */ \ - value /= weightFunc(num); /* Normalising sum of kernel component to unity */ \ - for (int j = 0; j < num; j++) { \ - int u = uCoords->data.S32[j], v = vCoords->data.S32[j]; /* Kernel coordinates */ \ - (TARGET)->kernel[v][u] += value * FUNC(poly->data.F32[j]); \ - } \ - /* The (0,0) kernel is subtracted from other kernels to preserve photometric scaling */ \ - if (kernels->spatialOrder > 0) { \ - kernel->kernel[0][0] += subValue; \ - } \ - break; \ - } \ - default: \ - psAbort("Should never get here."); \ - } \ - } \ -} -#else +// Apply the appropriate weighting to the kernel value +static inline float kernelWeighting(float value, // Value to apply weighting + bool variance // Do variance weighting? + ) +{ + return variance ? PS_SQR(value) : value; +} + + // Generate an image of the solved kernel static inline psKernel *solvedKernel(psKernel *kernel, // Kernel, to return @@ -207,6 +83,6 @@ const pmSubtractionKernels *kernels, // Kernel basis functions const psImage *polyValues, // Spatial polynomial values - double (*weightFunc)(double value) // Function for weighting - ) + bool varianceWeighting // Do variance weighting? + ) { int numKernels = kernels->num; // Number of kernel basis functions @@ -231,10 +107,12 @@ for (int i = 0; i < numKernels; i++) { - double value = 0.0; // Value of kernel coefficient + double sum = 0.0; // 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]; - } - } + sum += polyValues->data.F64[yOrder][xOrder] * solution->data.F64[index]; + } + } + float value = kernelWeighting(sum, varianceWeighting); // Value to be added to kernel + float subValue = kernelWeighting(- sum, varianceWeighting); // Value to be added for subtraction switch (kernels->type) { @@ -242,8 +120,8 @@ int u = kernels->u->data.S32[i]; // Offset in x int v = kernels->v->data.S32[i]; // Offset in y - kernel->kernel[v][u] += weightFunc(value); + 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] += weightFunc(- value); + kernel->kernel[0][0] += subValue; } break; @@ -258,12 +136,13 @@ // Normalising sum of kernel component to unity - float norm = 1.0 / weightFunc((uStop - uStart) * (vStop - vStart)); + float norm = kernelWeighting(1.0 / (float)((uStop - uStart) * (vStop - vStart)), + varianceWeighting); for (int v = vStart; v <= vStop; v++) { for (int u = uStart; u <= uStop; u++) { - kernel->kernel[v][u] += norm * weightFunc(value); + 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] += weightFunc(- value); + kernel->kernel[0][0] += subValue; } } @@ -278,5 +157,6 @@ for (int v = -size; v <= size; v++) { for (int u = -size; u <= size; u++) { - kernel->kernel[v][u] += weightFunc(preCalc->kernel[v][u] * value); + kernel->kernel[v][u] += kernelWeighting(preCalc->kernel[v][u], + varianceWeighting) * value; // Photometric scaling is built into the preCalc kernel --- no subtraction! } @@ -287,8 +167,8 @@ int u = kernels->u->data.S32[i]; // Offset in x int v = kernels->v->data.S32[i]; // Offset in y - kernel->kernel[v][u] += weightFunc(value); + kernel->kernel[v][u] += value; if (subtract) { // The (0,0) element is subtracted from most kernels to preserve photometric scaling - kernel->kernel[0][0] += weightFunc(- value); + kernel->kernel[0][0] += subValue; } } @@ -300,5 +180,6 @@ for (int v = -size; v <= size; v++) { for (int u = -size; u <= size; u++) { - kernel->kernel[v][u] += weightFunc(preCalc->kernel[v][u] * value); + kernel->kernel[v][u] += kernelWeighting(preCalc->kernel[v][u], + varianceWeighting) * value; // Photometric scaling is built into the preCalc kernel --- no subtraction! } @@ -308,5 +189,5 @@ case PM_SUBTRACTION_KERNEL_RINGS: { if (i == kernels->subIndex) { - kernel->kernel[0][0] += weightFunc(value); + kernel->kernel[0][0] += value; break; } @@ -319,7 +200,7 @@ for (int j = 0; j < num; j++) { int u = uCoords->data.S32[j], v = vCoords->data.S32[j]; // Kernel coordinates - kernel->kernel[v][u] += weightFunc(value * poly->data.F32[j]); + kernel->kernel[v][u] += kernelWeighting(poly->data.F32[j], varianceWeighting) * value; } - kernel->kernel[0][0] += weightFunc(- value * num); + kernel->kernel[0][0] += kernelWeighting(num, varianceWeighting) * subValue; break; } @@ -331,5 +212,4 @@ return kernel; } -#endif // Generate the convolved pixel value @@ -429,17 +309,4 @@ return NAN; } - -// Weighting function for use with solvedKernel: no weighting applied, suitable for combining image pixels -static inline double imageWeighting(double value) -{ - return value; -} - -// Weighting function for use with solvedKernel: weighting suitable for combining variances -static inline double varianceWeighting(double value) -{ - return PS_SQR(value); -} - // Mark a pixel as blank in the image, mask and weight @@ -830,9 +697,5 @@ psImage *polyValues = spatialPolyValues(kernels->spatialOrder, xNorm, yNorm); // Polynomial terms -#ifdef USE_FUNCTIONS_INSTEAD_OF_MACROS - kernelImage = solvedKernel(kernelImage, solution, kernels, polyValues, imageWeighting); -#else - SOLVED_KERNEL(kernelImage, solution, kernels, polyValues, imageWeighting); -#endif + kernelImage = solvedKernel(kernelImage, solution, kernels, polyValues, false); for (int y = yStamp - footprint, j = 0; y <= yStamp + footprint; y++, j++) { for (int x = xStamp - footprint, i = 0; x <= xStamp + footprint; x++, i++) { @@ -923,10 +786,5 @@ // The appropriate kernel -#ifdef USE_FUNCTIONS_INSTEAD_OF_MACROS - psKernel *kernel = solvedKernel(NULL, solution, kernels, polyValues, imageWeighting); -#else - psKernel *kernel = NULL; - SOLVED_KERNEL(kernel, solution, kernels, polyValues, imageWeighting); -#endif + psKernel *kernel = solvedKernel(NULL, solution, kernels, polyValues, false); psFree(polyValues); @@ -1077,16 +935,8 @@ 2.0 * (float)(i + size + 1 - numCols/2.0) / (float)numCols, 2.0 * (float)(j + size + 1 - numRows/2.0) / (float)numRows); -#ifdef USE_FUNCTIONS_INSTEAD_OF_MACROS - kernelImage = solvedKernel(kernelImage, solution, kernels, polyValues, imageWeighting); + kernelImage = solvedKernel(kernelImage, solution, kernels, polyValues, false); if (inWeight) { - kernelWeight = solvedKernel(kernelWeight, solution, kernels, polyValues, - varianceWeighting); - } -#else - SOLVED_KERNEL(kernelImage, solution, kernels, polyValues, imageWeighting); - if (inWeight) { - SOLVED_KERNEL(kernelWeight, solution, kernels, polyValues, varianceWeighting); - } -#endif + kernelWeight = solvedKernel(kernelWeight, solution, kernels, polyValues, true); + } for (int y = j; y < PS_MIN(j + fullSize, yMax); y++) {