Changeset 26502 for branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionKernels.c

Timestamp:

Jan 2, 2010, 6:13:37 PM (17 years ago)

Author:

eugene

Message:

adding window to dual convolution; fixing normalization for deconvolved-hermitians

File:

• branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionKernels.c (modified) (11 diffs)

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

@@ -132 +132 @@
             kernels->preCalc->data[index] = NULL;
             kernels->penalties->data.F32[index] = kernels->penalty * PS_SQR(PS_SQR(u) + PS_SQR(v));
+            if (!isfinite(kernels->penalties->data.F32[index])) {
+                psAbort ("invalid penalty");
+            }
 
             psTrace("psModules.imcombine", 7, "Kernel %d: %d %d\n", index, u, v);
@@ -146 +149 @@
 }
 
-bool pmSubtractionKernelPreCalcNormalize (pmSubtractionKernels *kernels, pmSubtractionKernelPreCalc *preCalc, int index, int size, int uOrder, int vOrder, float fwhm) {
+bool pmSubtractionKernelPreCalcNormalize (pmSubtractionKernels *kernels, pmSubtractionKernelPreCalc *preCalc, int index, int size, int uOrder, int vOrder, float fwhm, bool AlardLuptonStyle) {
 
     // Calculate moments
@@ -157 +160 @@
     }
 
-    // Normalize 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 = -size; v <= size; v++) {
-            for (int u = -size; u <= size; u++) {
-                sum += preCalc->kernel->kernel[v][u];
-            }
-        }
-        sum = 1.0 / sqrt(sum);
-        psBinaryOp(preCalc->xKernel, preCalc->xKernel, "*", psScalarAlloc(sum, PS_TYPE_F32));
-        psBinaryOp(preCalc->yKernel, preCalc->yKernel, "*", psScalarAlloc(sum, PS_TYPE_F32));
-        psBinaryOp(preCalc->kernel->image, preCalc->kernel->image, "*", psScalarAlloc(PS_SQR(sum), PS_TYPE_F32));
-
-#if 1
-        fprintf(stderr, "%d norm: %lf, null: %f\n", index, sum, preCalc->kernel->kernel[0][0]);
-#endif
-                    
-        preCalc->kernel->kernel[0][0] -= 1.0;
-        moment *= PS_SQR(sum);
-    }
-
-#if 1
+    // we have 4 cases here:
+    // 1) for odd functions, normalize the kernel by the maximum swing / Npix
+    // 2) for even functions, normalize the kernel to unity
+    // 3) for alard-lupton style normalization, subtract 1 from the 0,0 pixel for all even functions
+    // 4) for deconvolved hermitians, subtract 1 from the 0,0 pixel for the 0,0 function(s)
+
     double sum = 0.0;   // Sum of kernel component
+    double min = FLT_MAX;
+    double max = FLT_MIN;
+
     for (int v = -size; v <= size; v++) {
         for (int u = -size; u <= size; u++) {
             sum += preCalc->kernel->kernel[v][u];
+            min = PS_MIN(preCalc->kernel->kernel[v][u], min);
+            max = PS_MAX(preCalc->kernel->kernel[v][u], max);
         }
     }
-    fprintf(stderr, "%d sum: %lf\n", index, sum);
+#if 1
+    fprintf(stderr, "%d raw: %lf, null: %f, min: %lf, max: %lf, moment: %lf\n", index, sum, preCalc->kernel->kernel[0][0], min, max, moment);
+#endif
+
+    // only even terms have non-zero sums
+    if ((uOrder % 2 == 0) && (vOrder % 2 == 0)) {
+        moment /= sum;
+    } else {
+        moment = 0.0;
+    }
+
+    bool zeroNull = false;
+    float scale1D = 1.0 / sqrt(sum);
+    float scale2D = 1.0 / sum;
+
+    if (AlardLuptonStyle && (uOrder % 2 == 0 && vOrder % 2 == 0)) {
+        zeroNull = true;
+    }
+    if (!AlardLuptonStyle && (uOrder == 0 && vOrder == 0)) {
+        zeroNull = true;
+    }
+    if ((uOrder % 2) || (vOrder % 2)) {
+        // scale2D = 1.0 / (preCalc->kernel->image->numCols * preCalc->kernel->image->numRows * max);
+        scale2D = 1.0 / max;
+        scale1D = sqrt(scale2D);
+    }
+
+    psBinaryOp(preCalc->xKernel, preCalc->xKernel, "*", psScalarAlloc(scale1D, PS_TYPE_F32));
+    psBinaryOp(preCalc->yKernel, preCalc->yKernel, "*", psScalarAlloc(scale1D, PS_TYPE_F32));
+    psBinaryOp(preCalc->kernel->image, preCalc->kernel->image, "*", psScalarAlloc(scale2D, PS_TYPE_F32));
+    if (zeroNull) {
+        preCalc->kernel->kernel[0][0] -= 1.0;
+    }
+
+#if 1
+    sum = 0.0;   // Sum of kernel component
+    min = FLT_MAX;
+    max = FLT_MIN;
+    for (int v = -size; v <= size; v++) {
+        for (int u = -size; u <= size; u++) {
+            sum += preCalc->kernel->kernel[v][u];
+            min = PS_MIN(preCalc->kernel->kernel[v][u], min);
+            max = PS_MAX(preCalc->kernel->kernel[v][u], max);
+        }
+    }
+    fprintf(stderr, "%d mod: %lf, null: %f, min: %lf, max: %lf, scale: %f\n", index, sum, preCalc->kernel->kernel[0][0], min, max, scale2D);
 #endif
 
