Index: /trunk/ppMerge/src/Makefile.am =================================================================== --- /trunk/ppMerge/src/Makefile.am (revision 23409) +++ /trunk/ppMerge/src/Makefile.am (revision 23410) @@ -21,5 +21,5 @@ ppMergeFileGroup.c \ ppMergeReadChunk.c \ - ppMergeLoop_Threaded.c \ + ppMergeLoop.c \ ppMergeSetThreads.c \ ppMergeMask.c \ Index: /trunk/ppMerge/src/ppMergeLoop.c =================================================================== --- /trunk/ppMerge/src/ppMergeLoop.c (revision 23409) +++ /trunk/ppMerge/src/ppMergeLoop.c (revision 23410) @@ -1,3 +1,3 @@ -/** @file ppMergeLoop.c +/** @file ppMergeLoop_Threaded.c * * @brief @@ -6,10 +6,24 @@ * * @author IfA - * @version $Revision: 1.14 $ $Name: not supported by cvs2svn $ - * @date $Date: 2009-02-01 21:43:05 $ + * @version $Revision: 1.9 $ $Name: not supported by cvs2svn $ + * @date $Date: 2009-02-06 02:44:31 $ * Copyright 2009 Institute for Astronomy, University of Hawaii */ +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include + +#include +#include +#include + #include "ppMerge.h" + +// XXX this function is now sufficiently different for the major types, it would make sense to just +// split it into three: BASIC, SHUTTER, DARK bool ppMergeLoop(pmConfig *config) @@ -21,38 +35,39 @@ 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 + bool haveMasks = psMetadataLookupBool(&mdok, arguments, "INPUTS.MASKS"); // Do we have masks? + bool haveVariances = psMetadataLookupBool(&mdok, arguments, "INPUTS.VARIANCES"); // Do we have variances? + + psArray *inputs = ppMergeFileDataLevel(config, "PPMERGE.INPUT", PM_FPA_LEVEL_READOUT); // Input images + psArray *masks = NULL, *variances = NULL; // Input masks and variances 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 + if (haveVariances) { + variances = ppMergeFileDataLevel(config, "PPMERGE.INPUT.VARIANCE", PM_FPA_LEVEL_READOUT); + } + + int nThreads = psMetadataLookupS32 (&mdok, arguments, "NTHREADS"); + if (!mdok) nThreads = 0; + + // General combination parameters + 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 useVariances = psMetadataLookupBool(NULL, arguments, "VARIANCES"); // Use variances? + + // 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 - psImageMaskType maskVal; - psImageMaskType markVal; - if (!pmConfigMaskSetBits (&maskVal, &markVal, config)) { - psError (PS_ERR_UNKNOWN, true, "Unable to define the mask bit values"); - return false; + psImageMaskType maskVal, markVal; + if (!pmConfigMaskSetBits(&maskVal, &markVal, config)) { + psError (PS_ERR_UNKNOWN, true, "Unable to define the mask bit values"); + return false; } @@ -65,5 +80,5 @@ combination->iter = iter; combination->rej = rej; - combination->weights = useWeights; + combination->variances = useVariances; psMetadata *stats = NULL; ///< Statistics for output @@ -73,5 +88,5 @@ } - pmFPAview *view = pmFPAviewAlloc(0); ///< View to component of interest + pmFPAview *view = pmFPAviewAlloc(0); // View to component of interest // Retrieve data placed on analysis @@ -122,4 +137,5 @@ assert(output && output->fpa); pmFPA *outFPA = output->fpa; ///< Output FPA + pmHDU *lastHDU = NULL; // Last HDU that was updated int cellNum = 0; ///< Index of cell if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) { @@ -143,106 +159,179 @@ } + // Update the header + { + pmHDU *hdu = pmHDUGetHighest(outFPA, outChip, outCell); // HDU for file + if (hdu && hdu != lastHDU) { + ppMergeVersionHeader(hdu->header); + lastHDU = hdu; + } + } + + float shutterRef = NAN; ///< Reference shutter correction + pmReadout *pattern = NULL; + if (type == PPMERGE_TYPE_SHUTTER) { + shutterRef = pmShutterCorrectionReference(shutters->data[cellNum]); + pattern = pmReadoutAlloc(NULL); + } + pmReadout *outRO = pmReadoutAlloc(outCell); - psArray *readouts = psArrayAlloc(numFiles); ///< Input readouts + // open the input files (we need to do the work ourselves) 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; } - + } + + ppMergeFileGroup *fileGroup = NULL; + 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", j); + pmCell *inCell = pmFPAviewThisCell(view, input->fpa); ///< Input cell + readouts->data[j] = pmReadoutAlloc(inCell); + } + + fileGroup = ppMergeFileGroupAlloc(); + fileGroup->readouts = readouts; + fileGroup->read = false; + fileGroup->busy = false; + fileGroup->lastScan = 0; + fileGroup->firstScan = 0; + fileGroups->data[i] = fileGroup; + } + + // call the init functions + switch (type) { + case PPMERGE_TYPE_BIAS: + case PPMERGE_TYPE_FLAT: + case PPMERGE_TYPE_FRINGE: + psAssert (fileGroups->n > 0, "no valid file groups defined"); + fileGroup = fileGroups->data[0]; + if (!pmReadoutCombinePrepare(outRO, fileGroup->readouts, combination)) { + goto ERROR; + } + break; + case PPMERGE_TYPE_DARK: + psAssert (fileGroups->n > 0, "no valid file groups defined"); + fileGroup = fileGroups->data[0]; + if (!pmDarkCombinePrepare(outCell, fileGroup->readouts, darkOrdinates, darkNorm)) { + goto ERROR; + } + break; + case PPMERGE_TYPE_SHUTTER: + psAssert (fileGroups->n > 0, "no valid file groups defined"); + fileGroup = fileGroups->data[0]; + if (!pmShutterCorrectionGeneratePrepare(outRO, pattern, fileGroup->readouts, maskVal)) { + goto ERROR; + } + break; + default: + psAbort("Should never get here."); + } + + // Read input data by chunks + psTimerStart("ppMergeLoop"); + for (int numChunk = 0; true; numChunk++) { + + bool status = false; + fileGroup = ppMergeReadChunk(&status, fileGroups, config, numChunk); + if (!status) { + // Something went wrong + goto ERROR; + } + if (!fileGroup) { + // Nothing more to read + break; + } + + // Start a job + switch (type) { + case PPMERGE_TYPE_BIAS: + case PPMERGE_TYPE_FLAT: + case PPMERGE_TYPE_FRINGE: { + psThreadJob *job = psThreadJobAlloc("PPMERGE_READOUT_COMBINE"); ///< Job to start + + // 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); + + // call: pmReadoutCombine(outRO, fileGroup->readouts, zeros, scales, combination); + if (!psThreadJobAddPending(job)) { + goto ERROR; + } + break; + } + case PPMERGE_TYPE_DARK: { + psThreadJob *job = psThreadJobAlloc ("PPMERGE_DARK_COMBINE"); ///< Job to start + + // construct the arguments for this job + psArrayAdd(job->args, 1, outCell); + psArrayAdd(job->args, 1, fileGroup); + 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_IMAGE_MASK)); + + // call: pmDarkCombine(outCell, fileGroup->readouts, iter, rej, maskVal); + if (!psThreadJobAddPending(job)) { + goto ERROR; + } + break; + } + case PPMERGE_TYPE_SHUTTER: { + psThreadJob *job = psThreadJobAlloc ("PPMERGE_SHUTTER_CORRECTION"); + + // construct the arguments for this job + psArrayAdd(job->args, 1, outRO); + psArrayAdd(job->args, 1, pattern); + 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_IMAGE_MASK)); + + // call: pmShutterCorrectionGenerate(outRO, pattern, fileGroup->readouts, shutterRef, + // shutters->data[cellNum], iter, rej, maskVal); + if (!psThreadJobAddPending (job)) { + goto ERROR; + } + break; + } + default: + psAbort("Should never get here."); + } + } + + // Wait for the threads to finish and manage results + if (!psThreadPoolWait(false)) { + psError(PS_ERR_UNKNOWN, false, "Unable to combine images."); + return false; + } + + // we don't care about the results, just dump the done queue jobs + psThreadJob *job = NULL; // Job to dump + while ((job = psThreadJobGetDone())) { + psFree (job); + } + + psFree(fileGroups); + + // XXX eventually need to keep both the shutter and the pattern, as we do with dark + psFree(pattern); + + // Get list of cells for concepts averaging + psList *inCells = psListAlloc(NULL); // List of cells + for (int i = 0; i < numFiles; i++) { 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); + psListAdd(inCells, PS_LIST_TAIL, inCell); } if (!pmConceptsAverageCells(outCell, inCells, NULL, NULL, true)) { @@ -253,6 +342,4 @@ } psFree(inCells); - - psFree(readouts); // Plug supplementary images into their own FPAs @@ -266,11 +353,11 @@ psFree(countsRO); - // XXX EAM 2009.01.18 : sigmaCell and countsCell need to have their concepts copied from outCell. - // This was causing segfaults for VYSOS5; Why did this ever work for SIMTEST? - if (!pmConceptsCopyCell(countsCell, outCell)) { - psError(PS_ERR_UNKNOWN, false, "Unable to copy cell concepts."); - psFree(outRO); - goto ERROR; - } + // XXX EAM 2009.01.18 : sigmaCell and countsCell need to have their concepts copied from outCell. + // This was causing segfaults for VYSOS5; Why did this ever work for SIMTEST? + if (!pmConceptsCopyCell(countsCell, outCell)) { + psError(PS_ERR_UNKNOWN, false, "Unable to copy cell concepts."); + psFree(outRO); + goto ERROR; + } pmCell *sigmaCell = pmFPAfileThisCell(config->files, view, "PPMERGE.OUTPUT.SIGMA"); @@ -282,9 +369,9 @@ psFree(sigmaRO); - if (!pmConceptsCopyCell(sigmaCell, outCell)) { - psError(PS_ERR_UNKNOWN, false, "Unable to copy cell concepts."); - psFree(outRO); - goto ERROR; - } + if (!pmConceptsCopyCell(sigmaCell, outCell)) { + psError(PS_ERR_UNKNOWN, false, "Unable to copy cell concepts."); + psFree(outRO); + goto ERROR; + } } @@ -316,5 +403,5 @@ } - if (!ppStatsFPA(stats, outFPA, view, maskVal, config)) { + if (stats && !ppStatsFPA(stats, outFPA, view, maskVal, config)) { psError(PS_ERR_UNKNOWN, true, "Unable to generate stats for image."); goto ERROR; @@ -369,5 +456,5 @@ psFree(inputs); psFree(masks); - psFree(weights); + psFree(variances); psFree(stats); return true; @@ -378,5 +465,5 @@ psFree(inputs); psFree(masks); - psFree(weights); + psFree(variances); psFree(stats); return false; Index: unk/ppMerge/src/ppMergeLoop_Threaded.c =================================================================== --- /trunk/ppMerge/src/ppMergeLoop_Threaded.c (revision 23409) +++ (revision ) @@ -1,471 +1,0 @@ -/** @file ppMergeLoop_Threaded.c - * - * @brief - * - * @ingroup ppMerge - * - * @author IfA - * @version $Revision: 1.9 $ $Name: not supported by cvs2svn $ - * @date $Date: 2009-02-06 02:44:31 $ - * Copyright 2009 Institute for Astronomy, University of Hawaii - */ - -#ifdef HAVE_CONFIG_H -#include -#endif - -#include -#include - -#include -#include -#include - -#include "ppMerge.h" - -// XXX this function is now sufficiently different for the major types, it would make sense to just -// split it into three: BASIC, SHUTTER, DARK - -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 haveVariances = psMetadataLookupBool(&mdok, arguments, "INPUTS.VARIANCES"); // Do we have variances? - - psArray *inputs = ppMergeFileDataLevel(config, "PPMERGE.INPUT", PM_FPA_LEVEL_READOUT); // Input images - psArray *masks = NULL, *variances = NULL; // Input masks and variances - if (haveMasks) { - masks = ppMergeFileDataLevel(config, "PPMERGE.INPUT.MASK", PM_FPA_LEVEL_READOUT); - } - if (haveVariances) { - variances = ppMergeFileDataLevel(config, "PPMERGE.INPUT.VARIANCE", PM_FPA_LEVEL_READOUT); - } - - int nThreads = psMetadataLookupS32 (&mdok, arguments, "NTHREADS"); - if (!mdok) nThreads = 0; - - // General combination parameters - 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 useVariances = psMetadataLookupBool(NULL, arguments, "VARIANCES"); // Use variances? - - // 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 - psImageMaskType maskVal, 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->variances = useVariances; - - 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 - pmHDU *lastHDU = NULL; // Last HDU that was updated - 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; - } - - // Update the header - { - pmHDU *hdu = pmHDUGetHighest(outFPA, outChip, outCell); // HDU for file - if (hdu && hdu != lastHDU) { - ppMergeVersionHeader(hdu->header); - lastHDU = hdu; - } - } - - float shutterRef = NAN; ///< Reference shutter correction - pmReadout *pattern = NULL; - if (type == PPMERGE_TYPE_SHUTTER) { - shutterRef = pmShutterCorrectionReference(shutters->data[cellNum]); - pattern = pmReadoutAlloc(NULL); - } - - 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; - } - } - - ppMergeFileGroup *fileGroup = NULL; - 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", j); - pmCell *inCell = pmFPAviewThisCell(view, input->fpa); ///< Input cell - readouts->data[j] = pmReadoutAlloc(inCell); - } - - fileGroup = ppMergeFileGroupAlloc(); - fileGroup->readouts = readouts; - fileGroup->read = false; - fileGroup->busy = false; - fileGroup->lastScan = 0; - fileGroup->firstScan = 0; - fileGroups->data[i] = fileGroup; - } - - // call the init functions - switch (type) { - case PPMERGE_TYPE_BIAS: - case PPMERGE_TYPE_FLAT: - case PPMERGE_TYPE_FRINGE: - psAssert (fileGroups->n > 0, "no valid file groups defined"); - fileGroup = fileGroups->data[0]; - if (!pmReadoutCombinePrepare(outRO, fileGroup->readouts, combination)) { - goto ERROR; - } - break; - case PPMERGE_TYPE_DARK: - psAssert (fileGroups->n > 0, "no valid file groups defined"); - fileGroup = fileGroups->data[0]; - if (!pmDarkCombinePrepare(outCell, fileGroup->readouts, darkOrdinates, darkNorm)) { - goto ERROR; - } - break; - case PPMERGE_TYPE_SHUTTER: - psAssert (fileGroups->n > 0, "no valid file groups defined"); - fileGroup = fileGroups->data[0]; - if (!pmShutterCorrectionGeneratePrepare(outRO, pattern, fileGroup->readouts, maskVal)) { - goto ERROR; - } - break; - default: - psAbort("Should never get here."); - } - - // Read input data by chunks - psTimerStart("ppMergeLoop"); - for (int numChunk = 0; true; numChunk++) { - - bool status = false; - fileGroup = ppMergeReadChunk(&status, fileGroups, config, numChunk); - if (!status) { - // Something went wrong - goto ERROR; - } - if (!fileGroup) { - // Nothing more to read - break; - } - - // Start a job - switch (type) { - case PPMERGE_TYPE_BIAS: - case PPMERGE_TYPE_FLAT: - case PPMERGE_TYPE_FRINGE: { - psThreadJob *job = psThreadJobAlloc("PPMERGE_READOUT_COMBINE"); ///< Job to start - - // 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); - - // call: pmReadoutCombine(outRO, fileGroup->readouts, zeros, scales, combination); - if (!psThreadJobAddPending(job)) { - goto ERROR; - } - break; - } - case PPMERGE_TYPE_DARK: { - psThreadJob *job = psThreadJobAlloc ("PPMERGE_DARK_COMBINE"); ///< Job to start - - // construct the arguments for this job - psArrayAdd(job->args, 1, outCell); - psArrayAdd(job->args, 1, fileGroup); - 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_IMAGE_MASK)); - - // call: pmDarkCombine(outCell, fileGroup->readouts, iter, rej, maskVal); - if (!psThreadJobAddPending(job)) { - goto ERROR; - } - break; - } - case PPMERGE_TYPE_SHUTTER: { - psThreadJob *job = psThreadJobAlloc ("PPMERGE_SHUTTER_CORRECTION"); - - // construct the arguments for this job - psArrayAdd(job->args, 1, outRO); - psArrayAdd(job->args, 1, pattern); - 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_IMAGE_MASK)); - - // call: pmShutterCorrectionGenerate(outRO, pattern, fileGroup->readouts, shutterRef, - // shutters->data[cellNum], iter, rej, maskVal); - if (!psThreadJobAddPending (job)) { - goto ERROR; - } - break; - } - default: - psAbort("Should never get here."); - } - } - - // Wait for the threads to finish and manage results - if (!psThreadPoolWait(false)) { - psError(PS_ERR_UNKNOWN, false, "Unable to combine images."); - return false; - } - - // we don't care about the results, just dump the done queue jobs - psThreadJob *job = NULL; // Job to dump - while ((job = psThreadJobGetDone())) { - psFree (job); - } - - psFree(fileGroups); - - // XXX eventually need to keep both the shutter and the pattern, as we do with dark - psFree(pattern); - - // Get list of cells for concepts averaging - psList *inCells = psListAlloc(NULL); // List of cells - for (int i = 0; i < numFiles; i++) { - pmFPAfile *input = pmFPAfileSelectSingle(config->files, "PPMERGE.INPUT", i); - pmCell *inCell = pmFPAviewThisCell(view, input->fpa); ///< Input cell - psListAdd(inCells, PS_LIST_TAIL, inCell); - } - 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); - - // 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); - - // XXX EAM 2009.01.18 : sigmaCell and countsCell need to have their concepts copied from outCell. - // This was causing segfaults for VYSOS5; Why did this ever work for SIMTEST? - if (!pmConceptsCopyCell(countsCell, outCell)) { - psError(PS_ERR_UNKNOWN, false, "Unable to copy cell concepts."); - psFree(outRO); - goto ERROR; - } - - 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); - - if (!pmConceptsCopyCell(sigmaCell, outCell)) { - psError(PS_ERR_UNKNOWN, false, "Unable to copy cell concepts."); - psFree(outRO); - goto ERROR; - } - } - - // 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 (stats && !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(variances); - psFree(stats); - return true; - -ERROR: - psFree(view); - psFree(combination); - psFree(inputs); - psFree(masks); - psFree(variances); - psFree(stats); - return false; -} -