Index: trunk/ppMerge/src/Makefile.am =================================================================== --- trunk/ppMerge/src/Makefile.am (revision 9955) +++ trunk/ppMerge/src/Makefile.am (revision 9996) @@ -10,4 +10,5 @@ ppMergeConfig.c \ ppMergeData.c \ + ppMergeMask.c \ ppMergeOptions.c \ ppMergeScaleZero.c @@ -20,4 +21,5 @@ ppMergeConfig.h \ ppMergeData.h \ + ppMergeMask.h \ ppMergeOptions.h \ ppMergeScaleZero.h Index: trunk/ppMerge/src/ppMerge.c =================================================================== --- trunk/ppMerge/src/ppMerge.c (revision 9955) +++ trunk/ppMerge/src/ppMerge.c (revision 9996) @@ -14,4 +14,5 @@ #include "ppMergeCombine.h" #include "ppMergeScaleZero.h" +#include "ppMergeMask.h" //#include "ppMem.h" @@ -45,11 +46,17 @@ } - // Measure the background in each image psImage *scale = NULL; // The scalings psImage *zero = NULL; // The zeros - ppMergeScaleZero(&scale, &zero, data, options, config); - // Do the combination and write - ppMergeCombine(scale, zero, data, options, config); + if (options->mask) { + // Generate a mask + ppMergeMask(data, options, config); + } else { + // Measure the background in each image + ppMergeScaleZero(&scale, &zero, data, options, config); + + // Do the combination and write + ppMergeCombine(scale, zero, data, options, config); + } // Output the statistics Index: trunk/ppMerge/src/ppMergeCombine.c =================================================================== --- trunk/ppMerge/src/ppMergeCombine.c (revision 9955) +++ trunk/ppMerge/src/ppMergeCombine.c (revision 9996) @@ -78,31 +78,4 @@ #endif -#if 0 -// Write fringe extensions -static void writeFringes(psMetadata *fringes, // Table of extensions with fringe statistics to write - psFits *fits // FITS file to which to write - ) -{ - assert(fringes); - - // Write the fringe extensions - psMetadataIterator *fringesIter = psMetadataIteratorAlloc(fringes, PS_LIST_HEAD, false); // Iterator - psMetadataItem *fringeItem; // Item from iteration - while ((fringeItem = psMetadataGetAndIncrement(fringesIter))) { - const char *extname = fringeItem->name; // Extension name - psArray *fringe = fringeItem->data.V; // Fringe data (array of fringe components) - if (!pmFringesWriteFits(fits, NULL, fringe, extname)) { - psWarning("Unable to write fringe data to extension %s\n", extname); - } - } - - // Clear the list - while (psListLength(fringes->list) > 0) { - psMetadataRemoveIndex(fringes, PS_LIST_TAIL); - } - - return; -} -#endif // Combine the inputs @@ -133,4 +106,8 @@ pmFPAWrite(data->out, data->outFile, config->database, true, false); // Write header only pmChip *chip; // Chip of interest + psRandom *rng = NULL; // Random number generator; required for building a mask + if (options->mask) { + rng = psRandomAlloc(PS_RANDOM_TAUS, 0); + } while ((chip = pmFPAviewNextChip(view, fpa, 1))) { pmChipWrite(chip, data->outFile, config->database, true, false); // Write header only @@ -150,13 +127,10 @@ } - // Building a shutter correction - if (options->shutter) { - // Read in the pixels - int numRead = 0; // Number of inputs read - psVector *exptimes = psVectorAlloc(filenames->n, PS_TYPE_F32); // Reference measurements - psVectorInit(exptimes, 0); - psArray *images = psArrayAlloc(filenames->n); // Array of images - psArray *weights = psArrayAlloc(filenames->n); // Array of weights - psArray *masks = psArrayAlloc(filenames->n); // Array of masks + // Read bit by bit + int numRead; // Number of inputs read + int numScan = 0; + + do { + numRead = 0; for (int i = 0; i < filenames->n; i++) { if (! filenames->data[i] || strlen(filenames->data[i]) == 0) { @@ -167,126 +141,46 @@ continue; } - pmFPA *fpaIn = data->in->data[i]; // Input FPA - pmChip *chipIn = fpaIn->chips->data[view->chip]; // Input chip - pmCell *cellIn = chipIn->cells->data[view->cell]; // Input cell - if (!cellIn->file_exists) { - // || !cellIn->data_exists) { - continue; - } - if (!pmCellRead(cellIn, fits, config->database) || cellIn->readouts->n == 0) { - continue; - } - if (cellIn->readouts->n > 1) { - psWarning("Multiple readouts in file %d, chip %d, cell %d --- " - "only the first will be used.\n", i, view->chip, view->cell); - } - float exptime = psMetadataLookupF32(NULL, cellIn->concepts, "CELL.EXPOSURE"); - if (!isfinite(exptime)) { - psWarning("File %d, chip %d, cell %d has bad exposure time.\n", - i, view->chip, view->cell); - } - exptimes->data.F32[i] = exptime; - pmReadout *readoutIn = cellIn->readouts->data[0]; - pmReadoutSetWeight(readoutIn); - images->data[i] = psMemIncrRefCounter(readoutIn->image); - weights->data[i] = psMemIncrRefCounter(readoutIn->weight); - masks->data[i] = psMemIncrRefCounter(readoutIn->mask); - numRead++; + + if (!stack->data[i]) { + pmFPA *fpaIn = data->in->data[i]; // Input FPA + pmChip *chipIn = fpaIn->chips->data[view->chip]; // Input chip + pmCell *cellIn = chipIn->cells->data[view->cell]; // Input cell + stack->data[i] = pmReadoutAlloc(cellIn); // Input readout + } + + // Only reading and writing the first readout in each cell (plane 0) + if (pmReadoutReadNext(stack->data[i], fits, 0, options->rows)) { + // For anything except BIAS and DARK, we want a mask + if (options->scale || options->zero || options->shutter) { + pmReadoutSetMask(stack->data[i]); + } + if (options->shutter) { + pmReadoutSetWeight(stack->data[i]); + } + + numRead++; + } else { + psTrace("ppMerge", 3, "Unable to read from file %d for chip %d, " + "cell %d, scan %d\n", i, view->chip, view->cell, numScan); + } + } - if (numRead > 1) { - readout->image = pmShutterCorrectionMeasure(exptimes, images, weights, masks, - options->shutterSize, options->mean, - options->stdev, options->shutterIter, - options->shutterRej, - options->combine->maskVal); - psTrace("ppMerge", 5, "Chip %d, cell %d\n", view->chip, view->cell); - } else { - psError(PS_ERR_UNKNOWN, true, "Not enough images to combine: %d\n", numRead); - psFree(images); - psFree(exptimes); - return false; + + psTrace("ppMerge", 5, "Chip %d, cell %d, scan %d\n", view->chip, view->cell, numScan); + if (numRead > 0) { + + if (options->shutter) { + readout->image = pmShutterCorrectionMeasure(stack, options->shutterSize, + options->mean, options->stdev, + options->shutterIter, + options->shutterRej, + options->combine->maskVal); + } else { + pmReadoutCombine(readout, stack, cellZeros, cellScales, options->combine); + } } - psFree(images); - psFree(masks); - psFree(weights); - psFree(exptimes); - } else { -#if 1 - // Read bit by bit - - int numRead; // Number of inputs read - int numScan = 0; - do { - numRead = 0; - for (int i = 0; i < filenames->n; i++) { - if (! filenames->data[i] || strlen(filenames->data[i]) == 0) { - continue; - } - psFits *fits = data->files->data[i]; // FITS file handle - if (!fits) { - continue; - } - - if (!stack->data[i]) { - pmFPA *fpaIn = data->in->data[i]; // Input FPA - pmChip *chipIn = fpaIn->chips->data[view->chip]; // Input chip - pmCell *cellIn = chipIn->cells->data[view->cell]; // Input cell - stack->data[i] = pmReadoutAlloc(cellIn); // Input readout - } - - // Only reading and writing the first readout in each cell (plane 0) - if (pmReadoutReadNext(stack->data[i], fits, 0, options->rows)) { - // For anything except BIAS and DARK, we want a mask - if (options->scale || options->zero) { - pmReadoutSetMask(stack->data[i]); - } - numRead++; - } else { - psTrace("ppMerge", 3, "Unable to read from file %d for chip %d, " - "cell %d, scan %d\n", i, view->chip, view->cell, numScan); - } - } - - psTrace("ppMerge", 5, "Chip %d, cell %d, scan %d\n", view->chip, view->cell, numScan); - if (numRead > 0) { - pmReadoutCombine(readout, stack, cellZeros, cellScales, options->combine); - } - numScan++; - - } while (numRead > 0); - -#else - // Read whole cells at a time - - printf("Chip %d, cell %d\n", view->chip, view->cell); - int numRead = 0; // Number of inputs read - for (int i = 0; i < filenames->n; i++) { - if (! filenames->data[i] || strlen(filenames->data[i]) == 0) { - continue; - } - psFits *fits = data->files->data[i]; // FITS file handle - if (!fits) { - continue; - } - pmFPA *fpaIn = data->in->data[i]; // Input FPA - pmChip *chipIn = fpaIn->chips->data[view->chip]; // Input chip - pmCell *cellIn = chipIn->cells->data[view->cell]; // Input cell - if (!cellIn->file_exists || !cellIn->data_exists) { - continue; - } - if (!pmCellRead(cellIn, fits, config->database) || cellIn->readouts->n == 0) { - continue; - } - stack->data[i] = cellIn->readouts->data[0]; - numRead++; - } - if (numRead > 0) { - pmReadoutCombine(readout, stack, cellZeros, cellScales, options->combine); - psTrace("ppMerge", 5, "Chip %d, cell %d\n", view->chip, view->cell); - } - - memCheck(); -#endif - } + numScan++; + + } while (numRead > 0); psFree(stack); @@ -392,4 +286,5 @@ psFree(view); + psFree(rng); return true; Index: trunk/ppMerge/src/ppMergeMask.c =================================================================== --- trunk/ppMerge/src/ppMergeMask.c (revision 9996) +++ trunk/ppMerge/src/ppMergeMask.c (revision 9996) @@ -0,0 +1,126 @@ +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include +#include +#include +#include +#include + +#include "ppMerge.h" +#include "ppMergeData.h" +#include "ppMergeMask.h" + + +// Generate a mask +bool ppMergeMask(ppMergeData *data, // Data + ppMergeOptions *options, // Options + pmConfig *config // Configuration + ) +{ + psArray *filenames = psMetadataLookupPtr(NULL, config->arguments, "INPUT"); // The input file names + assert(filenames); // It should be here --- it's put here in ppMergeConfig + + psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS, 0); // Random number generator + pmFPA *fpaOut = data->out; // Output FPA + + // Iterate over each file + for (int i = 0; i < filenames->n; i++) { + if (! filenames->data[i] || strlen(filenames->data[i]) == 0) { + continue; + } + psFits *fits = data->files->data[i]; // FITS file handle + if (!fits) { + continue; + } + psTrace("ppMerge", 3, "File %d: %s\n", i, (const char*)filenames->data[i]); + + pmFPA *fpaIn = data->in->data[i]; // Input FPA + pmFPAview *view = pmFPAviewAlloc(0); // View of FPA, for iteration + pmChip *chipIn; // Input chip of interest + while ((chipIn = pmFPAviewNextChip(view, fpaIn, 1))) { + 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); + } + + pmCellFreeData(cellIn); + } + pmChipFreeData(chipIn); + } + pmFPAFreeData(fpaIn); + psFree(view); + } + psFree(rng); + + pmFPAview *view = pmFPAviewAlloc(0);// View of FPA, for iteration + pmFPAWrite(fpaOut, data->outFile, config->database, true, false); // Write header only + pmChip *chipOut; // Output chip of interest + while ((chipOut = pmFPAviewNextChip(view, fpaOut, 1))) { + pmChipWrite(chipOut, data->outFile, config->database, true, false); // Write header only + pmCell *cellOut; // Output cell of interest + while ((cellOut = pmFPAviewNextCell(view, fpaOut, 1))) { + pmCellWrite(cellOut, data->outFile, config->database, true); // Write header only + + psImage *suspect = psMetadataLookupPtr(NULL, cellOut->analysis, "MASK.SUSPECT"); + if (! suspect) { + continue; + } + + pmReadout *roOut = pmReadoutAlloc(cellOut); // Output readout + roOut->mask = pmMaskIdentifyBadPixels(suspect, options->maskBad, options->combine->maskVal); + psFree(roOut); // Drop reference + + if (cellOut->hdu && !cellOut->hdu->blankPHU) { + psTrace("ppMerge", 5, "Writing out cell HDU.\n"); + pmCellWriteMask(cellOut, data->outFile, config->database); + pmCellFreeData(cellOut); + } + } + + if (chipOut->hdu && !chipOut->hdu->blankPHU) { + psTrace("ppMerge", 5, "Writing out chip HDU.\n"); + pmChipWriteMask(chipOut, data->outFile, config->database); + pmChipFreeData(chipOut); + } + } + + if (fpaOut->hdu && !fpaOut->hdu->blankPHU) { + psTrace("ppMerge", 5, "Writing out FPA HDU.\n"); + pmFPAWriteMask(fpaOut, data->outFile, config->database); + } + pmFPAFreeData(fpaOut); + + psFree(view); + + return true; +} Index: trunk/ppMerge/src/ppMergeMask.h =================================================================== --- trunk/ppMerge/src/ppMergeMask.h (revision 9996) +++ trunk/ppMerge/src/ppMergeMask.h (revision 9996) @@ -0,0 +1,15 @@ +#ifndef PP_MERGE_MASK +#define PP_MERGE_MASK + +#include +#include "ppMergeData.h" +#include "ppMergeOptions.h" + +// Generate a mask +bool ppMergeMask(ppMergeData *data, // Data + ppMergeOptions *options, // Options + pmConfig *config // Configuration + ); + + +#endif Index: trunk/ppMerge/src/ppMergeOptions.c =================================================================== --- trunk/ppMerge/src/ppMergeOptions.c (revision 9955) +++ trunk/ppMerge/src/ppMergeOptions.c (revision 9996) @@ -37,4 +37,5 @@ options->fringe = false; options->shutter = false; + options->mask = false; options->sample = 1; options->mean = PS_STAT_SAMPLE_MEDIAN; @@ -47,4 +48,6 @@ options->shutterIter = 2; options->shutterRej = 3.0; + options->maskSuspect = 5.0; + options->maskBad = 10.0; options->onOff = 0; options->combine = pmCombineParamsAlloc(PS_STAT_SAMPLE_MEAN); @@ -143,4 +146,6 @@ OPTION_PARSE(options->shutterIter, recipe, "SHUTTER.ITER", S32); OPTION_PARSE(options->shutterRej, recipe, "SHUTTER.REJECT", F32); + OPTION_PARSE(options->maskSuspect, recipe, "MASK.SUSPECT", F32); + OPTION_PARSE(options->maskBad, recipe, "MASK.BAD", F32); options->combine->combine = parseStat(recipe, "COMBINE"); @@ -159,4 +164,5 @@ options->fringe = false; options->shutter = false; + options->mask = false; } else if (strcasecmp(type, "DARK") == 0) { options->zero = false; @@ -165,4 +171,5 @@ options->fringe = false; options->shutter = false; + options->mask = false; } else if (strcasecmp(type, "FLAT") == 0) { options->zero = false; @@ -171,4 +178,5 @@ options->fringe = false; options->shutter = false; + options->mask = false; } else if (strcasecmp(type, "FRINGE") == 0) { options->zero = true; @@ -177,4 +185,5 @@ options->fringe = true; options->shutter = false; + options->mask = false; } else if (strcasecmp(type, "SHUTTER") == 0) { options->zero = false; @@ -183,4 +192,13 @@ options->fringe = false; options->shutter = true; + options->mask = false; + options->rows = 0; // Read the whole image at once + } else if (strcasecmp(type, "MASK") == 0) { + options->zero = false; + options->scale = false; + options->darktime = false; + options->fringe = false; + options->shutter = false; + options->mask = true; } else { psLogMsg(__func__, PS_LOG_WARN, "Unrecognised image type: %s --- assuming BIAS.\n", type); @@ -190,4 +208,5 @@ options->fringe = false; options->shutter = false; + options->mask = false; } } else { @@ -198,4 +217,5 @@ options->fringe = false; options->shutter = false; + options->mask = false; } Index: trunk/ppMerge/src/ppMergeOptions.h =================================================================== --- trunk/ppMerge/src/ppMergeOptions.h (revision 9955) +++ trunk/ppMerge/src/ppMergeOptions.h (revision 9996) @@ -24,4 +24,5 @@ bool fringe; // Make fringe measurements? bool shutter; // Generate shutter correction? + bool mask; // Generate bad pixel mask? unsigned int sample; // Sampling factor for measuring the background psStatsOptions mean; // Statistic to use to measure the mean @@ -34,4 +35,6 @@ int shutterIter; // Number of iterations for shutter measurement float shutterRej; // Rejection limit for shutter measurement + float maskSuspect; // Threshold for identifying suspect pixels + float maskBad; // Threshold for identifying bad pixels ppOnOff onOff; // On/off pairs? pmCombineParams *combine; // Combination parameters