Changeset 25964 for branches/pap/psModules/src/imcombine/pmStack.c

Timestamp:

Oct 28, 2009, 5:35:05 PM (17 years ago)

Author:

Paul Price

Message:

Reworking stack combination because there are *three* modes for pixels going into the final stack (after rejection), not just two: tested and good, tested and rejected, and not tested. The code did not recognise the third, which is a distinct state because we don't want these pixels grown, as we do for rejected pixels. This cannot be fixed merely by using the 'safe' combination because that would discard the 'tested and good' pixels that have only a single unrejected input but are good because they have survived the testing process. Needed to add a new state into the combination process. Now I add these pixels straight into the 'reject' list. This requires a little bit more fiddling around in ppStack. Not sure it's working yet.

File:

• branches/pap/psModules/src/imcombine/pmStack.c (modified) (16 diffs)

branches/pap/psModules/src/imcombine/pmStack.c

@@ -290 +290 @@
 
 // Mark a pixel for inspection
+// Value in pixel doesn't seem to agree with the stack, so need to look closer
 static inline void combineInspect(const psArray *inputs, // Stack data
                                   int x, int y, // Pixel
@@ -301 +302 @@
     }
     data->inspect = psPixelsAdd(data->inspect, data->inspect->nalloc, x, y);
+    return;
+}
+
+// Mark a pixel for rejection
+// Cannot possibly inspect this pixel and confirm that it's good.
+// e.g., Only a single input
+static inline void combineReject(const psArray *inputs, // Stack data
+                                 int x, int y, // Pixel
+                                 int source // Source image index
+                                 )
+{
+    pmStackData *data = inputs->data[source]; // Stack data of interest
+    if (!data) {
+        psWarning("Can't find input data for source %d", source);
+        return;
+    }
+    data->reject = psPixelsAdd(data->reject, data->reject->nalloc, x, y);
     return;
 }
@@ -319 +337 @@
                           float rej, // Number of standard deviations at which to reject
                           float sys,    // Relative systematic error
-                          float discard,// Fraction of values to discard (Olympic weighted mean)
+                          float olympic,// Fraction of values to discard (Olympic weighted mean)
                           bool useVariance, // Use variance for rejection when combining?
                           bool safe,    // Combine safely?
-                          bool rejectInspect, // Reject values marked for inspection from combination?
+                          bool rejection, // Reject values marked for inspection from combination?
                           combineBuffer *buffer // Buffer for combination; to avoid multiple allocations
                          )
@@ -419 +437 @@
               }
               maskValue = 0;
-          }
-#ifdef TESTING
-          else {
+          } else {
+              if (!rejection) {
+                  combineReject(inputs, x, y, pixelSources->data.U16[0]);
+              }
+#ifdef TESTING
               numRejected = 1;
               if (x == TEST_X && y == TEST_Y) {
                   fprintf(stderr, "Single input to combine, safety on, pixel is bad.\n");
               }
+#endif
           }
-#endif
           break;
       }
