Changeset 14272 for trunk/psModules/src/imcombine/pmStack.c

Timestamp:

Jul 17, 2007, 1:33:09 PM (19 years ago)

Author:

Paul Price

Message:

Adding input/output for weight maps (variances). Using local
functions instead of psVectorStats because it is slow, and it assumes
that the input weights are errors, and I don't want to have to do
square-roots all over the place, etc.

File:

• trunk/psModules/src/imcombine/pmStack.c (modified) (19 diffs)

trunk/psModules/src/imcombine/pmStack.c

@@ -8 +8 @@
  *  @author GLG, MHPCC
  *
- *  @version $Revision: 1.7 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2007-07-11 01:32:41 $
+ *  @version $Revision: 1.8 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-07-17 23:33:09 $
  *
  *  Copyright 2004-2007 Institute for Astronomy, University of Hawaii
@@ -41 +41 @@
     psVector *pixels;                   // Pixel values
     psVector *masks;                    // Pixel masks
-    psStats *stats;                     // Statistics
+    psVector *weights;                  // Pixel weights
+    psVector *sort;                     // Buffer for sorting (to get a robust estimator of the standard dev)
 } combineBuffer;
 
@@ -48 +49 @@
     psFree(buffer->pixels);
     psFree(buffer->masks);
-    psFree(buffer->stats);
+    psFree(buffer->weights);
+    psFree(buffer->sort);
     return;
 }
@@ -61 +63 @@
     buffer->pixels = psVectorAlloc(numImages, PS_TYPE_F32);
     buffer->masks = psVectorAlloc(numImages, PS_TYPE_MASK);
-    buffer->stats = psStatsAlloc(stat);
+    buffer->weights = psVectorAlloc(numImages, PS_TYPE_F32);
+    buffer->sort = psVectorAlloc(numImages, PS_TYPE_F32);
 
     return buffer;
@@ -77 +80 @@
 
 
+// Generate a mean value for the combination
+// Not using psVectorStats because it assumes that the weights are errors, and weights by 1/error^2
+static bool combinationMean(float *mean, // Mean value, to return
+                            const psVector *values, // Values to combine
+                            const psVector *weights, // Weights to apply
+                            const psVector *masks, // Mask to apply
+                            psMaskType maskVal // Value to mask
+                            )
+{
+    assert(values && weights && masks);
+    assert(values->n == weights->n);
+    assert(values->n == masks->n);
+    assert(values->type.type == PS_TYPE_F32);
+    assert(weights->type.type == PS_TYPE_F32);
+    assert(masks->type.type == PS_TYPE_MASK);
+
+    float sumValueWeight = 0.0; // Sum of the value multiplied by the weight
+    float sumWeight = 0.0;      // Sum of the weights
+    for (int i = 0; i < values->n; i++) {
+        if (!(masks->data.PS_TYPE_MASK_DATA[i] & maskVal)) {
+            sumValueWeight += values->data.F32[i] * weights->data.F32[i];
+            sumWeight += weights->data.F32[i];
+        }
+    }
+    if (sumWeight <= 0) {
+        return false;
+    }
+
+    *mean = sumValueWeight / sumWeight;
+    return true;
+}
+
+// Generate a standard deviation for the combination
+// Not using psVectorStats because it has additional allocations which slow things down
+static bool combinationStdev(float *stdev, // Mean value, to return
+                             const psVector *values, // Values to combine
+                             const psVector *masks, // Mask to apply
+                             psMaskType maskVal, // Value to mask
+                             psVector *sortBuffer // Buffer for sorting
+                             )
+{
+    assert(values && masks);
+    assert(values->n == masks->n);
+    assert(values->type.type == PS_TYPE_F32);
+    assert(masks->type.type == PS_TYPE_MASK);
+    assert(sortBuffer && sortBuffer->nalloc >= values->n && sortBuffer->type.type == PS_TYPE_F32);
+
+    int num = 0;            // Number of valid values
+    for (int i = 0; i < values->n; i++) {
+        if (!(masks->data.PS_TYPE_MASK_DATA[i] & maskVal)) {
+            sortBuffer->data.F32[num++] = values->data.F32[i];
+        }
+    }
+    sortBuffer->n = num;
+    if (!psVectorSortInPlace(sortBuffer)) {
+        return false;
+    }
+    *stdev = 0.74 * (sortBuffer->data.F32[(int)(0.75 * num)] - sortBuffer->data.F32[(int)(0.25 * num)]);
+
+    return true;
+}
+
+// Generate a variance value for the combination
+static bool combinationVariance(float *variance, // Variance value, to return
+                                const psVector *variances, // Pixel variances to combine
+                                const psVector *weights, // Image weights to apply
+                                const psVector *masks, // Mask to apply
+                                psMaskType maskVal // Value to mask
+                                )
+{
+    assert(variances && weights && masks);
+    assert(variances->n == weights->n);
+    assert(variances->n == masks->n);
+    assert(variances->type.type == PS_TYPE_F32);
+    assert(weights->type.type == PS_TYPE_F32);
+    assert(masks->type.type == PS_TYPE_MASK);
+
+    // Get the variance in the combination.  We're not using the input pixel variances to generate a
+    // weighted average for the pixel flux (because that introduces systematic biases), so the variance
+    // of the output pixel value should simply be:
+    //     simga^2 = sum(weight_i^2 * sigma_i^2) / (sum(weight_i))^2
+    // This reduces, when the weights are all identically unity, to:
+    //     variance_combination = sum(variance_i) / N^2
+    // and if the variances are all equal:
+    //     variance_combination = variance_individual / N
+    // which makes sense --- the standard deviation of the combination is reduced by a factor of sqrt(N).
+
+    float sumWeights = 0.0;             // Sum of the global weights
+    float sumVarianceWeights = 0.0;     // Sum of the pixel variances multiplied by the global weights
+    for (int i = 0; i < variances->n; i++) {
+        if (!(masks->data.PS_TYPE_MASK_DATA[i] & maskVal)) {
+            sumVarianceWeights += variances->data.F32[i] * PS_SQR(weights->data.F32[i]);
+            sumWeights += weights->data.F32[i];
+        }
+    }
+
+    if (sumWeights <= 0) {
+        return false;
+    }
+
+    *variance = sumVarianceWeights / PS_SQR(sumWeights);
+    return true;
+}
+
 // Given a stack of images, combine with optional rejection.
 // Pixels in the stack that are rejected are marked for subsequent
 static bool combinePixels(psImage *image, // Combined image, for output
                           psImage *mask, // Combined mask, for output
+                          psImage *weight, // Combined variance map, for output
                           const psArray *inputs, // Stack data
-                          const psVector *weights, // Weights for data, or NULL
+                          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
@@ -104 +212 @@
     psVector *pixelData = buffer->pixels; // Values for the pixel of interest
     psVector *pixelMasks = buffer->masks; // Masks for the pixel of interest
-    psStats *stats = buffer->stats;     // Statistics
-
+    psVector *pixelWeights = buffer->weights; // Weights for the pixel of interest
+    psVector *sort = buffer->sort;      // Sort buffer
 
     // Extract the pixel and mask data
@@ -112 +220 @@
         pmStackData *data = inputs->data[i]; // Stack data of interest
         psImage *image = data->sky->image; // Image of interest
+        psImage *weight = data->sky->weight; // Weight map of interest
         psImage *mask = data->sky->mask; // Mask of interest
         pixelData->data.F32[i] = image->data.F32[y][x];
+        pixelWeights->data.F32[i] = weight->data.F32[y][x];
         pixelMasks->data.PS_TYPE_MASK_DATA[i] = mask->data.PS_TYPE_MASK_DATA[y][x];
         if (pixelMasks->data.PS_TYPE_MASK_DATA[i] & maskVal) {
@@ -122 +232 @@
         // Catch the case where everything is masked
         image->data.F32[y][x] = NAN;
+        if (weight) {
+            weight->data.F32[y][x] = NAN;
+        }
         mask->data.PS_TYPE_MASK_DATA[y][x] = bad;
         return false;
@@ -134 +247 @@
     // Record the value derived with no clipping, because pixels rejected using the harsh clipping applied in
     // the first pass might later be accepted.
-    stats->options |= PS_STAT_SAMPLE_MEAN;
-    if (!psVectorStats(stats, pixelData, weights, pixelMasks, maskVal)) {
-        // Can't do anything about an error except give it a NAN and mask
-        psErrorClear();
+    float mean, variance;               // Mean and variance of the combination
+    if (!combinationMean(&mean, pixelData, weights, pixelMasks, maskVal) ||
+        (weight && !combinationVariance(&variance, weights, pixelWeights, pixelMasks, maskVal))) {
         image->data.F32[y][x] = NAN;
         mask->data.PS_TYPE_MASK_DATA[y][x] = bad;
-        return false;
-    }
-    image->data.F32[y][x] = stats->sampleMean;
+        if (weight) {
+            weight->data.F32[y][x] = NAN;
+        }
+        return false;
+    }
+
+    image->data.F32[y][x] = mean;
+    if (weight) {
+        weight->data.F32[y][x] = variance;
+    }
     mask->data.PS_TYPE_MASK_DATA[y][x] = 0;
 
-    stats->options &= ~ PS_STAT_SAMPLE_MEAN;
-
+    // The clipping that follows is solely to identify suspect pixels.
+    // These suspect pixels will be inspected in more detail by other functions.
     long numClipped = LONG_MAX;         // Number of pixels clipped
     psMaskType ignore = maskVal | bad;  // Ignore values with this mask value
     for (int i = 0; i < numIter && numClipped > 0; i++) {
-#if 1
-        float mean = stats->sampleMedian;
-        float stdev = 0.74 * (stats->sampleUQ - stats->sampleLQ); // Rough estimate of the standard deviation
-#else
-        float mean = stats->robustMedian;
-        float stdev = stats->robustStdev;
-#endif
+        float stdev;                    // Standard deviation of the combination, for rejection
+        // Only get the mean if it's not the first iteration (because we got the mean earlier)
+        if ((i != 0 && !combinationMean(&mean, pixelData, weights, pixelMasks, maskVal)) ||
+            !combinationStdev(&stdev, pixelData, pixelMasks, maskVal, sort)) {
+            image->data.F32[y][x] = NAN;
+            mask->data.PS_TYPE_MASK_DATA[y][x] = bad;
+            weight->data.F32[y][x] = NAN;
+            return false;
+        }
+
         float limit = rej * stdev; // Rejection limit
         numClipped = 0;
@@ -165 +287 @@
             float diff = pixelData->data.F32[j] - mean;
             if (fabsf(diff) > limit) {
+                // Add the pixel as one to inspect
                 pmStackData *data = inputs->data[j]; // Stack data of interest
                 data->pixels = psPixelsAdd(data->pixels, PIXEL_LIST_BUFFER, x, y);
+                // Mask it so it's not considered in other iterations within this function
                 pixelMasks->data.PS_TYPE_MASK_DATA[j] |= bad;
                 numClipped++;
             }
         }
-        if (!psVectorStats(stats, pixelData, weights, pixelMasks, maskVal)) {
-            // Can't do anything about an error except give it a NAN and mask
-            psErrorClear();
-            image->data.F32[y][x] = NAN;
-            mask->data.PS_TYPE_MASK_DATA[y][x] = bad;
-            return false;
-        }
     }
 
@@ -186 +303 @@
 // Ensure the input array of pmStackData is valid, and get some details out of it
 static bool validateInputData(bool *haveDetector, // Do we have the detector images?
-                           bool *haveSky, // Do we have the sky images?
-                           bool *havePixels, // Do we have lists of pixels?
-                           int *num,    // Number of inputs
-                           int *numCols, int *numRows, // Size of (sky) images
-                           psArray *input // Input array of pmStackData to validate
+                              bool *haveSky, // Do we have the sky images?
+                              bool *haveSkyWeights, // Do we have weights in the sky images?
+                              bool *havePixels, // Do we have lists of pixels?
+                              int *num,    // Number of inputs
+                              int *numCols, int *numRows, // Size of (sky) images
+                              psArray *input // Input array of pmStackData to validate
     )
 {
@@ -203 +321 @@
             PS_ASSERT_IMAGE_NON_NULL(data->detector->image, false);
             PS_ASSERT_IMAGE_NON_NULL(data->detector->mask, false);
+            PS_ASSERT_IMAGES_SIZE_EQUAL(data->detector->image, data->detector->mask, false);
+            PS_ASSERT_IMAGE_TYPE(data->detector->image, PS_TYPE_F32, false);
+            PS_ASSERT_IMAGE_TYPE(data->detector->mask, PS_TYPE_MASK, false);
         }
         *haveSky = (data->sky != NULL);
+        *haveSkyWeights = false;
         if (*haveSky) {
             PS_ASSERT_IMAGE_NON_NULL(data->sky->image, false);
+            PS_ASSERT_IMAGE_TYPE(data->sky->image, PS_TYPE_F32, false);
             PS_ASSERT_IMAGE_NON_NULL(data->sky->mask, false);
+            PS_ASSERT_IMAGE_TYPE(data->sky->mask, PS_TYPE_MASK, false);
+            PS_ASSERT_IMAGES_SIZE_EQUAL(data->sky->image, data->sky->mask, false);
             *numCols = data->sky->image->numCols;
             *numRows = data->sky->image->numRows;
+            if (data->sky->weight) {
+                *haveSkyWeights = true;
+                PS_ASSERT_IMAGE_NON_NULL(data->sky->weight, false);
+                PS_ASSERT_IMAGES_SIZE_EQUAL(data->sky->image, data->sky->weight, false);
+                PS_ASSERT_IMAGE_TYPE(data->sky->weight, PS_TYPE_F32, false);
+            }
         }
         *havePixels = (data->pixels != NULL);
@@ -243 +374 @@
             PS_ASSERT_IMAGE_SIZE(data->sky->image, *numCols, *numRows, false);
             PS_ASSERT_IMAGES_SIZE_EQUAL(data->sky->image, data->sky->mask, false);
+            if (*haveSkyWeights) {
+                PS_ASSERT_IMAGE_NON_NULL(data->sky->weight, false);
+                PS_ASSERT_IMAGES_SIZE_EQUAL(data->sky->image, data->sky->weight, false);
+                PS_ASSERT_IMAGE_TYPE(data->sky->weight, PS_TYPE_F32, false);
+            }
         }
     }
