Index: /trunk/psModules/src/imcombine/pmSubtraction.c =================================================================== --- /trunk/psModules/src/imcombine/pmSubtraction.c (revision 14304) +++ /trunk/psModules/src/imcombine/pmSubtraction.c (revision 14305) @@ -4,6 +4,6 @@ * @author GLG, MHPCC * - * @version $Revision: 1.23 $ $Name: not supported by cvs2svn $ - * @date $Date: 2007-07-19 01:10:57 $ + * @version $Revision: 1.24 $ $Name: not supported by cvs2svn $ + * @date $Date: 2007-07-19 01:39:28 $ * * Copyright 2004-2007 Institute for Astronomy, University of Hawaii @@ -25,5 +25,5 @@ #include "pmSubtraction.h" -//#define USE_FUNCTIONS_INSTEAD_OF_MACROS // Can I pass the address of a static inline function? +#define USE_FUNCTIONS_INSTEAD_OF_MACROS // Can I pass the address of a static inline function? ////////////////////////////////////////////////////////////////////////////////////////////////////////////// @@ -209,4 +209,25 @@ 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++) { + kernel->kernel[v][u] += value * weightFunc(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 + kernel->kernel[v][u] += value; + } + // The (0,0) kernel is subtracted from other kernels to preserve photometric scaling + if (kernels->spatialOrder > 0 && i != kernels->subIndex) { + kernel->kernel[0][0] += subValue; + } + break; + } case PM_SUBTRACTION_KERNEL_ISIS: { psKernel *preCalc = kernels->preCalc->data[i];// Precalculated values @@ -355,4 +376,30 @@ return sum; } + 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 + double sub = 0.0; // Accumulated sum to subtract + 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]; + } + } + value = weightFunc(polyValue) * sum + weightFunc(-1.0) * sub; + } 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 + } + /* 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]; + } + return value; + } case PM_SUBTRACTION_KERNEL_ISIS: { psKernel *kernel = kernels->preCalc->data[index]; // The convolution kernel Index: /trunk/psModules/src/imcombine/pmSubtractionKernels.c =================================================================== --- /trunk/psModules/src/imcombine/pmSubtractionKernels.c (revision 14304) +++ /trunk/psModules/src/imcombine/pmSubtractionKernels.c (revision 14305) @@ -58,4 +58,5 @@ kernels->preCalc = NULL; kernels->size = size; + kernels->inner = 0; kernels->spatialOrder = spatialOrder; @@ -141,12 +142,4 @@ } -#if 0 - // Subtract a particular kernel in order to preserve photometric calibration across image - if (spatialOrder > 0 && index != kernels->subIndex) { - psKernel *subKernel = kernels->preCalc->data[kernels->subIndex]; // Kernel to subtract - (void)psBinaryOp(preCalc->image, preCalc->image, "-", subKernel->image); - } -#endif - // Iterate over spatial order. This loop creates the terms for // x^xOrder * y^yOrder such that (xOrder+yOrder) <= spatialOrder. @@ -390,4 +383,112 @@ } +// Grid United with Normal Kernel +pmSubtractionKernels *pmSubtractionKernelsGUNK(int size, int spatialOrder, const psVector *sigmas, + const psVector *orders, int inner) +{ + PS_ASSERT_INT_POSITIVE(size, NULL); + PS_ASSERT_INT_NONNEGATIVE(spatialOrder, NULL); + PS_ASSERT_VECTOR_NON_NULL(sigmas, NULL); + PS_ASSERT_VECTOR_TYPE(sigmas, PS_TYPE_F32, NULL); + PS_ASSERT_VECTOR_NON_NULL(orders, NULL); + PS_ASSERT_VECTOR_TYPE(orders, PS_TYPE_S32, NULL); + PS_ASSERT_VECTORS_SIZE_EQUAL(sigmas, orders, NULL); + PS_ASSERT_INT_NONNEGATIVE(inner, NULL); + PS_ASSERT_INT_LESS_THAN(inner, size, NULL); + + int numGaussians = sigmas->n; // Number of Gaussians + int numGaussianVars = 0; // Number of Gaussian variant functions in the kernel + for (int i = 0; i < numGaussians; i++) { + int gaussOrder = orders->data.S32[i]; // Polynomial order to apply to Gaussian + numGaussianVars += (gaussOrder + 1) * (gaussOrder + 2) / 2; + } + + int numInner = 2 * PS_SQR(inner) + 1; // Number of inner kernel elements + int numSpatial = (spatialOrder + 1) * (spatialOrder + 2) / 2; // Number of spatial variations of a kernel + + int num = (numGaussianVars + numInner) * numSpatial; // Total number of basis functions + + pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_GUNK, + size, spatialOrder); // The kernels + kernels->preCalc = psArrayAlloc(numGaussianVars * numSpatial); + kernels->inner = numGaussianVars * numSpatial; + + // Set the Gaussian kernel parameters + for (int i = 0, index = 0; i < numGaussians; i++) { + // Iterate over (u,v) order + for (int uOrder = 0; uOrder <= orders->data.S32[i]; uOrder++) { + for (int vOrder = 0; vOrder <= orders->data.S32[i] - uOrder; vOrder++) { + + + // Set the pre-calculated kernel + psKernel *preCalc = psKernelAlloc(-size, size, -size, size); + for (int v = -size; v <= size; v++) { + for (int u = -size; u <= size; u++) { + preCalc->kernel[v][u] = power(u, uOrder) * power(v, vOrder) * + expf(-0.5 * (PS_SQR(u) + PS_SQR(v)) / PS_SQR(sigmas->data.F32[i])); + } + } + + // Iterate over spatial order. This loop creates the terms for + // x^xOrder * y^yOrder such that (xOrder+yOrder) <= spatialOrder. + for (int xOrder = 0; xOrder <= spatialOrder; xOrder++) { + for (int yOrder = 0; yOrder <= spatialOrder - xOrder; yOrder++, index++) { + kernels->sigma->data.F32[index] = sigmas->data.F32[i]; + kernels->u->data.S32[index] = uOrder; + kernels->v->data.S32[index] = vOrder; + kernels->xOrder->data.S32[index] = xOrder; + kernels->yOrder->data.S32[index] = yOrder; + kernels->preCalc->data[index] = psMemIncrRefCounter(preCalc); + + psTrace("psModules.imcombine", 7, "Kernel %d: %f %d %d %d %d\n", index, + sigmas->data.F32[i], uOrder, vOrder, xOrder, yOrder); + } + } + + psFree(preCalc); // Drop reference + } + } + } + + + // Generate a grid set of kernels for each (u,v) + for (int v = - inner, index = kernels->inner; v <= inner; v++) { + for (int u = - inner; u <= inner; u++) { + // Iterate over spatial order. This loop creates the terms for + // x^xOrder * y^yOrder such that (xOrder+yOrder) <= spatialOrder. + for (int xOrder = 0; xOrder <= spatialOrder; xOrder++) { + for (int yOrder = 0; yOrder <= spatialOrder - xOrder; yOrder++, index++) { + kernels->u->data.S32[index] = u; + kernels->v->data.S32[index] = v; + kernels->xOrder->data.S32[index] = xOrder; + kernels->yOrder->data.S32[index] = yOrder; + + psTrace("psModules.imcombine", 7, "Kernel %d: %d %d %d %d\n", index, + u, v, xOrder, yOrder); + } + } + } + } + + kernels->subIndex = kernels->inner + numInner / 2; + assert(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); + + if (psTraceGetLevel("psModules.imcombine.kernel") >= 10) { + for (int i = 0; i < num; i++) { + psKernel *kernel = kernels->preCalc->data[i]; // Kernel of interest + psString kernelName = NULL; + psStringAppend(&kernelName, "kernel%03d.fits", i); + psFits *kernelFile = psFitsOpen(kernelName, "w"); + psFree(kernelName); + psFitsWriteImage(kernelFile, NULL, kernel->image, 0, NULL); + psFitsClose(kernelFile); + } + } + + return kernels; +} pmSubtractionKernels *pmSubtractionKernelsGenerate(pmSubtractionKernelsType type, int size, int spatialOrder, @@ -404,4 +505,6 @@ case PM_SUBTRACTION_KERNEL_FRIES: return pmSubtractionKernelsFRIES(size, spatialOrder, inner); + case PM_SUBTRACTION_KERNEL_GUNK: + return pmSubtractionKernelsGUNK(size, spatialOrder, sigmas, orders, inner); default: psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Unknown kernel type: %x", type); Index: /trunk/psModules/src/imcombine/pmSubtractionKernels.h =================================================================== --- /trunk/psModules/src/imcombine/pmSubtractionKernels.h (revision 14304) +++ /trunk/psModules/src/imcombine/pmSubtractionKernels.h (revision 14305) @@ -10,4 +10,5 @@ PM_SUBTRACTION_KERNEL_SPAM, ///< Summed Pixels for Advanced Matching --- summed delta functions PM_SUBTRACTION_KERNEL_FRIES, ///< Fibonacci Radius Increases Excellence of Subtraction + PM_SUBTRACTION_KERNEL_GUNK, ///< Grid United with Normal Kernel --- POIS and ISIS hybrid } pmSubtractionKernelsType; @@ -22,4 +23,5 @@ psArray *preCalc; ///< Array of images containing pre-calculated kernel (for ISIS) int size; ///< The half-size of the kernel + int inner; ///< The size of an inner region int spatialOrder; ///< The spatial order of the kernels } pmSubtractionKernels; @@ -60,4 +62,12 @@ ); +/// Generate GUNK kernels +pmSubtractionKernels *pmSubtractionKernelsGUNK(int size, ///< Half-size of the kernel + int spatialOrder, ///< Order of spatial variations + const psVector *sigmas, ///< Gaussian widths + const psVector *orders, ///< Polynomial order of gaussians + int inner ///< Inner radius containing grid of delta functions + ); + /// Generate a kernel of a specified type pmSubtractionKernels *pmSubtractionKernelsGenerate(pmSubtractionKernelsType type, ///< Kernel type