Changeset 26321 for branches/pap/psModules/src/imcombine/pmSubtractionEquation.c

Timestamp:

Dec 2, 2009, 12:04:51 PM (16 years ago)

Author:

Paul Price

Message:

Merging Gene's work on adding a window function to the stamps. I'd been thinking about something like this for a while, so glad someone's done it. It's got visualisation turned on though.

Location:

branches/pap/psModules/src/imcombine

Files:

• . (modified) (1 prop)
• pmSubtractionEquation.c (modified) (11 diffs)

branches/pap/psModules/src/imcombine/pmSubtractionEquation.c

@@ -18 +18 @@
 
 #define USE_WEIGHT                      // Include weight (1/variance) in equation?
+#define USE_WINDOW                      // Include weight (1/variance) in equation?
 
 
@@ -30 +31 @@
                                   const psKernel *reference, // Reference image (convolution source)
                                   const psKernel *weight,  // Weight image
+                                  const psKernel *window,  // Window image
                                   const psArray *convolutions,         // Convolutions for each kernel
                                   const pmSubtractionKernels *kernels, // Kernels
@@ -76 +78 @@
 #ifdef USE_WEIGHT
                     cc *= weight->kernel[y][x];
+#endif
+#ifdef USE_WINDOW
+                    if (window) {
+                        cc *= window->kernel[y][x];
+                    }
 #endif
                     sumCC += cc;
@@ -107 +114 @@
                 rc *= wtVal;
                 c *= wtVal;
+#endif
+#ifdef USE_WINDOW
+                if (window) {
+                    float winVal = window->kernel[y][x];
+                    ic *= winVal;
+                    rc *= winVal;
+                    c  *= winVal;
+                }
 #endif
                 sumIC += ic;
@@ -144 +159 @@
             ref *= wtVal;
             one *= wtVal;
+#endif
+#ifdef USE_WINDOW
+            if (window) {
+                float  winVal = window->kernel[y][x];
+                rr      *= winVal;
+                ir      *= winVal;
+                in      *= winVal;
+                ref *= winVal;
+                one *= winVal;
+            }
 #endif
             sumRR += rr;
@@ -603 +628 @@
 #endif
 
+    // XXX visualize the set of convolved stamps
+
     psImage *polyValues = p_pmSubtractionPolynomial(NULL, spatialOrder,
                                                     stamp->xNorm, stamp->yNorm); // Polynomial terms
@@ -620 +647 @@
       case PM_SUBTRACTION_MODE_1:
         status = calculateMatrixVector(stamp->matrix1, stamp->vector1, stamp->image2, stamp->image1,
-                                       stamp->weight, stamp->convolutions1, kernels, polyValues,
-                                       footprint);
+                                       stamp->weight, stamps->window, stamp->convolutions1, kernels,
+                                       polyValues, footprint);
         break;
       case PM_SUBTRACTION_MODE_2:
         status = calculateMatrixVector(stamp->matrix1, stamp->vector1, stamp->image1, stamp->image2,
-                                       stamp->weight, stamp->convolutions2, kernels, polyValues,
-                                       footprint);
+                                       stamp->weight, stamps->window, stamp->convolutions2, kernels,
+                                       polyValues, footprint);
         break;
       case PM_SUBTRACTION_MODE_DUAL:
@@ -1147 +1174 @@
     double devNorm = 1.0 / (double)numPixels; // Normalisation for deviations
     int numKernels = kernels->num;      // Number of kernels
+    double norm = NAN;
 
     psImage *polyValues = NULL;         // Polynomial values
     psKernel *residual = psKernelAlloc(-footprint, footprint, -footprint, footprint); // Residual image
+
+    // set up holding images for the visualization
+    pmSubtractionVisualShowFitInit (stamps);
 
     for (int i = 0; i < stamps->num; i++) {
@@ -1160 +1191 @@
         // Calculate coefficients of the kernel basis functions
         polyValues = p_pmSubtractionPolynomial(polyValues, kernels->spatialOrder, stamp->xNorm, stamp->yNorm);
-        double norm = p_pmSubtractionSolutionNorm(kernels); // Normalisation
+        norm = p_pmSubtractionSolutionNorm(kernels); // Normalisation
         double background = p_pmSubtractionSolutionBackground(kernels, polyValues);// Difference in background
 
@@ -1214 +1245 @@
                 }
             }
+
+            // XXX visualize the target, source, convolution and residual
+            pmSubtractionVisualShowFitAddStamp (target, source, residual, background, norm, i);
+
             for (int y = - footprint; y <= footprint; y++) {
                 for (int x = - footprint; x <= footprint; x++) {
@@ -1300 +1335 @@
         }
 #endif
-
+    
     }
     psFree(residual);
     psFree(polyValues);
 
+    psLogMsg("psModules.imcombine", PS_LOG_INFO, "normalization: %f", norm);
+    pmSubtractionVisualShowFit();
+    pmSubtractionVisualPlotFit(kernels);
+
     return deviations;
 }