Index: /trunk/ppMerge/src/Makefile.am =================================================================== --- /trunk/ppMerge/src/Makefile.am (revision 17226) +++ /trunk/ppMerge/src/Makefile.am (revision 17227) @@ -6,29 +6,14 @@ ppMerge_SOURCES = \ ppMerge.c \ - ppMergeCheckInputs.c \ - ppMergeCombine.c \ - ppMergeConfig.c \ - ppMergeData.c \ - ppMergeMask.c \ - ppMergeMaskSuspect.c \ - ppMergeMaskBad.c \ - ppMergeMaskWrite.c \ - ppMergeMaskGrow.c \ - ppMergeMaskAverageConcepts.c \ - ppMergeOptions.c \ + ppMergeArguments.c \ + ppMergeCamera.c \ + ppMergeFiles.c \ ppMergeScaleZero.c \ - ppMergeVersion.c + ppMergeLoop.c \ + ppMergeMask.c noinst_HEADERS = \ - ppMerge.h \ - ppMergeCheckInputs.h \ - ppMergeCombine.h \ - ppMergeConfig.h \ - ppMergeData.h \ - ppMergeMask.h \ - ppMergeOptions.h \ - ppMergeScaleZero.h \ - ppMergeVersion.h + ppMerge.h CLEANFILES = *~ Index: /trunk/ppMerge/src/ppMerge.c =================================================================== --- /trunk/ppMerge/src/ppMerge.c (revision 17226) +++ /trunk/ppMerge/src/ppMerge.c (revision 17227) @@ -4,15 +4,9 @@ #include +#include #include #include #include "ppMerge.h" -#include "ppMergeConfig.h" -#include "ppMergeData.h" -#include "ppMergeOptions.h" -#include "ppMergeCheckInputs.h" -#include "ppMergeCombine.h" -#include "ppMergeScaleZero.h" -#include "ppMergeMask.h" //#include "ppMem.h" @@ -27,17 +21,81 @@ { psLibInit(NULL); - psMemSetThreadSafety(false); // Turn off thread safety, for more + psMemSetThreadSafety(false); psTimerStart(TIMERNAME); - // Parse the configuration and arguments - // Open the input image(s) - // Determine camera, format from header if not already defined - // Construct camera in preparation for reading - pmConfig *config = ppMergeConfig(argc, argv); + psExit exitValue = PS_EXIT_SUCCESS; // Exit value for program + + pmConfig *config = pmConfigRead(&argc, argv, PPMERGE_RECIPE); // Configuration if (!config) { - psErrorStackPrint(stderr, "Unable to get configuration."); - exit(EXIT_FAILURE); + psErrorStackPrint(stderr, "Error reading configuration."); + exitValue = PS_EXIT_CONFIG_ERROR; + goto die; } + if (!ppMergeArguments(argc, argv, config)) { + psErrorStackPrint(stderr, "Error reading arguments."); + exitValue = PS_EXIT_CONFIG_ERROR; + goto die; + } + + ppMergeType type = psMetadataLookupS32(NULL, config->arguments, "TYPE"); // Type of frame + switch (type) { + case PPMERGE_TYPE_MASK: + if (!ppMergeMask(config)) { + psErrorStackPrint(stderr, "Error generating mask."); + exitValue = PS_EXIT_DATA_ERROR; + goto die; + } + break; + case PPMERGE_TYPE_BIAS: + case PPMERGE_TYPE_DARK: + case PPMERGE_TYPE_SHUTTER: + case PPMERGE_TYPE_FLAT: + case PPMERGE_TYPE_FRINGE: + if (!ppMergeScaleZero(config)) { + psErrorStackPrint(stderr, "Error getting scale and zero-points."); + exitValue = PS_EXIT_DATA_ERROR; + goto die; + } + if (!ppMergeLoop(config)) { + psErrorStackPrint(stderr, "Error performing merge."); + exitValue = PS_EXIT_PROG_ERROR; + goto die; + } + break; + default: + psAbort("Invalid frame type: %x", type); + } + + + // Output the statistics + bool mdok; // Status of MD lookup + psString statsName = psMetadataLookupStr(&mdok, config->arguments, "STATS.NAME"); // Statistics file name + if (mdok && statsName && strlen(statsName) > 0) { + psString resolved = pmConfigConvertFilename(statsName, config, true); // Resolved filename + FILE *statsFile = fopen(resolved, "w"); // Output statistics file + if (!statsFile) { + psError(PS_ERR_IO, true, "Unable to open statistics file %s for writing.", resolved); + psFree(resolved); + exitValue = PS_EXIT_CONFIG_ERROR; + goto die; + } + psFree(resolved); + psMetadata *stats = psMetadataLookupMetadata(&mdok, config->arguments, "STATS.DATA"); // Statistics + if (!stats) { + psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to find statistics"); + exitValue = PS_EXIT_PROG_ERROR; + goto die; + } + psString statsOut = psMetadataConfigFormat(stats); // String to write out + fprintf(statsFile, "%s", statsOut); + psFree(statsOut); + fclose(statsFile); + } + + + + +#if 0 // Set various tasks (define optional operations) ppMergeOptions *options = ppMergeOptionsParse(config); @@ -81,19 +139,17 @@ } -#if 0 - pmFPAPrint(stdout, data->out, true, true); -#endif - // Cleaning up psFree(data); psFree(options); +#endif + + die: + psTrace("ppSub", 1, "Finished at %f sec\n", psTimerMark("ppSub")); + psTimerStop(); + psFree(config); - - pmConceptsDone(); pmConfigDone(); psLibFinalize(); -// ppMemCheck(); - - return EXIT_SUCCESS; + exit(exitValue); } Index: /trunk/ppMerge/src/ppMerge.h =================================================================== --- /trunk/ppMerge/src/ppMerge.h (revision 17226) +++ /trunk/ppMerge/src/ppMerge.h (revision 17227) @@ -2,6 +2,81 @@ #define PP_MERGE_H -#define TIMERNAME "ppMerge" -#define PPMERGE_RECIPE "PPMERGE" +#define TIMERNAME "ppMerge" // Name for timer +#define PPMERGE_RECIPE "PPMERGE" // Recipe name + +// Type of frame to merge +typedef enum { + PPMERGE_TYPE_UNKNOWN, // Unknown type + PPMERGE_TYPE_BIAS, // Bias frame + PPMERGE_TYPE_DARK, // (Multi-)Dark frame + PPMERGE_TYPE_MASK, // Mask frame + PPMERGE_TYPE_SHUTTER, // Shutter frame + PPMERGE_TYPE_FLAT, // Flat-field frame (dome or sky) + PPMERGE_TYPE_FRINGE, // Fringe frame +} ppMergeType; + +// Files, for activation +typedef enum { + PPMERGE_FILES_ALL, // All files + PPMERGE_FILES_INPUT, // Input files + PPMERGE_FILES_OUTPUT // Output files +} ppMergeFiles; + +// Parse command-line arguments and recipe +bool ppMergeArguments(int argc, char *argv[], // Command-line arguments + pmConfig *config // Configuration + ); + +// Set up camera files +bool ppMergeCamera(pmConfig *config // Configuration + ); + +// Measure scale and zero-points +bool ppMergeScaleZero(pmConfig *config // Configuration + ); + +// Main loop to do the merging +bool ppMergeLoop(pmConfig *config // Configuration + ); + +// Main loop for masks +bool ppMergeMask(pmConfig *config // Configuration + ); + +// Read nominated input file +bool ppMergeFileReadInput(const pmConfig *config, // Configuration + pmReadout *readout, // Readout into which to read + int num, // Number of file in sequence + int rows // Number of rows to read at once + ); + +// Open nominated input file +bool ppMergeFileOpenInput(pmConfig *config, // Configuration + const pmFPAview *view, // View to open + int num // Number of file in sequence + ); + +// Set the data level for files specified by name; return an array of the files +psArray *ppMergeFileDataLevel(const pmConfig *config, // Configuration + const char *name, // Name of files + pmFPALevel level // Level for file data level + ); + +// Activate/deactivate a list of files +bool ppMergeFileActivate(const pmConfig *config, // Configuration + ppMergeFiles files, // Files to turn on/off + bool state // Activation state + ); + +// Activate/deactivate a single element for a list; return array of files +psArray *ppMergeFileActivateSingle(const pmConfig *config, // Configuration + ppMergeFiles files, // Files to turn on/off + bool state, // Activation state + int num // Number of file in sequence + ); + +// Return name of output pmFPAfile +psString ppMergeOutputFile(const pmConfig *config // Configuration + ); #endif Index: /trunk/ppMerge/src/ppMergeArguments.c =================================================================== --- /trunk/ppMerge/src/ppMergeArguments.c (revision 17227) +++ /trunk/ppMerge/src/ppMergeArguments.c (revision 17227) @@ -0,0 +1,343 @@ +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include +#include +#include +#include + +#include "ppMerge.h" + +// Print usage information and die +static void usage(const char *program, // Name of the program + psMetadata *arguments // Command-line arguments + ) +{ + fprintf(stderr, "\nPan-STARRS Detrend Merging\n\n"); + fprintf(stderr, "Usage: %s INPUT.mdc OUTPUT_ROOT\n" + "where INPUTS.mdc contains various METADATAs, each with:\n" + "\tIMAGE(STR): Image filename\n" + "\tMASK(STR): Mask filename\n" + "\tWEIGHT(STR) Weight map filename\n" + "where MASK and WEIGHT are optional.", + program); + fprintf(stderr, "\n"); + psArgumentHelp(arguments); +} + +// Get a float-point value from the command-line or recipe, and add it to the arguments +#define VALUE_ARG_RECIPE_FLOAT(ARGNAME, RECIPENAME, TYPE) { \ + ps##TYPE value = psMetadataLookup##TYPE(NULL, arguments, ARGNAME); \ + if (isnan(value)) { \ + bool mdok; \ + value = psMetadataLookup##TYPE(&mdok, recipe, RECIPENAME); \ + if (!mdok) { \ + psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Unable to find %s in recipe %s", \ + RECIPENAME, PPMERGE_RECIPE); \ + goto ERROR; \ + } \ + } \ + psMetadataAdd##TYPE(config->arguments, PS_LIST_TAIL, RECIPENAME, 0, NULL, value); \ +} + +// Get an integer value from the command-line or recipe, and add it to the arguments +#define VALUE_ARG_RECIPE_INT(ARGNAME, RECIPENAME, TYPE, UNSET) { \ + ps##TYPE value = psMetadataLookup##TYPE(NULL, arguments, ARGNAME); \ + if (value == UNSET) { \ + bool mdok; \ + value = psMetadataLookup##TYPE(&mdok, recipe, RECIPENAME); \ + if (!mdok) { \ + psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Unable to find %s in recipe %s", \ + RECIPENAME, PPMERGE_RECIPE); \ + goto ERROR; \ + } \ + } \ + psMetadataAdd##TYPE(config->arguments, PS_LIST_TAIL, RECIPENAME, 0, NULL, value); \ +} + +// Get a boolean from the command-line or recipe, and add it to the arguments if either is set +#define VALUE_ARG_RECIPE_BOOL(ARGNAME, RECIPENAME) { \ + bool value = (psMetadataLookupBool(NULL, arguments, ARGNAME) || \ + psMetadataLookupBool(NULL, recipe, RECIPENAME)); \ + psMetadataAddBool(config->arguments, PS_LIST_TAIL, RECIPENAME, 0, NULL, value); \ +} + +// Get a statistic name from the command-line or recipe, and add the enum to the arguments +#define VALUE_ARG_RECIPE_STAT(ARGNAME, RECIPENAME) { \ + const char *stat = psMetadataLookupStr(NULL, arguments, ARGNAME); \ + if (!stat) { \ + stat = psMetadataLookupStr(NULL, recipe, RECIPENAME); \ + if (!stat) { \ + psError(PS_ERR_BAD_PARAMETER_VALUE, false, "Unable to find %s in recipe %s", \ + RECIPENAME, PPMERGE_RECIPE); \ + goto ERROR; \ + } \ + } \ + psMetadataAddS32(config->arguments, PS_LIST_TAIL, RECIPENAME, 0, NULL, psStatsOptionFromString(stat)); \ +} + +// Get a string from the command-line or recipe, and add to the arguments +#define VALUE_ARG_RECIPE_STR(ARGNAME, RECIPENAME) { \ + const char *str = psMetadataLookupStr(NULL, arguments, ARGNAME); \ + if (!str) { \ + str = psMetadataLookupStr(NULL, recipe, RECIPENAME); \ + if (!str) { \ + psError(PS_ERR_BAD_PARAMETER_VALUE, false, "Unable to find %s in recipe %s", \ + RECIPENAME, PPMERGE_RECIPE); \ + goto ERROR; \ + } \ + } \ + psMetadataAddStr(config->arguments, PS_LIST_TAIL, RECIPENAME, 0, NULL, str); \ +} + +// Get a string from the command-line or recipe, and add the appropriate mask value to the arguments +#define VALUE_ARG_RECIPE_MASK(ARGNAME, RECIPENAME) { \ + const char *str = psMetadataLookupStr(NULL, arguments, ARGNAME); \ + if (!str) { \ + str = psMetadataLookupStr(NULL, recipe, RECIPENAME); \ + if (!str) { \ + psError(PS_ERR_BAD_PARAMETER_VALUE, false, "Unable to find %s in recipe %s", \ + RECIPENAME, PPMERGE_RECIPE); \ + goto ERROR; \ + } \ + } \ + psMaskType mask = pmConfigMask(str, config); \ + psMetadataAddU8(config->arguments, PS_LIST_TAIL, RECIPENAME, 0, NULL, mask); \ +} + +// Get a string value from the command-line and add it to the target +static bool valueArgStr(psMetadata *arguments, // Command-line arguments + const char *argName, // Argument name in the command-line arguments + const char *mdName, // Name for value in the metadata + psMetadata *target // Target metadata to which to add value + ) +{ + psString value = psMetadataLookupStr(NULL, arguments, argName); // Value of interest + if (value && strlen(value) > 0) { + return psMetadataAddStr(target, PS_LIST_TAIL, mdName, 0, NULL, value); + } + return false; +} + +bool ppMergeArguments(int argc, char *argv[], pmConfig *config) +{ + assert(config); + + psMetadata *arguments = psMetadataAlloc(); // Command-line arguments + psMetadataAddStr(arguments, PS_LIST_TAIL, "-type", 0, "Type of calibration frame", NULL); + psMetadataAddStr(arguments, PS_LIST_TAIL, "-stats", 0, "MDC file to hold statistics ", NULL); + // Standard combination parameters + psMetadataAddS32(arguments, PS_LIST_TAIL, "-rows", 0, "Rows to read per scan", 0); + psMetadataAddS32(arguments, PS_LIST_TAIL, "-sample", 0, "Sampling factor for background", 0); + psMetadataAddS32(arguments, PS_LIST_TAIL, "-iter", 0, "Number of rejection iterations", 0); + psMetadataAddF32(arguments, PS_LIST_TAIL, "-rej", 0, "Rejection threshold (sigma)", NAN); + psMetadataAddF32(arguments, PS_LIST_TAIL, "-fraclow", 0, "Fraction of low pixels to discard", NAN); + psMetadataAddF32(arguments, PS_LIST_TAIL, "-frachigh", 0, "Fraction of low pixels to discard", NAN); + psMetadataAddS32(arguments, PS_LIST_TAIL, "-nkeep", 0, "Minimum number of pixels in stack to keep", 0); + psMetadataAddBool(arguments, PS_LIST_TAIL, "-weights", 0, "Use image weights in combination?", false); + psMetadataAddStr(arguments, PS_LIST_TAIL, "-maskval", 0, "Mask value for input data", NULL); + psMetadataAddStr(arguments, PS_LIST_TAIL, "-combine", 0, "Statistic to use for combination", NULL); + psMetadataAddStr(arguments, PS_LIST_TAIL, "-mean", 0, "Statistic to use to measure the mean", NULL); + psMetadataAddStr(arguments, PS_LIST_TAIL, "-stdev", 0, "Statistic to use to measure the stdev", NULL); + + // Fringe construction parameters + psMetadataAddS32(arguments, PS_LIST_TAIL, "-fringe-num", 0, "Number of fringe regions", 0); + psMetadataAddS32(arguments, PS_LIST_TAIL, "-fringe-size", 0, "Half-size of fringe regions", 0); + psMetadataAddS32(arguments, PS_LIST_TAIL, "-fringe-xsmooth", 0, "Number of smoothing regions in x", 0); + psMetadataAddS32(arguments, PS_LIST_TAIL, "-fringe-ysmooth", 0, "Number of smoothing regions in y", 0); + + // Shutter construction parameters + psMetadataAddS32(arguments, PS_LIST_TAIL, "-shutter-size", 0, "Size for shutter measurement regions", 0); + psMetadataAddS32(arguments, PS_LIST_TAIL, "-shutter-iter", 0, "Number of iterations for shutter", 0); + psMetadataAddF32(arguments, PS_LIST_TAIL, "-shutter-rej", 0, "Rejection limit for shutter", NAN); + + // Mask construction parameters + psMetadataAddF32(arguments, PS_LIST_TAIL, "-mask-suspect", 0, "Threshold for suspect pixels (sigma)", NAN); + psMetadataAddF32(arguments, PS_LIST_TAIL, "-mask-bad", 0, "Threshold for bad pixels (sigma)", NAN); + psMetadataAddStr(arguments, PS_LIST_TAIL, "-mask-mode", 0, "Mode to identify bad pixels", NULL); + psMetadataAddS32(arguments, PS_LIST_TAIL, "-mask-grow", 0, "Number of pixels to grow final mask", 0); + psMetadataAddStr(arguments, PS_LIST_TAIL, "-mask-growval", 0, "Value to give grown mask pixels", NULL); + psMetadataAddBool(arguments, PS_LIST_TAIL, "-mask-chip", 0, "Measure mask statistics by chip?", false); + + if (argc == 1 || !psArgumentParse(arguments, &argc, argv) || argc != 3) { + usage(argv[0], arguments); + goto ERROR; + } + + unsigned int numBad = 0; // Number of bad lines + psMetadata *inputs = psMetadataConfigRead(NULL, &numBad, argv[1], false); // Information about inputs + if (!inputs || numBad > 0) { + psError(PS_ERR_BAD_PARAMETER_VALUE, false, "Unable to cleanly read MDC file with inputs."); + goto ERROR; + } + psMetadataAddMetadata(config->arguments, PS_LIST_TAIL, "INPUTS", 0, + "Metadata with input details", inputs); + psFree(inputs); + psMetadataAddStr(config->arguments, PS_LIST_TAIL, "OUTPUT", 0, + "Root name of the output image list", argv[2]); + + valueArgStr(arguments, "-stats", "STATS.NAME", config->arguments); + + + // Set the type of calibration frame + const char *typeStr = psMetadataLookupStr(NULL, arguments, "-type"); // Type of calibration + if (strlen(typeStr) <= 0) { + psError(PS_ERR_UNKNOWN, false, "No -type specified."); + goto ERROR; + } + ppMergeType type = PPMERGE_TYPE_UNKNOWN; // Enumerated type for frame type + if (strcasecmp(typeStr, "BIAS") == 0) { + type = PPMERGE_TYPE_BIAS; + } else if (strcasecmp(typeStr, "DARK") == 0) { + type = PPMERGE_TYPE_DARK; + } else if (strcasecmp(typeStr, "FLAT") == 0 || strcasecmp(typeStr, "SKYFLAT") == 0 || + strcasecmp(typeStr, "DOMEFLAT") == 0) { + type = PPMERGE_TYPE_FLAT; + } else if (strcasecmp(typeStr, "FRINGE") == 0) { + type = PPMERGE_TYPE_FRINGE; + } else if (strcasecmp(typeStr, "SHUTTER") == 0) { + type = PPMERGE_TYPE_SHUTTER; + } else if (strcasecmp(typeStr, "MASK") == 0 || + strcasecmp(typeStr, "DARKMASK") == 0 || + strcasecmp(typeStr, "FLATMASK") == 0) { + type = PPMERGE_TYPE_MASK; + } else { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Unrecognised image type: %s", typeStr); + goto ERROR; + } + psMetadataAddS32(config->arguments, PS_LIST_TAIL, "TYPE", 0, "Type of calibration frame", type); + + // Need to set the camera before the recipes can be read + if (!ppMergeCamera(config)) { + psError(PS_ERR_UNKNOWN, false, "Unable to setup cameras."); + goto ERROR; + } + + psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, PPMERGE_RECIPE); // Recipe for ppSim + if (!recipe) { + psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to find recipe %s", PPMERGE_RECIPE); + goto ERROR; + } + + // Standard combination parameters + VALUE_ARG_RECIPE_INT("-rows", "ROWS", S32, 0); + VALUE_ARG_RECIPE_INT("-sample", "SAMPLE", S32, 0); + VALUE_ARG_RECIPE_INT("-iter", "ITER", S32, 0); + VALUE_ARG_RECIPE_FLOAT("-rej", "REJ", F32); + VALUE_ARG_RECIPE_FLOAT("-fraclow", "FRACLOW", F32); + VALUE_ARG_RECIPE_FLOAT("-frachigh", "FRACHIGH", F32); + VALUE_ARG_RECIPE_INT("-nkeep", "NKEEP", S32, 0); + VALUE_ARG_RECIPE_BOOL("-weights", "WEIGHTS"); + VALUE_ARG_RECIPE_MASK("-maskval", "MASKVAL"); + VALUE_ARG_RECIPE_STAT("-combine", "COMBINE"); + VALUE_ARG_RECIPE_STAT("-mean", "MEAN"); + VALUE_ARG_RECIPE_STAT("-stdev", "STDEV"); + + // Fringe construction parameters + VALUE_ARG_RECIPE_INT("-fringe-num", "FRINGE.NUM", S32, 0); + VALUE_ARG_RECIPE_INT("-fringe-size", "FRINGE.SIZE", S32, 0); + VALUE_ARG_RECIPE_INT("-fringe-xsmooth", "FRINGE.XSMOOTH", S32, 0); + VALUE_ARG_RECIPE_INT("-fringe-ysmooth", "FRINGE.YSMOOTH", S32, 0); + + // Shutter construction parameters + VALUE_ARG_RECIPE_INT("-shutter-size", "SHUTTER.SIZE", S32, 0); + + // Mask construction parameters + VALUE_ARG_RECIPE_FLOAT("-mask-suspect", "MASK.SUSPECT", F32); + VALUE_ARG_RECIPE_FLOAT("-mask-bad", "MASK.BAD", F32); + VALUE_ARG_RECIPE_INT("-mask-grow", "MASK.GROW", S32, 0); + VALUE_ARG_RECIPE_MASK("-mask-growval", "MASK.GROWVAL"); + VALUE_ARG_RECIPE_BOOL("-mask-chip", "MASK.CHIPSTATS"); + + + const char *maskModeStr = psMetadataLookupStr(NULL, arguments, "-mask-mode"); // Mode to identify bad pix + if (!maskModeStr) { + maskModeStr = psMetadataLookupStr(NULL, recipe, "MASK.MODE"); + if (!maskModeStr) { + psError(PS_ERR_BAD_PARAMETER_VALUE, false, "No mask mode specified in recipe."); + goto ERROR; + } + } + pmMaskIdentifyMode maskMode = pmMaskIdentifyModeFromString(maskModeStr); + if (maskMode == PM_MASK_ID_NONE) { + psError(PS_ERR_BAD_PARAMETER_VALUE, false, "Invalid mask mode %s", maskModeStr); + goto ERROR; + } + psMetadataAddS32(config->arguments, PS_LIST_TAIL, "MASK.MODE", 0, "Mode for mask identification", + maskMode); + + + if (type == PPMERGE_TYPE_DARK) { + psMetadata *ordinates = psMetadataLookupMetadata(NULL, recipe, "DARK.ORDINATES"); // Ordinates info + psArray *translated = psArrayAllocEmpty(psListLength(ordinates->list)); // Translated version + + 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); + return false; + } + if (order <= 0) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, "ORDER not positive (%d) for DARK.ORDINATES %s", + order, item->name); + return false; + } + + pmDarkOrdinate *ord = pmDarkOrdinateAlloc(item->name, order); + ord->scale = scale; + ord->min = min; + ord->max = max; + psArrayAdd(translated, translated->n, ord); + psFree(ord); + } + psFree(iter); + + psMetadataAddArray(config->arguments, PS_LIST_TAIL, "DARK.ORDINATES", 0, + "Ordinates to fit for dark", translated); + psFree(translated); // Drop reference + + psString darkNorm = psMetadataLookupStr(NULL, recipe, "DARK.NORM"); // Normalisation concept + if (darkNorm && strcmp(darkNorm, "NONE") != 0) { + psMetadataAddStr(config->arguments, PS_LIST_TAIL, "DARK.NORM", 0, + "Normalisation concept for dark", darkNorm); + } + } + + +#if 0 + // Add concepts for scale and zero + // XXX These have never been used + 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); +#endif + + return true; + +ERROR: + psFree(arguments); + return false; +} + Index: /trunk/ppMerge/src/ppMergeCamera.c =================================================================== --- /trunk/ppMerge/src/ppMergeCamera.c (revision 17227) +++ /trunk/ppMerge/src/ppMergeCamera.c (revision 17227) @@ -0,0 +1,322 @@ +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include +#include +#include +#include + +#include "ppMerge.h" + +// Define an output file, with its own FPA +bool outputFile(pmConfig *config, // Configuration + const char *name, // Name of output file + pmFPAfileType type, // Type of file + const char *description, // Description of file + psMetadata *format, // Camera format + pmFPAview *view // View for PHU + ) +{ + assert(config); + assert(name && strlen(name) > 0); + assert(view); + + // Output image + pmFPA *fpa = pmFPAConstruct(config->camera); // FPA to contain the output + if (!fpa) { + psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to construct an FPA from camera configuration."); + return false; + } + + pmFPAfile *output = pmFPAfileDefineOutput(config, fpa, name); + psFree(fpa); // Drop reference + if (!output) { + psError(PS_ERR_IO, false, "Unable to generate output file from %s", name); + return false; + } + if (output->type != type) { + psError(PS_ERR_IO, true, "%s is not of type %s", name, pmFPAfileStringFromType(type)); + return false; + } + output->save = true; + + if (!pmFPAAddSourceFromView(fpa, description, view, format)) { + psError(PS_ERR_UNKNOWN, false, "Unable to generate output FPA."); + return false; + } + + return true; +} + + +bool ppMergeCamera(pmConfig *config) +{ + bool haveMasks = false; // Do we have masks? + bool haveWeights = false; // Do we have weight maps? + + ppMergeType type = psMetadataLookupS32(NULL, config->arguments, "TYPE"); // Type of frame + psMetadata *inputs = psMetadataLookupMetadata(NULL, config->arguments, "INPUTS"); // The inputs info + psMetadataIterator *iter = psMetadataIteratorAlloc(inputs, PS_LIST_HEAD, NULL); // Iterator + psMetadataItem *item; // Item from iteration + int numFiles = 0; // Number of files + while ((item = psMetadataGetAndIncrement(iter))) { + if (item->type != PS_DATA_METADATA) { + psError(PS_ERR_BAD_PARAMETER_TYPE, true, + "Component %s of the input metadata is not of type METADATA", item->name); + psFree(iter); + return false; + } + + psMetadata *input = item->data.md; // The input metadata of interest + + psString image = psMetadataLookupStr(NULL, input, "IMAGE"); // Name of image + if (!image || strlen(image) == 0) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Component %s lacks IMAGE of type STR", item->name); + psFree(iter); + return false; + } + + bool mdok; + psString mask = psMetadataLookupStr(&mdok, input, "MASK"); // Name of mask + psString weight = psMetadataLookupStr(&mdok, input, "WEIGHT"); // Name of weight map + + // Add the image file + psArray *imageFiles = psArrayAlloc(1); // Array of filenames for this FPA + imageFiles->data[0] = psMemIncrRefCounter(image); + psMetadataAddArray(config->arguments, PS_LIST_TAIL, "IMAGE.FILENAMES", PS_META_REPLACE, + "Filenames of image files", imageFiles); + psFree(imageFiles); + + bool found = false; // Found the file? + pmFPAfile *imageFile = pmFPAfileDefineFromArgs(&found, config, "PPMERGE.INPUT", "IMAGE.FILENAMES"); + if (!imageFile || !found) { + psError(PS_ERR_UNKNOWN, false, "Unable to define file from image %d (%s)", numFiles, image); + return false; + } + if (imageFile->type != PM_FPA_FILE_IMAGE) { + psError(PS_ERR_IO, true, "PPMERGE.INPUT is not of type IMAGE"); + return false; + } + + // Optionally add the mask file + if (mask && strlen(mask) > 0) { + psArray *maskFiles = psArrayAlloc(1); // Array of filenames for this FPA + maskFiles->data[0] = psMemIncrRefCounter(mask); + psMetadataAddArray(config->arguments, PS_LIST_TAIL, "MASK.FILENAMES", PS_META_REPLACE, + "Filenames of mask files", maskFiles); + psFree(maskFiles); + + bool status; + pmFPAfile *maskFile = pmFPAfileBindFromArgs(&status, imageFile, config, "PPMERGE.INPUT.MASK", + "MASK.FILENAMES"); + if (!status) { + psError(PS_ERR_UNKNOWN, false, "Unable to define file from mask %d (%s)", numFiles, mask); + return false; + } + if (maskFile->type != PM_FPA_FILE_MASK) { + psError(PS_ERR_IO, true, "PPMERGE.INPUT.MASK is not of type MASK"); + return false; + } + haveMasks = true; + } + + // Optionally add the weight file + if (weight && strlen(weight) > 0) { + haveWeights = true; + psArray *weightFiles = psArrayAlloc(1); // Array of filenames for this FPA + weightFiles->data[0] = psMemIncrRefCounter(weight); + psMetadataAddArray(config->arguments, PS_LIST_TAIL, "WEIGHT.FILENAMES", PS_META_REPLACE, + "Filenames of weight files", weightFiles); + psFree(weightFiles); + + bool status; + pmFPAfile *weightFile = pmFPAfileBindFromArgs(&status, imageFile, config, "PPMERGE.INPUT.WEIGHT", + "WEIGHT.FILENAMES"); + if (!status) { + psError(PS_ERR_UNKNOWN, false, "Unable to define file from weight %d (%s)", numFiles, weight); + return false; + } + if (weightFile->type != PM_FPA_FILE_WEIGHT) { + psError(PS_ERR_IO, true, "PPMERGE.INPUT.WEIGHT is not of type WEIGHT"); + return false; + } + haveWeights = true; + } + + numFiles++; + } + psFree(iter); + psMetadataRemoveKey(config->arguments, "IMAGE.FILENAMES"); + if (psMetadataLookup(config->arguments, "MASK.FILENAMES")) { + psMetadataRemoveKey(config->arguments, "MASK.FILENAMES"); + } + if (psMetadataLookup(config->arguments, "WEIGHT.FILENAMES")) { + psMetadataRemoveKey(config->arguments, "WEIGHT.FILENAMES"); + } + if (psMetadataLookup(config->arguments, "PSF.FILENAMES")) { + psMetadataRemoveKey(config->arguments, "PSF.FILENAMES"); + } + + psMetadataAddS32(config->arguments, PS_LIST_TAIL, "INPUTS.NUM", 0, "Number of input files", numFiles); + psMetadataAddBool(config->arguments, PS_LIST_TAIL, "INPUTS.MASKS", 0, "Got input masks?", haveMasks); + psMetadataAddBool(config->arguments, PS_LIST_TAIL, "INPUTS.WEIGHTS", 0, + "Got input weights?", haveWeights); + + + // Check that all the inputs are consistent + +// Check an FPA level to ensure the camera format and PHU are consistent across all input files +#define CHECK_LEVEL(HDU, CHIP, CELL) { \ + if (HDU) { \ + if (!phuView) { \ + phuView = pmFPAviewAlloc(0); \ + phuView->chip = CHIP; \ + phuView->cell = CELL; \ + } else if ((phuView->chip != (CHIP)) && (phuView->cell != (CELL))) { \ + psError(PS_ERR_UNKNOWN, true, "Differing PHU for input %d", i); \ + psFree(phuView); \ + return false; \ + } \ + if (!format) { \ + format = (HDU)->format; \ + } else if (format != (HDU)->format) { \ + psError(PS_ERR_UNKNOWN, true, "Camera format %d doesn't match: %p vs %p", \ + i, format, (HDU)->format); \ + psFree(phuView); \ + return false; \ + } \ + continue; \ + } \ +} + + psMetadata *format = NULL; // Camera format + pmFPAview *phuView = NULL; // View to PHU + for (int i = 0; i < numFiles; i++) { + pmFPAfile *input = pmFPAfileSelectSingle(config->files, "PPMERGE.INPUT", i); // File of interest + pmFPA *fpa = input->fpa; // FPA of interest + CHECK_LEVEL(fpa->hdu, -1, -1); + psArray *chips = fpa->chips; // Array of chips + for (int j = 0; j < chips->n; j++) { + pmChip *chip = chips->data[j]; // Chip of interest + CHECK_LEVEL(chip->hdu, j, -1); + psArray *cells = chip->cells; // Array of cells + for (int k = 0; k < cells->n; k++) { + pmCell *cell = cells->data[k]; // Cell of interest + CHECK_LEVEL(cell->hdu, j, k); + } + } + } + if (!phuView || !format) { + psError(PS_ERR_UNEXPECTED_NULL, true, "Unable to find PHU for input files."); + psFree(phuView); + return false; + } + + // Cull chips and cells that don't have data + // Otherwise the abundance of metadata in the concepts (esp. for GPC) can overload the memory + for (int i = 0; i < numFiles; i++) { + pmFPAfile *input = pmFPAfileSelectSingle(config->files, "PPMERGE.INPUT", i); // File of interest + pmFPA *fpa = input->fpa; // FPA of interest + 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; + cell->data_exists = false; + cell->file_exists = false; + culled++; + } + } + if (culled == cells->n) { + psFree(chip->concepts); + chip->concepts = NULL; + chip->data_exists = false; + chip->file_exists = false; + } + } + } + + // Count the cells + { + int numCells = 0; // Number of cells + pmFPAfile *input = pmFPAfileSelectSingle(config->files, "PPMERGE.INPUT", 0); // Representative file + pmFPA *fpa = input->fpa; // FPA for file + psArray *chips = fpa->chips; // Array of chips + for (int i = 0; i < chips->n; i++) { + pmChip *chip = chips->data[i]; // Chip of interest + psArray *cells = chip->cells; // Array of cells + for (int j = 0; j < cells->n; j++) { + pmCell *cell = cells->data[j]; // Cell of interest + pmHDU *hdu = pmHDUGetLowest(fpa, chip, cell); // HDU that would have data + if (hdu && !hdu->blankPHU) { + numCells++; + } + } + } + + psMetadataAddS32(config->arguments, PS_LIST_TAIL, "INPUTS.CELLS", 0, "Number of cells in input", + numCells); + } + + // Output image + pmFPA *fpa = pmFPAConstruct(config->camera); // FPA to contain the output + if (!fpa) { + psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to construct an FPA from camera configuration."); + psFree(phuView); + return false; + } + + psString outName = ppMergeOutputFile(config); // Name of output file + + pmFPAfileType fileType = PM_FPA_FILE_NONE; // Type of output file + switch (type) { + case PPMERGE_TYPE_BIAS: + fileType = PM_FPA_FILE_IMAGE; + break; + case PPMERGE_TYPE_DARK: + fileType = PM_FPA_FILE_DARK; + break; + case PPMERGE_TYPE_MASK: + fileType = PM_FPA_FILE_MASK; + break; + case PPMERGE_TYPE_SHUTTER: + fileType = PM_FPA_FILE_IMAGE; + break; + case PPMERGE_TYPE_FLAT: + fileType = PM_FPA_FILE_IMAGE; + break; + case PPMERGE_TYPE_FRINGE: + fileType = PM_FPA_FILE_FRINGE; + default: + psAbort("Unknown frame type: %x", type); + } + + if (!outputFile(config, outName, fileType, "Merged detrend", format, phuView)) { + psFree(outName); + psFree(phuView); + return false; + } + psFree(outName); + + if (!outputFile(config, "PPMERGE.OUTPUT.SIGMA", PM_FPA_FILE_IMAGE, "Merge sigma", format, phuView)) { + psFree(phuView); + return false; + } + + if (!outputFile(config, "PPMERGE.OUTPUT.COUNT", PM_FPA_FILE_IMAGE, "Merged count", format, phuView)) { + psFree(phuView); + return false; + } + + psFree(phuView); + + return true; +} Index: /trunk/ppMerge/src/ppMergeFiles.c =================================================================== --- /trunk/ppMerge/src/ppMergeFiles.c (revision 17227) +++ /trunk/ppMerge/src/ppMergeFiles.c (revision 17227) @@ -0,0 +1,214 @@ +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include +#include +#include + +#include "ppMerge.h" + +const char *allFiles[] = { "PPMERGE.INPUT", "PPMERGE.INPUT.MASK", "PPMERGE.INPUT.WEIGHT", + "PPMERGE.OUTPUT", "PPMERGE.OUTPUT.COUNT", "PPMERGE.OUTPUT.SIGMA", + NULL }; // All files +const char *inputFiles[] = { "PPMERGE.INPUT", "PPMERGE.INPUT.MASK", "PPMERGE.INPUT.WEIGHT", + NULL }; // Input files +const char *outputFiles[] = { "PPMERGE.OUTPUT", "PPMERGE.OUTPUT.COUNT", "PPMERGE.OUTPUT.SIGMA", + NULL }; // Output files + +// Select file list based on enum +static const char **selectFiles(ppMergeFiles files) +{ + switch (files) { + case PPMERGE_FILES_ALL: return allFiles; + case PPMERGE_FILES_INPUT: return inputFiles; + case PPMERGE_FILES_OUTPUT: return outputFiles; + default: + psAbort("Invalid file option"); + } + return NULL; +} + +bool ppMergeFileReadInput(const pmConfig *config, pmReadout *readout, int num, int rows) +{ + pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PPMERGE.INPUT", num); + if (!pmReadoutReadChunk(readout, file->fits, 0, rows, 0)) { + psError(PS_ERR_UNKNOWN, false, "Unable to read readout."); + return false; + } + bool mdok; // Status of MD lookup + if (psMetadataLookupBool(&mdok, config->arguments, "INPUTS.MASKS")) { + pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PPMERGE.INPUT.MASK", num); + if (!pmReadoutReadChunkMask(readout, file->fits, 0, rows, 0)) { + psError(PS_ERR_UNKNOWN, false, "Unable to read readout mask."); + return false; + } + } + if (psMetadataLookupBool(&mdok, config->arguments, "INPUTS.WEIGHTS")) { + pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PPMERGE.INPUT.WEIGHT", num); + if (!pmReadoutReadChunkWeight(readout, file->fits, 0, rows, 0)) { + psError(PS_ERR_UNKNOWN, false, "Unable to read readout weight."); + return false; + } + } + return true; +} + +bool ppMergeFileOpenInput(pmConfig *config, const pmFPAview *view, int num) +{ + { + pmFPAfile *input = pmFPAfileSelectSingle(config->files, "PPMERGE.INPUT", num); + pmFPAview *fileView = pmFPAviewForLevel(input->fileLevel, view); + if (!pmFPAfileOpen(input, fileView, config)) { + psError(PS_ERR_UNKNOWN, false, "Unable to open image file %d", num); + psFree(fileView); + return false; + } + psFree(fileView); + } + bool mdok; // Status of MD lookup + if (psMetadataLookupBool(&mdok, config->arguments, "INPUTS.MASKS")) { + pmFPAfile *mask = pmFPAfileSelectSingle(config->files, "PPMERGE.INPUT.MASK", num); // Mask file + pmFPAview *fileView = pmFPAviewForLevel(mask->fileLevel, view); + if (!pmFPAfileOpen(mask, fileView, config)) { + psError(PS_ERR_UNKNOWN, false, "Unable to open mask file %d", num); + psFree(fileView); + return false; + } + psFree(fileView); + } + if (psMetadataLookupBool(&mdok, config->arguments, "INPUTS.WEIGHTS")) { + pmFPAfile *weight = pmFPAfileSelectSingle(config->files, "PPMERGE.INPUT.WEIGHT", num); // Weight file + pmFPAview *fileView = pmFPAviewForLevel(weight->fileLevel, view); + if (!pmFPAfileOpen(weight, fileView, config)) { + psError(PS_ERR_UNKNOWN, false, "Unable to open weight file %d", num); + psFree(fileView); + return false; + } + psFree(fileView); + } + return true; +} + +psArray *ppMergeFileDataLevel(const pmConfig *config, const char *name, pmFPALevel level) +{ + assert(config); + assert(name); + + int numFiles = psMetadataLookupS32(NULL, config->arguments, "INPUTS.NUM"); // Number of input files + assert(numFiles > 0); + psArray *files = psArrayAlloc(numFiles); // Files of interest + for (int i = 0; i < numFiles; i++) { + pmFPAfile *file = pmFPAfileSelectSingle(config->files, name, i); // Image file + file->dataLevel = level; + file->freeLevel = level; + files->data[i] = psMemIncrRefCounter(file); + } + return files; +} + +bool ppMergeFileActivate(const pmConfig *config, ppMergeFiles files, bool state) +{ + assert(config); + + bool mdok; // Status of MD lookup + bool haveMasks = psMetadataLookupBool(&mdok, config->arguments, "INPUTS.MASKS"); // Do we have masks? + bool haveWeights = psMetadataLookupBool(&mdok, config->arguments, "INPUTS.WEIGHTS"); // Got weights? + + const char **fileList = selectFiles(files); // Files to activate + for (int i = 0; fileList[i] != NULL; i++) { + if (!haveMasks && strcmp(fileList[i], "PPMERGE.INPUT.MASK") == 0) { + continue; + } + if (!haveWeights && strcmp(fileList[i], "PPMERGE.INPUT.WEIGHT") == 0) { + continue; + } + psString name = NULL; // Name of file + if (strcmp(fileList[i], "PPMERGE.OUTPUT") == 0) { + name = ppMergeOutputFile(config); + } else { + name = psStringCopy(fileList[i]); + } + + if (!pmFPAfileActivate(config->files, state, name)) { + psError(PS_ERR_UNKNOWN, false, "Unable to activate file %s", name); + psFree(name); + return false; + } + psFree(name); + } + + return true; +} + + + +psArray *ppMergeFileActivateSingle(const pmConfig *config, ppMergeFiles files, bool state, int num) +{ + assert(config); + + bool mdok; // Status of MD lookup + bool haveMasks = psMetadataLookupBool(&mdok, config->arguments, "INPUTS.MASKS"); // Do we have masks? + bool haveWeights = psMetadataLookupBool(&mdok, config->arguments, "INPUTS.WEIGHTS"); // Got weights? + + psList *list = psListAlloc(NULL); // List of files + const char **fileList = selectFiles(files); // Files to activate + for (int i = 0; fileList[i] != NULL; i++) { + if (!haveMasks && strcmp(fileList[i], "PPMERGE.INPUT.MASK") == 0) { + continue; + } + if (!haveWeights && strcmp(fileList[i], "PPMERGE.INPUT.WEIGHT") == 0) { + continue; + } + + psString name = NULL; // Name of file + if (strcmp(fileList[i], "PPMERGE.OUTPUT") == 0) { + name = ppMergeOutputFile(config); + } else { + name = psStringCopy(fileList[i]); + } + + pmFPAfile *file = pmFPAfileActivateSingle(config->files, state, name, num); // Activated file + psFree(name); + psListAdd(list, PS_LIST_TAIL, file); + } + + psArray *array = psListToArray(list); + psFree(list); + + return array; +} + + +psString ppMergeOutputFile(const pmConfig *config) +{ + ppMergeType type = psMetadataLookupS32(NULL, config->arguments, "TYPE"); // Type of frame + const char *outSuffix = NULL; // Suffix for output file + switch (type) { + case PPMERGE_TYPE_BIAS: + outSuffix = "BIAS"; + break; + case PPMERGE_TYPE_DARK: + outSuffix = "DARK"; + break; + case PPMERGE_TYPE_MASK: + outSuffix = "MASK"; + break; + case PPMERGE_TYPE_SHUTTER: + outSuffix = "SHUTTER"; + break; + case PPMERGE_TYPE_FLAT: + outSuffix = "FLAT"; + break; + case PPMERGE_TYPE_FRINGE: + outSuffix = "FRINGE"; + default: + psAbort("Unknown frame type: %x", type); + } + + psString outName = NULL; // Name of output file + psStringAppend(&outName, "PPMERGE.OUTPUT.%s", outSuffix); + + return outName; +} Index: /trunk/ppMerge/src/ppMergeLoop.c =================================================================== --- /trunk/ppMerge/src/ppMergeLoop.c (revision 17227) +++ /trunk/ppMerge/src/ppMergeLoop.c (revision 17227) @@ -0,0 +1,346 @@ +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include + +#include +#include +#include + +#include "ppMerge.h" + +bool ppMergeLoop(pmConfig *config) +{ + assert(config); + + psMetadata *arguments = config->arguments; // Arguments + ppMergeType type = psMetadataLookupS32(NULL, config->arguments, "TYPE"); // Type of frame + int numFiles = psMetadataLookupS32(NULL, arguments, "INPUTS.NUM"); // Number of input files + bool mdok; // Status of MD lookup + bool haveMasks = psMetadataLookupBool(&mdok, arguments, "INPUTS.MASKS"); // Do we have masks? + bool haveWeights = psMetadataLookupBool(&mdok, arguments, "INPUTS.WEIGHTS"); // Do we have weights? + + psArray *inputs = ppMergeFileDataLevel(config, "PPMERGE.INPUT", PM_FPA_LEVEL_READOUT); // Input images + psArray *masks = NULL, *weights = NULL; // Input masks and weights + if (haveMasks) { + masks = ppMergeFileDataLevel(config, "PPMERGE.INPUT.MASK", PM_FPA_LEVEL_READOUT); + } + if (haveWeights) { + weights = ppMergeFileDataLevel(config, "PPMERGE.INPUT.WEIGHT", PM_FPA_LEVEL_READOUT); + } + + // General combination parameters + int rows = psMetadataLookupS32(NULL, arguments, "ROWS"); // Number of rows to read per chunk + int iter = psMetadataLookupS32(NULL, arguments, "ITER"); // Number of rejection iterations + float rej = psMetadataLookupF32(NULL, arguments, "REJ"); // Rejection level + float fraclow = psMetadataLookupF32(NULL, arguments, "FRACLOW"); // Reject fraction of low pixels + float frachigh = psMetadataLookupF32(NULL, arguments, "FRACHIGH"); // Reject fraction of hi pixels + int nKeep = psMetadataLookupS32(NULL, arguments, "NKEEP"); // Minimum number of values to keep + psMaskType maskVal = psMetadataLookupU8(NULL, arguments, "MASKVAL"); // Value to mask + psStatsOptions combineStat = psMetadataLookupS32(NULL, arguments, "COMBINE"); // Combination statistic + bool useWeights = psMetadataLookupBool(NULL, arguments, "WEIGHTS"); // Use weights? + + // Fringe parameters + int fringeNum = psMetadataLookupS32(NULL, arguments, "FRINGE.NUM"); // Number of fringe points + int fringeSize = psMetadataLookupS32(NULL, arguments, "FRINGE.SIZE"); // Size of fringe regions + int fringeSmoothX = psMetadataLookupS32(NULL, arguments, "FRINGE.XSMOOTH"); // Smoothing regions in x + int fringeSmoothY = psMetadataLookupS32(NULL, arguments, "FRINGE.YSMOOTH"); // Smoothing regions in y + + pmCombineParams *combination = pmCombineParamsAlloc(combineStat); // Combination parameters + combination->maskVal = maskVal; + combination->blank = pmConfigMask("BLANK", config); + combination->nKeep = nKeep; + combination->fracHigh = frachigh; + combination->fracLow = fraclow; + combination->iter = iter; + combination->rej = rej; + combination->weights = useWeights; + + // Retrieve data placed on analysis + psVector *scales = NULL, *zeros = NULL; // Scale and zeroes for combination + psArray *shutters = NULL; // Shutter correction data + switch (type) { + case PPMERGE_TYPE_FRINGE: + zeros = psMetadataLookupPtr(NULL, arguments, "ZEROS"); + if (!zeros) { + psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to find ZEROS"); + goto ERROR; + } + // Flow through + case PPMERGE_TYPE_FLAT: + scales = psMetadataLookupPtr(NULL, arguments, "SCALES"); + if (!scales) { + psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to find SCALES"); + goto ERROR; + } + break; + case PPMERGE_TYPE_SHUTTER: + shutters = psMetadataLookupPtr(NULL, arguments, "SHUTTER"); + if (!shutters) { + psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to find SHUTTER"); + goto ERROR; + } + break; + case PPMERGE_TYPE_MASK: + case PPMERGE_TYPE_BIAS: + case PPMERGE_TYPE_DARK: + break; + default: + psAbort("Should never get here."); + } + + // Dark parameters + psArray *darkOrdinates = psMetadataLookupPtr(NULL, arguments, "DARK.ORDINATES"); // Dark info + psString darkNorm = psMetadataLookupStr(&mdok, arguments, "DARK.NORM"); // Dark normalisation + + psMetadata *stats = NULL; // Statistics for output + if (psMetadataLookup(config->arguments, "STATS.NAME")) { + stats = psMetadataAlloc(); + psMetadataAddMetadata(config->arguments, PS_LIST_TAIL, "STATS.DATA", 0, "Statistics output", stats); + } + + if (!ppMergeFileActivate(config, PPMERGE_FILES_ALL, true)) { + psError(PS_ERR_UNKNOWN, false, "Unable to activate files."); + goto ERROR; + } + psString outName = ppMergeOutputFile(config); // Name of output file + pmFPAfile *output = psMetadataLookupPtr(NULL, config->files, outName); // Output file + psFree(outName); + assert(output && output->fpa); + pmFPA *outFPA = output->fpa; // Output FPA + pmFPAview *view = pmFPAviewAlloc(0); // View to component of interest + int cellNum = 0; // Index of cell + if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) { + goto ERROR; + } + pmChip *outChip; // Chip of interest + while ((outChip = pmFPAviewNextChip(view, outFPA, 1))) { + if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) { + goto ERROR; + } + pmCell *outCell; // Cell of interest + while ((outCell = pmFPAviewNextCell(view, outFPA, 1))) { + if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) { + goto ERROR; + } + + pmHDU *hdu = pmHDUGetLowest(outFPA, outChip, outCell); // HDU for cell + if (!hdu || hdu->blankPHU) { + // No data here + continue; + } + + pmReadout *outRO = pmReadoutAlloc(outCell); + + psArray *readouts = psArrayAlloc(numFiles); // Input readouts + for (int i = 0; i < numFiles; i++) { + // We need to do some of the opening ourselves + if (!ppMergeFileOpenInput(config, view, i)) { + psError(PS_ERR_UNKNOWN, false, "Unable to open file %d", i); + goto ERROR; + } + + pmFPAfile *input = pmFPAfileSelectSingle(config->files, "PPMERGE.INPUT", i); + pmCell *inCell = pmFPAviewThisCell(view, input->fpa); // Input cell + readouts->data[i] = pmReadoutAlloc(inCell); + } + + float shutterRef = NAN; // Reference shutter correction + if (type == PPMERGE_TYPE_SHUTTER) { + shutterRef = pmShutterCorrectionReference(shutters->data[cellNum]); + } + + // Read convolutions by chunks + bool more = true; // More to read? + for (int numChunk = 0; more; numChunk++) { + psTrace("ppStack", 2, "Initial stack of chunk %d....\n", numChunk); + for (int i = 0; i < numFiles; i++) { + pmReadout *inRO = readouts->data[i]; // Input readout + + // Read a chunk from a file + #define READ_CHUNK(NAME,TYPE) { \ + pmFPAfile *file = pmFPAfileSelectSingle(config->files, NAME, i); \ + if (!pmReadoutReadChunk##TYPE(inRO, file->fits, 0, rows, 0)) { \ + psError(PS_ERR_IO, false, "Unable to read chunk %d for file %s %d", \ + numChunk, NAME, i); \ + psFree(readouts); \ + psFree(outRO); \ + goto ERROR; \ + } \ + } + + READ_CHUNK("PPMERGE.INPUT", /* Blank */); + if (haveMasks) { + READ_CHUNK("PPMERGE.INPUT.MASK", Mask); + } + if (haveWeights) { + READ_CHUNK("PPMERGE.INPUT.WEIGHT", Weight); + } + } + + switch (type) { + case PPMERGE_TYPE_SHUTTER: + if (!pmShutterCorrectionGenerate(outRO, NULL, readouts, shutterRef, + shutters->data[cellNum], iter, rej, maskVal)) { + psFree(readouts); + psFree(outRO); + goto ERROR; + } + break; + case PPMERGE_TYPE_DARK: + if (!pmDarkCombine(outCell, readouts, darkOrdinates, darkNorm, iter, rej, maskVal)) { + psFree(readouts); + psFree(outRO); + goto ERROR; + } + break; + case PPMERGE_TYPE_BIAS: + case PPMERGE_TYPE_FLAT: + case PPMERGE_TYPE_FRINGE: + if (!pmReadoutCombine(outRO, readouts, zeros, scales, combination)) { + psFree(readouts); + psFree(outRO); + goto ERROR; + } + break; + default: + psAbort("Should never get here."); + } + + + for (int i = 0; i < numFiles && more; i++) { + pmReadout *inRO = readouts->data[i]; + + // Check to see if there's more chunks to read + #define MORE_CHUNK(NAME,TYPE) { \ + pmFPAfile *file = pmFPAfileSelectSingle(config->files, NAME, i); \ + more &= pmReadoutMore##TYPE(inRO, file->fits, 0, rows); \ + } + + MORE_CHUNK("PPMERGE.INPUT", /* Blank */); + if (haveMasks) { + MORE_CHUNK("PPMERGE.INPUT.MASK", Mask); + } + if (haveWeights) { + MORE_CHUNK("PPMERGE.INPUT.WEIGHT", Weight); + } + } + } + + // Get list of cells for concepts averaging + psList *inCells = psListAlloc(NULL); // List of cells + for (int i = 0; i < numFiles; i++) { + pmReadout *readout = readouts->data[i]; // Readout of interest + psListAdd(inCells, PS_LIST_TAIL, readout->parent); + } + if (!pmConceptsAverageCells(outCell, inCells, NULL, NULL, true)) { + psError(PS_ERR_UNKNOWN, false, "Unable to average cell concepts."); + psFree(inCells); + psFree(outRO); + goto ERROR; + } + psFree(inCells); + + psFree(readouts); + + // Plug supplementary images into their own FPAs + { + pmCell *countsCell = pmFPAfileThisCell(config->files, view, "PPMERGE.OUTPUT.COUNT"); + pmReadout *countsRO = pmReadoutAlloc(countsCell); // Readout with count of inputs per pixel + psImage *counts = psMetadataLookupPtr(NULL, outRO->analysis, PM_READOUT_STACK_ANALYSIS_COUNT); + countsRO->image = psImageCopy(countsRO->image, counts, PS_TYPE_F32); + psMetadataRemoveKey(outRO->analysis, PM_READOUT_STACK_ANALYSIS_COUNT); + countsRO->data_exists = countsCell->data_exists = countsCell->parent->data_exists = true; + psFree(countsRO); + + pmCell *sigmaCell = pmFPAfileThisCell(config->files, view, "PPMERGE.OUTPUT.SIGMA"); + pmReadout *sigmaRO = pmReadoutAlloc(sigmaCell); // Readout with stdev per pixel + psImage *sigma = psMetadataLookupPtr(NULL, outRO->analysis, PM_READOUT_STACK_ANALYSIS_SIGMA); + sigmaRO->image = psImageCopy(sigmaRO->image, sigma, PS_TYPE_F32); + psMetadataRemoveKey(outRO->analysis, PM_READOUT_STACK_ANALYSIS_SIGMA); + sigmaRO->data_exists = sigmaCell->data_exists = sigmaCell->parent->data_exists = true; + psFree(sigmaRO); + } + + // 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 (type == PPMERGE_TYPE_FRINGE && outRO) { + pmFringeRegions *regions = pmFringeRegionsAlloc(fringeNum, fringeSize, fringeSize, + fringeSmoothX, fringeSmoothY); + pmFringeStats *fringe = pmFringeStatsMeasure(regions, outRO, maskVal); + psFree(regions); + if (!fringe) { + psError(PS_ERR_UNKNOWN, false, "Unable to measure fringe statistics.\n"); + psFree(outRO); + goto ERROR; + } + + psArray *fringes = psArrayAlloc(1); // Array of fringes + fringes->data[0] = fringe; + + pmFringesFormat(outCell, NULL, fringes); + psFree(fringes); // Drop reference + } + + if (!ppStatsFPA(stats, outFPA, view, maskVal | pmConfigMask("BLANK", config), config)) { + psError(PS_ERR_UNKNOWN, true, "Unable to generate stats for image."); + goto ERROR; + } + + psFree(outRO); + cellNum++; + + if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) { + goto ERROR; + } + } + if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) { + goto ERROR; + } + } + + // Get list of FPAs for concepts averaging + psList *inFPAs = psListAlloc(NULL); // List of FPAs + for (int i = 0; i < numFiles; i++) { + pmFPAfile *input = inputs->data[i]; // Input file + psListAdd(inFPAs, PS_LIST_TAIL, input->fpa); + } + if (!pmConceptsAverageFPAs(output->fpa, inFPAs)) { + psError(PS_ERR_UNKNOWN, false, "Unable to average FPA concepts."); + psFree(inFPAs); + goto ERROR; + } + psFree(inFPAs); + + + if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) { + goto ERROR; + } + + psFree(view); + psFree(combination); + psFree(inputs); + psFree(masks); + psFree(weights); + psFree(stats); + return true; + +ERROR: + psFree(view); + psFree(combination); + psFree(inputs); + psFree(masks); + psFree(weights); + psFree(stats); + return false; +} + Index: /trunk/ppMerge/src/ppMergeMask.c =================================================================== --- /trunk/ppMerge/src/ppMergeMask.c (revision 17226) +++ /trunk/ppMerge/src/ppMergeMask.c (revision 17227) @@ -5,6 +5,6 @@ #include #include -#include #include + #include #include @@ -12,203 +12,394 @@ #include "ppMerge.h" -#include "ppMergeData.h" -#include "ppMergeMask.h" - - -// Generate a mask -bool ppMergeMask(ppMergeData *data, // Data - ppMergeOptions *options, // Options - pmConfig *config // Configuration - ) + +static bool mergeMask(pmConfig *config, // Configuration + const pmFPAview *view, // View to chip + bool writeOut, // Write output? + psRandom *rng, // Random number generator + psMetadata *stats // Statistics output + ) { - for (int i = 0; i < 2; i++) { - if (!ppMergeMaskSuspect (data, options, config)) { - psError(PS_ERR_UNKNOWN, true, "failed on mask suspect.\n"); - return false; - } - - if (!ppMergeMaskBad (data, options, config)) { - psError(PS_ERR_UNKNOWN, true, "failed on mask bad.\n"); - return false; - } - } - - if (!ppMergeMaskAverageConcepts (data, options, config)) { - psError(PS_ERR_UNKNOWN, true, "failed on average concepts.\n"); - return false; - } - - if (!ppMergeMaskGrow (data, options, config)) { - psError(PS_ERR_UNKNOWN, true, "failed on mask grow.\n"); - return false; - } - - if (!ppMergeMaskWrite (data, options, config)) { - psError(PS_ERR_UNKNOWN, true, "failed on mask write.\n"); - return false; - } + assert(config); + assert(view); + assert(view->chip != -1 && view->cell == -1 && view->readout == -1); + + bool mdok; // Status of MD lookup + int numFiles = psMetadataLookupS32(NULL, config->arguments, "INPUTS.NUM"); // Number of input files + psStatsOptions meanStat = psMetadataLookupS32(NULL, config->arguments, "MEAN"); // Statistic for mean + psStatsOptions stdevStat = psMetadataLookupS32(NULL, config->arguments, "STDEV"); // Statistic for stdev + psMaskType maskVal = psMetadataLookupU8(NULL, config->arguments, "MASKVAL"); // Value to mask + int sample = psMetadataLookupS32(NULL, config->arguments, "SAMPLE"); // Size of sample for statistics + bool chipStats = psMetadataLookupBool(&mdok, config->arguments, "MASK.CHIPSTATS"); // Statistics on chip? + float maskSuspect = psMetadataLookupF32(NULL, config->arguments, "MASK.SUSPECT"); // Threshold for suspect pixels + float maskBad = psMetadataLookupF32(NULL, config->arguments, "MASK.BAD"); // Threshold for bad pixels + pmMaskIdentifyMode maskMode = psMetadataLookupS32(NULL, config->arguments, "MASK.MODE"); // Mode for identifying bad pixels + int maskGrow = psMetadataLookupS32(NULL, config->arguments, "MASK.GROW"); // Radius to grow mask + psMaskType maskGrowVal = psMetadataLookupU8(NULL, config->arguments, "MASK.GROWVAL"); // Value for grown mask + + psStats *statistics = psStatsAlloc(meanStat | stdevStat); // Statistics for background + + psString outName = ppMergeOutputFile(config); // Name of output file + pmChip *outChip = pmFPAfileThisChip(config->files, view, outName); // Output chip + psFree(outName); + + // For each input file, get the statistics, which can be calculated at the chip or cell levels + psVector *values = psVectorAlloc(sample, PS_TYPE_F32); // Pixel values for statistics + pmFPAview *inView = pmFPAviewAlloc(0); // View for input + for (int i = 0; i < numFiles; i++) { + pmFPAfileActivate(config->files, false, NULL); + psArray *files = ppMergeFileActivateSingle(config, PPMERGE_FILES_INPUT, true, i); // Input files + pmFPAfile *input = files->data[0]; // Input file + psFree(files); + pmFPA *inFPA = input->fpa; // Input FPA + *inView = *view; + + int valueIndex = 0; // Index for vector of pixel values + + pmCell *inCell; // Input cell + while ((inCell = pmFPAviewNextCell(inView, inFPA, 1))) { + pmHDU *hdu = pmHDUFromCell(inCell); // HDU for cell + if (!hdu || hdu->blankPHU) { + // No data here + continue; + } + psTrace("ppMerge", 1, "Getting suspect pixels for file %d chip %d cell %d", + i, inView->chip, inView->cell); + + if (!pmFPAfileIOChecks(config, inView, PM_FPA_BEFORE)) { + psFree(inView); + goto MERGE_MASK_ERROR; + } + + if (!ppMergeFileOpenInput(config, view, i)) { + psError(PS_ERR_UNKNOWN, false, "Unable to open file %d", i); + psFree(inView); + goto MERGE_MASK_ERROR; + } + + if (inCell->readouts->n > 1) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, + "File %d chip %d cell %d contains more than one readout (%ld)", + i, inView->chip, inView->cell, inCell->readouts->n); + psFree(inView); + goto MERGE_MASK_ERROR; + } + + pmReadout *readout; + if (inCell->readouts && inCell->readouts->n == 1) { + readout = psMemIncrRefCounter(inCell->readouts->data[0]); // Input readout + } else { + readout = pmReadoutAlloc(inCell); + } + + if (!ppMergeFileReadInput(config, readout, i, 0)) { + psError(PS_ERR_UNKNOWN, false, "Unable to read readout %d", i); + psFree(inView); + psFree(readout); + goto MERGE_MASK_ERROR; + } + + pmCell *outCell = pmFPAfileThisCell(config->files, inView, "PPMERGE.OUTPUT.MASK"); // Output cell + pmReadout *outRO = NULL; // Output readout + if (outCell->readouts && outCell->readouts->n == 1) { + outRO = psMemIncrRefCounter(outCell->readouts->data[0]); + } else { + outRO = pmReadoutAlloc(outCell); + } + psImage *outMask = outRO->mask; // Output mask image (for iterative generation of mask) + + psImage *image = readout->image, *mask = readout->mask; // Image and mask + int numCols = readout->image->numCols, numRows = readout->image->numRows; // Image size + int numPix = numCols * numRows; // Number of pixels + int numCells = chipStats ? readout->parent->parent->cells->n : 1; // Number of cells + int num = PS_MIN(numPix, sample / numCells); // Number of values to add + if (!chipStats) { + valueIndex = 0; + } + for (int i = 0; i < num; i++) { + int pixel = numPix * psRandomUniform(rng); + int x = pixel % numCols; + int y = pixel / numCols; + if ((mask && (mask->data.PS_TYPE_MASK_DATA[y][x] & maskVal)) || + !isfinite(image->data.F32[y][x]) || + (outMask && (outMask->data.PS_TYPE_MASK_DATA[y][x] & maskVal))) { + continue; + } + + values->data.F32[valueIndex++] = image->data.F32[y][x]; + } + + if (!chipStats) { + values->n = valueIndex; + if (!psVectorStats(statistics, values, NULL, NULL, 0)) { + psError(PS_ERR_UNKNOWN, false, "Unable to do statistics on readout."); + psFree(inView); + psFree(readout); + goto MERGE_MASK_ERROR; + } + + if (!pmMaskFlagSuspectPixels(outRO, readout, psStatsGetValue(statistics, meanStat), + psStatsGetValue(statistics, stdevStat), maskSuspect, maskVal)) { + psError(PS_ERR_UNKNOWN, false, "Unable to find suspect values in file %d", i); + psFree(inView); + psFree(readout); + goto MERGE_MASK_ERROR; + } + pmCellFreeData(inCell); + + if (!pmFPAfileIOChecks(config, inView, PM_FPA_AFTER)) { + psFree(inView); + psFree(readout); + goto MERGE_MASK_ERROR; + } + } + psFree(readout); + psFree(outRO); + } + + // Additional run through cells if we want chip-level statistics + if (chipStats && valueIndex > 0) { + values->n = valueIndex; + if (!psVectorStats(statistics, values, NULL, NULL, 0)) { + psError(PS_ERR_UNKNOWN, false, "Unable to do statistics on chip."); + goto MERGE_MASK_ERROR; + } + inView->cell = -1; + while ((inCell = pmFPAviewNextCell(inView, inFPA, 1))) { + pmHDU *hdu = pmHDUFromCell(inCell); // HDU for cell + if (!hdu || hdu->blankPHU) { + // No data here + continue; + } + pmReadout *readout = inCell->readouts->data[0]; // Readout of interest + + inView->readout = 0; + pmReadout *outRO = pmFPAfileThisReadout(config->files, inView, "PPMERGE.OUTPUT.MASK"); + + if (!pmMaskFlagSuspectPixels(outRO, readout, psStatsGetValue(statistics, meanStat), + psStatsGetValue(statistics, stdevStat), maskSuspect, maskVal)) { + psError(PS_ERR_UNKNOWN, false, "Unable to find suspect values in file %d", i); + goto MERGE_MASK_ERROR; + } + pmCellFreeData(inCell); + + inView->readout = -1; + if (!pmFPAfileIOChecks(config, inView, PM_FPA_AFTER)) { + psFree(inView); + goto MERGE_MASK_ERROR; + } + } + } + } + psFree(inView); + psFree(statistics); statistics = NULL; + psFree(values); values = NULL; + + + // Another run through the chip to threshold on the suspects + pmFPAfileActivate(config->files, false, NULL); + if (writeOut) { + ppMergeFileActivate(config, PPMERGE_FILES_OUTPUT, true); + } + pmFPA *outFPA = outChip->parent; // Output FPA + pmCell *outCell; // Output cell + pmFPAview *outView = pmFPAviewAlloc(0); // View into output FPA + *outView = *view; + while ((outCell = pmFPAviewNextCell(outView, outFPA, 1))) { + pmHDU *hdu = pmHDUFromCell(outCell); // HDU for cell + if (!hdu || hdu->blankPHU) { + // No data here + continue; + } + + psTrace("ppMerge", 1, "Getting bad pixels for chip %d cell %d", outView->chip, outView->cell); + + assert(outCell->readouts && outCell->readouts->n == 1); + pmReadout *outRO = outCell->readouts->data[0]; // Output readout + if (!pmMaskIdentifyBadPixels(outRO, maskVal, maskBad, maskMode)) { + psError(PS_ERR_UNKNOWN, false, "Unable to mask bad pixels"); + goto MERGE_MASK_ERROR; + } + + // Supplementary outputs + { + // The counts image is fairly useless, but it preserves the model + pmCell *countsCell = pmFPAfileThisCell(config->files, outView, "PPMERGE.OUTPUT.COUNT"); + pmReadout *countsRO = pmReadoutAlloc(countsCell); // Readout with count of inputs per pixel + countsRO->image = psImageAlloc(outRO->mask->numCols, outRO->mask->numRows, PS_TYPE_F32); + psImageInit(countsRO->image, numFiles); + countsRO->data_exists = countsCell->data_exists = countsCell->parent->data_exists = true; + psFree(countsRO); + + pmCell *sigmaCell = pmFPAfileThisCell(config->files, outView, "PPMERGE.OUTPUT.SIGMA"); + pmReadout *sigmaRO = pmReadoutAlloc(sigmaCell); // Readout with suspect image + psImage *suspect = psMetadataLookupPtr(NULL, outRO->analysis, PM_MASK_ANALYSIS_SUSPECT); + sigmaRO->image = psImageCopy(sigmaRO->image, suspect, PS_TYPE_F32); + psMetadataRemoveKey(outRO->analysis, PM_MASK_ANALYSIS_SUSPECT); + sigmaRO->data_exists = sigmaCell->data_exists = sigmaCell->parent->data_exists = true; + psFree(sigmaRO); + } + + if (maskGrowVal > 0) { + psImage *grown = psImageGrowMask(NULL, outRO->mask, maskVal, maskGrow, maskGrowVal); // Grown mask + psFree(outRO->mask); + outRO->mask = grown; + } + + if (writeOut) { + if (!pmFPAfileIOChecks(config, outView, PM_FPA_BEFORE)) { + psFree(outView); + goto MERGE_MASK_ERROR; + } + + outRO->data_exists = outCell->data_exists = outChip->data_exists = true; + + // Average concepts + psList *cells = psListAlloc(NULL); // List of cells, for concept averaging + for (int i = 0; i < numFiles; i++) { + pmFPAfile *inFile = pmFPAfileSelectSingle(config->files, "PPMERGE.INPUT", i); // Input file + pmCell *inCell = pmFPAviewThisCell(outView, inFile->fpa); // Input cell + psListAdd(cells, PS_LIST_TAIL, inCell); + } + if (!pmConceptsAverageCells(outCell, cells, NULL, NULL, true)) { + psError(PS_ERR_UNKNOWN, false, "Unable to average cell concepts."); + psFree(cells); + psFree(outRO); + psFree(outView); + return false; + } + psFree(cells); + + // Statistics on the merged cell using a fake image + outRO->image = psImageAlloc(outRO->mask->numCols, outRO->mask->numRows, PS_TYPE_F32); + psImageInit(outRO->image, 1.0); + if (!ppStatsFPA(stats, outRO->parent->parent->parent, outView, + maskVal | pmConfigMask("BLANK", config), config)) { + psError(PS_ERR_UNEXPECTED_NULL, true, "Unable to generate stats for image."); + psFree(outRO); + psFree(outView); + return false; + } + psFree(outRO->image); + outRO->image = NULL; + + // Write + if (!pmFPAfileIOChecks(config, outView, PM_FPA_AFTER)) { + psFree(outView); + goto MERGE_MASK_ERROR; + } + } + } + psFree(outView); + + ppMergeFileActivate(config, PPMERGE_FILES_ALL, true); return true; + + +MERGE_MASK_ERROR: + psFree(statistics); + psFree(values); + return false; } -bool ppMergeMaskReadoutStats(const pmReadout *readout, - const pmReadout *output, - ppMergeOptions *options, // Options - psRandom *rng) + + + + +bool ppMergeMask(pmConfig *config) { - PS_ASSERT_PTR_NON_NULL(readout, false); - PS_ASSERT_IMAGE_NON_NULL(readout->image, false); - PS_ASSERT_IMAGE_NON_EMPTY(readout->image, false); - PS_ASSERT_IMAGE_TYPE(readout->image, PS_TYPE_F32, false); - if (readout->mask) { - PS_ASSERT_IMAGE_NON_EMPTY(readout->mask, false); - PS_ASSERT_IMAGES_SIZE_EQUAL(readout->image, readout->mask, false); - PS_ASSERT_IMAGE_TYPE(readout->mask, PS_TYPE_MASK, false); - } - - psImage *mask = NULL; - psImage *image = readout->image; // Image of interest - - if (output) { - mask = output->mask; // Corresponding mask - } - - if (rng) { - psMemIncrRefCounter(rng); - } else { - rng = psRandomAlloc(PS_RANDOM_TAUS, 0); - } - - // XXX note that this now will accept any of several stats options - psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV); - stats->nSubsample = options->sample; - - if (!psImageBackground(stats, NULL, image, mask, options->combine->maskVal, rng)) { - psError(PS_ERR_UNKNOWN, false, "Failure to measure image statistics.\n"); - psFree(stats); - psFree(rng); - return false; - } - if (!isfinite(stats->robustMedian) || !isfinite(stats->robustUQ) || !isfinite(stats->robustLQ)) { - psError(PS_ERR_UNKNOWN, false, "invalide image statistics (nan).\n"); - psFree(stats); - psFree(rng); - return false; - } - - psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "READOUT.MEDIAN", PS_META_REPLACE, "image stats", stats->robustMedian); - psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "READOUT.STDEV", PS_META_REPLACE, "image stats", stats->robustStdev); - + assert(config); + + bool mdok; // Status of MD lookup + int numFiles = psMetadataLookupS32(NULL, config->arguments, "INPUTS.NUM"); // Number of inputs + bool haveMasks = psMetadataLookupBool(&mdok, config->arguments, "INPUTS.MASKS"); // Do we have masks? + bool haveWeights = psMetadataLookupBool(&mdok, config->arguments, "INPUTS.WEIGHTS"); // Do we have weights? + int iter = psMetadataLookupS32(NULL, config->arguments, "ITER"); // Number of rejection iterations + + PS_ASSERT_INT_POSITIVE(iter, false); + + pmFPAview *view = pmFPAviewAlloc(0); // View to component of interest + psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS, 0); // Random number generator + + psMetadata *stats = NULL; // Statistics for output + if (psMetadataLookup(config->arguments, "STATS.NAME")) { + stats = psMetadataAlloc(); + psMetadataAddMetadata(config->arguments, PS_LIST_TAIL, "STATS.DATA", 0, "Statistics output", stats); + } + + psArray *inputs = ppMergeFileDataLevel(config, "PPMERGE.INPUT", PM_FPA_LEVEL_READOUT); // Input images + psFree(inputs); + if (haveMasks) { + psArray *masks = ppMergeFileDataLevel(config, "PPMERGE.INPUT.MASK", PM_FPA_LEVEL_READOUT); + psFree(masks); + } + if (haveWeights) { + psArray *weights = ppMergeFileDataLevel(config, "PPMERGE.INPUT.WEIGHT", PM_FPA_LEVEL_READOUT); + psFree(weights); + } + + if (!ppMergeFileActivate(config, PPMERGE_FILES_INPUT, true)) { + psError(PS_ERR_UNKNOWN, false, "Unable to activate files."); + goto PPMERGE_MASK_ERROR; + } + if (!ppMergeFileActivate(config, PPMERGE_FILES_OUTPUT, true)) { + psError(PS_ERR_UNKNOWN, false, "Unable to activate files."); + goto PPMERGE_MASK_ERROR; + } + + psString outName = ppMergeOutputFile(config); // Name of output file + pmFPAfile *output = psMetadataLookupPtr(NULL, config->files, outName); // Output file + psFree(outName); + assert(output && output->fpa); + pmFPA *outFPA = output->fpa; // Output FPA + + if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) { + goto PPMERGE_MASK_ERROR; + } + pmChip *outChip; // Chip of interest + while ((outChip = pmFPAviewNextChip(view, outFPA, 1))) { + if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) { + goto PPMERGE_MASK_ERROR; + } + + for (int i = 0; i < iter; i++) { + if (!mergeMask(config, view, (i == iter - 1), rng, stats)) { + psError(PS_ERR_UNKNOWN, false, "Unable to merge chip %d", view->chip); + goto PPMERGE_MASK_ERROR; + } + } + + if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) { + goto PPMERGE_MASK_ERROR; + } + } + + psList *fpaList = psListAlloc(NULL);// List of FPAs for concept averaging + for (int i = 0; i < numFiles; i++) { + pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PPMERGE.INPUT", i); // Input file + psListAdd(fpaList, PS_LIST_TAIL, file->fpa); + } + if (!pmConceptsAverageFPAs(outFPA, fpaList)) { + psError(PS_ERR_UNKNOWN, false, "Unable to average FPA concepts."); + psFree(fpaList); + goto PPMERGE_MASK_ERROR; + } + psFree(fpaList); + + if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) { + goto PPMERGE_MASK_ERROR; + } + + psFree(stats); + psFree(view); psFree(rng); + return true; + +PPMERGE_MASK_ERROR: + psFree(stats); + psFree(view); + psFree(rng); + return false; } -bool ppMergeMaskChipStats (const pmChip *chip, - const pmChip *output, - ppMergeOptions *options, - psRandom *rng) -{ - PS_ASSERT_PTR_NON_NULL(chip, false); - - // XXX note that this now will accept any of several stats options - psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV); - - if (rng) { - psMemIncrRefCounter(rng); - } else { - rng = psRandomAlloc(PS_RANDOM_TAUS, 0); - } - - // accumulate a vector of data values using options->sample per readout - psVector *values = psVectorAllocEmpty(options->sample, PS_TYPE_F32); // Vector containing subsample - - for (int nCell = 0; nCell < chip->cells->n; nCell++) { - pmCell *cell = chip->cells->data[nCell]; - if (cell->readouts->n == 0) continue; - - for (int nReadout = 0; nReadout < cell->readouts->n; nReadout++) { - pmReadout *readout = cell->readouts->data[nReadout]; - if (!readout->image) continue; - - psTrace("ppMerge", 4, "Measure statistics for cell %d, readout %d\n", nCell, nReadout); - - psImage *image = readout->image; // Image of interest - - pmCell *cellOutput = output->cells->data[nCell]; - pmReadout *readoutOutput = NULL; - psImage *mask = NULL; - if (cellOutput->readouts->n > 0) { - readoutOutput = cellOutput->readouts->data[0]; - mask = readoutOutput->mask; // Corresponding mask - } - - // Size of image - long nx = image->numCols; - long ny = image->numRows; - const int Npixels = nx*ny; // Total number of pixels - const int Nsubset = PS_MIN(options->sample, Npixels); // Number of pixels in subset - - values = psVectorRealloc (values, values->n + Nsubset); // make sure we have enough space - - // select a subset of the image pixels to measure the stats - for (long i = 0; i < Nsubset; i++) { - double frnd = psRandomUniform(rng); - int pixel = Npixels * frnd; - int ix = pixel % nx; - int iy = pixel / nx; - - if (!isfinite(image->data.F32[iy][ix])) continue; - if (mask && (mask->data.U8[iy][ix] & options->combine->maskVal)) continue; - - float value = image->data.F32[iy][ix]; - values->data.F32[values->n] = value; - values->n ++; - } - } - } - - // no valid data, skip the chip - if (!values->n) { - psFree(values); - psFree(stats); - psFree(rng); - return true; - } - - // calculate the statistics - if (!psVectorStats (stats, values, NULL, NULL, 0)) { - psError(PS_ERR_UNKNOWN, false, "Unable to measure statistics for chip"); - psFree(values); - psFree(stats); - psFree(rng); - return false; - } - if (!isfinite(stats->robustMedian) || !isfinite(stats->robustUQ) || !isfinite(stats->robustLQ)) { - psError(PS_ERR_UNKNOWN, false, "invalid image statistics (nan).\n"); - psFree(values); - psFree(stats); - psFree(rng); - return false; - } - - // supply the stats to the readout analysis metadata - for (int nCell = 0; nCell < chip->cells->n; nCell++) { - pmCell *cell = chip->cells->data[nCell]; - if (cell->readouts->n == 0) continue; - - for (int nReadout = 0; nReadout < cell->readouts->n; nReadout++) { - pmReadout *readout = cell->readouts->data[nReadout]; - if (!readout->image) continue; - - psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "READOUT.MEDIAN", PS_META_REPLACE, "image stats", stats->robustMedian); - psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "READOUT.STDEV", PS_META_REPLACE, "image stats", stats->robustStdev); - } - } - - psLogMsg ("ppMerge", PS_LOG_INFO, "statistics for chip: %f +/- %f\n", stats->robustMedian, stats->robustStdev); - - psFree(values); - psFree(stats); - psFree(rng); - return true; -} Index: /trunk/ppMerge/src/ppMergeScaleZero.c =================================================================== --- /trunk/ppMerge/src/ppMergeScaleZero.c (revision 17226) +++ /trunk/ppMerge/src/ppMergeScaleZero.c (revision 17227) @@ -10,339 +10,207 @@ #include "ppMerge.h" -#include "ppMergeScaleZero.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 - psArray **shutters, // The shutter correction data for each cell - ppMergeData *data,// The data - const ppMergeOptions *options, // The options - const pmConfig *config // The configuration - ) +bool ppMergeScaleZero(pmConfig *config) { - assert(data); - assert(options); assert(config); - if (!options->scale && !options->zero && !options->shutter) { - return true; // We did everything we were asked for + ppMergeType type = psMetadataLookupS32(NULL, config->arguments, "TYPE"); // Type of frame + int numInputs = psMetadataLookupS32(NULL, config->arguments, "INPUTS.NUM"); // Number of inputs + int numCells = psMetadataLookupS32(NULL, config->arguments, "INPUTS.CELLS"); // Number of cells + psStatsOptions meanStat = psMetadataLookupS32(NULL, config->arguments, "MEAN"); // Statistic for mean + psStatsOptions stdevStat = psMetadataLookupS32(NULL, config->arguments, "STDEV"); // Statistic for stdev + psMaskType maskVal = psMetadataLookupU8(NULL, config->arguments, "MASKVAL"); // Value to mask + int shutterSize = psMetadataLookupS32(NULL, config->arguments, "SHUTTER.SIZE"); // Size of shutter region + + psVector *gains = NULL; // Gains for each cell + psArray *shutters = NULL; // Shutter data for each cell + psStats *stats = NULL; // Statistics for background + psImage *background = NULL; // Background measurements per cell per file + + switch (type) { + case PPMERGE_TYPE_BIAS: + case PPMERGE_TYPE_DARK: + // Nothing to measure + return true; + case PPMERGE_TYPE_FLAT: + case PPMERGE_TYPE_FRINGE: + gains = psVectorAlloc(numCells, PS_TYPE_F32); + background = psImageAlloc(numCells, numInputs, PS_TYPE_F32); + psImageInit(background, NAN); + stats = psStatsAlloc(meanStat); + break; + case PPMERGE_TYPE_SHUTTER: + shutters = psArrayAlloc(numCells); + break; + case PPMERGE_TYPE_MASK: + default: + break; } - - assert(config->camera); // Need the camera configuration - assert(config->arguments); // Need the list of files - - 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)); - assert(!options->shutter || shutters); - assert(!shutters || !*shutters || (*shutters)->n == data->numCells); - - // Allocate the outputs - if (options->scale) { - if (*scales) { - psFree(*scales); - } - *scales = psImageAlloc(data->numCells, filenames->n, PS_TYPE_F32); + psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS, 0); // Random number generator + pmFPAview *view = NULL; // View into FPA + + for (int i = 0; i < numInputs; i++) { + pmFPAfileActivate(config->files, false, NULL); + psArray *files = ppMergeFileActivateSingle(config, PPMERGE_FILES_INPUT, true, i); // Activated files + pmFPAfile *input = files->data[0]; // Representative file; should be the image (not mask or weight) + pmFPA *fpa = input->fpa; // FPA of interest + view = pmFPAviewAlloc(0); // View to component of interest + int cellNum = 0; // Index for cell + if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) { + goto ERROR; + } + pmChip *chip; // Chip of interest + while ((chip = pmFPAviewNextChip(view, fpa, 1))) { + if (!chip->file_exists) { + continue; + } + if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) { + goto ERROR; + } + + pmCell *cell; // Cell of interest + while ((cell = pmFPAviewNextCell(view, fpa, 1))) { + if (!cell->file_exists) { + continue; + } + if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) { + goto ERROR; + } + + if (!cell->data_exists) { + continue; + } + + if (cell->readouts->n > 1) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, + "File %d chip %d cell %d contains more than one readout (%ld)", + i, view->chip, view->cell, cell->readouts->n); + goto ERROR; + } + pmReadout *readout = cell->readouts->data[0]; // Readout of interest + + switch (type) { + case PPMERGE_TYPE_FLAT: + case PPMERGE_TYPE_FRINGE: { + // Extract the gain + float gain = psMetadataLookupF32(NULL, cell->concepts, "CELL.GAIN"); // Cell gain + if (!isfinite(gain)) { + psError(PS_ERR_BAD_PARAMETER_VALUE, false, + "CELL.GAIN for file %d chip %d cell %d is not set.", + i, view->chip, view->cell); + goto ERROR; + } + gains->data.F32[cellNum] = gain; + + // Measure the background + if (!psImageBackground(stats, NULL, readout->image, readout->mask, maskVal, rng)) { + psError(PS_ERR_UNKNOWN, false, + "Unable to get statistics for file %d chip %d cell %d", + i, view->chip, view->cell); + goto ERROR; + } + background->data.F32[i][cellNum] = psStatsGetValue(stats, meanStat); + break; + } + case PPMERGE_TYPE_SHUTTER: { + pmShutterCorrectionData *shutter = shutters->data[cellNum]; // Shutter correction data + if (!shutter) { + shutter = pmShutterCorrectionDataAlloc(readout->image->numCols, + readout->image->numRows, + shutterSize); + shutters->data[cellNum] = shutter; + } + if (!pmShutterCorrectionAddReadout(shutter, readout, meanStat, stdevStat, + maskVal, rng)) { + psError(PS_ERR_UNKNOWN, false, + "Can't add file %d chip %d cell %d to shutter correction.", + i, view->chip, view->cell); + goto ERROR; + } + break; + } + case PPMERGE_TYPE_MASK: + default: + psAbort("Should never get here."); + } + + cellNum++; + + if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) { + goto ERROR; + } + } + + if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) { + goto ERROR; + } + } + + if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) { + goto ERROR; + } + + psFree(view); + +#if 0 + // Reset files for reading again + for (int i = 0; i < files->n; i++) { + pmFPAfile *file = files->data[i]; // File of interest + } +#endif + psFree(files); } - if (options->zero) { - if (*zeros) { - psFree(*zeros); - } - *zeros = psImageAlloc(data->numCells, filenames->n, PS_TYPE_F32); + + psFree(rng); rng = NULL; + psFree(stats); stats = NULL; + + // Store results + switch (type) { + case PPMERGE_TYPE_FRINGE: + psMetadataAddImage(config->arguments, PS_LIST_TAIL, "ZEROS", 0, + "Zero to subtract from each input cell", background); + // Flow through + case PPMERGE_TYPE_FLAT: { + // Need to normalize over the focal plane + if (psTraceGetLevel("ppMerge") > 9) { + for (int i = 0; i < gains->n; i++) { + psTrace("ppMerge", 10, "Gain for cell %d is %f\n", i, gains->data.F32[i]); + } + } + psVector *fluxes = NULL; // Solution to fluxes + if (!pmFlatNormalize(&fluxes, &gains, background)) { + psError(PS_ERR_UNKNOWN, false, "Normalisation failed to converge --- continuing anyway."); + psFree(fluxes); + goto ERROR; + } + + psMetadataAddVector(config->arguments, PS_LIST_TAIL, "SCALES", 0, + "Scale to divide into each input file", fluxes); + psFree(fluxes); // Drop reference + break; + } + case PPMERGE_TYPE_SHUTTER: + psMetadataAddArray(config->arguments, PS_LIST_TAIL, "SHUTTER", 0, + "Shutter data", shutters); + break; + case PPMERGE_TYPE_MASK: + default: + psAbort("Should never get here."); } - psRandom *rng = NULL; // Random number generator - if (options->shutter) { - if (*shutters) { - psFree(*shutters); - } - *shutters = psArrayAlloc(data->numCells); - rng = psRandomAlloc(PS_RANDOM_TAUS, 0); - } - - bool fromConcepts = false; // Do we get the scale and zero points from the concepts? - bool first = true; // Are we on the first cell (that sets the standard for the rest)? - bool done = false; // Are we done going through the list? - bool mdok = true; // Status of MD lookup - for (long i = 0; i < data->in->n && !done; i++) { - pmFPA *fpa = data->in->data[i]; // The FPA - if (!fpa) { - continue; - } - long cellNum = -1; // Number of the cell - psArray *chips = fpa->chips; // The array of chips - for (long j = 0; j < chips->n && !done; j++) { - pmChip *chip = chips->data[j]; // The chip - if (!chip) { - continue; - } - psArray *cells = chip->cells; // The array of cells - for (long k = 0; k < cells->n && !done; k++) { - pmCell *cell = cells->data[k]; // The cell - if (!cell) { - continue; - } - cellNum++; - - if (options->scale) { - float scale = psMetadataLookupF32(&mdok, cell->concepts, "PPMERGE.SCALE"); // The scale - if (mdok && !isnan(scale)) { - if (!first && !fromConcepts) { - psWarning("PPMERGE.SCALE and PPMERGE.ZERO have been set " - "for some, but not all cells --- we will re-measure it for all cells."); - done = true; - continue; - } - fromConcepts = true; - (*scales)->data.F32[i][cellNum] = scale; - psTrace("ppMerge", 9, "Scale for input %ld, chip %ld, cell %ld: %f\n", - i, j, k, scale); - } else if (!first && fromConcepts) { - psWarning("PPMERGE.SCALE and PPMERGE.ZERO have been set " - "for some, but not all cells --- we will re-measure it for all cells."); - fromConcepts = false; - done = true; - continue; - } - } - - if (options->zero) { - float zero = psMetadataLookupF32(&mdok, cell->concepts, "PPMERGE.ZERO"); // The zero - if (mdok && !isnan(zero)) { - if (!first && !fromConcepts) { - psWarning("PPMERGE.SCALE and PPMERGE.ZERO have been set " - "for some, but not all cells --- we will re-measure it for all cells."); - done = true; - continue; - } - fromConcepts = true; - (*zeros)->data.F32[i][cellNum] = zero; - psTrace("ppMerge", 9, "Zero for input %ld, chip %ld, cell %ld: %f\n", i, j, k, zero); - } else if (!first && fromConcepts) { - psWarning("PPMERGE.SCALE and PPMERGE.ZERO have been set " - "for some, but not all cells --- we will re-measure it for all cells."); - fromConcepts = false; - done = true; - continue; - } - } - - first = false; - } - } - } - - if (fromConcepts) { - // We've already done everything we need to - psFree(rng); - return true; - } - - psImage *background = psImageAlloc(data->numCells, filenames->n, PS_TYPE_F32); // Background measurements - psImageInit(background, NAN); - psStats *bgStats = psStatsAlloc(options->mean); // Statistic to measure the background - psVector *gains = NULL; // The gains for each cell - if (options->scale) { - gains = psVectorAlloc(data->numCells, PS_TYPE_F32); - } - - pmFPAview *view = pmFPAviewAlloc(0); - - 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("ppMerge", 9, "Opening %s to get background...\n", name); - psFits *inFile = data->files->data[i]; // The FITS file to read - 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)) { - continue; - } - - // Scaling by the background - if (options->scale || options->zero) { - if (!pmCellRead(cell, inFile, config->database)) { - // Nothing here - pmCellFreeData(cell); - continue; - } - - if (cell->readouts->n > 1) { - psWarning("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) { - pmCellFreeData(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])) { - psWarning("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 - psVector *sampleMask = NULL; // Mask for sample - if (readout->mask) { - sampleMask = psVectorAlloc(sampleSize, PS_TYPE_U8); - } - 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]; - } - } - status = psVectorStats(bgStats, sample, sampleMask, NULL, options->combine->maskVal); - if (!status) { - psTrace("ppMerge", 3, "failed to get stats for for %s, cell %d is %f\n", name, cellNum, - background->data.F32[i][cellNum]); - psErrorClear(); - } - psFree(sample); - psFree(sampleMask); - background->data.F32[i][cellNum] = psStatsGetValue(bgStats, options->mean); - psTrace("ppMerge", 3, "Background for %s, cell %d is %f\n", name, cellNum, - background->data.F32[i][cellNum]); - } - - // Shutter correction - if (options->shutter) { - if (!pmCellRead(cell, inFile, config->database)) { - // Nothing here - pmCellFreeData(cell); - continue; - } - - if (cell->readouts->n > 1) { - psWarning("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 - - pmShutterCorrectionData *shutter = (*shutters)->data[cellNum]; // Shutter correction data - if (!shutter) { - shutter = pmShutterCorrectionDataAlloc(readout->image->numCols, - readout->image->numRows, - options->shutterSize); - (*shutters)->data[cellNum] = shutter; - } - if (!pmShutterCorrectionAddReadout(shutter, readout, options->mean, options->stdev, - options->combine->maskVal, rng)) { - psWarning("Can't add file %s chip %d cell %d to shutter correction --- ignored.", - name, j, k); - status = false; - } - } - - - pmCellFreeData(cell); - } - pmChipFreeData(chip); - } - pmFPAFreeData(fpa); - } + + psFree(background); + psFree(shutters); + return true; + +ERROR: + // Common path for errors + psFree(gains); + psFree(background); + psFree(shutters); psFree(rng); - psFree(bgStats); + psFree(stats); psFree(view); - - if (options->scale) { - // Need to normalize over the focal plane - if (psTraceGetLevel("ppMerge") > 9) { - for (int i = 0; i < gains->n; i++) { - psTrace("ppMerge", 10, "Gain for cell %d is %f\n", i, gains->data.F32[i]); - } - } - psVector *fluxes = NULL; // Solution to fluxes - if (!pmFlatNormalize(&fluxes, &gains, background)) { - psWarning("Normalisation failed to converge --- continuing anyway.\n"); - status = false; - } - psFree(gains); - - 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->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("ppMerge") > 9) { - psImage *scalesDeref = *scales; // Dereference the pointer - for (int i = 0; i < scalesDeref->numRows; i++) { - for (int j = 0; j < scalesDeref->numCols; j++) { - psTrace("ppMerge", 9, "Scale for exposure %d, cell %d is: %f\n", i, j, - scalesDeref->data.F32[i][j]); - } - } - } - if (zeros && *zeros && psTraceGetLevel("ppMerge") > 9) { - psImage *zerosDeref = *zeros; // Dereference the pointer - for (int i = 0; i < zerosDeref->numRows; i++) { - for (int j = 0; j < zerosDeref->numCols; j++) { - psTrace("ppMerge", 9, "Zero for exposure %d, cell %d is: %f\n", i, j, - zerosDeref->data.F32[i][j]); - } - } - } - - psFree(background); - - return status; + return false; }