Changeset 16423 for trunk/psModules/src/imcombine/pmStackReject.c

Timestamp:

Feb 13, 2008, 2:13:09 PM (18 years ago)

Author:

Paul Price

Message:

Starting to attempt to support stacking of images that have been read in piecemeal (scan by scan). Added parameters to pmSubtractionKernels to specify size of the (original) image, so that normalisation can be done (during convolution) when the image size is not the original image size. Doing some algebraic gymnastics in the pixel rejection to account for the image being smaller than the original.

File:

• trunk/psModules/src/imcombine/pmStackReject.c (modified) (7 diffs)

trunk/psModules/src/imcombine/pmStackReject.c

@@ -11 +11 @@
 #define PIXEL_LIST_BUFFER 100           // Number of pixels to add to list at a time
 
-psPixels *pmStackReject(const psPixels *in, float threshold, const psArray *regions, const psArray *kernels)
+psPixels *pmStackReject(const psPixels *in, const psRegion *valid, float threshold,
+                        const psArray *subRegions, const psArray *kernels)
 {
     PS_ASSERT_PIXELS_NON_NULL(in, NULL);
     PS_ASSERT_FLOAT_LARGER_THAN_OR_EQUAL(threshold, 0.0, NULL);
     PS_ASSERT_FLOAT_LESS_THAN_OR_EQUAL(threshold, 1.0, NULL);
-    PS_ASSERT_ARRAY_NON_NULL(regions, NULL);
+    PS_ASSERT_ARRAY_NON_NULL(subRegions, NULL);
     PS_ASSERT_ARRAY_NON_NULL(kernels, NULL);
-    PS_ASSERT_ARRAYS_SIZE_EQUAL(regions, kernels, NULL);
+    PS_ASSERT_ARRAYS_SIZE_EQUAL(subRegions, kernels, NULL);
 
     // Get the original image size
-    int numRegions = regions->n;        // Number of regions
+    int numRegions = subRegions->n;        // Number of regions
     int numCols = 0, numRows = 0;       // Size of original image
     int minCols = INT_MAX, minRows = INT_MAX; // Minimum coordinate for image --- should be 0,0
     int size = 0;                       // Size of kernel
     for (int i = 0; i < numRegions; i++) {
-        psRegion *region = regions->data[i]; // Region of interest
-        if (region->x0 < minCols) {
-            minCols = region->x0;
+        psRegion *subRegion = subRegions->data[i]; // Region of interest
+        if (subRegion->x0 < minCols) {
+            minCols = subRegion->x0;
         }
-        if (region->y0 < minRows) {
-            minRows = region->y0;
+        if (subRegion->y0 < minRows) {
+            minRows = subRegion->y0;
         }
-        if (region->x1 > numCols) {
-            numCols = region->x1;
+        if (subRegion->x1 > numCols) {
+            numCols = subRegion->x1;
         }
-        if (region->y1 > numRows) {
-            numRows = region->y1;
+        if (subRegion->y1 > numRows) {
+            numRows = subRegion->y1;
         }
 
@@ -49 +50 @@
         }
     }
-    if (minCols != 0 || minRows != 0) {
-        psError(PS_ERR_BAD_PARAMETER_VALUE, true,
-                "Some error with image regions --- minimum coordinate is not 0,0");
-        return NULL;
+
+    // Adjust the size for the size of the subimage
+    if (valid) {
+        minCols = PS_MAX(valid->x0, minCols);
+        minRows = PS_MAX(valid->y0, minRows);
+        numCols = PS_MIN(valid->x1, numCols);
+        numRows = PS_MIN(valid->y1, numRows);
     }
 
-    psImage *mask = psPixelsToMask(NULL, in, psRegionSet(0, numCols, 0, numRows), 0x01); // Mask image
+    psImage *mask = psPixelsToMask(NULL, in, psRegionSet(minCols, numCols, minRows, numRows), 0x01); // Mask
     psImage *image = psImageCopy(NULL, mask, PS_TYPE_F32); // Floating-point version, so we can convolve
     psFree(mask);
@@ -63 +67 @@
     pmReadout *inRO = pmReadoutAlloc(NULL); // Readout with input image
     inRO->image = image;
+    inRO->col0 = minCols;
+    inRO->row0 = minRows;
     for (int i = 0; i < numRegions; i++) {
-        psRegion *region = regions->data[i]; // Region of interest
+        psRegion *subRegion = subRegions->data[i]; // Region of interest
+        if (valid && (subRegion->x0 > valid->x1 || subRegion->x1 < valid->x0 ||
+                      subRegion->y0 > valid->y1 || subRegion->y1 < valid->y0)) {
+            // Region is outside of our sub-image
+            continue;
+        }
         pmSubtractionKernels *kernel = kernels->data[i]; // Kernel of interest
-        if (!pmSubtractionConvolve(convRO, NULL, inRO, NULL, NULL, 0, region, kernel, true)) {
+        if (!pmSubtractionConvolve(convRO, NULL, inRO, NULL, NULL, 0, subRegion, kernel, true)) {
             psError(PS_ERR_UNKNOWN, false, "Unable to convolve mask image in region %d.", i);
             psFree(convRO);
@@ -77 +88 @@
 
         // 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;
+        float xNorm = (subRegion->x0 + 0.5 * (subRegion->x1 - subRegion->x0) - kernel->numCols/2.0) /
+            (float)kernel->numCols;
+        float yNorm = (subRegion->y0 + 0.5 * (subRegion->y1 - subRegion->y0) - kernel->numRows/2.0) /
+            (float)kernel->numRows;
         psImage *image = pmSubtractionKernelImage(kernel, xNorm, yNorm, false);
         if (!image) {
@@ -94 +107 @@
         psFree(image);
 
-        psImage *subConv = psImageSubset(convRO->image, *region); // Sub-image of convolved image
+        psImage *subConv = psImageSubset(convRO->image, *subRegion); // Sub-image of convolved image
         psBinaryOp(subConv, subConv, "*", psScalarAlloc(1.0 / sum, PS_TYPE_F32));
     }
@@ -114 +127 @@
     // Now, grow the mask to include everything that touches a bad pixel in the convolution
     mask = psPixelsToMask(NULL, bad, psRegionSet(0, numCols - 1, 0, numRows - 1), 0xff);
+    assert(mask->numCols == numCols && mask->numRows == numRows);
     for (int i = 0; i < bad->n; i++) {
         int xPix = bad->data[i].x, yPix = bad->data[i].y; // Coordinates of interest
@@ -123 +137 @@
         for (int y = yMin; y <= yMax; y++) {
             for (int x = xMin; x <= xMax; x++) {
+                assert(x < mask->numCols && y < mask->numRows);
                 mask->data.PS_TYPE_MASK_DATA[y][x] = 0xff;
             }