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

Timestamp:

Oct 15, 2009, 4:58:42 PM (17 years ago)

Author:

Paul Price

Message:

Fixed implementation of chi^{2 penalties. Using penalties going as r}4 (following Yuan & Akerlof 2008).

File:

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

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

@@ -379 +379 @@
 
 // Add in penalty term to least-squares vector
-static bool calculatePenalty(psVector *vector, // Vector to which to add in penalty term
-                             const pmSubtractionKernels *kernels // Kernel parameters
+static bool calculatePenalty(psImage *matrix,                     // Matrix to which to add in penalty term
+                             psVector *vector,                    // Vector to which to add in penalty term
+                             const pmSubtractionKernels *kernels, // Kernel parameters
+                             float norm                           // Normalisation
     )
 {
@@ -393 +395 @@
         for (int yOrder = 0, index = i; yOrder <= spatialOrder; yOrder++) {
             for (int xOrder = 0; xOrder <= spatialOrder - yOrder; xOrder++, index += numKernels) {
-                vector->data.F64[index] -= penalties->data.F32[i];
+                matrix->data.F64[index][index] -= norm * penalties->data.F32[i];
             }
         }
@@ -786 +788 @@
 #endif
 
-        calculatePenalty(sumVector, kernels);
+        int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index for background
+        calculatePenalty(sumMatrix, sumVector, kernels, sumMatrix->data.F64[bgIndex][bgIndex]);
 
 #ifdef TESTING
@@ -870 +873 @@
             }
         }
-        calculatePenalty(sumVector1, kernels);
-        calculatePenalty(sumVector2, kernels);
+
+        int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index for background
+        calculatePenalty(sumMatrix1, sumVector1, kernels, sumMatrix1->data.F64[bgIndex][bgIndex]);
+        calculatePenalty(sumMatrix2, sumVector2, kernels, -sumMatrix1->data.F64[bgIndex][bgIndex]);
 
         // Pure matrix operations