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

Timestamp:

Mar 27, 2009, 3:08:33 PM (17 years ago)

Author:

Paul Price

Message:

Ensure we can reject pixels when combining images without convolution.

File:

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

trunk/psModules/src/imcombine/pmStack.c

@@ -46 +46 @@
     psVector *weights;                  // Pixel weightings
     psVector *sources;                  // Pixel sources (which image did they come from?)
+    psVector *limits;                   // Rejection limits
     psVector *sort;                     // Buffer for sorting (to get a robust estimator of the standard dev)
 } combineBuffer;
@@ -56 +57 @@
     psFree(buffer->weights);
     psFree(buffer->sources);
+    psFree(buffer->limits);
     psFree(buffer->sort);
     return;
@@ -71 +73 @@
     buffer->weights = psVectorAlloc(numImages, PS_TYPE_F32);
     buffer->sources = psVectorAlloc(numImages, PS_TYPE_U16);
+    buffer->limits = psVectorAlloc(numImages, PS_TYPE_F32);
     buffer->sort = psVectorAlloc(numImages, PS_TYPE_F32);
     return buffer;
@@ -319 +322 @@
                           bool useVariance, // Use variance for rejection when combining?
                           bool safe,    // Combine safely?
+                          bool rejectInspect, // Reject values marked for inspection from combination?
                           combineBuffer *buffer // Buffer for combination; to avoid multiple allocations
                          )
@@ -336 +340 @@
     psVector *pixelWeights = buffer->weights; // Image weights for the pixel of interest
     psVector *pixelSources = buffer->sources; // Sources for the pixel of interest
+    psVector *pixelLimits = buffer->limits; // Limits for the pixel of interest
     psVector *sort = buffer->sort;      // Sort buffer
 
@@ -375 +380 @@
     pixelWeights->n = num;
     pixelSources->n = num;
+    pixelLimits->n = num;
 
 #ifdef TESTING
@@ -389 +395 @@
     // Default option is that the pixel is bad
     float imageValue = NAN, varianceValue = NAN; // Value for combined image and variance map
-    psImageMaskType maskValue = bad;         // Value for combined mask
+    psImageMaskType maskValue = bad;    // Value for combined mask
     switch (num) {
       case 0:
@@ -449 +455 @@
               if (PS_SQR(diff) > PS_SQR(rej) * sigma2) {
                   // Not consistent: mark both for inspection
-                  combineInspect(inputs, x, y, pixelSources->data.U16[0]);
-                  combineInspect(inputs, x, y, pixelSources->data.U16[1]);
+                  if (rejectInspect) {
+                      imageValue = NAN;
+                      varianceValue = NAN;
+                      maskValue = bad;
+                  } else {
+                      combineInspect(inputs, x, y, pixelSources->data.U16[0]);
+                      combineInspect(inputs, x, y, pixelSources->data.U16[1]);
+                  }
 #ifdef TESTING
                   if (x == TEST_X && y == TEST_Y) {
@@ -501 +513 @@
 #endif
 
-                      pixelVariances->data.F32[i] = rej2 * (pixelVariances->data.F32[i] + sysVar);
+                      pixelLimits->data.F32[i] = rej2 * (pixelVariances->data.F32[i] + sysVar);
                   }
               }
@@ -541 +553 @@
 #define MASK_PIXEL_FOR_INSPECTION() \
     pixelMasks->data.PS_TYPE_VECTOR_MASK_DATA[j] = 0xff; \
-    combineInspect(inputs, x, y, pixelSources->data.U16[j]); \
+    if (!rejectInspect) { \
+        combineInspect(inputs, x, y, pixelSources->data.U16[j]); \
+    } \
     numClipped++; \
     totalClipped++;
@@ -553 +567 @@
                       // Comparing squares --- cheaper than lots of sqrts
                       // pixelVariances includes the rejection limit, from above
-                      if (PS_SQR(diff) > pixelVariances->data.F32[j]) {
+                      if (PS_SQR(diff) > pixelLimits->data.F32[j]) {
                           MASK_PIXEL_FOR_INSPECTION();
 #ifdef TESTING
                           if (x == TEST_X && y == TEST_Y) {
                               fprintf(stderr, "Rejecting input %d based on variance: %f > %f\n",
-                                      j, diff, sqrtf(pixelVariances->data.F32[j]));
+                                      j, diff, sqrtf(pixelLimits->data.F32[j]));
                           }
 #endif
@@ -573 +587 @@
               }
           }
+
+          if (rejectInspect && totalClipped > 0) {
+              // Get rid of the masked values
+              // The alternative to this is to make combinationMeanVariance() accept a mask
+              int good = 0;            // Index of good value
+              for (int i = 0; i < num; i++) {
+                  if (pixelMasks->data.PS_TYPE_VECTOR_MASK_DATA[i]) {
+                      continue;
+                  }
+                  if (i != good) {
+                      pixelData->data.F32[good] = pixelData->data.F32[i];
+                      pixelVariances->data.F32[good] = pixelVariances->data.F32[i];
+                      pixelWeights->data.F32[good] = pixelWeights->data.F32[i];
+                      pixelData->data.F32[good] = pixelData->data.F32[i];
+                  }
+                  good++;
+              }
+              pixelData->n = good;
+              pixelVariances->n = good;
+              pixelWeights->n = good;
+              if (combinationMeanVariance(&mean, &var, pixelData, pixelVariances, pixelWeights)) {
+                  imageValue = mean;
+                  varianceValue = var;
+                  maskValue = 0;
+              } else {
+                  imageValue = NAN;
+                  varianceValue = NAN;
+                  maskValue = bad;
+              }
+          }
+
           break;
       }
@@ -734 +779 @@
 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 entire, bool useVariance, bool safe, bool rejectInspect)
 {
     PS_ASSERT_PTR_NON_NULL(combined, false);
@@ -838 +883 @@
                     combinePixels(combinedImage, combinedMask, combinedVariance, input, weights,
                                   addVariance, reject, x, y, maskVal, bad, numIter, rej, sys, discard,
-                                  useVariance, safe, buffer);
+                                  useVariance, safe, rejectInspect, buffer);
                 }
             }
@@ -853 +898 @@
                 combinePixels(combinedImage, combinedMask, combinedVariance, input, weights,
                               addVariance, reject, x, y, maskVal, bad, numIter, rej, sys, discard,
-                              useVariance, safe, buffer);
+                              useVariance, safe, rejectInspect, buffer);
             }
         }
@@ -881 +926 @@
                 combinePixels(combinedImage, combinedMask, combinedVariance, input, weights,
                               addVariance, NULL, x, y, maskVal, bad, numIter, rej, sys, discard,
-                              useVariance, safe, buffer);
+                              useVariance, safe, rejectInspect, buffer);
             }
         }