@@ -412 +548 @@
     bool haveDetector;                  // Do we have the detector images?
     bool haveSky;                       // Do we have the sky images?
+    bool haveSkyWeights;                // Do we have the sky weight maps?
     bool havePixels;                    // Do we have lists of pixels?
     int num;                            // Number of inputs
     int numCols, numRows;               // Size of (sky) images
-    if (!validateInputData(&haveDetector, &haveSky, &havePixels, &num, &numCols, &numRows, input)) {
+    if (!validateInputData(&haveDetector, &haveSky, &haveSkyWeights, &havePixels, &num,
+                           &numCols, &numRows, input)) {
         return false;
     }
@@ -461 +599 @@
         psImage *combinedImage = combined->image; // Combined image
         psImage *combinedMask = combined->mask; // Combined mask
+        psImage *combinedWeight = combined->weight; // Combined mask
 
         psArray *pixelMap = pixelMapGenerate(input, numCols, numRows); // Map of pixels to source
@@ -472 +611 @@
             int x = pixels->data[i].x, y = pixels->data[i].y; // Coordinates of interest
             psVector *reject = pixelMapQuery(pixelMap, x, y); // Inspect these images closely
-            combinePixels(combinedImage, combinedMask, input, weights, reject, x, y,
+            combinePixels(combinedImage, combinedMask, combinedWeight, input, weights, reject, x, y,
                           maskVal, bad, numIter, rej, buffer);
         }
@@ -491 +630 @@
         }
 
