Changeset 41872 for branches/eam_branches/psModules.20211028/src/imcombine/pmStack.c

Timestamp:

Oct 29, 2021, 5:17:26 AM (5 years ago)

Author:

eugene

Message:

combine-by-percent must select a symmetric subset

File:

• branches/eam_branches/psModules.20211028/src/imcombine/pmStack.c (modified) (3 diffs)

branches/eam_branches/psModules.20211028/src/imcombine/pmStack.c

@@ -1514 +1514 @@
 
 # define MIN_GOOD_PERCENTILE 2
+# define SUSPECT_FRACTION 0.65
+
+// Comparison and swap functions for sorting values directly
+#define SORT_VV_COMPARE(A,B) (pixelData->data.F32[A] < pixelData->data.F32[B])
+#define SORT_VV_SWAP(TYPE,A,B) {                                        \
+    if (A != B) {                                                       \
+      psF32 tempVal = pixelData->data.F32[A];                           \
+      pixelData->data.F32[A] = pixelData->data.F32[B];                  \
+      pixelData->data.F32[B] = tempVal;                                 \
+      psF32 tempVar = pixelVariances->data.F32[A];                      \
+      pixelVariances->data.F32[A] = pixelVariances->data.F32[B];                \
+      pixelVariances->data.F32[B] = tempVar;                            \
+    }                                                                   \
+}
+
+#define SORT_VALUE_AND_VAR(A,B) { PSSORT(A->n, SORT_VV_COMPARE, SORT_VV_SWAP, F32); }
+  
+
 bool pmStackCombineByPercentile(
     pmReadout *combined,
     psArray *stackData,
-    psF64 minRange,
+    psF64 rejectFraction,
     psF64 maxRange,
     psImageMaskType badMaskBits,        // treat these bits as 'bad'
@@ -1541 +1559 @@
     psFree(stackReadouts);
 
-    psVector *pixelData = psVectorAlloc(stackData->n, PS_TYPE_F32);
-    psVector *pixelMask = psVectorAlloc(stackData->n, PS_TYPE_VECTOR_MASK);
+    // make sure the output readout matches the inputs, and set to blank by default
+    pmReadoutUpdateSize(combined, minInputCols, minInputRows, xSize, ySize, true, true, blankMaskBits);
+
+    psVector *pixelData      = psVectorAlloc(stackData->n, PS_TYPE_F32);
+    psVector *pixelVariances = psVectorAlloc(stackData->n, PS_TYPE_F32);
+    psVector *pixelMask      = psVectorAlloc(stackData->n, PS_TYPE_VECTOR_MASK);
+
+    int nGoodBits[16]; // accumulate the good pixel bits here for fuzzy logic
+    psAssert (sizeof(psImageMaskType) == 2, "psImageMaskType is not the expected size");
+    memset (nGoodBits, 0, 16*sizeof(int));
 
     for (int y = minInputRows; y < maxInputRows; y++) {
@@ -1552 +1578 @@
                 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;
+
+                psImage *image    = ro->image;
+                psImage *variance = ro->variance;
+                psImage *mask     = ro->mask;
 
                 int xIn = x - data->readout->col0;
                 int yIn = y - data->readout->row0; // Coordinates on input readout
 
+                // skip obviously bad input data
                 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)
+                if (fabs(image->data.F32[yIn][xIn]) > 1e5) continue; // XXX this cut is a bit arbitrary..
+                if (mask->data.PS_TYPE_IMAGE_MASK_DATA[yIn][xIn] & badMaskBits) continue;
+
                 // 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;
+                psImageMaskType value = 0x0001;
+                for (int nbit = 0; nbit < 16; nbit ++) {
+                  if (mask->data.PS_TYPE_IMAGE_MASK_DATA[yIn][xIn] & value) {
+                    nGoodBits[nbit] ++;
                   }
+                  value <<= 1;
                 }
-# endif
+
+                // pixelWeights->data.F32[nGood] = data->weight;
+                // pixelExps->data.F32[nGood] = data->exp;
+                // pixelSources->data.U16[nGood] = i;
 
                 pixelData->data.F32[nGood] = image->data.F32[yIn][xIn];
+                if (variance) {
+                  // pixelVariances->data.F32[numGood] = variance->data.F32[yIn][xIn] * addVariance->data.F32[i];
+                  pixelVariances->data.F32[nGood] = variance->data.F32[yIn][xIn];
+                }
+
                 nGood ++;
             }
             pixelData->n = nGood;
-            psVectorSortInPlace (pixelData);
-
-            if (nGood < MIN_GOOD_PERCENTILE) {
+
+
+            // rather than define a min and max value,
+            // what we really want is a symmetric selection about the middle,
+            // or the output is biased.  If N % 2 = 1, then 
+            
+            // we are going to exclude rejectFraction*nGood measurements.  But the
+            // rejection needs to be symmetric
+
+# define AT_LEAST 0
+# if (AT_LEAST)
+            int Ns = MIN(MAX(1, 0.5*rejectFraction * nGood), nGood);
+            int Ne = nGood - Ns;
+            int Npt = Ne - Ns;
+# else
+            int Ns = MIN(MAX(0, 0.5*rejectFraction * nGood), nGood);
+            int Ne = nGood - Ns;
+            int Npt = Ne - Ns;
+# endif
+            if ((Npt < 1) || (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;
+                combined->mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] = blankMaskBits; // probably not needed since it was set above.
                 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;
+            // sort both pixelData and pixelVariance
+            SORT_VALUE_AND_VAR (pixelData, pixelVariances);
+
+            // set the given (suspect) mask bit if nGoodBits[i] > f*nGood
+            // in other words if more than 65% of the good inputs had one of
+            // these bits set, then we should set that bit in the output mask
+            psImageMaskType value = 0x0001;
+            psImageMaskType outputMask = 0x0000;
+            for (int nbit = 0; nbit < 16; nbit ++) {
+              if (nGoodBits[nbit] > SUSPECT_FRACTION*nGood) {
+                outputMask |= 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;
+            float varSum = 0.0;
             for (int n = Ns; n < Ne; n++) {
                 sum += pixelData->data.F32[n];
+                varSum += pixelVariances->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);
-            }
+            // alternative: calculate the stdev of the pixel values
+            // 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)
+            float varValue = varSum / (Npt*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;
+            combined->mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] = outputMask;
         }
     }