Changeset 20444 for trunk/ppImage/src/ppImageReplaceBackground.c

Timestamp:

Oct 28, 2008, 1:49:32 PM (18 years ago)

Author:

Paul Price

Message:

Implementing background subtraction in the place of replacing bad pixels with the background. Bad pixels are automatically identified and set to 0 in the process, so this is effectively the same thing as before, except it subtracts the background as well. The background (from psphot) is built if it doesn't exist.

File:

• trunk/ppImage/src/ppImageReplaceBackground.c (modified) (1 diff)

trunk/ppImage/src/ppImageReplaceBackground.c

@@ -3 +3 @@
 #endif
 
-# include "ppImage.h"
-
-enum {MODE_NONE, MODE_MODEL, MODE_VALUE};
+#include "ppImage.h"
 
 // In this function, we perform the psphot analysis routine for the chip-mosaicked images
-bool ppImageReplaceBackground (pmConfig *config, pmFPAview *view, ppImageOptions *options) {
+bool ppImageSubtractBackground(pmConfig *config, const pmFPAview *view, const ppImageOptions *options)
+{
+    psAssert(config, "Need configuration");
+    psAssert(view, "Need view to chip");
+    psAssert(options, "Need options");
 
-    bool status;
-    pmCell *cell;
-    pmReadout *readout;
+    if (!options->doBG) {
+        return true;
+    }
 
-    // find the reference source image
-    pmFPAfile *input = psMetadataLookupPtr (&status, config->files, "PPIMAGE.CHIP");
+    bool status;                        // Status of MD lookup
+    pmFPAfile *input = psMetadataLookupPtr(&status, config->files, "PPIMAGE.CHIP"); // File to correct
     if (!status) {
-        psError(PSPHOT_ERR_CONFIG, false, "PSPHOT.INPUT I/O file is not defined");
+        psError(PS_ERR_UNEXPECTED_NULL, true, "PPIMAGE.CHIP file is not defined");
         return false;
     }
 
-    // select the appropriate recipe information
-    psMetadata *ppImageRecipe  = psMetadataLookupPtr (&status, config->recipes, "PPIMAGE");
-
-    // select the appropriate recipe information (needed by psphotModelBackground)
-    psMetadata *psphotRecipe  = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE);
-
-    char *replaceSkyModeWord = psMetadataLookupStr (&status, ppImageRecipe, "REPLACE.MODE");
-    float replaceSkyValue = 0.0;
-
-    int replaceSkyMode = MODE_NONE;
-    if (!strcasecmp (replaceSkyModeWord, "MODEL")) {
-        replaceSkyMode = MODE_MODEL;
-    }
-    if (!strcasecmp (replaceSkyModeWord, "VALUE")) {
-        replaceSkyMode = MODE_VALUE;
-        replaceSkyValue = psMetadataLookupF32 (&status, ppImageRecipe, "REPLACE.VALUE");
-    }
-    psAssert (replaceSkyMode != MODE_NONE, "replace sky mode not defined");
-
-    // XXX Should this be options->maskValue or & ~options->satMask? the latter will leave saturated pixels high
-    // psMaskType maskVal  = options->maskValue & ~options->satMask;
-    psMaskType maskVal  = options->maskValue;
-
-    pmFPAfile *modelFile = NULL;
-    if (replaceSkyMode == MODE_MODEL) {
-        // find the model data, if it exists yet (if not, we build it below)
-        modelFile = psMetadataLookupPtr(&status, config->files, "PSPHOT.BACKMDL");
-
-        // user-defined masks to test for good/bad pixels (build from recipe list if not yet set)
-        psMetadataAddU8 (psphotRecipe, PS_LIST_TAIL, "MASK.PSPHOT", PS_META_REPLACE, "user-defined mask", maskVal);
+    // The background model file may be defined, though the model hasn't been built
+    // If so, we will build it below.
+    pmFPAfile *modelFile = psMetadataLookupPtr(&status, config->files, "PSPHOT.BACKMDL"); // Background model
+    if (!status) {
+        psError(PS_ERR_UNEXPECTED_NULL, true, "PSPHOT.BACKMDL file is not defined");
+        return false;
     }
 
-    // iterate over the cells and readout for this chip
-    while ((cell = pmFPAviewNextCell (view, input->fpa, 1)) != NULL) {
-        psLogMsg ("ppImagePhotom", 5, "Cell %d: %x %x\n", view->cell, cell->file_exists, cell->process);
-        if (! cell->process || ! cell->file_exists) { continue; }
+    psMetadata *ppImageRecipe = psMetadataLookupPtr(NULL, config->recipes, RECIPE_NAME);
+    psAssert(ppImageRecipe, "Need PPIMAGE recipe");
+    psMetadata *psphotRecipe = psMetadataLookupPtr(NULL, config->recipes, PSPHOT_RECIPE);
+    psAssert(psphotRecipe, "Need PSPHOT recipe");
 
-        // process each of the readouts
-        while ((readout = pmFPAviewNextReadout (view, input->fpa, 1)) != NULL) {
-            if (! readout->data_exists) { continue; }
-            if (! readout->mask) { continue; }
+    // XXX Should this be options->maskValue or options->maskValue & ~options->satMask?
+    //     The latter will leave saturated pixels high
+    psMaskType maskVal = options->maskValue;
 
-            // replace masked pixels with values from model (unbinning not needed)
-            pmReadout *modelReadout = NULL;
-            psImageBinning *binning = NULL;
-            psImage *model = NULL;
-            if (replaceSkyMode == MODE_MODEL) {
-                // we are using this pmFPAfile as an I/O file: select readout or create
-                if (modelFile) {
-                    modelReadout = pmFPAviewThisReadout (view, modelFile->fpa);
-                }
-                if (!modelReadout) {
-                    psphotModelBackground (config, view, "PPIMAGE.CHIP");
-                    modelFile = psMetadataLookupPtr(&status, config->files, "PSPHOT.BACKMDL");
-                    assert (modelFile);
-                    if (modelFile->mode == PM_FPA_MODE_INTERNAL) {
-                        modelReadout = modelFile->readout;
-                    } else {
-                        modelReadout = pmFPAviewThisReadout (view, modelFile->fpa);
-                    }
-                    assert (modelReadout);
-                }
-                binning = psMetadataLookupPtr(&status, psphotRecipe, "PSPHOT.BACKGROUND.BINNING");
-                model = modelReadout->image;
-            }
+    // user-defined masks to test for good/bad pixels (build from recipe list if not yet set)
+    psMetadataAddU8(psphotRecipe, PS_LIST_TAIL, "MASK.PSPHOT", PS_META_REPLACE, "user-defined mask", maskVal);
 
-            psImage *mask = readout->mask;
-            psImage *image = readout->image;
+    // Since we are working on a chip-mosaicked image, there should only be a single cell and readout
+    pmChip *chip = pmFPAviewThisChip(view, input->fpa); // Chip of interest
+    if (!chip) {
+        psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to find chip");
+        return false;
+    }
+    if (chip->cells->n == 0) {
+        psWarning("Chip has no cells");
+        return true;
+    }
+    if (chip->cells->n > 1) {
+        psWarning("Chip has %ld cells; only the first will be processed", chip->cells->n);
+    }
+    pmCell *cell = chip->cells->data[0]; // Cell of interest
+    if (!cell || !cell->process || !cell->file_exists) {
+        // Nothing to process
+        return true;
+    }
+    if (cell->readouts->n == 0) {
+        psWarning("Cell has no readouts");
+        return true;
+    }
+    if (cell->readouts->n > 1) {
+        psWarning("Cell has %ld readouts; only the first will be processed", cell->readouts->n);
+    }
+    pmReadout *ro = cell->readouts->data[0]; // Readout of interest
+    if (!ro || !ro->data_exists) {
+        // Nothing to process
+        return true;
+    }
+    psImage *image = ro->image, *mask = ro->mask; // Image and mask of interest
 
-            for (int iy = 0; iy < image->numRows; iy++) {
-                for (int ix = 0; ix < image->numCols; ix++) {
-                    if (!(mask->data.U8[iy][ix] && maskVal)) continue;
-                    if (replaceSkyMode == MODE_MODEL) {
-                        image->data.F32[iy][ix] = psImageUnbinPixel_V2(ix, iy, model, binning);
-                    } else {
-                        image->data.F32[iy][ix] = replaceSkyValue;
-                    }
-                }
-            }
-        }
+
+    pmFPAview roView = *view;           // View to readout
+    roView.cell = roView.readout = 0;
+    pmReadout *modelRO = pmFPAviewThisReadout(&roView, modelFile->fpa); // Background model
+    if (!modelRO) {
+        // Create the background model
+        if (!psphotModelBackground(config, &roView, "PPIMAGE.CHIP")) {
+            psError(PS_ERR_UNKNOWN, false, "Unable to model background");
+            return false;
+        }
+        modelRO = (modelFile->mode == PM_FPA_MODE_INTERNAL) ? modelFile->readout :
+                   pmFPAviewThisReadout(&roView, modelFile->fpa);
+        if (!modelRO) {
+            psError(PS_ERR_UNEXPECTED_NULL, true, "Unable to find background model");
+            return false;
+        }
+    }
+    psImageBinning *binning = psMetadataLookupPtr(&status, psphotRecipe,
+                                                  "PSPHOT.BACKGROUND.BINNING"); // Binning for model
+    psImage *modelImage = modelRO->image; // Background model
+
+    // Do the background subtraction
+    int numCols = image->numCols, numRows = image->numRows; // Size of image
+    for (int y = 0; y < numRows; y++) {
+        for (int x = 0; x < numCols; x++) {
+            if (mask && mask->data.PS_TYPE_MASK_DATA[y][x] & maskVal) {
+                image->data.F32[y][x] = 0.0;
+            } else {
+                float value = psImageUnbinPixel_V2(x, y, modelImage, binning); // Background value
+                if (!isfinite(value)) {
+                    image->data.F32[y][x] = NAN;
+                    mask->data.PS_TYPE_MASK_DATA[y][x] |= options->badMask;
+                } else {
+                    image->data.F32[y][x] -= value;
+                }
+            }
+        }
     }
 
-    if (replaceSkyMode == MODE_MODEL) {
-        pmFPAfileDropInternal (config->files, "PSPHOT.BACKMDL");
-        pmFPAfileDropInternal (config->files, "PSPHOT.BACKMDL.STDEV");
-    }
+    pmFPAfileDropInternal(config->files, "PSPHOT.BACKMDL");
+    pmFPAfileDropInternal(config->files, "PSPHOT.BACKMDL.STDEV");
 
     return true;