+        psImage *combinedWeights = combined->weight; // Combined weight map
+        if (haveSkyWeights && !combinedWeights) {
+            combined->weight = psImageAlloc(numCols, numRows, PS_TYPE_F32);
+            combinedWeights = combined->weight;
+        }
+
         // Generate the pixel lists in which to place the rejected pixels
         if (numIter != 0) {
@@ -501 +646 @@
         for (int y = 0; y < numRows; y++) {
             for (int x = 0; x < numCols; x++) {
-                combinePixels(combinedImage, combinedMask, input, weights, NULL, x, y,
+                combinePixels(combinedImage, combinedMask, combinedWeights, input, weights, NULL, x, y,
                               maskVal, bad, numIter, rej, buffer);
             }
@@ -538 +683 @@
     bool haveDetector;                  // Do we have the detector images?
     bool haveSky;                       // Do we have the sky images?
+    bool haveSkyWeights;                // Do we have the sky weight maps?
     bool havePixels;                    // Do we have lists of pixels?
     int num;                            // Number of inputs
     int numCols, numRows;               // Size of (sky) images
-    if (!validateInputData(&haveDetector, &haveSky, &havePixels, &num, &numCols, &numRows, input)) {
+    if (!validateInputData(&haveDetector, &haveSky, &haveSkyWeights, &havePixels,
+                           &num, &numCols, &numRows, input)) {
         return false;
     }