Changeset 29506

Timestamp:

Oct 21, 2010, 6:14:32 AM (16 years ago)

Author:

eugene

Message:

subtract kernel(0) instead of delta function

Location:

branches/eam_branches/ipp-20100823/psModules/src/imcombine

Files:

• pmSubtractionEquation.c (modified) (5 diffs)
• pmSubtractionKernels.c (modified) (2 diffs)

branches/eam_branches/ipp-20100823/psModules/src/imcombine/pmSubtractionEquation.c

@@ -21 +21 @@
 
 // XXX TEST:
-# define USE_WINDOW                      // window to avoid neighbor contamination
+//# define USE_WINDOW                      // window to avoid neighbor contamination
 
 # define PENALTY false
@@ -1091 +1091 @@
         }
 
-#if 1
+#if 0
         psImage *save = psImageCopy(NULL, sumMatrix, PS_TYPE_F32);
         psFitsWriteImageSimple ("sumMatrix.fits", save, NULL);
@@ -1115 +1115 @@
         }
 
+# if (0)
         for (int i = 0; i < solution->n; i++) {
             psLogMsg("psModules.imcombine", PS_LOG_DETAIL, "Single solution %d: %lf\n", i, solution->data.F64[i]);
         }
+# endif
 
         if (!kernels->solution1) {
@@ -1169 +1171 @@
         }
 
-#if 1
+#if 0
         psImage *save = psImageCopy(NULL, sumMatrix, PS_TYPE_F32);
         psFitsWriteImageSimple ("sumMatrix.fits", save, NULL);
@@ -1187 +1189 @@
         solution = psMatrixSolveSVD(solution, sumMatrix, sumVector, NAN);
 
-#if (1)
+#if (0)
         for (int i = 0; i < solution->n; i++) {
             fprintf(stderr, "Dual solution %d: %lf\n", i, solution->data.F64[i]);

branches/eam_branches/ipp-20100823/psModules/src/imcombine/pmSubtractionKernels.c

@@ -183 +183 @@
 }
 
+# define CENTRAL_DELTA 0
+
+# if (CENTRAL_DELTA)
+
 // XXX *** this code used the central pixel to force zero net flux,
 // Alard actually uses kernel(0) for this purpose (for even order, kernel[i] = kernel'[i] - kernel[0])
-
 static bool pmSubtractionKernelPreCalcNormalize(pmSubtractionKernels *kernels, pmSubtractionKernelPreCalc *preCalc,
                                                 int index, int uOrder, int vOrder, float fwhm,
@@ -290 +293 @@
     return true;
 }
+
+# else /* CENTRAL_DELTA */
+
+static bool zeroIsNormal = false;
+
+// this code uses kernel(0) to force zero flux, and is invalid for other kinds of normalizations
+static bool pmSubtractionKernelPreCalcNormalize(pmSubtractionKernels *kernels, pmSubtractionKernelPreCalc *preCalc,
+                                                int index, int uOrder, int vOrder, float fwhm,
+                                                bool AlardLuptonStyle, bool forceZeroNull)
+{
+    // 1) for odd functions: no renormalization 
+    // 2) for even functions, normalize the kernel to unity
+    // 3) for even functions & index > 0, subtract kernel(0)
+
+    // Calculate moments
+    double sum = 0.0, sum2 = 0.0;           // Sum of kernel component
+    float min = INFINITY, max = -INFINITY;  // Minimum and maximum kernel value
+
+    for (int v = preCalc->kernel->yMin; v <= preCalc->kernel->yMax; v++) {
+        for (int u = preCalc->kernel->xMin; u <= preCalc->kernel->xMax; u++) {
+            double value = preCalc->kernel->kernel[v][u];
+            double value2 = PS_SQR(value);
+            sum += value;
+            sum2 += value2;
+            min = PS_MIN(value, min);
+            max = PS_MAX(value, max);
+        }
+    }
+
+#if 0
+    fprintf(stderr, "%d raw: %lf, null: %f, min: %lf, max: %lf\n", index, sum, preCalc->kernel->kernel[0][0], min, max);
+#endif
+
+    float scale2D = 1.0;                // Scaling for 2-D kernels
+    float scale1D = 1.0;                // Scaling for 1-D kernels
+
+    if (uOrder % 2 == 0 && vOrder % 2 == 0) {
+
+        scale2D = 1.0 / sum;            // Scaling for 2-D kernels
+        scale1D = sqrtf(scale2D);               // Scaling for 1-D kernels
+
+        for (int v = preCalc->kernel->yMin; v <= preCalc->kernel->yMax; v++) {
+            for (int u = preCalc->kernel->xMin; u <= preCalc->kernel->xMax; u++) {
+                preCalc->kernel->kernel[v][u] *= scale2D;
+            }
+        }
+        if (index == 0) {
+            zeroIsNormal = true;
+        } else {
+            psAssert(zeroIsNormal, "failed to normalize zero kernel first");
+            pmSubtractionKernelPreCalc *zeroKernel = kernels->preCalc->data[0];
+            psAssert(zeroKernel, "failed to supply zero kernel");
+            for (int v = preCalc->kernel->yMin; v <= preCalc->kernel->yMax; v++) {
+                for (int u = preCalc->kernel->xMin; u <= preCalc->kernel->xMax; u++) {
+                    preCalc->kernel->kernel[v][u] -= zeroKernel->kernel->kernel[v][u];
+                }
+            }
+        }
+
+        // XXX why do we bother renormlizing the 1D kernels?  I don't think we use them again...
+        if (preCalc->xKernel) {
+            psBinaryOp(preCalc->xKernel, preCalc->xKernel, "*", psScalarAlloc(scale1D, PS_TYPE_F32));
+        }
+        if (preCalc->yKernel) {
+            psBinaryOp(preCalc->yKernel, preCalc->yKernel, "*", psScalarAlloc(scale1D, PS_TYPE_F32));
+        }
+    }
+
+#if 0
+    {
+        double Sum = 0.0;   // Sum of kernel component
+        double Sum2 = 0.0;   // Sum of kernel component
+        float min = INFINITY, max = -INFINITY;  // Minimum and maximum kernel value
+        for (int v = preCalc->kernel->yMin; v <= preCalc->kernel->yMax; v++) {
+            for (int u = preCalc->kernel->xMin; u <= preCalc->kernel->xMax; u++) {
+                double value = preCalc->kernel->kernel[v][u];
+                Sum += value;
+                Sum2 += PS_SQR(value);
+                min = PS_MIN(preCalc->kernel->kernel[v][u], min);
+                max = PS_MAX(preCalc->kernel->kernel[v][u], max);
+            }
+        }
+        fprintf(stderr, "%d sum: %lf, sum2: %lf, null: %f, min: %lf, max: %lf, scale: %f\n", index, Sum, Sum2, preCalc->kernel->kernel[0][0], min, max, scale2D);
+    }
+#endif
+
+    kernels->widths->data.F32[index] = fwhm;
+    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;
+    psTrace("psModules.imcombine", 7, "Kernel %d: %f %d %d\n", index, fwhm, uOrder, vOrder);
+
+    return true;
+}
+# endif /* Central Delta */
 
 pmSubtractionKernels *p_pmSubtractionKernelsRawISIS(int size, int spatialOrder,