Index: trunk/psModules/src/imcombine/pmSubtraction.c =================================================================== --- trunk/psModules/src/imcombine/pmSubtraction.c (revision 13398) +++ trunk/psModules/src/imcombine/pmSubtraction.c (revision 13406) @@ -4,6 +4,6 @@ * @author GLG, MHPCC * - * @version $Revision: 1.10 $ $Name: not supported by cvs2svn $ - * @date $Date: 2007-05-16 21:13:51 $ + * @version $Revision: 1.11 $ $Name: not supported by cvs2svn $ + * @date $Date: 2007-05-17 04:13:19 $ * * Copyright 2004-2007 Institute for Astronomy, University of Hawaii @@ -53,4 +53,75 @@ } +// 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 if the weightFunc can't be "static inline", they still appear 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)->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; \ + } \ + case PM_SUBTRACTION_KERNEL_ISIS: { \ + psKernel *preCalc = (KERNELS)->preCalc->data[i]; /* Precalculated values */ \ + psKernel *subKernel = (KERNELS)->preCalc->data[(KERNELS)->subIndex]; /* Kernel to subtract */ \ + for (int v = -size; v <= size; v++) { \ + for (int u = -size; u <= size; u++) { \ + (TARGET)->kernel[v][u] += value * FUNC(preCalc->kernel[v][u]); \ + /* Subtract (0,0) kernel from other kernels to preserve photometric scaling */ \ + if ((KERNELS)->spatialOrder > 0 && i != (KERNELS)->subIndex) { \ + (TARGET)->kernel[v][u] += subValue * FUNC(subKernel->kernel[v][u]); \ + } \ + } \ + } \ + break; \ + } \ + default: \ + psAbort("Should never get here."); \ + } \ + } \ +} + +#if 0 // Generate an image of the solved kernel static inline psKernel *solvedKernel(psKernel *kernel, // Kernel, to return @@ -117,5 +188,65 @@ return kernel; } - +#endif + +// 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 if the weightFunc can't be "static inline", they still appear 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; \ + } \ + 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; v <= size; v++) { \ + for (int u = -size; u <= size; u++) { \ + sum += FUNC(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 += FUNC(subKernel->kernel[v][u]) * (IMAGE)->data.F32[(Y) + v][(X) + u]; \ + } \ + } \ + } \ + TARGET = FUNC(polyValue) * sum + FUNC(-1.0) * sub; \ + break; \ + } \ + default: \ + psAbort("Should never get here."); \ + } \ +} + +#if 0 // Generate the convolved pixel value static inline double convolvePixel(const pmSubtractionKernels *kernels, // Kernel basis functions @@ -165,4 +296,5 @@ return NAN; } +#endif // Weighting function for use with convolvePixel: no weighting applied, suitable for combining image pixels @@ -269,6 +401,11 @@ // Generate the convolutions for (int i = 0; i < numKernels; i++) { +#if 0 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 } @@ -450,7 +587,12 @@ for (int x = footprint; x <= footprint; x++) { for (int k = 0; k < numKernels; k++) { +#if 0 convolvedStamp->data.F32[y + footprint][x + footprint] += convolvePixel(kernels, k, xStamp + x, yStamp + y, refImage, polyValues, imageWeighting); +#else + CONVOLVE_PIXEL(convolvedStamp->data.F32[y + footprint][x + footprint], kernels, + k, xStamp + x, yStamp + y, refImage, polyValues, imageWeighting); +#endif } } @@ -532,5 +674,10 @@ // The appropriate kernel +#if 0 psKernel *kernel = solvedKernel(NULL, solution, kernels, polyValues, imageWeighting); +#else + psKernel *kernel = NULL; + SOLVED_KERNEL(kernel, solution, kernels, polyValues, imageWeighting); +#endif psFree(polyValues); @@ -630,4 +777,5 @@ 2.0 * (float)(i + size - numCols/2.0) / (float)numCols, 2.0 * (float)(j + size - numRows/2.0) / (float)numRows); +#if 0 kernelImage = solvedKernel(kernelImage, solution, kernels, polyValues, imageWeighting); if (inWeight) { @@ -635,4 +783,10 @@ varianceWeighting); } +#else + SOLVED_KERNEL(kernelImage, solution, kernels, polyValues, imageWeighting); + if (inWeight) { + SOLVED_KERNEL(kernelWeight, solution, kernels, polyValues, varianceWeighting); + } +#endif for (int y = j; y < j + fullSize; y++) {