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

Timestamp:

Jan 28, 2010, 5:32:15 PM (16 years ago)

Author:

Paul Price

Message:

Attempting to use a wide function for the normalisation component, rather than delta-function.

File:

• branches/pap/psModules/src/imcombine/pmSubtractionEquation.c (modified) (15 diffs)

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

@@ -229 +229 @@
                                       const psArray *convolutions1, // Convolutions of image 1 for each kernel
                                       const psArray *convolutions2, // Convolutions of image 2 for each kernel
+                                      const psKernel *normConv,     // Normalisation, convolved
                                       const pmSubtractionKernels *kernels, // Kernels
                                       const psImage *polyValues, // Spatial polynomial values
@@ -370 +371 @@
         double sumB = 0.0;              // Sum of B products (for matrix, background)
         double sumI2 = 0.0;             // Sum of I_2 (for vector, background)
+
+        double sumAN = 0.0;             // Matrix: A*norm
+        double sumBN = 0.0;             // Matrix: B*norm
+
         for (int y = - footprint; y <= footprint; y++) {
             for (int x = - footprint; x <= footprint; x++) {
                 double a = iConv1->kernel[y][x];
                 double b = iConv2->kernel[y][x];
+                double n = normConv->kernel[y][x];
                 float i1 = image1->kernel[y][x];
                 float i2 = image2->kernel[y][x];
@@ -381 +387 @@
                 double ai1 = a * i1;
                 double bi1 = b * i1;
+
+                double an = a * n;
+                double bn = b * n;
 
                 if (weight) {
@@ -391 +400 @@
                     b *= wtVal;
                     i2 *= wtVal;
+
+                    an *= wtVal;
+                    bn *= wtVal;
                 }
                 if (window) {
@@ -401 +413 @@
                     b *= wtVal;
                     i2 *= wtVal;
+
+                    an *= wtVal;
+                    bn *= wtVal;
                 }
                 sumAI2 += ai2;
@@ -409 +424 @@
                 sumB += b;
                 sumI2 += i2;
+
+                sumAN += an;
+                sumBN += bn;
             }
         }
@@ -420 +438 @@
             double b   = sumB * poly[iTerm];
 
+            double an  = sumAN * poly[iTerm];
+            double bn  = sumBN * poly[iTerm];
+
             if ((mode & PM_SUBTRACTION_EQUATION_NORM) && (mode & PM_SUBTRACTION_EQUATION_KERNELS)) {
-                matrix->data.F64[iIndex][normIndex] = ai1;
-                matrix->data.F64[normIndex][iIndex] = ai1;
-                matrix->data.F64[iIndex + numParams][normIndex] = bi1;
-                matrix->data.F64[normIndex][iIndex + numParams] = bi1;
+//                matrix->data.F64[iIndex][normIndex] = ai1;
+//                matrix->data.F64[normIndex][iIndex] = ai1;
+//                matrix->data.F64[iIndex + numParams][normIndex] = bi1;
+//                matrix->data.F64[normIndex][iIndex + numParams] = bi1;
+                matrix->data.F64[iIndex][normIndex] = an;
+                matrix->data.F64[normIndex][iIndex] = an;
+                matrix->data.F64[iIndex + numParams][normIndex] = bn;
+                matrix->data.F64[normIndex][iIndex + numParams] = bn;
             }
             if ((mode & PM_SUBTRACTION_EQUATION_BG) && (mode & PM_SUBTRACTION_EQUATION_KERNELS)) {
@@ -453 +478 @@
     double sumI1I2 = 0.0;               // Sum of I_1.I_2 (for vector, normalisation)
     double normI1 = 0.0, normI2 = 0.0;  // Sum of I_1 and I_2 within the normalisation window
+
+    double sumN = 0.0;              // Matrix: bg*norm
+    double sumN2 = 0.0;             // Matrix: norm*norm
+    double sumNI2 = 0.0;            // Vector: I2*norm
+
     for (int y = - footprint; y <= footprint; y++) {
         for (int x = - footprint; x <= footprint; x++) {
@@ -461 +491 @@
             double one = 1.0;
             double i1i2 = i1 * i2;
+
+            double n = normConv->kernel[y][x];
+            double n2 = n * n;
+            double ni2 = n * i2;
 
             if (PS_SQR(x) + PS_SQR(y) <= PS_SQR(normWindow)) {
@@ -474 +508 @@
                 i2 *= wtVal;
                 i1i2 *= wtVal;
+
+                n *= wtVal;
+                n2 *= wtVal;
+                ni2 *= wtVal;
             }
             if (window) {
@@ -482 +520 @@
                 i2 *= wtVal;
                 i1i2 *= wtVal;
+
+                n *= wtVal;
+                n2 *= wtVal;
+                ni2 *= wtVal;
             }
             sumI1 += i1;
@@ -488 +530 @@
             sumI2 += i2;
             sumI1I2 += i1i2;
+
+            sumN += n;
+            sumN2 += n2;
+            sumNI2 += ni2;
         }
     }
@@ -496 +542 @@
 
     if (mode & PM_SUBTRACTION_EQUATION_NORM) {
-        matrix->data.F64[normIndex][normIndex] = sumI1I1;
-        vector->data.F64[normIndex] = sumI1I2;
+//        matrix->data.F64[normIndex][normIndex] = sumI1I1;
+//        vector->data.F64[normIndex] = sumI1I2;
+        matrix->data.F64[normIndex][normIndex] = sumN2;
+        vector->data.F64[normIndex] = sumNI2;
+
     }
     if (mode & PM_SUBTRACTION_EQUATION_BG) {
@@ -504 +553 @@
     }
     if ((mode & PM_SUBTRACTION_EQUATION_NORM) && (mode & PM_SUBTRACTION_EQUATION_BG)) {
-        matrix->data.F64[bgIndex][normIndex] = sumI1;
-        matrix->data.F64[normIndex][bgIndex] = sumI1;
+//        matrix->data.F64[bgIndex][normIndex] = sumI1;
+//        matrix->data.F64[normIndex][bgIndex] = sumI1;
+        matrix->data.F64[bgIndex][normIndex] = sumN;
+        matrix->data.F64[normIndex][bgIndex] = sumN;
     }
     return true;
@@ -729 +780 @@
                                            stamp->image1, stamp->image2,
                                            weight, window, stamp->convolutions1, stamp->convolutions2,
-                                           kernels, polyValues, footprint, stamps->normWindow, mode);
+                                           stamp->normConv, kernels, polyValues, footprint,
+                                           stamps->normWindow, mode);
         break;
       default: