Changeset 16604

Timestamp:

Feb 22, 2008, 9:24:42 AM (18 years ago)

Author:

Paul Price

Message:

Merging in ppStack development branch --- ppStack now works with incremental reads. Small conflicts (due to stupid patch application by CVS) resolved.

Location:

trunk/psModules/src/imcombine

Files:

• pmReadoutCombine.c (modified) (1 diff)
• pmStack.c (modified) (12 diffs)
• pmStackReject.c (modified) (6 diffs)
• pmSubtraction.c (modified) (4 diffs)
• pmSubtraction.h (modified) (2 diffs)
• pmSubtractionMatch.c (modified) (1 diff)

trunk/psModules/src/imcombine/pmReadoutCombine.c

@@ -113 +113 @@
     }
 
-    pmReadoutUpdateSize(output, minInputCols, minInputRows, xSize, ySize, true);
+    pmReadoutUpdateSize(output, minInputCols, minInputRows, xSize, ySize, true, params->weights,
+                        params->blank);
     psTrace("psModules.imcombine", 7, "Output minimum: %d,%d\n", output->col0, output->row0);
 
     psStatsOptions combineStdev = 0; // Statistics option for weights
     if (params->weights) {
-
-        if (!output->weight) {
-            output->weight = psImageAlloc(xSize, ySize, PS_TYPE_F32);
-        }
-        if (output->weight->numCols < xSize || output->weight->numRows < ySize) {
-            psImage *newWeight = psImageAlloc(xSize, ySize, PS_TYPE_F32);
-            psImageInit(newWeight, 0.0);
-            psImageOverlaySection(newWeight, output->weight, output->col0, output->row0, "=");
-            psFree(output->weight);
-            output->weight = newWeight;
-        }
-
         if (first) {
             psMetadataAddStr(hdu->header, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK,

trunk/psModules/src/imcombine/pmStack.c

@@ -8 +8 @@
  *  @author GLG, MHPCC
  *
- *  @version $Revision: 1.17 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2008-02-14 23:33:09 $
  *  Copyright 2004-2007 Institute for Astronomy, University of Hawaii
  *
@@ -23 +21 @@
 #include "pmHDU.h"
 #include "pmFPA.h"
+#include "pmReadoutStack.h"
 #include "pmConceptsAverage.h"
 
@@ -42 +41 @@
     psVector *weights;                  // Pixel weights
     psVector *sort;                     // Buffer for sorting (to get a robust estimator of the standard dev)
+#if 0
+    int x0, y0;                         // Offset from original image to combination region
+    int nx, ny;                         // Number of pixels to combine
+#endif
 } combineBuffer;
 
@@ -64 +67 @@
     buffer->weights = psVectorAlloc(numImages, PS_TYPE_F32);
     buffer->sort = psVectorAlloc(numImages, PS_TYPE_F32);
+
+#if 0
+    buffer->x0 = 0;
+    buffer->y0 = 0;
+    buffer->nx = 0;
+    buffer->ny = 0;
+#endif
 
     return buffer;
@@ -136 +146 @@
         return false;
     }
-    *median = num % 2 ? (sortBuffer->data.F32[num / 2] + sortBuffer->data.F32[num / 2 + 1]) / 2.0 :
-        sortBuffer->data.F32[num / 2];
-    *stdev = 0.74 * (sortBuffer->data.F32[(int)(0.75 * num)] - sortBuffer->data.F32[(int)(0.25 * num)]);
+    *median = num % 2 ? sortBuffer->data.F32[num / 2] :
+        (sortBuffer->data.F32[num / 2] + sortBuffer->data.F32[num / 2 + 1]) / 2.0 ;
+#if 0
+    if (num < NUM_DIRECT_STDEV) {
+#endif
+        // If there are not many values, the direct standard deviation is better
+        double sum = 0.0;
+        for (int i = 0; i < num; i++) {
+            sum += PS_SQR(sortBuffer->data.F32[i] - *median);
+        }
+        *stdev = sqrt(sum / (float)(num - 1));
+#if 0
+    } else {
+        *stdev = 0.74 * (sortBuffer->data.F32[(int)(0.75 * num)] - sortBuffer->data.F32[(int)(0.25 * num)]);
+    }
+#endif
 
     return true;
@@ -186 +209 @@
 
 // Given a stack of images, combine with optional rejection.
-// Pixels in the stack that are rejected are marked for subsequent
+// Pixels in the stack that are rejected are marked for subsequent inspection
 static bool combinePixels(psImage *image, // Combined image, for output
                           psImage *mask, // Combined mask, for output
@@ -193 +216 @@
                           const psVector *weights, // Global (single value) weights for data, or NULL
                           const psVector *reject, // Indices of pixels to reject, or NULL
-                          int x, int y, // Coordinates of interest
+                          int x, int y, // Coordinates of interest; frame of output image
                           psMaskType maskVal, // Value to mask
                           psMaskType bad, // Value to give bad pixels
@@ -223 +246 @@
             continue;
         }
+        int xIn = x - data->readout->col0, yIn = y - data->readout->row0; // Coordinates on input readout
         psImage *image = data->readout->image; // Image of interest
         psImage *weight = data->readout->weight; // Weight map of interest
         psImage *mask = data->readout->mask; // Mask of interest
-        pixelData->data.F32[i] = image->data.F32[y][x];
+        pixelData->data.F32[i] = image->data.F32[yIn][xIn];
         if (weight) {
-            pixelWeights->data.F32[i] = weight->data.F32[y][x];
-        }
-        pixelMasks->data.PS_TYPE_MASK_DATA[i] = mask->data.PS_TYPE_MASK_DATA[y][x];
+            pixelWeights->data.F32[i] = weight->data.F32[yIn][xIn];
+        }
+        pixelMasks->data.PS_TYPE_MASK_DATA[i] = mask->data.PS_TYPE_MASK_DATA[yIn][xIn];
         if (pixelMasks->data.PS_TYPE_MASK_DATA[i] & maskVal) {
             numBad++;
@@ -563 +587 @@
     }
 
+    // Get the sizes
+    psArray *stack = psArrayAlloc(input->n);
+    for (int i = 0; i < stack->n; i++) {
+        pmStackData *data = input->data[i]; // Stack data
+        stack->data[i] = psMemIncrRefCounter(data->readout);
+    }
+    int minInputCols, maxInputCols, minInputRows, maxInputRows; // Smallest and largest values to combine
+    int xSize, ySize;                   // Size of the output image
+    if (!pmReadoutStackValidate(&minInputCols, &maxInputCols, &minInputRows, &maxInputRows, &xSize, &ySize,
+                                stack)) {
+        psError(PS_ERR_UNKNOWN, false, "Input stack is not valid.");
+        psFree(stack);
+        return false;
+    }
+    psFree(stack);
+    pmReadoutUpdateSize(combined, minInputCols, minInputRows, xSize, ySize, true, true, bad);
+    psTrace("psModules.imcombine", 1, "Combining [%d:%d,%d:%d] (%dx%d)\n",
+            minInputCols, maxInputCols, minInputRows, maxInputRows, xSize, ySize);
+
+
     // Buffer for combination
     combineBuffer *buffer = combineBufferAlloc(num, numIter == 0 ? PS_STAT_SAMPLE_MEAN :
@@ -573 +617 @@
         psImage *combinedWeight = combined->weight; // Combined mask
 
-        psArray *pixelMap = pixelMapGenerate(input, numCols, numRows); // Map of pixels to source
+        psArray *pixelMap = pixelMapGenerate(input, maxInputCols, maxInputRows); // Map of pixels to source
         psPixels *pixels = NULL;            // Total list of pixels, with no duplicates
         for (int i = 0; i < num; i++) {
@@ -584 +628 @@
         }
         for (int i = 0; i < pixels->n; i++) {
+            // Pixel coordinates are in the frame of the original image
             int x = pixels->data[i].x, y = pixels->data[i].y; // Coordinates of interest
+            if (x < minInputCols || x >= maxInputCols || y < minInputRows || y >= maxInputRows) {
+                continue;
+            }
             psVector *reject = pixelMapQuery(pixelMap, x, y); // Inspect these images closely
             combinePixels(combinedImage, combinedMask, combinedWeight, input, weights, reject, x, y,
@@ -622 +670 @@
         }
 
-        for (int y = 0; y < numRows; y++) {
-            for (int x = 0; x < numCols; x++) {
+        for (int y = minInputRows; y < maxInputRows; y++) {
+            for (int x = minInputCols; x < maxInputCols; x++) {
                 combinePixels(combinedImage, combinedMask, combinedWeights, input, weights, NULL, x, y,
                               maskVal, bad, numIter, rej, buffer);

trunk/psModules/src/imcombine/pmStackReject.c

@@ -25 +25 @@
     int numRegions = regions->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 minCols = INT_MAX, minRows = INT_MAX; // Minimum coordinate for image
+    int size = 0;                       // Size of kernel
     for (int i = 0; i < numRegions; i++) {
         psRegion *region = regions->data[i]; // Region of interest
@@ -46 +47 @@
         return NULL;
     }
+    psTrace("psModules.imcombine", 1, "Rejecting [%d:%d,%d:%d]\n", minCols, numCols, minRows, numRows);
 
     psImage *mask = psPixelsToMask(NULL, in, psRegionSet(0, numCols, 0, numRows), 0x01); // Mask image
@@ -56 +58 @@
     inRO->image = image;
     for (int i = 0; i < numRegions; i++) {
-        psRegion *region = regions->data[i]; // Region of interest
-        if (!pmSubtractionConvolve(convRO, NULL, inRO, NULL, NULL, 0, region, kernels, true)) {
+        psRegion *region = subRegions->data[i]; // Region of interest
+        if (valid && (region->x0 > valid->x1 || region->x1 < valid->x0 ||
+                      region->y0 > valid->y1 || region->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, false, true)) {
             psError(PS_ERR_UNKNOWN, false, "Unable to convolve mask image in region %d.", i);
             psFree(convRO);
@@ -94 +102 @@
     // Threshold the convolved image
     psPixels *bad = psPixelsAllocEmpty(PIXEL_LIST_BUFFER); // List of pixels that should be masked
-    for (int y = 0; y < numRows; y++) {
-        for (int x = 0; x < numCols; x++) {
+    for (int y = size; y < convolved->numRows - size; y++) {
+        for (int x = size; x < convolved->numCols - size; x++) {
             if (convolved->data.F32[y][x] > threshold) {
-                bad = psPixelsAdd(bad, PIXEL_LIST_BUFFER, x, y);
+                // Pixel coordinates in "bad" correspond to the full image
+                bad = psPixelsAdd(bad, PIXEL_LIST_BUFFER, x + minCols, y + minRows);
             }
         }
@@ -103 +112 @@
     psFree(convolved);
 
-    // Now, we want to convolve the original pixels properly
-    mask = psPixelsToMask(NULL, bad, psRegionSet(0, numCols, 0, numRows), 0xff);
-    int size = kernels->size;           // Size of kernels
+    // Now, grow the mask to include everything that touches a bad pixel in the convolution
+    int x0 = minCols, y0 = minRows;     // Offset for mask image
+    mask = psPixelsToMask(NULL, bad, psRegionSet(x0, numCols - 1, y0, numRows - 1), 0xff);
     for (int i = 0; i < bad->n; i++) {
-        int xPix = bad->data[i].x, yPix = bad->data[i].y; // Coordinates of interest
+        int xPix = bad->data[i].x - x0, yPix = bad->data[i].y - y0; // Coordinates in frame of mask image
         // Convolution limits
         int xMin = PS_MAX(xPix - size, 0);
-        int xMax = PS_MIN(xPix + size, numCols - 1);
+        int xMax = PS_MIN(xPix + size, mask->numCols - 1);
         int yMin = PS_MAX(yPix - size, 0);
-        int yMax = PS_MIN(yPix + size, numRows - 1);
+        int yMax = PS_MIN(yPix + size, mask->numRows - 1);
         for (int y = yMin; y <= yMax; y++) {
             for (int x = xMin; x <= xMax; x++) {
@@ -119 +128 @@
         }
     }
-
     bad = psPixelsFromMask(bad, mask, 0xff);
     psFree(mask);
 
+    // Convert coordinates to frame of original image
+    for (int i = 0; i < bad->n; i++) {
+        int x = bad->data[i].x + x0;
+        int y = bad->data[i].y + y0;
+        if (x < 0 || x >= numCols || y < 0 || y >= numRows) {
+            psWarning("Bad pixel coordinate %d: %d,%d --- ignored.",
+                      i, x, y);
+            continue;
+        }
+        bad->data[i].x = x;
+        bad->data[i].y = y;
+    }
+
     return bad;
 }

trunk/psModules/src/imcombine/pmSubtraction.c

@@ -3 +3 @@
  *  @author Paul Price, IfA
  *  @author GLG, MHPCC
- *
- *  @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
@@ -645 +642 @@
 bool pmSubtractionConvolve(pmReadout *out1, pmReadout *out2, const pmReadout *ro1, const pmReadout *ro2,
                            const psImage *subMask, psMaskType blank, const psRegion *region,
-                           const pmSubtractionKernels *kernels, bool useFFT)
+                           const pmSubtractionKernels *kernels, bool doBG, bool useFFT)
 {
     PS_ASSERT_PTR_NON_NULL(out1, false);
@@ -707 +704 @@
     }
 
-    psImage *image, *weight;            // Image and weight map to convolve
+    const pmReadout *source;            // Source for image parameters
     switch (kernels->mode) {
       case PM_SUBTRACTION_MODE_TARGET:
       case PM_SUBTRACTION_MODE_1:
       case PM_SUBTRACTION_MODE_DUAL:
-        image = ro1->image;
-        weight = ro1->weight;
+        source = ro1;
         break;
       case PM_SUBTRACTION_MODE_2:
-        image = ro2->image;
-        weight = ro2->weight;
+        source = ro2;
         break;
       default:
         psAbort("Unsupported subtraction mode: %x", kernels->mode);
     }
-
-    int numCols = ro1->image->numCols, numRows = ro1->image->numRows; // Image dimensions
+    psImage *image = source->image, *weight = source->weight; // Image and weight map to convolve
+    int numCols = image->numCols, numRows = image->numRows; // Image dimensions
+    int x0 = source->col0, y0 = source->row0; // Image offset
 
     psImage *convImage1 = out1->image;   // Convolved image
@@ -815 +811 @@
             // (with the stamps) that it does not vary rapidly on this scale.
             polyValues = p_pmSubtractionPolynomial(polyValues, kernels->spatialOrder, xNorm, yNorm);
-            float background = p_pmSubtractionSolutionBackground(kernels, polyValues); // Background term
-
+            float background = doBG ? p_pmSubtractionSolutionBackground(kernels, polyValues) :
+                0.0; // Background term
             psRegion subRegion = psRegionSet(i, xSubMax, j, ySubMax); // Sub-region to convolve
 

trunk/psModules/src/imcombine/pmSubtraction.h

@@ -6 +6 @@
  * @author GLG, MHPCC
  *
- * @version $Revision: 1.22 $ $Name: not supported by cvs2svn $
- * @date $Date: 2007-12-07 01:57:15 $
+ * @version $Revision: 1.23 $ $Name: not supported by cvs2svn $
+ * @date $Date: 2008-02-22 19:24:42 $
  * Copyright 2004-207 Institute for Astronomy, University of Hawaii
  */
@@ -92 +92 @@
                            const psRegion *region, ///< Region to convolve (or NULL)
                            const pmSubtractionKernels *kernels, ///< Kernel parameters
+                           bool doBG,   ///< Apply background term?
                            bool useFFT  ///< Use Fast Fourier Transform for the convolution?
     );

trunk/psModules/src/imcombine/pmSubtractionMatch.c

@@ -434 +434 @@
 
             psTrace("psModules.imcombine", 2, "Convolving...\n");
-            if (!pmSubtractionConvolve(conv1, conv2, ro1, ro2, subMask, maskBlank, region, kernels, useFFT)) {
+            if (!pmSubtractionConvolve(conv1, conv2, ro1, ro2, subMask, maskBlank, region, kernels,
+                                       true, useFFT)) {
                 psError(PS_ERR_UNKNOWN, false, "Unable to convolve image.");
                 goto MATCH_ERROR;