Changeset 15913 for trunk/ppMerge/src/ppMergeMask.c

Timestamp:

Dec 24, 2007, 11:12:34 AM (18 years ago)

Author:

eugene

Message:

moved stats calculation out of pmMaskFlagSuspectPixels; added options for per-readout vs per-chip stats

File:

• trunk/ppMerge/src/ppMergeMask.c (modified) (2 diffs)

trunk/ppMerge/src/ppMergeMask.c

@@ -43 +43 @@
         pmChip *chipIn;                 // Input chip of interest
         while ((chipIn = pmFPAviewNextChip(view, fpaIn, 1))) {
+
+            // handle chip vs cell statistics & avoid reading the data twice
+
+            // load the data of all cells 
             pmCell *cellIn;             // Input cell of interest
             while ((cellIn = pmFPAviewNextCell(view, fpaIn, 1))) {
-                if (!pmCellRead(cellIn, fits, config->database)) {
-                    continue;
-                }
-                if (cellIn->readouts->n == 0) {
-                    continue;
-                }
-
-                pmCell *cellOut = pmFPAviewThisCell(view, fpaOut); // Output cell
-                // Suspect pixels image
-                psImage *suspect = psMetadataLookupPtr(NULL, cellOut->analysis, "MASK.SUSPECT");
-                bool first = suspect ? false : true;
-
-                pmReadout *roIn;        // Input readout of interest
-                while ((roIn = pmFPAviewNextReadout(view, fpaIn, 1))) {
-                    if (!roIn->image) {
-                        continue;
-                    }
-                    psTrace("ppMerge", 4, "Flagging suspect pixels in chip %d, cell %d, ro %d\n",
-                            view->chip, view->cell, view->readout);
-                    float frac = options->sample / (float)(roIn->image->numCols * roIn->image->numRows);
-                    suspect = pmMaskFlagSuspectPixels(suspect, roIn, options->maskSuspect,
-                                                      options->combine->maskVal, frac, rng);
-                }
-
-                if (first) {
-                    psMetadataAddImage(cellOut->analysis, PS_LIST_TAIL, "MASK.SUSPECT", 0,
-                                       "Suspect pixels", suspect);
-                    psFree(suspect);
-                }
-
+                if (!pmCellRead(cellIn, fits, config->database)) continue;
+                if (cellIn->readouts->n == 0) continue;
+            }
+
+            // calculate the readout statistics either for each readout, or across the entire chip
+            if (options->statsByChip) {
+                ppMergeMaskChipStats (chipIn, options, rng);
+            } else {
+                // calculate the stats for each cell independently
+                while ((cellIn = pmFPAviewNextCell(view, fpaIn, 1))) {
+                    if (cellIn->readouts->n == 0) continue;
+                    pmReadout *roIn;        // Input readout of interest
+                    while ((roIn = pmFPAviewNextReadout(view, fpaIn, 1))) {
+                        if (!roIn->image) continue;
+                        psTrace("ppMerge", 4, "Measure statistics for chip %d, cell %d, ro %d\n",
+                                view->chip, view->cell, view->readout);
+                        ppMergeMaskReadoutStats (roIn, options, rng);
+                    }
+                }
+            }
+
+            // apply the measured statistics to determine the outliers to be masked
+            while ((cellIn = pmFPAviewNextCell(view, fpaIn, 1))) {
+                if (cellIn->readouts->n == 0) continue;
+
+                pmCell *cellOut = pmFPAviewThisCell(view, fpaOut); // Output cell
+                // Suspect pixels image
+                psImage *suspect = psMetadataLookupPtr(NULL, cellOut->analysis, "MASK.SUSPECT");
+                bool first = suspect ? false : true;
+
+                pmReadout *roIn;        // Input readout of interest
+                while ((roIn = pmFPAviewNextReadout(view, fpaIn, 1))) {
+                    if (!roIn->image) {
+                        continue;
+                    }
+                    psTrace("ppMerge", 4, "Flagging suspect pixels in chip %d, cell %d, ro %d\n",
+                            view->chip, view->cell, view->readout);
+                    suspect = pmMaskFlagSuspectPixels(suspect, roIn, options->maskSuspect, options->combine->maskVal);
+                }
+
+                if (first) {
+                    psMetadataAddImage(cellOut->analysis, PS_LIST_TAIL, "MASK.SUSPECT", 0,
+                                       "Suspect pixels", suspect);
+                    psFree(suspect);
+                }
                 pmCellFreeData(cellIn);
             }
@@ -195 +214 @@
     return true;
 }
