Changeset 14106 for trunk/psModules/src/imcombine/pmSubtractionKernels.c

Timestamp:

Jul 10, 2007, 1:54:26 PM (19 years ago)

Author:

Paul Price

Message:

Adding additional convolution kernels: SPAM (binned pixels around a central single-pixel core) and FRIES (binning changes as a (Fibonacci-like) function of radius).

File:

• trunk/psModules/src/imcombine/pmSubtractionKernels.c (modified) (4 diffs)

trunk/psModules/src/imcombine/pmSubtractionKernels.c

@@ -15 +15 @@
     psFree(kernels->v);
     psFree(kernels->sigma);
+    psFree(kernels->uStop);
+    psFree(kernels->vStop);
     psFree(kernels->xOrder);
     psFree(kernels->yOrder);
@@ -48 +50 @@
     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 +191 @@
 }
 
+/// 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 +401 @@
       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);