Changeset 14305

Timestamp:

Jul 18, 2007, 3:39:28 PM (19 years ago)

Author:

Paul Price

Message:

Adding GUNK (Grid United with Normal Kernels --- a POIS/ISIS hybrid) kernels.

Location:

trunk/psModules/src/imcombine

Files:

• pmSubtraction.c (modified) (4 diffs)
• pmSubtractionKernels.c (modified) (4 diffs)
• pmSubtractionKernels.h (modified) (3 diffs)

trunk/psModules/src/imcombine/pmSubtraction.c

@@ -4 +4 @@
  *  @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 +25 @@
 #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 +209 @@
               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 +376 @@
           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

trunk/psModules/src/imcombine/pmSubtractionKernels.c

@@ -58 +58 @@
     kernels->preCalc = NULL;
     kernels->size = size;
+    kernels->inner = 0;
     kernels->spatialOrder = spatialOrder;
 
@@ -141 +142 @@
                 }
 
-#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 +383 @@
 }
 
+// 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 +505 @@
       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);

trunk/psModules/src/imcombine/pmSubtractionKernels.h

@@ -10 +10 @@
     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 +23 @@
     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 +62 @@
     );
 
+/// 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