Changeset 26502

Timestamp:

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

Author:

eugene

Message:

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

Location:

branches/eam_branches/20091201/psModules/src/imcombine

Files:

• pmSubtractionEquation.c (modified) (11 diffs)
• pmSubtractionKernels.c (modified) (11 diffs)
• pmSubtractionStamps.c (modified) (7 diffs)
• pmSubtractionVisual.c (modified) (8 diffs)

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

@@ -15 +15 @@
 #include "pmSubtractionVisual.h"
 
-//#define TESTING                         // TESTING output for debugging; may not work with threads!
+// #define TESTING                         // TESTING output for debugging; may not work with threads!
 
 #define USE_WEIGHT                      // Include weight (1/variance) in equation?
-#define USE_WINDOW                      // Include weight (1/variance) in equation?
+// #define USE_WINDOW                      // Include weight (1/variance) in equation?
 
 
@@ -157 +157 @@
             if (mode & PM_SUBTRACTION_EQUATION_KERNELS) {
                 vector->data.F64[iIndex] = sumIC * poly[iTerm];
-                if (!(mode & PM_SUBTRACTION_EQUATION_NORM)) {
+                // XXX TEST for Hermitians: do not calculate A - norm*B - \sum(k x B),
+                // instead, calculate A - \sum(k x B), with full hermitians
+                if (0 && !(mode & PM_SUBTRACTION_EQUATION_NORM)) {
                     // subtract norm * sumRC * poly[iTerm]
                     psAssert (kernels->solution1, "programming error: define solution first!");
@@ -228 +230 @@
                                       const psKernel *image2, // Image 2
                                       const psKernel *weight,  // Weight image
+                                      const psKernel *window,  // Window image
                                       const psArray *convolutions1, // Convolutions of image 1 for each kernel
                                       const psArray *convolutions2, // Convolutions of image 2 for each kernel
@@ -291 +294 @@
                         ab *= wtVal;
                     }
+                    if (window) {
+                        float wtVal = window->kernel[y][x];
+                        aa *= wtVal;
+                        bb *= wtVal;
+                        ab *= wtVal;
+                    }
                     sumAA += aa;
                     sumBB += bb;
@@ -319 +328 @@
                     if (weight) {
                         ab *= weight->kernel[y][x];
+                    }
+                    if (window) {
+                        ab *= window->kernel[y][x];
                     }
                     sumAB += ab;
@@ -365 +377 @@
                     i2 *= wtVal;
                 }
+                if (window) {
+                    float wtVal = window->kernel[y][x];
+                    ai2 *= wtVal;
+                    bi2 *= wtVal;
+                    ai1 *= wtVal;
+                    bi1 *= wtVal;
+                    i1i2 *= wtVal;
+                    a *= wtVal;
+                    b *= wtVal;
+                    i2 *= wtVal;
+                }
                 sumAI2 += ai2;
                 sumBI2 += bi2;
@@ -411 +434 @@
             if (weight) {
                 float wtVal = weight->kernel[y][x];
+                i1 *= wtVal;
+                i1i1 *= wtVal;
+                one *= wtVal;
+                i2 *= wtVal;
+                i1i2 *= wtVal;
+            }
+            if (window) {
+                float wtVal = window->kernel[y][x];
                 i1 *= wtVal;
                 i1i1 *= wtVal;
@@ -516 +547 @@
             for (int xOrder = 0; xOrder <= spatialOrder - yOrder; xOrder++, index += numKernels) {
                 // Contribution to chi^2: a_i^2 P_i
+                if (!isfinite(penalties->data.F32[i])) {
+                    psAbort ("invalid penalty");
+                }
                 matrix->data.F64[index][index] -= norm * penalties->data.F32[i];
             }
@@ -702 +736 @@
         break;
       case PM_SUBTRACTION_MODE_DUAL:
-        psAbort ("dual is disabled: need to add window and calculation mode");
         if (new) {
             stamp->matrix2 = psImageAlloc(numKernels * numSpatial, numKernels * numSpatial, PS_TYPE_F64);
@@ -714 +747 @@
 #endif
         status = calculateDualMatrixVector(stamp->matrix1, stamp->vector1, stamp->matrix2, stamp->vector2,
-                                           stamp->matrixX, stamp->image1, stamp->image2, weight,
+                                           stamp->matrixX, stamp->image1, stamp->image2, weight, window,
                                            stamp->convolutions1, stamp->convolutions2, kernels, polyValues,
                                            footprint);
@@ -1587 +1620 @@
                 }
             }
+
+            // XXX visualize the target, source, convolution and residual
+            pmSubtractionVisualShowFitAddStamp (image2, image1, residual, background, norm, i);
+
             for (int y = - footprint; y <= footprint; y++) {
                 for (int x = - footprint; x <= footprint; x++) {

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,

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

@@ -81 +81 @@
 
 // Search a region for a suitable stamp
-bool stampSearch(int *xStamp, int *yStamp, // Coordinates of stamp, to return
+bool stampSearch(float *xStamp, float *yStamp, // Coordinates of stamp, to return
                  float *fluxStamp, // Flux of stamp, to return
                  const psImage *image1, const psImage *image2, // Images to search
@@ -116 +116 @@
                 ((image1) ? image1->data.F32[y][x] : image2->data.F32[y][x]); // Flux at pixel
             if (flux > *fluxStamp) {
-                *xStamp = x;
-                *yStamp = y;
+                *xStamp = x + 0.5;
+                *yStamp = y + 0.5;
                 *fluxStamp = flux;
                 found = true;
@@ -405 +405 @@
             numSearch++;
 
-            int xStamp = 0, yStamp = 0; // Coordinates of stamp
+            float xStamp = 0.0, yStamp = 0.0; // Coordinates of stamp
             float fluxStamp = -INFINITY; // Flux of stamp
             bool goodStamp = false;     // Found a good stamp?
@@ -437 +437 @@
                     goodStamp = stampSearch(&xStamp, &yStamp, &fluxStamp, image1, image2, thresh1, thresh2,
                                             subMask, xMin, xMax, yMin, yMax, numCols, numRows, border);
+
+                    // XXX reset for a test:
+                    xStamp = xList->data.F32[j];
+                    yStamp = yList->data.F32[j];
+                    fprintf (stderr, "find: %d %d ==> %5.1f %5.1f\n", xCentre, yCentre, xStamp, yStamp);
                 }
             } else {
@@ -543 +548 @@
         }
 
+        fprintf (stderr, "stamp: %5.1f %5.1f == %d %d\n", xStamp, yStamp, xPix, yPix);
+
         bool found = false;
         for (int j = 0; j < numStamps && !found; j++) {
@@ -735 +742 @@
             return false;
         }
+        fprintf (stderr, "stamp: %5.1f %5.1f == %d %d\n", stamp->x, stamp->y, x, y);
 
         // Catch memory leaks --- these should have been freed and NULLed before
@@ -811 +819 @@
             y->data.F32[numOut] = source->peak->yf;
         }
+        fprintf (stderr, "stamp: %5.1f %5.1f\n", x->data.F32[numOut], y->data.F32[numOut]);
         numOut++;
     }

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

@@ -251 +251 @@
         int yPix = ySub * (2*footprint + 1 + 3) + footprint;
 
+        double sum = 0.0;
         for (int y = -footprint; y <= footprint; y++) {
             for (int x = -footprint; x <= footprint; x++) {
                 output->data.F32[y + yPix][x + xPix] = kernel->kernel[y][x];
+                sum += kernel->kernel[y][x];
             }
         }
+        fprintf (stderr, "kernel %d, sum %f\n", i, sum);
     }                                                    
         
@@ -274 +277 @@
 
     // choose the brightest stamp
-    pmSubtractionStamp *maxStamp = stamps->stamps->data[0];
-    float maxFlux = maxStamp->flux; 
-    for (int i = 1; i < stamps->num; i++) {
+    pmSubtractionStamp *maxStamp = NULL;
+    float maxFlux = NAN;
+    for (int i = 0; i < stamps->num; i++) {
         pmSubtractionStamp *stamp = stamps->stamps->data[i];
+        if (!isfinite(stamp->flux)) continue;
+        if (!stamp->convolutions1 && !stamp->convolutions2) continue;
+        if (!maxStamp) {
+            maxFlux = stamp->flux;
+            maxStamp = stamp;
+            continue;
+        }
         if (stamp->flux > maxFlux) {
             maxFlux = stamp->flux;
             maxStamp = stamp;
         }
+    }
+
+    if (!isfinite(maxStamp->flux)) {
+        fprintf (stderr, "no valid stamps?\n");
     }
 
@@ -308 +322 @@
             int yPix = ySub * (2*footprint + 1 + 3) + footprint;
             
+            double sum = 0.0;
             for (int y = -footprint; y <= footprint; y++) {
                 for (int x = -footprint; x <= footprint; x++) {
                     output->data.F32[y + yPix][x + xPix] = kernel->kernel[y][x];
+                    sum += kernel->kernel[y][x];
                 }
             }
+            fprintf (stderr, "kernel %d, sum %f\n", i, sum);
         }               
         pmVisualScaleImage(kapa2, output, "Image", 0, true);
@@ -339 +356 @@
             int yPix = ySub * (2*footprint + 1 + 3) + footprint;
             
+            double sum = 0.0;
             for (int y = -footprint; y <= footprint; y++) {
                 for (int x = -footprint; x <= footprint; x++) {
                     output->data.F32[y + yPix][x + xPix] = kernel->kernel[y][x];
+                    sum += kernel->kernel[y][x];
                 }
             }
+            fprintf (stderr, "kernel %d, sum %f\n", i, sum);
         }               
         pmVisualScaleImage(kapa2, output, "Image", 1, true);
