Index: /trunk/psModules/src/imcombine/pmSubtraction.c =================================================================== --- /trunk/psModules/src/imcombine/pmSubtraction.c (revision 14359) +++ /trunk/psModules/src/imcombine/pmSubtraction.c (revision 14360) @@ -4,6 +4,6 @@ * @author GLG, MHPCC * - * @version $Revision: 1.31 $ $Name: not supported by cvs2svn $ - * @date $Date: 2007-07-20 21:18:50 $ + * @version $Revision: 1.32 $ $Name: not supported by cvs2svn $ + * @date $Date: 2007-07-21 04:30:57 $ * * Copyright 2004-2007 Institute for Astronomy, University of Hawaii @@ -170,4 +170,24 @@ 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 */ \ + 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."); \ @@ -274,4 +294,26 @@ break; } + case PM_SUBTRACTION_KERNEL_RINGS: { + if (i == kernels->subIndex) { + kernel->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 + + for (int j = 0; j < num; j++) { + int u = uCoords->data.S32[j], v = vCoords->data.S32[j]; // Kernel coordinates + kernel->kernel[v][u] += value * weightFunc(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."); @@ -367,4 +409,21 @@ } return polyValue * sum - sub; + } + case PM_SUBTRACTION_KERNEL_RINGS: { + if (index == kernels->subIndex) { + return image->data.F32[0][0]; + } + psArray *preCalc = kernels->preCalc->data[index]; // 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 + double sum = 0.0; // Accumulated sum from convolution + 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]; + } + // The (0,0) kernel is subtracted from other kernels to preserve photometric scaling + return polyValue * sum - image->data.F32[0][0]; } default: Index: /trunk/psModules/src/imcombine/pmSubtractionKernels.c =================================================================== --- /trunk/psModules/src/imcombine/pmSubtractionKernels.c (revision 14359) +++ /trunk/psModules/src/imcombine/pmSubtractionKernels.c (revision 14360) @@ -8,4 +8,6 @@ #include "pmSubtractionKernels.h" +#define RINGS_BUFFER 10 // Buffer size for RINGS data + // Free function for pmSubtractionKernels @@ -14,5 +16,5 @@ psFree(kernels->u); psFree(kernels->v); - psFree(kernels->sigma); + psFree(kernels->widths); psFree(kernels->uStop); psFree(kernels->vStop); @@ -50,5 +52,5 @@ kernels->u = psVectorAlloc(numBasisFunctions, PS_TYPE_S32); kernels->v = psVectorAlloc(numBasisFunctions, PS_TYPE_S32); - kernels->sigma = NULL; + kernels->widths = NULL; kernels->uStop = NULL; kernels->vStop = NULL; @@ -134,5 +136,5 @@ psFree(params); - kernels->sigma = psVectorAlloc(num, PS_TYPE_F32); + kernels->widths = psVectorAlloc(num, PS_TYPE_F32); kernels->preCalc = psArrayAlloc(num); @@ -156,5 +158,5 @@ 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->widths->data.F32[index] = sigmas->data.F32[i]; kernels->u->data.S32[index] = uOrder; kernels->v->data.S32[index] = vOrder; @@ -435,5 +437,5 @@ psFree(params); - kernels->sigma = psVectorAlloc(numGaussianVars * numSpatial, PS_TYPE_F32); + kernels->widths = psVectorAlloc(numGaussianVars * numSpatial, PS_TYPE_F32); kernels->preCalc = psArrayAlloc(numGaussianVars * numSpatial); kernels->inner = numGaussianVars * numSpatial; @@ -459,5 +461,5 @@ 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->widths->data.F32[index] = sigmas->data.F32[i]; kernels->u->data.S32[index] = uOrder; kernels->v->data.S32[index] = vOrder; @@ -505,7 +507,130 @@ } +// RINGS --- just what it says +pmSubtractionKernels *pmSubtractionKernelsRINGS(int size, int spatialOrder, int ringsOrder) +{ + PS_ASSERT_INT_POSITIVE(size, NULL); + PS_ASSERT_INT_NONNEGATIVE(spatialOrder, NULL); + PS_ASSERT_INT_NONNEGATIVE(ringsOrder, NULL); + + int numRings = size + 1; // Number of rings; 0 --> size, inclusive + int numPoly = (ringsOrder + 1) * (ringsOrder + 2) / 2; // Number of polynomial variants of each ring + int numSpatial = (spatialOrder + 1) * (spatialOrder + 2) / 2; // Number of spatial variations of a kernel + + int num = numRings * numPoly * numSpatial; // Total number of basis functions + + pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_GUNK, + size, spatialOrder); // The kernels + + psLogMsg("psModules.imcombine", PS_LOG_INFO, "RINGS kernel: %d,%d,%d --> %d elements", + size, ringsOrder, spatialOrder, num); + + kernels->widths = psVectorAlloc(num, PS_TYPE_S32); + kernels->preCalc = psArrayAlloc(num); + + // Set the Gaussian kernel parameters + for (int i = 0, index = 0; i < size; i++) { + // Iterate over (u,v) order + for (int uOrder = 0; uOrder <= (i == 0 ? 0 : ringsOrder); uOrder++) { + for (int vOrder = 0; vOrder <= (i == 0 ? 0 : ringsOrder) - uOrder; vOrder++) { + + psArray *data = psArrayAlloc(3); // Container for data + psVector *uCoords = data->data[0] = psVectorAllocEmpty(RINGS_BUFFER, PS_TYPE_S32); // u coords + psVector *vCoords = data->data[1] = psVectorAllocEmpty(RINGS_BUFFER, PS_TYPE_S32); // v coords + psVector *poly = data->data[2] = psVectorAllocEmpty(RINGS_BUFFER, PS_TYPE_F32); // Polynomial + if (i == 0) { + uCoords->data.S32[0] = vCoords->data.S32[0] = 0; + poly->data.F32[0] = 0; + uCoords->n = vCoords->n = poly->n = 1; + } else { + float radius = i; // Radius of ring + float lower2 = PS_SQR(radius - 0.5); // Lower limit of radius^2 + float upper2 = PS_SQR(radius + 0.5); // Upper limit of radius^2 + + int j = 0; // Index for data + for (int v = 1; v <= size; v++) { + int v2 = PS_SQR(v); // Square of v + float vPolyPlus = power(v, vOrder); // Value of (+v)^vOrder + float vPolyMinus = power(-v, vOrder); // Value of (-v)^vOrder + + // u = 0 + uCoords->data.S32[j] = 0; + vCoords->data.S32[j] = v; + poly->data.F32[j] = (uOrder == 0 ? vPolyPlus : 0.0); + j++; + + uCoords->data.S32[j] = 0; + vCoords->data.S32[j] = -v; + poly->data.F32[j] = (uOrder == 0 ? vPolyMinus : 0.0); + j++; + + psVectorExtend(uCoords, RINGS_BUFFER, 2); + psVectorExtend(vCoords, RINGS_BUFFER, 2); + psVectorExtend(poly, RINGS_BUFFER, 2); + + for (int u = 1; u <= v; u++) { + int u2 = PS_SQR(u); // Square of u + int distance2 = u2 + v2; // Distance from the centre + if (distance2 > lower2 && distance2 < upper2) { + float uPolyPlus = power(u, uOrder); // Value of (+u)^uOrder + float uPolyMinus = power(-u, uOrder); // Value of (-u)^uOrder + + uCoords->data.S32[j] = u; + vCoords->data.S32[j] = v; + poly->data.F32[j] = uPolyPlus * vPolyPlus; + j++; + + uCoords->data.S32[j] = u; + vCoords->data.S32[j] = -v; + poly->data.F32[j] = uPolyPlus * vPolyMinus; + j++; + + uCoords->data.S32[j] = -u; + vCoords->data.S32[j] = v; + poly->data.F32[j] = uPolyMinus * vPolyPlus; + j++; + + uCoords->data.S32[j] = -u; + vCoords->data.S32[j] = -v; + poly->data.F32[j] = uPolyMinus * vPolyMinus; + j++; + + psVectorExtend(uCoords, RINGS_BUFFER, 4); + psVectorExtend(vCoords, RINGS_BUFFER, 4); + psVectorExtend(poly, RINGS_BUFFER, 4); + } + } + } + } + + // 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->preCalc->data[index] = data; + kernels->widths->data.S32[index] = PS_SQR(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; + + psTrace("psModules.imcombine", 7, "Kernel %d: %d %d %d %d %d\n", index, + i, uOrder, vOrder, xOrder, yOrder); + } + } + } + } + } + + kernels->subIndex = 0; + assert(kernels->xOrder->data.S32[kernels->subIndex] == 0 && + kernels->yOrder->data.S32[kernels->subIndex] == 0); + + return kernels; +} + pmSubtractionKernels *pmSubtractionKernelsGenerate(pmSubtractionKernelsType type, int size, int spatialOrder, const psVector *sigmas, const psVector *orders, int inner, - int binning) + int binning, int ringsOrder) { switch (type) { @@ -520,4 +645,6 @@ case PM_SUBTRACTION_KERNEL_GUNK: return pmSubtractionKernelsGUNK(size, spatialOrder, sigmas, orders, inner); + case PM_SUBTRACTION_KERNEL_RINGS: + return pmSubtractionKernelsRINGS(size, ringsOrder, spatialOrder); 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 14359) +++ /trunk/psModules/src/imcombine/pmSubtractionKernels.h (revision 14360) @@ -11,4 +11,5 @@ PM_SUBTRACTION_KERNEL_FRIES, ///< Fibonacci Radius Increases Excellence of Subtraction PM_SUBTRACTION_KERNEL_GUNK, ///< Grid United with Normal Kernel --- POIS and ISIS hybrid + PM_SUBTRACTION_KERNEL_RINGS, ///< Rings Instead of the Normal Gaussian Subtraction } pmSubtractionKernelsType; @@ -17,5 +18,5 @@ pmSubtractionKernelsType type; ///< Type of kernels --- allowing the use of multiple kernels psVector *u, *v; ///< Offset (for POIS) or polynomial order (for ISIS) - psVector *sigma; ///< Gaussian widths (ISIS only) + psVector *widths; ///< Gaussian widths (ISIS) or ring radius (RINGS) psVector *uStop, *vStop; ///< Width of kernel element (SPAM,FRIES only) psVector *xOrder, *yOrder; ///< Spatial Polynomial order (for all) @@ -70,4 +71,11 @@ ); +/// Generate RINGS kernels +pmSubtractionKernels *pmSubtractionKernelsRINGS(int size, ///< Half-size of the kernel + int spatialOrder, ///< Order of spatial variations + int ringsOrder ///< Polynomial order + ); + + /// Generate a kernel of a specified type pmSubtractionKernels *pmSubtractionKernelsGenerate(pmSubtractionKernelsType type, ///< Kernel type @@ -77,5 +85,6 @@ const psVector *orders, ///< Polynomial order of gaussians int inner, ///< Inner radius to preserve unbinned - int binning ///< Kernel binning factor + int binning, ///< Kernel binning factor + int ringsOrder ///< Polynomial order for RINGS );