+
+bool ppMergeMaskReadoutStats(const pmReadout *readout, 
+                             ppMergeOptions *options, // Options
+                             psRandom *rng)
+{
+    PS_ASSERT_PTR_NON_NULL(readout, false);
+    PS_ASSERT_IMAGE_NON_NULL(readout->image, false);
+    PS_ASSERT_IMAGE_NON_EMPTY(readout->image, false);
+    PS_ASSERT_IMAGE_TYPE(readout->image, PS_TYPE_F32, false);
+    if (readout->mask) {
+        PS_ASSERT_IMAGE_NON_EMPTY(readout->mask, false);
+        PS_ASSERT_IMAGES_SIZE_EQUAL(readout->image, readout->mask, false);
+        PS_ASSERT_IMAGE_TYPE(readout->mask, PS_TYPE_MASK, false);
+    }
+
+    psImage *image = readout->image;    // Image of interest
+    psImage *mask = readout->mask;      // Corresponding mask
+
+    if (rng) {
+        psMemIncrRefCounter(rng);
+    } else {
+        rng = psRandomAlloc(PS_RANDOM_TAUS, 0);
+    }
+
+    // XXX note that this now will accept any of several stats options
+    psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV);
+    stats->nSubsample = options->sample;
+
+    if (!psImageBackground(stats, NULL, image, mask, options->combine->maskVal, rng)) {
+        psError(PS_ERR_UNKNOWN, false, "Failure to measure image statistics.\n");
+        psFree(stats);
+        psFree(rng);
+        return false;
+    }
+    if (!isfinite(stats->robustMedian) || !isfinite(stats->robustUQ) || !isfinite(stats->robustLQ)) {
+        psError(PS_ERR_UNKNOWN, false, "invalide image statistics (nan).\n");
+        psFree(stats);
+        psFree(rng);
+        return false;
+    }
+
+    psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "READOUT.MEDIAN", PS_META_REPLACE, "image stats", stats->robustMedian);
+    psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "READOUT.STDEV",  PS_META_REPLACE, "image stats", stats->robustStdev);
+
+    psFree(rng);
+    return true;
+}
+
+bool ppMergeMaskChipStats (const pmChip *chip,
+                           ppMergeOptions *options,
+                           psRandom *rng)
+{
+    PS_ASSERT_PTR_NON_NULL(chip, false);
+
+    // XXX note that this now will accept any of several stats options
+    psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV);
+
+    if (rng) {
+        psMemIncrRefCounter(rng);
+    } else {
+        rng = psRandomAlloc(PS_RANDOM_TAUS, 0);
+    }
+
+    // accumulate a vector of data values using options->sample per readout
+    psVector *values = psVectorAllocEmpty(options->sample, PS_TYPE_F32); // Vector containing subsample
+
+    for (int nCell = 0; nCell < chip->cells->n; nCell++) {
+        pmCell *cell = chip->cells->data[nCell];
+        if (cell->readouts->n == 0) continue;
+
+        for (int nReadout = 0; nReadout < cell->readouts->n; nReadout++) {
+            pmReadout *readout = cell->readouts->data[nReadout];
+            if (!readout->image) continue;
+
+            psTrace("ppMerge", 4, "Measure statistics for cell %d, readout %d\n", nCell, nReadout);
+
+            psImage *image = readout->image;    // Image of interest
+            psImage *mask = readout->mask;      // Corresponding mask
+
+            // Size of image
+            long nx = image->numCols;
+            long ny = image->numRows;
+            const int Npixels = nx*ny;  // Total number of pixels
+            const int Nsubset = PS_MIN(options->sample, Npixels); // Number of pixels in subset
+            
+            values = psVectorRealloc (values, values->n + Nsubset); // make sure we have enough space
+
+            // select a subset of the image pixels to measure the stats
+            for (long i = 0; i < Nsubset; i++) {
+                double frnd = psRandomUniform(rng);
+                int pixel = Npixels * frnd;
+                int ix = pixel % nx;
+                int iy = pixel / nx;
+
+                if (!isfinite(image->data.F32[iy][ix])) continue;
+                if (mask && (mask->data.U8[iy][ix] & options->combine->maskVal)) continue;
+
+                float value = image->data.F32[iy][ix];
+                values->data.F32[values->n] = value;
+                values->n ++;
+            }
+        }
+    }
+
+    // calculate the statistics
+    if (!psVectorStats (stats, values, NULL, NULL, 0)) {
+        psError(PS_ERR_UNKNOWN, false, "Unable to measure statistics for chip");
+        psFree(values);
+        psFree(stats);
+        psFree(rng);
+        return false;
+    }
+    if (!isfinite(stats->robustMedian) || !isfinite(stats->robustUQ) || !isfinite(stats->robustLQ)) {
+        psError(PS_ERR_UNKNOWN, false, "invalide image statistics (nan).\n");
+        psFree(values);
+        psFree(stats);
+        psFree(rng);
+        return false;
+    }
+
+    // supply the stats to the readout analysis metadata
+    for (int nCell = 0; nCell < chip->cells->n; nCell++) {
+        pmCell *cell = chip->cells->data[nCell];
+        if (cell->readouts->n == 0) continue;
+
+        for (int nReadout = 0; nReadout < cell->readouts->n; nReadout++) {
+            pmReadout *readout = cell->readouts->data[nReadout];
+            if (!readout->image) continue;
+
+            psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "READOUT.MEDIAN", PS_META_REPLACE, "image stats", stats->robustMedian);
+            psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "READOUT.STDEV",  PS_META_REPLACE, "image stats", stats->robustStdev);
+        }
+    }
+
+    psFree(values);
+    psFree(stats);
+    psFree(rng);
+    return true;
+}