Index: /branches/eam_branch_20080719/ppMerge/src/Makefile.am =================================================================== --- /branches/eam_branch_20080719/ppMerge/src/Makefile.am (revision 18759) +++ /branches/eam_branch_20080719/ppMerge/src/Makefile.am (revision 18759) @@ -0,0 +1,29 @@ +bin_PROGRAMS = ppMerge + +ppMerge_CFLAGS = $(PPSTATS_CFLAGS) $(PSMODULE_CFLAGS) $(PSLIB_CFLAGS) +ppMerge_LDFLAGS = $(PPSTATS_LIBS) $(PSMODULE_LIBS) $(PSLIB_LIBS) + +ppMerge_SOURCES = \ + ppMerge.c \ + ppMergeArguments.c \ + ppMergeCamera.c \ + ppMergeFiles.c \ + ppMergeScaleZero.c \ + ppMergeLoop.c \ + ppMergeFileGroup.c \ + ppMergeMask.c + +# ppMergeLoop_Threaded.c +# ppMergeThreadLauncher.c + + +noinst_HEADERS = \ + ppMerge.h + +CLEANFILES = *~ + +clean-local: + -rm -f TAGS +# Tags for emacs +tags: + etags `find . -name \*.[ch] -print` Index: /branches/eam_branch_20080719/ppMerge/src/ppMerge.c =================================================================== --- /branches/eam_branch_20080719/ppMerge/src/ppMerge.c (revision 18759) +++ /branches/eam_branch_20080719/ppMerge/src/ppMerge.c (revision 18759) @@ -0,0 +1,144 @@ +#include "ppMerge.h" + +// Yet to do: +// +// 1. Mask pixels with less than the minimum number of electrons +// 2. Sampling factor for background measurement +// 3. On/off pairs + +int main(int argc, char **argv) +{ + psLibInit(NULL); + psMemSetThreadSafety(false); + psTimerStart(TIMERNAME); + + psExit exitValue = PS_EXIT_SUCCESS; // Exit value for program + + pmConfig *config = pmConfigRead(&argc, argv, PPMERGE_RECIPE); // Configuration + if (!config) { + 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, 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); + if (!options) { + psErrorStackPrint(stderr, "Unable to parse options."); + exit(EXIT_FAILURE); + } + + // Check the inputs + ppMergeData *data = ppMergeCheckInputs(options, config); + if (!data) { + psErrorStackPrint(stderr, "Not enough valid input files."); + exit(EXIT_FAILURE); + } + + if (options->mask) { + // Generate a mask + ppMergeMask(data, options, config); + } else { + + psImage *scale = NULL; // The scalings + psImage *zero = NULL; // The zeros + psArray *shutters = NULL; // The shutter correction data + + // Measure the background in each image + ppMergeScaleZero(&scale, &zero, &shutters, data, options, config); + + // Do the combination and write + ppMergeCombine(scale, zero, shutters, data, options, config); + + psFree(scale); + psFree(zero); + psFree(shutters); + } + + // Output the statistics + if (data->statsFile && data->stats) { + psString statsOut = psMetadataConfigFormat(data->stats); // String to write out + fprintf(data->statsFile, "%s", statsOut); + psFree(statsOut); + } + + // Cleaning up + psFree(data); + psFree(options); +#endif + + die: + psTrace("ppSub", 1, "Finished at %f sec\n", psTimerMark("ppSub")); + psTimerStop(); + + psFree(config); + pmConfigDone(); + psLibFinalize(); + + exit(exitValue); +} Index: /branches/eam_branch_20080719/ppMerge/src/ppMerge.h =================================================================== --- /branches/eam_branch_20080719/ppMerge/src/ppMerge.h (revision 18759) +++ /branches/eam_branch_20080719/ppMerge/src/ppMerge.h (revision 18759) @@ -0,0 +1,105 @@ +#ifndef PP_MERGE_H +#define PP_MERGE_H + +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include +#include +#include +#include +#include + +#define TIMERNAME "ppMerge" // Name for timer +#define PPMERGE_RECIPE "PPMERGE" // Recipe name +#define THREADED 0 + +// 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; + +typedef struct { + psArray *readouts; + bool read; + bool busy; +} ppMergeFileGroup; + +// 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(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 + ); + + +ppMergeFileGroup *ppMergeFileGroupAlloc(); +ppMergeFileGroup *ppMergeReadChunk (psArray *fileGroups, pmConfig *config, int numChunk); +void *ppMergeThreadLauncher (void *data); + +#endif Index: /branches/eam_branch_20080719/ppMerge/src/ppMergeArguments.c =================================================================== --- /branches/eam_branch_20080719/ppMerge/src/ppMergeArguments.c (revision 18759) +++ /branches/eam_branch_20080719/ppMerge/src/ppMergeArguments.c (revision 18759) @@ -0,0 +1,388 @@ +#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); + + // XXX EAM : not clear this should be allowed on the command line. + // 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-set-value",0, "Value to set for output mask pixels", NULL); + psMetadataAddBool(arguments, PS_LIST_TAIL, "-mask-chip", 0, "Measure mask statistics by chip?", false); + + psMetadataAddBool(arguments, PS_LIST_TAIL, "-mask-smooth-suspect", 0, "smooth suspect pixels?", false); + psMetadataAddF32(arguments, PS_LIST_TAIL, "-mask-smooth-scale", 0, "smoothing sigma for suspect pixels", NAN); + + // chip & cell selections are used to limit chips to be processed + int argnum; + if ((argnum = psArgumentGet (argc, argv, "-chip"))) { + psArgumentRemove (argnum, &argc, argv); + psMetadataAddStr (config->arguments, PS_LIST_TAIL, "CHIP_SELECTIONS", PS_DATA_STRING, "", argv[argnum]); + psArgumentRemove (argnum, &argc, argv); + } + if ((argnum = psArgumentGet (argc, argv, "-cell"))) { + psArgumentRemove (argnum, &argc, argv); + psMetadataAddStr (config->arguments, PS_LIST_TAIL, "CELL_SELECTIONS", PS_DATA_STRING, "", argv[argnum]); + psArgumentRemove (argnum, &argc, argv); + } + +# if (THREADED) + // Number of threads + if ((argnum = psArgumentGet(argc, argv, "-threads"))) { + psArgumentRemove(argnum, &argc, argv); + int nThreads = atoi(argv[argnum]); + psMetadataAddS32(config->arguments, PS_LIST_TAIL, "NTHREADS", 0, "number of warp threads", nThreads); + psArgumentRemove(argnum, &argc, argv); + + // create the thread pool with number of desired threads, supplying our thread launcher function + // XXX need to determine the number of threads from the config data + psThreadPoolInit (nThreads, &ppMergeThreadLauncher); + } +# endif + + 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; + goto VALID; + } + if (strcasecmp(typeStr, "DARK") == 0 || + strcasecmp(typeStr, "DARK_PREMASK") == 0) { + type = PPMERGE_TYPE_DARK; + goto VALID; + } + if (strcasecmp(typeStr, "FLAT") == 0 || + strcasecmp(typeStr, "DOMEFLAT") == 0 || + strcasecmp(typeStr, "SKYFLAT") == 0 || + strcasecmp(typeStr, "FLAT_PREMASK") == 0 || + strcasecmp(typeStr, "SKYFLAT_PREMASK") == 0 || + strcasecmp(typeStr, "DOMEFLAT_PREMASK") == 0) { + type = PPMERGE_TYPE_FLAT; + goto VALID; + } + if (strcasecmp(typeStr, "FRINGE") == 0) { + type = PPMERGE_TYPE_FRINGE; + goto VALID; + } + if (strcasecmp(typeStr, "SHUTTER") == 0) { + type = PPMERGE_TYPE_SHUTTER; + goto VALID; + } + if (strcasecmp(typeStr, "MASK") == 0 || + strcasecmp(typeStr, "DARKMASK") == 0 || + strcasecmp(typeStr, "FLATMASK") == 0) { + type = PPMERGE_TYPE_MASK; + goto VALID; + } + psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Unrecognised image type: %s", typeStr); + goto ERROR; + + VALID: + 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"); + + // XXX we do not supply this on the command line + // VALUE_ARG_RECIPE_MASK("-maskval", "MASKVAL"); + + VALUE_ARG_RECIPE_STAT("-combine", "COMBINE"); + VALUE_ARG_RECIPE_STAT("-mean", "MEAN"); + VALUE_ARG_RECIPE_STAT("-stdev", "STDEV"); + + // Fringe construction parameters + VALUE_ARG_RECIPE_INT("-fringe-num", "FRINGE.NUM", S32, 0); + VALUE_ARG_RECIPE_INT("-fringe-size", "FRINGE.SIZE", S32, 0); + VALUE_ARG_RECIPE_INT("-fringe-xsmooth", "FRINGE.XSMOOTH", S32, 0); + VALUE_ARG_RECIPE_INT("-fringe-ysmooth", "FRINGE.YSMOOTH", S32, 0); + + // Shutter construction parameters + VALUE_ARG_RECIPE_INT("-shutter-size", "SHUTTER.SIZE", S32, 0); + + // Mask construction parameters + VALUE_ARG_RECIPE_FLOAT("-mask-suspect", "MASK.SUSPECT", F32); + VALUE_ARG_RECIPE_FLOAT("-mask-bad", "MASK.BAD", F32); + VALUE_ARG_RECIPE_INT("-mask-grow", "MASK.GROW", S32, 0); + VALUE_ARG_RECIPE_STR("-mask-set-value", "MASK.SET.VALUE"); + VALUE_ARG_RECIPE_BOOL("-mask-chip", "MASK.CHIPSTATS"); + + VALUE_ARG_RECIPE_BOOL("-mask-smooth-suspect", "MASK.SMOOTH.SUSPECT"); + VALUE_ARG_RECIPE_FLOAT("-mask-smooth-scale", "MASK.SMOOTH.SCALE", F32); + + 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/eam_branch_20080719/ppMerge/src/ppMergeCamera.c =================================================================== --- /branches/eam_branch_20080719/ppMerge/src/ppMergeCamera.c (revision 18759) +++ /branches/eam_branch_20080719/ppMerge/src/ppMergeCamera.c (revision 18759) @@ -0,0 +1,353 @@ +#include "ppMerge.h" + +// Define an output file, with its own FPA +bool outputFile(pmConfig *config, // Configuration + const char *name, // Name of output file + pmFPAfileType type, // Type of file + const char *description, // Description of file + psMetadata *format, // Camera format + pmFPAview *view // View for PHU + ) +{ + assert(config); + assert(name && strlen(name) > 0); + assert(view); + + // Output image + pmFPA *fpa = pmFPAConstruct(config->camera, config->cameraName); // FPA to contain the output + if (!fpa) { + psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to construct an FPA from camera configuration."); + return false; + } + + // Chip selection: turn on only the chips specified. + bool status; + char *chipLine = psMetadataLookupStr(&status, config->arguments, "CHIP_SELECTIONS"); + psArray *chips = psStringSplitArray (chipLine, ",", false); + if (chips->n > 0) { + pmFPASelectChip (fpa, -1, true); // deselect all chips + for (int i = 0; i < chips->n; i++) { + int chipNum = atoi(chips->data[i]); + if (! pmFPASelectChip(fpa, chipNum, false)) { + psError(PS_ERR_IO, false, "Chip number %d doesn't exist in camera.\n", chipNum); + psFree (chips); + return false; + } + } + } + psFree (chips); + + // Cell selection: turn on only the cells specified. Note that this affects the same cell + // for all chips + char *cellLine = psMetadataLookupStr(&status, config->arguments, "CELL_SELECTIONS"); + psArray *cells = psStringSplitArray (cellLine, ",", false); + if (cells->n > 0) { + for (int i = 0; i < fpa->chips->n; i++) { + pmChip *chip = fpa->chips->data[i]; + pmChipSelectCell (chip, -1, true); // deselect all cells + for (int j = 0; j < cells->n; j++) { + int cellNum = atoi(cells->data[j]); + if (! pmChipSelectCell(chip, cellNum, false)) { + psError(PS_ERR_IO, false, "Cell number %d doesn't exist in camera.\n", cellNum); + psFree (cells); + return false; + } + } + } + } + psFree (cells); + + pmFPAfile *output = pmFPAfileDefineOutput(config, fpa, name); + psFree(fpa); // Drop reference + if (!output) { + psError(PS_ERR_IO, false, "Unable to generate output file from %s", name); + return false; + } + if (output->type != type) { + psError(PS_ERR_IO, true, "%s is not of type %s", name, pmFPAfileStringFromType(type)); + return false; + } + output->save = true; + + if (!pmFPAAddSourceFromView(fpa, description, view, format)) { + psError(PS_ERR_UNKNOWN, false, "Unable to generate output FPA."); + return false; + } + + return true; +} + + +bool ppMergeCamera(pmConfig *config) +{ + bool haveMasks = false; // Do we have masks? + bool haveWeights = false; // Do we have weight maps? + + ppMergeType type = psMetadataLookupS32(NULL, config->arguments, "TYPE"); // Type of frame + psMetadata *inputs = psMetadataLookupMetadata(NULL, config->arguments, "INPUTS"); // The inputs info + psMetadataIterator *iter = psMetadataIteratorAlloc(inputs, PS_LIST_HEAD, NULL); // Iterator + psMetadataItem *item; // Item from iteration + int numFiles = 0; // Number of files + while ((item = psMetadataGetAndIncrement(iter))) { + if (item->type != PS_DATA_METADATA) { + psError(PS_ERR_BAD_PARAMETER_TYPE, true, + "Component %s of the input metadata is not of type METADATA", item->name); + psFree(iter); + return false; + } + + psMetadata *input = item->data.md; // The input metadata of interest + + psString image = psMetadataLookupStr(NULL, input, "IMAGE"); // Name of image + if (!image || strlen(image) == 0) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Component %s lacks IMAGE of type STR", item->name); + psFree(iter); + return false; + } + + bool mdok; + psString mask = psMetadataLookupStr(&mdok, input, "MASK"); // Name of mask + psString weight = psMetadataLookupStr(&mdok, input, "WEIGHT"); // Name of weight map + + // Add the image file + psArray *imageFiles = psArrayAlloc(1); // Array of filenames for this FPA + imageFiles->data[0] = psMemIncrRefCounter(image); + psMetadataAddArray(config->arguments, PS_LIST_TAIL, "IMAGE.FILENAMES", PS_META_REPLACE, + "Filenames of image files", imageFiles); + psFree(imageFiles); + + bool found = false; // Found the file? + pmFPAfile *imageFile = pmFPAfileDefineFromArgs(&found, config, "PPMERGE.INPUT", "IMAGE.FILENAMES"); + if (!imageFile || !found) { + psError(PS_ERR_UNKNOWN, false, "Unable to define file from image %d (%s)", numFiles, image); + return false; + } + if (imageFile->type != PM_FPA_FILE_IMAGE) { + psError(PS_ERR_IO, true, "PPMERGE.INPUT is not of type IMAGE"); + return false; + } + + // Optionally add the mask file + if (mask && strlen(mask) > 0) { + psArray *maskFiles = psArrayAlloc(1); // Array of filenames for this FPA + maskFiles->data[0] = psMemIncrRefCounter(mask); + psMetadataAddArray(config->arguments, PS_LIST_TAIL, "MASK.FILENAMES", PS_META_REPLACE, + "Filenames of mask files", maskFiles); + psFree(maskFiles); + + bool status; + pmFPAfile *maskFile = pmFPAfileBindFromArgs(&status, imageFile, config, "PPMERGE.INPUT.MASK", + "MASK.FILENAMES"); + if (!status) { + psError(PS_ERR_UNKNOWN, false, "Unable to define file from mask %d (%s)", numFiles, mask); + return false; + } + if (maskFile->type != PM_FPA_FILE_MASK) { + psError(PS_ERR_IO, true, "PPMERGE.INPUT.MASK is not of type MASK"); + return false; + } + haveMasks = true; + } + + // Optionally add the weight file + if (weight && strlen(weight) > 0) { + haveWeights = true; + psArray *weightFiles = psArrayAlloc(1); // Array of filenames for this FPA + weightFiles->data[0] = psMemIncrRefCounter(weight); + psMetadataAddArray(config->arguments, PS_LIST_TAIL, "WEIGHT.FILENAMES", PS_META_REPLACE, + "Filenames of weight files", weightFiles); + psFree(weightFiles); + + bool status; + pmFPAfile *weightFile = pmFPAfileBindFromArgs(&status, imageFile, config, "PPMERGE.INPUT.WEIGHT", + "WEIGHT.FILENAMES"); + if (!status) { + psError(PS_ERR_UNKNOWN, false, "Unable to define file from weight %d (%s)", numFiles, weight); + return false; + } + if (weightFile->type != PM_FPA_FILE_WEIGHT) { + psError(PS_ERR_IO, true, "PPMERGE.INPUT.WEIGHT is not of type WEIGHT"); + return false; + } + haveWeights = true; + } + + numFiles++; + } + psFree(iter); + psMetadataRemoveKey(config->arguments, "IMAGE.FILENAMES"); + if (psMetadataLookup(config->arguments, "MASK.FILENAMES")) { + psMetadataRemoveKey(config->arguments, "MASK.FILENAMES"); + } + if (psMetadataLookup(config->arguments, "WEIGHT.FILENAMES")) { + psMetadataRemoveKey(config->arguments, "WEIGHT.FILENAMES"); + } + if (psMetadataLookup(config->arguments, "PSF.FILENAMES")) { + psMetadataRemoveKey(config->arguments, "PSF.FILENAMES"); + } + + psMetadataAddS32(config->arguments, PS_LIST_TAIL, "INPUTS.NUM", 0, "Number of input files", numFiles); + psMetadataAddBool(config->arguments, PS_LIST_TAIL, "INPUTS.MASKS", 0, "Got input masks?", haveMasks); + psMetadataAddBool(config->arguments, PS_LIST_TAIL, "INPUTS.WEIGHTS", 0, + "Got input weights?", haveWeights); + + + // Check that all the inputs are consistent + +// Check an FPA level to ensure the camera format and PHU are consistent across all input files +#define CHECK_LEVEL(HDU, CHIP, CELL) { \ + if (HDU) { \ + if (!phuView) { \ + phuView = pmFPAviewAlloc(0); \ + phuView->chip = CHIP; \ + phuView->cell = CELL; \ + } else if ((phuView->chip != (CHIP)) && (phuView->cell != (CELL))) { \ + psError(PS_ERR_UNKNOWN, true, "Differing PHU for input %d", i); \ + psFree(phuView); \ + return false; \ + } \ + if (!format) { \ + format = (HDU)->format; \ + } else if (format != (HDU)->format) { \ + psError(PS_ERR_UNKNOWN, true, "Camera format %d doesn't match: %p vs %p", \ + i, format, (HDU)->format); \ + psFree(phuView); \ + return false; \ + } \ + continue; \ + } \ +} + + psMetadata *format = NULL; // Camera format + pmFPAview *phuView = NULL; // View to PHU + for (int i = 0; i < numFiles; i++) { + pmFPAfile *input = pmFPAfileSelectSingle(config->files, "PPMERGE.INPUT", i); // File of interest + pmFPA *fpa = input->fpa; // FPA of interest + CHECK_LEVEL(fpa->hdu, -1, -1); + psArray *chips = fpa->chips; // Array of chips + for (int j = 0; j < chips->n; j++) { + pmChip *chip = chips->data[j]; // Chip of interest + CHECK_LEVEL(chip->hdu, j, -1); + psArray *cells = chip->cells; // Array of cells + for (int k = 0; k < cells->n; k++) { + pmCell *cell = cells->data[k]; // Cell of interest + CHECK_LEVEL(cell->hdu, j, k); + } + } + } + if (!phuView || !format) { + psError(PS_ERR_UNEXPECTED_NULL, true, "Unable to find PHU for input files."); + psFree(phuView); + return false; + } + + // Cull chips and cells that don't have data + // Otherwise the abundance of metadata in the concepts (esp. for GPC) can overload the memory + for (int i = 0; i < numFiles; i++) { + pmFPAfile *input = pmFPAfileSelectSingle(config->files, "PPMERGE.INPUT", i); // File of interest + pmFPA *fpa = input->fpa; // FPA of interest + psArray *chips = fpa->chips; // Array of chips in output + for (int i = 0; i < chips->n; i++) { + pmChip *chip = chips->data[i]; // Chip of interest + psArray *cells = chip->cells; // Array of cells + int culled = 0; // Number of culled cells + for (int j = 0; j < cells->n; j++) { + pmCell *cell = cells->data[j]; + pmHDU *hdu = pmHDUGetLowest(fpa, chip, cell); // HDU for cell + if (!hdu || hdu->blankPHU) { + cell->data_exists = false; + cell->file_exists = false; + culled++; + if (cell->concepts) { + psFree(cell->concepts); + cell->concepts = NULL; + } + } + } + if (culled == cells->n) { + chip->data_exists = false; + chip->file_exists = false; + if (chip->concepts) { + psFree(chip->concepts); + chip->concepts = NULL; + } + } + } + } + + // 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, config->cameraName); // FPA to contain the output + if (!fpa) { + psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to construct an FPA from camera configuration."); + psFree(phuView); + return false; + } + + psString outName = ppMergeOutputFile(config); // Name of output file + + pmFPAfileType fileType = PM_FPA_FILE_NONE; // Type of output file + switch (type) { + case PPMERGE_TYPE_BIAS: + fileType = PM_FPA_FILE_IMAGE; + break; + case PPMERGE_TYPE_DARK: + fileType = PM_FPA_FILE_DARK; + break; + case PPMERGE_TYPE_MASK: + fileType = PM_FPA_FILE_MASK; + break; + case PPMERGE_TYPE_SHUTTER: + fileType = PM_FPA_FILE_IMAGE; + break; + case PPMERGE_TYPE_FLAT: + fileType = PM_FPA_FILE_IMAGE; + break; + case PPMERGE_TYPE_FRINGE: + fileType = PM_FPA_FILE_FRINGE; + default: + psAbort("Unknown frame type: %x", type); + } + + if (!outputFile(config, outName, fileType, "Merged detrend", format, phuView)) { + psFree(outName); + psFree(phuView); + return false; + } + psFree(outName); + + if (!outputFile(config, "PPMERGE.OUTPUT.SIGMA", PM_FPA_FILE_IMAGE, "Merge sigma", format, phuView)) { + psFree(phuView); + return false; + } + + if (!outputFile(config, "PPMERGE.OUTPUT.COUNT", PM_FPA_FILE_IMAGE, "Merged count", format, phuView)) { + psFree(phuView); + return false; + } + + psFree(phuView); + + return true; +} Index: /branches/eam_branch_20080719/ppMerge/src/ppMergeFileGroup.c =================================================================== --- /branches/eam_branch_20080719/ppMerge/src/ppMergeFileGroup.c (revision 18759) +++ /branches/eam_branch_20080719/ppMerge/src/ppMergeFileGroup.c (revision 18759) @@ -0,0 +1,20 @@ +# include "ppMerge.h" + +void ppMergeFileGroupFree (ppMergeFileGroup *group) { + + if (!group) return; + + psFree (group->readouts); + return; +} + +ppMergeFileGroup *ppMergeFileGroupAlloc() { + + ppMergeFileGroup *group = (ppMergeFileGroup *) psAlloc(sizeof(ppMergeFileGroup)); + psMemSetDeallocator(group, (psFreeFunc) ppMergeFileGroupFree); + + group->readouts = NULL; + group->read = false; + group->busy = false; + return group; +} Index: /branches/eam_branch_20080719/ppMerge/src/ppMergeFiles.c =================================================================== --- /branches/eam_branch_20080719/ppMerge/src/ppMergeFiles.c (revision 18759) +++ /branches/eam_branch_20080719/ppMerge/src/ppMergeFiles.c (revision 18759) @@ -0,0 +1,205 @@ +#include "ppMerge.h" + +const char *allFiles[] = { "PPMERGE.INPUT", "PPMERGE.INPUT.MASK", "PPMERGE.INPUT.WEIGHT", + "PPMERGE.OUTPUT", "PPMERGE.OUTPUT.COUNT", "PPMERGE.OUTPUT.SIGMA", + NULL }; // All files +const char *inputFiles[] = { "PPMERGE.INPUT", "PPMERGE.INPUT.MASK", "PPMERGE.INPUT.WEIGHT", + NULL }; // Input files +const char *outputFiles[] = { "PPMERGE.OUTPUT", "PPMERGE.OUTPUT.COUNT", "PPMERGE.OUTPUT.SIGMA", + NULL }; // Output files + +// Select file list based on enum +static const char **selectFiles(ppMergeFiles files) +{ + switch (files) { + case PPMERGE_FILES_ALL: return allFiles; + case PPMERGE_FILES_INPUT: return inputFiles; + case PPMERGE_FILES_OUTPUT: return outputFiles; + default: + psAbort("Invalid file option"); + } + return NULL; +} + +bool ppMergeFileReadInput(pmConfig *config, pmReadout *readout, int num, int rows) +{ + pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PPMERGE.INPUT", num); + if (!pmReadoutReadChunk(readout, file->fits, 0, rows, 0, config)) { + 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, config)) { + 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, config)) { + 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/eam_branch_20080719/ppMerge/src/ppMergeLoop.c =================================================================== --- /branches/eam_branch_20080719/ppMerge/src/ppMergeLoop.c (revision 18759) +++ /branches/eam_branch_20080719/ppMerge/src/ppMergeLoop.c (revision 18759) @@ -0,0 +1,358 @@ +#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 + 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 + + // set the mask and mark bit values based on the named masks + psMaskType maskVal; + psMaskType markVal; + if (!pmConfigMaskSetBits (&maskVal, &markVal, config)) { + psError (PS_ERR_UNKNOWN, true, "Unable to define the mask bit values"); + return false; + } + + pmCombineParams *combination = pmCombineParamsAlloc(combineStat); // Combination parameters + combination->maskVal = maskVal; + combination->blank = pmConfigMaskGet("BLANK", config); + combination->nKeep = nKeep; + combination->fracHigh = frachigh; + combination->fracLow = fraclow; + combination->iter = iter; + combination->rej = rej; + combination->weights = useWeights; + + 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); + } + + pmFPAview *view = pmFPAviewAlloc(0); // View to component of interest + + // 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 + + + 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 + int cellNum = 0; // Index of cell + if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) { + goto ERROR; + } + pmChip *outChip; // Chip of interest + while ((outChip = pmFPAviewNextChip(view, outFPA, 1))) { + if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) { + goto ERROR; + } + pmCell *outCell; // Cell of interest + while ((outCell = pmFPAviewNextCell(view, outFPA, 1))) { + if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) { + goto ERROR; + } + + pmHDU *hdu = pmHDUGetLowest(outFPA, outChip, outCell); // HDU for cell + if (!hdu || hdu->blankPHU) { + // No data here + continue; + } + + pmReadout *outRO = pmReadoutAlloc(outCell); + + psArray *readouts = psArrayAlloc(numFiles); // Input readouts + for (int i = 0; i < numFiles; i++) { + // We need to do some of the opening ourselves + if (!ppMergeFileOpenInput(config, view, i)) { + psError(PS_ERR_UNKNOWN, false, "Unable to open file %d", i); + goto ERROR; + } + + pmFPAfile *input = pmFPAfileSelectSingle(config->files, "PPMERGE.INPUT", i); + pmCell *inCell = pmFPAviewThisCell(view, input->fpa); // Input cell + readouts->data[i] = pmReadoutAlloc(inCell); + } + + float shutterRef = NAN; // Reference shutter correction + if (type == PPMERGE_TYPE_SHUTTER) { + shutterRef = pmShutterCorrectionReference(shutters->data[cellNum]); + } + + // Read convolutions by chunks + bool more = true; // More to read? + for (int numChunk = 0; more; numChunk++) { + psTrace("ppStack", 2, "Initial stack of chunk %d....\n", numChunk); + for (int i = 0; i < numFiles; i++) { + pmReadout *inRO = readouts->data[i]; // Input readout + + // Read a chunk from a file + #define READ_CHUNK(NAME,TYPE) { \ + pmFPAfile *file = pmFPAfileSelectSingle(config->files, NAME, i); \ + if (!pmReadoutReadChunk##TYPE(inRO, file->fits, 0, rows, 0, config)) { \ + psError(PS_ERR_IO, false, "Unable to read chunk %d for file %s %d", \ + numChunk, NAME, i); \ + psFree(readouts); \ + psFree(outRO); \ + goto ERROR; \ + } \ + } + + READ_CHUNK("PPMERGE.INPUT", /* Blank */); + if (haveMasks) { + READ_CHUNK("PPMERGE.INPUT.MASK", Mask); + } + if (haveWeights) { + READ_CHUNK("PPMERGE.INPUT.WEIGHT", Weight); + } + } + + switch (type) { + case PPMERGE_TYPE_SHUTTER: + if (!pmShutterCorrectionGenerate(outRO, NULL, readouts, shutterRef, + shutters->data[cellNum], iter, rej, maskVal)) { + psFree(readouts); + psFree(outRO); + goto ERROR; + } + break; + case PPMERGE_TYPE_DARK: + if (!pmDarkCombine(outCell, readouts, darkOrdinates, darkNorm, iter, rej, maskVal)) { + psFree(readouts); + psFree(outRO); + goto ERROR; + } + break; + case PPMERGE_TYPE_BIAS: + case PPMERGE_TYPE_FLAT: + case PPMERGE_TYPE_FRINGE: + if (!pmReadoutCombine(outRO, readouts, zeros, scales, combination)) { + psFree(readouts); + psFree(outRO); + goto ERROR; + } + break; + default: + psAbort("Should never get here."); + } + + + for (int i = 0; i < numFiles && more; i++) { + pmReadout *inRO = readouts->data[i]; + + // Check to see if there's more chunks to read + #define MORE_CHUNK(NAME,TYPE) { \ + pmFPAfile *file = pmFPAfileSelectSingle(config->files, NAME, i); \ + more &= pmReadoutMore##TYPE(inRO, file->fits, 0, rows, config); \ + } + + MORE_CHUNK("PPMERGE.INPUT", /* Blank */); + if (haveMasks) { + MORE_CHUNK("PPMERGE.INPUT.MASK", Mask); + } + if (haveWeights) { + MORE_CHUNK("PPMERGE.INPUT.WEIGHT", Weight); + } + } + } + + // Get list of cells for concepts averaging + psList *inCells = psListAlloc(NULL); // List of cells + for (int i = 0; i < numFiles; i++) { + pmReadout *readout = readouts->data[i]; // Readout of interest + psListAdd(inCells, PS_LIST_TAIL, readout->parent); + } + if (!pmConceptsAverageCells(outCell, inCells, NULL, NULL, true)) { + psError(PS_ERR_UNKNOWN, false, "Unable to average cell concepts."); + psFree(inCells); + psFree(outRO); + goto ERROR; + } + psFree(inCells); + + psFree(readouts); + + // Plug supplementary images into their own FPAs + { + pmCell *countsCell = pmFPAfileThisCell(config->files, view, "PPMERGE.OUTPUT.COUNT"); + pmReadout *countsRO = pmReadoutAlloc(countsCell); // Readout with count of inputs per pixel + psImage *counts = psMetadataLookupPtr(NULL, outRO->analysis, PM_READOUT_STACK_ANALYSIS_COUNT); + countsRO->image = psImageCopy(countsRO->image, counts, PS_TYPE_F32); + psMetadataRemoveKey(outRO->analysis, PM_READOUT_STACK_ANALYSIS_COUNT); + countsRO->data_exists = countsCell->data_exists = countsCell->parent->data_exists = true; + psFree(countsRO); + + pmCell *sigmaCell = pmFPAfileThisCell(config->files, view, "PPMERGE.OUTPUT.SIGMA"); + pmReadout *sigmaRO = pmReadoutAlloc(sigmaCell); // Readout with stdev per pixel + psImage *sigma = psMetadataLookupPtr(NULL, outRO->analysis, PM_READOUT_STACK_ANALYSIS_SIGMA); + sigmaRO->image = psImageCopy(sigmaRO->image, sigma, PS_TYPE_F32); + psMetadataRemoveKey(outRO->analysis, PM_READOUT_STACK_ANALYSIS_SIGMA); + sigmaRO->data_exists = sigmaCell->data_exists = sigmaCell->parent->data_exists = true; + psFree(sigmaRO); + } + + // Measure the fringes for this cell + // + // XXX Need to deal with multiple components: we will do this by building up the components + // Read the existing fringe measurements + // Use existing regions to measure fringe statistics + // Add the new fringe measurements to the existing fringe measurements. + // Write the appended fringe measurements. + // Read in the "output" file to get the existing components. + // Put the new readout into the cell after the existing readouts. + if (type == PPMERGE_TYPE_FRINGE && outRO) { + pmFringeRegions *regions = pmFringeRegionsAlloc(fringeNum, fringeSize, fringeSize, + fringeSmoothX, fringeSmoothY); + pmFringeStats *fringe = pmFringeStatsMeasure(regions, outRO, maskVal); + psFree(regions); + if (!fringe) { + psError(PS_ERR_UNKNOWN, false, "Unable to measure fringe statistics.\n"); + psFree(outRO); + goto ERROR; + } + + psArray *fringes = psArrayAlloc(1); // Array of fringes + fringes->data[0] = fringe; + + pmFringesFormat(outCell, NULL, fringes); + psFree(fringes); // Drop reference + } + + if (!ppStatsFPA(stats, outFPA, view, maskVal, 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 (outChip->data_exists) { + psList *inChips = psListAlloc(NULL); + for (int i=0; i < numFiles; i++) { + pmChip *chip = pmFPAviewThisChip(view, ((pmFPAfile *)inputs->data[i])->fpa); + psListAdd(inChips, PS_LIST_TAIL, chip); + } + if (!pmConceptsAverageChips(outChip, inChips, true)) { + psError(PS_ERR_UNKNOWN, false, "Unable to average Chip concepts."); + psFree(inChips); + goto ERROR; + } + psFree(inChips); + } + if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) { + goto ERROR; + } + } + + // Get list of FPAs for concepts averaging + psList *inFPAs = psListAlloc(NULL); // List of FPAs + for (int i = 0; i < numFiles; i++) { + pmFPAfile *input = inputs->data[i]; // Input file + psListAdd(inFPAs, PS_LIST_TAIL, input->fpa); + } + if (!pmConceptsAverageFPAs(output->fpa, inFPAs)) { + psError(PS_ERR_UNKNOWN, false, "Unable to average FPA concepts."); + psFree(inFPAs); + goto ERROR; + } + psFree(inFPAs); + + + if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) { + goto ERROR; + } + + psFree(view); + psFree(combination); + psFree(inputs); + psFree(masks); + psFree(weights); + psFree(stats); + return true; + +ERROR: + psFree(view); + psFree(combination); + psFree(inputs); + psFree(masks); + psFree(weights); + psFree(stats); + return false; +} + Index: /branches/eam_branch_20080719/ppMerge/src/ppMergeLoop_Threaded.c =================================================================== --- /branches/eam_branch_20080719/ppMerge/src/ppMergeLoop_Threaded.c (revision 18759) +++ /branches/eam_branch_20080719/ppMerge/src/ppMergeLoop_Threaded.c (revision 18759) @@ -0,0 +1,403 @@ +#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); + } + + int nThreads = psMetadataLookupS32 (&mdok, arguments, "NTHREADS"); + if (!mdok) nThreads = 0; + + // 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 + 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 + + // set the mask and mark bit values based on the named masks + psMaskType maskVal; + psMaskType markVal; + if (!pmConfigMaskSetBits (&maskVal, &markVal, config)) { + psError (PS_ERR_UNKNOWN, true, "Unable to define the mask bit values"); + return false; + } + + pmCombineParams *combination = pmCombineParamsAlloc(combineStat); // Combination parameters + combination->maskVal = maskVal; + combination->blank = pmConfigMaskGet("BLANK", config); + combination->nKeep = nKeep; + combination->fracHigh = frachigh; + combination->fracLow = fraclow; + combination->iter = iter; + combination->rej = rej; + combination->weights = useWeights; + + 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); + } + + pmFPAview *view = pmFPAviewAlloc(0); // View to component of interest + + // 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 + + + 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 + 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; + } + + float shutterRef = NAN; // Reference shutter correction + if (type == PPMERGE_TYPE_SHUTTER) { + shutterRef = pmShutterCorrectionReference(shutters->data[cellNum]); + } + + pmReadout *outRO = pmReadoutAlloc(outCell); + + // open the input files (we need to do the work ourselves) + for (int i = 0; i < numFiles; i++) { + if (!ppMergeFileOpenInput(config, view, i)) { + psError(PS_ERR_UNKNOWN, false, "Unable to open file %d", i); + goto ERROR; + } + } + + psArray *fileGroups = psArrayAlloc (nThreads + 1); + + // generate readouts for each input file in each file group + for (int i = 0; i < fileGroups->n; i++) { + psArray *readouts = psArrayAlloc(numFiles); // Input readouts + for (int j = 0; j < numFiles; j++) { + pmFPAfile *input = pmFPAfileSelectSingle(config->files, "PPMERGE.INPUT", i); + pmCell *inCell = pmFPAviewThisCell(view, input->fpa); // Input cell + readouts->data[j] = pmReadoutAlloc(inCell); + } + + ppMergeFileGroup *fileGroup = ppMergeFileGroupAlloc(); + fileGroup->readouts = readouts; + fileGroup->read = false; + fileGroup->busy = false; + fileGroups->data[i] = fileGroup; + } + + // Read input data by chunks + for (int numChunk = 0; true; numChunk++) { + + ppMergeFileGroup *fileGroup = ppMergeReadChunk (&status, fileGroups, config, numChunk); + if (!status) goto ERROR; + if (!fileGroup) break; + + switch (type) { + case PPMERGE_TYPE_SHUTTER: + if (nThreads) { + // allocate a job + psThreadJob *job = psThreadJobAlloc ("PPMERGE_SHUTTER_CORRECTION", 0); + + // construct the arguments for this job + psArrayAdd (job->args, 1, outRO); + psArrayAdd (job->args, 1, fileGroup); + psArrayAdd (job->args, 1, psScalarAlloc(shutterRef, PS_TYPE_F32)) + psArrayAdd (job->args, 1, shutters->data[cellNum]); + psArrayAdd (job->args, 1, psScalarAlloc(iter, PS_TYPE_S32)); + psArrayAdd (job->args, 1, psScalarAlloc(rej, PS_TYPE_F32)); + psArrayAdd (job->args, 1, psScalarAlloc(maskVal, PS_TYPE_U8)); + + psThreadJobAddPending (job); + } else { + if (!pmShutterCorrectionGenerate(outRO, NULL, fileGroup->readouts, shutterRef, shutters->data[cellNum], iter, rej, maskVal)) { + goto ERROR; + } + fileGroup->busy = false; + } + break; + case PPMERGE_TYPE_DARK: + if (nThreads) { + // allocate a job + psThreadJob *job = psThreadJobAlloc ("PPMERGE_DARK_COMBINE", 0); + + // construct the arguments for this job + psArrayAdd (job->args, 1, outCell); + psArrayAdd (job->args, 1, fileGroup); + psArrayAdd (job->args, 1, darkOrdinates); + psArrayAdd (job->args, 1, darkNorm); + psArrayAdd (job->args, 1, psScalarAlloc(iter, PS_TYPE_S32)); + psArrayAdd (job->args, 1, psScalarAlloc(rej, PS_TYPE_F32)); + psArrayAdd (job->args, 1, psScalarAlloc(maskVal, PS_TYPE_U8)); + + psThreadJobAddPending (job); + } else { + if (!pmDarkCombine(outCell, fileGroup->readouts, darkOrdinates, darkNorm, iter, rej, maskVal)) { + goto ERROR; + } + fileGroup->busy = false; + } + break; + case PPMERGE_TYPE_BIAS: + case PPMERGE_TYPE_FLAT: + case PPMERGE_TYPE_FRINGE: + if (nThreads) { + // allocate a job + psThreadJob *job = psThreadJobAlloc ("PPMERGE_READOUT_COMBINE", 0); + + // construct the arguments for this job + psArrayAdd (job->args, 1, outRO); + psArrayAdd (job->args, 1, fileGroup); + psArrayAdd (job->args, 1, zeros); + psArrayAdd (job->args, 1, scales); + psArrayAdd (job->args, 1, combination); + + psThreadJobAddPending (job); + } else { + if (!pmReadoutCombine(outRO, fileGroup->readouts, zeros, scales, combination)) { + goto ERROR; + } + fileGroup->busy = false; + } + break; + default: + psAbort("Should never get here."); + } + } + + // wait for the threads to finish and manage results + if (nThreads) { + // wait here for the threaded jobs to finish + if (!psThreadPoolWait ()) { + psError(PS_ERR_UNKNOWN, false, "Unable to interpolate image."); + return false; + } + fprintf (stderr, "success for threaded jobs\n"); + + // we don't care about the results, just dump the done queue jobs + psThreadJob *job = NULL; + while ((job = psThreadJobGetDone()) != NULL) { + psFree (job); + } + } + + // Get list of cells for concepts averaging + psList *inCells = psListAlloc(NULL); // List of cells + for (int i = 0; i < numFiles; i++) { + pmReadout *readout = readouts->data[i]; // Readout of interest + psListAdd(inCells, PS_LIST_TAIL, readout->parent); + } + if (!pmConceptsAverageCells(outCell, inCells, NULL, NULL, true)) { + psError(PS_ERR_UNKNOWN, false, "Unable to average cell concepts."); + psFree(inCells); + psFree(outRO); + goto ERROR; + } + psFree(inCells); + + psFree(fileGroups); + + // Plug supplementary images into their own FPAs + { + pmCell *countsCell = pmFPAfileThisCell(config->files, view, "PPMERGE.OUTPUT.COUNT"); + pmReadout *countsRO = pmReadoutAlloc(countsCell); // Readout with count of inputs per pixel + psImage *counts = psMetadataLookupPtr(NULL, outRO->analysis, PM_READOUT_STACK_ANALYSIS_COUNT); + countsRO->image = psImageCopy(countsRO->image, counts, PS_TYPE_F32); + psMetadataRemoveKey(outRO->analysis, PM_READOUT_STACK_ANALYSIS_COUNT); + countsRO->data_exists = countsCell->data_exists = countsCell->parent->data_exists = true; + psFree(countsRO); + + pmCell *sigmaCell = pmFPAfileThisCell(config->files, view, "PPMERGE.OUTPUT.SIGMA"); + pmReadout *sigmaRO = pmReadoutAlloc(sigmaCell); // Readout with stdev per pixel + psImage *sigma = psMetadataLookupPtr(NULL, outRO->analysis, PM_READOUT_STACK_ANALYSIS_SIGMA); + sigmaRO->image = psImageCopy(sigmaRO->image, sigma, PS_TYPE_F32); + psMetadataRemoveKey(outRO->analysis, PM_READOUT_STACK_ANALYSIS_SIGMA); + sigmaRO->data_exists = sigmaCell->data_exists = sigmaCell->parent->data_exists = true; + psFree(sigmaRO); + } + + // Measure the fringes for this cell + // + // XXX Need to deal with multiple components: we will do this by building up the components + // Read the existing fringe measurements + // Use existing regions to measure fringe statistics + // Add the new fringe measurements to the existing fringe measurements. + // Write the appended fringe measurements. + // Read in the "output" file to get the existing components. + // Put the new readout into the cell after the existing readouts. + if (type == PPMERGE_TYPE_FRINGE && outRO) { + pmFringeRegions *regions = pmFringeRegionsAlloc(fringeNum, fringeSize, fringeSize, + fringeSmoothX, fringeSmoothY); + pmFringeStats *fringe = pmFringeStatsMeasure(regions, outRO, maskVal); + psFree(regions); + if (!fringe) { + psError(PS_ERR_UNKNOWN, false, "Unable to measure fringe statistics.\n"); + psFree(outRO); + goto ERROR; + } + + psArray *fringes = psArrayAlloc(1); // Array of fringes + fringes->data[0] = fringe; + + pmFringesFormat(outCell, NULL, fringes); + psFree(fringes); // Drop reference + } + + if (!ppStatsFPA(stats, outFPA, view, maskVal, 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 (outChip->data_exists) { + psList *inChips = psListAlloc(NULL); + for (int i=0; i < numFiles; i++) { + pmChip *chip = pmFPAviewThisChip(view, ((pmFPAfile *)inputs->data[i])->fpa); + psListAdd(inChips, PS_LIST_TAIL, chip); + } + if (!pmConceptsAverageChips(outChip, inChips, true)) { + psError(PS_ERR_UNKNOWN, false, "Unable to average Chip concepts."); + psFree(inChips); + goto ERROR; + } + psFree(inChips); + } + if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) { + goto ERROR; + } + } + + // Get list of FPAs for concepts averaging + psList *inFPAs = psListAlloc(NULL); // List of FPAs + for (int i = 0; i < numFiles; i++) { + pmFPAfile *input = inputs->data[i]; // Input file + psListAdd(inFPAs, PS_LIST_TAIL, input->fpa); + } + if (!pmConceptsAverageFPAs(output->fpa, inFPAs)) { + psError(PS_ERR_UNKNOWN, false, "Unable to average FPA concepts."); + psFree(inFPAs); + goto ERROR; + } + psFree(inFPAs); + + + if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) { + goto ERROR; + } + + psFree(view); + psFree(combination); + psFree(inputs); + psFree(masks); + psFree(weights); + psFree(stats); + return true; + +ERROR: + psFree(view); + psFree(combination); + psFree(inputs); + psFree(masks); + psFree(weights); + psFree(stats); + return false; +} + Index: /branches/eam_branch_20080719/ppMerge/src/ppMergeMask.c =================================================================== --- /branches/eam_branch_20080719/ppMerge/src/ppMergeMask.c (revision 18759) +++ /branches/eam_branch_20080719/ppMerge/src/ppMergeMask.c (revision 18759) @@ -0,0 +1,457 @@ +#include "ppMerge.h" + +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 + ) +{ + 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 + 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 + + bool smoothSuspect = psMetadataLookupBool(&mdok, config->arguments, "MASK.SMOOTH.SUSPECT"); // Radius to grow mask + float smoothScale = psMetadataLookupF32(&mdok, config->arguments, "MASK.SMOOTH.SCALE"); // Radius to grow mask + + psMaskType markVal; + psMaskType maskValRaw; + if (!pmConfigMaskSetBits (&maskValRaw, &markVal, config)) { + psError (PS_ERR_UNKNOWN, true, "Unable to define the mask bit values"); + return false; + } + + char *maskOutName = psMetadataLookupStr (&mdok, config->arguments, "MASK.SET.VALUE"); + psMaskType maskValOut = pmConfigMaskGet (maskOutName, config); + if (!maskValOut) { + psError (PS_ERR_UNKNOWN, true, "Undefined output mask bit value"); + return false; + } + + 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); + + int numCells = 1; + if (chipStats) { + // count the number of active cells for this chip: + numCells = 0; + for (int i = 0; i < outChip->cells->n; i++) { + pmCell *cell = outChip->cells->data[i]; + if (!cell->process) continue; + numCells ++; + } + } + + // For each input file, get the statistics, which can be calculated at the chip or cell levels + psVector *values = psVectorAlloc(sample, PS_TYPE_F32); // Pixel values for statistics + pmFPAview *inView = pmFPAviewAlloc(0); // View for input + for (int i = 0; i < numFiles; i++) { + pmFPAfileActivate(config->files, false, NULL); + psArray *files = ppMergeFileActivateSingle(config, PPMERGE_FILES_INPUT, true, i); // Input files + pmFPAfile *input = files->data[0]; // Input file + psFree(files); + pmFPA *inFPA = input->fpa; // Input FPA + *inView = *view; + + int valueIndex = 0; // Index for vector of pixel values + + pmCell *inCell; // Input cell + while ((inCell = pmFPAviewNextCell(inView, inFPA, 1))) { + + // the output FPA structure carries the information about which cells to process + pmCell *outCell = pmFPAfileThisCell(config->files, inView, "PPMERGE.OUTPUT.MASK"); // Output cell + if (!outCell->process) continue; + + 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; + } + + pmReadout *outRO = NULL; // Output readout + if (outCell->readouts && outCell->readouts->n == 1) { + outRO = psMemIncrRefCounter(outCell->readouts->data[0]); + } else { + outRO = pmReadoutAlloc(outCell); + } + psImage *outMask = outRO->mask; // Output mask image (for iterative generation of mask) + + psImage *image = readout->image, *mask = readout->mask; // Image and mask + int numCols = readout->image->numCols, numRows = readout->image->numRows; // Image size + int numPix = numCols * numRows; // Number of pixels + int 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] & maskValRaw)) continue; + if (outMask && (outMask->data.PS_TYPE_MASK_DATA[y][x] & maskValOut)) continue; + if (!isfinite(image->data.F32[y][x])) 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; + } + + float mean = psStatsGetValue(statistics, meanStat); + float stdev = psStatsGetValue(statistics, stdevStat); + + // this function increments the count for each suspect pixel in each input plane + // maskValRaw is used to test for valid input pixels + if (!pmMaskFlagSuspectPixels(outRO, readout, mean, stdev, maskSuspect, maskValRaw)) { + psError(PS_ERR_UNKNOWN, false, "Unable to find suspect values in file %d", i); + psFree(inView); + psFree(readout); + goto MERGE_MASK_ERROR; + } + pmCellFreeData(inCell); + + if (!pmFPAfileIOChecks(config, inView, PM_FPA_AFTER)) { + psFree(inView); + psFree(readout); + goto MERGE_MASK_ERROR; + } + } + psFree(readout); + psFree(outRO); + } + + // Additional run through cells if we want chip-level statistics + if (chipStats && valueIndex > 0) { + values->n = valueIndex; + if (!psVectorStats(statistics, values, NULL, NULL, 0)) { + psError(PS_ERR_UNKNOWN, false, "Unable to do statistics on chip."); + goto MERGE_MASK_ERROR; + } + inView->cell = -1; + while ((inCell = pmFPAviewNextCell(inView, inFPA, 1))) { + + // the output FPA structure carries the information about which cells to process + pmCell *outCell = pmFPAfileThisCell(config->files, inView, "PPMERGE.OUTPUT.MASK"); // Output cell + if (!outCell->process) continue; + + pmHDU *hdu = pmHDUFromCell(inCell); // HDU for cell + if (!hdu || hdu->blankPHU) continue; + + pmReadout *readout = inCell->readouts->data[0]; // Readout of interest + + inView->readout = 0; + pmReadout *outRO = pmFPAfileThisReadout(config->files, inView, "PPMERGE.OUTPUT.MASK"); + + float mean = psStatsGetValue(statistics, meanStat); + float stdev = psStatsGetValue(statistics, stdevStat); + + if (!pmMaskFlagSuspectPixels(outRO, readout, mean, stdev, maskSuspect, maskValRaw)) { + psError(PS_ERR_UNKNOWN, false, "Unable to find suspect values in file %d", i); + goto MERGE_MASK_ERROR; + } + pmCellFreeData(inCell); + + inView->readout = -1; + if (!pmFPAfileIOChecks(config, inView, PM_FPA_AFTER)) { + psFree(inView); + goto MERGE_MASK_ERROR; + } + } + } + } + psFree(inView); + psFree(statistics); statistics = NULL; + psFree(values); values = NULL; + + + // Another run through the chip to threshold on the suspects + pmFPAfileActivate(config->files, false, NULL); + if (writeOut) { + ppMergeFileActivate(config, PPMERGE_FILES_OUTPUT, true); + } + pmFPA *outFPA = outChip->parent; // Output FPA + pmCell *outCell; // Output cell + pmFPAview *outView = pmFPAviewAlloc(0); // View into output FPA + *outView = *view; + while ((outCell = pmFPAviewNextCell(outView, outFPA, 1))) { + + // skip inactive cells + if (!outCell->process) continue; + + pmHDU *hdu = pmHDUFromCell(outCell); // HDU for cell + if (!hdu || hdu->blankPHU) continue; + + psTrace("ppMerge", 1, "Getting bad pixels for chip %d cell %d", outView->chip, outView->cell); + + assert(outCell->readouts && outCell->readouts->n == 1); + pmReadout *outRO = outCell->readouts->data[0]; // Output readout + + if (smoothSuspect) { + // XXX test output of suspect pixel image + psImage *suspects = psMetadataLookupPtr(NULL, outRO->analysis, PM_MASK_ANALYSIS_SUSPECT); // Suspect img + assert (suspects); + psImageSmooth (suspects, smoothScale, 3); // extend smoothing region to 3-sigma + } + + // set the bad pixels to the value 'maskVal' + if (!pmMaskIdentifyBadPixels(outRO, maskValOut, maskBad, maskMode)) { + psError(PS_ERR_UNKNOWN, false, "Unable to mask bad pixels"); + goto MERGE_MASK_ERROR; + } + + // Supplementary outputs + { + // The counts image is fairly useless, but it preserves the model + pmCell *countsCell = pmFPAfileThisCell(config->files, outView, "PPMERGE.OUTPUT.COUNT"); + pmReadout *countsRO = pmReadoutAlloc(countsCell); // Readout with count of inputs per pixel + countsRO->image = psImageAlloc(outRO->mask->numCols, outRO->mask->numRows, PS_TYPE_F32); + psImageInit(countsRO->image, numFiles); + countsRO->data_exists = countsCell->data_exists = countsCell->parent->data_exists = true; + psFree(countsRO); + + pmCell *sigmaCell = pmFPAfileThisCell(config->files, outView, "PPMERGE.OUTPUT.SIGMA"); + pmReadout *sigmaRO = pmReadoutAlloc(sigmaCell); // Readout with suspect image + psImage *suspect = psMetadataLookupPtr(NULL, outRO->analysis, PM_MASK_ANALYSIS_SUSPECT); + sigmaRO->image = psImageCopy(sigmaRO->image, suspect, PS_TYPE_F32); + psMetadataRemoveKey(outRO->analysis, PM_MASK_ANALYSIS_SUSPECT); + sigmaRO->data_exists = sigmaCell->data_exists = sigmaCell->parent->data_exists = true; + psFree(sigmaRO); + } + + if (maskGrow > 0) { + psImage *grown = psImageGrowMask(NULL, outRO->mask, maskValOut, maskGrow, maskValOut); // Grown mask + psFree(outRO->mask); + outRO->mask = grown; + } + + if (writeOut) { + if (!pmFPAfileIOChecks(config, outView, PM_FPA_BEFORE)) { + psFree(outView); + goto MERGE_MASK_ERROR; + } + + outRO->data_exists = outCell->data_exists = outChip->data_exists = true; + + // Average concepts + psList *cells = psListAlloc(NULL); // List of cells, for concept averaging + for (int i = 0; i < numFiles; i++) { + pmFPAfile *inFile = pmFPAfileSelectSingle(config->files, "PPMERGE.INPUT", i); // Input file + pmCell *inCell = pmFPAviewThisCell(outView, inFile->fpa); // Input cell + psListAdd(cells, PS_LIST_TAIL, inCell); + } + if (!pmConceptsAverageCells(outCell, cells, NULL, NULL, true)) { + psError(PS_ERR_UNKNOWN, false, "Unable to average cell concepts."); + psFree(cells); + psFree(outRO); + psFree(outView); + return false; + } + psFree(cells); + + // Statistics on the merged cell using a fake image + outRO->image = psImageAlloc(outRO->mask->numCols, outRO->mask->numRows, PS_TYPE_F32); + psImageInit(outRO->image, 1.0); + if (!ppStatsFPA(stats, outRO->parent->parent->parent, outView, maskValOut, config)) { + psError(PS_ERR_UNEXPECTED_NULL, true, "Unable to generate stats for image."); + psFree(outRO); + psFree(outView); + return false; + } + psFree(outRO->image); + outRO->image = NULL; + + // Write + if (!pmFPAfileIOChecks(config, outView, PM_FPA_AFTER)) { + psFree(outView); + goto MERGE_MASK_ERROR; + } + } + } + psFree(outView); + + ppMergeFileActivate(config, PPMERGE_FILES_ALL, true); + + return true; + + +MERGE_MASK_ERROR: + psFree(statistics); + psFree(values); + return false; +} + +bool ppMergeMask(pmConfig *config) +{ + assert(config); + + bool mdok; // Status of MD lookup + int numFiles = psMetadataLookupS32(NULL, config->arguments, "INPUTS.NUM"); // Number of inputs + bool haveMasks = psMetadataLookupBool(&mdok, config->arguments, "INPUTS.MASKS"); // Do we have masks? + bool haveWeights = psMetadataLookupBool(&mdok, config->arguments, "INPUTS.WEIGHTS"); // Do we have weights? + int iter = psMetadataLookupS32(NULL, config->arguments, "ITER"); // Number of rejection iterations + + PS_ASSERT_INT_POSITIVE(iter, false); + + pmFPAview *view = pmFPAviewAlloc(0); // View to component of interest + psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS, 0); // Random number generator + + psMetadata *stats = NULL; // Statistics for output + if (psMetadataLookup(config->arguments, "STATS.NAME")) { + stats = psMetadataAlloc(); + psMetadataAddMetadata(config->arguments, PS_LIST_TAIL, "STATS.DATA", 0, "Statistics output", stats); + } + + psArray *inputs = ppMergeFileDataLevel(config, "PPMERGE.INPUT", PM_FPA_LEVEL_READOUT); // Input images + psFree(inputs); + if (haveMasks) { + psArray *masks = ppMergeFileDataLevel(config, "PPMERGE.INPUT.MASK", PM_FPA_LEVEL_READOUT); + psFree(masks); + } + if (haveWeights) { + psArray *weights = ppMergeFileDataLevel(config, "PPMERGE.INPUT.WEIGHT", PM_FPA_LEVEL_READOUT); + psFree(weights); + } + + if (!ppMergeFileActivate(config, PPMERGE_FILES_INPUT, true)) { + psError(PS_ERR_UNKNOWN, false, "Unable to activate files."); + goto PPMERGE_MASK_ERROR; + } + if (!ppMergeFileActivate(config, PPMERGE_FILES_OUTPUT, true)) { + psError(PS_ERR_UNKNOWN, false, "Unable to activate files."); + goto PPMERGE_MASK_ERROR; + } + + // XXX this function should use pmConfigMaskReadHeader () to get the named values defined + // for the input masks. + + 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 (!outChip->process) continue; + + 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 (outChip->data_exists) { + psList *inChips = psListAlloc(NULL); + for (int i=0; i < numFiles; i++) { + pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PPMERGE.INPUT", i); // Input file + pmChip *chip = pmFPAviewThisChip(view, file->fpa); + psListAdd(inChips, PS_LIST_TAIL, chip); + } + + // XXX I need to call pmConfigMaskWriteHeader for the PHU somewhere, after it is created! + + if (!pmConceptsAverageChips(outChip, inChips, true)) { + psError(PS_ERR_UNKNOWN, false, "Unable to average Chip concepts."); + psFree(inChips); + goto PPMERGE_MASK_ERROR; + } + psFree(inChips); + } + 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; +} + Index: /branches/eam_branch_20080719/ppMerge/src/ppMergeReadChunk.c =================================================================== --- /branches/eam_branch_20080719/ppMerge/src/ppMergeReadChunk.c (revision 18759) +++ /branches/eam_branch_20080719/ppMerge/src/ppMergeReadChunk.c (revision 18759) @@ -0,0 +1,81 @@ +# include "ppMerge.h" + +ppMergeFileGroup *ppMergeReadChunk (psArray *fileGroups, pmConfig *config, int numChunk) { + + // select an available fileGroup + + bool haveMasks = psMetadataLookupBool(&mdok, arguments, "INPUTS.MASKS"); // Do we have masks? + bool haveWeights = psMetadataLookupBool(&mdok, arguments, "INPUTS.WEIGHTS"); // Do we have weights? + + while (1) { + // check for any fileGroups which can read data + for (int j = 0; j < fileGroups->n; j++) { + ppMergeFileGroup *fileGroup = fileGroups->data[j]; + if (fileGroup->read) continue; + + psArray *readouts = fileGroup->readouts; + + psTrace("ppStack", 2, "Reading data for chunk %d into fileGroup %d....n", numChunk, j); + for (int i = 0; i < readouts->n; i++) { + pmReadout *inRO = readouts->data[i]; // Input readout + + // Read a chunk from a file + pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PPMERGE.INPUT", i); + if (!pmReadoutReadChunk(inRO, file->fits, 0, rows, 0, config)) { + psError(PS_ERR_IO, false, "Unable to read chunk %d for file PPMERGE.INPUT %d", numChunk, i); + return NULL; + } + + if (haveMasks) { + pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PPMERGE.INPUT.MASK", i); + if (!pmReadoutReadChunkMask(inRO, file->fits, 0, rows, 0, config)) { + psError(PS_ERR_IO, false, "Unable to read chunk %d for file PPMERGE.INPUT.MASK %d", numChunk, NAME, i); + return NULL; + } + } + + if (haveWeights) { + pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PPMERGE.INPUT.WEIGHT", i); + if (!pmReadoutReadChunkWeight(inRO, file->fits, 0, rows, 0, config)) { + psError(PS_ERR_IO, false, "Unable to read chunk %d for file PPMERGE.INPUT.WEIGHT %d", numChunk, NAME, i); + return NULL; + } + } + } + fileGroup->read = fileGroup->busy = true; + return fileGroup; + } + + // check for any fileGroups which are done processing + bool wait = true; + bool more = true; + for (int j = 0; j < fileGroups->n; j++) { + ppMergeFileGroup *fileGroup = fileGroups->data[j]; + if (!fileGroup->read || fileGroup->busy) continue; + + wait = false; + psArray *readouts = fileGroup->readouts; + // any more data to be read? + for (int i = 0; i < readouts->n && more; i++) { + pmReadout *inRO = readouts->data[i]; + + pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PPMERGE.INPUT", i); + more &= pmReadoutMore(inRO, file->fits, 0, rows, config); + + if (haveMasks) { + pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PPMERGE.INPUT.MASK", i); + more &= pmReadoutMoreMask(inRO, file->fits, 0, rows, config); + } + if (haveWeights) { + pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PPMERGE.INPUT.WEIGHT", i); + more &= pmReadoutMoreWeight(inRO, file->fits, 0, rows, config); + } + } + fileGroup->read = false; + } + if (!more) return NULL; + + if (wait) usleep (10000); + } + return NULL; +} Index: /branches/eam_branch_20080719/ppMerge/src/ppMergeScaleZero.c =================================================================== --- /branches/eam_branch_20080719/ppMerge/src/ppMergeScaleZero.c (revision 18759) +++ /branches/eam_branch_20080719/ppMerge/src/ppMergeScaleZero.c (revision 18759) @@ -0,0 +1,211 @@ +#include "ppMerge.h" + +// Get the scale and zero for each chip of each input +bool ppMergeScaleZero(pmConfig *config) +{ + assert(config); + + 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 + int shutterSize = psMetadataLookupS32(NULL, config->arguments, "SHUTTER.SIZE"); // Size of shutter region + + psMaskType maskVal = pmConfigMaskGet("MASK.VALUE", config); // Value to mask + + psVector *gains = NULL; // Gains for each cell + psArray *shutters = NULL; // Shutter data for each cell + psStats *stats = NULL; // Statistics for background + psImage *background = NULL; // Background measurements per cell per file + + switch (type) { + case PPMERGE_TYPE_BIAS: + case PPMERGE_TYPE_DARK: + // Nothing to measure + return true; + case PPMERGE_TYPE_FLAT: + case PPMERGE_TYPE_FRINGE: + gains = psVectorAlloc(numCells, PS_TYPE_F32); + background = psImageAlloc(numCells, numInputs, PS_TYPE_F32); + psImageInit(background, NAN); + stats = psStatsAlloc(meanStat); + break; + case PPMERGE_TYPE_SHUTTER: + shutters = psArrayAlloc(numCells); + break; + case PPMERGE_TYPE_MASK: + default: + break; + } + psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS, 0); // Random number generator + pmFPAview *view = NULL; // View into FPA + + for (int i = 0; i < numInputs; i++) { + pmFPAfileActivate(config->files, false, NULL); + psArray *files = ppMergeFileActivateSingle(config, PPMERGE_FILES_INPUT, true, i); // Activated files + pmFPAfile *input = files->data[0]; // Representative file; should be the image (not mask or weight) + pmFPA *fpa = input->fpa; // FPA of interest + view = pmFPAviewAlloc(0); // View to component of interest + int cellNum = 0; // Index for cell + if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) { + goto ERROR; + } + pmChip *chip; // Chip of interest + while ((chip = pmFPAviewNextChip(view, fpa, 1))) { + if (!chip->file_exists) { + continue; + } + if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) { + goto ERROR; + } + + pmCell *cell; // Cell of interest + while ((cell = pmFPAviewNextCell(view, fpa, 1))) { + if (!cell->file_exists) { + continue; + } + if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) { + goto ERROR; + } + + if (!cell->data_exists) { + continue; + } + + if (cell->readouts->n > 1) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, + "File %d chip %d cell %d contains more than one readout (%ld)", + i, view->chip, view->cell, cell->readouts->n); + goto ERROR; + } + pmReadout *readout = cell->readouts->data[0]; // Readout of interest + + switch (type) { + case PPMERGE_TYPE_FLAT: + case PPMERGE_TYPE_FRINGE: { + // Extract the gain + float gain = psMetadataLookupF32(NULL, cell->concepts, "CELL.GAIN"); // Cell gain + if (!isfinite(gain)) { + // psError(PS_ERR_BAD_PARAMETER_VALUE, false, + // "CELL.GAIN for file %d chip %d cell %d is not set.", + // i, view->chip, view->cell); + // goto ERROR; + psWarning ("CELL.GAIN for file %d chip %d cell %d is NaN", i, view->chip, view->cell); + } + 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; + psWarning ("Unable to get statistics for file %d chip %d cell %d", i, view->chip, view->cell); + background->data.F32[i][cellNum] = NAN; + } else { + background->data.F32[i][cellNum] = psStatsGetValue(stats, meanStat); + } + break; + } + case PPMERGE_TYPE_SHUTTER: { + pmShutterCorrectionData *shutter = shutters->data[cellNum]; // Shutter correction data + if (!shutter) { + shutter = pmShutterCorrectionDataAlloc(readout->image->numCols, + readout->image->numRows, + shutterSize); + shutters->data[cellNum] = shutter; + } + if (!pmShutterCorrectionAddReadout(shutter, readout, meanStat, stdevStat, + maskVal, rng)) { + psError(PS_ERR_UNKNOWN, false, + "Can't add file %d chip %d cell %d to shutter correction.", + i, view->chip, view->cell); + goto ERROR; + } + break; + } + case PPMERGE_TYPE_MASK: + default: + psAbort("Should never get here."); + } + + cellNum++; + + if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) { + goto ERROR; + } + } + + if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) { + goto ERROR; + } + } + + if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) { + goto ERROR; + } + + psFree(view); + +#if 0 + // Reset files for reading again + for (int i = 0; i < files->n; i++) { + pmFPAfile *file = files->data[i]; // File of interest + } +#endif + psFree(files); + } + + psFree(rng); rng = NULL; + psFree(stats); stats = NULL; + + // Store results + switch (type) { + case PPMERGE_TYPE_FRINGE: + psMetadataAddImage(config->arguments, PS_LIST_TAIL, "ZEROS", 0, + "Zero to subtract from each input cell", background); + // Flow through + case PPMERGE_TYPE_FLAT: { + // Need to normalize over the focal plane + if (psTraceGetLevel("ppMerge") > 9) { + for (int i = 0; i < gains->n; i++) { + psTrace("ppMerge", 10, "Gain for cell %d is %f\n", i, gains->data.F32[i]); + } + } + psVector *fluxes = NULL; // Solution to fluxes + if (!pmFlatNormalize(&fluxes, &gains, background)) { + psError(PS_ERR_UNKNOWN, false, "Normalisation failed to converge --- continuing anyway."); + psFree(fluxes); + goto ERROR; + } + + psMetadataAddVector(config->arguments, PS_LIST_TAIL, "SCALES", 0, + "Scale to divide into each input file", fluxes); + psFree(fluxes); // Drop reference + break; + } + case PPMERGE_TYPE_SHUTTER: + psMetadataAddArray(config->arguments, PS_LIST_TAIL, "SHUTTER", 0, + "Shutter data", shutters); + break; + case PPMERGE_TYPE_MASK: + default: + psAbort("Should never get here."); + } + + psFree(background); + psFree(shutters); + return true; + +ERROR: + // Common path for errors + psFree(gains); + psFree(background); + psFree(shutters); + psFree(rng); + psFree(stats); + psFree(view); + return false; +} + Index: /branches/eam_branch_20080719/ppMerge/src/ppMergeThreadLauncher.c =================================================================== --- /branches/eam_branch_20080719/ppMerge/src/ppMergeThreadLauncher.c (revision 18759) +++ /branches/eam_branch_20080719/ppMerge/src/ppMergeThreadLauncher.c (revision 18759) @@ -0,0 +1,93 @@ +# include "ppMerge.h" + +// each thread runs this function, starting a new job when it finished with an old one +// it is called with a (void *) pointer to its own thread pointer +void *ppMergeThreadLauncher (void *data) { + + psThread *self = data; + psThreadJob *job = NULL; + + while (1) { + + // if we get an error, just wait until we are cleared or killed + while (self->fault) { + usleep (10000); + } + + // request a new job, if there are none available, sleep a bit + // we have to lock here so the job queue cannot be empty yet no threads busy + psThreadLock(); + while ((job = psThreadJobGetPending ()) == NULL) { + psThreadUnlock(); + usleep (10000); + } + self->busy = true; + psThreadUnlock(); + + // list all allowed job types here + + // pmReadoutCombine + if (!strcmp (job->type, "PPMERGE_READOUT_COMBINE")) { + psAssert (job->args->n == 5, "invalid number of arguments to pmReadoutCombine"); + + pmReadout *output = job->args->data[0]; + ppMergeFileGroup *fileGroup = job->args->data[1]; + psVector *zero = job->args->data[2]; + psVector *scale = job->args->data[3]; + pmCombineParams *params = job->args->data[4]; + + bool status = pmReadoutCombine (output, fileGroup->readouts, zero, scale, params); + if (!status) { + self->fault = true; + } + // we do not have to lock here because this transition is not tied to the job queue + fileGroup->busy = false; + self->busy = false; + continue; + } + + // pmDarkCombine + if (!strcmp (job->type, "PPMERGE_DARK_COMBINE")) { + psAssert (job->args->n == 7, "invalid number of arguments to pmDarkCombine"); + + pmCell *outCell = job->args->data[0]; + ppMergeFileGroup *fileGroup = job->args->data[1]; + psArray *darkOrdinates = job->args->data[2]; + psString darkNorm = job->args->data[3]; + psScalar *iter = job->args->data[4]; + psScalar *rej = job->args->data[5]; + psScalar *maskVal = job->args->data[6]; + + bool status = pmDarkCombine(outCell, fileGroup->readouts, darkOrdinates, darkNorm, iter->data.S32, rej->data.F32, maskVal->data.U8); + if (!status) { + self->fault = true; + } + // we do not have to lock here because this transition is not tied to the job queue + fileGroup->busy = false; + self->busy = false; + continue; + } + + // pmShutterCorrectionGenerate + if (!strcmp (job->type, "PPMERGE_SHUTTER_CORRECTION")) { + psAssert (job->args->n == 7, "invalid number of arguments to pmDarkCombine"); + + pmReadout *output = job->args->data[0]; + ppMergeFileGroup *fileGroup = job->args->data[1]; + psScalar *shutterRef = job->args->data[2]; + pmShutterCorrectionData *data = job->args->data[3]; + psScalar *iter = job->args->data[4]; + psScalar *rej = job->args->data[5]; + psScalar *maskVal = job->args->data[6]; + + bool status = pmShutterCorrectionGenerate(output, NULL, fileGroup->readouts, shutterRef->data.F32, data, iter->data.S32, rej->data.F32, maskVal->data.U8); + if (!status) { + self->fault = true; + } + // we do not have to lock here because this transition is not tied to the job queue + fileGroup->busy = false; + self->busy = false; + continue; + } + } +} Index: /branches/eam_branch_20080719/psLib/src/types/psMetadata.c =================================================================== --- /branches/eam_branch_20080719/psLib/src/types/psMetadata.c (revision 18759) +++ /branches/eam_branch_20080719/psLib/src/types/psMetadata.c (revision 18759) @@ -0,0 +1,1585 @@ +/** @file psMetadata.c + * + * + * @brief Contains metadata structures, enumerations and functions prototypes. + * + * This file defines metadata item, metadata type, metadata flags, metadata containers, and function + * prototypes necessary creating psLib metadata APIs + * + * @ingroup Metadata + * + * @author Robert DeSonia, MHPCC + * @author Ross Harman, MHPCC + * + * @version $Revision: 1.171 $ $Name: not supported by cvs2svn $ + * @date $Date: 2008-07-24 22:42:28 $ + * + * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii + */ + +#ifdef HAVE_CONFIG_H +# include "config.h" +#endif + +/******************************************************************************/ +/* INCLUDE FILES */ +/******************************************************************************/ +#include +#include +#include + +#include "fitsio.h" +#include "psType.h" +#include "psMemory.h" +#include "psError.h" +#include "psAbort.h" +#include "psList.h" +#include "psHash.h" +#include "psVector.h" +#include "psMetadata.h" +#include "psLookupTable.h" +#include "psString.h" + +#include "psAssert.h" +#include "psLogMsg.h" +#include "psTrace.h" + +/******************************************************************************/ +/* DEFINE STATEMENTS */ +/******************************************************************************/ + +#define MAX_STRING_SIZE 1024 // Maximum string size for psMetadatItem + +/******************************************************************************/ +/* TYPE DEFINITIONS */ +/******************************************************************************/ + +// None + +/*****************************************************************************/ +/* GLOBAL VARIABLES */ +/*****************************************************************************/ + +// None + +/*****************************************************************************/ +/* FILE STATIC VARIABLES */ +/*****************************************************************************/ + +static psS32 metadataId = 0; + +/*****************************************************************************/ +/* FUNCTION IMPLEMENTATION - LOCAL */ +/*****************************************************************************/ + +static psMetadataItem* makeMetaMulti(psHash* table, + const char* key, + psMetadataItem* existing) +{ + + if (existing != NULL && existing->type == PS_DATA_METADATA_MULTI) { + return existing; + } + + // move any existing entry into a psList + psList* newList = psListAlloc(existing); + + psMetadataItem* item = psMetadataItemAlloc(key, + PS_DATA_METADATA_MULTI, + "", + newList); + + if (existing != NULL) { + psHashRemove(table,key); // take out the old entry + } + + psHashAdd(table, key, item); // put in the new entry + psMemDecrRefCounter(item); // get rid of extra reference + + // free local references of newly allocated items. + psFree(newList); + + return item; +} + +static void metadataItemFree(psMetadataItem* metadataItem) +{ + psDataType type; + + type = metadataItem->type; + + + psMemDecrRefCounter(metadataItem->name); + psMemDecrRefCounter(metadataItem->comment); + + if(!PS_DATA_IS_PRIMITIVE(type)) { + psFree(metadataItem->data.V); + } +} + +static void metadataIteratorFree(psMetadataIterator* iter) +{ + psFree(iter->iter); + + if (iter->regex != NULL) { + regfree(iter->regex); + psFree(iter->regex); + } +} + +static void metadataFree(psMetadata* metadata) +{ + psFree(metadata->list); + psFree(metadata->hash); +} + +/*****************************************************************************/ +/* FUNCTION IMPLEMENTATION - PUBLIC */ +/*****************************************************************************/ + +psMetadataItem* p_psMetadataItemAlloc(const char *file, + unsigned int lineno, + const char *func, + const char *name, + psDataType type, + const char *comment, + ...) +{ + va_list argPtr; + psMetadataItem* metadataItem = NULL; + + // Get the variable list parameters to pass to allocation function + va_start(argPtr, comment); + + // Call metadata item allocation + metadataItem = p_psMetadataItemAllocV(file, lineno, func, name, type, comment, argPtr); + + // Clean up stack after variable arguement has been used + va_end(argPtr); + + return metadataItem; +} + +bool psMemCheckMetadataItem(psPtr ptr) +{ + PS_ASSERT_PTR(ptr, false); + return ( psMemGetDeallocator(ptr) == (psFreeFunc)metadataItemFree ); +} + + +#define METADATAITEM_ALLOC_TYPE(NAME,TYPE,METATYPE) \ +psMetadataItem* psMetadataItemAlloc##NAME(const char* name, \ + const char* comment, \ + TYPE value) \ +{ \ + return psMetadataItemAlloc(name, METATYPE, comment, value); \ +} + +METADATAITEM_ALLOC_TYPE(Str,const char*,PS_DATA_STRING) +METADATAITEM_ALLOC_TYPE(F32,psF32,PS_DATA_F32) +METADATAITEM_ALLOC_TYPE(F64,psF64,PS_DATA_F64) +METADATAITEM_ALLOC_TYPE(S8,psS8,PS_DATA_S8) +METADATAITEM_ALLOC_TYPE(S16,psS16,PS_DATA_S16) +METADATAITEM_ALLOC_TYPE(S32,psS32,PS_DATA_S32) +METADATAITEM_ALLOC_TYPE(S64,psS64,PS_DATA_S64) +METADATAITEM_ALLOC_TYPE(U8,psU8,PS_DATA_U8) +METADATAITEM_ALLOC_TYPE(U16,psU16,PS_DATA_U16) +METADATAITEM_ALLOC_TYPE(U32,psU32,PS_DATA_U32) +METADATAITEM_ALLOC_TYPE(U64,psU64,PS_DATA_U64) +METADATAITEM_ALLOC_TYPE(Bool,bool,PS_DATA_BOOL) + +psMetadataItem* psMetadataItemAllocPtr(const char *name, + psDataType type, + const char *comment, + psPtr value) +{ + return (psMetadataItemAlloc(name, type, comment, value)); +} + + +psMetadataItem* p_psMetadataItemAllocV(const char *file, + unsigned int lineno, + const char *func, + const char *name, + psDataType type, + const char *comment, + va_list argPtr) +{ + PS_ASSERT_STRING_NON_EMPTY(name, NULL); + PS_ASSERT_INT_POSITIVE(type, NULL); + + // Allocate metadata item + psMetadataItem *metadataItem = (psMetadataItem*) p_psAlloc(file, lineno, func, sizeof(psMetadataItem)); + metadataItem->data.V = NULL; + + // Set deallocator + psMemSetDeallocator(metadataItem, (psFreeFunc) metadataItemFree); + + // XXX Does it make sense allocate a NULL string here? It seems like this + // is just a waste of memory but something, somewhere, is probably + // dependant on this behavior so I'm leaving it in place for the time + // being. -JH + + // Allocate and set metadata item comment + if (comment == NULL) { + // Per SDRS, null isn't allowed, must use "" instead + metadataItem->comment = p_psStringCopy(file, lineno, func, ""); + } else { + metadataItem->comment = p_psStringCopy(file, lineno, func, comment); + } + + // Set metadata item unique id + *(psS32 *)(&metadataItem->id) = ++metadataId; + + // Set metadata item type + metadataItem->type = type & PS_METADATA_TYPE_MASK; + + // Allocate and set metadata item name + metadataItem->name = p_psStringCopy(file, lineno, func, name); + + // Set metadata item value + switch(metadataItem->type) { + case PS_DATA_BOOL: + metadataItem->data.B = (bool)va_arg(argPtr, psS32); + break; + case PS_DATA_S8: + metadataItem->data.S8 = (psS8)va_arg(argPtr, psS32); + break; + case PS_DATA_S16: + metadataItem->data.S16 = (psS16)va_arg(argPtr, psS32); + break; + case PS_DATA_S32: + metadataItem->data.S32 = (psS32)va_arg(argPtr, psS32); + break; + case PS_DATA_S64: + metadataItem->data.S64 = (psS64)va_arg(argPtr, psS64); + break; + case PS_DATA_U8: + metadataItem->data.U8 = (psU8)va_arg(argPtr, psU32); + break; + case PS_DATA_U16: + metadataItem->data.U16 = (psU16)va_arg(argPtr, psU32); + break; + case PS_DATA_U32: + metadataItem->data.U32 = (psU32)va_arg(argPtr, psU32); + break; + case PS_DATA_U64: + metadataItem->data.U64 = (psU64)va_arg(argPtr, psU64); + break; + case PS_DATA_F32: + metadataItem->data.F32 = (psF32)va_arg(argPtr, psF64); + break; + case PS_DATA_F64: + metadataItem->data.F64 = (psF64)va_arg(argPtr, psF64); + break; + case PS_DATA_STRING: { + // Copy input strings, so they can be messed with + char *string = va_arg(argPtr, char*); + metadataItem->data.str = (string ? p_psStringCopy(file, lineno, func, string) : NULL); + break; + } + case PS_DATA_ARRAY: // psArray + case PS_DATA_BITSET: // psBitSet + case PS_DATA_CUBE: // psCube + case PS_DATA_FITS: // psFits + case PS_DATA_HASH: // psHash + case PS_DATA_HISTOGRAM: // psHistogram + case PS_DATA_IMAGE: // psImage + case PS_DATA_KERNEL: // psKernel + case PS_DATA_LIST: // psList + case PS_DATA_LOOKUPTABLE: // psLookupTable + case PS_DATA_METADATA: // psMetadata + case PS_DATA_METADATAITEM: // psMetadataItem + case PS_DATA_MINIMIZATION: // psMinimization + case PS_DATA_PIXELS: // psPixels + case PS_DATA_PLANE: // psPlane + case PS_DATA_PLANEDISTORT: // psPlaneDistort + case PS_DATA_PLANETRANSFORM: // psPlaneTransform + case PS_DATA_POLYNOMIAL1D: // psPolynomial1D + case PS_DATA_POLYNOMIAL2D: // psPolynomial2D + case PS_DATA_POLYNOMIAL3D: // psPolynomial3D + case PS_DATA_POLYNOMIAL4D: // psPolynomial4D + case PS_DATA_PROJECTION: // psProjection + case PS_DATA_REGION: // psRegion + case PS_DATA_SCALAR: // psScalar + case PS_DATA_SPHERE: // psSphere + case PS_DATA_SPHEREROT: // psSphereTransform + case PS_DATA_SPLINE1D: // psSpline1D + case PS_DATA_STATS: // psStats + case PS_DATA_TIME: // psTime + case PS_DATA_VECTOR: // psVector + case PS_DATA_UNKNOWN: // "other" + case PS_DATA_METADATA_MULTI: + // Copy of input data not performed due to variability of data types + metadataItem->data.V = p_psMemIncrRefCounter(file, lineno, func, va_arg(argPtr, psPtr)); + break; + default: + psError(PS_ERR_BAD_PARAMETER_VALUE, true, _("Specified psDataType, %d, is not supported."), type); + psFree(metadataItem); + metadataItem = NULL; + break; + } + + return metadataItem; +} + +psMetadata* p_psMetadataAlloc(const char *file, + unsigned int lineno, + const char *func) +{ + psList* list = NULL; + psHash* hash = NULL; + psMetadata* metadata = NULL; + + // Allocate metadata + metadata = (psMetadata*) p_psAlloc(file, lineno, func, sizeof(psMetadata)); + // Set deallocator + psMemSetDeallocator(metadata, (psFreeFunc) metadataFree); + + // Allocate metadata's internal containers + list = (psList*) psListAlloc(NULL); + hash = (psHash*) psHashAlloc(10); + + metadata->list = list; + metadata->hash = hash; + + return metadata; +} + + +bool psMemCheckMetadata(psPtr ptr) +{ + PS_ASSERT_PTR(ptr, false); + return ( psMemGetDeallocator(ptr) == (psFreeFunc)metadataFree ); +} + + +psMetadataItem *p_psMetadataItemCopy(const char *file, + unsigned int lineno, + const char *func, + const psMetadataItem *in) +{ + PS_ASSERT_METADATA_ITEM_NON_NULL(in,NULL); + + psMetadataItem *newItem = NULL; // New metadata item, to be returned + + #define PS_METADATA_ITEM_COPY_CASE(NAME,TYPE) \ +case PS_DATA_##NAME: \ + newItem = p_psMetadataItemAlloc(file, lineno, func, in->name, PS_DATA_##NAME, in->comment, in->data.TYPE); \ + break; \ + + switch (in->type) { + // Simple types + PS_METADATA_ITEM_COPY_CASE(BOOL,B); + PS_METADATA_ITEM_COPY_CASE(S8,S8); + PS_METADATA_ITEM_COPY_CASE(S16,S16); + PS_METADATA_ITEM_COPY_CASE(S32,S32); + PS_METADATA_ITEM_COPY_CASE(S64,S64); + PS_METADATA_ITEM_COPY_CASE(U8,U8); + PS_METADATA_ITEM_COPY_CASE(U16,U16); + PS_METADATA_ITEM_COPY_CASE(U32,U32); + PS_METADATA_ITEM_COPY_CASE(U64,U64); + PS_METADATA_ITEM_COPY_CASE(F32,F32); + PS_METADATA_ITEM_COPY_CASE(F64,F64); + PS_METADATA_ITEM_COPY_CASE(STRING,V); // This will copy the string, not point at it. + case PS_DATA_VECTOR: { + PS_ASSERT_PTR_NON_NULL(in->data.V, NULL); + psVector *vecCopy = psVectorCopy(NULL, (psVector*)(in->data.V), + ((psVector*)(in->data.V))->type.type); + newItem = p_psMetadataItemAlloc(file, lineno, func, in->name, PS_DATA_VECTOR, in->comment, vecCopy); + psFree(vecCopy); // Drop reference + break; + } + case PS_DATA_TIME: { + psTime *timeCopy = NULL; + // pass through NULL values + if (in->data.V) { + timeCopy = psTimeCopy( (psTime*)(in->data.V) ); + if (timeCopy == NULL) { + psError(PS_ERR_BAD_PARAMETER_NULL, false, + "Error copying time. Time skipped.\n"); + } + } + newItem = p_psMetadataItemAlloc(file, lineno, func, in->name, PS_DATA_TIME, in->comment, timeCopy); + if (timeCopy) { + psFree(timeCopy); // Drop reference + } + break; + } + case PS_DATA_METADATA: { + // Metadata: copy the next level and stuff that in too + psMetadata *metadata = psMetadataCopy(NULL, in->data.md); + newItem = p_psMetadataItemAlloc(file, lineno, func, in->name, PS_DATA_METADATA, in->comment, metadata); + psFree(metadata); // Drop reference + break; + } + case PS_DATA_REGION: { + psRegion *region = in->data.V; // The region + psRegion *new = psRegionAlloc(region->x0, region->x1, region->y0, region->y1); // Copy of the region + newItem = p_psMetadataItemAlloc(file, lineno, func, in->name, PS_DATA_REGION, in->comment, new); + psFree(new); // Drop reference + break; + } + default: + // Other kinds of pointers + // XXX EAM : why was this a warning?? + // psWarning("Copying a pointer in the metadata item: %s (%x)\n", in->name, in->type); + newItem = p_psMetadataItemAlloc(file, lineno, func, in->name, in->type, in->comment, in->data.V); + break; + } + + return newItem; +} + +psMetadata *p_psMetadataCopy(const char *file, + unsigned int lineno, + const char *func, + psMetadata *out, + const psMetadata *in) +{ + PS_ASSERT_METADATA_NON_NULL(in, NULL); + + bool errorOccurred = false; + bool outAlloced = false; + if (out == NULL) { + out = p_psMetadataAlloc(file, lineno, func); + outAlloced = true; + } + psMetadataIterator *iter = psMetadataIteratorAlloc(in, PS_LIST_HEAD, NULL); + psMetadataItem *inItem = NULL; + while ((inItem = psMetadataGetAndIncrement(iter))) { + // Need to look for MULTI, which won't be picked up using the iterator. + psMetadataItem *multiCheckItem = psMetadataLookup(in, inItem->name); + unsigned int flag = PS_META_REPLACE; // Flag to indicate MULTI; otherwise, replace + if (multiCheckItem->type == PS_DATA_METADATA_MULTI) { + psTrace("psLib.types", 10, "MULTI: %s (%s)\n", inItem->name, inItem->comment); + flag = PS_META_DUPLICATE_OK; + } + psTrace("psLib.types", 5, "Copying %s (%s)...\n", inItem->name, inItem->comment); + + // Copy the item and add it on + psMetadataItem *newItem = psMetadataItemCopy(inItem); // Copied item + if (!psMetadataAddItem(out, newItem, PS_LIST_TAIL, flag)) { + psError(PS_ERR_UNKNOWN, false, "Error copying %s (%s) in the metadata\n", + inItem->name, inItem->comment); + if (outAlloced) { + psFree(out); + psFree(newItem); + psFree(iter); + return NULL; + } else { + errorOccurred = true; + } + } + psFree(newItem); // Drop reference + } + psFree(iter); + + if (errorOccurred) { + return NULL; + } + + return out; +} + +// this function copies the input metadata to the output, raising an error +// if any entries in 'in' a) do not exist in 'out' or b) have a different type +bool p_psMetadataUpdate(const char *file, + unsigned int lineno, + const char *func, + psMetadata *out, + const psMetadata *in) +{ + PS_ASSERT_METADATA_NON_NULL(in, NULL); + PS_ASSERT_METADATA_NON_NULL(out, NULL); + + bool result = true; + + // we loop over the metadata container by item, not by name. this pushes us directly into + // the elements of the MULTI items; we miss the containers. We need to check the + // (possible) MULTI containers to see how to set the flags. If we encounter a MULTI which + // is requesting REPLACE, then the first item on the input list from that MULTI should be + // added with the REPLACE flag turned on; the rest should have the REPLACE flag turned off. + // save a list of the MULTI entries which has REPLACE turned off? + + psArray *multiItems = psArrayAllocEmpty (128); + + psMetadataIterator *iter = psMetadataIteratorAlloc(in, PS_LIST_HEAD, NULL); + psMetadataItem *inItem = NULL; + while ((inItem = psMetadataGetAndIncrement(iter))) { + + // it is not possible to have RESET && UPDATE + // is it possible to have !RESET && !UPDATE? + // are these mutually exclusive conditions? + // input MULTI & RESET : replace an existing MULTI or ITEM of same name + // input MULTI & UPDATE : supplement an existing MULTI or ITEM of same name + // in both cases the name must exist in 'out' + + // Need to look for MULTI, which won't be picked up using the iterator. + psMetadataItem *multiCheckItem = psMetadataLookup(in, inItem->name); + + // default operation for the new metadata item (replace existing entry, require existence & type) + unsigned int flag = PS_META_REPLACE | PS_META_REQUIRE_ENTRY | PS_META_REQUIRE_TYPE; + + // for MULTI items, the mode is carried on the multi container + if (PS_METADATA_ITEM_GET_TYPE(multiCheckItem) == PS_DATA_METADATA_MULTI) { + psTrace("psLib.types", 10, "MULTI: %s (%s)\n", inItem->name, inItem->comment); + if (multiCheckItem->type & PS_META_UPDATE_FOLDER) { + psTrace("psLib.types", 10, "supplement MULTI with new entries\n"); + flag = PS_META_DUPLICATE_OK | PS_META_REQUIRE_ENTRY | PS_META_REQUIRE_TYPE; + } else { + multiCheckItem->type |= PS_META_UPDATE_FOLDER; + // save this multi so we can reset their flags below + psArrayAdd (multiItems, 128, multiCheckItem); + } + } + psTrace("psLib.types", 5, "Copying %s (%s)...\n", inItem->name, inItem->comment); + + // for METADATA folders, do something different? + if (PS_METADATA_ITEM_GET_TYPE(inItem) == PS_DATA_METADATA) { + if (inItem->type & PS_META_UPDATE_FOLDER) { + flag = PS_META_UPDATE_FOLDER | PS_META_REQUIRE_ENTRY | PS_META_REQUIRE_TYPE; + } + } + + // Copy the item and add it on. report all errors so we get a listing + psMetadataItem *newItem = psMetadataItemCopy(inItem); // Copied item + if (!psMetadataAddItem(out, newItem, PS_LIST_TAIL, flag)) { + psError(PS_ERR_UNKNOWN, false, "Error copying %s\n", inItem->name); + result = false; + } + psFree(newItem); // Drop reference + } + psFree(iter); + + // remove the UPDATE_FOLDER flag from the following MULTI items + for (int i = 0; i < multiItems->n; i++) { + psMetadataItem *item = multiItems->data[i]; + item->type &= ~PS_META_UPDATE_FOLDER; + } + psFree (multiItems); + + if (!result) { + psError(PS_ERR_UNKNOWN, false, "failed to update metadata\n"); + } + + return result; +} + +// this function copies the input metadata to the output, supplementing existing metadata +// folders with the contents from corresponding folders in the output +bool p_psMetadataOverlay(const char *file, + unsigned int lineno, + const char *func, + psMetadata *out, + const psMetadata *in) +{ + PS_ASSERT_METADATA_NON_NULL(in, NULL); + PS_ASSERT_METADATA_NON_NULL(out, NULL); + + bool result = true; + + psMetadataIterator *iter = psMetadataIteratorAlloc(in, PS_LIST_HEAD, NULL); + psMetadataItem *inItem = NULL; + while ((inItem = psMetadataGetAndIncrement(iter))) { + // Need to look for MULTI, which won't be picked up using the iterator. + psMetadataItem *multiCheckItem = psMetadataLookup(in, inItem->name); + unsigned int flag = PS_META_REPLACE; // Flag to indicate MULTI; otherwise, replace + if (multiCheckItem->type == PS_DATA_METADATA_MULTI) { + psTrace("psLib.types", 10, "MULTI: %s (%s)\n", inItem->name, inItem->comment); + flag = PS_META_DUPLICATE_OK; + } + psTrace("psLib.types", 5, "Copying %s (%s)...\n", inItem->name, inItem->comment); + + // if this is a metadata, and it has a corresponding match, overlay them + if (inItem->type == PS_DATA_METADATA) { + bool status; + psMetadata *outFolder = psMetadataLookupMetadata (&status, out, inItem->name); + if (outFolder) { + if (!psMetadataOverlay (outFolder, inItem->data.md)) { + fprintf (stderr, "Error overlaying metadata folder\n"); + result = false; + } + continue; + } + } + + // Copy the item and add it on. report all errors so we get a listing + psMetadataItem *newItem = psMetadataItemCopy(inItem); // Copied item + if (!psMetadataAddItem(out, newItem, PS_LIST_TAIL, flag)) { + fprintf (stderr, "Error copying %s\n", inItem->name); + result = false; + } + psFree(newItem); // Drop reference + } + psFree(iter); + + if (!result) { + psError(PS_ERR_UNKNOWN, false, "failed to update metadata\n"); + } + return result; +} + +// may need to extend this to change the keyname in the copy +bool psMetadataItemSupplement(psMetadata *out, + const psMetadata *in, + const char *key) +{ + PS_ASSERT_METADATA_NON_NULL(in, false); + PS_ASSERT_METADATA_NON_NULL(out, false); + PS_ASSERT_STRING_NON_EMPTY(key, false); + + psMetadataItem *item = psMetadataLookup(in, key); + if (!item) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Could not find '%s' in metadata.\n", key); + return false; + } + if (!psMetadataAddItem(out, item, PS_LIST_TAIL, PS_META_REPLACE) ) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Could not add %s to metadata.\n", key); + return false; + } + + return true; +} + +// XXX is it sensible that item is a 'const' here? +bool psMetadataAddItem(psMetadata *md, + const psMetadataItem *item, + int location, + psS32 flags) +{ + PS_ASSERT_METADATA_NON_NULL(md,false); + PS_ASSERT_METADATA_ITEM_NON_NULL(item,false); + + psHash *mdTable = md->hash; + psList *mdList = md->list; + char *key = item->name; + + // See if key is already in table + psMetadataItem *existingEntry = psHashLookup(mdTable, key); + + // this block handles cases for the MULTI items + if (item->type == PS_DATA_METADATA_MULTI) { + // the incoming entry is PS_DATA_METADATA_MULTI + + // Shouldn't have a second reference to the same MULTI in a single Metadata! + // XXX not sure I understand this case + if (item == existingEntry) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Cannot have 2 references to the same MULTI in a single Metadata!"); + return false; + } + + if (flags & PS_META_REPLACE) { + // drop the existing entry or entries + psMetadataRemoveKey(md, key); + } else { + // elevate the existing hash entry to be PS_DATA_METADATA_MULTI + existingEntry = makeMetaMulti(mdTable,key,existingEntry); + } + + // add all the items in the incoming entry to metadata + psList* list = item->data.list; + if (list != NULL) { + psListIterator* iter = psListIteratorAlloc(list,PS_LIST_HEAD,true); + psMetadataItem* listItem; + while ((listItem=(psMetadataItem*)psListGetAndIncrement(iter)) != NULL) { + psMetadataAddItem(md,listItem,location,flags); + } + psFree(iter); + } + + return true; // all done. + } + + // special block for items which are METADATA folders + if (existingEntry && (item->type == PS_DATA_METADATA)) { + // REPLACE and UPDATE must be mutually exclusive + psAssert ((!((flags & PS_META_REPLACE) && (flags & PS_META_UPDATE_FOLDER))), "cannot have both REPLACE and UPDATE"); + + # if (0) + // handle the case of replace below + if (flags & PS_META_REPLACE) { + // drop the existing entry (skip if we are replacing with same pointer) + // XXX what if existingEntry is a MULTI? drop all? + // XXX this segment does not check for matched types + if (item != existingEntry) { + psMetadataRemoveKey(md, key); + } + existingEntry = NULL; + } + # endif + + if (flags & PS_META_UPDATE_FOLDER) { + if (existingEntry->type != PS_DATA_METADATA) { + psError(PS_ERR_UNKNOWN, false, "invalid to request UPDATE for metadata which matches another type"); + return false; + } + + // merge the existing entry : this completes the insert + if (!psMetadataCopy ((psMetadata *)existingEntry->data.V, (psMetadata *)item->data.V)) { + psError(PS_ERR_UNKNOWN, false, "failed to copy new metadata on existing"); + return false; + } + return true; + } + } + + // how the item is added to the hash depends on prior existence, flags, etc. + // XXX i think these cases are overloaded - are all combinations possible? + if (existingEntry) { // prior existence + + // duplicate entries allowed - add another entry. + if ((existingEntry->type == PS_DATA_METADATA_MULTI) || (flags & PS_META_DUPLICATE_OK)) { + + // make sure the existing entry is PS_DATA_METADATA_MULTI + existingEntry = makeMetaMulti(mdTable,key,existingEntry); + + // add item to the existing hash's list of duplicate entries + if (!psListAdd(existingEntry->data.list, PS_LIST_TAIL, (psMetadataItem *) item) ) { + psError(PS_ERR_UNKNOWN, false, _("Failed to add metadata item, %s, to metadata collection list."), key); + return false; + } + // add to the metadata list of entries + if (!psListAdd(mdList, location, (psMetadataItem *) item)) { + psError(PS_ERR_UNKNOWN, false, _("Failed to add metadata item, %s, to metadata collection list."), key); + return false; + } + return true; + } + + // replace entry instead of creating a duplicate entry. + if (flags & PS_META_REPLACE) { + + if ((flags & PS_META_REQUIRE_TYPE) && (existingEntry->type != item->type)) { + psError (PS_ERR_UNKNOWN, true, _("Existing item %s does not match type (%x) for item requiring matching type (%x)"), key, existingEntry->type, item->type); + return false; + } + if (item != existingEntry) { + // Replacing an item with itself can lead to memory corruption + // when you blow away what you're trying to add + psMetadataRemoveKey(md, key); + // add to the metadata has of entries + if (!psHashAdd(mdTable, key, (psMetadataItem *) item)) { + psError(PS_ERR_UNKNOWN, false, _("Failed to add metadata item, %s, to metadata collection list."), key); + return false; + } + // add to the metadata list of entries + if (!psListAdd(mdList, location, (psMetadataItem *) item)) { + psError(PS_ERR_UNKNOWN, false, _("Failed to add metadata item, %s, to metadata collection list."), key); + return false; + } + return true; + } else { + // adding myself back to the list; this is a NOP + return true; + } + } + + // if specified, keep the existing entry + if (flags & PS_META_NO_REPLACE) { + return true; + } + + // default is to error on duplicate entry. + psError(PS_ERR_BAD_PARAMETER_VALUE, true, _("Duplicate metadata item name: %s is not allowed. Use a psMetadataFlags option to allow such action."), item->name); + return false; + } + + // OK, this is a new item. + if (flags & PS_META_REQUIRE_ENTRY) { + psError (PS_ERR_UNKNOWN, true, _("No matching item found for item requiring existing entry (%s)"), key); + return false; + } + + // Node doesn't exist - Add new metadata item to metadata collection's hash + psHashAdd(mdTable, key, (psMetadataItem *) item); + + // Create a multi, if requested + if (flags & PS_META_DUPLICATE_OK) { + makeMetaMulti(mdTable, key, (psMetadataItem *) item); // Casting away const! + } + + if (!psListAdd(mdList, location, (psPtr)item)) { + psError(PS_ERR_UNKNOWN, false, _("Failed to add metadata item, %s, to metadata collection list."), key); + return false; + } + + return true; +} + +bool psMetadataAdd(psMetadata *md, + long location, + const char *name, + int format, + const char *comment,...) +{ + PS_ASSERT_METADATA_NON_NULL(md, false); + + va_list argPtr; + va_start(argPtr, comment); + bool result = psMetadataAddV(md,location,name,format,comment,argPtr); + va_end(argPtr); + + return result; +} + +bool psMetadataAddV(psMetadata *md, + long location, + const char *name, + int format, + const char *comment, + va_list list) +{ + PS_ASSERT_METADATA_NON_NULL(md,false); + + psMetadataItem* metadataItem = psMetadataItemAllocV(name, format & PS_METADATA_TYPE_MASK, comment, list); + + if (!psMetadataAddItem(md, metadataItem, location, format & PS_METADATA_FLAGS_MASK)) { + psError(PS_ERR_UNKNOWN,false,_("Failed to add metadata item to metadata collection list.")); + psFree(metadataItem); + return false; + } + // Decrement reference count, since the metadata item is now in metadata collection and no longer needed + psFree(metadataItem); + + return true; +} + +#define METADATA_ADD_TYPE(NAME,TYPE,METATYPE) \ +bool psMetadataAdd##NAME(psMetadata* md, long where, const char* name, \ + int format, const char* comment, TYPE value) { \ + return psMetadataAdd(md, where, name, format | METATYPE, comment, value); \ +} + +METADATA_ADD_TYPE(Bool,bool,PS_DATA_BOOL) +METADATA_ADD_TYPE(S8,psS8,PS_DATA_S8) +METADATA_ADD_TYPE(S16,psS16,PS_DATA_S16) +METADATA_ADD_TYPE(S32,psS32,PS_DATA_S32) +METADATA_ADD_TYPE(S64,psS64,PS_DATA_S64) +METADATA_ADD_TYPE(U8,psU8,PS_DATA_U8) +METADATA_ADD_TYPE(U16,psU16,PS_DATA_U16) +METADATA_ADD_TYPE(U32,psU32,PS_DATA_U32) +METADATA_ADD_TYPE(U64,psU64,PS_DATA_U64) +METADATA_ADD_TYPE(F32,psF32,PS_DATA_F32) +METADATA_ADD_TYPE(F64,psF64,PS_DATA_F64) +METADATA_ADD_TYPE(List,psList*,PS_DATA_LIST) +METADATA_ADD_TYPE(Str,const char*,PS_DATA_STRING) +METADATA_ADD_TYPE(Vector,psVector*,PS_DATA_VECTOR) +METADATA_ADD_TYPE(Image,psImage*,PS_DATA_IMAGE) +METADATA_ADD_TYPE(Hash,psHash*,PS_DATA_HASH) +METADATA_ADD_TYPE(LookupTable,psLookupTable*,PS_DATA_LOOKUPTABLE) +METADATA_ADD_TYPE(Metadata,psMetadata*,PS_DATA_METADATA) +METADATA_ADD_TYPE(Array,psArray*,PS_DATA_ARRAY) +METADATA_ADD_TYPE(Time,psTime*,PS_DATA_TIME) +METADATA_ADD_TYPE(Unknown,psPtr,PS_DATA_UNKNOWN) + +bool psMetadataAddPtr(psMetadata* md, long where, const char* name, + psDataType type, const char* comment, psPtr value) +{ + return psMetadataAdd(md, where, name, type, comment, value); +} + + + +// Remove by key name +bool psMetadataRemoveKey(psMetadata *md, + const char *key) +{ + PS_ASSERT_METADATA_NON_NULL(md, false); + PS_ASSERT_STRING_NON_EMPTY(key, false); + + psList* mdList = md->list; + psHash* mdTable = md->hash; + + psMetadataItem* entry = psHashLookup(mdTable,key); + if (entry == NULL) { + psError(PS_ERR_UNKNOWN, false, + _("Failed to remove metadata item, %s, from metadata table."), key); + return false; + } + if (entry->type == PS_DATA_METADATA_MULTI) { + psMetadataItem* listItem; + psListIterator* iter = psListIteratorAlloc( + entry->data.list, + PS_LIST_HEAD,true); + while ((listItem=psListGetAndIncrement(iter)) != NULL) { + psListRemoveData(mdList, listItem); + } + psFree(iter); + psHashRemove(mdTable,key); + } else { + psListRemoveData(mdList, entry); + psHashRemove(mdTable, key); + } + + return true; +} + + +// Remove by index +bool psMetadataRemoveIndex(psMetadata *md, + long location) +{ + PS_ASSERT_METADATA_NON_NULL(md, false); + + psList* mdList = md->list; + psHash* mdTable = md->hash; + + psMetadataItem* entry = psListGet(mdList, location); + if (entry == NULL) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, + _("Could not find metadata item at index %ld."), location); + return false; + } + const char *key = entry->name; + + if (key == NULL) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, + _("Failed to remove metadata item, at index %ld, from metadata list."), + location); + return false; + } + + psMetadataItem* tableItem = psHashLookup(mdTable, key); + if (tableItem == NULL) { + psError(PS_ERR_UNKNOWN, false, + _("Failed to remove metadata item, %s, from metadata table."), key); + return false; + } + + if (tableItem->type == PS_DATA_METADATA_MULTI) { + // multiple entries with same key, remove just the specified one + psListRemoveData(tableItem->data.list, entry); + } else { + //Tested below. psHashRemove can't return false here. + psHashRemove(mdTable, key); + } + psListRemove(mdList, location); + + return true; +} + +psMetadataItem* psMetadataLookup(const psMetadata *md, + const char *key) +{ + PS_ASSERT_METADATA_NON_NULL(md,NULL); + PS_ASSERT_STRING_NON_EMPTY(key,NULL); + + return (psMetadataItem*)psHashLookup(md->hash, key); +} + +void* psMetadataLookupPtr(bool *status, + const psMetadata *md, + const char *key) +{ + PS_ASSERT_METADATA_NON_NULL(md,NULL); + + if (status) { + *status = true; + } + + psMetadataItem *metadataItem = psMetadataLookup(md, key); + if (metadataItem == NULL) { + if (status) { + *status = false; + } + return NULL; + } + if (metadataItem->type == PS_DATA_METADATA_MULTI) { + // if multiple keys found, use the first. + // metadataItem = (psMetadataItem*)((metadataItem->data.list)->head); + metadataItem = (psMetadataItem*)(metadataItem->data.list->head->data); + if (status) { + *status = true; + } + } + + if(PS_DATA_IS_PRIMITIVE(metadataItem->type)) { + if (status) { + *status = false; + } + psError(PS_ERR_BAD_PARAMETER_VALUE, true, + _("Specified psDataType, %d, is not supported."), + metadataItem->type); + return NULL; + } else { + return metadataItem->data.V; + } +} + +#define psMetadataLookupNumTYPE(TYPE) \ +ps##TYPE psMetadataLookup##TYPE(bool *status, const psMetadata *md, const char *key) \ +{ \ + psMetadataItem *metadataItem = NULL; \ + ps##TYPE value = 0; \ + \ + if (status) { \ + *status = true; \ + } \ + \ + metadataItem = psMetadataLookup(md, key); \ + if(metadataItem == NULL) { \ + if (status) { \ + *status = false; \ + } \ + return 0; \ + } \ + if (metadataItem->type == PS_DATA_METADATA_MULTI) { \ + /* if multiple keys found, use the first. */ \ + metadataItem = (psMetadataItem*)((metadataItem->data.list)->head); \ + } \ + \ + switch (metadataItem->type) { \ + case PS_DATA_S8: \ + value = (ps##TYPE)metadataItem->data.S8; \ + break; \ + case PS_DATA_S16: \ + value = (ps##TYPE)metadataItem->data.S16; \ + break; \ + case PS_DATA_S32: \ + value = (ps##TYPE)metadataItem->data.S32; \ + break; \ + case PS_DATA_S64: \ + value = (ps##TYPE)metadataItem->data.S64; \ + break; \ + case PS_DATA_U8: \ + value = (ps##TYPE)metadataItem->data.U8; \ + break; \ + case PS_DATA_U16: \ + value = (ps##TYPE)metadataItem->data.U16; \ + break; \ + case PS_DATA_U32: \ + value = (ps##TYPE)metadataItem->data.U32; \ + break; \ + case PS_DATA_U64: \ + value = (ps##TYPE)metadataItem->data.U64; \ + break; \ + case PS_DATA_F32: \ + value = (ps##TYPE)metadataItem->data.F32; \ + break; \ + case PS_DATA_F64: \ + value = (ps##TYPE)metadataItem->data.F64; \ + break; \ + case PS_DATA_BOOL: \ + if (metadataItem->data.B) { \ + value = 1; \ + } \ + break; \ + default: \ + /* if you get to this point, the value is not a number. */ \ + if (status) { \ + *status = false; \ + } \ + break; \ + } \ + \ + /* psFree(metadataItem); currently, the lookup doesn't increment the ref count */ \ + return value; \ +} + +psMetadataLookupNumTYPE(F32) +psMetadataLookupNumTYPE(F64) +psMetadataLookupNumTYPE(S8) +psMetadataLookupNumTYPE(S16) +psMetadataLookupNumTYPE(S32) +psMetadataLookupNumTYPE(S64) +psMetadataLookupNumTYPE(U8) +psMetadataLookupNumTYPE(U16) +psMetadataLookupNumTYPE(U32) +psMetadataLookupNumTYPE(U64) +psMetadataLookupNumTYPE(Bool) + +psMetadataItem* psMetadataGet(const psMetadata *md, + int location) +{ + psMetadataItem* entry = NULL; + + PS_ASSERT_METADATA_NON_NULL(md,NULL); + + entry = (psMetadataItem*) psListGet(md->list, location); + if (entry == NULL) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, _("Could not find metadata item at index %d."), location); + return NULL; + } + + return entry; +} + +// XXX should md be const? +psMetadataIterator* p_psMetadataIteratorAlloc(const char *file, + unsigned int lineno, + const char *func, + const psMetadata* md, + long location, + const char* regex) +{ + PS_ASSERT_METADATA_NON_NULL(md,NULL); + + psMetadataIterator* newIter = p_psAlloc(file, lineno, func, sizeof(psMetadataIterator)); + psMemSetDeallocator(newIter, (psFreeFunc) metadataIteratorFree); + + newIter->iter = psListIteratorAlloc(md->list, location, false); + + if (regex) { + newIter->regex = psAlloc(sizeof(regex_t)); + int regRtn = regcomp(newIter->regex, regex, 0); + if (regRtn != 0) { + char errMsg[256]; + regerror(regRtn, newIter->regex, errMsg, 256); + psError(PS_ERR_BAD_PARAMETER_VALUE, true, + _("Specified regular expression is invalid. %s."), + errMsg); + psFree(newIter); + return NULL; + } + if (!psMetadataIteratorSet(newIter, location)) { + psError(PS_ERR_UNKNOWN, false, "Unable to set location=%ld for metadata iterator.\n", location); + psFree(newIter); + return NULL; + } + } else { + newIter->regex = NULL; + } + + return newIter; +} + +bool psMetadataIteratorSet(psMetadataIterator* iterator, + long location) +{ + PS_ASSERT_METADATA_ITERATOR_NON_NULL(iterator, NULL); + psListIterator *listIter = iterator->iter; + + regex_t *regex = iterator->regex; + + // handle trivial case where no regex subsetting is required. + if (regex == NULL) { + return psListIteratorSet(listIter, location); + } + + // If there's a regex, then we need to count into the list. + // We count ONLY those entries that match the regex + + // Here we count in from the tail + if (location < 0) { + psListIteratorSet(listIter, PS_LIST_TAIL); + psMetadataItem *item; // Item from iteration + int match = 0; // Match number + while (listIter->cursor && (item = listIter->cursor->data)) { + if (regexec(regex, item->name, 0, NULL, 0) == 0) { + // this key is a match + match--; + if (match == location) { + break; + } + } + (void)psListGetAndDecrement(listIter); + } + + // We return "true" even if we didn't find the Nth match, because the iterator is set appropriately + // (off the front), and subsequent calls to psMetadataGetAnd{In,De}crement will return NULL. + return true; + } + + // Here we count in from the head + psListIteratorSet(listIter, PS_LIST_HEAD); + psMetadataItem *item; // Item from iteration + int match = -1; // Match number + while (listIter->cursor && (item = listIter->cursor->data)) { + if (regexec(regex, item->name, 0, NULL, 0) == 0) { + // this key is a match + match++; + if (match == location) { + break; + } + } + (void)psListGetAndIncrement(listIter); + } + // We return "true" even if we didn't find the Nth match, because the iterator is set appropriately + // (off the end), and subsequent calls to psMetadataGetAnd{In,De}crement will return NULL. + return true; +} + +psMetadataItem* psMetadataGetAndIncrement(psMetadataIterator* iterator) +{ + PS_ASSERT_METADATA_ITERATOR_NON_NULL(iterator,NULL); + psListIterator* listIter = iterator->iter; + regex_t* regex = iterator->regex; + + // handle trivial case where no regex subsetting is required. + if (regex == NULL) { + return (psMetadataItem*)psListGetAndIncrement(listIter); + } + + // Iterate until we find something matching the regex + psMetadataItem *newItem; // New MD item from iteration + while ((newItem = psListGetAndIncrement(listIter))) { + if (regexec(regex, newItem->name, 0, NULL, 0) == 0) { + // this key is a match + break; + } + } + return newItem; +} + +psMetadataItem* psMetadataGetAndDecrement(psMetadataIterator* iterator) +{ + PS_ASSERT_METADATA_ITERATOR_NON_NULL(iterator, NULL); + psListIterator* listIter = iterator->iter; + regex_t* regex = iterator->regex; + + // handle trivial case where no regex subsetting is required. + if (regex == NULL) { + return (psMetadataItem*)psListGetAndDecrement(listIter); + } + + // Iterate until we find something matching the regex + psMetadataItem *newItem; // New MD item from iteration + while ((newItem = psListGetAndDecrement(listIter))) { + if (regexec(regex, newItem->name, 0, NULL, 0) == 0) { + // this key is a match + break; + } + } + return newItem; +} + +psMetadata *psMetadataLookupMetadata(bool *status, + const psMetadata *md, + const char *key) +{ + PS_ASSERT_METADATA_NON_NULL(md,NULL); + psMetadataItem *item = psMetadataLookup((psMetadata*)md, key); // The metadata with instruments + psMetadata *value = NULL; // The value to return + if (!item) { + // The given key isn't in the metadata + if (status) { + *status = false; + } else { + psError(PS_ERR_IO, true, "Couldn't find %s in the metadata.\n", key); + } + } else if (item->type != PS_DATA_METADATA) { + // The value at the key isn't metadata + if (status) { + *status = false; + } else { + psLogMsg(__func__, PS_LOG_DETAIL, "%s isn't of type PS_DATA_META, as expected.\n", key); + } + // value = NULL; + } else { + // We have the requested metadata + if (status) { + *status = true; + } + value = item->data.md; // The requested metadata + } + return value; +} + +psTime *psMetadataLookupTime(bool *status, + const psMetadata *md, + const char *key) +{ + PS_ASSERT_METADATA_NON_NULL(md,NULL); + + psMetadataItem *item = psMetadataLookup((psMetadata*)md, key); + if (!item) { + // The given key isn't in the metadata + if (status) { + *status = false; + } else { + psError(PS_ERR_IO, true, "Couldn't find %s in the metadata.\n", key); + } + return NULL; + } + + if (item->type != PS_DATA_TIME) { + // The value at the key isn't metadata + if (status) { + *status = false; + } else { + psLogMsg(__func__, PS_LOG_DETAIL, "%s isn't of type PS_DATA_TIME, as expected.\n", key); + } + return NULL; + } + + // We have the requested metadata + if (status) { + *status = true; + } + + return item->data.V; +} + + +psString psMetadataLookupStr(bool *status, + const psMetadata *md, + const char *key) +{ + PS_ASSERT_METADATA_NON_NULL(md,NULL); + + psMetadataItem *item = psMetadataLookup((psMetadata*)md, key); // The metadata with instruments + // char *value = NULL; // The value to return + psString value = NULL; + if (!item) { + // The given key isn't in the metadata + if (status) { + *status = false; + } else { + psError(PS_ERR_IO, true, "Couldn't find %s in the metadata.\n", key); + } + } else if (item->type != PS_DATA_STRING) { + // The value at the key isn't of the desired type + if (status) { + *status = false; + } else { + psLogMsg(__func__, PS_LOG_DETAIL, "%s isn't of type PS_DATA_STRING, as expected.\n", key); + } + // value = NULL; + } else { + // We have the requested metadata + if (status) { + *status = true; + } + value = item->data.V; + } + return value; +} + +psList *psMetadataKeys(psMetadata *md) +{ + PS_ASSERT_METADATA_NON_NULL(md, NULL); + return psHashKeyList(md->hash); +} + + +bool psMetadataPrint(FILE *fd, + const psMetadata *md, + int level) +{ + PS_ASSERT_METADATA_NON_NULL(md,false); + if (fd == NULL) { + fd = stdout; + } else { + if ( fprintf(fd, "\n") < 0 ) { + psError(PS_ERR_IO, true, + "Invalid file pointer in psMetadataPrint. Could not write to fd.\n"); + return false; + } + } + + psErrorClear(); // we're going to _append_ errors to the stack + bool noErrors = true; // have see seen any errors? + + // Casting away const --- the addition of an iterator should not be considered an invasion of "const". + psMetadataIterator *iter = psMetadataIteratorAlloc(md, PS_LIST_HEAD, NULL); // Iterator + psMetadataItem *item = NULL; // Item from metadata + while ( (item = psMetadataGetAndIncrement(iter)) ) { + //So that unsupported types are not printed here, adding a switch statement + switch(item->type) { + case PS_DATA_STRING: + case PS_DATA_BOOL: + case PS_DATA_S8: + case PS_DATA_S16: + case PS_DATA_S32: + case PS_DATA_S64: + case PS_DATA_U8: + case PS_DATA_U16: + case PS_DATA_U32: + case PS_DATA_U64: + case PS_DATA_F32: + case PS_DATA_F64: + case PS_DATA_METADATA: + case PS_DATA_REGION: + case PS_DATA_LIST: + case PS_DATA_TIME: + case PS_DATA_VECTOR: + // Indent... + // fprintf(fd, "\n"); + for (int i = 0; i < level; i++) { + fprintf(fd, " "); + } + fprintf(fd, "%s", item->name); + if (item->comment && strlen(item->comment) > 0) { + fprintf(fd, " (%s)", item->comment); + } + fprintf(fd, ": "); + break; + case PS_DATA_UNKNOWN: + default: + break; + } + switch (item->type) { + case PS_DATA_STRING: + fprintf(fd, "%s\n", item->data.str); + break; + case PS_DATA_BOOL: + if (item->data.B) { + fprintf(fd, "True\n"); + } else { + fprintf(fd, "False\n"); + } + break; + case PS_DATA_S8: + fprintf(fd, "%d\n", item->data.S8); + break; + case PS_DATA_S16: + fprintf(fd, "%d\n", item->data.S16); + break; + case PS_DATA_S32: + fprintf(fd, "%d\n", item->data.S32); + break; + case PS_DATA_S64: + fprintf(fd, "%jd\n", item->data.S64); + break; + case PS_DATA_U8: + fprintf(fd, "%u\n", item->data.U8); + break; + case PS_DATA_U16: + fprintf(fd, "%u\n", item->data.U16); + break; + case PS_DATA_U32: + fprintf(fd, "%u\n", item->data.U32); + break; + case PS_DATA_U64: + fprintf(fd, "%ju\n", item->data.U64); + break; + case PS_DATA_F32: + fprintf(fd, "%f\n", item->data.F32); + break; + case PS_DATA_F64: + fprintf(fd, "%f\n", item->data.F64); + break; + case PS_DATA_METADATA: + if (!item->data.V) { + fprintf(fd, "(null)\n"); + break; + } + psMetadataPrint(fd, item->data.V, level + 1); + for (int i = 0; i < level; i++) { + fprintf(fd, " "); + } + fprintf(fd, "%s -- END\n", item->name); + break; + case PS_DATA_REGION: + if (!item->data.V) { + fprintf(fd, "(null)\n"); + break; + } + psString region = psRegionToString(*(psRegion*)item->data.V); + fprintf(fd, "%s\n", region); + psFree(region); + break; + case PS_DATA_LIST: + if (!item->data.V) { + fprintf(fd, "(null)\n"); + break; + } + fprintf(fd, "\n", ((psList*)item->data.V)->n); + break; + case PS_DATA_TIME: + if (!item->data.V) { + fprintf(fd, "(null)\n"); + break; + } + psString time = psTimeToISO(item->data.V); + fprintf(fd, "%s\n", time); + psFree(time); + break; + case PS_DATA_VECTOR: + if (!item->data.V) { + fprintf(fd, "(null)\n"); + break; + } + psVector *vector = item->data.V; + switch (vector->type.type) { + case PS_DATA_U8: + fprintf(fd, "U8 "); + for (int i = 0; i < vector->n; i++) { + fprintf(fd, "%u ", vector->data.U8[i]); + } + fprintf(fd, "\n"); + break; + case PS_DATA_U16: + fprintf(fd, "U16 "); + for (int i = 0; i < vector->n; i++) { + fprintf(fd, "%u ", vector->data.U16[i]); + } + fprintf(fd, "\n"); + break; + case PS_DATA_U32: + fprintf(fd, "U32 "); + for (int i = 0; i < vector->n; i++) { + fprintf(fd, "%u ", vector->data.U32[i]); + } + fprintf(fd, "\n"); + break; + case PS_DATA_U64: + fprintf(fd, "U64 "); + for (int i = 0; i < vector->n; i++) { + fprintf(fd, "%ju ", vector->data.U64[i]); + } + fprintf(fd, "\n"); + break; + case PS_DATA_S8: + fprintf(fd, "S8 "); + for (int i = 0; i < vector->n; i++) { + fprintf(fd, "%d ", vector->data.S8[i]); + } + fprintf(fd, "\n"); + break; + case PS_DATA_S16: + fprintf(fd, "S16 "); + for (int i = 0; i < vector->n; i++) { + fprintf(fd, "%d ", vector->data.S16[i]); + } + fprintf(fd, "\n"); + break; + case PS_DATA_S32: + fprintf(fd, "S32 "); + for (int i = 0; i < vector->n; i++) { + fprintf(fd, "%d ", vector->data.S32[i]); + } + fprintf(fd, "\n"); + break; + case PS_DATA_S64: + fprintf(fd, "S64 "); + for (int i = 0; i < vector->n; i++) { + fprintf(fd, "%jd ", vector->data.S64[i]); + } + fprintf(fd, "\n"); + break; + case PS_DATA_F32: + fprintf(fd, "F32 "); + for (int i = 0; i < vector->n; i++) { + fprintf(fd, "%f ", vector->data.F32[i]); + } + fprintf(fd, "\n"); + break; + case PS_DATA_F64: + fprintf(fd, "F64 "); + for (int i = 0; i < vector->n; i++) { + fprintf(fd, "%f ", vector->data.F64[i]); + } + fprintf(fd, "\n"); + break; + case PS_DATA_UNKNOWN: + default: + psError(PS_ERR_BAD_PARAMETER_VALUE, true, + _("Specified psDataType, %d, is not supported."), vector->type.type); + psFree(iter); + fprintf(fd, "\n"); + return false; + } + break; + default: + psError(PS_ERR_IO, false, "Non-printable metadata type for \"%s\": %x\n", + item->name, item->type); + noErrors = false; + break; + } + } + psFree(iter); + + return noErrors; +}