Index: unk/ppMerge/src/ppMergeBackground.c =================================================================== --- /trunk/ppMerge/src/ppMergeBackground.c (revision 7260) +++ (revision ) @@ -1,258 +1,0 @@ -#include -#include -#include -#include - -#include "ppMerge.h" -#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 -static inline double getStat(const psStats *stats, // Statistics structure - psStatsOptions option // Statistic option to return - ) -{ - switch (option) { - case PS_STAT_SAMPLE_MEAN: - return stats->sampleMean; - case PS_STAT_SAMPLE_MEDIAN: - return stats->sampleMedian; - case PS_STAT_ROBUST_MEDIAN: - return stats->robustMedian; - case PS_STAT_FITTED_MEAN: - return stats->fittedMean; - case PS_STAT_CLIPPED_MEAN: - return stats->clippedMean; - default: - psAbort(__func__, "Invalid statistics option: %x\n", option); - } - return NAN; // Token return -} - - -// Get the scale and zero for each chip of each input -bool ppMergeScaleZero(psImage **scales, // The scales for each integration/cell - psImage **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 - - 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 || ((*scales)->type.type == PS_TYPE_F32 && - (*scales)->numCols == data->numCells && - (*scales)->numRows == filenames->n)); - assert(!options->zero || zeros); - assert(!zeros || !*zeros || ((*zeros)->type.type == PS_TYPE_F32 && - (*zeros)->numCols == data->numCells && - (*zeros)->numRows == filenames->n)); - - psImage *background = psImageAlloc(data->numCells, filenames->n, PS_TYPE_F32); // 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_F32); - gains->n = data->numCells; - } - if (options->exptime) { - exptime = psImageAlloc(filenames->n, data->numCells, PS_TYPE_F32); - } - - 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 - if (!name || strlen(name) == 0) { - continue; - } - psTrace(__func__, 9, "Opening %s to get background...\n", name); - psFits *inFile = psFitsOpen(filenames->data[i], "r"); // The FITS file to read - if (!inFile) { - psLogMsg(__func__, PS_LOG_WARN, "Unable to open input file %s --- ignored.\n", name); - status = false; - continue; - } - - pmFPA *fpa = data->in->data[i]; // The FPA for this input - int cellNum = -1; // Number of the cell - psArray *chips = fpa->chips; // Array of chips - 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 k = 0; k < cells->n; k++) { - pmCell *cell = cells->data[k]; // The cell of interest - if (!cell) { - continue; - } - 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.F32[i][cellNum] = psMetadataLookupF32(&mdok, cell->concepts, - "CELL.EXPOSURE"); - if (!mdok || isnan(exptime->data.F32[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.F32[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 - gains->data.F32[cellNum] = psMetadataLookupF32(&mdok, cell->concepts, "CELL.GAIN"); - if (!mdok || isnan(gains->data.F32[cellNum])) { - psLogMsg(__func__, PS_LOG_WARN, "CELL.GAIN for file %s chip %d cell %d is not " - "set.\n", name, j, k); - gains->data.F32[cellNum] = NAN; - status = false; - } - } - - // Get the background - int sampleSize = (image->numCols * image->numRows) / options->sample; // Size of sample - psVector *sample = psVectorAlloc(sampleSize, PS_TYPE_F32); // Sample of the image - sample->n = sampleSize; - psVector *sampleMask = NULL; // Mask for sample - if (readout->mask) { - sampleMask = psVectorAlloc(sampleSize, PS_TYPE_U8); - sampleMask->n = sampleSize; - } - psImage *mask = readout->mask; // The mask image - for (long i = 0; i < sampleSize; i++) { - int j = i * options->sample; // Index into image - int x = j % image->numCols; // x index - int y = j / image->numCols; // y index - sample->data.F32[i] = image->data.F32[y][x]; - if (readout->mask) { - sampleMask->data.U8[i] = mask->data.U8[y][x]; - } - } - psVectorStats(bgStats, sample, sampleMask, NULL, options->combine->maskVal); - psFree(sample); - psFree(sampleMask); - background->data.F32[i][cellNum] = getStat(bgStats, options->background); - psTrace(__func__, 3, "Background for %s, cell %d is %f\n", name, cellNum, - background->data.F32[i][cellNum]); - } - - pmCellFreeData(cell); - } - pmChipFreeData(chip); - } - pmFPAFreeData(fpa); - psFitsClose(inFile); - } - psFree(bgStats); - - if (options->scale) { - // Need to normalize over the focal plane - bool converge = true; // Did we converge? - if (psTraceGetLevel(__func__) > 9) { - for (int i = 0; i < gains->n; i++) { - psTrace(__func__, 10, "Gain for cell %d is %f\n", i, gains->data.F32[i]); - } - } - 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_F32); - psImage *scalesDeref = *scales; // Dereference the pointer - - for (int i = 0; i < scalesDeref->numRows; i++) { - psF32 bg = fluxes->data.F32[i]; - for (int j = 0; j < scalesDeref->numCols; j++) { - scalesDeref->data.F32[i][j] = bg; - } - } - } - - if (options->exptime) { - if (!options->scale) { - // Copy over the exposure times - psImageCopy(*scales, exptime, PS_TYPE_F32); - } 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_F32); - } - } - - // Diagnostic stuff - if (scales && *scales && psTraceGetLevel(__func__) > 9) { - psImage *scalesDeref = *scales; // Dereference the pointer - for (int i = 0; i < scalesDeref->numRows; i++) { - for (int j = 0; j < scalesDeref->numCols; j++) { - psTrace(__func__, 9, "Scale for exposure %d, cell %d is: %f\n", i, j, - scalesDeref->data.F32[i][j]); - } - } - } - if (zeros && *zeros && psTraceGetLevel(__func__) > 9) { - psImage *zerosDeref = *zeros; // Dereference the pointer - for (int i = 0; i < zerosDeref->numRows; i++) { - for (int j = 0; j < zerosDeref->numCols; j++) { - psTrace(__func__, 9, "Zero for exposure %d, cell %d is: %f\n", i, j, - zerosDeref->data.F32[i][j]); - } - } - } - - psFree(background); - - return status; -} - Index: unk/ppMerge/src/ppMergeBackground.h =================================================================== --- /trunk/ppMerge/src/ppMergeBackground.h (revision 7260) +++ (revision ) @@ -1,18 +1,0 @@ -#ifndef PP_MERGE_BACKGROUND_H -#define PP_MERGE_BACKGROUND_H - -#include -#include - -#include "ppMergeData.h" -#include "ppMergeOptions.h" - -// Get the scale and zero for each chip of each input -bool ppMergeScaleZero(psImage **scales, // The scales for each integration/cell - psImage **zeros, // The zeroes for each integration/cell - ppMergeData *data,// The data - const ppMergeOptions *options, // The options - const pmConfig *config // The configuration - ); - -#endif Index: /trunk/ppMerge/src/ppMergeScaleZero.c =================================================================== --- /trunk/ppMerge/src/ppMergeScaleZero.c (revision 7261) +++ /trunk/ppMerge/src/ppMergeScaleZero.c (revision 7261) @@ -0,0 +1,258 @@ +#include +#include +#include +#include + +#include "ppMerge.h" +#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 +static inline double getStat(const psStats *stats, // Statistics structure + psStatsOptions option // Statistic option to return + ) +{ + switch (option) { + case PS_STAT_SAMPLE_MEAN: + return stats->sampleMean; + case PS_STAT_SAMPLE_MEDIAN: + return stats->sampleMedian; + case PS_STAT_ROBUST_MEDIAN: + return stats->robustMedian; + case PS_STAT_FITTED_MEAN: + return stats->fittedMean; + case PS_STAT_CLIPPED_MEAN: + return stats->clippedMean; + default: + psAbort(__func__, "Invalid statistics option: %x\n", option); + } + return NAN; // Token return +} + + +// Get the scale and zero for each chip of each input +bool ppMergeScaleZero(psImage **scales, // The scales for each integration/cell + psImage **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 + + 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 || ((*scales)->type.type == PS_TYPE_F32 && + (*scales)->numCols == data->numCells && + (*scales)->numRows == filenames->n)); + assert(!options->zero || zeros); + assert(!zeros || !*zeros || ((*zeros)->type.type == PS_TYPE_F32 && + (*zeros)->numCols == data->numCells && + (*zeros)->numRows == filenames->n)); + + psImage *background = psImageAlloc(data->numCells, filenames->n, PS_TYPE_F32); // 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_F32); + gains->n = data->numCells; + } + if (options->exptime) { + exptime = psImageAlloc(filenames->n, data->numCells, PS_TYPE_F32); + } + + 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 + if (!name || strlen(name) == 0) { + continue; + } + psTrace(__func__, 9, "Opening %s to get background...\n", name); + psFits *inFile = psFitsOpen(filenames->data[i], "r"); // The FITS file to read + if (!inFile) { + psLogMsg(__func__, PS_LOG_WARN, "Unable to open input file %s --- ignored.\n", name); + status = false; + continue; + } + + pmFPA *fpa = data->in->data[i]; // The FPA for this input + int cellNum = -1; // Number of the cell + psArray *chips = fpa->chips; // Array of chips + 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 k = 0; k < cells->n; k++) { + pmCell *cell = cells->data[k]; // The cell of interest + if (!cell) { + continue; + } + 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.F32[i][cellNum] = psMetadataLookupF32(&mdok, cell->concepts, + "CELL.EXPOSURE"); + if (!mdok || isnan(exptime->data.F32[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.F32[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 + gains->data.F32[cellNum] = psMetadataLookupF32(&mdok, cell->concepts, "CELL.GAIN"); + if (!mdok || isnan(gains->data.F32[cellNum])) { + psLogMsg(__func__, PS_LOG_WARN, "CELL.GAIN for file %s chip %d cell %d is not " + "set.\n", name, j, k); + gains->data.F32[cellNum] = NAN; + status = false; + } + } + + // Get the background + int sampleSize = (image->numCols * image->numRows) / options->sample; // Size of sample + psVector *sample = psVectorAlloc(sampleSize, PS_TYPE_F32); // Sample of the image + sample->n = sampleSize; + psVector *sampleMask = NULL; // Mask for sample + if (readout->mask) { + sampleMask = psVectorAlloc(sampleSize, PS_TYPE_U8); + sampleMask->n = sampleSize; + } + psImage *mask = readout->mask; // The mask image + for (long i = 0; i < sampleSize; i++) { + int j = i * options->sample; // Index into image + int x = j % image->numCols; // x index + int y = j / image->numCols; // y index + sample->data.F32[i] = image->data.F32[y][x]; + if (readout->mask) { + sampleMask->data.U8[i] = mask->data.U8[y][x]; + } + } + psVectorStats(bgStats, sample, sampleMask, NULL, options->combine->maskVal); + psFree(sample); + psFree(sampleMask); + background->data.F32[i][cellNum] = getStat(bgStats, options->background); + psTrace(__func__, 3, "Background for %s, cell %d is %f\n", name, cellNum, + background->data.F32[i][cellNum]); + } + + pmCellFreeData(cell); + } + pmChipFreeData(chip); + } + pmFPAFreeData(fpa); + psFitsClose(inFile); + } + psFree(bgStats); + + if (options->scale) { + // Need to normalize over the focal plane + bool converge = true; // Did we converge? + if (psTraceGetLevel(__func__) > 9) { + for (int i = 0; i < gains->n; i++) { + psTrace(__func__, 10, "Gain for cell %d is %f\n", i, gains->data.F32[i]); + } + } + 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_F32); + psImage *scalesDeref = *scales; // Dereference the pointer + + for (int i = 0; i < scalesDeref->numRows; i++) { + psF32 bg = fluxes->data.F32[i]; + for (int j = 0; j < scalesDeref->numCols; j++) { + scalesDeref->data.F32[i][j] = bg; + } + } + } + + if (options->exptime) { + if (!options->scale) { + // Copy over the exposure times + psImageCopy(*scales, exptime, PS_TYPE_F32); + } 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_F32); + } + } + + // Diagnostic stuff + if (scales && *scales && psTraceGetLevel(__func__) > 9) { + psImage *scalesDeref = *scales; // Dereference the pointer + for (int i = 0; i < scalesDeref->numRows; i++) { + for (int j = 0; j < scalesDeref->numCols; j++) { + psTrace(__func__, 9, "Scale for exposure %d, cell %d is: %f\n", i, j, + scalesDeref->data.F32[i][j]); + } + } + } + if (zeros && *zeros && psTraceGetLevel(__func__) > 9) { + psImage *zerosDeref = *zeros; // Dereference the pointer + for (int i = 0; i < zerosDeref->numRows; i++) { + for (int j = 0; j < zerosDeref->numCols; j++) { + psTrace(__func__, 9, "Zero for exposure %d, cell %d is: %f\n", i, j, + zerosDeref->data.F32[i][j]); + } + } + } + + psFree(background); + + return status; +} + Index: /trunk/ppMerge/src/ppMergeScaleZero.h =================================================================== --- /trunk/ppMerge/src/ppMergeScaleZero.h (revision 7261) +++ /trunk/ppMerge/src/ppMergeScaleZero.h (revision 7261) @@ -0,0 +1,18 @@ +#ifndef PP_MERGE_BACKGROUND_H +#define PP_MERGE_BACKGROUND_H + +#include +#include + +#include "ppMergeData.h" +#include "ppMergeOptions.h" + +// Get the scale and zero for each chip of each input +bool ppMergeScaleZero(psImage **scales, // The scales for each integration/cell + psImage **zeros, // The zeroes for each integration/cell + ppMergeData *data,// The data + const ppMergeOptions *options, // The options + const pmConfig *config // The configuration + ); + +#endif