Index: unk/ppMerge/src/ppMergeCheckInputs.c =================================================================== --- /trunk/ppMerge/src/ppMergeCheckInputs.c (revision 17232) +++ (revision ) @@ -1,176 +1,0 @@ -#ifdef HAVE_CONFIG_H -#include -#endif - -#include -#include -#include -#include -#include - -#include "ppMerge.h" -#include "ppMergeCheckInputs.h" -#include "ppMergeData.h" - -// Check input files to make sure everything's consistent -ppMergeData *ppMergeCheckInputs(ppMergeOptions *options, // Options - pmConfig *config // Configuration - ) -{ - ppMergeData *data = ppMergeDataAlloc(); // The data, to return - - // Output file - psString outName = psMetadataLookupStr(NULL, config->arguments, "OUTPUT"); // The output file name - assert(outName); // It should be there! - outName = pmConfigConvertFilename(outName, config, true); - data->outFile = psFitsOpen(outName, "w"); // Output FITS file - if (!data->outFile) { - // There's no point in continuing if we can't open the output - psErrorStackPrint(stderr, "Can't open output image: %s\n", outName); - psFree(outName); - exit(EXIT_FAILURE); - } - psFree(outName); - - // Statistics file - psString statsName = psMetadataLookupStr(NULL, config->arguments, "-stats"); // Name for statistics file - if (statsName && strlen(statsName) > 0) { - psString resolved = pmConfigConvertFilename(statsName, config, true); // Resolved filename - data->statsFile = fopen(resolved, "w"); - if (!data->statsFile) { - psError(PS_ERR_IO, true, "Unable to open statistics file %s for writing.\n", resolved); - psFree(resolved); - psFree(data); - return NULL; - } - psFree(resolved); - } - - psArray *filenames = psMetadataLookupPtr(NULL, config->arguments, "INPUT"); // The input file names - assert(filenames); - if (!data->in) { - data->in = psArrayAlloc(filenames->n); - } - if (!data->files) { - data->files = psArrayAlloc(filenames->n); - } - int numGood = 0; // Number of good files - for (int i = 0; i < filenames->n; i++) { - psString name = filenames->data[i]; // The name of the file - if (!name || strlen(name) == 0) { - continue; - } - psTrace("ppMerge", 1, "Checking input file %s....\n", name); - psString resolved = pmConfigConvertFilename(name, config, false); // Resolved file name - psFits *inFile = psFitsOpen(resolved, "r"); // The FITS file to read - if (!inFile) { - psLogMsg(__func__, PS_LOG_WARN, "Unable to open input file %s --- ignored.\n", resolved); - // Kick it out - psFree(filenames->data[i]); - filenames->data[i] = NULL; - continue; - } - psFree(resolved); - psMetadata *header = psFitsReadHeader(NULL, inFile); // The FITS (primary) header - data->files->data[i] = inFile; - - // The formats must be identical. The chief reason for this is so that we know what output format to - // use. I guess one could specify a different output format on the command line, but how do we - // generate a PHU for that? Perhaps we could revisit this restriction in the future (construct an - // FPAview from the specified camera format configuration, and use pmFPAAddSourceFromView), but for - // now it's less hassle just to limit the output format to be the input format. - if (!options->format) { - options->format = pmConfigCameraFormatFromHeader(config, header, true); - if (!options->format) { - psLogMsg(__func__, PS_LOG_WARN, "Unable to identify camera format for input file %s --- " - "ignored.\n", name); - // Kick it out - psFree(header); - data->in->data[i] = NULL; - continue; - } - } else { - bool valid = false; - if (!pmConfigValidateCameraFormat(&valid, options->format, header)) { - psError (PS_ERR_UNKNOWN, false, "Error in config scripts\n"); - exit (PS_EXIT_CONFIG_ERROR); - } - if (!valid) { - psLogMsg(__func__, PS_LOG_WARN, "Input file %s doesn't match camera format --- ignored.\n", name); - // Kick it out - psFree(header); - data->in->data[i] = NULL; - continue; - } - } - - pmFPA *fpa = pmFPAConstruct(config->camera); - pmFPAview *view = pmFPAAddSourceFromHeader(fpa, header, options->format); - psFree(view); - - // Cull chips and cells that don't have data - // Otherwise the abundance of metadata in the concepts (esp. for GPC) can overload the memory - psArray *chips = fpa->chips; // Array of chips in output - for (int i = 0; i < chips->n; i++) { - pmChip *chip = chips->data[i]; // Chip of interest - psArray *cells = chip->cells; // Array of cells - int culled = 0; // Number of culled cells - for (int j = 0; j < cells->n; j++) { - pmCell *cell = cells->data[j]; - pmHDU *hdu = pmHDUGetLowest(fpa, chip, cell); // HDU for cell - if (!hdu || hdu->blankPHU) { - psFree(cell->concepts); - cell->concepts = NULL; - culled++; - } - } - if (culled == cells->n) { - psFree(chip->concepts); - chip->concepts = NULL; - } - } - data->in->data[i] = fpa; - - - // Use the first valid input as the basis for the output --- including the header - if (!data->out) { - psTrace("ppMerge", 5, "Constructing output using %s as a template.\n", name); - data->out = pmFPAConstruct(config->camera); - pmFPAview *view = pmFPAAddSourceFromHeader(data->out, header, options->format); - psFree(view); - } - psFree(header); - - psTrace("ppMerge", 3, "%s checks out.\n", name); - numGood++; - } - - // Count the cells - int numCells = 0; // Number of cells in the output FPA - psArray *chips = data->out->chips; // Array of chips in output - for (int i = 0; i < chips->n; i++) { - pmChip *chip = chips->data[i]; // Chip of interest - if (!chip) { - continue; - } - psArray *cells = chip->cells; // Array of cells - for (int j = 0; j < cells->n; j++) { - pmCell *cell = cells->data[j]; - if (cell) { - numCells++; - } - } - } - data->numCells = numCells; - psTrace("ppMerge", 3, "Output has %d cells.\n", numCells); - - psTrace("ppMerge", 3, "We have %d good inputs.\n", numGood); - if (numGood > 1) { - return data; - } - - psFree(data); - return NULL; -} - - Index: unk/ppMerge/src/ppMergeCheckInputs.h =================================================================== --- /trunk/ppMerge/src/ppMergeCheckInputs.h (revision 17232) +++ (revision ) @@ -1,13 +1,0 @@ -#ifndef PP_MERGE_CHECK_INPUTS_H -#define PP_MERGE_CHECK_INPUTS_H - -#include -#include "ppMergeOptions.h" -#include "ppMergeData.h" - -// Check input files to make sure everything's consistent -ppMergeData *ppMergeCheckInputs(ppMergeOptions *options, // Options - pmConfig *config // Configuration - ); - -#endif Index: unk/ppMerge/src/ppMergeCombine.c =================================================================== --- /trunk/ppMerge/src/ppMergeCombine.c (revision 17232) +++ (revision ) @@ -1,434 +1,0 @@ -#ifdef HAVE_CONFIG_H -#include -#endif - -#include -#include -#include -#include -#include -#include -#include - -#include "ppMerge.h" -#include "ppMergeData.h" -#include "ppMergeCombine.h" -#include "ppMergeVersion.h" - -#define TESTING - - -#if 0 -static FILE *dumpFile = NULL; - -static psMemId mbAlloc(psMemBlock *mb) -{ - if (!dumpFile) { - dumpFile = fopen("memBlocks.dat", "w"); - } - fprintf(dumpFile, "Alloc: %12lu\t%12zd\t%s:%d\n", mb->id, mb->userMemorySize, - mb->file, mb->lineno); - return 1; -} - -static void dumpDone(void) -{ - fclose(dumpFile); - exit(EXIT_FAILURE); -} -#endif - -#if 0 -static psMemId memId = 0; -static void memDump(void) -{ - psMemBlock **leaks = NULL; - int numLeaks = psMemCheckLeaks(memId, &leaks, NULL, true); - FILE *memFile = fopen("mem.dat", "w"); - fprintf(memFile, "# MemBlock Size Source\n"); - for (int i = 0; i < numLeaks; i++) { - psMemBlock *mb = leaks[i]; - fprintf(memFile, "%12lu\t%12zd\t%s:%d\n", mb->id, mb->userMemorySize, - mb->file, mb->lineno); - } - fclose(memFile); - psFree(leaks); -} -#endif - -#if 0 -static void memCheck(void) -{ - return; - if (psTraceGetLevel("ppMerge") > 9) { - psMemBlock **leaks = NULL; - int numLeaks = psMemCheckLeaks(0, &leaks, NULL, true); - size_t largestSize = 0; - psMemId largest = 0; - size_t totalSize = 0; - for (int i = 0; i < numLeaks; i++) { - psMemBlock *mb = leaks[i]; - totalSize += mb->userMemorySize; - if (mb->userMemorySize > largestSize) { - largestSize = mb->userMemorySize; - largest = mb->id; - } - } - psFree(leaks); - psTrace("ppMerge", 0, "Memory in use: %zd\n", totalSize); - psTrace("ppMerge", 0, "Largest block: %ld\n", largest); - psTrace("ppMerge", 0, "sbrk(): %p\n", sbrk(0)); - } - return; -} -#endif - - -// Combine the inputs -bool ppMergeCombine(psImage *scales, // Scales for each cell of each integration, or NULL - psImage *zeros, // Zeros for each cell of each integration, or NULL - psArray *shutters, // Shutter correction data for each cell, or NULL - 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 - - // Sanity checks - assert(!options->scale || scales); - assert(!scales || (scales->type.type == PS_TYPE_F32 && - scales->numCols == data->numCells && - scales->numRows == filenames->n)); - assert(!options->zero || zeros); - assert(!zeros || (zeros->type.type == PS_TYPE_F32 && - zeros->numCols == data->numCells && - zeros->numRows == filenames->n)); - assert(!options->shutter || shutters); - assert(!shutters || (shutters->n == data->numCells)); - - // Iterate over the FPA - pmFPA *fpa = data->out; // Output FPA - pmFPAview *view = pmFPAviewAlloc(0);// View of FPA, for iteration - int cellNum = -1; // Cell number in the whole FPA - if (data->out->hdu) { - pmFPAUpdateNames(data->out, NULL, NULL); - } - 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 - pmHDU *lastHDU = NULL; // Last HDU to be updated - if (options->mask) { - rng = psRandomAlloc(PS_RANDOM_TAUS, 0); - } - while ((chip = pmFPAviewNextChip(view, fpa, 1))) { - if (chip->hdu) { - // Data will exist soon - pmFPAUpdateNames(data->out, chip, NULL); - chip->data_exists = true; - } - pmChipWrite(chip, data->outFile, config->database, true, false); // Write header only - pmCell *cell; // Cell of interest - while ((cell = pmFPAviewNextCell(view, fpa, 1))) { - cellNum++; - - pmHDU *hdu = pmHDUGetLowest(data->out, chip, cell); // HDU for cell - if (!hdu || hdu->blankPHU) { - pmCellWrite(cell, data->outFile, config->database, true); // Write header only - continue; - } - if (cell->hdu) { - // Data will exist soon - pmFPAUpdateNames(data->out, chip, cell); - chip->data_exists = cell->data_exists = true; - } - pmReadout *readout = pmReadoutAlloc(cell); // Output readout of interest - psArray *stack = psArrayAlloc(filenames->n); // Stack of readouts to combine - psVector *cellScales = NULL; // Scales for this cell - if (scales) { - cellScales = psImageCol(NULL, scales, cellNum); - } - psVector *cellZeros = NULL; // Zeros for this cell - if (zeros) { - cellZeros = psImageCol(NULL, zeros, cellNum); - } - - // Read bit by bit - int numRead; // Number of inputs read - int numScan = 0; - - // Put version metadata into header - if (hdu && hdu != lastHDU) { - if (!hdu->header) { - hdu->header = psMetadataAlloc(); - } - ppMergeVersionMetadata(hdu->header); - lastHDU = hdu; - } - - float shutterRef = NAN; // Reference shutter correction - if (options->shutter) { - shutterRef = pmShutterCorrectionReference(shutters->data[cellNum]); - } - - 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) - bool readOK; - if (pmReadoutReadNext(&readOK, stack->data[i], fits, 0, options->rows)) { - if (!readOK) { - psError(PS_ERR_IO, false, "Failed to read concepts for cell.\n"); - psErrorStackPrint(stderr, "trouble reading data!\n"); - exit (1); - } - // If we're creating a bias or a dark, we don't want to generate a mask - if ((options->zero || options->scale || options->shutter || options->mask) && - options->combine->maskVal) { - pmReadoutSetMask(stack->data[i], options->satMask, options->badMask); - } - - // If we're combining with weights, we want to generate weights. - if (options->combine->weights && !options->dark) { - - // If it's a bias or dark, set the gain to zero: noise only contributed by read - if (!options->zero && !options->scale) { - pmReadoutSetWeight(stack->data[i], false); - } else { - pmReadoutSetWeight(stack->data[i], true); - } - } - - numRead++; - } else { - psTrace("ppMerge", 3, "Finished reading file %d, 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) { - if (options->shutter) { - pmShutterCorrectionGenerate(readout, NULL, stack, shutterRef, shutters->data[cellNum], - options->shutterIter, options->shutterRej, - options->combine->maskVal); - } else if (options->dark) { - pmDarkCombine(cell, stack, options->darkOrdinates, options->darkNorm, - options->combine->iter, options->combine->rej, - options->combine->maskVal); - } else { - pmReadoutCombine(readout, stack, cellZeros, cellScales, options->combine); - } - } - numScan++; - - } while (numRead > 0); - - // Get list of cells for concepts averaging - psList *inCells = psListAlloc(NULL); // List of cells - for (int i = 0; i < filenames->n; i++) { - if (! filenames->data[i] || strlen(filenames->data[i]) == 0) { - continue; - } - pmCell *cellIn = pmFPAviewThisCell(view, data->in->data[i]); // Input cell - psListAdd(inCells, PS_LIST_TAIL, cellIn); - } - if (!pmConceptsAverageCells(cell, inCells, NULL, NULL, true)) { - psError(PS_ERR_UNKNOWN, false, "Unable to average cell concepts."); - psFree(inCells); - return false; - } - psFree(inCells); - - psFree(stack); - -#if 0 - // Set the dark time for the output image, since we normalised - if (options->darktime) { - psMetadataItem *darkItem = psMetadataLookup(cell->concepts, "CELL.DARKTIME"); - darkItem->data.F32 = 1.0; - psMetadataItem *expItem = psMetadataLookup(cell->concepts, "CELL.EXPOSURE"); - expItem->data.F32 = 1.0; - } -#endif - - // Measure the fringes for this cell - // - // XXX Need to deal with multiple components: we will do this by building up the components - // Read the existing fringe measurements - // Use existing regions to measure fringe statistics - // Add the new fringe measurements to the existing fringe measurements. - // Write the appended fringe measurements. - // Read in the "output" file to get the existing components. - // Put the new readout into the cell after the existing readouts. - if (options->fringe && readout->image) { - pmFringeRegions *regions = pmFringeRegionsAlloc(options->fringeNum, options->fringeSize, - options->fringeSize, options->fringeSmoothX, - options->fringeSmoothY); // Fringe regions - pmFringeStats *fringe = pmFringeStatsMeasure(regions, readout, options->combine->maskVal); - psFree(regions); - if (!fringe) { - psError(PS_ERR_UNKNOWN, false, "Unable to measure fringe statistics.\n"); - psFree(readout); - return false; - } - - psArray *fringes = psArrayAlloc(1); // Array of fringes - fringes->data[0] = fringe; - - pmFringesFormat(cell, NULL, fringes); - psFree(fringes); // Drop reference - } - - if (readout->image) { - // Add MD5 information for cell - pmHDU *hdu = pmHDUFromCell(cell); // HDU that owns the cell - const char *chipName = psMetadataLookupStr(NULL, chip->concepts, "CHIP.NAME"); - const char *cellName = psMetadataLookupStr(NULL, cell->concepts, "CELL.NAME"); - - psString headerName = NULL; // Header name for MD5 - psStringAppend(&headerName, "MD5_%s_%s", chipName, cellName); - - psVector *md5 = psImageMD5(readout->image); // md5 hash - psString md5string = psMD5toString(md5); // String - psFree(md5); - psMetadataAddStr(hdu->header, PS_LIST_TAIL, headerName, PS_META_REPLACE, - "Image MD5", md5string); - psFree(md5string); - psFree(headerName); - - - // Statistics on the merged cell - if (data->statsFile) { - if (!data->stats) { - data->stats = psMetadataAlloc(); - } - if (!ppStatsFPA(data->stats, data->out, view, - options->combine->maskVal | pmConfigMask("BLANK", config), - config)) { - psError(PS_ERR_UNEXPECTED_NULL, true, "Unable to generate stats for image.\n"); - return false; - } - } - } - - psFree(readout); // Drop reference - - - // We threw away the bias sections --- record this - psMetadataItem *biassecItem = psMetadataLookup(cell->concepts, "CELL.BIASSEC"); // Item of BIASSEC - psList *biassecList = biassecItem->data.V; // List of BIASSECs - while (psListRemove(biassecList, PS_LIST_TAIL)); // Removing all entries - - // Blow away the cell data - for (int i = 0; i < filenames->n; i++) { - pmFPA *fpaIn = data->in->data[i]; // Input FPA - if (!fpaIn) { continue; } // was not a valid input file - pmChip *chipIn = fpaIn->chips->data[view->chip]; // Input chip - pmCell *cellIn = chipIn->cells->data[view->cell]; // Input cell - pmCellFreeData(cellIn); - } - - // Write the pixels - if (cell->hdu && !cell->hdu->blankPHU) { - psTrace("ppMerge", 5, "Writing out cell HDU.\n"); - if (options->dark) { - pmCellWriteDark(cell, data->outFile, config->database, false); - } else { - pmCellWrite(cell, data->outFile, config->database, false); - if (options->fringe) { - pmCellWriteTable(data->outFile, cell, "FRINGE"); - } - } - pmCellFreeData(cell); - } - } - - // Blow away the chip data - for (int i = 0; i < filenames->n; i++) { - pmFPA *fpaIn = data->in->data[i]; // Input FPA - if (!fpaIn) { continue; } // was not a valid input file - pmChip *chipIn = fpaIn->chips->data[view->chip]; // Input chip - pmChipFreeData(chipIn); - } - - // Write the pixels - if (chip->hdu && !chip->hdu->blankPHU) { - psTrace("ppMerge", 5, "Writing out chip HDU.\n"); - if (options->dark) { - pmChipWriteDark(chip, data->outFile, config->database, false, false); - } else { - pmChipWrite(chip, data->outFile, config->database, false, false); - if (options->fringe) { - pmChipWriteTable(data->outFile, chip, "FRINGE"); - } - } - - pmChipFreeData(chip); - } - } - - // Blow away the FPA data - for (int i = 0; i < filenames->n; i++) { - pmFPA *fpaIn = data->in->data[i]; // Input FPA - if (!fpaIn) { continue; } // was not a valid input file - pmFPAFreeData(fpaIn); - } - - if (data->out->hdu && !data->out->hdu->blankPHU) { - // Write the pixels - psTrace("ppMerge", 5, "Writing out FPA HDU.\n"); - if (options->dark) { - pmFPAWriteDark(data->out, data->outFile, config->database, false, false); - } else { - pmFPAWrite(data->out, data->outFile, config->database, false, false); - if (options->fringe) { - pmFPAWriteTable(data->outFile, fpa, "FRINGE"); - } - } - } - - pmFPAFreeData(data->out); - - psFree(view); - psFree(rng); - - - // Get list of FPAs for concepts averaging - psList *inFPAs = psListAlloc(NULL); // List of FPAs - for (int i = 0; i < filenames->n; i++) { - if (! filenames->data[i] || strlen(filenames->data[i]) == 0) { - continue; - } - pmFPA *fpaIn = data->in->data[i]; // Input FPA - psListAdd(inFPAs, PS_LIST_TAIL, fpaIn); - } - - if (!pmConceptsAverageFPAs(fpa, inFPAs)) { - psError(PS_ERR_UNKNOWN, false, "Unable to average FPA concepts."); - psFree(inFPAs); - return false; - } - psFree(inFPAs); - - - return true; - -} Index: unk/ppMerge/src/ppMergeCombine.h =================================================================== --- /trunk/ppMerge/src/ppMergeCombine.h (revision 17232) +++ (revision ) @@ -1,19 +1,0 @@ -#ifndef PP_MERGE_COMBINE_H -#define PP_MERGE_COMBINE_H - -#include -#include - -#include "ppMergeData.h" -#include "ppMergeOptions.h" - -// Combine readouts -bool ppMergeCombine(psImage *scales, // Scales for each cell of each integration, or NULL - psImage *zeros, // Zeros for each cell of each integration, or NULL - psArray *shutters, // Shutter correction data for each cell, or NULL - ppMergeData *data, // Data - ppMergeOptions *options, // Options - pmConfig *config // Configuration - ); - -#endif Index: unk/ppMerge/src/ppMergeConfig.c =================================================================== --- /trunk/ppMerge/src/ppMergeConfig.c (revision 17232) +++ (revision ) @@ -1,93 +1,0 @@ -#ifdef HAVE_CONFIG_H -#include -#endif - -#include -#include -#include - -#include "ppMerge.h" -#include "ppMergeConfig.h" - -#define DONT_USE_DB - -// Output usage information -static void usage(const char *programName, // Name of the program - psMetadata *arguments // Arguments list - ) -{ - printf("Merge multiple calibration frames into a master frame by stacking.\n\n" - "Usage:\n" - "\t%s OUTPUT.fits INPUT1.fits INPUT2.fits ...\n" - "\n", programName); - psArgumentHelp(arguments); - psFree(arguments); - exit(EXIT_FAILURE); -} - -pmConfig *ppMergeConfig(int argc, char **argv) -{ - pmConfig *config = pmConfigRead(&argc, argv, PPMERGE_RECIPE); - // Load the site-wide configuration information - if (! config) { - psErrorStackPrint(stderr, "Can't find site configuration!\n"); - exit(EXIT_FAILURE); - } - - psMetadata *arguments = psMetadataAlloc(); // Command-line arguments - psMetadataAddStr(arguments, PS_LIST_TAIL, "-type", 0, "Type of calibration frame", ""); - psMetadataAddBool(arguments, PS_LIST_TAIL, "-zero", 0, "Subtract background?", false); - psMetadataAddBool(arguments, PS_LIST_TAIL, "-scale", 0, "Scale by background?", false); - psMetadataAddBool(arguments, PS_LIST_TAIL, "-exptime", 0, "Scale by the exposure time?", false); - psMetadataAddS32(arguments, PS_LIST_TAIL, "-onoff", 0, "Number of on/off pairs", 0); - psMetadataAddStr(arguments, PS_LIST_TAIL, "-stats", 0, "MDC file to hold statistics ", NULL); - - if (argc == 1) { - psFree(config); - usage(argv[0], arguments); - } - - // Parse the arguments - if (!psArgumentParse(arguments, &argc, argv) || argc < 3) { - psErrorStackPrint(stderr, "Unable to parse arguments."); - psFree(config); - usage(argv[0], arguments); - } - - // Add the output image to the arguments list - psMetadataAddStr(config->arguments, PS_LIST_TAIL, "OUTPUT", 0, "Name of the output image", - argv[1]); - - psMetadataCopy(config->arguments, arguments); - psFree(arguments); - - // Everything remaining must be input files - if (argc - 2 <= 1) { - psErrorStackPrint(stderr, "No files to combine.\n"); - exit(EXIT_FAILURE); - } - psArray *files = psArrayAlloc(argc - 2); - for (int i = 2; i < argc; i++) { - files->data[i - 2] = psStringCopy(argv[i]); - } - psMetadataAddPtr(config->arguments, PS_LIST_TAIL, "INPUT", PS_DATA_ARRAY, - "Array of inputs images", files); - psFree(files); // Drop reference - -#ifndef DONT_USE_DB -#ifdef HAVE_PSDB - // Define database handle, if required - config->database = pmConfigDB(config); -#endif -#endif - - // Add concepts for scale and zero - psMetadataItem *scaleItem = psMetadataItemAllocF32("PPMERGE.SCALE", "Scaling for ppMerge", NAN); - psMetadataItem *zeroItem = psMetadataItemAllocF32("PPMERGE.ZERO", "Zero offset for ppMerge", NAN); - pmConceptRegister(scaleItem, NULL, NULL, false, PM_FPA_LEVEL_CELL); - pmConceptRegister(zeroItem, NULL, NULL, false, PM_FPA_LEVEL_CELL); - psFree(scaleItem); - psFree(zeroItem); - - return config; -} Index: unk/ppMerge/src/ppMergeConfig.h =================================================================== --- /trunk/ppMerge/src/ppMergeConfig.h (revision 17232) +++ (revision ) @@ -1,10 +1,0 @@ -#ifndef PP_MERGE_CONFIG_H -#define PP_MERGE_CONFIG_H - -#include - -// Get the configuration information -pmConfig *ppMergeConfig(int argc, char **argv // The standard command-line parameters (but pointer to number) - ); - -#endif Index: unk/ppMerge/src/ppMergeData.c =================================================================== --- /trunk/ppMerge/src/ppMergeData.c (revision 17232) +++ (revision ) @@ -1,50 +1,0 @@ -#ifdef HAVE_CONFIG_H -#include -#endif - -#include -#include - -#include "ppMerge.h" -#include "ppMergeData.h" - -// Free function for ppMergeData -static void mergeDataFree(ppMergeData *data // Data to free - ) -{ - if (data->files) { - for (long i = 0; i < data->files->n; i++) { - psFitsClose(data->files->data[i]); - data->files->data[i] = NULL; - } - psFree(data->files); - } - psFree(data->in); - psFree(data->out); - if (data->outFile) { - psFitsClose(data->outFile); - data->outFile = NULL; - } - if (data->statsFile) { - fclose(data->statsFile); - data->statsFile = NULL; - } - psFree(data->stats); -} - -// Allocator for ppMergeData -ppMergeData *ppMergeDataAlloc(void) -{ - ppMergeData *data = psAlloc(sizeof(ppMergeData)); // The data, to return - psMemSetDeallocator(data, (psFreeFunc)mergeDataFree); - - data->numCells = 0; - data->files = NULL; - data->in = NULL; - data->out = NULL; - data->outFile = NULL; - data->stats = NULL; - data->statsFile = NULL; - - return data; -} Index: unk/ppMerge/src/ppMergeData.h =================================================================== --- /trunk/ppMerge/src/ppMergeData.h (revision 17232) +++ (revision ) @@ -1,23 +1,0 @@ -#ifndef PP_MERGE_DATA_H -#define PP_MERGE_DATA_H - -#include -#include - -// Container for the data -typedef struct { - int numCells; // Number of (valid) cells in the FPA - psArray *files; // Input file pointers - psArray *in; // Input FPA structures - pmFPA *out; // Output FPA structure - psFits *outFile; // FITS file handle for output - psMetadata *stats; // Statistics on the combined image - FILE *statsFile; // File stream for statistics output -} ppMergeData; - - -// Allocator -ppMergeData *ppMergeDataAlloc(void); - - -#endif Index: unk/ppMerge/src/ppMergeMask.h =================================================================== --- /trunk/ppMerge/src/ppMergeMask.h (revision 17232) +++ (revision ) @@ -1,50 +1,0 @@ -#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 - ); - -bool ppMergeMaskSuspect(ppMergeData *data, // Data - ppMergeOptions *options, // Options - pmConfig *config // Configuration - ); - -bool ppMergeMaskBad(ppMergeData *data, // Data - ppMergeOptions *options, // Options - pmConfig *config // Configuration - ); - -bool ppMergeMaskAverageConcepts(ppMergeData *data, // Data - ppMergeOptions *options, // Options - pmConfig *config // Configuration - ); - -bool ppMergeMaskWrite(ppMergeData *data, // Data - ppMergeOptions *options, // Options - pmConfig *config // Configuration - ); - -bool ppMergeMaskGrow(ppMergeData *data, // Data - ppMergeOptions *options, // Options - pmConfig *config // Configuration - ); - -bool ppMergeMaskGrowReadout (pmReadout *readout, psMaskType maskVal, int nPixels); - -bool ppMergeMaskReadoutStats(const pmReadout *readout, - const pmReadout *output, - ppMergeOptions *options, // Options - psRandom *rng); - -bool ppMergeMaskChipStats (const pmChip *chip, - const pmChip *output, - ppMergeOptions *options, - psRandom *rng); -#endif Index: unk/ppMerge/src/ppMergeMaskAverageConcepts.c =================================================================== --- /trunk/ppMerge/src/ppMergeMaskAverageConcepts.c (revision 17232) +++ (revision ) @@ -1,83 +1,0 @@ -#ifdef HAVE_CONFIG_H -#include -#endif - -#include -#include -#include -#include -#include -#include -#include - -#include "ppMerge.h" -#include "ppMergeData.h" -#include "ppMergeMask.h" - - -// Generate a mask -bool ppMergeMaskAverageConcepts(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 - - pmFPA *fpaOut = data->out; // Output FPA - - pmFPAview *view = pmFPAviewAlloc(0);// View of FPA, for iteration - if (fpaOut->hdu) { - pmFPAUpdateNames(fpaOut, NULL, NULL); - } - - pmChip *chipOut; // Output chip of interest - while ((chipOut = pmFPAviewNextChip(view, fpaOut, 1))) { - - pmCell *cellOut; // Output cell of interest - while ((cellOut = pmFPAviewNextCell(view, fpaOut, 1))) { - - pmHDU *hdu = pmHDUGetLowest(fpaOut, chipOut, cellOut); - if (!hdu || hdu->blankPHU) { - continue; - } - - // Get list of cells for concepts averaging - psList *inCells = psListAlloc(NULL); // List of cells - for (int i = 0; i < filenames->n; i++) { - if (! filenames->data[i] || strlen(filenames->data[i]) == 0) { - continue; - } - pmCell *cellIn = pmFPAviewThisCell(view, data->in->data[i]); // Input cell - psListAdd(inCells, PS_LIST_TAIL, cellIn); - } - if (!pmConceptsAverageCells(cellOut, inCells, NULL, NULL, true)) { - psError(PS_ERR_UNKNOWN, false, "Unable to average cell concepts."); - psFree(inCells); - return false; - } - psFree(inCells); - } - } - - // Get list of FPAs for concepts averaging - psList *inFPAs = psListAlloc(NULL); // List of FPAs - for (int i = 0; i < filenames->n; i++) { - if (! filenames->data[i] || strlen(filenames->data[i]) == 0) { - continue; - } - pmFPA *fpaIn = data->in->data[i]; // Input FPA - psListAdd(inFPAs, PS_LIST_TAIL, fpaIn); - } - - if (!pmConceptsAverageFPAs(fpaOut, inFPAs)) { - psError(PS_ERR_UNKNOWN, false, "Unable to average FPA concepts."); - psFree(inFPAs); - return false; - } - psFree(inFPAs); - - psFree(view); - - return true; -} Index: unk/ppMerge/src/ppMergeMaskBad.c =================================================================== --- /trunk/ppMerge/src/ppMergeMaskBad.c (revision 17232) +++ (revision ) @@ -1,83 +1,0 @@ -#ifdef HAVE_CONFIG_H -#include -#endif - -#include -#include -#include -#include -#include -#include -#include - -#include "ppMerge.h" -#include "ppMergeData.h" -#include "ppMergeMask.h" - - -// Generate a mask -bool ppMergeMaskBad(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 - - pmFPA *fpaOut = data->out; // Output FPA - - pmFPAview *view = pmFPAviewAlloc(0);// View of FPA, for iteration - - pmChip *chipOut; // Output chip of interest - while ((chipOut = pmFPAviewNextChip(view, fpaOut, 1))) { - - pmCell *cellOut; // Output cell of interest - while ((cellOut = pmFPAviewNextCell(view, fpaOut, 1))) { - - psImage *suspect = psMetadataLookupPtr(NULL, cellOut->analysis, "MASK.SUSPECT"); - if (! suspect) { - continue; - } - - pmReadout *roOut = NULL; - if (cellOut->readouts->n > 0) { - roOut = psMemIncrRefCounter (cellOut->readouts->data[0]); - psFree (roOut->mask); - psTrace("ppMerge", 5, "save results with mask from previous pass\n"); - } else { - roOut = pmReadoutAlloc(cellOut); // Output readout - } - - roOut->mask = pmMaskIdentifyBadPixels(suspect, options->combine->maskVal, filenames->n, options->maskBad, options->maskMode); - roOut->data_exists = cellOut->data_exists = chipOut->data_exists = true; - - // Statistics on the merged cell - if (data->statsFile) { - if (!data->stats) { - data->stats = psMetadataAlloc(); - } - - // Build a fake image to do statistics - roOut->image = psImageAlloc(roOut->mask->numCols, roOut->mask->numRows, PS_TYPE_F32); - psImageInit(roOut->image, 1.0); - if (!ppStatsFPA(data->stats, data->out, view, - options->combine->maskVal | pmConfigMask("BLANK", config), - config)) { - psError(PS_ERR_UNEXPECTED_NULL, true, "Unable to generate stats for image.\n"); - return false; - } - psFree(roOut->image); - roOut->image = NULL; - } - psFree(roOut); // Drop reference - - // drop this reference (unless??) - psMetadataRemoveKey (cellOut->analysis, "MASK.SUSPECT"); - } - } - - psFree(view); - - return true; -} - Index: unk/ppMerge/src/ppMergeMaskByImageStats.c =================================================================== --- /trunk/ppMerge/src/ppMergeMaskByImageStats.c (revision 17232) +++ (revision ) @@ -1,87 +1,0 @@ -#ifdef HAVE_CONFIG_H -#include -#endif - -#include -#include -#include -#include -#include -#include -#include - -#include "ppMerge.h" -#include "ppMergeData.h" -#include "ppMergeMask.h" - - -// Generate a mask image consisting of pixels which are outliers from the image mean -bool ppMergeMaskByImageStats(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); - - return true; -} Index: unk/ppMerge/src/ppMergeMaskByPixelStats.c =================================================================== --- /trunk/ppMerge/src/ppMergeMaskByPixelStats.c (revision 17232) +++ (revision ) @@ -1,177 +1,0 @@ -#ifdef HAVE_CONFIG_H -#include -#endif - -#include -#include -#include -#include -#include -#include -#include - -#include "ppMerge.h" -#include "ppMergeData.h" -#include "ppMergeCombine.h" -#include "ppMergeVersion.h" - -// XXX for the moment, this function and ppMergeMaskByImageStats both load (and free) the input -// data - -// Combine the inputs -bool ppMergeMaskByPixelStats(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 - - // Iterate over the FPA - pmFPA *fpa = data->out; // Output FPA - pmFPAview *view = pmFPAviewAlloc(0);// View of FPA, for iteration - int cellNum = -1; // Cell number in the whole FPA - if (data->out->hdu) { - pmFPAUpdateNames(data->out, NULL, NULL); - } - - pmChip *chip; // Chip of interest - pmHDU *lastHDU = NULL; // Last HDU to be updated - - while ((chip = pmFPAviewNextChip(view, fpa, 1))) { - if (chip->hdu) { - // Data will exist soon - pmFPAUpdateNames(data->out, chip, NULL); - chip->data_exists = true; - } - pmCell *cell; // Cell of interest - while ((cell = pmFPAviewNextCell(view, fpa, 1))) { - cellNum++; - if (cell->hdu) { - // Data will exist soon - pmFPAUpdateNames(data->out, chip, cell); - chip->data_exists = cell->data_exists = true; - } - pmReadout *readout = pmReadoutAlloc(cell); // Output readout of interest - psArray *stack = psArrayAlloc(filenames->n); // Stack of readouts to combine - - // Read bit by bit - int numRead; // Number of inputs read - int numScan = 0; - - // Put version metadata into header - pmHDU *hdu = pmHDUFromCell(cell); - if (hdu && hdu != lastHDU) { - if (!hdu->header) { - hdu->header = psMetadataAlloc(); - } - ppMergeVersionMetadata(hdu->header); - lastHDU = hdu; - } - - while (1) { - 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) - bool readOK; - if (pmReadoutReadNext(&readOK, stack->data[i], fits, 0, options->rows)) { - if (!readOK) { - psError(PS_ERR_IO, false, "Failed to read concepts for cell.\n"); - psErrorStackPrint(stderr, "trouble reading data!\n"); - exit (1); - } - 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) { - break; - } - ppMergeMaskReadoutByPixelStats(readout, stack, options->combine); - numScan++; - } - - // Get list of cells for concepts averaging - // XXX only do this operation once between Image and Pixel Stats - psList *inCells = psListAlloc(NULL); // List of cells - for (int i = 0; i < filenames->n; i++) { - if (! filenames->data[i] || strlen(filenames->data[i]) == 0) { - continue; - } - pmCell *cellIn = pmFPAviewThisCell(view, data->in->data[i]); // Input cell - psListAdd(inCells, PS_LIST_TAIL, cellIn); - } - if (!pmConceptsAverageCells(cell, inCells, NULL, NULL, true)) { - psError(PS_ERR_UNKNOWN, false, "Unable to average cell concepts."); - psFree(inCells); - return false; - } - psFree(inCells); - psFree(stack); - psFree(readout); // Drop reference - - // Blow away the cell data - for (int i = 0; i < filenames->n; i++) { - pmFPA *fpaIn = data->in->data[i]; // Input FPA - if (!fpaIn) { continue; } // was not a valid input file - pmChip *chipIn = fpaIn->chips->data[view->chip]; // Input chip - pmCell *cellIn = chipIn->cells->data[view->cell]; // Input cell - pmCellFreeData(cellIn); - } - } - - // Blow away the chip data - for (int i = 0; i < filenames->n; i++) { - pmFPA *fpaIn = data->in->data[i]; // Input FPA - if (!fpaIn) { continue; } // was not a valid input file - pmChip *chipIn = fpaIn->chips->data[view->chip]; // Input chip - pmChipFreeData(chipIn); - } - } - - // Blow away the FPA data - for (int i = 0; i < filenames->n; i++) { - pmFPA *fpaIn = data->in->data[i]; // Input FPA - if (!fpaIn) { continue; } // was not a valid input file - pmFPAFreeData(fpaIn); - } - - psFree(view); - - // Get list of FPAs for concepts averaging - psList *inFPAs = psListAlloc(NULL); // List of FPAs - for (int i = 0; i < filenames->n; i++) { - if (! filenames->data[i] || strlen(filenames->data[i]) == 0) { - continue; - } - pmFPA *fpaIn = data->in->data[i]; // Input FPA - psListAdd(inFPAs, PS_LIST_TAIL, fpaIn); - } - if (!pmConceptsAverageFPAs(fpa, inFPAs)) { - psError(PS_ERR_UNKNOWN, false, "Unable to average FPA concepts."); - psFree(inFPAs); - return false; - } - psFree(inFPAs); - - return true; -} Index: unk/ppMerge/src/ppMergeMaskGrow.c =================================================================== --- /trunk/ppMerge/src/ppMergeMaskGrow.c (revision 17232) +++ (revision ) @@ -1,89 +1,0 @@ -#ifdef HAVE_CONFIG_H -#include -#endif - -#include -#include -#include -#include -#include -#include -#include - -#include "ppMerge.h" -#include "ppMergeData.h" -#include "ppMergeMask.h" - - -// Generate a mask -bool ppMergeMaskGrow(ppMergeData *data, // Data - ppMergeOptions *options, // Options - pmConfig *config // Configuration - ) -{ - pmFPA *fpaOut = data->out; // Output FPA - - pmFPAview *view = pmFPAviewAlloc(0);// View of FPA, for iteration - - pmChip *chipOut; // Output chip of interest - while ((chipOut = pmFPAviewNextChip(view, fpaOut, 1))) { - - pmCell *cellOut; // Output cell of interest - while ((cellOut = pmFPAviewNextCell(view, fpaOut, 1))) { - - pmReadout *roOut; - while ((roOut = pmFPAviewNextReadout(view, fpaOut, 1))) { - - // do this N iterations (XXX add to options) - ppMergeMaskGrowReadout(roOut, options->combine->maskVal, options->growPixels); - ppMergeMaskGrowReadout(roOut, options->combine->maskVal, options->growPixels); - } - } - } - - psFree(view); - - return true; -} - -bool ppMergeMaskGrowReadout (pmReadout *readout, psMaskType maskVal, int nPixels) { - - if (!readout->mask) return true; - - psImage *oldMask = readout->mask; - psImage *newMask = psImageCopy (NULL, oldMask, PS_TYPE_U8); - - psImageInit (newMask, 0); - - int Nx = oldMask->numCols; - int Ny = oldMask->numRows; - - for (int iy = 0; iy < Ny; iy++) { - for (int ix = 0; ix < Nx; ix++) { - int nSum = 0; - for (int jy = -1; jy <= +1; jy++) { - int ny = iy + jy; - if (ny < 0) continue; - if (ny >= Ny) continue; - for (int jx = -1; jx <= +1; jx++) { - if (!jx && !jy) continue; - int nx = ix + jx; - if (nx < 0) continue; - if (nx >= Nx) continue; - if (oldMask->data.U8[ny][nx] & maskVal) { - nSum ++; - } - } - } - if (nSum >= nPixels) { - newMask->data.U8[iy][ix] = maskVal; - } else { - newMask->data.U8[iy][ix] = oldMask->data.U8[iy][ix]; - } - } - } - psFree (readout->mask); - readout->mask = newMask; - - return true; -} Index: unk/ppMerge/src/ppMergeMaskOutput.c =================================================================== --- /trunk/ppMerge/src/ppMergeMaskOutput.c (revision 17232) +++ (revision ) @@ -1,157 +1,0 @@ -#ifdef HAVE_CONFIG_H -#include -#endif - -#include -#include -#include -#include -#include -#include -#include - -#include "ppMerge.h" -#include "ppMergeData.h" -#include "ppMergeMask.h" - - -// write out the mask: assumes it has been saved on the cell->analysis as MASK.UNION -bool ppMergeMaskOutput(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 - - pmFPA *fpaOut = data->out; // Output FPA - - pmFPAview *view = pmFPAviewAlloc(0);// View of FPA, for iteration - if (fpaOut->hdu) { - pmFPAUpdateNames(fpaOut, NULL, NULL); - } - pmFPAWriteMask(fpaOut, data->outFile, config->database, true, false); // Write header only - pmChip *chipOut; // Output chip of interest - while ((chipOut = pmFPAviewNextChip(view, fpaOut, 1))) { - if (chipOut->hdu) { - chipOut->data_exists = true; - pmFPAUpdateNames(fpaOut, chipOut, NULL); - } - pmChipWriteMask(chipOut, data->outFile, config->database, true, false); // Write header only - pmCell *cellOut; // Output cell of interest - while ((cellOut = pmFPAviewNextCell(view, fpaOut, 1))) { - if (cellOut->hdu) { - chipOut->data_exists = cellOut->data_exists = true; - pmFPAUpdateNames(fpaOut, chipOut, cellOut); - } - pmCellWriteMask(cellOut, data->outFile, config->database, true); // Write header only - - psImage *maskUnion = psMetadataLookupPtr(NULL, cellOut->analysis, "MASK.UNION"); - if (! maskUnion) { - continue; - } - - pmReadout *roOut = pmReadoutAlloc(cellOut); // Output readout - roOut->mask = maskUnion; - roOut->data_exists = cellOut->data_exists = chipOut->data_exists = true; - - // XXX skip this? Get list of cells for concepts averaging - { - psList *inCells = psListAlloc(NULL); // List of cells - for (int i = 0; i < filenames->n; i++) { - if (! filenames->data[i] || strlen(filenames->data[i]) == 0) { - continue; - } - pmCell *cellIn = pmFPAviewThisCell(view, data->in->data[i]); // Input cell - psListAdd(inCells, PS_LIST_TAIL, cellIn); - } - if (!pmConceptsAverageCells(cellOut, inCells, NULL, NULL, true)) { - psError(PS_ERR_UNKNOWN, false, "Unable to average cell concepts."); - psFree(inCells); - return false; - } - psFree(inCells); - } - - { - // Add MD5 information for cell - pmHDU *hdu = pmHDUFromCell(cell); // HDU that owns the cell - const char *chipName = psMetadataLookupStr(NULL, chip->concepts, "CHIP.NAME"); - const char *cellName = psMetadataLookupStr(NULL, cell->concepts, "CELL.NAME"); - - psString headerName = NULL; // Header name for MD5 - psStringAppend(&headerName, "MD5_%s_%s", chipName, cellName); - - psVector *md5 = psImageMD5(readout->image); // md5 hash - psString md5string = psMD5toString(md5); // String - psFree(md5); - psMetadataAddStr(hdu->header, PS_LIST_TAIL, headerName, PS_META_REPLACE, - "Image MD5", md5string); - psFree(md5string); - psFree(headerName); - } - - // Statistics on the merged cell - if (data->statsFile) { - if (!data->stats) { - data->stats = psMetadataAlloc(); - } - - // Build a fake image to do statistics - roOut->image = psImageAlloc(roOut->mask->numCols, roOut->mask->numRows, PS_TYPE_F32); - psImageInit(roOut->image, 1.0); - if (!ppStatsFPA(data->stats, data->out, view, - options->combine->maskVal | pmConfigMask("BLANK", config), - config)) { - psError(PS_ERR_UNEXPECTED_NULL, true, "Unable to generate stats for image.\n"); - return false; - } - psFree(roOut->image); - roOut->image = NULL; - } - - psFree(roOut); // Drop reference - - if (cellOut->hdu && !cellOut->hdu->blankPHU) { - psTrace("ppMerge", 5, "Writing out cell HDU.\n"); - pmCellWriteMask(cellOut, data->outFile, config->database, false); - pmCellFreeData(cellOut); - } - } - - if (chipOut->hdu && !chipOut->hdu->blankPHU) { - psTrace("ppMerge", 5, "Writing out chip HDU.\n"); - pmChipWriteMask(chipOut, data->outFile, config->database, false, false); - pmChipFreeData(chipOut); - } - } - - // Get list of FPAs for concepts averaging - { - psList *inFPAs = psListAlloc(NULL); // List of FPAs - for (int i = 0; i < filenames->n; i++) { - if (! filenames->data[i] || strlen(filenames->data[i]) == 0) { - continue; - } - pmFPA *fpaIn = data->in->data[i]; // Input FPA - psListAdd(inFPAs, PS_LIST_TAIL, fpaIn); - } - - if (!pmConceptsAverageFPAs(fpaOut, inFPAs)) { - psError(PS_ERR_UNKNOWN, false, "Unable to average FPA concepts."); - psFree(inFPAs); - return false; - } - psFree(inFPAs); - } - - if (fpaOut->hdu && !fpaOut->hdu->blankPHU) { - psTrace("ppMerge", 5, "Writing out FPA HDU.\n"); - pmFPAWriteMask(fpaOut, data->outFile, config->database, false, false); - } - pmFPAFreeData(fpaOut); - - psFree(view); - - return true; -} Index: unk/ppMerge/src/ppMergeMaskReadoutByPixelStats.c =================================================================== --- /trunk/ppMerge/src/ppMergeMaskReadoutByPixelStats.c (revision 17232) +++ (revision ) @@ -1,327 +1,0 @@ -#ifdef HAVE_CONFIG_H -#include -#endif - -#include -#include -#include -#include -#include -#include - -// XXX: Maybe add support for S16 and S32 types. Currently, only F32 supported. -bool ppMergeMaskReadoutByPixelStats(pmReadout *output, const psArray *inputs, const pmCombineParams *params) -{ - // Check inputs - PS_ASSERT_PTR_NON_NULL(output, false); - PS_ASSERT_ARRAY_NON_NULL(inputs, false); - PS_ASSERT_PTR_NON_NULL(params, false); - PS_ASSERT_FLOAT_WITHIN_RANGE(params->fracLow, 0.0, 1.0, false); - PS_ASSERT_FLOAT_WITHIN_RANGE(params->fracHigh, 0.0, 1.0, false); - if (params->combine != PS_STAT_SAMPLE_MEAN && params->combine != PS_STAT_SAMPLE_MEDIAN && - params->combine != PS_STAT_ROBUST_MEDIAN && params->combine != PS_STAT_FITTED_MEAN && - params->combine != PS_STAT_CLIPPED_MEAN) { - psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Combination method is not SAMPLE_MEAN, SAMPLE_MEDIAN, " - "ROBUST_MEDIAN, FITTED_MEAN or CLIPPED_MEAN.\n"); - return false; - } - - // XXX is it possible / desired to use the weight in this analysis? - for (int i = 0; i < inputs->n; i++) { - pmReadout *readout = inputs->data[i]; // Readout of interest - if (params->weights && !readout->weight) { - psError(PS_ERR_UNEXPECTED_NULL, true, - "Rejection based on weights requested, but no weights supplied for image %d.\n", i); - return false; - } - } - - bool first = !output->image; // First pass through? - - pmHDU *hdu = pmHDUFromReadout(output); // Output HDU - if (!hdu) { - psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to find HDU for readout.\n"); - return false; - } - - if (first) { - psString comment = NULL; // Comment to add to header - psStringAppend(&comment, "Combining using statistic: %x", params->combine); - psMetadataAddStr(hdu->header, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK, comment, ""); - psFree(comment); - } - - psStats *stats = psStatsAlloc(params->combine); // The statistics to use in the combination - if (params->combine == PS_STAT_CLIPPED_MEAN) { - stats->clipSigma = params->rej; - stats->clipIter = params->iter; - - if (first) { - psString comment = NULL; // Comment to add to header - psStringAppend(&comment, "Combination clipping: %d iterations, rejection at %f sigma", - params->iter, params->rej); - psMetadataAddStr(hdu->header, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK, comment, ""); - psFree(comment); - } - } - - int minInputCols, maxInputCols, minInputRows, maxInputRows; // Smallest and largest values to combine - int xSize, ySize; // Size of the output image - if (!pmReadoutStackValidate(&minInputCols, &maxInputCols, &minInputRows, &maxInputRows, &xSize, &ySize, - inputs)) { - psError(PS_ERR_UNKNOWN, false, "No valid input readouts."); - return false; - } - - pmReadoutUpdateSize(output, minInputCols, minInputRows, xSize, ySize, true); - psTrace("psModules.imcombine", 7, "Output minimum: %d,%d\n", output->col0, output->row0); - - psStatsOptions combineStdev = 0; // Statistics option for weights - if (params->weights) { - - if (!output->weight) { - output->weight = psImageAlloc(xSize, ySize, PS_TYPE_F32); - } - if (output->weight->numCols < xSize || output->weight->numRows < ySize) { - psImage *newWeight = psImageAlloc(xSize, ySize, PS_TYPE_F32); - psImageInit(newWeight, 0.0); - psImageOverlaySection(newWeight, output->weight, output->col0, output->row0, "="); - psFree(output->weight); - output->weight = newWeight; - } - - if (first) { - psMetadataAddStr(hdu->header, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK, - "Using input weights to combine images", ""); - } - - // Get the correct statistics option for weights - switch (params->combine) { - case PS_STAT_SAMPLE_MEAN: - case PS_STAT_SAMPLE_MEDIAN: - combineStdev = PS_STAT_SAMPLE_STDEV; - break; - case PS_STAT_ROBUST_MEDIAN: - combineStdev = PS_STAT_ROBUST_STDEV; - break; - case PS_STAT_FITTED_MEAN: - combineStdev = PS_STAT_FITTED_STDEV; - break; - case PS_STAT_CLIPPED_MEAN: - combineStdev = PS_STAT_CLIPPED_STDEV; - break; - default: - psAbort("Should never get here --- checked params->combine before.\n"); - } - stats->options |= combineStdev; - } - - // We loop through each pixel in the output image. We loop through each input readout. We determine if - // that output pixel is contained in the image from that readout. If so, we save it in psVector pixels. - // If not, we set a mask for that element in pixels. Then, we mask off pixels not between fracLow and - // fracHigh. Then we call the vector stats routine on those pixels/mask. Then we set the output pixel - // value to the result of the stats call. - - psVector *pixels = psVectorAlloc(inputs->n, PS_TYPE_F32); // Stack of pixels - psF32 *pixelsData = pixels->data.F32; // Dereference pixels - - psVector *mask = psVectorAlloc(inputs->n, PS_TYPE_U8); // Mask for stack - psU8 *maskData = mask->data.U8; // Dereference mask - - psVector *weights = NULL; // Stack of weights - psVector *errors = NULL; // Stack of errors (sqrt of variance/weights), for psVectorStats - psF32 *weightsData = NULL; // Dereference weights - if (params->weights) { - weights = psVectorAlloc(inputs->n, PS_TYPE_F32); // Stack of weights - weightsData = weights->data.F32; - } - psVector *index = NULL; // The indices to sort the pixels - - float keepFrac = 1.0 - params->fracLow - params->fracHigh; // Fraction of pixels to keep - if (keepFrac != 1.0 && first) { - psString comment = NULL; // Comment to add to header - psStringAppend(&comment, "Min/max rejection: %f high, %f low, keep %d", - params->fracHigh, params->fracLow, params->nKeep); - psMetadataAddStr(hdu->header, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK, comment, ""); - psFree(comment); - } - - psMaskType maskVal = params->maskVal; // The mask value - if (maskVal && first) { - psString comment = NULL; // Comment to add to header - psStringAppend(&comment, "Mask for combination: %x", maskVal); - psMetadataAddStr(hdu->header, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK, comment, ""); - psFree(comment); - } - - #ifndef PS_NO_TRACE - psTrace("psModules.imcombine", 3, "Iterating output: %d --> %d, %d --> %d\n", - minInputCols - output->col0, maxInputCols - output->col0, - minInputRows - output->row0, maxInputRows - output->row0); - if (psTraceGetLevel("psModules.imcombine") >= 3) { - for (int r = 0; r < inputs->n; r++) { - pmReadout *readout = inputs->data[r]; // Input readout - psTrace("psModules.imcombine", 3, "Iterating input %d: %d --> %d, %d --> %d\n", r, - minInputCols - readout->col0, maxInputCols - readout->col0, - minInputRows - readout->row0, maxInputRows - readout->row0); - } - } - #endif - - // Dereference output products - psF32 **outputImage = output->image->data.F32; // Output image - psU8 **outputMask = output->mask->data.U8; // Output mask - psF32 **outputWeight = NULL; // Output weight map - if (output->weight) { - outputWeight = output->weight->data.F32; - } - - psVector *invScale = NULL; // Inverse scale; pre-calculated for efficiency - if (scale) { - invScale = (psVector*)psBinaryOp(NULL, psScalarAlloc(1.0, PS_TYPE_F32), "/", (const psPtr)scale); - } - - for (int i = minInputRows; i < maxInputRows; i++) { - int yOut = i - output->row0; // y position on output readout - #ifdef SHOW_BUSY - - if (psTraceGetLevel("psModules.imcombine") > 9) { - printf("Processing row %d\r", i); - fflush(stdout); - } - #endif - for (int j = minInputCols; j < maxInputCols; j++) { - int xOut = j - output->col0; // x position on output readout - - int numValid = 0; // Number of valid pixels in the stack - memset(maskData, 0, mask->n * sizeof(psU8)); // Reset the mask - for (int r = 0; r < inputs->n; r++) { - pmReadout *readout = inputs->data[r]; // Input readout - int yIn = i - readout->row0; // y position on input readout - int xIn = j - readout->col0; // x position on input readout - psImage *image = readout->image; // The readout image - - #if 0 // This should have been taken care of already: - // Check bounds - if (xIn < 0 || xIn >= image->numCols || yIn < 0 || yIn >= image->numRows) { - continue; - } - #endif - - pixelsData[r] = image->data.F32[yIn][xIn]; - if (!isfinite(pixelsData[r])) { - maskData[r] = 1; - continue; - } - - // Check mask - psImage *roMask = readout->mask; // The mask image - if (roMask && roMask->data.U8[yIn][xIn] & maskVal) { - maskData[r] = 1; - continue; - } - - if (params->weights) { - weightsData[r] = readout->weight->data.F32[yIn][xIn]; - } - - if (zero) { - pixelsData[r] -= zero->data.F32[r]; - } - if (scale) { - pixelsData[r] *= invScale->data.F32[r]; - if (params->weights) { - weightsData[r] *= invScale->data.F32[r] * invScale->data.F32[r]; - } - } - - numValid++; - } - - if (numValid == 0) { - outputMask[yOut][xOut] = params->blank; - outputImage[yOut][xOut] = NAN; - continue; - } - - // Apply fracLow,fracHigh if there are enough pixels - if (numValid * keepFrac >= params->nKeep && keepFrac != 1.0) { - index = psVectorSortIndex(index, pixels); - int numLow = numValid * params->fracLow; // Number of low pixels to clip - int numHigh = numValid * params->fracHigh; // Number of high pixels to clip - // Low pixels - psS32 *indexData = index->data.S32; // Dereference index - for (int k = 0, numMasked = 0; numMasked < numLow && k < index->n; k++) { - // Don't count the ones that are already masked - if (!maskData[indexData[k]]) { - maskData[indexData[k]] = 1; - numMasked++; - } - } - // High pixels - for (int k = pixels->n - 1, numMasked = 0; numMasked < numHigh && k >= 0; k--) { - // Don't count the ones that are already masked - if (! maskData[indexData[k]]) { - maskData[indexData[k]] = 1; - numMasked++; - } - } - } - - // XXXXX this step probably is very expensive : convert errors to variance everywhere? - if (params->weights) { - errors = (psVector*)psUnaryOp(errors, weights, "sqrt"); - } - - // Combination - if (!psVectorStats(stats, pixels, errors, mask, 1)) { - // Can't do much about it, but it's not worth worrying about - psErrorClear(); - outputImage[yOut][xOut] = NAN; - outputMask[yOut][xOut] = params->blank; - if (params->weights) { - outputWeight[yOut][xOut] = NAN; - } - } else { - outputImage[yOut][xOut] = psStatsGetValue(stats, params->combine); - if (!isfinite(outputImage[yOut][xOut])) { - outputMask[yOut][xOut] = params->blank; - } - if (params->weights) { - float stdev = psStatsGetValue(stats, combineStdev); - outputWeight[yOut][xOut] = PS_SQR(stdev); // Variance - // XXXX this is not the correct formal error. - // also, the weighted mean is not obviously the correct thing here - } - } - } - } - #ifdef SHOW_BUSY - if (psTraceGetLevel("psModules.imcombine") > 9) { - printf("\n"); - } - #endif - psFree(index); - psFree(pixels); - psFree(mask); - psFree(weights); - psFree(errors); - psFree(stats); - psFree(invScale); - - // Update the "concepts" - psList *inputCells = psListAlloc(NULL); // List of cells - for (long i = 0; i < inputs->n; i++) { - pmReadout *readout = inputs->data[i]; // Readout of interest - psListAdd(inputCells, PS_LIST_TAIL, readout->parent); - } - bool success = pmConceptsAverageCells(output->parent, inputCells, NULL, NULL, true); - psFree(inputCells); - - output->data_exists = true; - output->parent->data_exists = true; - output->parent->parent->data_exists = true; - - return success; -} - Index: unk/ppMerge/src/ppMergeMaskSuspect.c =================================================================== --- /trunk/ppMerge/src/ppMergeMaskSuspect.c (revision 17232) +++ (revision ) @@ -1,117 +1,0 @@ -#ifdef HAVE_CONFIG_H -#include -#endif - -#include -#include -#include -#include -#include -#include -#include - -#include "ppMerge.h" -#include "ppMergeData.h" -#include "ppMergeMask.h" - - -// Generate a mask -bool ppMergeMaskSuspect(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))) { - - // 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; - } - - // calculate the readout statistics either for each readout, or across the entire chip - if (options->statsByChip) { - pmChip *chipOut = pmFPAviewThisChip(view, fpaOut); // Output cell - if (!ppMergeMaskChipStats (chipIn, chipOut, options, rng)) { - psAbort ("stats problem"); - } - } else { - // calculate the stats for each cell independently - while ((cellIn = pmFPAviewNextCell(view, fpaIn, 1))) { - if (cellIn->readouts->n == 0) continue; - pmCell *cellOut = pmFPAviewThisCell(view, fpaOut); // Output cell - - pmReadout *roOut = NULL; - if (cellOut->readouts->n > 0) { - roOut = cellOut->readouts->data[0]; - psTrace("ppMerge", 5, "masking with results from previous pass\n"); - } - - 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, roOut, 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); - } - pmChipFreeData(chipIn); - } - pmFPAFreeData(fpaIn); - psFree(view); - } - psFree(rng); - - return true; -} Index: unk/ppMerge/src/ppMergeMaskWrite.c =================================================================== --- /trunk/ppMerge/src/ppMergeMaskWrite.c (revision 17232) +++ (revision ) @@ -1,69 +1,0 @@ -#ifdef HAVE_CONFIG_H -#include -#endif - -#include -#include -#include -#include -#include -#include -#include - -#include "ppMerge.h" -#include "ppMergeData.h" -#include "ppMergeMask.h" - - -// Generate a mask -bool ppMergeMaskWrite(ppMergeData *data, // Data - ppMergeOptions *options, // Options - pmConfig *config // Configuration - ) -{ - pmFPA *fpaOut = data->out; // Output FPA - - pmFPAview *view = pmFPAviewAlloc(0);// View of FPA, for iteration - if (fpaOut->hdu) { - pmFPAUpdateNames(fpaOut, NULL, NULL); - } - pmFPAWriteMask(fpaOut, data->outFile, config->database, true, false); // Write header only - pmChip *chipOut; // Output chip of interest - while ((chipOut = pmFPAviewNextChip(view, fpaOut, 1))) { - if (chipOut->hdu) { - chipOut->data_exists = true; - pmFPAUpdateNames(fpaOut, chipOut, NULL); - } - pmChipWriteMask(chipOut, data->outFile, config->database, true, false); // Write header only - pmCell *cellOut; // Output cell of interest - while ((cellOut = pmFPAviewNextCell(view, fpaOut, 1))) { - if (cellOut->hdu) { - chipOut->data_exists = cellOut->data_exists = true; - pmFPAUpdateNames(fpaOut, chipOut, cellOut); - } - pmCellWriteMask(cellOut, data->outFile, config->database, true); // Write header only - - if (cellOut->hdu && !cellOut->hdu->blankPHU) { - psTrace("ppMerge", 5, "Writing out cell HDU.\n"); - pmCellWriteMask(cellOut, data->outFile, config->database, false); - pmCellFreeData(cellOut); - } - } - - if (chipOut->hdu && !chipOut->hdu->blankPHU) { - psTrace("ppMerge", 5, "Writing out chip HDU.\n"); - pmChipWriteMask(chipOut, data->outFile, config->database, false, false); - pmChipFreeData(chipOut); - } - } - - if (fpaOut->hdu && !fpaOut->hdu->blankPHU) { - psTrace("ppMerge", 5, "Writing out FPA HDU.\n"); - pmFPAWriteMask(fpaOut, data->outFile, config->database, false, false); - } - pmFPAFreeData(fpaOut); - - psFree(view); - - return true; -} Index: unk/ppMerge/src/ppMergeOptions.c =================================================================== --- /trunk/ppMerge/src/ppMergeOptions.c (revision 17232) +++ (revision ) @@ -1,323 +1,0 @@ -#ifdef HAVE_CONFIG_H -#include -#endif - -#include -#include -#include -#include -#include - -#include "ppMerge.h" -#include "ppMergeOptions.h" - -#define ARRAY_BUFFER 4 // Number of values to add at once - -////////////////////////////////////////////////////////////////////////////////////////////////////////////// -// ppMergeOptions -////////////////////////////////////////////////////////////////////////////////////////////////////////////// - -// Free function -static void mergeOptionsFree(ppMergeOptions *options // Options to free - ) -{ - psFree(options->format); - psFree(options->combine); - psFree(options->darkOrdinates); - psFree(options->darkNorm); -} - -// Allocator -ppMergeOptions *ppMergeOptionsAlloc(const pmConfig *config) -{ - ppMergeOptions *options = psAlloc(sizeof(ppMergeOptions)); // The options, to return - psMemSetDeallocator(options, (psFreeFunc)mergeOptionsFree); - - options->format = NULL; - options->rows = 0; - options->minElectrons = NAN; - options->zero = false; - options->scale = false; - options->dark = false; - options->fringe = false; - options->shutter = false; - options->mask = false; - options->sample = 1; - options->mean = PS_STAT_SAMPLE_MEDIAN; - options->stdev = PS_STAT_SAMPLE_STDEV; - options->fringeNum = 100; - options->fringeSize = 10; - options->fringeSmoothX = 5; - options->fringeSmoothY = 5; - options->shutterSize = 10; - options->shutterIter = 2; - options->shutterRej = 3.0; - options->maskSuspect = 5.0; - options->maskBad = 10.0; - options->maskMode = PM_MASK_ID_VALUE; - options->statsByChip = true; - options->onOff = 0; - options->combine = pmCombineParamsAlloc(PS_STAT_SAMPLE_MEAN); - options->combine->rej = 3.0; - options->combine->iter = 1; - options->combine->fracHigh = 0.0; - options->combine->fracLow = 0.0; - options->combine->nKeep = 1; - options->combine->maskVal = 0xff; - options->combine->weights = false; - options->combine->blank = pmConfigMask("BLANK", config); - options->satMask = 0x00; - options->badMask = 0x00; - options->growPixels = 0; - options->darkOrdinates = NULL; - options->darkNorm = NULL; - return options; -} - -////////////////////////////////////////////////////////////////////////////////////////////////////////////// -// ppMergeOptionsParse -////////////////////////////////////////////////////////////////////////////////////////////////////////////// - -// Parse a recipe option according to its type -#define OPTION_PARSE(OPTION,MD,NAME,TYPE) \ -{ \ - psMetadataItem *item = psMetadataLookup(MD, NAME); \ - if (item) { \ - OPTION = psMetadataItemParse##TYPE(item); \ - } else { \ - psWarning("Recipe option %s isn't specified; using default.\n", NAME); \ - } \ -} - -// Parse a statistic -static psStatsOptions parseStat(psMetadata *source, // Source of the statistics option - const char *name // Name of the statistics option - ) -{ - bool mdok = true; // Status of MD lookup - const char *stat = psMetadataLookupStr(&mdok, source, name); // The statistic string - if (!mdok || !stat || strlen(stat) == 0) { - return 0; - } - return psStatsOptionFromString(stat); -} - -// Parse the options -ppMergeOptions *ppMergeOptionsParse(pmConfig *config // Configuration - ) -{ - ppMergeOptions *options = ppMergeOptionsAlloc(config); // The merge options - - // We need to work out the camera before we can get the recipe. Take the first input and inspect it. - if (!config->camera) { - psArray *filenames = psMetadataLookupPtr(NULL, config->arguments, "INPUT"); // The input file names - psString resolved = pmConfigConvertFilename(filenames->data[0], config, false); // Resolved file name - psFits *inFile = psFitsOpen(resolved, "r"); // The FITS file to read - if (!inFile) { - psError(PS_ERR_IO, false, "Unable to open input file %s to determine camera.\n", resolved); - psFree(resolved); - psFree(config); - exit(EXIT_FAILURE); - } - psFree(resolved); - psMetadata *header = psFitsReadHeader(NULL, inFile); // The FITS (primary) header - psFitsClose(inFile); - - options->format = pmConfigCameraFormatFromHeader(config, header, true); - psFree(header); - if (!options->format) { - psLogMsg(__func__, PS_LOG_WARN, "Unable to identify camera format for input file %s\n", - (char *)filenames->data[0]); - exit(EXIT_FAILURE); - } - } - - // we must have the recipes by this point - assert (config->recipes); - - // Now we can read the recipe - bool mdok = true; // Status of MD lookup - psMetadata *recipe = psMetadataLookupMetadata(&mdok, config->recipes, PPMERGE_RECIPE); // Recipe information - if (!mdok || !recipe) { - psError(PS_ERR_IO, true, "Unable to find recipe %s", PPMERGE_RECIPE); - exit(EXIT_FAILURE); - } - - // First, deal with the recipe. These are parameters that will typically be constant for a camera. - OPTION_PARSE(options->rows, recipe, "ROWS", S32); - OPTION_PARSE(options->minElectrons, recipe, "ELECTRONS", F32); - OPTION_PARSE(options->sample, recipe, "SAMPLE", S32); - OPTION_PARSE(options->combine->rej, recipe, "REJ", F32); - OPTION_PARSE(options->combine->iter, recipe, "ITER", S32); - OPTION_PARSE(options->combine->fracHigh, recipe, "FRACHIGH", F32); - OPTION_PARSE(options->combine->fracLow, recipe, "FRACLOW", F32); - OPTION_PARSE(options->combine->nKeep, recipe, "NKEEP", S32); - OPTION_PARSE(options->combine->weights, recipe, "WEIGHTS", Bool); - OPTION_PARSE(options->fringeNum, recipe, "FRINGE.NUM", S32); - OPTION_PARSE(options->fringeSize, recipe, "FRINGE.SIZE", S32); - OPTION_PARSE(options->fringeSmoothX, recipe, "FRINGE.XSMOOTH", S32); - OPTION_PARSE(options->fringeSmoothY, recipe, "FRINGE.YSMOOTH", S32); - OPTION_PARSE(options->shutterSize, recipe, "SHUTTER.SIZE", S32); - 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); - OPTION_PARSE(options->statsByChip, recipe, "STATS.BY.CHIP", Bool); - OPTION_PARSE(options->growPixels, recipe, "MASK.GROW.NPIX", S32); - - const char *masks = psMetadataLookupStr(NULL, recipe, "MASKVAL"); - options->combine->maskVal = pmConfigMask(masks, config); - - options->combine->combine = parseStat(recipe, "COMBINE"); - options->mean = parseStat(recipe, "MEAN"); - options->stdev = parseStat(recipe, "STDEV"); - - // Now the command-line options. These are parameters that depend on what type of frame is being combined - - // Set options based on the type of calibration frame - const char *type = psMetadataLookupStr(NULL, config->arguments, "-type"); // The type of calibration frame - if (strlen(type) > 0) { - if (strcasecmp(type, "BIAS") == 0) { - options->zero = false; - options->scale = false; - options->dark = false; - options->fringe = false; - options->shutter = false; - options->mask = false; - } else if (strcasecmp(type, "DARK") == 0) { - options->zero = false; - options->scale = false; - options->dark = true; - options->fringe = false; - options->shutter = false; - options->mask = false; - } else if (strcasecmp(type, "FLAT") == 0) { - options->zero = false; - options->scale = true; - options->dark = false; - options->fringe = false; - options->shutter = false; - options->mask = false; - } else if (strcasecmp(type, "SKYFLAT") == 0) { - options->zero = false; - options->scale = true; - options->dark = false; - options->fringe = false; - options->shutter = false; - options->mask = false; - } else if (strcasecmp(type, "DOMEFLAT") == 0) { - options->zero = false; - options->scale = true; - options->dark = false; - options->fringe = false; - options->shutter = false; - options->mask = false; - } else if (strcasecmp(type, "FRINGE") == 0) { - options->zero = true; - options->scale = true; - options->dark = false; - options->fringe = true; - options->shutter = false; - options->mask = false; - } else if (strcasecmp(type, "SHUTTER") == 0) { - options->zero = false; - options->scale = false; - options->dark = false; - options->fringe = false; - options->shutter = true; - options->mask = false; - } else if (strcasecmp(type, "MASK") == 0 || - strcasecmp(type, "DARKMASK") == 0 || - strcasecmp(type, "FLATMASK") == 0) { - options->zero = false; - options->scale = false; - options->dark = false; - options->fringe = false; - options->shutter = false; - options->mask = true; - } else { - psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Unrecognised image type: %s", type); - psFree(options); - return NULL; - } - } else { - psError(PS_ERR_BAD_PARAMETER_VALUE, true, "No image type specified."); - psFree(options); - return NULL; - } - - const char *maskMode = psMetadataLookupStr(NULL, recipe, "MASK.MODE"); // The type of calibration frame - if (!maskMode) { - psError(PS_ERR_BAD_PARAMETER_VALUE, true, "No mask mode specified in recipe."); - psFree(options); - return NULL; - } - - options->maskMode = pmMaskIdentifyModeFromString (maskMode); - if (options->maskMode == PM_MASK_ID_NONE) { - psError(PS_ERR_BAD_PARAMETER_VALUE, true, "invalid mask mode %s.", maskMode); - psFree(options); - return NULL; - } - -#if 0 - // Number of on/off images - OPTION_PARSE(options->onOff, config->arguments, "-onoff", S32); -#endif - - // Masking options - options->satMask = pmConfigMask("SAT", config); - options->badMask = pmConfigMask("BAD", config); - - if (options->dark) { - psMetadata *ordinates = psMetadataLookupMetadata(NULL, recipe, "DARK.ORDINATES"); // Ordinates info - options->darkOrdinates = psArrayAllocEmpty(psListLength(ordinates->list)); - - psMetadataIterator *iter = psMetadataIteratorAlloc(ordinates, PS_LIST_HEAD, NULL); // Iterator - psMetadataItem *item; // Item from iteration - while ((item = psMetadataGetAndIncrement(iter))) { - int order = 0; // Polynomial order - bool scale = false; // Scale values? - float min = NAN, max = NAN; // Minimum and maximum values for scaling - switch (item->type) { - case PS_TYPE_S32: - order = item->data.S32; - break; - case PS_DATA_METADATA: - order = psMetadataLookupS32(NULL, item->data.md, "ORDER"); - bool mdok; // Status of MD lookup - scale = psMetadataLookupBool(&mdok, item->data.md, "SCALE"); - min = psMetadataLookupF32(&mdok, item->data.md, "MIN"); - max = psMetadataLookupF32(&mdok, item->data.md, "MAX"); - break; - default: - psError(PS_ERR_BAD_PARAMETER_TYPE, true, - "Type of DARK.ORDINATES entry %s (%x) is not METADATA or S32", - item->name, item->type); - psFree(options); - return false; - } - if (order <= 0) { - psError(PS_ERR_BAD_PARAMETER_VALUE, true, "ORDER not positive (%d) for DARK.ORDINATES %s", - order, item->name); - psFree(options); - return false; - } - - pmDarkOrdinate *ord = pmDarkOrdinateAlloc(item->name, order); - ord->scale = scale; - ord->min = min; - ord->max = max; - psArrayAdd(options->darkOrdinates, options->darkOrdinates->n, ord); - psFree(ord); - } - psFree(iter); - - psString darkNorm = psMetadataLookupStr(NULL, recipe, "DARK.NORM"); // Normalisation concept - if (darkNorm && strcmp(darkNorm, "NONE") != 0) { - options->darkNorm = psStringCopy(darkNorm); - } - } - - return options; -} Index: unk/ppMerge/src/ppMergeOptions.h =================================================================== --- /trunk/ppMerge/src/ppMergeOptions.h (revision 17232) +++ (revision ) @@ -1,59 +1,0 @@ -#ifndef PP_MERGE_OPTIONS_H -#define PP_MERGE_OPTIONS_H - -#include -#include - -// Mode of on/off pairs; the value corresponds to how many images are in each set -typedef enum { - PP_ONOFF_ABBA = -1, // On/off pairs in the ABBA mode - PP_ONOFF_NONE = 0, // No on/off pairs - PP_ONOFF_ABAB = 1, // On/off pairs in the ABAB mode (one image each) - PP_ONOFF_AABB = 2, // On/off pairs, two images each - // And so on and so forth... just use a number beyond this -} ppOnOff; - -// Options for ppMerge -typedef struct { - psMetadata *format; // Camera format configuration - unsigned int rows; // Number of rows to read at once - float minElectrons; // Minimum number of electrons for useful signal - bool zero; // Subtract background before combining? - bool scale; // Scale by the background before combining? - bool dark; // Generate dark? - 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 - psStatsOptions stdev; // Statistic to use to measure the stdev - int fringeNum; // Number of fringe regions per cell - int fringeSize; // Size of fringe regions - int fringeSmoothX; // Number of smoothing regions per cell, in x - int fringeSmoothY; // Number of smoothing regions per cell, in y - int shutterSize; // Size for shutter measurement regions - 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 - bool statsByChip; // measure statistics for masking by chip or readout? - pmMaskIdentifyMode maskMode; // how to set the limit based on the threshold value above? - ppOnOff onOff; // On/off pairs? - pmCombineParams *combine; // Combination parameters - psMaskType satMask; // Mask value for saturated pixels - psMaskType badMask; // Mask value for bad (low) pixels - int growPixels; - psArray *darkOrdinates; // Ordinates for dark combination - psString darkNorm; // Dark normalisation concept -} ppMergeOptions; - -// Allocator -ppMergeOptions *ppMergeOptionsAlloc(const pmConfig *config); - - -// Parse the options for ppMerge -ppMergeOptions *ppMergeOptionsParse(pmConfig *config // Configuration - ); - - -#endif