Changeset 26502 for branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionEquation.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/pmSubtractionEquation.c (modified) (11 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++) {