Changeset 24951 for branches/czw_branch/cleanup/psModules/src/imcombine/pmSubtractionEquation.c

Timestamp:

Jul 30, 2009, 5:20:29 PM (17 years ago)

Author:

watersc1

Message:

Finished up my edits to the detrend cleanup, and some changes to my
copy of burntool and the pslib astrometry. Detrend cleanup has not
been tested yet. That's up next.

Location:

branches/czw_branch/cleanup

Files:

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

branches/czw_branch/cleanup/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: