Changeset 17005

Timestamp:

Mar 17, 2008, 11:38:43 AM (18 years ago)

Author:

Paul Price

Message:

Allow an entire image to be rejected, in which case the entire combined image needs to be recalculated. Also split the pmStackData->pixels into separate inspect and reject parts, to minimise confusion.

Location:

trunk/psModules/src/imcombine

Files:

• pmStack.c (modified) (17 diffs)
• pmStack.h (modified) (3 diffs)

trunk/psModules/src/imcombine/pmStack.c

@@ -8 +8 @@
  *  @author GLG, MHPCC
  *
- *  @version $Revision: 1.29 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2008-03-06 21:38:39 $
+ *  @version $Revision: 1.30 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2008-03-17 21:38:43 $
  *  Copyright 2004-2007 Institute for Astronomy, University of Hawaii
  *
@@ -80 +80 @@
 {
     psFree(data->readout);
-    psFree(data->pixels);
+    psFree(data->reject);
+    psFree(data->inspect);
     return;
 }
@@ -172 +173 @@
     } else {
         *median = num % 2 ? sortBuffer->data.F32[num / 2] :
-            (sortBuffer->data.F32[num / 2] + sortBuffer->data.F32[num / 2 + 1]) / 2.0;
+            (sortBuffer->data.F32[num / 2 - 1] + sortBuffer->data.F32[num / 2]) / 2.0;
         if (stdev) {
             if (num <= NUM_DIRECT_STDEV) {
@@ -180 +181 @@
                     sum += PS_SQR(sortBuffer->data.F32[i] - *median);
                 }
-                *stdev = sqrt(sum / (float)(num - 1));
+                *stdev = sqrt(sum / (double)(num - 1));
             } else {
                 // Standard deviation from the interquartile range
@@ -201 +202 @@
     pmStackData *data = inputs->data[source]; // Stack data of interest
     if (!data) {
+        psWarning("Can't find input data for source %d", source);
         return;
     }
-    data->pixels = psPixelsAdd(data->pixels, PIXEL_LIST_BUFFER, x, y);
+    data->inspect = psPixelsAdd(data->inspect, PIXEL_LIST_BUFFER, x, y);
     return;
 }
@@ -391 +393 @@
 // Ensure the input array of pmStackData is valid, and get some details out of it
 static bool validateInputData(bool *haveVariances, // Do we have variance maps in the sky images?
-                              bool *havePixels, // Do we have lists of pixels?
+                              bool *haveRejects, // Do we have lists of rejected pixels?
                               int *num,    // Number of inputs
                               int *numCols, int *numRows, // Size of (sky) images
@@ -420 +422 @@
         PS_ASSERT_IMAGE_TYPE(data->readout->weight, PS_TYPE_F32, false);
     }
-    *havePixels = (data->pixels != NULL);
+    *haveRejects = (data->reject != NULL);
 
     // Make sure the rest correspond with the first
@@ -433 +435 @@
             return false;
         }
-        if ((*havePixels && !data->pixels) || (data->pixels && !*havePixels)) {
-            psError(PS_ERR_UNEXPECTED_NULL, true, "The pixels are specified in some but not all inputs.");
+        if ((*haveRejects && !data->reject) || (data->reject && !*haveRejects)) {
+            psError(PS_ERR_UNEXPECTED_NULL, true,
+                    "The rejected pixels are specified in some but not all inputs.");
             return false;
         }
@@ -473 +476 @@
             continue;
         }