@@ -196 +233 @@
     kernels->preCalc->data[index] = preCalc;
     kernels->penalties->data.F32[index] = kernels->penalty * fabsf(moment);
+    if (!isfinite(kernels->penalties->data.F32[index])) {
+        psAbort ("invalid penalty");
+    }
 
     psTrace("psModules.imcombine", 7, "Kernel %d: %f %d %d %f\n", index,
@@ -250 +290 @@
 
                 pmSubtractionKernelPreCalc *preCalc = pmSubtractionKernelPreCalcAlloc(PM_SUBTRACTION_KERNEL_ISIS, uOrder, vOrder, size, sigma); // structure to hold precalculated values
-                pmSubtractionKernelPreCalcNormalize (kernels, preCalc, index, size, uOrder, vOrder, fwhms->data.F32[i]);
+                pmSubtractionKernelPreCalcNormalize (kernels, preCalc, index, size, uOrder, vOrder, fwhms->data.F32[i], true);
             }
         }
@@ -302 +342 @@
             for (int vOrder = 0; vOrder <= orders->data.S32[i] - uOrder; vOrder++, index++) {
                 pmSubtractionKernelPreCalc *preCalc = pmSubtractionKernelPreCalcAlloc(PM_SUBTRACTION_KERNEL_HERM, uOrder, vOrder, size, sigma); // structure to hold precalculated values
-                pmSubtractionKernelPreCalcNormalize (kernels, preCalc, index, size, uOrder, vOrder, fwhms->data.F32[i]);
+                pmSubtractionKernelPreCalcNormalize (kernels, preCalc, index, size, uOrder, vOrder, fwhms->data.F32[i], true);
             }
         }
@@ -338 +378 @@
         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;
+        num += PS_SQR(gaussOrder + 1);
     }
 
@@ -349 +389 @@
     // XXXXX hard-wired reference sigma for now of 1.7 pix (== 4.0 pix fwhm == 1.0 arcsec in simtest)
     // generate the Gaussian deconvolution kernel
-    # define DECONV_SIGMA 1.7
+    # define DECONV_SIGMA 1.6
     psKernel *kernelGauss = pmSubtractionDeconvolveGauss (size, DECONV_SIGMA);
+
+# if 1
+    psArray *deconKernels = psArrayAllocEmpty(100);
+# endif
 
     // Set the kernel parameters
@@ -357 +401 @@
         // 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++) {
+            for (int vOrder = 0; vOrder <= orders->data.S32[i]; vOrder++, index++) {
 
                 pmSubtractionKernelPreCalc *preCalc = pmSubtractionKernelPreCalcAlloc(PM_SUBTRACTION_KERNEL_HERM, uOrder, vOrder, size, sigma); // structure to hold precalculated values
@@ -365 +409 @@
                 preCalc->kernel = pmSubtractionDeconvolveKernel(kernelTarget, kernelGauss); // Kernel
 
+                pmSubtractionKernelPreCalcNormalize (kernels, preCalc, index, size, uOrder, vOrder, fwhms->data.F32[i], false);
+
                 // XXXX test demo that deconvolved kernel is valid
-                // psImage *kernelConv = psImageConvolveFFT(NULL, preCalc->kernel->image, NULL, 0, kernelGauss);
-                // pmSubtractionVisualShowSubtraction (kernelTarget->image, preCalc->kernel->image, kernelConv);
-
-                pmSubtractionKernelPreCalcNormalize (kernels, preCalc, index, size, uOrder, vOrder, fwhms->data.F32[i]);
+# if 1
+                psImage *kernelConv = psImageConvolveFFT(NULL, preCalc->kernel->image, NULL, 0, kernelGauss);
+                psArrayAdd (deconKernels, 100, kernelConv);
+                psFree (kernelConv);
+
+                if (!uOrder && !vOrder){
+                    pmSubtractionVisualShowSubtraction (kernelTarget->image, preCalc->kernel->image, kernelConv);
+                }
+# endif
             }
         }
     }
+
+# if 1
+    psImage *dot = psImageAlloc(deconKernels->n, deconKernels->n, PS_TYPE_F32);
+    for (int i = 0; i < deconKernels->n; i++) {
+        for (int j = 0; j <= i; j++) {
+            psImage *t1 = deconKernels->data[i];
+            psImage *t2 = deconKernels->data[j];
+
+            double sum = 0.0;
+            for (int iy = 0; iy < t1->numRows; iy++) {
+                for (int ix = 0; ix < t1->numCols; ix++) {
+                    sum += t1->data.F32[iy][ix] * t2->data.F32[iy][ix];
+                }
+            }
+            dot->data.F32[j][i] = sum;
+            dot->data.F32[i][j] = sum;
+        }
+    }
+    pmSubtractionVisualShowSubtraction (dot, NULL, NULL);
+    psFree (dot);
+    psFree (deconKernels);
+# endif
 
     return kernels;
@@ -857 +930 @@
                 kernels->v->data.S32[index] = vOrder;
                 kernels->penalties->data.F32[index] = kernels->penalty * fabsf(moment);
+                if (!isfinite(kernels->penalties->data.F32[index])) {
+                    psAbort ("invalid penalty");
+                }
 
                 psTrace("psModules.imcombine", 7, "Kernel %d: %d %d %d\n", index,