@@ -459 +479 @@
               if (PS_SQR(diff) > PS_SQR(rej) * sigma2) {
                   // Not consistent: mark both for inspection
-                  if (rejectInspect) {
+                  if (rejection) {
                       imageValue = NAN;
                       varianceValue = NAN;
@@ -550 +570 @@
                   }
 #endif
-                  median = combinationWeightedOlympic(pixelData, pixelWeights, pixelMasks, discard, sort);
+                  median = combinationWeightedOlympic(pixelData, pixelWeights, pixelMasks, olympic, sort);
               }
 
@@ -563 +583 @@
 #define MASK_PIXEL_FOR_INSPECTION() \
     pixelMasks->data.PS_TYPE_VECTOR_MASK_DATA[j] = 0xff; \
-    if (!rejectInspect) { \
+    if (!rejection) { \
         combineInspect(inputs, x, y, pixelSources->data.U16[j]); \
     } \
@@ -605 +625 @@
           }
 
-          if (rejectInspect && totalClipped > 0) {
+          if (rejection && totalClipped > 0) {
               // Get rid of the masked values
               // The alternative to this is to make combinationMeanVariance() accept a mask
@@ -658 +678 @@
 // 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 *haveRejects, // Do we have lists of rejected pixels?
                               int *num,    // Number of inputs
                               int *numCols, int *numRows, // Size of (sky) images
-                              psArray *input // Input array of pmStackData to validate
+                              const psArray *input, // Input array of pmStackData to validate
+                              const pmReadout *output // Output readout
     )
 {
@@ -687 +707 @@
         PS_ASSERT_IMAGE_TYPE(data->readout->variance, PS_TYPE_F32, false);
     }
-    *haveRejects = (data->reject != NULL);
+    bool haveRejects = (data->reject != NULL); // Do we have rejected pixels?
 
     // Make sure the rest correspond with the first
@@ -700 +720 @@
             return false;
         }
-        if ((*haveRejects && !data->reject) || (data->reject && !*haveRejects)) {
+        if ((haveRejects && !data->reject) || (data->reject && !haveRejects)) {
             psError(PS_ERR_UNEXPECTED_NULL, true,
                     "The rejected pixels are specified in some but not all inputs.");
@@ -716 +736 @@
             PS_ASSERT_IMAGE_TYPE(data->readout->variance, PS_TYPE_F32, false);
         }
+    }
+
+    PM_ASSERT_READOUT_NON_NULL(output, false);
+    if (output->image) {
+        PS_ASSERT_IMAGE_NON_NULL(output->image, false);
+        PS_ASSERT_IMAGE_TYPE(output->image, PS_TYPE_F32, false);
+        PS_ASSERT_IMAGE_NON_NULL(output->mask, false);
+        PS_ASSERT_IMAGE_TYPE(output->mask, PS_TYPE_IMAGE_MASK, false);
+        PS_ASSERT_IMAGES_SIZE_EQUAL(output->image, output->mask, false);
     }
 
@@ -802 +831 @@
 /// Stack input images
 bool pmStackCombine(pmReadout *combined, psArray *input, psImageMaskType maskVal, psImageMaskType bad,
-                    int kernelSize, int numIter, float rej, float sys, float discard,
-                    bool entire, bool useVariance, bool safe, bool rejectInspect)
-{
-    PS_ASSERT_PTR_NON_NULL(combined, false);
+                    int kernelSize, int numIter, float rej, float sys, float olympic,
+                    bool useVariance, bool safe, bool rejection)
+{
     bool haveVariances;                 // Do we have the variance maps?
-    bool haveRejects;                   // Do we have lists of rejected pixels?
     int num;                            // Number of inputs
     int numCols, numRows;               // Size of (sky) images
-    if (!validateInputData(&haveVariances, &haveRejects, &num, &numCols, &numRows, input)) {
+    if (!validateInputData(&haveVariances, &num, &numCols, &numRows, input, combined)) {
         return false;
     }
@@ -821 +848 @@
         PS_ASSERT_FLOAT_LARGER_THAN(rej, 0.0, false);
     }
-    if (haveRejects) {
-        // This is a subsequent combination, so expect that the image and mask already exist
-        PS_ASSERT_IMAGE_NON_NULL(combined->image, false);
-        PS_ASSERT_IMAGE_TYPE(combined->image, PS_TYPE_F32, false);
-        PS_ASSERT_IMAGE_NON_NULL(combined->mask, false);
-        PS_ASSERT_IMAGE_TYPE(combined->mask, PS_TYPE_IMAGE_MASK, false);
-        PS_ASSERT_IMAGES_SIZE_EQUAL(combined->image, combined->mask, false);
-    }
     if (useVariance && !haveVariances) {
         psWarning("Unable to use variance in rejection if no variance maps supplied --- option turned off");
@@ -852 +871 @@
         }
 #endif
-        if (!haveRejects && !data->inspect) {
-            data->inspect = psPixelsAllocEmpty(PIXEL_LIST_BUFFER);
+        if (!rejection) {
+            // Ensure pixels can be put on the appropriate list
+            if (!data->inspect) {
+                data->inspect = psPixelsAllocEmpty(PIXEL_LIST_BUFFER);
+            }
+            if (!data->reject) {
+                data->reject = psPixelsAllocEmpty(PIXEL_LIST_BUFFER);
+            }
         }
     }
@@ -882 +907 @@
     combineBuffer *buffer = combineBufferAlloc(num);
 
-    if (haveRejects) {
-        psImage *combinedImage = combined->image; // Combined image
-        psImage *combinedMask = combined->mask; // Combined mask
-        psImage *combinedVariance = combined->variance; // Combined variance map
-
-        psArray *pixelMap = pixelMapGenerate(input, minInputCols, maxInputCols,
-                                             minInputRows, maxInputRows); // Map of pixels to source
-        psPixels *pixels = NULL;            // Total list of pixels, with no duplicates
+    // Pull the products out, allocate if necessary
+    psImage *combinedImage = combined->image; // Combined image
+    if (!combinedImage) {
+        combined->image = psImageAlloc(numCols, numRows, PS_TYPE_F32);
+        combinedImage = combined->image;
+    }
+    psImage *combinedMask = combined->mask; // Combined mask
+    if (!combinedMask) {
+        combined->mask = psImageAlloc(numCols, numRows, PS_TYPE_IMAGE_MASK);
+        combinedMask = combined->mask;
+    }
+
+    psImage *combinedVariance = combined->variance; // Combined variance map
+    if (haveVariances && !combinedVariance) {
+        combined->variance = psImageAlloc(numCols, numRows, PS_TYPE_F32);
+        combinedVariance = combined->variance;
+    }
+
+    // Set up rejection list
+    psArray *pixelMap = NULL;           // Map of pixels to source
+    if (rejection) {
+        pixelMap = pixelMapGenerate(input, minInputCols, maxInputCols, minInputRows, maxInputRows);
+    }
+
+    // Combine each pixel
+    for (int y = minInputRows; y < maxInputRows; y++) {
+        for (int x = minInputCols; x < maxInputCols; x++) {
+            psVector *reject = NULL; // Images to reject for this pixel
+            if (rejection) {
+                reject = pixelMapQuery(pixelMap, minInputCols, minInputRows, x, y);
+            }
+            combinePixels(combinedImage, combinedMask, combinedVariance, input, weights,
+                          addVariance, reject, x, y, maskVal, bad, numIter, rej, sys, olympic,
+                          useVariance, safe, rejection, buffer);
+        }
+    }
+
+    psFree(pixelMap);
+    psFree(weights);
+    psFree(buffer);
+
+#ifndef PS_NO_TRACE
+    if (!rejection && psTraceGetLevel("psModules.imcombine") >= 5) {
         for (int i = 0; i < num; i++) {
-            pmStackData *data = input->data[i]; // Stacking data; contains the list of pixels
-            if (!data) {
+            pmStackData *data = input->data[i]; // Stack data for this input
+            if (!data || !data->inspect) {
                 continue;
             }
-            pixels = psPixelsConcatenate(pixels, data->reject);
-        }
-        pixels = psPixelsDuplicates(pixels, pixels);
-
-        if (entire) {
-            // Combine entire image
-            for (int y = minInputRows; y < maxInputRows; y++) {
-                for (int x = minInputCols; x < maxInputCols; x++) {
-                    psVector *reject = pixelMapQuery(pixelMap, minInputCols, minInputRows,
-                                                     x, y); // Reject these images
-                    combinePixels(combinedImage, combinedMask, combinedVariance, input, weights,
-                                  addVariance, reject, x, y, maskVal, bad, numIter, rej, sys, discard,
-                                  useVariance, safe, rejectInspect, 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, minInputCols, minInputRows,
-                                                 x, y); // Reject these images
-                combinePixels(combinedImage, combinedMask, combinedVariance, input, weights,
-                              addVariance, reject, x, y, maskVal, bad, numIter, rej, sys, discard,
-                              useVariance, safe, rejectInspect, buffer);
-            }
-        }
-        psFree(pixels);
-        psFree(pixelMap);
-    } else {
-        // Pull the products out, allocate if necessary
-        psImage *combinedImage = combined->image; // Combined image
-        if (!combinedImage) {
-            combined->image = psImageAlloc(numCols, numRows, PS_TYPE_F32);
-            combinedImage = combined->image;
-        }
-        psImage *combinedMask = combined->mask; // Combined mask
-        if (!combinedMask) {
-            combined->mask = psImageAlloc(numCols, numRows, PS_TYPE_IMAGE_MASK);
-            combinedMask = combined->mask;
-        }
-
-        psImage *combinedVariance = combined->variance; // Combined variance map
-        if (haveVariances && !combinedVariance) {
-            combined->variance = psImageAlloc(numCols, numRows, PS_TYPE_F32);
-            combinedVariance = combined->variance;
-        }
-
-        for (int y = minInputRows; y < maxInputRows; y++) {
-            for (int x = minInputCols; x < maxInputCols; x++) {
-                combinePixels(combinedImage, combinedMask, combinedVariance, input, weights,
-                              addVariance, NULL, x, y, maskVal, bad, numIter, rej, sys, discard,
-                              useVariance, safe, rejectInspect, buffer);
-            }
-        }
-
-#ifndef PS_NO_TRACE
-        if (psTraceGetLevel("psModules.imcombine") >= 5) {
-            for (int i = 0; i < num; i++) {
-                pmStackData *data = input->data[i]; // Stack data for this input
-                if (!data || !data->inspect) {
-                    continue;
-                }
-                psTrace("psModules.imcombine", 5, "Image %d: %ld pixels to inspect.\n", i, data->inspect->n);
-            }
-        }
-#endif
-    }
-
-    psFree(weights);
-    psFree(buffer);
+            psTrace("psModules.imcombine", 5, "Image %d: %ld pixels to inspect.\n", i, data->inspect->n);
+        }
+    }
+#endif
 
     return true;