Index: unk/ppNorm/src/ppNorm.c =================================================================== --- /trunk/ppNorm/src/ppNorm.c (revision 9507) +++ (revision ) @@ -1,33 +1,0 @@ -#ifdef HAVE_CONFIG_H -#include -#endif - -#include -#include -#include - -#include "ppNorm.h" - -int main(int argc, char *argv[]) -{ - - psLibInit(NULL); - psTimerStart(TIMERNAME); - - // Parse the configuration and arguments - pmConfig *config = pmConfigRead(&argc, argv); - - // Get the options, open the files - ppNormData *data = ppNormSetup(config); - - // Go through the FPA and do the hard work - ppNormLoop(data, config); - - psFree(data); - psFree(config); - pmConceptsDone(); - pmConfigDone(); - psLibFinalize(); - - return EXIT_SUCCESS; -} Index: unk/ppNorm/src/ppNorm.h =================================================================== --- /trunk/ppNorm/src/ppNorm.h (revision 9507) +++ (revision ) @@ -1,21 +1,0 @@ -#ifndef PP_NORM_H -#define PP_NORM_H - -#define RECIPENAME "PPNORM" -#define TIMERNAME "PPNORM" -#define INPUTNAME "PPNORM.INPUT" -#define OUTPUTNAME "PPNORM.OUTPUT" - -#include -#include "ppNormData.h" - -// Set up the options and input/output files -ppNormData *ppNormSetup(pmConfig *config // Configuration - ); - -// Loop over the input image and do all the hard work -void ppNormLoop(ppNormData *data, // The data - pmConfig *config // Configuration - ); - -#endif Index: unk/ppNorm/src/ppNormData.c =================================================================== --- /trunk/ppNorm/src/ppNormData.c (revision 9507) +++ (revision ) @@ -1,31 +1,0 @@ -#ifdef HAVE_CONFIG_H -#include -#endif - -#include -#include - -#include "ppNormData.h" - -static void normDataFree(ppNormData *data // Data to free - ) -{ - // Nothing to free; this function is just for identification purposes - return; -} - - -ppNormData *ppNormDataAlloc(void) -{ - ppNormData *data = psAlloc(sizeof(ppNormData)); // Newly allocated data - psMemSetDeallocator(data, (psFreeFunc)normDataFree); - - data->sample = 0; - data->maskVal = 0; - data->stat = 0; - data->rej = 0.0; - data->iter = 0; - - return data; -} - Index: unk/ppNorm/src/ppNormData.h =================================================================== --- /trunk/ppNorm/src/ppNormData.h (revision 9507) +++ (revision ) @@ -1,19 +1,0 @@ -#ifndef PP_NORM_DATA_H -#define PP_NORM_DATA_H - -#include -#include - -// Options for processing -typedef struct { - float sample; // Fraction of cell to sample for statistics - psMaskType maskVal; // Mask value for images - psStatsOptions stat; // Statistic to use - float rej; // Rejection limit - int iter; // Rejection iterations -} ppNormData; - -// Allocator -ppNormData *ppNormDataAlloc(void); - -#endif Index: unk/ppNorm/src/ppNormLoop.c =================================================================== --- /trunk/ppNorm/src/ppNormLoop.c (revision 9507) +++ (revision ) @@ -1,249 +1,0 @@ -#ifdef HAVE_CONFIG_H -#include -#endif - -#include -#include -#include -#include -#include - -#include "ppNorm.h" - -#define MAX_ITER 10 // Maximum number of iterations for pmFlatNormalize -#define TOLERANCE 1e-3 // Tolerance to reach for pmFlatNormalize - - -// Count the number of cells in an FPA -static int countCells(pmFPAfile *file, // The FPA file for which to count cells - pmConfig *config // Configuration - ) -{ - pmFPAview *view = pmFPAviewAlloc(0);// View of FPA - long numCells = 0; // Number of cells - - pmChip *chip; // Chip, from iteration - while ((chip = pmFPAviewNextChip(view, file->fpa, 1))) { - if (!chip->file_exists) { - continue; - } - - pmCell *cell; // Cell, from iteration - while ((cell = pmFPAviewNextCell(view, file->fpa, 1))) { - if (!cell->file_exists) { - continue; - } - - numCells++; - } - } - - psFree(view); - return numCells; -} - - - - - -void ppNormLoop(ppNormData *data, // The data - pmConfig *config // Configuration - ) -{ - assert(psStatsSingleOption(data->stat)); - psStats *bgStats = psStatsAlloc(data->stat); // Statistics for the background - bgStats->clipIter = data->iter; - bgStats->clipSigma = data->rej; - - bool mdok; // Status of MD lookup - pmFPAfile *input = psMetadataLookupPtr(&mdok, config->files, INPUTNAME); // Input FPA file - if (!input || !mdok) { - psAbort(__func__, "Should never get here: I'm sure I saved it.\n"); - } - - long numCells = countCells(input, config); // Number cells - - psImage *backgrounds = psImageAlloc(numCells, 1, PS_TYPE_F32); // Matrix of backgrounds. Only single exp. - psVector *gains = psVectorAlloc(numCells, PS_TYPE_F32); // Vector of gains for the cells - gains->n = numCells; - - // Iterate through the FPA, getting the background - pmFPAview *view = pmFPAviewAlloc(0);// View of FPA - long cellNum = -1; // Cell number, for indexing - if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) { - psError(PS_ERR_IO, false, "Unable to check file I/O.\n"); - goto done; - } - pmChip *chip; // Chip, from iteration - while ((chip = pmFPAviewNextChip(view, input->fpa, 1))) { - if (!chip->file_exists) { - continue; - } - const char *chipName = psMetadataLookupStr(NULL, chip->concepts, "CHIP.NAME"); // Name of chip - if (!pmFPAfileIOChecks (config, view, PM_FPA_BEFORE)) { - psError(PS_ERR_IO, false, "Unable to check file I/O at chip %s.\n", chipName); - goto done; - } - pmCell *cell; // Cell, from iteration - while ((cell = pmFPAviewNextCell(view, input->fpa, 1))) { - if (!cell->file_exists) { - continue; - } - cellNum++; - const char *cellName = psMetadataLookupStr(NULL, cell->concepts, "CELL.NAME"); // Name of cell - if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) { - psError(PS_ERR_IO, false, "Unable to check file I/O at chip %s, cell %s.\n", - chipName, cellName); - goto done; - } - - if (cell->readouts->n > 1) { - psLogMsg(__func__, PS_LOG_WARN, "Chip %s, cell %s contains %d readouts --- " - "only the first will be used.\n", chipName, cellName, cell->readouts->n); - } - pmReadout *readout = cell->readouts->data[0]; - - if (data->sample <= 0.0) { - psImageStats(bgStats, readout->image, readout->mask, data->maskVal); - } else { - // Apply sampling - psImage *image = readout->image; // The image of interest - psImage *mask = readout->mask; // The mask image - int numSamples = data->sample * image->numCols * image->numRows; // Number of samples - int sampleSpace = 1.0 / data->sample; // Space between samples - psVector *sampleValues = psVectorAlloc(numSamples, PS_TYPE_F32); // Vector of samples - sampleValues->n = numSamples; - psVector *sampleMask = NULL; // Corresponding mask - if (mask) { - sampleMask = psVectorAlloc(numSamples, PS_TYPE_U8); - sampleMask->n = numSamples; - } - for (int i = 0; i < numSamples; i++) { - int j = i * sampleSpace; - int y = j / image->numRows; - int x = j % image->numRows; - sampleValues->data.F32[i] = image->data.F32[y][x]; - if (mask) { - sampleMask->data.U8[i] = mask->data.U8[y][x]; - } - } - psVectorStats(bgStats, sampleValues, sampleMask, NULL, data->maskVal); - psFree(sampleValues); - psFree(sampleMask); - } - - backgrounds->data.F32[0][cellNum] = psStatsGetValue(bgStats, data->stat); - psTrace(__func__, 3, "Background for chip %s, cell %s (%d) is %f\n", chipName, cellName, - cellNum, backgrounds->data.F32[0][cellNum]); - - gains->data.F32[cellNum] = psMetadataLookupF32(&mdok, cell->concepts, "CELL.GAIN"); - if (!mdok || isnan(gains->data.F32[cellNum])) { - psLogMsg(__func__, PS_LOG_WARN, "Unable to access CELL.GAIN for cell --- ignored.\n"); - } - - if (!pmFPAfileIOChecks (config, view, PM_FPA_AFTER)) { - psError(PS_ERR_IO, false, "Unable to check file I/O at chip %s, cell %s.\n", - chipName, cellName); - goto done; - } - pmCellFreeData(cell); - } - if (!pmFPAfileIOChecks (config, view, PM_FPA_AFTER)) { - psError(PS_ERR_IO, false, "Unable to check file I/O at chip %s.\n", chipName); - goto done; - } - pmChipFreeData(chip); - } - if (!pmFPAfileIOChecks (config, view, PM_FPA_AFTER)) { - psError(PS_ERR_IO, false, "Unable to check file I/O.\n"); - goto done; - } - pmFPAFreeData(input->fpa); - - bool converge = false; // Did the normalisation converge? - psVector *fluxes = pmFlatNormalize(&converge, gains, backgrounds, MAX_ITER, TOLERANCE); - if (!converge || !fluxes) { - psError(PS_ERR_UNKNOWN, false, "Normalisation failed to converge.\n"); - goto done; - } - - assert(fluxes->n == 1); // Only a single integration - if (fluxes->data.F32[0] <= 0.0) { - psError(PS_ERR_UNKNOWN, false, "Non-positive value (%f) returned for normalisation.\n", - fluxes->data.F32[0]); - goto done; - } - float invFlux = 1.0 / fluxes->data.F32[0]; // The inverse flux; multiply by this to normalise - - // Go through the image, normalising the fluxes - pmFPAfileActivate(config->files, true, OUTPUTNAME); - pmFPAviewReset(view); - - { - long cellNum = -1; // Cell number, for indexing - if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) { - psError(PS_ERR_IO, false, "Unable to check file I/O.\n"); - goto done; - } - pmChip *chip; // Chip, from iteration - while ((chip = pmFPAviewNextChip(view, input->fpa, 1))) { - if (!chip->file_exists) { - continue; - } - const char *chipName = psMetadataLookupStr(NULL, chip->concepts, "CHIP.NAME"); // Name of chip - if (!pmFPAfileIOChecks (config, view, PM_FPA_BEFORE)) { - psError(PS_ERR_IO, false, "Unable to check file I/O at chip %s.\n", chipName); - goto done; - } - pmCell *cell; // Cell, from iteration - while ((cell = pmFPAviewNextCell(view, input->fpa, 1))) { - if (!cell->file_exists) { - continue; - } - cellNum++; - const char *cellName = psMetadataLookupStr(NULL, cell->concepts, "CELL.NAME"); // Name of cell - if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) { - psError(PS_ERR_IO, false, "Unable to check file I/O at chip %s, cell %s.\n", - chipName, cellName); - goto done; - } - - if (cell->readouts->n > 1) { - psLogMsg(__func__, PS_LOG_WARN, "Chip %s, cell %s contains %d readouts --- " - "only the first will be used.\n", chipName, cellName, cell->readouts->n); - } - pmReadout *readout = cell->readouts->data[0]; - - psTrace(__func__, 3, "Dividing chip %s, cell %s (%d) by %f\n", chipName, cellName, - cellNum, fluxes->data.F32[0]); - - psBinaryOp(readout->image, readout->image, "*", psScalarAlloc(invFlux, PS_TYPE_F32)); - - if (!pmFPAfileIOChecks (config, view, PM_FPA_AFTER)) { - psError(PS_ERR_IO, false, "Unable to check file I/O at chip %s, cell %s.\n", - chipName, cellName); - goto done; - } - pmCellFreeData(cell); - } - if (!pmFPAfileIOChecks (config, view, PM_FPA_AFTER)) { - psError(PS_ERR_IO, false, "Unable to check file I/O at chip %s.\n", chipName); - goto done; - } - pmChipFreeData(chip); - } - if (!pmFPAfileIOChecks (config, view, PM_FPA_AFTER)) { - psError(PS_ERR_IO, false, "Unable to check file I/O.\n"); - goto done; - } - pmFPAFreeData(input->fpa); - } - -done: - psFree(view); - psFree(bgStats); - psFree(backgrounds); - psFree(gains); - - return; -} Index: unk/ppNorm/src/ppNormSetup.c =================================================================== --- /trunk/ppNorm/src/ppNormSetup.c (revision 9507) +++ (revision ) @@ -1,165 +1,0 @@ -#ifdef HAVE_CONFIG_H -#include -#endif - -#ifdef HAVE_CONFIG_H -#include -#endif - -#include -#include -#include -#include - -#include "ppNorm.h" -#include "ppNormData.h" - -// Print usage information and die -static void usageAndDie(pmConfig *config // Configuration (contains the arguments list) - ) -{ - printf("Normalize the input FPA.\n\n" - "Usage:\n" - "\t%s [-file INPUT.fits] [-list LIST] OUTPUT\n" - "\n", config->argv[0]); - psArgumentHelp(config->arguments); - psFree(config); - psLibFinalize(); - pmConceptsDone(); - pmConfigDone(); - exit(EXIT_FAILURE); -} - -ppNormData *ppNormSetup(pmConfig *config // Configuration - ) -{ - // Setup and parse command-line arguments - psMetadata *arguments = config->arguments; // Arguments - psMetadataAddStr(arguments, PS_LIST_TAIL, "-stat", 0, "Statistic to use", NULL); - psMetadataAddF32(arguments, PS_LIST_TAIL, "-sample", 0, "Sampling fraction", 0.0); - psMetadataAddF32(arguments, PS_LIST_TAIL, "-rej", 0, "Rejection level", 0.0); - psMetadataAddS32(arguments, PS_LIST_TAIL, "-iter", 0, "Rejection iterations", 0); - - // No command-line arguments: print the help - if (*config->argc == 1) { - usageAndDie(config); - } - - // Get input files - if (!pmConfigFileSetsMD(config->arguments, config, "INPUT", "-file", "-list")) { - psError(PS_ERR_IO, false, "Unable to open input files.\n"); - usageAndDie(config); - } - - if (!psArgumentParse(arguments, config->argc, config->argv)) { - printf("Unable to parse command-line arguments.\n\n"); - usageAndDie(config); - } - - // The output (root) name is the final thing left on the command-line - if (*config->argc != 2) { - usageAndDie(config); - } - // Parse the command-line options - ppNormData *data = ppNormDataAlloc(); // The data - psMetadataAddStr(config->arguments, PS_LIST_TAIL, "OUTPUT", 0, - "Name of the output image", config->argv[1]); - - // Read the camera - bool status = false; // Status of function calls - pmFPAfile *input = pmFPAfileDefineFromArgs(&status, config, INPUTNAME, "INPUT"); - if (!status) { - psError(PS_ERR_IO, false, "Failed to build FPA from %s.INPUT.\n", INPUTNAME); - goto die; - } - - pmFPAfile *output = pmFPAfileDefineOutput(config, input->fpa, OUTPUTNAME); - if (!output) { - psError(PS_ERR_UNKNOWN, false, "Failed to build output from %s.OUTPUT.\n", OUTPUTNAME); - goto die; - } - pmFPAfileActivate(config->files, false, OUTPUTNAME); // We don't want to use this until later - -#if 0 - psFree(input); // Drop reference - psFree(output); // Drop reference -#endif - - // Get the command-line arguments - const char *statString = psMetadataLookupStr(NULL, arguments, "-stat"); // Statistic of choice - if (statString && strlen(statString) > 0) { - data->stat = psStatsOptionFromString(statString); - if (data->stat == 0) { - psError(PS_ERR_BAD_PARAMETER_VALUE, false, "Unable to interpret statistic: %s\n", statString); - goto die; - } - } - data->sample = psMetadataLookupF32(NULL, arguments, "-sample"); - - // Determine recipe parameters - bool mdok; // Status of MD lookup - psMetadata *recipe = psMetadataLookupMD(&mdok, config->recipes, RECIPENAME); - if (!mdok || !recipe) { - psLogMsg(__func__, PS_LOG_WARN, "Unable to find recipe %s --- attempting to proceed.\n", RECIPENAME); - return data; - } - - if (data->sample == 0) { - float sample = psMetadataLookupF32(&mdok, recipe, "SAMPLE"); // Sample fraction - if (mdok && sample > 0) { - data->sample = sample; - } else { - psLogMsg(__func__, PS_LOG_WARN, "SAMPLE in recipe is not of type F32 and positive --- " - "retaining default.\n"); - } - } - - psMaskType maskVal = psMetadataLookupU8(&mdok, recipe, "MASKVAL"); // Mask value - if (mdok) { - data->maskVal = maskVal; - } else { - psLogMsg(__func__, PS_LOG_WARN, "ITER in recipe is not of type U8 --- retaining default.\n"); - } - - if (data->stat == 0) { - const char *statName = psMetadataLookupStr(&mdok, recipe, "STAT"); // Statistics name - if (!mdok || !statName || strlen(statName) == 0) { - psError(PS_ERR_UNKNOWN, true, "Unable to find STAT of type STR in the recipe.\n"); - goto die; - } - data->stat = psStatsOptionFromString(statName); - if (data->stat == 0) { - psError(PS_ERR_BAD_PARAMETER_VALUE, false, "Unable to parse STAT in recipe: %s\n", statName); - goto die; - } - } - - // Clipping options - if (data->stat == PS_STAT_CLIPPED_MEAN && data->iter == 0 && isnan(data->rej)) { - int iter = psMetadataLookupS32(&mdok, recipe, "ITER"); // Number of rejection iterations - if (mdok && iter > 0) { - data->iter = iter; - } else { - psError(PS_ERR_BAD_PARAMETER_VALUE, true, "ITER in recipe is not of type S32 and positive.\n"); - goto die; - } - float rej = psMetadataLookupF32(&mdok, recipe, "REJ"); // Rejection level - if (mdok && rej > 0) { - data->rej = rej; - } else { - psError(PS_ERR_BAD_PARAMETER_VALUE, true, "REJ in recipe is not of type F32 and positive.\n"); - goto die; - } - } - - return data; - - // Common path for error conditions: clean up and exit. -die: - psFree(config); - psFree(data); - pmConceptsDone(); - pmConfigDone(); - psLibFinalize(); - exit(EXIT_FAILURE); -}