Index: trunk/psModules/src/imcombine/pmSubtractionKernels.c =================================================================== --- trunk/psModules/src/imcombine/pmSubtractionKernels.c (revision 13390) +++ trunk/psModules/src/imcombine/pmSubtractionKernels.c (revision 14106) @@ -15,4 +15,6 @@ psFree(kernels->v); psFree(kernels->sigma); + psFree(kernels->uStop); + psFree(kernels->vStop); psFree(kernels->xOrder); psFree(kernels->yOrder); @@ -48,7 +50,9 @@ kernels->u = psVectorAlloc(numBasisFunctions, PS_TYPE_S32); kernels->v = psVectorAlloc(numBasisFunctions, PS_TYPE_S32); + kernels->sigma = NULL; + kernels->uStop = NULL; + kernels->vStop = NULL; kernels->xOrder = psVectorAlloc(numBasisFunctions, PS_TYPE_S32); kernels->yOrder = psVectorAlloc(numBasisFunctions, PS_TYPE_S32); - kernels->sigma = NULL; kernels->subIndex = 0; kernels->preCalc = NULL; @@ -187,6 +191,208 @@ } +/// 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 + ) +{ + PS_ASSERT_INT_POSITIVE(size, NULL); + PS_ASSERT_INT_NONNEGATIVE(spatialOrder, NULL); + PS_ASSERT_INT_NONNEGATIVE(inner, NULL); + PS_ASSERT_INT_LARGER_THAN(size, inner, NULL); + PS_ASSERT_INT_POSITIVE(binning, NULL); + + // The outer region should be divisible by the "binning"; otherwise allocate remainder to the inner region + int numOuter = (size - inner) / binning; // Number of summed pixels in the outer region + int numInner = inner + (size - inner) % binning; // Number of pixels in the inner region + assert(numOuter * binning + numInner == size); + int numTotal = numOuter + numInner; // Total number of summed pixels + + psTrace("psModules.imcombine", 3, "Inner: %d Outer: %d\n", numInner, numOuter); + + int num = PS_SQR(2 * numTotal + 1) * + (spatialOrder + 1) * (spatialOrder + 2) / 2; // 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 + + kernels->uStop = psVectorAlloc(num, PS_TYPE_F32); + kernels->vStop = psVectorAlloc(num, PS_TYPE_F32); + + psVector *locations = psVectorAlloc(2 * numTotal + 1, PS_TYPE_S32); // Locations for each kernel element + psVector *widths = psVectorAlloc(2 * numTotal + 1, PS_TYPE_S32); // Widths for each kernel element + locations->data.S32[numTotal] = 0; + widths->data.S32[numTotal] = 0; + for (int i = 1; i <= numInner; i++) { + locations->data.S32[numTotal + i] = i; + widths->data.S32[numTotal + i] = 0; + locations->data.S32[numTotal - i] = - i; + widths->data.S32[numTotal - i] = 0; + } + for (int i = numInner + 1; i <= numTotal; i++) { + locations->data.S32[numTotal + i] = locations->data.S32[numTotal + i - 1] + + widths->data.S32[numTotal + i - 1] + 1; + widths->data.S32[numTotal + i] = binning - 1; + locations->data.S32[numTotal - i] = locations->data.S32[numTotal - i + 1] - + widths->data.S32[numTotal - i + 1] - binning; + widths->data.S32[numTotal - i] = binning - 1; + } + + if (psTraceGetLevel("psModules.imcombine") >= 10) { + for (int i = 0; i < 2 * numTotal + 1; i++) { + psTrace("psModules.imcombine", 10, "%d: %d -> %d\n", i, locations->data.S32[i], + locations->data.S32[i] + widths->data.S32[i]); + } + } + + // Set the kernel parameters + for (int i = - numTotal, index = 0; i <= numTotal; i++) { + int u = locations->data.S32[numTotal + i]; // Location of pixel + int uStop = u + widths->data.S32[numTotal + i]; // Width of pixel + + for (int j = - numTotal; j <= numTotal; j++) { + int v = locations->data.S32[numTotal + j]; // Location of pixel + int vStop = v + widths->data.S32[numTotal + j]; // Width of pixel + + // 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->uStop->data.S32[index] = uStop; + kernels->vStop->data.S32[index] = vStop; + kernels->xOrder->data.S32[index] = xOrder; + kernels->yOrder->data.S32[index] = yOrder; + + psTrace("psModules.imcombine", 7, "Kernel %d: %d %d %d %d %d %d\n", index, + u, uStop, v, vStop, xOrder, yOrder); + } + } + } + } + + kernels->subIndex = (num - (spatialOrder + 1) * (spatialOrder + 2) / 2) / 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 && + kernels->xOrder->data.S32[kernels->subIndex] == 0 && + kernels->yOrder->data.S32[kernels->subIndex] == 0); + + psFree(locations); + psFree(widths); + + return kernels; +} + + +/// 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 + ) +{ + PS_ASSERT_INT_POSITIVE(size, NULL); + PS_ASSERT_INT_NONNEGATIVE(spatialOrder, NULL); + PS_ASSERT_INT_NONNEGATIVE(inner, NULL); + PS_ASSERT_INT_LARGER_THAN(size, inner, NULL); + + int fibNum = 0; // Number of Fibonacci values + int fibLast = 1, fibTotal = 2; // Fibonacci sequence + while (fibTotal < size - inner) { + int temp = fibTotal; + fibTotal += fibLast; + fibLast = temp; + fibNum++; + } + + int numInner = inner; // Number of pixels in the inner region + int numOuter = fibNum; // Number of summed pixels in the outer region + int numTotal = numOuter + numInner; // Total number of summed pixels + + psTrace("psModules.imcombine", 3, "Inner: %d Outer: %d\n", numInner, numOuter); + + int num = PS_SQR(2 * numTotal + 1) * + (spatialOrder + 1) * (spatialOrder + 2) / 2; // 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 + kernels->uStop = psVectorAlloc(num, PS_TYPE_F32); + kernels->vStop = psVectorAlloc(num, PS_TYPE_F32); + + psVector *start = psVectorAlloc(2 * numTotal + 1, PS_TYPE_S32); + psVector *stop = psVectorAlloc(2 * numTotal + 1, PS_TYPE_S32); + start->data.S32[numTotal] = 0; + stop->data.S32[numTotal] = 0; + for (int i = 1; i <= numInner; i++) { + start->data.S32[numTotal + i] = i; + stop->data.S32[numTotal + i] = i; + start->data.S32[numTotal - i] = -i; + stop->data.S32[numTotal - i] = -i; + } + for (int i = numInner + 1, fibLast = 1, fib = 2, temp; i <= numTotal; + i++, fib = (temp = fib) + fibLast, fibLast = temp) { + start->data.S32[numTotal + i] = stop->data.S32[numTotal + i - 1] + 1; + stop->data.S32[numTotal + i] = PS_MIN(start->data.S32[numTotal + i] + fib - 1, size); + start->data.S32[numTotal - i] = - stop->data.S32[numTotal + i]; + stop->data.S32[numTotal - i] = - start->data.S32[numTotal + i]; + } + + if (psTraceGetLevel("psModules.imcombine") >= 10) { + for (int i = 0; i < 2 * numTotal + 1; i++) { + psTrace("psModules.imcombine", 10, "%d: %d -> %d\n", i, start->data.S32[i], stop->data.S32[i]); + } + } + + // Set the kernel parameters + for (int i = - numTotal, index = 0; i <= numTotal; i++) { + int u = start->data.S32[numTotal + i]; // Location of pixel + int uStop = stop->data.S32[numTotal + i]; // Width of pixel + for (int j = - numTotal; j <= numTotal; j++) { + int v = start->data.S32[numTotal + j]; // Location of pixel + int vStop = stop->data.S32[numTotal + j]; // Width of pixel + + // 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->uStop->data.S32[index] = uStop; + kernels->vStop->data.S32[index] = vStop; + kernels->xOrder->data.S32[index] = xOrder; + kernels->yOrder->data.S32[index] = yOrder; + + psTrace("psModules.imcombine", 7, "Kernel %d: %d %d %d %d %d %d\n", index, + u, uStop, v, vStop, xOrder, yOrder); + } + } + } + } + + kernels->subIndex = (num - (spatialOrder + 1) * (spatialOrder + 2) / 2) / 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 && + kernels->xOrder->data.S32[kernels->subIndex] == 0 && + kernels->yOrder->data.S32[kernels->subIndex] == 0); + + psFree(start); + psFree(stop); + + return kernels; +} + + pmSubtractionKernels *pmSubtractionKernelsGenerate(pmSubtractionKernelsType type, int size, int spatialOrder, - const psVector *sigmas, const psVector *orders) + const psVector *sigmas, const psVector *orders, int inner, + int binning) { switch (type) { @@ -195,4 +401,8 @@ case PM_SUBTRACTION_KERNEL_ISIS: return pmSubtractionKernelsISIS(size, spatialOrder, sigmas, orders); + case PM_SUBTRACTION_KERNEL_SPAM: + return pmSubtractionKernelsSPAM(size, spatialOrder, inner, binning); + case PM_SUBTRACTION_KERNEL_FRIES: + return pmSubtractionKernelsFRIES(size, spatialOrder, inner); default: psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Unknown kernel type: %x", type);