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