-        assert(data->pixels);
-        psPixels *pixels = data->pixels;// The pixels of interest
+        assert(data->reject);
+        psPixels *pixels = data->reject; // The rejected pixels
         for (int j = 0; j < pixels->n; j++) {
             int x = pixels->data[j].x, y = pixels->data[j].y; // Coordinates of interest
@@ -518 +521 @@
 
     data->readout = psMemIncrRefCounter(readout);
-    data->pixels = NULL;
+    data->reject = NULL;
+    data->inspect = NULL;
     data->weight = weight;
 
@@ -526 +530 @@
 /// Stack input images
 bool pmStackCombine(pmReadout *combined, psArray *input, psMaskType maskVal, psMaskType bad,
-                    int kernelSize, int numIter, float rej, bool useVariance, bool safe)
+                    int kernelSize, int numIter, float rej, bool entire, bool useVariance, bool safe)
 {
     PS_ASSERT_PTR_NON_NULL(combined, false);
     bool haveVariances;                 // Do we have the variance maps?
-    bool havePixels;                    // Do we have lists of pixels?
+    bool haveRejects;                   // Do we have lists of rejected pixels?
     int num;                            // Number of inputs
     int numCols, numRows;               // Size of (sky) images
-    if (!validateInputData(&haveVariances, &havePixels, &num, &numCols, &numRows, input)) {
+    if (!validateInputData(&haveVariances, &haveRejects, &num, &numCols, &numRows, input)) {
         return false;
     }
@@ -544 +548 @@
         PS_ASSERT_FLOAT_LARGER_THAN(rej, 0.0, false);
     }
-    if (havePixels) {
+    if (haveRejects) {
         // This is a subsequent combination, so expect that the image and mask already exist
         PS_ASSERT_IMAGE_NON_NULL(combined->image, false);
@@ -557 +561 @@
     }
 
-    // Pull out the image weightings
-    psVector *weights = psVectorAlloc(num, PS_TYPE_F32);
+    psVector *weights = psVectorAlloc(num, PS_TYPE_F32); // Relative weighting for each image
+    psArray *stack = psArrayAlloc(num); // Stack of readouts
     for (int i = 0; i < num; i++) {
         pmStackData *data = input->data[i]; // Stack data for this input
         if (!data) {
             weights->data.F32[i] = 0.0;
+            continue;
         }
         weights->data.F32[i] = data->weight;
-    }
-
-    // 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
@@ -598 +598 @@
     combineBuffer *buffer = combineBufferAlloc(num);
 
-    if (havePixels) {
-        // Only combine select pixels
+    if (haveRejects) {
         psImage *combinedImage = combined->image; // Combined image
         psImage *combinedMask = combined->mask; // Combined mask
@@ -611 +610 @@
                 continue;
             }
-            pixels = psPixelsConcatenate(pixels, data->pixels);
-        }
-        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, combinedVariance, input, weights, reject, x, y,
-                          maskVal, bad, numIter, rej, useVariance, safe, buffer);
+            pixels = psPixelsConcatenate(pixels, data->reject);
+        }
+
+        if (entire) {
+            // Combine entire image
+            for (int y = minInputRows; y < maxInputRows; y++) {
+                for (int x = minInputCols; x < maxInputCols; x++) {
+                    psVector *reject = pixelMapQuery(pixelMap, x, y); // Inspect these images closely
+                    combinePixels(combinedImage, combinedMask, combinedVariance, input, weights, reject, x, y,
+                                  maskVal, bad, numIter, rej, useVariance, safe, buffer);
+                }
+            }
+        } else {
+            // Only combine previously rejected pixels
+            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, combinedVariance, input, weights, reject, x, y,
+                              maskVal, bad, numIter, rej, useVariance, safe, buffer);
+            }
         }
         psFree(pixels);
@@ -644 +656 @@
         }
 
-        // Generate the pixel lists in which to place the rejected pixels
-        if (numIter != 0) {
-            for (int i = 0; i < num; i++) {
-                pmStackData *data = input->data[i]; // Stack data for this input
-                if (!data) {
-                    continue;
-                }
-                data->pixels = psPixelsAllocEmpty(PIXEL_LIST_BUFFER);
-            }
-        }
-
         for (int y = minInputRows; y < maxInputRows; y++) {
             for (int x = minInputCols; x < maxInputCols; x++) {
@@ -669 +670 @@
                     continue;
                 }
-                psTrace("psModules.imcombine", 5, "Image %d: %ld pixels to inspect.\n", i, data->pixels->n);
+                psTrace("psModules.imcombine", 5, "Image %d: %ld pixels to inspect.\n", i, data->inspect->n);
             }
         }

trunk/psModules/src/imcombine/pmStack.h

@@ -8 +8 @@
  * @author GLG, MHPCC
  *
- * @version $Revision: 1.6 $ $Name: not supported by cvs2svn $
- * @date $Date: 2008-02-27 21:16:57 $
+ * @version $Revision: 1.7 $ $Name: not supported by cvs2svn $
+ * @date $Date: 2008-03-17 21:38:43 $
  *
  * Copyright 2004-2007 Institute for Astronomy, University of Hawaii
@@ -28 +28 @@
 typedef struct {
     pmReadout *readout;                 ///< Warped readout (sky cell)
-    psPixels *pixels;                   ///< Pixels to inspect or reject
-    float weight;                       ///< Weight to apply
+    psPixels *reject;                   ///< Pixels to reject
+    psPixels *inspect;                  ///< Pixels to inspect
+    float weight;                       ///< Relative weighting for image
 } pmStackData;
 
@@ -45 +46 @@
                     int numIter,        ///< Number of iterations
                     float rej,          ///< Rejection limit (standard deviations)
+                    bool entire,        ///< Combine entire image even if rejection lists provided?
                     bool useVariance,   ///< Use variance values for rejection?
                     bool safe           ///< Play safe with small numbers of input pixels (mask if N <= 2)?