Changeset 41868

Timestamp:

Oct 28, 2021, 7:54:55 AM (5 years ago)

Author:

eugene

Message:

updates from laptop: fix conflict

File:

• branches/eam_branches/psModules.20211015/src/imcombine/pmStack.c (modified) (1 diff)

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

@@ -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;
+                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) continue;
+
+                // 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);
+}
+
+# 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(