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

Timestamp:

Aug 22, 2007, 3:59:31 PM (19 years ago)

Author:

Paul Price

Message:

Adding check on return of pmSubtractionKernelImage.
Fixing coordinates (supposed to be -1..1).
Moving thresholding out of loop (if we need to keep the
entire convolved image, we may as well threshold once).

File:

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

trunk/psModules/src/imcombine/pmSubtraction.c

@@ -4 +4 @@
  *  @author GLG, MHPCC
  *
- *  @version $Revision: 1.44 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2007-08-22 02:56:53 $
+ *  @version $Revision: 1.45 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-08-23 01:59:31 $
  *
  *  Copyright 2004-2007 Institute for Astronomy, University of Hawaii
@@ -1064 +1064 @@
 
     // Convolve the image with the kernel --- we're basically applying a matched filter and then thresholding
-    psImage *convolved = NULL;          // Convolved mask image
     psPixels *out = NULL;               // List of pixels that should be masked
     for (int i = 0; i < numRegions; i++) {
@@ -1079 +1078 @@
         // the kernel may scale the unit level that we've inserted.
 
-        // Image of the kernel at the centre
-        psImage *kernel = pmSubtractionKernelImage(solution, kernels,
-                                                   region->x0 + 0.5 * (region->x1 - region->x0),
-                                                   region->y0 + 0.5 * (region->y1 - region->y0));
+        // Image of the kernel at the centre of the region
+        float xNorm = (region->x0 + 0.5 * (region->x1 - region->x0) - numCols/2.0) / (float)numCols;
+        float yNorm = (region->y0 + 0.5 * (region->y1 - region->y0) - numRows/2.0) / (float)numRows;
+        psImage *kernel = pmSubtractionKernelImage(solution, kernels, xNorm, yNorm);
+        if (!kernel) {
+            psError(PS_ERR_UNKNOWN, false, "Unable to generate kernel image.");
+            psFree(convolved);
+            psFree(image);
+            return NULL;
+        }
         float sum = 0.0;
         for (int y = 0; y < kernel->numRows; y++) {
@@ -1091 +1096 @@
         psFree(kernel);
 
-        // Threshold the convolved image
-        float adjustedThreshold = sum * threshold; // Threshold adjusted for the kernel normalisation
-        for (int y = region->y0; y < region->y1; y++) {
-            for (int x = region->x0; x < region->x1; x++) {
-                if (image->data.F32[y][x] > adjustedThreshold) {
-                    out = psPixelsAdd(out, PIXEL_LIST_BUFFER, x, y);
-                }
-            }
-        }
+        psImage *subConv = psImageSubset(convolved, *region); // Sub-image of convolved image
+        psBinaryOp(subConv, subConv, "*", psScalarAlloc(1.0 / sum, PS_TYPE_F32));
     }
     psFree(image);
 
+    // Threshold the convolved image
+    for (int y = 0; y < numRows; y++) {
+        for (int x = 0; x < numCols; x++) {
+            if (convolved->data.F32[y][x] > threshold) {
+                out = psPixelsAdd(out, PIXEL_LIST_BUFFER, x, y);
+            }
+        }
+    }
+    psFree(convolved);
+
     return out;
 }