Changeset 26473

Timestamp:

Dec 22, 2009, 9:55:18 AM (16 years ago)

Author:

eugene

Message:

adding hermitian kernels

Location:

branches/eam_branches/20091201/psModules/src/imcombine

Files:

• pmSubtractionKernels.c (modified) (7 diffs)
• pmSubtractionKernels.h (modified) (4 diffs)

branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionKernels.c

@@ -10 +10 @@
 #include "pmSubtraction.h"
 #include "pmSubtractionKernels.h"
+#include "pmSubtractionHermitian.h"
 
 #define RINGS_BUFFER 10                 // Buffer size for RINGS data
-
 
 // Free function for pmSubtractionKernels
@@ -57 +57 @@
     for (int i = 0, x = -size; x <= size; i++, x++) {
         kernel->data.F32[i] = norm * power(x, order) * expf(expNorm * PS_SQR(x));
+    }
+
+    return kernel;
+}
+
+// Generate 1D convolution kernel for HERM (normalized for 2D)
+static psVector *subtractionKernelHERM(float sigma, // Gaussian width
+                                       int order, // Polynomial order
+                                       int size // Kernel half-size
+    )
+{
+    int fullSize = 2 * size + 1;        // Full size of kernel
+    psVector *kernel = psVectorAlloc(fullSize, PS_TYPE_F32); // Kernel to return
+
+    // for now, we are only allowing equal orders and sigmas in X and Y
+    float nf = exp(lgamma(order + 1));
+    float norm = 1.0 / sqrt(nf*sigma*sqrt(M_2_PI));
+
+    for (int i = 0, x = -size; x <= size; i++, x++) {
+        float xf = x / sigma;
+        float z = -0.25*xf*xf;
+        kernel->data.F32[i] = norm * p_pmSubtractionHermitianPolynomial(xf, order) * exp(z);
     }
 
@@ -189 +211 @@
 
 #if 0
+                double sum = 0.0;   // Sum of kernel component
+                for (int v = -size; v <= size; v++) {
+                    for (int u = -size; u <= size; u++) {
+                        sum += kernel->kernel[v][u];
+                    }
+                }
+                fprintf(stderr, "%d sum: %lf\n", index, sum);
+#endif
+
+                kernels->widths->data.F32[index] = fwhms->data.F32[i];
+                kernels->u->data.S32[index] = uOrder;
+                kernels->v->data.S32[index] = vOrder;
+                if (kernels->preCalc->data[index]) {
+                    psFree(kernels->preCalc->data[index]);
+                }
+                kernels->preCalc->data[index] = preCalc;
+                kernels->penalties->data.F32[index] = kernels->penalty * fabsf(moment);
+
+                psTrace("psModules.imcombine", 7, "Kernel %d: %f %d %d %f\n", index,
+                        fwhms->data.F32[i], uOrder, vOrder, fabsf(moment));
+            }
+        }
+    }
+
+    return kernels;
+}
+
+pmSubtractionKernels *pmSubtractionKernelsHERM(int size, int spatialOrder,
+                                               const psVector *fwhmsIN, const psVector *ordersIN,
+                                               float penalty, pmSubtractionMode mode)
+{
+    PS_ASSERT_VECTOR_NON_NULL(fwhmsIN, NULL);
+    PS_ASSERT_VECTOR_TYPE(fwhmsIN, PS_TYPE_F32, NULL);
+    PS_ASSERT_VECTOR_NON_NULL(ordersIN, NULL);
+    PS_ASSERT_VECTOR_TYPE(ordersIN, PS_TYPE_S32, NULL);
+    PS_ASSERT_VECTORS_SIZE_EQUAL(fwhmsIN, ordersIN, NULL);
+    PS_ASSERT_INT_POSITIVE(size, NULL);
+    PS_ASSERT_INT_NONNEGATIVE(spatialOrder, NULL);
+
+    // check the requested fwhm values: any values <= 0.0 should be dropped
+    psVector *fwhms  = psVectorAllocEmpty (fwhmsIN->n, PS_TYPE_F32);
+    psVector *orders = psVectorAllocEmpty (ordersIN->n, PS_TYPE_S32);
+    for (int i = 0; i < fwhmsIN->n; i++) {
+        if (fwhmsIN->data.F32[i] <= FLT_EPSILON) continue;
+        psVectorAppend(fwhms, fwhmsIN->data.F32[i]);
+        psVectorAppend(orders, ordersIN->data.S32[i]);
+    }
+
+    int numGaussians = fwhms->n;       // Number of Gaussians
+
+    int num = 0;                        // Number of basis functions
+    psString params = NULL;             // List of parameters
+    for (int i = 0; i < numGaussians; i++) {
+        int gaussOrder = orders->data.S32[i]; // Polynomial order to apply to Gaussian
+        psStringAppend(&params, "(%.1f,%d)", fwhms->data.F32[i], orders->data.S32[i]);
+        num += (gaussOrder + 1) * (gaussOrder + 2) / 2;
+    }
+
+    pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_HERM, size, spatialOrder, penalty, mode); // The kernels
+    psStringAppend(&kernels->description, "HERM(%d,%s,%d,%.2e)", size, params, spatialOrder, penalty);
+
+    psLogMsg("psModules.imcombine", PS_LOG_INFO, "HERM kernel: %s,%d --> %d elements", params, spatialOrder, num);
+    psFree(params);
+
+    // Set the kernel parameters
+    int fullSize = 2 * size + 1;        // Full size of kernels
+    for (int i = 0, index = 0; i < numGaussians; i++) {
+        float sigma = fwhms->data.F32[i] / (2.0 * sqrtf(2.0 * logf(2.0))); // Gaussian sigma
+        // 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++, index++) {
+                psArray *preCalc = psArrayAlloc(3); // Array to hold precalculated values
+                psVector *xKernel = preCalc->data[0] = subtractionKernelHERM(sigma, uOrder, size); // x Kernel
+                psVector *yKernel = preCalc->data[1] = subtractionKernelHERM(sigma, vOrder, size); // y Kernel
+                psKernel *kernel = preCalc->data[2] = psKernelAlloc(-size, size, -size, size);      // Kernel
+
+                // Calculate moments
+                double moment = 0.0;    // Moment, for penalty
+                for (int v = -size, y = 0; v <= size; v++, y++) {
+                    for (int u = -size, x = 0; u <= size; u++, x++) {
+                        double value = xKernel->data.F32[x] * yKernel->data.F32[y]; // Value of kernel
+                        kernel->kernel[v][u] = value;
+                        moment += value * PS_SQR((PS_SQR(u) + PS_SQR(v)));
+                    }
+                }
+
+                // Normalise sum of kernel component to unity for even functions
+                if (uOrder % 2 == 0 && vOrder % 2 == 0) {
+                    double sum = 0.0;   // Sum of kernel component
+                    for (int v = 0; v < fullSize; v++) {
+                        for (int u = 0; u < fullSize; u++) {
+                            sum += xKernel->data.F32[u] * yKernel->data.F32[v];
+                        }
+                    }
+                    sum = 1.0 / sqrt(sum);
+                    psBinaryOp(xKernel, xKernel, "*", psScalarAlloc(sum, PS_TYPE_F32));
+                    psBinaryOp(yKernel, yKernel, "*", psScalarAlloc(sum, PS_TYPE_F32));
+                    psBinaryOp(kernel->image, kernel->image, "*", psScalarAlloc(PS_SQR(sum), PS_TYPE_F32));
+
+#if 1
+                    fprintf(stderr, "%d norm: %lf, null: %f\n", index, sum,kernel->kernel[0][0]);
+#endif
+                    
+                    kernel->kernel[0][0] -= 1.0;
+                    moment *= PS_SQR(sum);
+                }
+
+
+#if 1
                 double sum = 0.0;   // Sum of kernel component
                 for (int v = -size; v <= size; v++) {
@@ -662 +793 @@
       case PM_SUBTRACTION_KERNEL_ISIS:
         return pmSubtractionKernelsISIS(size, spatialOrder, fwhms, orders, penalty, mode);
+      case PM_SUBTRACTION_KERNEL_HERM:
+        return pmSubtractionKernelsHERM(size, spatialOrder, fwhms, orders, penalty, mode);
       case PM_SUBTRACTION_KERNEL_SPAM:
         return pmSubtractionKernelsSPAM(size, spatialOrder, inner, binning, penalty, mode);
@@ -724 +857 @@
     float penalty = 0.0;                // Penalty for wideness
 
-    if (strncmp(description, "ISIS", 4) == 0) {
-        // XXX Support for GUNK
+    // ISIS and HERM have the same description layout
+    if (!strncmp(description, "ISIS", 4) || !strcmp (description, "HERM")) {
+        // XXX Support for GUNK (not yet supported)
         if (strstr(description, "+POIS")) {
             type = PM_SUBTRACTION_KERNEL_GUNK;
             psAbort("Deciphering GUNK kernels (%s) is not currently supported.", description);
-        } else {
-            type = PM_SUBTRACTION_KERNEL_ISIS;
-            char *ptr = (char*)description + 5;    // Eat "ISIS("
-            PARSE_STRING_NUMBER(size, ptr, ',', parseStringInt);
-
-            // Count the number of Gaussians
-            int numGauss = 0;
-            for (char *string = ptr; string; string = strchr(string + 1, '(')) {
-                numGauss++;
-            }
-
-            fwhms = psVectorAlloc(numGauss, PS_TYPE_F32);
-            orders = psVectorAlloc(numGauss, PS_TYPE_S32);
-
-            for (int i = 0; i < numGauss; i++) {
-                ptr++;                  // Eat the '('
-                PARSE_STRING_NUMBER(fwhms->data.F32[i], ptr, ',', parseStringFloat); // Eat "1.234,"
-                PARSE_STRING_NUMBER(orders->data.S32[i], ptr, ')', parseStringInt); // Eat "3)"
-            }
-
-            ptr++;                      // Eat ','
-            PARSE_STRING_NUMBER(spatialOrder, ptr, ',', parseStringInt);
-            penalty = parseStringFloat(ptr);
-        }
-    } else if (strncmp(description, "RINGS", 5) == 0) {
+        } 
+
+        type = pmSubtractionKernelsTypeFromString (description);
+        psAssert (type != PM_SUBTRACTION_KERNEL_NONE, "must  be ISIS or HERM");
+
+        char *ptr = (char*)description + 5;    // Eat "ISIS(" or "HERM("
+        PARSE_STRING_NUMBER(size, ptr, ',', parseStringInt);
+
+        // Count the number of Gaussians
+        int numGauss = 0;
+        for (char *string = ptr; string; string = strchr(string + 1, '(')) {
+            numGauss++;
+        }
+
+        fwhms = psVectorAlloc(numGauss, PS_TYPE_F32);
+        orders = psVectorAlloc(numGauss, PS_TYPE_S32);
+
+        for (int i = 0; i < numGauss; i++) {
+            ptr++;                  // Eat the '('
+            PARSE_STRING_NUMBER(fwhms->data.F32[i], ptr, ',', parseStringFloat); // Eat "1.234,"
+            PARSE_STRING_NUMBER(orders->data.S32[i], ptr, ')', parseStringInt); // Eat "3)"
+        }
+
+        ptr++;                      // Eat ','
+        PARSE_STRING_NUMBER(spatialOrder, ptr, ',', parseStringInt);
+        penalty = parseStringFloat(ptr);
+
+        return pmSubtractionKernelsGenerate(type, size, spatialOrder, fwhms, orders, inner, binning, ringsOrder, penalty, mode);
+    } 
+
+    if (strncmp(description, "RINGS", 5) == 0) {
         type = PM_SUBTRACTION_KERNEL_RINGS;
         char *ptr = (char*)description + 6;
@@ -761 +901 @@
         PARSE_STRING_NUMBER(spatialOrder, ptr, ',', parseStringInt);
         PARSE_STRING_NUMBER(penalty, ptr, ')', parseStringInt);
-    } else {
-        psAbort("Deciphering kernels other than ISIS and RINGS is not currently supported.");
-    }
-
+        return pmSubtractionKernelsGenerate(type, size, spatialOrder, fwhms, orders, inner, binning, ringsOrder, penalty, mode);
+    } 
+
+    psAbort("Deciphering kernels other than ISIS, HERM and RINGS is not currently supported.");
 
     return pmSubtractionKernelsGenerate(type, size, spatialOrder, fwhms, orders,
@@ -778 +918 @@
     if (strcasecmp(type, "ISIS") == 0) {
         return PM_SUBTRACTION_KERNEL_ISIS;
+    }
+    if (strcasecmp(type, "HERM") == 0) {
+        return PM_SUBTRACTION_KERNEL_HERM;
     }
     if (strcasecmp(type, "SPAM") == 0) {

branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionKernels.h

@@ -10 +10 @@
     PM_SUBTRACTION_KERNEL_POIS,         ///< Pan-STARRS Optimal Image Subtraction --- delta functions
     PM_SUBTRACTION_KERNEL_ISIS,         ///< Traditional kernel --- gaussians modified by polynomials
+    PM_SUBTRACTION_KERNEL_HERM,         ///< Hermitian polynomial kernels
     PM_SUBTRACTION_KERNEL_SPAM,         ///< Summed Pixels for Advanced Matching --- summed delta functions
     PM_SUBTRACTION_KERNEL_FRIES,        ///< Fibonacci Radius Increases Excellence of Subtraction
@@ -30 +31 @@
     psString description;               ///< Description of the kernel parameters
     long num;                           ///< Number of kernel components (not including the spatial ones)
-    psVector *u, *v;                    ///< Offset (for POIS) or polynomial order (for ISIS)
-    psVector *widths;                   ///< Gaussian FWHMs (ISIS)
+    psVector *u, *v;                    ///< Offset (for POIS) or polynomial order (for ISIS or HERM)
+    psVector *widths;                   ///< Gaussian FWHMs (ISIS or HERM)
     psVector *uStop, *vStop;            ///< Width of kernel element (SPAM,FRIES only)
-    psArray *preCalc;                   ///< Array of images containing pre-calculated kernel (for ISIS)
+    psArray *preCalc;                   ///< Array of images containing pre-calculated kernel (for ISIS or HERM)
     float penalty;                      ///< Penalty for wideness
     psVector *penalties;                ///< Penalty for each kernel component
@@ -64 +65 @@
         PS_ASSERT_VECTOR_SIZE((KERNELS)->widths, (KERNELS)->num, RETURNVALUE); \
     } \
+    if ((KERNELS)->type == PM_SUBTRACTION_KERNEL_HERM) { \
+        PS_ASSERT_VECTOR_NON_NULL((KERNELS)->widths, RETURNVALUE); \
+        PS_ASSERT_VECTOR_TYPE((KERNELS)->widths, PS_TYPE_F32, RETURNVALUE); \
+        PS_ASSERT_VECTOR_SIZE((KERNELS)->widths, (KERNELS)->num, RETURNVALUE); \
+    } \
     if ((KERNELS)->uStop || (KERNELS)->vStop) { \
         PS_ASSERT_VECTOR_NON_NULL((KERNELS)->uStop, RETURNVALUE); \
@@ -142 +148 @@
                                                );
 
+/// Generate HERM kernels
+pmSubtractionKernels *pmSubtractionKernelsHERM(int size, ///< Half-size of the kernel
+                                               int spatialOrder, ///< Order of spatial variations
+                                               const psVector *fwhms, ///< Gaussian FWHMs
+                                               const psVector *orders, ///< order of hermitian polynomials
+                                               float penalty, ///< Penalty for wideness
+                                               pmSubtractionMode mode ///< Mode for subtraction
+                                               );
+
 /// Generate SPAM kernels
 pmSubtractionKernels *pmSubtractionKernelsSPAM(int size, ///< Half-size of the kernel