Index: trunk/ppMerge/src/ppMergeMask.c =================================================================== --- trunk/ppMerge/src/ppMergeMask.c (revision 15937) +++ trunk/ppMerge/src/ppMergeMask.c (revision 17227) @@ -5,6 +5,6 @@ #include #include -#include #include + #include #include @@ -12,203 +12,394 @@ #include "ppMerge.h" -#include "ppMergeData.h" -#include "ppMergeMask.h" - - -// Generate a mask -bool ppMergeMask(ppMergeData *data, // Data - ppMergeOptions *options, // Options - pmConfig *config // Configuration - ) + +static bool mergeMask(pmConfig *config, // Configuration + const pmFPAview *view, // View to chip + bool writeOut, // Write output? + psRandom *rng, // Random number generator + psMetadata *stats // Statistics output + ) { - for (int i = 0; i < 2; i++) { - if (!ppMergeMaskSuspect (data, options, config)) { - psError(PS_ERR_UNKNOWN, true, "failed on mask suspect.\n"); - return false; - } - - if (!ppMergeMaskBad (data, options, config)) { - psError(PS_ERR_UNKNOWN, true, "failed on mask bad.\n"); - return false; - } - } - - if (!ppMergeMaskAverageConcepts (data, options, config)) { - psError(PS_ERR_UNKNOWN, true, "failed on average concepts.\n"); - return false; - } - - if (!ppMergeMaskGrow (data, options, config)) { - psError(PS_ERR_UNKNOWN, true, "failed on mask grow.\n"); - return false; - } - - if (!ppMergeMaskWrite (data, options, config)) { - psError(PS_ERR_UNKNOWN, true, "failed on mask write.\n"); - return false; - } + assert(config); + assert(view); + assert(view->chip != -1 && view->cell == -1 && view->readout == -1); + + bool mdok; // Status of MD lookup + int numFiles = psMetadataLookupS32(NULL, config->arguments, "INPUTS.NUM"); // Number of input files + psStatsOptions meanStat = psMetadataLookupS32(NULL, config->arguments, "MEAN"); // Statistic for mean + psStatsOptions stdevStat = psMetadataLookupS32(NULL, config->arguments, "STDEV"); // Statistic for stdev + psMaskType maskVal = psMetadataLookupU8(NULL, config->arguments, "MASKVAL"); // Value to mask + int sample = psMetadataLookupS32(NULL, config->arguments, "SAMPLE"); // Size of sample for statistics + bool chipStats = psMetadataLookupBool(&mdok, config->arguments, "MASK.CHIPSTATS"); // Statistics on chip? + float maskSuspect = psMetadataLookupF32(NULL, config->arguments, "MASK.SUSPECT"); // Threshold for suspect pixels + float maskBad = psMetadataLookupF32(NULL, config->arguments, "MASK.BAD"); // Threshold for bad pixels + pmMaskIdentifyMode maskMode = psMetadataLookupS32(NULL, config->arguments, "MASK.MODE"); // Mode for identifying bad pixels + int maskGrow = psMetadataLookupS32(NULL, config->arguments, "MASK.GROW"); // Radius to grow mask + psMaskType maskGrowVal = psMetadataLookupU8(NULL, config->arguments, "MASK.GROWVAL"); // Value for grown mask + + psStats *statistics = psStatsAlloc(meanStat | stdevStat); // Statistics for background + + psString outName = ppMergeOutputFile(config); // Name of output file + pmChip *outChip = pmFPAfileThisChip(config->files, view, outName); // Output chip + psFree(outName); + + // For each input file, get the statistics, which can be calculated at the chip or cell levels + psVector *values = psVectorAlloc(sample, PS_TYPE_F32); // Pixel values for statistics + pmFPAview *inView = pmFPAviewAlloc(0); // View for input + for (int i = 0; i < numFiles; i++) { + pmFPAfileActivate(config->files, false, NULL); + psArray *files = ppMergeFileActivateSingle(config, PPMERGE_FILES_INPUT, true, i); // Input files + pmFPAfile *input = files->data[0]; // Input file + psFree(files); + pmFPA *inFPA = input->fpa; // Input FPA + *inView = *view; + + int valueIndex = 0; // Index for vector of pixel values + + pmCell *inCell; // Input cell + while ((inCell = pmFPAviewNextCell(inView, inFPA, 1))) { + pmHDU *hdu = pmHDUFromCell(inCell); // HDU for cell + if (!hdu || hdu->blankPHU) { + // No data here + continue; + } + psTrace("ppMerge", 1, "Getting suspect pixels for file %d chip %d cell %d", + i, inView->chip, inView->cell); + + if (!pmFPAfileIOChecks(config, inView, PM_FPA_BEFORE)) { + psFree(inView); + goto MERGE_MASK_ERROR; + } + + if (!ppMergeFileOpenInput(config, view, i)) { + psError(PS_ERR_UNKNOWN, false, "Unable to open file %d", i); + psFree(inView); + goto MERGE_MASK_ERROR; + } + + if (inCell->readouts->n > 1) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, + "File %d chip %d cell %d contains more than one readout (%ld)", + i, inView->chip, inView->cell, inCell->readouts->n); + psFree(inView); + goto MERGE_MASK_ERROR; + } + + pmReadout *readout; + if (inCell->readouts && inCell->readouts->n == 1) { + readout = psMemIncrRefCounter(inCell->readouts->data[0]); // Input readout + } else { + readout = pmReadoutAlloc(inCell); + } + + if (!ppMergeFileReadInput(config, readout, i, 0)) { + psError(PS_ERR_UNKNOWN, false, "Unable to read readout %d", i); + psFree(inView); + psFree(readout); + goto MERGE_MASK_ERROR; + } + + pmCell *outCell = pmFPAfileThisCell(config->files, inView, "PPMERGE.OUTPUT.MASK"); // Output cell + pmReadout *outRO = NULL; // Output readout + if (outCell->readouts && outCell->readouts->n == 1) { + outRO = psMemIncrRefCounter(outCell->readouts->data[0]); + } else { + outRO = pmReadoutAlloc(outCell); + } + psImage *outMask = outRO->mask; // Output mask image (for iterative generation of mask) + + psImage *image = readout->image, *mask = readout->mask; // Image and mask + int numCols = readout->image->numCols, numRows = readout->image->numRows; // Image size + int numPix = numCols * numRows; // Number of pixels + int numCells = chipStats ? readout->parent->parent->cells->n : 1; // Number of cells + int num = PS_MIN(numPix, sample / numCells); // Number of values to add + if (!chipStats) { + valueIndex = 0; + } + for (int i = 0; i < num; i++) { + int pixel = numPix * psRandomUniform(rng); + int x = pixel % numCols; + int y = pixel / numCols; + if ((mask && (mask->data.PS_TYPE_MASK_DATA[y][x] & maskVal)) || + !isfinite(image->data.F32[y][x]) || + (outMask && (outMask->data.PS_TYPE_MASK_DATA[y][x] & maskVal))) { + continue; + } + + values->data.F32[valueIndex++] = image->data.F32[y][x]; + } + + if (!chipStats) { + values->n = valueIndex; + if (!psVectorStats(statistics, values, NULL, NULL, 0)) { + psError(PS_ERR_UNKNOWN, false, "Unable to do statistics on readout."); + psFree(inView); + psFree(readout); + goto MERGE_MASK_ERROR; + } + + if (!pmMaskFlagSuspectPixels(outRO, readout, psStatsGetValue(statistics, meanStat), + psStatsGetValue(statistics, stdevStat), maskSuspect, maskVal)) { + psError(PS_ERR_UNKNOWN, false, "Unable to find suspect values in file %d", i); + psFree(inView); + psFree(readout); + goto MERGE_MASK_ERROR; + } + pmCellFreeData(inCell); + + if (!pmFPAfileIOChecks(config, inView, PM_FPA_AFTER)) { + psFree(inView); + psFree(readout); + goto MERGE_MASK_ERROR; + } + } + psFree(readout); + psFree(outRO); + } + + // Additional run through cells if we want chip-level statistics + if (chipStats && valueIndex > 0) { + values->n = valueIndex; + if (!psVectorStats(statistics, values, NULL, NULL, 0)) { + psError(PS_ERR_UNKNOWN, false, "Unable to do statistics on chip."); + goto MERGE_MASK_ERROR; + } + inView->cell = -1; + while ((inCell = pmFPAviewNextCell(inView, inFPA, 1))) { + pmHDU *hdu = pmHDUFromCell(inCell); // HDU for cell + if (!hdu || hdu->blankPHU) { + // No data here + continue; + } + pmReadout *readout = inCell->readouts->data[0]; // Readout of interest + + inView->readout = 0; + pmReadout *outRO = pmFPAfileThisReadout(config->files, inView, "PPMERGE.OUTPUT.MASK"); + + if (!pmMaskFlagSuspectPixels(outRO, readout, psStatsGetValue(statistics, meanStat), + psStatsGetValue(statistics, stdevStat), maskSuspect, maskVal)) { + psError(PS_ERR_UNKNOWN, false, "Unable to find suspect values in file %d", i); + goto MERGE_MASK_ERROR; + } + pmCellFreeData(inCell); + + inView->readout = -1; + if (!pmFPAfileIOChecks(config, inView, PM_FPA_AFTER)) { + psFree(inView); + goto MERGE_MASK_ERROR; + } + } + } + } + psFree(inView); + psFree(statistics); statistics = NULL; + psFree(values); values = NULL; + + + // Another run through the chip to threshold on the suspects + pmFPAfileActivate(config->files, false, NULL); + if (writeOut) { + ppMergeFileActivate(config, PPMERGE_FILES_OUTPUT, true); + } + pmFPA *outFPA = outChip->parent; // Output FPA + pmCell *outCell; // Output cell + pmFPAview *outView = pmFPAviewAlloc(0); // View into output FPA + *outView = *view; + while ((outCell = pmFPAviewNextCell(outView, outFPA, 1))) { + pmHDU *hdu = pmHDUFromCell(outCell); // HDU for cell + if (!hdu || hdu->blankPHU) { + // No data here + continue; + } + + psTrace("ppMerge", 1, "Getting bad pixels for chip %d cell %d", outView->chip, outView->cell); + + assert(outCell->readouts && outCell->readouts->n == 1); + pmReadout *outRO = outCell->readouts->data[0]; // Output readout + if (!pmMaskIdentifyBadPixels(outRO, maskVal, maskBad, maskMode)) { + psError(PS_ERR_UNKNOWN, false, "Unable to mask bad pixels"); + goto MERGE_MASK_ERROR; + } + + // Supplementary outputs + { + // The counts image is fairly useless, but it preserves the model + pmCell *countsCell = pmFPAfileThisCell(config->files, outView, "PPMERGE.OUTPUT.COUNT"); + pmReadout *countsRO = pmReadoutAlloc(countsCell); // Readout with count of inputs per pixel + countsRO->image = psImageAlloc(outRO->mask->numCols, outRO->mask->numRows, PS_TYPE_F32); + psImageInit(countsRO->image, numFiles); + countsRO->data_exists = countsCell->data_exists = countsCell->parent->data_exists = true; + psFree(countsRO); + + pmCell *sigmaCell = pmFPAfileThisCell(config->files, outView, "PPMERGE.OUTPUT.SIGMA"); + pmReadout *sigmaRO = pmReadoutAlloc(sigmaCell); // Readout with suspect image + psImage *suspect = psMetadataLookupPtr(NULL, outRO->analysis, PM_MASK_ANALYSIS_SUSPECT); + sigmaRO->image = psImageCopy(sigmaRO->image, suspect, PS_TYPE_F32); + psMetadataRemoveKey(outRO->analysis, PM_MASK_ANALYSIS_SUSPECT); + sigmaRO->data_exists = sigmaCell->data_exists = sigmaCell->parent->data_exists = true; + psFree(sigmaRO); + } + + if (maskGrowVal > 0) { + psImage *grown = psImageGrowMask(NULL, outRO->mask, maskVal, maskGrow, maskGrowVal); // Grown mask + psFree(outRO->mask); + outRO->mask = grown; + } + + if (writeOut) { + if (!pmFPAfileIOChecks(config, outView, PM_FPA_BEFORE)) { + psFree(outView); + goto MERGE_MASK_ERROR; + } + + outRO->data_exists = outCell->data_exists = outChip->data_exists = true; + + // Average concepts + psList *cells = psListAlloc(NULL); // List of cells, for concept averaging + for (int i = 0; i < numFiles; i++) { + pmFPAfile *inFile = pmFPAfileSelectSingle(config->files, "PPMERGE.INPUT", i); // Input file + pmCell *inCell = pmFPAviewThisCell(outView, inFile->fpa); // Input cell + psListAdd(cells, PS_LIST_TAIL, inCell); + } + if (!pmConceptsAverageCells(outCell, cells, NULL, NULL, true)) { + psError(PS_ERR_UNKNOWN, false, "Unable to average cell concepts."); + psFree(cells); + psFree(outRO); + psFree(outView); + return false; + } + psFree(cells); + + // Statistics on the merged cell using a fake image + outRO->image = psImageAlloc(outRO->mask->numCols, outRO->mask->numRows, PS_TYPE_F32); + psImageInit(outRO->image, 1.0); + if (!ppStatsFPA(stats, outRO->parent->parent->parent, outView, + maskVal | pmConfigMask("BLANK", config), config)) { + psError(PS_ERR_UNEXPECTED_NULL, true, "Unable to generate stats for image."); + psFree(outRO); + psFree(outView); + return false; + } + psFree(outRO->image); + outRO->image = NULL; + + // Write + if (!pmFPAfileIOChecks(config, outView, PM_FPA_AFTER)) { + psFree(outView); + goto MERGE_MASK_ERROR; + } + } + } + psFree(outView); + + ppMergeFileActivate(config, PPMERGE_FILES_ALL, true); return true; + + +MERGE_MASK_ERROR: + psFree(statistics); + psFree(values); + return false; } -bool ppMergeMaskReadoutStats(const pmReadout *readout, - const pmReadout *output, - ppMergeOptions *options, // Options - psRandom *rng) + + + + +bool ppMergeMask(pmConfig *config) { - 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 *mask = NULL; - psImage *image = readout->image; // Image of interest - - if (output) { - mask = output->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); - + assert(config); + + bool mdok; // Status of MD lookup + int numFiles = psMetadataLookupS32(NULL, config->arguments, "INPUTS.NUM"); // Number of inputs + bool haveMasks = psMetadataLookupBool(&mdok, config->arguments, "INPUTS.MASKS"); // Do we have masks? + bool haveWeights = psMetadataLookupBool(&mdok, config->arguments, "INPUTS.WEIGHTS"); // Do we have weights? + int iter = psMetadataLookupS32(NULL, config->arguments, "ITER"); // Number of rejection iterations + + PS_ASSERT_INT_POSITIVE(iter, false); + + pmFPAview *view = pmFPAviewAlloc(0); // View to component of interest + psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS, 0); // Random number generator + + psMetadata *stats = NULL; // Statistics for output + if (psMetadataLookup(config->arguments, "STATS.NAME")) { + stats = psMetadataAlloc(); + psMetadataAddMetadata(config->arguments, PS_LIST_TAIL, "STATS.DATA", 0, "Statistics output", stats); + } + + psArray *inputs = ppMergeFileDataLevel(config, "PPMERGE.INPUT", PM_FPA_LEVEL_READOUT); // Input images + psFree(inputs); + if (haveMasks) { + psArray *masks = ppMergeFileDataLevel(config, "PPMERGE.INPUT.MASK", PM_FPA_LEVEL_READOUT); + psFree(masks); + } + if (haveWeights) { + psArray *weights = ppMergeFileDataLevel(config, "PPMERGE.INPUT.WEIGHT", PM_FPA_LEVEL_READOUT); + psFree(weights); + } + + if (!ppMergeFileActivate(config, PPMERGE_FILES_INPUT, true)) { + psError(PS_ERR_UNKNOWN, false, "Unable to activate files."); + goto PPMERGE_MASK_ERROR; + } + if (!ppMergeFileActivate(config, PPMERGE_FILES_OUTPUT, true)) { + psError(PS_ERR_UNKNOWN, false, "Unable to activate files."); + goto PPMERGE_MASK_ERROR; + } + + psString outName = ppMergeOutputFile(config); // Name of output file + pmFPAfile *output = psMetadataLookupPtr(NULL, config->files, outName); // Output file + psFree(outName); + assert(output && output->fpa); + pmFPA *outFPA = output->fpa; // Output FPA + + if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) { + goto PPMERGE_MASK_ERROR; + } + pmChip *outChip; // Chip of interest + while ((outChip = pmFPAviewNextChip(view, outFPA, 1))) { + if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) { + goto PPMERGE_MASK_ERROR; + } + + for (int i = 0; i < iter; i++) { + if (!mergeMask(config, view, (i == iter - 1), rng, stats)) { + psError(PS_ERR_UNKNOWN, false, "Unable to merge chip %d", view->chip); + goto PPMERGE_MASK_ERROR; + } + } + + if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) { + goto PPMERGE_MASK_ERROR; + } + } + + psList *fpaList = psListAlloc(NULL);// List of FPAs for concept averaging + for (int i = 0; i < numFiles; i++) { + pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PPMERGE.INPUT", i); // Input file + psListAdd(fpaList, PS_LIST_TAIL, file->fpa); + } + if (!pmConceptsAverageFPAs(outFPA, fpaList)) { + psError(PS_ERR_UNKNOWN, false, "Unable to average FPA concepts."); + psFree(fpaList); + goto PPMERGE_MASK_ERROR; + } + psFree(fpaList); + + if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) { + goto PPMERGE_MASK_ERROR; + } + + psFree(stats); + psFree(view); psFree(rng); + return true; + +PPMERGE_MASK_ERROR: + psFree(stats); + psFree(view); + psFree(rng); + return false; } -bool ppMergeMaskChipStats (const pmChip *chip, - const pmChip *output, - 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 - - pmCell *cellOutput = output->cells->data[nCell]; - pmReadout *readoutOutput = NULL; - psImage *mask = NULL; - if (cellOutput->readouts->n > 0) { - readoutOutput = cellOutput->readouts->data[0]; - mask = readoutOutput->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 ++; - } - } - } - - // no valid data, skip the chip - if (!values->n) { - psFree(values); - psFree(stats); - psFree(rng); - return true; - } - - // 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, "invalid 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); - } - } - - psLogMsg ("ppMerge", PS_LOG_INFO, "statistics for chip: %f +/- %f\n", stats->robustMedian, stats->robustStdev); - - psFree(values); - psFree(stats); - psFree(rng); - return true; -}