Changeset 41843 for branches/eam_branches/ipp-dev-20210817/psModules/src/imcombine/pmStack.c

Timestamp:

Oct 16, 2021, 3:54:37 PM (5 years ago)

Author:

eugene

Message:

add utility command to report on pmFPAfiles; fix typo; add quick combine-by-percentile; remove some old and excessively-verbose test blocks

File:

• branches/eam_branches/ipp-dev-20210817/psModules/src/imcombine/pmStack.c (modified) (4 diffs)

branches/eam_branches/ipp-dev-20210817/psModules/src/imcombine/pmStack.c

@@ -961 +961 @@
 
     // KMM values;
-    float Punimodal = 1.0,KMMmean = NAN,KMMsigma = NAN,KMMpi = NAN;
+    float Punimodal = 1.0, KMMmean = NAN, KMMsigma = NAN, KMMpi = NAN;
     int KMM_MINIMUM_INPUTS = 6;
     bool useKMM = false;
@@ -1510 +1510 @@
 }
 
+# define MIN_GOOD_PERCENTILE 2
+bool pmStackCombineByPercentile(
+    pmReadout *combined,
+    psArray *stackData,
+    psF64 minRange,
+    psF64 maxRange,
+    psImageMaskType badMaskBits,        // treat these bits as 'bad'
+    psImageMaskType suspectMaskBits,    // treat these bits as 'suspect'
+    psImageMaskType blankMaskBits       // use this mask value for pixels missing input data (distinguish between Ninput = 0 and Ngood = 0?)
+) {
+    int minInputCols, maxInputCols, minInputRows, maxInputRows; // Smallest and largest values to combine
+    int xSize, ySize;                   // Size of the output image
+
+    // we need to copy the readouts to their own array for validation
+    psArray *stackReadouts = psArrayAlloc(stackData->n);
+    for (int i = 0; i < stackData->n; i++) {
+        pmStackData *data = stackData->data[i]; // Stack data for this input
+        pmReadout *ro = data->readout;  // Readout of interest
+        stackReadouts->data[i] = psMemIncrRefCounter(ro);
+    }
+
+    if (!pmReadoutStackValidate(&minInputCols, &maxInputCols, &minInputRows, &maxInputRows, &xSize, &ySize, stackReadouts)) {
+        psError(psErrorCodeLast(), false, "Input stack is not valid.");
+        psFree(stackReadouts);
+        return false;
+    }
+    psFree(stackReadouts);
+
+    psVector *pixelData = psVectorAlloc(stackData->n, PS_TYPE_F32);
+    psVector *pixelMask = psVectorAlloc(stackData->n, PS_TYPE_VECTOR_MASK);
+
+    for (int y = minInputRows; y < maxInputRows; y++) {
+        for (int x = minInputCols; x < maxInputCols; x++) {
+
+            int nGood = 0; 
+            for (int i = 0; i < stackData->n; i++) {
+
+                pmStackData *data = stackData->data[i]; // Stack data for this input
+                pmReadout *ro = data->readout;  // Readout of interest
+                psImage *image = ro->image;
+
+                int xIn = x - data->readout->col0;
+                int yIn = y - data->readout->row0; // Coordinates on input readout
+
+                if (!isfinite(image->data.F32[yIn][xIn])) continue;
+                if (fabs(image->data.F32[yIn][xIn]) > 1e5) continue;
+                // XXX need to test the input mask as well here
+
+                pixelData->data.F32[nGood] = image->data.F32[yIn][xIn];
+                nGood ++;
+            }
+            pixelData->n = nGood;
+            psVectorSortInPlace (pixelData);
+
+            if (nGood < MIN_GOOD_PERCENTILE) {
+                combined->image->data.F32[y][x] = NAN;
+                combined->variance->data.F32[y][x] = NAN;
+                combined->mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] = 1;
+                continue;
+            }
+
+            int Ns = MIN(MAX(0, minRange * nGood),nGood); // e.g., 0.1 * 50 = 5
+            int Ne = MIN(MAX(0, maxRange * nGood),nGood); // e.g., 0.9 * 50 = 45 
+            int Npt = Ne - Ns;
+
+            float sum = 0.0;
+            for (int n = Ns; n < Ne; n++) {
+                sum += pixelData->data.F32[n];
+            }
+            float mean = sum / (float) Npt;
+
+            float varSum = 0.0;
+            for (int n = Ns; n < Ne; n++) {
+                varSum += SQ(pixelData->data.F32[n] - mean);
+            }
+            // variance on the mean (stdev / sqrt(N))^2
+            float varValue = varSum / (Npt - 1) / Npt;
+
+            combined->image->data.F32[y][x] = mean;
+            combined->variance->data.F32[y][x] = varValue;
+            combined->mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] = 0;
+        }
+    }
+  
+    psFree(pixelData);
+    psFree(pixelMask);
+
+    return (true);
+}
+
 /// Stack input images
 bool pmStackCombine(
@@ -1644 +1734 @@
         expweight = expmaps->variance;
     }
-
-/*     // Pre-reject inputs using KMM bimodality test. */
-/*     if (1)  { */
-/* /\*       KMMRejectUnpopular(input,x,y); *\/ */
-/*       rejection = true; */
-/*     } */
    
     // Set up rejection list
@@ -1660 +1744 @@
     for (int y = minInputRows; y < maxInputRows; y++) {
         for (int x = minInputCols; x < maxInputCols; x++) {
+            // this is only for TESTING:
             CHECKPIX(x, y, "Combining pixel %d,%d: %x %x %f %f %f %f %d %d %d\n", x, y, badMaskBits, blankMaskBits, iter, rej, sys, olympic, useVariance, safe, rejection);
-
-/*          // Pre-reject inputs using KMM bimodality test. */
-/*        if (1)  { */
-/*          KMMRejectUnpopular(input,x,y); */
-/* /\*      rejection = true; *\/ */
-/*        } */
-/*        else { */
-/*          KMMRejectBright(input,x,y); */
-/*        } */
-
-#ifdef TESTING
-            if (PS_SQR(x - TEST_X) + PS_SQR(y - TEST_Y) <= PS_SQR(TEST_RADIUS)) {
-                fprintf (stderr, "combine\n");
-            }
-# endif
 
             psVector *reject = NULL; // Images to reject for this pixel