Changeset 24844 for trunk/psModules/src/imcombine/pmSubtractionEquation.c

Timestamp:

Jul 17, 2009, 2:38:57 PM (17 years ago)

Author:

Paul Price

Message:

Correct deviations for normalisation if convolving image 1. This is required to make deviations from convolving different images comparable.

File:

• trunk/psModules/src/imcombine/pmSubtractionEquation.c (modified) (5 diffs)

trunk/psModules/src/imcombine/pmSubtractionEquation.c

@@ -793 +793 @@
 
 #ifdef TESTING
-              // XXX double-check for NAN in data:
-                for (int iy = 0; iy < stamp->matrix1->numRows; iy++) {
-                    for (int ix = 0; ix < stamp->matrix1->numCols; ix++) {
-                        if (!isfinite(stamp->matrix1->data.F64[iy][ix])) {
-                            fprintf (stderr, "WARNING: NAN in matrix1\n");
-                        }
-                    }
-                }
-                for (int ix = 0; ix < stamp->vector1->n; ix++) {
-                    if (!isfinite(stamp->vector1->data.F64[ix])) {
-                        fprintf (stderr, "WARNING: NAN in vector1\n");
-                    }
-                }
+              // XXX double-check for NAN in data:
+                for (int iy = 0; iy < stamp->matrix1->numRows; iy++) {
+                    for (int ix = 0; ix < stamp->matrix1->numCols; ix++) {
+                        if (!isfinite(stamp->matrix1->data.F64[iy][ix])) {
+                            fprintf (stderr, "WARNING: NAN in matrix1\n");
+                        }
+                    }
+                }
+                for (int ix = 0; ix < stamp->vector1->n; ix++) {
+                    if (!isfinite(stamp->vector1->data.F64[ix])) {
+                        fprintf (stderr, "WARNING: NAN in vector1\n");
+                    }
+                }
 #endif
 
@@ -818 +818 @@
 
 #ifdef TESTING
-        for (int ix = 0; ix < sumVector->n; ix++) {
-            if (!isfinite(sumVector->data.F64[ix])) {
-                fprintf (stderr, "WARNING: NAN in vector1\n");
-            }
-        }
+        for (int ix = 0; ix < sumVector->n; ix++) {
+            if (!isfinite(sumVector->data.F64[ix])) {
+                fprintf (stderr, "WARNING: NAN in vector1\n");
+            }
+        }
 #endif
 
@@ -828 +828 @@
 
 #ifdef TESTING
-        for (int ix = 0; ix < sumVector->n; ix++) {
-            if (!isfinite(sumVector->data.F64[ix])) {
-                fprintf (stderr, "WARNING: NAN in vector1\n");
-            }
-        }
+        for (int ix = 0; ix < sumVector->n; ix++) {
+            if (!isfinite(sumVector->data.F64[ix])) {
+                fprintf (stderr, "WARNING: NAN in vector1\n");
+            }
+        }
         {
             psImage *inverse = psMatrixInvert(NULL, sumMatrix, NULL);
@@ -866 +866 @@
 
 #ifdef TESTING
-        // XXX double-check for NAN in data:
-        for (int ix = 0; ix < kernels->solution1->n; ix++) {
-            if (!isfinite(kernels->solution1->data.F64[ix])) {
-                fprintf (stderr, "WARNING: NAN in vector1\n");
-            }
-        }
+        // XXX double-check for NAN in data:
+        for (int ix = 0; ix < kernels->solution1->n; ix++) {
+            if (!isfinite(kernels->solution1->data.F64[ix])) {
+                fprintf (stderr, "WARNING: NAN in vector1\n");
+            }
+        }
 #endif
 
@@ -1127 +1127 @@
                 source = stamp->image1;
                 convolutions = stamp->convolutions1;
+
+                // Having convolved image1 and changed its normalisation, we need to renormalise the residual
+                // so that it is on the scale of image1.
+                psImage *image = pmSubtractionKernelImage(kernels, stamp->xNorm, stamp->yNorm,
+                                                          false); // Kernel image
+                if (!image) {
+                    psError(PS_ERR_UNKNOWN, false, "Unable to generate image of kernel.");
+                    return false;
+                }
+                double sumKernel = 0;   // Sum of kernel, for normalising residual
+                int size = kernels->size; // Half-size of kernel
+                int fullSize = 2 * size + 1; // Full size of kernel
+                for (int y = 0; y < fullSize; y++) {
+                    for (int x = 0; x < fullSize; x++) {
+                        sumKernel += image->data.F32[y][x];
+                    }
+                }
+                psFree(image);
+                devNorm = 1.0 / sumKernel / numPixels;
                 break;
               case PM_SUBTRACTION_MODE_2: