Index: trunk/psModules/src/imcombine/pmSubtractionKernels.c =================================================================== --- trunk/psModules/src/imcombine/pmSubtractionKernels.c (revision 15432) +++ trunk/psModules/src/imcombine/pmSubtractionKernels.c (revision 15756) @@ -4,7 +4,9 @@ #include +#include #include #include +#include "pmSubtraction.h" #include "pmSubtractionKernels.h" @@ -22,7 +24,9 @@ psFree(kernels->vStop); psFree(kernels->preCalc); -} - -// Raise a number to an integer power + psFree(kernels->solution1); + psFree(kernels->solution2); +} + +// Raise an integer to an integer power static inline long power(int value, // Value int exp // Exponent @@ -45,16 +49,9 @@ bool p_pmSubtractionKernelsAddGrid(pmSubtractionKernels *kernels, int start, int size) { - PS_ASSERT_PTR_NON_NULL(kernels, false); - PS_ASSERT_VECTOR_NON_NULL(kernels->u, false); - PS_ASSERT_VECTOR_NON_NULL(kernels->v, false); - PS_ASSERT_VECTOR_NON_NULL(kernels->widths, false); - PS_ASSERT_VECTOR_TYPE(kernels->u, PS_TYPE_S32, false); - PS_ASSERT_VECTOR_TYPE(kernels->v, PS_TYPE_S32, false); - PS_ASSERT_VECTOR_TYPE(kernels->widths, PS_TYPE_F32, false); - PS_ASSERT_ARRAY_NON_NULL(kernels->preCalc, false); + PM_ASSERT_SUBTRACTION_KERNELS_NON_NULL(kernels, false); PS_ASSERT_INT_NONNEGATIVE(start, false); PS_ASSERT_INT_NONNEGATIVE(size, false); - int numNew = PS_SQR(2 * size + 1); // Number of new kernel parameters to add + int numNew = PS_SQR(2 * size + 1) - 1; // Number of new kernel parameters to add // Ensure the sizes match @@ -68,4 +65,9 @@ for (int v = - size, index = start; v <= size; v++) { for (int u = - size; u <= size; u++, index++) { + if (v == 0 && u == 0) { + // Skip normalisation component: added explicitly + index--; + continue; + } kernels->widths->data.F32[index] = NAN; kernels->u->data.S32[index] = u; @@ -81,5 +83,6 @@ pmSubtractionKernels *p_pmSubtractionKernelsRawISIS(int size, int spatialOrder, - const psVector *fwhms, const psVector *orders) + const psVector *fwhms, const psVector *orders, + pmSubtractionMode mode) { PS_ASSERT_VECTOR_NON_NULL(fwhms, NULL); @@ -97,10 +100,10 @@ for (int i = 0; i < numGaussians; i++) { int gaussOrder = orders->data.S32[i]; // Polynomial order to apply to Gaussian - psStringAppend(¶ms, "(%.2f,%d)", fwhms->data.F32[i], orders->data.S32[i]); + psStringAppend(¶ms, "(%.1f,%d)", fwhms->data.F32[i], orders->data.S32[i]); num += (gaussOrder + 1) * (gaussOrder + 2) / 2; } pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_ISIS, - size, spatialOrder); // The kernels + size, spatialOrder, mode); // The kernels psStringAppend(&kernels->description, "ISIS(%d,%s,%d)", size, params, spatialOrder); @@ -149,6 +152,4 @@ } - kernels->subIndex = -1; - return kernels; } @@ -159,5 +160,5 @@ pmSubtractionKernels *pmSubtractionKernelsAlloc(int numBasisFunctions, pmSubtractionKernelsType type, - int size, int spatialOrder) + int size, int spatialOrder, pmSubtractionMode mode) { pmSubtractionKernels *kernels = psAlloc(sizeof(pmSubtractionKernels)); // Kernels, to return @@ -173,22 +174,24 @@ kernels->uStop = NULL; kernels->vStop = NULL; - kernels->subIndex = 0; kernels->size = size; kernels->inner = 0; kernels->spatialOrder = spatialOrder; kernels->bgOrder = 0; - - return kernels; -} - -pmSubtractionKernels *pmSubtractionKernelsPOIS(int size, int spatialOrder) + kernels->mode = mode; + kernels->solution1 = NULL; + kernels->solution2 = NULL; + + return kernels; +} + +pmSubtractionKernels *pmSubtractionKernelsPOIS(int size, int spatialOrder, pmSubtractionMode mode) { PS_ASSERT_INT_POSITIVE(size, NULL); PS_ASSERT_INT_NONNEGATIVE(spatialOrder, NULL); - int num = PS_SQR(2 * size + 1); // Number of basis functions + int num = PS_SQR(2 * size + 1) - 1; // Number of basis functions pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_POIS, - size, spatialOrder); // The kernels + size, spatialOrder, mode); // The kernels psStringAppend(&kernels->description, "POIS(%d,%d)", size, spatialOrder); psLogMsg("psModules.imcombine", PS_LOG_INFO, "POIS kernel: %d,%d --> %d elements", @@ -199,8 +202,4 @@ } - kernels->subIndex = num / 2; - assert(kernels->u->data.S32[kernels->subIndex] == 0 && - kernels->v->data.S32[kernels->subIndex] == 0); - return kernels; } @@ -208,23 +207,22 @@ pmSubtractionKernels *pmSubtractionKernelsISIS(int size, int spatialOrder, - const psVector *fwhms, const psVector *orders) + const psVector *fwhms, const psVector *orders, + pmSubtractionMode mode) { pmSubtractionKernels *kernels = p_pmSubtractionKernelsRawISIS(size, spatialOrder, - fwhms, orders); // Kernels + fwhms, orders, mode); // Kernels if (!kernels) { return NULL; } - // Subtract a reference kernel from all the others to maintain flux scaling + // Subtract normalisation from all the others to maintain flux scaling if (spatialOrder > 0) { - int subIndex = 0; // Index of kernel to subtract from others - psKernel *subtract = kernels->preCalc->data[subIndex]; // Kernel to subtract from others int numGaussians = fwhms->n; // Number of Gaussians for (int i = 0, index = 0; i < numGaussians; i++) { for (int uOrder = 0; uOrder <= orders->data.S32[i]; uOrder++) { for (int vOrder = 0; vOrder <= orders->data.S32[i] - uOrder; vOrder++, index++) { - if (uOrder % 2 == 0 && vOrder % 2 == 0 && index != subIndex) { + if (uOrder % 2 == 0 && vOrder % 2 == 0) { psKernel *kernel = kernels->preCalc->data[index]; // Kernel to subtract - psBinaryOp(kernel->image, kernel->image, "-", subtract->image); + kernel->kernel[0][0] -= 1.0; } } @@ -242,16 +240,19 @@ psFitsWriteImage(kernelFile, NULL, kernel->image, 0, NULL); psFitsClose(kernelFile); - } - } - - return kernels; -} - -/// Generate SPAM kernels -pmSubtractionKernels *pmSubtractionKernelsSPAM(int size, ///< Half-size of the kernel - int spatialOrder, ///< Order of spatial variations - int inner, ///< Inner radius to preserve unbinned - int binning ///< Kernel binning factor - ) + double sum = 0.0; + for (int y = 0; y < kernel->image->numRows; y++) { + for (int x = 0; x < kernel->image->numCols; x++) { + sum += kernel->image->data.F32[y][x]; + } + } + psTrace("psModules.imcombine.kernel", 10, "Kernel %d sum: %lf\n", i, sum); + } + } + + return kernels; +} + +pmSubtractionKernels *pmSubtractionKernelsSPAM(int size, int spatialOrder, int inner, int binning, + pmSubtractionMode mode) { PS_ASSERT_INT_POSITIVE(size, NULL); @@ -269,10 +270,10 @@ psTrace("psModules.imcombine", 3, "Inner: %d Outer: %d\n", numInner, numOuter); - int num = PS_SQR(2 * numTotal + 1); // Number of basis functions + int num = PS_SQR(2 * numTotal + 1) - 1; // Number of basis functions psTrace("psModules.imcombine", 3, "Number of basis functions: %d\n", num); pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_SPAM, - size, spatialOrder); // The kernels + size, spatialOrder, mode); // The kernels psStringAppend(&kernels->description, "SPAM(%d,%d,%d,%d)", size, inner, binning, spatialOrder); @@ -280,6 +281,6 @@ size, inner, binning, spatialOrder, num); - kernels->uStop = psVectorAlloc(num, PS_TYPE_F32); - kernels->vStop = psVectorAlloc(num, PS_TYPE_F32); + kernels->uStop = psVectorAlloc(num, PS_TYPE_S32); + kernels->vStop = psVectorAlloc(num, PS_TYPE_S32); psVector *locations = psVectorAlloc(2 * numTotal + 1, PS_TYPE_S32); // Locations for each kernel element @@ -314,4 +315,9 @@ for (int j = - numTotal; j <= numTotal; j++, index++) { + if (i == 0 && j == 0) { + // Skip normalisation component: added explicitly + index--; + continue; + } int v = locations->data.S32[numTotal + j]; // Location of pixel int vStop = v + widths->data.S32[numTotal + j]; // Width of pixel @@ -327,10 +333,4 @@ } - kernels->subIndex = num / 2; - assert(kernels->u->data.S32[kernels->subIndex] == 0 && - kernels->v->data.S32[kernels->subIndex] == 0 && - kernels->uStop->data.S32[kernels->subIndex] == 0 && - kernels->vStop->data.S32[kernels->subIndex] == 0); - psFree(locations); psFree(widths); @@ -340,9 +340,5 @@ -/// Generate FRIES kernels -pmSubtractionKernels *pmSubtractionKernelsFRIES(int size, ///< Half-size of the kernel - int spatialOrder, ///< Order of spatial variations - int inner ///< Inner radius to preserve unbinned - ) +pmSubtractionKernels *pmSubtractionKernelsFRIES(int size, int spatialOrder, int inner, pmSubtractionMode mode) { PS_ASSERT_INT_POSITIVE(size, NULL); @@ -366,10 +362,10 @@ psTrace("psModules.imcombine", 3, "Inner: %d Outer: %d\n", numInner, numOuter); - int num = PS_SQR(2 * numTotal + 1); // Number of basis functions + int num = PS_SQR(2 * numTotal + 1) - 1; // Number of basis functions psTrace("psModules.imcombine", 3, "Number of basis functions: %d\n", num); pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_FRIES, - size, spatialOrder); // The kernels + size, spatialOrder, mode); // The kernels psStringAppend(&kernels->description, "FRIES(%d,%d,%d)", size, inner, spatialOrder); @@ -377,6 +373,6 @@ size, inner, spatialOrder, num); - kernels->uStop = psVectorAlloc(num, PS_TYPE_F32); - kernels->vStop = psVectorAlloc(num, PS_TYPE_F32); + kernels->uStop = psVectorAlloc(num, PS_TYPE_S32); + kernels->vStop = psVectorAlloc(num, PS_TYPE_S32); psVector *start = psVectorAlloc(2 * numTotal + 1, PS_TYPE_S32); @@ -409,4 +405,9 @@ int uStop = stop->data.S32[numTotal + i]; // Width of pixel for (int j = - numTotal; j <= numTotal; j++, index++) { + if (i == 0 && j == 0) { + // Skip normalisation component: added explicitly + index--; + continue; + } int v = start->data.S32[numTotal + j]; // Location of pixel int vStop = stop->data.S32[numTotal + j]; // Width of pixel @@ -422,10 +423,4 @@ } - kernels->subIndex = num / 2; - assert(kernels->u->data.S32[kernels->subIndex] == 0 && - kernels->v->data.S32[kernels->subIndex] == 0 && - kernels->uStop->data.S32[kernels->subIndex] == 0 && - kernels->vStop->data.S32[kernels->subIndex] == 0); - psFree(start); psFree(stop); @@ -436,5 +431,5 @@ // Grid United with Normal Kernel pmSubtractionKernels *pmSubtractionKernelsGUNK(int size, int spatialOrder, const psVector *fwhms, - const psVector *orders, int inner) + const psVector *orders, int inner, pmSubtractionMode mode) { PS_ASSERT_INT_POSITIVE(size, NULL); @@ -449,5 +444,5 @@ pmSubtractionKernels *kernels = p_pmSubtractionKernelsRawISIS(size, spatialOrder, - fwhms, orders); // Kernels + fwhms, orders, mode); // Kernels psStringPrepend(&kernels->description, "GUNK="); psStringAppend(&kernels->description, "+POIS(%d,%d)", inner, spatialOrder); @@ -469,5 +464,4 @@ int numISIS = kernels->num; // Number of ISIS kernels - int numInner = PS_SQR(2 * inner + 1); // Number of inner kernel elements if (!p_pmSubtractionKernelsAddGrid(kernels, numISIS, inner)) { @@ -475,13 +469,10 @@ } - kernels->subIndex = numInner/2 + numISIS; - assert(kernels->u->data.S32[kernels->subIndex] == 0 && - kernels->v->data.S32[kernels->subIndex] == 0); - return kernels; } // RINGS --- just what it says -pmSubtractionKernels *pmSubtractionKernelsRINGS(int size, int spatialOrder, int inner, int ringsOrder) +pmSubtractionKernels *pmSubtractionKernelsRINGS(int size, int spatialOrder, int inner, int ringsOrder, + pmSubtractionMode mode) { PS_ASSERT_INT_POSITIVE(size, NULL); @@ -505,14 +496,14 @@ } - int numInner = inner; // Number of pixels in the inner region + int numInner = inner - 1; // Number of pixels in the inner region int numOuter = fibNum; // Number of summed pixels in the outer region int numRings = numOuter + numInner; // Number of rings (not including the central pixel) - int numPoly = (ringsOrder + 1) * (ringsOrder + 2) / 2; // Number of polynomial variants of each ring - - int num = numRings * numPoly + 1; // Total number of basis functions + int numPoly = PM_SUBTRACTION_POLYTERMS(ringsOrder); // Number of polynomial variants of each ring + + int num = numRings * numPoly; // Total number of basis functions pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_RINGS, - size, spatialOrder); // The kernels + size, spatialOrder, mode); // The kernels psStringAppend(&kernels->description, "RINGS(%d,%d,%d,%d)", size, inner, ringsOrder, spatialOrder); @@ -522,13 +513,9 @@ // Set the Gaussian kernel parameters int fibIndex = 1, fibIndexMinus1 = 0; // Fibonnacci parameters - int radiusLast = 0; // Last radius - for (int i = 0, index = 0; i < numRings + 1; i++) { - + int radiusLast = 1; // Last radius + for (int i = 1, index = 0; i < numRings + 1; i++) { float lower2; // Lower limit of radius^2 float upper2; // Upper limit of radius^2 - if (i == 0) { - lower2 = 0; - upper2 = PS_SQR(0.5); - } else if (i <= numInner) { + if (i <= inner) { // A ring every pixel width float radius = i; @@ -572,5 +559,4 @@ for (int v = -size; v <= size; v++) { int v2 = PS_SQR(v); // Square of v -// float vPoly = power(v, vOrder); // Value of v^vOrder float vPoly = powf(v/(float)size, vOrder); // Value of v^vOrder @@ -579,5 +565,4 @@ int distance2 = u2 + v2; // Distance from the centre if (distance2 > lower2 && distance2 < upper2) { -// float uPoly = power(u, uOrder); // Value of u^uOrder float uPoly = powf(u/(float)size, uOrder); // Value of u^uOrder @@ -627,6 +612,4 @@ } - kernels->subIndex = 0; - return kernels; } @@ -634,19 +617,19 @@ pmSubtractionKernels *pmSubtractionKernelsGenerate(pmSubtractionKernelsType type, int size, int spatialOrder, const psVector *fwhms, const psVector *orders, int inner, - int binning, int ringsOrder) + int binning, int ringsOrder, pmSubtractionMode mode) { switch (type) { case PM_SUBTRACTION_KERNEL_POIS: - return pmSubtractionKernelsPOIS(size, spatialOrder); + return pmSubtractionKernelsPOIS(size, spatialOrder, mode); case PM_SUBTRACTION_KERNEL_ISIS: - return pmSubtractionKernelsISIS(size, spatialOrder, fwhms, orders); + return pmSubtractionKernelsISIS(size, spatialOrder, fwhms, orders, mode); case PM_SUBTRACTION_KERNEL_SPAM: - return pmSubtractionKernelsSPAM(size, spatialOrder, inner, binning); + return pmSubtractionKernelsSPAM(size, spatialOrder, inner, binning, mode); case PM_SUBTRACTION_KERNEL_FRIES: - return pmSubtractionKernelsFRIES(size, spatialOrder, inner); + return pmSubtractionKernelsFRIES(size, spatialOrder, inner, mode); case PM_SUBTRACTION_KERNEL_GUNK: - return pmSubtractionKernelsGUNK(size, spatialOrder, fwhms, orders, inner); + return pmSubtractionKernelsGUNK(size, spatialOrder, fwhms, orders, inner, mode); case PM_SUBTRACTION_KERNEL_RINGS: - return pmSubtractionKernelsRINGS(size, spatialOrder, inner, ringsOrder); + return pmSubtractionKernelsRINGS(size, spatialOrder, inner, ringsOrder, mode); default: psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Unknown kernel type: %x", type);