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

Timestamp:

Nov 4, 2021, 6:05:18 PM (5 years ago)

Author:

eugene

Message:

merge changes from eam_branches/ipp-dev-20210817 (pmPattern updates, new inverse transform extra orders api, forward transform uses ORD, pmSourceIO_CMF.c.in conversion to pmFitsTableNew)

Location:

trunk/psModules

Files:

• . (modified) (1 prop)
• src/imcombine/pmStack.c (modified) (5 diffs)

trunk/psModules/src/imcombine/pmStack.c

@@ -756 +756 @@
         // *goodMask &= mask->data.PS_TYPE_IMAGE_MASK_DATA[yIn][xIn]; // save the mask bits still used
         // check for set bits and increment counter as appropriate
+        // count the number of times a given mask bit is set in the input pixels.
+        // NOTE: since we have explicitly skipped the pixels with any bad bits, these are only
+        // the suspect bits (nGoodBits is a bit of a misnomer: it is more like 'nSuspectBitsForGoodInputs'
         { 
             psImageMaskType value = 0x0001;
@@ -961 +964 @@
 
     // 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 +1513 @@
 }
 
+# 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;
+                psImage *mask = ro->mask;
+
+                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
+                if (mask->data.PS_TYPE_IMAGE_MASK_DATA[yIn][xIn] & badMaskBits);
+
+                // XXX save the count of suspect bits
+
+# if (0)
+                // count the number of times a given mask bit is set in the input pixels.
+                // NOTE: since we have explicitly skipped the pixels with any bad bits, these are only
+                // the suspect bits (nGoodBits is a bit of a misnomer: it is more like 'nSuspectBitsForGoodInputs'
+                if (1) { 
+                  psImageMaskType value = 0x0001;
+                  for (int nbit = 0; nbit < 16; nbit ++) {
+                    if (mask->data.PS_TYPE_IMAGE_MASK_DATA[yIn][xIn] & value) {
+                      // XXX nGoodBits[nbit] ++;
+                    }
+                    value <<= 1;
+                  }
+                }
+# endif
+
+                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;
+            }
+
+# if (0)
+# define SUSPECT_FRACTION 0.65
+            // set the mask bits if nGoodBits[i] > f*numGood
+            if (0) {
+              psImageMaskType value = 0x0001;
+              for (int nbit = 0; nbit < 16; nbit ++) {
+                if (nGoodBits[nbit] > SUSPECT_FRACTION*numGood) {
+                  *goodMask |= value;
+                }
+                value <<= 1;
+              }
+            }
+# endif
+
+            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;
+
+            // XXX this is probably not the correct variance to save
+            // the predicted variance should be the 1 / sum (1/var)
+
+            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 +1773 @@
         expweight = expmaps->variance;
     }
-
-/*     // Pre-reject inputs using KMM bimodality test. */
-/*     if (1)  { */
-/* /\*       KMMRejectUnpopular(input,x,y); *\/ */
-/*       rejection = true; */
-/*     } */
    
     // Set up rejection list
@@ -1660 +1783 @@
     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