Changeset 18248 for branches/pap_branch_080617/psModules/src/imcombine/pmSubtractionEquation.c

Timestamp:

Jun 20, 2008, 3:27:33 PM (18 years ago)

Author:

Paul Price

Message:

Dual convolution seems to be working. Removed division by 'weight' when generating equation --- seems to bias the fits.

File:

• branches/pap_branch_080617/psModules/src/imcombine/pmSubtractionEquation.c (modified) (12 diffs)

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

@@ -30 +30 @@
     for (int y = - footprint; y <= footprint; y++) {
         for (int x = - footprint; x <= footprint; x++) {
-            sum += image1->kernel[y][x] * image2->kernel[y][x] / weight->kernel[y][x];
+            sum += image1->kernel[y][x] * image2->kernel[y][x] / 1.0; // weight->kernel[y][x];
         }
     }
@@ -194 +194 @@
             for (int y = - footprint; y <= footprint; y++) {
                 for (int x = - footprint; x <= footprint; x++) {
-                    sumC += conv->kernel[y][x] / weight->kernel[y][x];
+                    sumC += conv->kernel[y][x] / 1.0; // weight->kernel[y][x];
                 }
             }
@@ -217 +217 @@
         for (int y = - footprint; y <= footprint; y++) {
             for (int x = - footprint; x <= footprint; x++) {
-                double invNoise2 = 1.0 / weight->kernel[y][x];
+                double invNoise2 = 1.0 / 1.0; // weight->kernel[y][x];
                 double value = input->kernel[y][x] * invNoise2;
                 sumI += value;
@@ -276 +276 @@
         for (int y = - footprint; y <= footprint; y++) {
             for (int x = - footprint; x <= footprint; x++) {
-                    sumTC += target->kernel[y][x] * conv->kernel[y][x] / weight->kernel[y][x];
+                sumTC += target->kernel[y][x] * conv->kernel[y][x] / 1.0; // weight->kernel[y][x];
             }
         }
@@ -296 +296 @@
         for (int y = - footprint; y <= footprint; y++) {
             for (int x = - footprint; x <= footprint; x++) {
-                float value = target->kernel[y][x] / weight->kernel[y][x];
+                float value = target->kernel[y][x] / 1.0; // weight->kernel[y][x];
                 sumIT += value * input->kernel[y][x];
                 sumT += value;
@@ -365 +365 @@
         for (int y = - footprint; y <= footprint; y++) {
             for (int x = - footprint; x <= footprint; x++) {
-                sumC += conv->kernel[y][x] / weight->kernel[y][x];
+                sumC += conv->kernel[y][x] / 1.0; // weight->kernel[y][x];
             }
         }
@@ -384 +384 @@
 
 // Add in penalty term to least-squares vector
-static bool calculatePenalty(int numPixels, // Number of pixels; for normalisation
-                             psVector *vector, // Vector to which to add in penalty term
+static bool calculatePenalty(psVector *vector, // Vector to which to add in penalty term
                              const pmSubtractionKernels *kernels // Kernel parameters
     )
@@ -393 +392 @@
     }
 
-    float penalty = numPixels * kernels->penalty; // Penalty value
     psVector *penalties = kernels->penalties; // Penalties for each kernel component
     int spatialOrder = kernels->spatialOrder; // Order of spatial variations
@@ -400 +398 @@
         for (int yOrder = 0, index = i; yOrder <= spatialOrder; yOrder++) {
             for (int xOrder = 0; xOrder <= spatialOrder - yOrder; xOrder++, index += numKernels) {
-                vector->data.F64[index] -= penalty * penalties->data.F32[i];
+                vector->data.F64[index] -= penalties->data.F32[i];
             }
         }
@@ -711 +709 @@
             }
         }
-        calculatePenalty(numStamps * PS_SQR(2 * stamps->footprint + 1), sumVector, kernels);
+        calculatePenalty(sumVector, kernels);
 
         psVector *permutation = NULL;       // Permutation vector, required for LU decomposition
@@ -758 +756 @@
             }
         }
-        calculatePenalty(numStamps * PS_SQR(2 * stamps->footprint + 1), sumVector1, kernels);
-        calculatePenalty(numStamps * PS_SQR(2 * stamps->footprint + 1), sumVector2, kernels);
+        calculatePenalty(sumVector1, kernels);
+        calculatePenalty(sumVector2, kernels);
 
 #if 0
@@ -812 +810 @@
         psImage *F = (psImage*)psBinaryOp(NULL, CtBiC, "-", A);
         assert(F->numRows == numParams && F->numCols == numParams);
-        float det = NAN;
+        float det = 0.0;
         psImage *Fi = psMatrixInvert(NULL, F, &det);
         assert(Fi->numRows == numParams && Fi->numCols == numParams);