@@ -375 +395 @@
         overlay[Noverlay].type = KII_OVERLAY_BOX;
         if ((stamp->x < 1.0) && (stamp->y < 1.0)) {
-            fprintf (stderr, "stamp zero: %f %f\n", stamp->x, stamp->y);
+            // fprintf (stderr, "stamp zero: %f %f\n", stamp->x, stamp->y);
             continue;
         }
         if (!isfinite(stamp->x) && !isfinite(stamp->y)) {
-            fprintf (stderr, "stamp nan: %f %f\n", stamp->x, stamp->y);
+            // fprintf (stderr, "stamp nan: %f %f\n", stamp->x, stamp->y);
             continue;
         }
@@ -409 +429 @@
 bool pmSubtractionVisualShowFitInit(pmSubtractionStampList *stamps) {
 
-    if (!pmVisualIsVisual()) return true;
+    // if (!pmVisualIsVisual()) return true;
 
     // generate 4 storage images large enough to hold the N stamps:
@@ -446 +466 @@
 bool pmSubtractionVisualShowFitAddStamp(psKernel *target, psKernel *source, psKernel *convolution, double background, double norm, int index) {
 
-    if (!pmVisualIsVisual()) return true;
+    // if (!pmVisualIsVisual()) return true;
 
     double sum;
@@ -517 +537 @@
 bool pmSubtractionVisualShowFit() {
 
-    if (!pmVisualIsVisual()) return true;
+    // if (!pmVisualIsVisual()) return true;
     if (!pmVisualInitWindow(&kapa1, "ppSub:Images")) return false;
     if (!pmVisualInitWindow(&kapa2, "ppSub:Misc")) return false;