IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 20444


Ignore:
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.

Location:
trunk/ppImage/src
Files:
5 edited

Legend:

Unmodified
Added
Removed

  • trunk/ppImage/src/ppImage.h

    r20419 r20444  
    3535    bool doFringe;                      // Fringe subtraction
    3636    bool doPhotom;                      // Source identification and photometry
     37    bool doBG;                          // Background subtraction
    3738    bool doAstromChip;                  // per-chip Astrometry
    3839    bool doAstromMosaic;                // full-mosaic Astrometry
    3940    bool doStats;                       // call ppStats on the image
    40     bool replaceMasked;                 // fill in masked values with background model
    4141
    4242    // output files requested
    … …  
    134134
    135135bool ppImageRebinChip (pmConfig *config, pmFPAview *view, ppImageOptions *options, char *outName);
    136 bool ppImageRebinReadout (pmReadout *output, pmReadout *input, pmFPAfile *outFile, ppImageOptions *options);
    137136
    138137bool ppImagePhotom (pmConfig *config, pmFPAview *view);
    … …  
    140139bool ppImageAddstar (pmConfig *config);
    141140
    142 bool ppImageReplaceBackground (pmConfig *config, pmFPAview *view, ppImageOptions *options);
     141// Subtract background from the chip-mosaicked image
     142bool ppImageSubtractBackground(
     143    pmConfig *config,                   // Configuration
     144    const pmFPAview *view,              // View to chip of interest
     145    const ppImageOptions *options       // Processing options
     146    );
    143147
    144148bool ppImageMosaicChip (pmConfig *config, const ppImageOptions *options, const pmFPAview *view,

  • trunk/ppImage/src/ppImageLoop.c

    r20419 r20444  
    129129
    130130        // replace the masked pixels with the background level
    131         if (options->replaceMasked) {
    132             if (!ppImageReplaceBackground(config, view, options)) {
    133                 ESCAPE("Unable to replace masked pixels with background level");
     131        if (options->doBG) {
     132            if (!ppImageSubtractBackground(config, view, options)) {
     133                ESCAPE("Unable to subtract background");
    134134            }
    135135        }

  • trunk/ppImage/src/ppImageOptions.c

    r18556 r20444  
    247247    options->doAstromChip   = psMetadataLookupBool(NULL, recipe, "ASTROM.CHIP");
    248248    options->doAstromMosaic = psMetadataLookupBool(NULL, recipe, "ASTROM.MOSAIC");
     249    options->doBG           = psMetadataLookupBool(NULL, recipe, "BACKGROUND");
    249250
    250251    // even if not requested explicitly, if any of these are set, build an internal mask and weight:
    … …  
    253254        options->doMaskBuild = true;
    254255        options->doWeightBuild = true;
    255     } else if (options->doMask) {
     256    } else if (options->doMask || options->doBG) {
    256257        options->doMaskBuild = true;
    257258    }
    … …  
    267268    options->fringeKeep = psMetadataLookupF32(NULL, recipe, "FRINGE.KEEP");
    268269
    269     options->replaceMasked  = psMetadataLookupBool(NULL, recipe, "REPLACE.MASKED");
    270270
    271271    return options;

  • trunk/ppImage/src/ppImageParseCamera.c

    r18896 r20444  
    7474            return NULL;
    7575        }
    76         // XXX have ppImageDefineFile return the pmFPAfile?
    77         pmFPAfile *mask = psMetadataLookupPtr(&status, config->files, "PPIMAGE.MASK");
    78         psAssert (mask, "mask not defined?  not possible!");
    79 
    80         // Need to read the names of bit masks from the mask header and set them in the
    81         // recipe.  If we are loading this from the detrend db, this action will happen
    82         // when the file is resolved.
    83         if (!mask->detrend) {
    84             // XXX need to load the mask bit names from one of the headers
    85             // this grabs the first available hdu : no guarantee that it will be valid, though
    86             pmHDU *hdu = pmHDUGetFirst (mask->fpa);
    87             if (!hdu) {
    88                 psError(PS_ERR_IO, true, "no valid HDU for PPIMAGE.INPUT.MASK");
    89                 return NULL;
    90             }
    91             // XXX is this consistent with the pmConfigMaskReadHeader call above?
    92             if (!pmConfigMaskReadHeader (config, hdu->header)) {
    93                 psError(PS_ERR_IO, false, "error in mask bits");
    94                 return NULL;
    95             }
    96         }
     76        // XXX have ppImageDefineFile return the pmFPAfile?
     77        pmFPAfile *mask = psMetadataLookupPtr(&status, config->files, "PPIMAGE.MASK");
     78        psAssert (mask, "mask not defined?  not possible!");
     79
     80        // Need to read the names of bit masks from the mask header and set them in the
     81        // recipe.  If we are loading this from the detrend db, this action will happen
     82        // when the file is resolved.
     83        if (!mask->detrend) {
     84            // XXX need to load the mask bit names from one of the headers
     85            // this grabs the first available hdu : no guarantee that it will be valid, though
     86            pmHDU *hdu = pmHDUGetFirst (mask->fpa);
     87            if (!hdu) {
     88                psError(PS_ERR_IO, true, "no valid HDU for PPIMAGE.INPUT.MASK");
     89                return NULL;
     90            }
     91            // XXX is this consistent with the pmConfigMaskReadHeader call above?
     92            if (!pmConfigMaskReadHeader (config, hdu->header)) {
     93                psError(PS_ERR_IO, false, "error in mask bits");
     94                return NULL;
     95            }
     96        }
    9797    }
    9898    if (options->doShutter) {
    … …  
    114114    int nThreads = psMetadataLookupS32 (&status, config->arguments, "NTHREADS");
    115115    if (nThreads > 0) {
    116         int nScanRows = psMetadataLookupS32 (&status, recipe, "SCAN.ROWS");
    117         pmDetrendSetThreadTasks (nScanRows);
     116        int nScanRows = psMetadataLookupS32 (&status, recipe, "SCAN.ROWS");
     117        pmDetrendSetThreadTasks (nScanRows);
    118118    }
    119119
    … …  
    274274    // For photometry, we operate on the chip-mosaicked image
    275275    // we create a copy of the mosaicked image for psphot so we can write out a clean image
    276     if (options->doPhotom) {
     276    if (options->doPhotom || options->doBG) {
    277277
    278278        // this file is just used as a carrier; output files (eg, PSPHOT.RESID) are defined by

  • trunk/ppImage/src/ppImageReplaceBackground.c

    r18556 r20444  
    33#endif
    44
    5 # include "ppImage.h"
    6 
    7 enum {MODE_NONE, MODE_MODEL, MODE_VALUE};
     5#include "ppImage.h"
    86
    97// In this function, we perform the psphot analysis routine for the chip-mosaicked images
    10 bool ppImageReplaceBackground (pmConfig *config, pmFPAview *view, ppImageOptions *options) {
     8bool ppImageSubtractBackground(pmConfig *config, const pmFPAview *view, const ppImageOptions *options)
     9{
     10    psAssert(config, "Need configuration");
     11    psAssert(view, "Need view to chip");
     12    psAssert(options, "Need options");
    1113
    12     bool status;
    13     pmCell *cell;
    14     pmReadout *readout;
     14    if (!options->doBG) {
     15        return true;
     16    }
    1517
    16     // find the reference source image
    17     pmFPAfile *input = psMetadataLookupPtr (&status, config->files, "PPIMAGE.CHIP");
     18    bool status;                        // Status of MD lookup
     19    pmFPAfile *input = psMetadataLookupPtr(&status, config->files, "PPIMAGE.CHIP"); // File to correct
    1820    if (!status) {
    19         psError(PSPHOT_ERR_CONFIG, false, "PSPHOT.INPUT I/O file is not defined");
     21        psError(PS_ERR_UNEXPECTED_NULL, true, "PPIMAGE.CHIP file is not defined");
    2022        return false;
    2123    }
    2224
    23     // select the appropriate recipe information
    24     psMetadata *ppImageRecipe  = psMetadataLookupPtr (&status, config->recipes, "PPIMAGE");
    25 
    26     // select the appropriate recipe information (needed by psphotModelBackground)
    27     psMetadata *psphotRecipe  = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE);
    28 
    29     char *replaceSkyModeWord = psMetadataLookupStr (&status, ppImageRecipe, "REPLACE.MODE");
    30     float replaceSkyValue = 0.0;
    31 
    32     int replaceSkyMode = MODE_NONE;
    33     if (!strcasecmp (replaceSkyModeWord, "MODEL")) {
    34         replaceSkyMode = MODE_MODEL;
    35     }
    36     if (!strcasecmp (replaceSkyModeWord, "VALUE")) {
    37         replaceSkyMode = MODE_VALUE;
    38         replaceSkyValue = psMetadataLookupF32 (&status, ppImageRecipe, "REPLACE.VALUE");
    39     }
    40     psAssert (replaceSkyMode != MODE_NONE, "replace sky mode not defined");
    41 
    42     // XXX Should this be options->maskValue or & ~options->satMask? the latter will leave saturated pixels high
    43     // psMaskType maskVal  = options->maskValue & ~options->satMask;
    44     psMaskType maskVal  = options->maskValue;
    45 
    46     pmFPAfile *modelFile = NULL;
    47     if (replaceSkyMode == MODE_MODEL) {
    48         // find the model data, if it exists yet (if not, we build it below)
    49         modelFile = psMetadataLookupPtr(&status, config->files, "PSPHOT.BACKMDL");
    50 
    51         // user-defined masks to test for good/bad pixels (build from recipe list if not yet set)
    52         psMetadataAddU8 (psphotRecipe, PS_LIST_TAIL, "MASK.PSPHOT", PS_META_REPLACE, "user-defined mask", maskVal);
     25    // The background model file may be defined, though the model hasn't been built
     26    // If so, we will build it below.
     27    pmFPAfile *modelFile = psMetadataLookupPtr(&status, config->files, "PSPHOT.BACKMDL"); // Background model
     28    if (!status) {
     29        psError(PS_ERR_UNEXPECTED_NULL, true, "PSPHOT.BACKMDL file is not defined");
     30        return false;
    5331    }
    5432
    55     // iterate over the cells and readout for this chip
    56     while ((cell = pmFPAviewNextCell (view, input->fpa, 1)) != NULL) {
    57         psLogMsg ("ppImagePhotom", 5, "Cell %d: %x %x\n", view->cell, cell->file_exists, cell->process);
    58         if (! cell->process || ! cell->file_exists) { continue; }
     33    psMetadata *ppImageRecipe = psMetadataLookupPtr(NULL, config->recipes, RECIPE_NAME);
     34    psAssert(ppImageRecipe, "Need PPIMAGE recipe");
     35    psMetadata *psphotRecipe = psMetadataLookupPtr(NULL, config->recipes, PSPHOT_RECIPE);
     36    psAssert(psphotRecipe, "Need PSPHOT recipe");
    5937
    60         // process each of the readouts
    61         while ((readout = pmFPAviewNextReadout (view, input->fpa, 1)) != NULL) {
    62             if (! readout->data_exists) { continue; }
    63             if (! readout->mask) { continue; }
     38    // XXX Should this be options->maskValue or options->maskValue & ~options->satMask?
     39    //     The latter will leave saturated pixels high
     40    psMaskType maskVal = options->maskValue;
    6441
    65             // replace masked pixels with values from model (unbinning not needed)
    66             pmReadout *modelReadout = NULL;
    67             psImageBinning *binning = NULL;
    68             psImage *model = NULL;
    69             if (replaceSkyMode == MODE_MODEL) {
    70                 // we are using this pmFPAfile as an I/O file: select readout or create
    71                 if (modelFile) {
    72                     modelReadout = pmFPAviewThisReadout (view, modelFile->fpa);
    73                 }
    74                 if (!modelReadout) {
    75                     psphotModelBackground (config, view, "PPIMAGE.CHIP");
    76                     modelFile = psMetadataLookupPtr(&status, config->files, "PSPHOT.BACKMDL");
    77                     assert (modelFile);
    78                     if (modelFile->mode == PM_FPA_MODE_INTERNAL) {
    79                         modelReadout = modelFile->readout;
    80                     } else {
    81                         modelReadout = pmFPAviewThisReadout (view, modelFile->fpa);
    82                     }
    83                     assert (modelReadout);
    84                 }
    85                 binning = psMetadataLookupPtr(&status, psphotRecipe, "PSPHOT.BACKGROUND.BINNING");
    86                 model = modelReadout->image;
    87             }
     42    // user-defined masks to test for good/bad pixels (build from recipe list if not yet set)
     43    psMetadataAddU8(psphotRecipe, PS_LIST_TAIL, "MASK.PSPHOT", PS_META_REPLACE, "user-defined mask", maskVal);
    8844
    89             psImage *mask = readout->mask;
    90             psImage *image = readout->image;
     45    // Since we are working on a chip-mosaicked image, there should only be a single cell and readout
     46    pmChip *chip = pmFPAviewThisChip(view, input->fpa); // Chip of interest
     47    if (!chip) {
     48        psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to find chip");
     49        return false;
     50    }
     51    if (chip->cells->n == 0) {
     52        psWarning("Chip has no cells");
     53        return true;
     54    }
     55    if (chip->cells->n > 1) {
     56        psWarning("Chip has %ld cells; only the first will be processed", chip->cells->n);
     57    }
     58    pmCell *cell = chip->cells->data[0]; // Cell of interest
     59    if (!cell || !cell->process || !cell->file_exists) {
     60        // Nothing to process
     61        return true;
     62    }
     63    if (cell->readouts->n == 0) {
     64        psWarning("Cell has no readouts");
     65        return true;
     66    }
     67    if (cell->readouts->n > 1) {
     68        psWarning("Cell has %ld readouts; only the first will be processed", cell->readouts->n);
     69    }
     70    pmReadout *ro = cell->readouts->data[0]; // Readout of interest
     71    if (!ro || !ro->data_exists) {
     72        // Nothing to process
     73        return true;
     74    }
     75    psImage *image = ro->image, *mask = ro->mask; // Image and mask of interest
    9176
    92             for (int iy = 0; iy < image->numRows; iy++) {
    93                 for (int ix = 0; ix < image->numCols; ix++) {
    94                     if (!(mask->data.U8[iy][ix] && maskVal)) continue;
    95                     if (replaceSkyMode == MODE_MODEL) {
    96                         image->data.F32[iy][ix] = psImageUnbinPixel_V2(ix, iy, model, binning);
    97                     } else {
    98                         image->data.F32[iy][ix] = replaceSkyValue;
    99                     }
    100                 }
    101             }
    102         }
     77
     78    pmFPAview roView = *view;           // View to readout
     79    roView.cell = roView.readout = 0;
     80    pmReadout *modelRO = pmFPAviewThisReadout(&roView, modelFile->fpa); // Background model
     81    if (!modelRO) {
     82        // Create the background model
     83        if (!psphotModelBackground(config, &roView, "PPIMAGE.CHIP")) {
     84            psError(PS_ERR_UNKNOWN, false, "Unable to model background");
     85            return false;
     86        }
     87        modelRO = (modelFile->mode == PM_FPA_MODE_INTERNAL) ? modelFile->readout :
     88                   pmFPAviewThisReadout(&roView, modelFile->fpa);
     89        if (!modelRO) {
     90            psError(PS_ERR_UNEXPECTED_NULL, true, "Unable to find background model");
     91            return false;
     92        }
     93    }
     94    psImageBinning *binning = psMetadataLookupPtr(&status, psphotRecipe,
     95                                                  "PSPHOT.BACKGROUND.BINNING"); // Binning for model
     96    psImage *modelImage = modelRO->image; // Background model
     97
     98    // Do the background subtraction
     99    int numCols = image->numCols, numRows = image->numRows; // Size of image
     100    for (int y = 0; y < numRows; y++) {
     101        for (int x = 0; x < numCols; x++) {
     102            if (mask && mask->data.PS_TYPE_MASK_DATA[y][x] & maskVal) {
     103                image->data.F32[y][x] = 0.0;
     104            } else {
     105                float value = psImageUnbinPixel_V2(x, y, modelImage, binning); // Background value
     106                if (!isfinite(value)) {
     107                    image->data.F32[y][x] = NAN;
     108                    mask->data.PS_TYPE_MASK_DATA[y][x] |= options->badMask;
     109                } else {
     110                    image->data.F32[y][x] -= value;
     111                }
     112            }
     113        }
    103114    }
    104115
    105     if (replaceSkyMode == MODE_MODEL) {
    106         pmFPAfileDropInternal (config->files, "PSPHOT.BACKMDL");
    107         pmFPAfileDropInternal (config->files, "PSPHOT.BACKMDL.STDEV");
    108     }
     116    pmFPAfileDropInternal(config->files, "PSPHOT.BACKMDL");
     117    pmFPAfileDropInternal(config->files, "PSPHOT.BACKMDL.STDEV");
    109118
    110119    return true;
Note: See TracChangeset for help on using the changeset viewer.