Changeset 27840 for branches/simtest_nebulous_branches/psModules/src/imcombine/pmStackReject.c

Timestamp:

May 3, 2010, 8:50:52 AM (16 years ago)

Author:

eugene

Message:

updates from trunk

Location:

branches/simtest_nebulous_branches

Files:

• . (modified) (1 prop)
• psModules (modified) (1 prop)
• psModules/src/imcombine/pmStackReject.c (modified) (15 diffs)

branches/simtest_nebulous_branches/psModules/src/imcombine/pmStackReject.c

@@ -10 +10 @@
 #include "pmSubtractionThreads.h"
 #include "pmSubtractionKernels.h"
+
+#include "pmStackReject.h"
 
 #define PIXEL_LIST_BUFFER 100           // Number of pixels to add to list at a time
@@ -35 +37 @@
 {
     int size = kernels->size;           // Half-size of convolution kernel
-    psImage *polyValues = p_pmSubtractionPolynomialFromCoords(NULL, kernels, numCols, numRows,
-                                                              xMin + size + 1, yMin + size + 1); // Polynomial
+    int x = PS_MIN(xMin + size + 1, kernels->xMax); // x coordinate of interest
+    int y = PS_MIN(yMin + size + 1, kernels->yMax); // y coordinate of interest
+
+    psImage *polyValues = p_pmSubtractionPolynomialFromCoords(NULL, kernels, x, y); // Polynomial
     int box = p_pmSubtractionBadRadius(NULL, kernels, polyValues, false, poorFrac); // Radius of bad box
     psTrace("psModules.imcombine", 10, "Growing by %d", box);
@@ -43 +47 @@
     if (box > 0) {
         // Convolve a subimage, then stick it in the target
-        // XXX if (threaded) {
-        // XXX     psMutexLock(source);
-        // XXX }
         psImage *mask = psImageSubset(source, psRegionSet(xMin - box, xMax + box,
                                                           yMin - box, yMax + box)); // Mask to convolve
-        // XXX if (threaded) {
-        // XXX     psMutexUnlock(source);
-        // XXX }
         psImage *convolved = psImageConvolveMask(NULL, mask, PM_STACK_MASK_BAD, PM_STACK_MASK_CONVOLVE,
                                                  -box, box, -box, box); // Convolved mask
-        // XXX if (threaded) {
-        // XXX     psMutexLock(source);
-        // XXX }
         psFree(mask);
-        // XXX if (threaded) {
-        // XXX     psMutexUnlock(source);
-        // XXX }
 
         int numBytes = (xMax - xMin) * PSELEMTYPE_SIZEOF(PS_TYPE_IMAGE_MASK); // Number of bytes to copy
@@ -111 +103 @@
 
     if (!pmSubtractionThreaded()) {
-        pmSubtractionThreadsInit(NULL, NULL);
+        pmSubtractionThreadsInit();
     }
 
@@ -125 +117 @@
 
 
-psPixels *pmStackReject(const psPixels *in, int numCols, int numRows, float threshold, float poorFrac,
-                        int stride, const psArray *subRegions, const psArray *subKernels)
+psPixels *pmStackReject(const psPixels *in, int numCols, int numRows, float threshold, int stride,
+                        const psArray *subRegions, const psArray *subKernels)
 {
     PS_ASSERT_PIXELS_NON_NULL(in, NULL);
@@ -162 +154 @@
     pmReadout *inRO = pmReadoutAlloc(NULL); // Readout with input image
     inRO->image = image;
-    // XXX if (threaded) {
-    // XXX     psMutexInit(image);
-    // XXX }
+    convRO->image = psImageAlloc(image->numCols, image->numRows, PS_TYPE_F32);
     for (int i = 0; i < numRegions; i++) {
         psRegion *region = subRegions->data[i]; // Region of interest
         pmSubtractionKernels *kernels = subKernels->data[i]; // Kernel of interest
-        if (!pmSubtractionConvolve(convRO, NULL, inRO, NULL, NULL, stride, 0, 0, 1.0, 0.0,
+        if (!pmSubtractionConvolve(NULL, convRO, NULL, inRO, NULL, stride, 0, 0, 1.0, 0.0, 0.0,
                                    region, kernels, false, true)) {
-            psError(PS_ERR_UNKNOWN, false, "Unable to convolve mask image in region %d.", i);
+            psError(psErrorCodeLast(), false, "Unable to convolve mask image in region %d.", i);
             psFree(convRO);
             psFree(inRO);
@@ -180 +170 @@
 
         // Image of the kernel at the centre of the region
-        float xNorm = (region->x0 + 0.5 * (region->x1 - region->x0) - kernels->numCols/2.0) /
-            (float)kernels->numCols;
-        float yNorm = (region->y0 + 0.5 * (region->y1 - region->y0) - kernels->numRows/2.0) /
-            (float)kernels->numRows;
-        psImage *kernel = pmSubtractionKernelImage(kernels, xNorm, yNorm, false);
+        psImage *kernel = pmSubtractionKernelImage(kernels, 0.5, 0.5, false);
         if (!kernel) {
-            psError(PS_ERR_UNKNOWN, false, "Unable to generate kernel image.");
+            psError(psErrorCodeLast(), false, "Unable to generate kernel image.");
             psFree(convRO);
             psFree(inRO);
@@ -211 +197 @@
         }
     }
-    // XXX if (threaded) {
-    // XXX     psMutexDestroy(image);
-    // XXX }
     psFree(inRO);
     psImage *convolved = psMemIncrRefCounter(convRO->image);
@@ -241 +224 @@
     }
     psTrace("psModules.imcombine", 7, "Found %ld bad pixels", bad->n);
-
-    // Now, grow the mask to include everything that touches a bad pixel in the convolution
-    psImage *source = psPixelsToMask(NULL, bad, psRegionSet(0, numCols - 1, 0, numRows - 1),
+    psFree(convolved);
+
+    return bad;
+}
+
+
+psPixels *pmStackRejectGrow(const psPixels *in, int numCols, int numRows, float poorFrac,
+                            const psArray *subRegions, const psArray *subKernels)
+{
+    PS_ASSERT_PIXELS_NON_NULL(in, NULL);
+    PS_ASSERT_ARRAY_NON_NULL(subRegions, NULL);
+    PS_ASSERT_ARRAY_NON_NULL(subKernels, NULL);
+    PS_ASSERT_ARRAYS_SIZE_EQUAL(subRegions, subKernels, NULL);
+
+    psImage *source = psPixelsToMask(NULL, in, psRegionSet(0, numCols - 1, 0, numRows - 1),
                                      PM_STACK_MASK_BAD); // Mask image to grow
 
@@ -262 +257 @@
     bool oldThreads = psImageConvolveSetThreads(false); // Old value of threading for psImageColvolve
 
-    psImage *target = psImageRecycle(convolved, numCols, numRows, PS_TYPE_IMAGE_MASK); // Grown image
+    psImage *target = psImageAlloc(numCols, numRows, PS_TYPE_IMAGE_MASK); // Grown image
     psImageInit(target, 0);
-    // XXX if (threaded) {
-    // XXX     psMutexInit(source);
-    // XXX }
     for (int i = 0; i < subRegions->n; i++) {
         psRegion *region = subRegions->data[i]; // Subtraction region
@@ -287 +279 @@
                     psArray *args = job->args; // Job arguments
                     psArrayAdd(args, 1, target);
-                    // XXX psMutexLock(source);
                     psArrayAdd(args, 1, source);
-                    // XXX psMutexUnlock(source);
                     psArrayAdd(args, 1, kernels);
                     PS_ARRAY_ADD_SCALAR(args, numCols, PS_TYPE_S32);
@@ -307 +297 @@
                 } else if (!stackRejectGrow(target, source, kernels, numCols, numRows,
                                             i, xSubMax, j, ySubMax, poorFrac)) {
-                    psError(PS_ERR_UNKNOWN, false, "Unable to grow bad pixels.");
+                    psError(psErrorCodeLast(), false, "Unable to grow bad pixels.");
                     psFree(source);
                     psFree(target);
@@ -317 +307 @@
 
     if (!psThreadPoolWait(false)) {
-        psError(PS_ERR_UNKNOWN, false, "Unable to grow bad pixels.");
+        psError(psErrorCodeLast(), false, "Unable to grow bad pixels.");
         psFree(source);
         psFree(target);
@@ -332 +322 @@
         }
 
-        // XXX psMutexDestroy(source);
     }
 
@@ -351 +340 @@
 
     psFree(source);
-    bad = psPixelsFromMask(bad, target, PM_STACK_MASK_ALL);
+    psPixels *bad = psPixelsFromMask(NULL, target, PM_STACK_MASK_ALL); // All bad pixels
+    psFree(target);
     psTrace("psModules.imcombine", 7, "Total %ld bad pixels", bad->n);