Changeset 16479 for trunk/psModules/src/imcombine/pmSubtraction.c

Timestamp:

Feb 14, 2008, 1:33:09 PM (18 years ago)

Author:

Paul Price

Message:

Rewinding state of psModules/src/imcombine to the last merge (pap_merge_080122), so that ppStack on the mainline will build (I had branched ppStack but made changes to psModules/src/imcombine on the mainline; this check in is restoring to the state of the mainline while the development proceeds on branch pap_branch_080214)

File:

• trunk/psModules/src/imcombine/pmSubtraction.c (modified) (6 diffs)

trunk/psModules/src/imcombine/pmSubtraction.c

@@ -4 +4 @@
  *  @author GLG, MHPCC
  *
- *  @version $Revision: 1.79 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2008-02-14 23:18:34 $
+ *  @version $Revision: 1.80 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2008-02-14 23:33:09 $
  *
  *  Copyright 2004-2007 Institute for Astronomy, University of Hawaii
@@ -50 +50 @@
 
     // Take the square of the normal kernel
+    double sumNormal = 0.0, sumVariance = 0.0; // Sum of the normal and variance kernels
     for (int v = yMin; v <= yMax; v++) {
         for (int u = xMin; u <= xMax; u++) {
-            out->kernel[v][u] = PS_SQR(normalKernel->kernel[v][u]);
-        }
-    }
+            float value = normalKernel->kernel[v][u]; // Value of interest
+            float value2 = PS_SQR(value); // Value squared
+            sumNormal += value;
+            sumVariance += value2;
+            out->kernel[v][u] = value2;
+        }
+    }
+
+    // Normalise so that the sum of the variance kernel is the square of the sum of the normal kernel
+    // This is required to keep the relative scaling between the image and the weight map
+    psBinaryOp(out->image, out->image, "*", psScalarAlloc(PS_SQR(sumNormal) / sumVariance, PS_TYPE_F32));
 
     return out;
@@ -681 +690 @@
         PS_ASSERT_IMAGE_TYPE(out1->weight, PS_TYPE_F32, false);
     }
-    if (region && psRegionIsNaN(*region)) {
-        psString string = psRegionToString(*region);
-        psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Input region (%s) contains NAN values", string);
-        psFree(string);
-        return false;
+    if (region) {
+        if (psRegionIsNaN(*region)) {
+            psString string = psRegionToString(*region);
+            psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Input region (%s) contains NAN values", string);
+            psFree(string);
+            return false;
+        }
+        if (region->x0 < 0 || region->x1 > ro1->image->numCols ||
+            region->y0 < 0 || region->y1 > ro1->image->numRows) {
+            psString string = psRegionToString(*region);
+            psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Input region (%s) does not fit in image (%dx%d)",
+                    string, ro1->image->numCols, ro1->image->numRows);
+            psFree(string);
+            return false;
+        }
     }
 
@@ -705 +724 @@
 
     int numCols = ro1->image->numCols, numRows = ro1->image->numRows; // Image dimensions
-    int x0 = ro1->col0, y0 = ro1->row0; // Image offset
 
     psImage *convImage1 = out1->image;   // Convolved image
@@ -762 +780 @@
         yMax = PS_MIN(region->y1, yMax);
     }
-
-    // Size to use when calculating normalised coordinates (different from actual size when convolving
-    // subimage)
-    int xNormSize = (kernels->numCols > 0 ? kernels->numCols : numCols);
-    int yNormSize = (kernels->numRows > 0 ? kernels->numRows : numRows);
 
     psMaskType maskSource;              // Mask these pixels when convolving
@@ -794 +807 @@
     for (int j = yMin; j < yMax; j += fullSize) {
         int ySubMax = PS_MIN(j + fullSize, yMax); // Range for subregion of interest
-        float yNorm = 2.0 * (float)(j + y0 + size + 1 - yNormSize/2.0) /
-            (float)yNormSize; // Normalised coordinate
+        float yNorm = 2.0 * (float)(j + size + 1 - numRows/2.0) / (float)numRows; // Normalised coordinate
         for (int i = xMin; i < xMax; i += fullSize) {
             int xSubMax = PS_MIN(i + fullSize, xMax); // Range for subregion of interest
-            float xNorm = 2.0 * (float)(i + x0 + size + 1 - xNormSize/2.0) /
-                (float)xNormSize; // Normalised coordinate
+            float xNorm = 2.0 * (float)(i + size + 1 - numCols/2.0) / (float)numCols; // Normalised coordinate
 
             // Only generate polynomial values every kernel footprint, since we have already assumed