Index: trunk/ppMerge/src/ppMergeBackground.c =================================================================== --- trunk/ppMerge/src/ppMergeBackground.c (revision 7001) +++ trunk/ppMerge/src/ppMergeBackground.c (revision 7061) @@ -5,4 +5,8 @@ #include "ppMergeBackground.h" + +#define MAX_ITER 10 // Maximum number of iterations for pmFlatNormalize +#define TOLERANCE 1e-3 // Tolerance to reach for pmFlatNormalize + // Extract a particular statistic from the stats @@ -29,37 +33,40 @@ -// Get the background for each chip of each input -psImage *ppMergeBackground(const ppMergeOptions *options, // The options - const pmConfig *config // The configuration +// Get the scale and zero for each chip of each input +bool ppMergeScaleZero(ppImage **scales, // The scales for each integration/cell + ppImage **zeros, // The zeroes for each integration/cell + ppMergeData *data,// The data + const ppMergeOptions *options, // The options + const pmConfig *config // The configuration ) { assert(config->camera); // Need the camera configuration - int numCells = 0; // Number of cells in the model FPA - { - // Count the cells - pmFPA *fpa = pmFPAConstruct(config->camera); // The FPA - psArray *chips = fpa->chips; // Array of chips - 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++; - } - } - } - psFree(fpa); - } + psArray *filenames = psMetadataLookupPtr(NULL, config->arguments, "INPUT"); // The input file names assert(filenames); // It should be here --- it's put here in ppMergeConfig - psImage *background = psImageAlloc(filenames->n, numCells, PS_TYPE_F64); // Background measurements + // Sanity checks + assert(!options->scale || scales); + assert(!scales || !*scales || ((*scales)->type.type == PS_TYPE_F64 && + (*scales)->numCols == data->numCells && + (*scales)->numRows == filenames->n)); + assert(!options->zero || zeros); + assert(!zeros || !*zeros || ((*zeros)->type.type == PS_TYPE_F64 && + (*zeros)->numCols == data->numCells && + (*zeros)->numRows == filenames->n)); + + psImage *background = psImageAlloc(filenames->n, data->numCells, PS_TYPE_F64); // Background measurements psImageInit(background, NAN); psStats *bgStats = psStatsAlloc(options->background); // Statistic to measure the background - + psVector *gains = NULL; // The gains for each cell + psImage *exptime = NULL; // The exposure time for each integration of each cell + if (options->scale) { + gains = psVectorAlloc(data->numCells, PS_TYPE_F64); + } + if (options->exptime) { + exptime = psImageAlloc(filenames->n, data->numCells, PS_TYPE_F64); + } + + bool status = true; // Status of getting the scale and zero --- did everything go right? for (int i = 0; i < filenames->n; i++) { psString name = filenames->data[i]; // The name of the file @@ -70,63 +77,131 @@ if (!inFile) { psErrorPrint(PS_ERR_IO, false, "Unable to open input file %s --- ignored.\n", name); + status = false; continue; } - psMetadata *header = psFitsReadHeader(NULL, inFile); // The FITS (primary) header - psMetadata *format = pmConfigCameraFormatFromHeader(config, header); // The camera format - psFree(header); - if (!format) { - psErrorPrint(PS_ERR_IO, false, "Unable to identify camera format for input file %s --- " - "ignored.\n", name); - continue; - } - - pmFPA *fpa = pmFPAConstruct(config->camera); // The FPA - pmFPAview *view = pmFPAAddSourceFromHeader(fpa, phu, format); // View of the PHU - int cellNum = 0; // Number of the cell + + pmFPA *fpa = in->data[i]; // The FPA for this input + int cellNum = -1; // Number of the cell psArray *chips = fpa->chips; // Array of chips - for (int i = 0; i < chips->n; i++) { - pmChip *chip = chips->data[i]; // The chip of interest + for (int j = 0; j < chips->n; j++) { + pmChip *chip = chips->data[j]; // The 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]; // The cell of interest + for (int k = 0; k < cells->n; k++) { + pmCell *cell = cells->data[k]; // The cell of interest if (!cell) { continue; } - - if (!pmCellRead(cell, inFile, config->database)) { - // Nothing here - psFree(cell); - cellNum++; - continue; - } - - if (cell->readouts->n > 1) { - psLogMsg(__func__, PS_LOG_WARN, "File %s chip %d cell %d contains more than one readout " - "--- ignoring all but the first.\n", name, i, j); - } - - pmReadout *readout = cell->readouts->data[0]; // The readout of interest - psImage *image = readout->image; // The pixels of interest - if (!image) { - psFree(cell); - cellNum++; - continue; - } - - psImageStats(bgStats, image, readout->mask, options->maskVal); - background[i][cellNum++] = getStat(bgStats, options->background); - psFree(cell); + cellNum++; + + if (!pmCellReadHeader(cell, inFile)) { + psLogMsg(__func__, PS_LOG_WARN, "Unable to read header for file %s chip %d cell %d --- " + "ignored.\n", name, j, k); + status = false; + } + + // Normalising by the exposure time + if (options->exptime) { + bool mdok = true; // Status of MD lookup + exptime->data.F64[i][cellNum] = psMetadataLookup(&mdok, cell->concepts, "CELL.EXPOSURE"); + if (!mdok || isnan(exptime->data.F64[i][cellNum])) { + psLogMsg(__func__, PS_LOG_WARN, "CELL.EXPOSURE for file %s chip %d cell %d is not " + "set.\n", name, j, k); + exptime->data.F64[i][cellNum] = NAN; + status = false; + } + } + + // Scaling by the background + if (options->scale || options->zero) { + if (!pmCellRead(cell, inFile, config->database)) { + // Nothing here + psFree(cell); + continue; + } + + if (cell->readouts->n > 1) { + psLogMsg(__func__, PS_LOG_WARN, "File %s chip %d cell %d contains more than one " + "readout --- ignoring all but the first.\n", name, j, k); + status = false; + } + + pmReadout *readout = cell->readouts->data[0]; // The readout of interest + psImage *image = readout->image; // The pixels of interest + if (!image) { + psFree(cell); + continue; + } + + // Get the gain + if (options->scale) { + bool mdok = true; // Status of MD lookup + gain->data.F32[cellNum] = psMetadataLookupF32(&mdok, cell->concepts, "CELL.GAIN"); + if (!mdok || isnan(gain->data.F32[cellNum])) { + psLogMsg(__func__, PS_LOG_WARN, "CELL.GAIN for file %s chip %d cell %d is not " + "set.\n", name, j, k); + gain->data.F32[cellNum] = NAN; + status = false; + } + } + + // Get the background + psImageStats(bgStats, image, readout->mask, options->maskVal); + background[i][cellNum] = getStat(bgStats, options->background); + } + + psCellFreeData(cell); } - psFree(chip); - } - psFree(fpa); - psFree(view); + psChipFreeData(chip); + } + psFPAFreeData(fpa); psFitsClose(inFile); } psFree(bgStats); - return background; + if (options->scale) { + // Need to normalize over the focal plane + bool converge = true; // Did we converge? + psVector *fluxes = pmFlatNormalize(&converge, gains, background, MAX_ITER, TOLERANCE); + if (!converge || !fluxes) { + psLogMsg(__func__, PS_LOG_WARN, "Normalisation failed to converge --- continuing anyway.\n"); + status = false; + } + psFree(gains); + + *scales = psImageCopy(*scales, background, PS_TYPE_F64); + psImage *scalesDeref = *scales; // Dereference the pointer + + for (int i = 0; i < scalesDeref->numRows; i++) { + psF64 bg = fluxes->data.F64[i]; + for (int j = 0; j < scalesDeref->numCols; j++) { + scalesDeref->data.F64[i][j] = bg; + } + } + } + + if (options->exptime) { + if (!options->scale) { + // Copy over the exposure times + psImageCopy(*scales, exptime, PS_TYPE_F64); + } else { + psBinaryOp(*scales, *scales, "*", exptime); + } + } + + if (options->zero) { + if (!*zeros) { + // This is much faster than copying! + *zeros = psMemIncrRefCounter(background); + } else { + *zeros = psImageCopy(*zeros, background, PS_TYPE_F64); + } + } + + psFree(background); + + return status; } +