Index: trunk/ppStack/src/Makefile.am =================================================================== --- trunk/ppStack/src/Makefile.am (revision 19283) +++ trunk/ppStack/src/Makefile.am (revision 19337) @@ -14,5 +14,6 @@ ppStackVersion.c \ ppStackMatch.c \ - ppStackSources.c + ppStackSources.c \ + ppStackThread.c noinst_HEADERS = \ Index: trunk/ppStack/src/ppStack.h =================================================================== --- trunk/ppStack/src/ppStack.h (revision 19283) +++ trunk/ppStack/src/ppStack.h (revision 19337) @@ -5,5 +5,57 @@ #define PPSTACK_INSPECT_PIXELS "PPSTACK.PIXELS" // Name of rejected pixels metadata items +#include #include + +// Mask values for inputs +typedef enum { + PPSTACK_MASK_MATCH = 0x01, // PSF-matching failed + PPSTACK_MASK_REJECT = 0x02, // Rejection failed + PPSTACK_MASK_BAD = 0x04, // Bad image (too many pixels rejected) + PPSTACK_MASK_ALL = 0xff // All errors +} ppStackMask; + +// Thread for stacking chunks +// +// Each input file contributes a readout, into which is read a chunk from that file +typedef struct { + psArray *readouts; // Input readouts to read and stack + bool read; // Has the scan been read? + bool busy; // Is the scan being processed? + int firstScan; // First row of the chunk to be read for this group + int lastScan; // Last row of the chunk to be read for this group +} ppStackThread; + +// Allocator +ppStackThread *ppStackThreadAlloc(psArray *readouts // Inputs readouts to read and stack + ); + +// Data for threads +typedef struct { + psArray *threads; // Threads doing stacking + int lastScan; // Last row that's been read + psArray *imageFits; // FITS file pointers for images + psArray *maskFits; // FITS file pointers for masks + psArray *weightFits; // FITS file pointers for weights +} ppStackThreadData; + +// Set up thread data +ppStackThreadData *ppStackThreadDataSetup(const psArray *cells, // Array of input cells + const psArray *imageNames, // Names of images to read + const psArray *maskNames, // Names of masks to read + const psArray *weightNames, // Names of weight maps to read + const pmConfig *config // Configuration + ); + +// Read chunk into the first available file thread +ppStackThread *ppStackThreadRead(bool *status, // Status of read + ppStackThreadData *stack, // Stacks available for reading + pmConfig *config, // Configuration + int numChunk, // Chunk number (only for interest) + int overlap // Overlap between subsequent scans + ); + +// Initialise the threads +void ppStackThreadInit(void); // Setup command-line arguments @@ -31,9 +83,15 @@ // Perform stacking on a readout -bool ppStackReadoutInitial(const pmConfig *config, // Configuration - pmReadout *outRO, // Output readout - const psArray *readouts, // Input readouts - const psArray *regions, // Array with array of regions used in each PSF matching - const psArray *kernels // Array with array of kernels used in each PSF matching +// +// Returns an array of pixels to inspect for each input image +psArray *ppStackReadoutInitial(const pmConfig *config, // Configuration + pmReadout *outRO, // Output readout + const psArray *readouts, // Input readouts + const psArray *regions, // Array with array of regions used in each PSF match + const psArray *kernels // Array with array of kernels used in each PSF match + ); + +// Thread entry point for ppStackReadoutInitial +bool ppStackReadoutInitialThread(psThreadJob *job // Job to process ); @@ -43,4 +101,8 @@ const psArray *readouts, // Input readouts const psArray *rejected // Array with pixels rejected in each image + ); + +// Thread entry point for ppStackReadoutFinal +bool ppStackReadoutFinalThread(psThreadJob *job // Job to process ); @@ -62,11 +124,11 @@ /// Convolve image to match specified seeing -bool ppStackMatch(pmReadout *readout, ///< Readout to be convolved; replaced with output +bool ppStackMatch(pmReadout *readout, // Readout to be convolved; replaced with output psArray **regions, // Array of regions used in each PSF matching, returned psArray **kernels, // Array of kernels used in each PSF matching, returned - const psArray *sources, ///< Array of sources - const pmPSF *psf, ///< Target PSF - psRandom *rng, ///< Random number generator - const pmConfig *config ///< Configuration + const psArray *sources, // Array of sources + const pmPSF *psf, // Target PSF + psRandom *rng, // Random number generator + const pmConfig *config // Configuration ); @@ -75,7 +137,7 @@ /// /// Corrects the sources to have consistent magnitudes. Returns the source lists. -psArray *ppStackSourceListAdd(psArray *lists, ///< List to which to add, or NULL - psArray *sources, ///< Sources to add - const pmConfig *config ///< Configuration +psArray *ppStackSourceListAdd(psArray *lists, // List to which to add, or NULL + psArray *sources, // Sources to add + const pmConfig *config // Configuration ); @@ -83,10 +145,10 @@ /// /// Corrects the sources to have consistent magnitudes where possible. Returns the sources -psArray *ppStackSourceListCombine(psArray *lists, ///< Source lists - const pmConfig *config ///< Configuration +psArray *ppStackSourceListCombine(psArray *lists, // Source lists + const pmConfig *config // Configuration ); /// Print sources into a ds9 regions file -void ppStackSourcesPrint(const psArray *sources ///< Sources to print +void ppStackSourcesPrint(const psArray *sources // Sources to print ); Index: trunk/ppStack/src/ppStackLoop.c =================================================================== --- trunk/ppStack/src/ppStackLoop.c (revision 19283) +++ trunk/ppStack/src/ppStackLoop.c (revision 19337) @@ -12,5 +12,5 @@ #include "ppStack.h" -#define TESTING +//#define TESTING #define WCS_TOLERANCE 0.001 // Tolerance for WCS @@ -95,27 +95,4 @@ } -#if 0 -// Set the data level for a list of files -static void fileSetDataLevel(pmConfig *config, // Configuration - char **files, // Files for which to set level - pmFPALevel level // Level to set - ) -{ - assert(config); - assert(files); - - for (int i = 0; files[i] != NULL; i++) { - psArray *selected = pmFPAfileSelect(config->files, files[i]); // Selected files of interest - for (int j = 0; j < selected->n; j++) { - pmFPAfile *file = selected->data[j]; - assert(file); - file->dataLevel = level; - } - psFree(selected); - } - return; -} -#endif - // Iterate down the hierarchy, loading files; we can get away with this because we're working on skycells static pmFPAview *filesIterateDown(pmConfig *config // Configuration @@ -244,5 +221,4 @@ } int num = psMetadataLookupS32(NULL, config->arguments, "INPUTS.NUM"); // Number of inputs - int numScans = psMetadataLookupS32(NULL, recipe, "ROWS"); // Number of scans to read at once psMetadata *ppsub = psMetadataLookupMetadata(NULL, config->recipes, "PPSUB"); // PPSUB recipe @@ -412,4 +388,6 @@ int numGood = 0; // Number of good frames int numCols = 0, numRows = 0; // Size of image + psVector *inputMask = psVectorAlloc(num, PS_TYPE_U8); // Mask for inputs + psVectorInit(inputMask, 0); for (int i = 0; i < num; i++) { psTrace("ppStack", 2, "Convolving input %d of %d to target PSF....\n", i, num); @@ -423,4 +401,5 @@ psFree(targetPSF); psFree(rng); + psFree(inputMask); return false; } @@ -439,4 +418,5 @@ psFree(targetPSF); psFree(rng); + psFree(inputMask); return false; } @@ -447,4 +427,5 @@ if (!ppStackMatch(readout, ®ions, &kernels, sources, targetPSF, rng, config)) { psErrorStackPrint(stderr, "Unable to match image %d --- ignoring.", i); + inputMask->data.U8[i] = PPSTACK_MASK_MATCH; psErrorClear(); continue; @@ -474,7 +455,10 @@ if (numGood == 0) { psError(PS_ERR_UNKNOWN, false, "No good images to combine."); - return false; - } - + psFree(subKernels); + psFree(subRegions); + psFree(cells); + psFree(inputMask); + return false; + } #ifdef TESTING @@ -482,155 +466,140 @@ #endif + ppStackThreadInit(); + ppStackThreadData *stack = ppStackThreadDataSetup(cells, imageNames, maskNames, weightNames, + config); // Data for stacking + psFree(cells); + memDump("preinitial"); - // Stack the convolved files + // Stack the convolvd files psTrace("ppStack", 1, "Initial stack of convolved images....\n"); + pmReadout *outRO = NULL; // Output readout + pmFPAview *view = NULL; // View to readout { + int row0, col0; // Offset for readout + int numCols, numRows; // Size of readout + ppStackThread *thread = stack->threads->data[0]; // Representative thread + if (!pmReadoutStackSetOutputSize(&col0, &row0, &numCols, &numRows, thread->readouts)) { + psError(PS_ERR_UNKNOWN, false, "problem setting output readout size."); + psFree(subKernels); + psFree(subRegions); + psFree(stack); + psFree(inputMask); + return false; + } + pmFPAfileActivate(config->files, false, NULL); fileActivation(config, combineFiles, true); - pmFPAview *view = filesIterateDown(config); + view = filesIterateDown(config); if (!view) { - psFree(cells); psFree(subKernels); psFree(subRegions); - return false; - } + psFree(stack); + psFree(inputMask); + return false; + } + pmCell *outCell = pmFPAfileThisCell(config->files, view, "PPSTACK.OUTPUT"); // Output cell - pmReadout *outRO = pmReadoutAlloc(outCell); // Output readout - - psArray *readouts = psArrayAlloc(num); // Readouts for convolved images - for (int i = 0; i < num; i++) { - pmCell *cell = cells->data[i]; // Cell of interest - if (!cell) { - // Bad image - continue; - } - pmReadout *ro = cell->readouts->data[0]; // Readout of interest - if (!ro) { - ro = pmReadoutAlloc(cell); - } - readouts->data[i] = ro; // Readout into which to read - } - psFree(cells); - - // FITS files - psArray *imageFits = psArrayAlloc(num); - psArray *maskFits = psArrayAlloc(num); - psArray *weightFits = psArrayAlloc(num); - for (int i = 0; i < num; i++) { - if (!readouts->data[i]) { - // Bad image - continue; - } - // Resolved names - psString imageResolved = pmConfigConvertFilename(imageNames->data[i], config, false, false); - psString maskResolved = pmConfigConvertFilename(maskNames->data[i], config, false, false); - psString weightResolved = pmConfigConvertFilename(weightNames->data[i], config, false, false); - imageFits->data[i] = psFitsOpen(imageResolved, "r"); - maskFits->data[i] = psFitsOpen(maskResolved, "r"); - weightFits->data[i] = psFitsOpen(weightResolved, "r"); - psFree(imageResolved); - psFree(maskResolved); - psFree(weightResolved); - if (!imageFits->data[i] || !maskFits->data[i] || !weightFits->data[i]) { - psError(PS_ERR_UNKNOWN, false, "Unable to open convolved files %s, %s, %s", - (char*)imageNames->data[i], (char*)maskNames->data[i], (char*)weightNames->data[i]); + outRO = pmReadoutAlloc(outCell); // Output readout + + psString maskBadStr = psMetadataLookupStr(NULL, recipe, "MASK.BAD"); // Name of bits to mask for bad + psMaskType maskBad = pmConfigMaskGet(maskBadStr, config); // Bits to mask for bad pixels + if (!pmReadoutStackDefineOutput(outRO, col0, row0, numCols, numRows, true, true, maskBad)) { + psError(PS_ERR_UNKNOWN, false, "Unable to prepare output."); + psFree(subKernels); + psFree(subRegions); + psFree(stack); + psFree(inputMask); + psFree(view); + psFree(outRO); + return false; + } + + bool status; // Status of read + for (int numChunk = 0; true; numChunk++) { + ppStackThread *thread = ppStackThreadRead(&status, stack, config, numChunk, overlap); + if (!status) { + // Something went wrong + psError(PS_ERR_UNKNOWN, false, "Unable to read chunk %d", numChunk); psFree(subKernels); psFree(subRegions); + psFree(stack); + psFree(inputMask); + psFree(view); + psFree(outRO); return false; } - } - - // 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); + if (!thread) { + // Nothing more to read + break; + } + + // call: ppStackReadoutInitial(config, outRO, readouts, subRegions, subKernels) + psThreadJob *job = psThreadJobAlloc("PPSTACK_INITIAL_COMBINE"); // Job to start + psArrayAdd(job->args, 1, thread); + psArrayAdd(job->args, 1, config); + psArrayAdd(job->args, 1, outRO); + psArrayAdd(job->args, 1, subRegions); + psArrayAdd(job->args, 1, subKernels); + if (!psThreadJobAddPending(job)) { + psFree(job); + psFree(subKernels); + psFree(subRegions); + psFree(stack); + psFree(inputMask); + psFree(view); + psFree(outRO); + return false; + } + psFree(job); + } + + if (!psThreadPoolWait(false)) { + psError(PS_ERR_UNKNOWN, false, "Unable to do initial combination."); + psFree(subKernels); + psFree(subRegions); + psFree(stack); + psFree(inputMask); + psFree(view); + psFree(outRO); + return false; + } + + memDump("initial"); + } + + // Pixel rejection + psArray *rejected = psArrayAlloc(num); + { + psArray *inspect = psArrayAlloc(num); // Pixels to inspect + int numRejected = 0; // Number of inputs rejected completely + + // Count images rejected out of hand + for (int i = 0; i < num; i++) { + if (inputMask->data.U8[i]) { + numRejected++; + } + inspect->data[i] = psPixelsAllocEmpty(PIXELS_BUFFER); + } + + // Harvest the jobs, combining the inspection lists + psThreadJob *job; // Completed job + while ((job = psThreadJobGetDone())) { + psArray *results = job->results; // Results of job + psAssert(results && results->n == num, "Results array."); for (int i = 0; i < num; i++) { - pmReadout *readout = readouts->data[i]; - if (!readout) { - // Bad image + if (inputMask->data.U8[i]) { continue; } - - if (!pmReadoutReadChunk(readout, imageFits->data[i], 0, numScans, overlap, - config) || - !pmReadoutReadChunkMask(readout, maskFits->data[i], 0, numScans, overlap, - config) || - !pmReadoutReadChunkWeight(readout, weightFits->data[i], 0, numScans, overlap, - config)) { - psError(PS_ERR_IO, false, "Unable to read chunk %d for file %d", numChunk, i); - psFree(readouts); - psFree(subKernels); - psFree(subRegions); - psFree(outRO); - psFree(view); - return false; - } - } - -#ifndef PS_NO_TRACE - { - pmReadout *ro = NULL; // Representative readout - for (int i = 0; i < num && !ro; i++) { - ro = readouts->data[i]; - } - psAssert(ro, "There should be a readout here."); - psTrace("ppStack", 4, "Stack: [%d:%d,%d:%d]\n", ro->col0, ro->col0 + ro->image->numCols, - ro->row0, ro->row0 + ro->image->numRows); - } -#endif - - if (!ppStackReadoutInitial(config, outRO, readouts, subRegions, subKernels)) { - psError(PS_ERR_UNKNOWN, false, "Unable to stack images.\n"); - psFree(readouts); - psFree(subKernels); - psFree(subRegions); - psFree(outRO); - psFree(view); - return false; - } - - for (int i = 0; i < num && more; i++) { - pmReadout *readout = readouts->data[i]; - if (!readout) { - // Bad image - continue; - } - more &= pmReadoutMore(readout, imageFits->data[i], 0, numScans, config); - more &= pmReadoutMoreMask(readout, maskFits->data[i], 0, numScans, config); - more &= pmReadoutMoreWeight(readout, weightFits->data[i], 0, numScans, config); - } - } - - memDump("initial"); - - // Reset for the second read - // Extract the rejection lists - psArray *rejected = psArrayAlloc(num); // Pixels to inspect - int numRejected = 0; // Number of inputs rejected completely + inspect->data[i] = psPixelsConcatenate(inspect->data[i], results->data[i]); + } + psFree(job); + } + for (int i = 0; i < num; i++) { - pmReadout *ro = readouts->data[i]; // Readout of interest - if (!ro) { - // Bad image - numRejected++; + if (inputMask->data.U8[i]) { continue; } - psPixels *inspect = psPixelsAllocEmpty(PIXELS_BUFFER); // Inspection list for this readout - psMetadataIterator *iter = psMetadataIteratorAlloc(ro->analysis, PS_LIST_HEAD, - "^" PPSTACK_INSPECT_PIXELS "$"); // Iterator - psMetadataItem *item; - while ((item = psMetadataGetAndIncrement(iter))) { - psPixels *pixels = item->data.V; // Rejected pixels - if (!pixels) { - continue; - } - psTrace("ppStack", 5, "Adding %ld pixels to inspect to image %d", pixels->n, i); - inspect = psPixelsConcatenate(inspect, pixels); - } - psFree(iter); - psMetadataRemoveKey(ro->analysis, PPSTACK_INSPECT_PIXELS); - pmReadoutFreeData(ro); - - psLogMsg("ppStack", PS_LOG_INFO, "%ld total pixels to inspect from image %d", inspect->n, i); #ifdef TESTING @@ -648,6 +617,6 @@ #endif - psPixels *reject = pmStackReject(inspect, numCols, numRows, threshold, poorFrac, - subRegions->data[i], subKernels->data[i]); // Pixels to reject + psPixels *reject = pmStackReject(inspect->data[i], numCols, numRows, threshold, poorFrac, + subRegions->data[i], subKernels->data[i]); // Rejected pixels #ifdef TESTING @@ -665,12 +634,15 @@ #endif - psFree(inspect); + psFree(inspect->data[i]); + inspect->data[i] = NULL; + if (!reject) { psWarning("Rejection on image %d didn't work --- reject entire image.", i); numRejected++; + inputMask->data.U8[i] = PPSTACK_MASK_REJECT; } else { - float frac = reject->n / (float)(outRO->image->numCols * outRO->image->numRows); // Pixel frac + float frac = reject->n / (float)(numCols * numRows); // Pixel fraction psLogMsg("ppStack", PS_LOG_INFO, "%ld pixels rejected from image %d (%.1f%%)", - reject->n, i, frac * 100.0); + reject->n, i, frac * 100.0); if (frac > imageRej) { psWarning("Image %d rejected completely because rejection fraction (%.3f) " @@ -679,12 +651,12 @@ // reject == NULL means reject image completely reject = NULL; + inputMask->data.U8[i] = PPSTACK_MASK_BAD; numRejected++; } } rejected->data[i] = reject; - - memDump("reject"); - - } + } + + psFree(inspect); psFree(subKernels); psFree(subRegions); @@ -692,204 +664,184 @@ if (numRejected == num) { psError(PS_ERR_UNKNOWN, true, "All inputs completely rejected; unable to proceed."); - psFree(readouts); + psFree(stack); psFree(rejected); + psFree(inputMask); + psFree(view); psFree(outRO); + return false; + } + } + + memDump("reject"); + + // Final combination + psTrace("ppStack", 2, "Final stack of convolved images....\n"); + { + stack->lastScan = 0; // Reset read + bool status; // Status of read + for (int numChunk = 0; true; numChunk++) { + ppStackThread *thread = ppStackThreadRead(&status, stack, config, numChunk, 0); + if (!status) { + // Something went wrong + psError(PS_ERR_UNKNOWN, false, "Unable to read chunk %d", numChunk); + psFree(stack); + psFree(rejected); + psFree(inputMask); + psFree(view); + psFree(outRO); + return false; + } + if (!thread) { + // Nothing more to read + break; + } + + // call: ppStackReadoutFinal(config, outRO, readouts, rejected) + psThreadJob *job = psThreadJobAlloc("PPSTACK_FINAL_COMBINE"); // Job to start + psArrayAdd(job->args, 1, thread); + psArrayAdd(job->args, 1, config); + psArrayAdd(job->args, 1, outRO); + psArrayAdd(job->args, 1, rejected); + if (!psThreadJobAddPending(job)) { + psFree(job); + psFree(stack); + psFree(rejected); + psFree(inputMask); + psFree(view); + psFree(outRO); + return false; + } + psFree(job); + } + + if (!psThreadPoolWait(true)) { + psError(PS_ERR_UNKNOWN, false, "Unable to do final combination."); + psFree(stack); + psFree(inputMask); + psFree(rejected); psFree(view); - return false; - } - - memDump("prefinal"); - - // Read convolutions by chunks - psTrace("ppStack", 2, "Final stack of convolved images....\n"); - more = true; - for (int numChunk = 0; more; numChunk++) { - psTrace("ppStack", 2, "Final stack of chunk %d....\n", numChunk); - for (int i = 0; i < num; i++) { - if (!rejected->data[i]) { - continue; - } - pmReadout *readout = readouts->data[i]; - assert(readout); - - if (!pmReadoutReadChunk(readout, imageFits->data[i], 0, numScans, 0, config) || - !pmReadoutReadChunkMask(readout, maskFits->data[i], 0, numScans, 0, config) || - !pmReadoutReadChunkWeight(readout, weightFits->data[i], 0, numScans, 0, - config)) { - psError(PS_ERR_IO, false, "Unable to read chunk %d for file %d", numChunk, i); - psFree(readouts); - psFree(rejected); - psFree(outRO); - psFree(view); - return false; - } - } - -#ifndef PS_NO_TRACE - { - pmReadout *ro = NULL; - for (int i = 0; i < num && !ro; i++) { - if (!rejected->data[i]) { - continue; + psFree(outRO); + return false; + } + } + + memDump("final"); + + psTrace("ppStack", 2, "Cleaning up after combination....\n"); + + // Ensure masked regions really look masked + { + psString maskBadStr = psMetadataLookupStr(NULL, recipe, "MASK.BAD"); // Name of bits for bad + psMaskType maskBad = pmConfigMaskGet(maskBadStr, config); // Bits to mask for bad pixels + if (!pmReadoutMaskApply(outRO, maskBad)) { + psWarning("Unable to apply mask"); + } + } + + // Close up + bool wcsDone = false; // Have we done the WCS? + for (int i = 0; i < num; i++) { + if (!inputMask->data.U8[i]) { + continue; + } + if (!wcsDone) { + // Copy astrometry over + wcsDone = true; + ppStackThread *thread = stack->threads->data[0]; // Representative stack + pmReadout *inRO = thread->readouts->data[i]; // Template readout + pmHDU *inHDU = pmHDUFromCell(inRO->parent); // Template HDU + pmHDU *outHDU = pmHDUFromCell(outRO->parent); // Output HDU + pmChip *outChip = outRO->parent->parent; // Output chip + pmFPA *outFPA = outChip->parent; // Output FPA + if (!outHDU || !inHDU) { + psWarning("Unable to find HDU at FPA level to copy astrometry."); + } else { + if (!pmAstromReadWCS(outFPA, outChip, inHDU->header, 1.0)) { + psErrorClear(); + psWarning("Unable to read WCS astrometry from input FPA."); + wcsDone = false; + } else { + if (!outHDU->header) { + outHDU->header = psMetadataAlloc(); } - ro = readouts->data[i]; - } - if (ro && ro->image) { - psTrace("ppStack", 4, "Stack: [%d:%d,%d:%d]\n", ro->col0, ro->col0 + ro->image->numCols, - ro->row0, ro->row0 + ro->image->numRows); - } - } -#endif - - if (!ppStackReadoutFinal(config, outRO, readouts, rejected)) { - psError(PS_ERR_UNKNOWN, false, "Unable to stack images.\n"); - psFree(readouts); - psFree(rejected); - psFree(outRO); - psFree(view); - return false; - } - - for (int i = 0; i < num && more; i++) { - pmReadout *readout = readouts->data[i]; - if (!readout) { - continue; - } - more &= pmReadoutMore(readout, imageFits->data[i], 0, numScans, config); - more &= pmReadoutMoreMask(readout, maskFits->data[i], 0, numScans, config); - more &= pmReadoutMoreWeight(readout, weightFits->data[i], 0, numScans, config); - } - } - - psFree(rejected); - for (int i = 0; i < readouts->n; i++) { - pmReadout *ro = readouts->data[i]; // Readout of interest - if (!ro) { - continue; - } - pmReadoutFreeData(ro); - } - psFree(readouts); - - // Ensure masked regions really look masked - { - psString maskBadStr = psMetadataLookupStr(NULL, recipe, "MASK.BAD"); // Name of bits for bad - psMaskType maskBad = pmConfigMaskGet(maskBadStr, config); // Bits to mask for bad pixels - if (!pmReadoutMaskApply(outRO, maskBad)) { - psWarning("Unable to apply mask"); - } - } - - // Close up - bool wcsDone = false; // Have we done the WCS? - for (int i = 0; i < num; i++) { - if (!readouts->data[i]) { - continue; - } - if (!wcsDone) { - // Copy astrometry over - wcsDone = true; - pmReadout *inRO = readouts->data[i]; // Template readout - pmHDU *inHDU = pmHDUFromCell(inRO->parent); // Template HDU - pmHDU *outHDU = pmHDUFromCell(outRO->parent); // Output HDU - pmChip *outChip = outRO->parent->parent; // Output chip - pmFPA *outFPA = outChip->parent; // Output FPA - if (!outHDU || !inHDU) { - psWarning("Unable to find HDU at FPA level to copy astrometry."); - } else { - if (!pmAstromReadWCS(outFPA, outChip, inHDU->header, 1.0)) { + if (!pmAstromWriteWCS(outHDU->header, outFPA, outChip, WCS_TOLERANCE)) { psErrorClear(); - psWarning("Unable to read WCS astrometry from input FPA."); + psWarning("Unable to write WCS astrometry to output FPA."); wcsDone = false; - } else { - if (!outHDU->header) { - outHDU->header = psMetadataAlloc(); - } - if (!pmAstromWriteWCS(outHDU->header, outFPA, outChip, WCS_TOLERANCE)) { - psErrorClear(); - psWarning("Unable to write WCS astrometry to output FPA."); - wcsDone = false; - } } } } - - psFitsClose(imageFits->data[i]); - psFitsClose(maskFits->data[i]); - psFitsClose(weightFits->data[i]); - imageFits->data[i] = NULL; - maskFits->data[i] = NULL; - weightFits->data[i] = NULL; - if (tempDelete) { - psString imageResolved = pmConfigConvertFilename(imageNames->data[i], config, false, false); - psString maskResolved = pmConfigConvertFilename(maskNames->data[i], config, false, false); - psString weightResolved = pmConfigConvertFilename(weightNames->data[i], config, false, false); - if (unlink(imageResolved) == -1 || unlink(maskResolved) == -1 || - unlink(weightResolved) == -1) { - psWarning("Unable to delete temporary files for image %d", i); - } - psFree(imageResolved); - psFree(maskResolved); - psFree(weightResolved); - } - } - psFree(imageNames); - psFree(maskNames); - psFree(weightNames); - psFree(imageFits); - psFree(maskFits); - psFree(weightFits); - - memDump("final"); - - if (psMetadataLookupBool(&mdok, recipe, "PHOTOMETRY")) { - - fileActivation(config, combineFiles, false); - fileActivation(config, photFiles, true); - pmFPAview *photView = filesIterateDown(config); - - psTrace("ppStack", 1, "Photometering stacked image....\n"); - if (!ppStackPhotometry(config, outRO, view)) { - psError(PS_ERR_UNKNOWN, false, "Unable to perform photometry on output."); - psFree(outRO); - psFree(photView); - psFree(view); - return false; - } + } + + if (tempDelete) { + psString imageResolved = pmConfigConvertFilename(imageNames->data[i], config, false, false); + psString maskResolved = pmConfigConvertFilename(maskNames->data[i], config, false, false); + psString weightResolved = pmConfigConvertFilename(weightNames->data[i], config, false, false); + if (unlink(imageResolved) == -1 || unlink(maskResolved) == -1 || + unlink(weightResolved) == -1) { + psWarning("Unable to delete temporary files for image %d", i); + } + psFree(imageResolved); + psFree(maskResolved); + psFree(weightResolved); + } + } + psFree(imageNames); + psFree(maskNames); + psFree(weightNames); + + psFree(inputMask); + psFree(stack); + + memDump("cleanup"); + + if (psMetadataLookupBool(&mdok, recipe, "PHOTOMETRY")) { + psTrace("ppStack", 1, "Photometering stacked image....\n"); + + fileActivation(config, combineFiles, false); + fileActivation(config, photFiles, true); + pmFPAview *photView = filesIterateDown(config); + if (!ppStackPhotometry(config, outRO, photView)) { + psError(PS_ERR_UNKNOWN, false, "Unable to perform photometry on output."); + psFree(outRO); psFree(photView); - - fileActivation(config, combineFiles, true); - } - - memDump("phot"); - - // Statistics on output - if (stats) { - psTrace("ppStack", 1, "Gathering statistics on stacked image....\n"); - psString maskBadStr = psMetadataLookupStr(NULL, recipe, "MASK.BAD"); // Name of bits for bad - psMaskType maskBad = pmConfigMaskGet(maskBadStr, config); // Bits to mask for bad pixels - - ppStatsFPA(stats, outCell->parent->parent, view, maskBad, config); - } - - memDump("stats"); - - // Put version information into the header - pmHDU *hdu = pmHDUFromCell(outRO->parent); - if (!hdu) { - psError(PS_ERR_UNKNOWN, false, "Unable to find HDU for output."); - return false; - } - if (!hdu->header) { - hdu->header = psMetadataAlloc(); - } - ppStackVersionMetadata(hdu->header); - - // Write out the output files - filesIterateUp(config); - + return false; + } + psFree(photView); + + fileActivation(config, combineFiles, true); + } + + memDump("phot"); + + // Statistics on output + if (stats) { + psTrace("ppStack", 1, "Gathering statistics on stacked image....\n"); + psString maskBadStr = psMetadataLookupStr(NULL, recipe, "MASK.BAD"); // Name of bits for bad + psMaskType maskBad = pmConfigMaskGet(maskBadStr, config); // Bits to mask for bad pixels + + ppStatsFPA(stats, outRO->parent->parent->parent, view, maskBad, config); + } + + memDump("stats"); + + // Put version information into the header + pmHDU *hdu = pmHDUFromCell(outRO->parent); + if (!hdu) { + psError(PS_ERR_UNKNOWN, false, "Unable to find HDU for output."); psFree(outRO); psFree(view); - } + return false; + } + if (!hdu->header) { + hdu->header = psMetadataAlloc(); + } + ppStackVersionMetadata(hdu->header); + + // Write out the output files + filesIterateUp(config); + + psFree(outRO); + psFree(view); // Write out summary statistics Index: trunk/ppStack/src/ppStackReadout.c =================================================================== --- trunk/ppStack/src/ppStackReadout.c (revision 19283) +++ trunk/ppStack/src/ppStackReadout.c (revision 19337) @@ -12,7 +12,43 @@ //#define TESTING // Write debugging output? - -bool ppStackReadoutInitial(const pmConfig *config, pmReadout *outRO, const psArray *readouts, - const psArray *regions, const psArray *kernels) +bool ppStackReadoutInitialThread(psThreadJob *job) +{ + PS_ASSERT_THREAD_JOB_NON_NULL(job, false); + + psArray *args = job->args; // Arguments + ppStackThread *thread = args->data[0]; // Thread + pmConfig *config = args->data[1]; // Configuration + pmReadout *outRO = args->data[2]; // Output readout + psArray *subRegions = args->data[3]; // Regions for PSF-matching + psArray *subKernels = args->data[4]; // Kernels for PSF-matching + + psArray *inspect = ppStackReadoutInitial(config, outRO, thread->readouts, subRegions, subKernels); + + job->results = inspect; + thread->busy = false; + + return inspect ? true : false; +} + +bool ppStackReadoutFinalThread(psThreadJob *job) +{ + PS_ASSERT_THREAD_JOB_NON_NULL(job, false); + + psArray *args = job->args; // Arguments + ppStackThread *thread = args->data[0]; // Thread + pmConfig *config = args->data[1]; // Configuration + pmReadout *outRO = args->data[2]; // Output readout + psArray *rejected = args->data[3]; // Rejected pixels + + bool status = ppStackReadoutFinal(config, outRO, thread->readouts, rejected); // Status of operation + + thread->busy = false; + + return status; +} + + +psArray *ppStackReadoutInitial(const pmConfig *config, pmReadout *outRO, const psArray *readouts, + const psArray *regions, const psArray *kernels) { assert(config); @@ -88,4 +124,5 @@ // Save list of pixels to inspect + psArray *inspect = psArrayAlloc(num); // List of pixels to inspect for (int i = 0; i < num; i++) { pmStackData *data = stack->data[i]; // Data for this image @@ -97,7 +134,5 @@ continue; } - psMetadataAddPtr(readout->analysis, PS_LIST_TAIL, PPSTACK_INSPECT_PIXELS, - PS_DATA_PIXELS | PS_META_DUPLICATE_OK, "Pixels to inspect from initial combination", - data->inspect); + inspect->data[i] = psMemIncrRefCounter(data->inspect); } psFree(stack); @@ -105,5 +140,5 @@ sectionNum++; - return true; + return inspect; } Index: trunk/ppStack/src/ppStackThread.c =================================================================== --- trunk/ppStack/src/ppStackThread.c (revision 19337) +++ trunk/ppStack/src/ppStackThread.c (revision 19337) @@ -0,0 +1,251 @@ +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include +#include +#include + +#include "ppStack.h" + +#define THREAD_WAIT 10000 // Microseconds to wait if thread is not available + + +static void stackThreadFree(ppStackThread *thread) +{ + psFree(thread->readouts); + return; +} + +ppStackThread *ppStackThreadAlloc(psArray *readouts) +{ + ppStackThread *thread = psAlloc(sizeof(ppStackThread)); + psMemSetDeallocator(thread, (psFreeFunc)stackThreadFree); + + thread->readouts = psMemIncrRefCounter(readouts); + thread->read = false; + thread->busy = false; + thread->firstScan = 0; + thread->lastScan = 0; + return thread; +} + +static void stackThreadDataFree(ppStackThreadData *stack) +{ + psFree(stack->threads); + for (int i = 0; i < stack->imageFits->n; i++) { + psFitsClose(stack->imageFits->data[i]); + psFitsClose(stack->maskFits->data[i]); + psFitsClose(stack->weightFits->data[i]); + stack->imageFits->data[i] = stack->maskFits->data[i] = stack->weightFits->data[i] = NULL; + } + psFree(stack->imageFits); + psFree(stack->maskFits); + psFree(stack->weightFits); + return; +} + +ppStackThreadData *ppStackThreadDataSetup(const psArray *cells, const psArray *imageNames, + const psArray *maskNames, const psArray *weightNames, + const pmConfig *config) +{ + PS_ASSERT_ARRAY_NON_NULL(cells, NULL); + PS_ASSERT_ARRAYS_SIZE_EQUAL(cells, imageNames, NULL); + PS_ASSERT_ARRAYS_SIZE_EQUAL(cells, maskNames, NULL); + PS_ASSERT_ARRAYS_SIZE_EQUAL(cells, weightNames, NULL); + + ppStackThreadData *stack = psAlloc(sizeof(ppStackThreadData)); // Thread data, to return + psMemSetDeallocator(stack, (psFreeFunc)stackThreadDataFree); + + stack->lastScan = 0; + + int numInputs = cells->n; // Number of inputs + stack->imageFits = psArrayAlloc(numInputs); + stack->maskFits = psArrayAlloc(numInputs); + stack->weightFits = psArrayAlloc(numInputs); + for (int i = 0; i < numInputs; i++) { + if (!cells->data[i]) { + // Bad image + continue; + } + // Resolved names + psString imageResolved = pmConfigConvertFilename(imageNames->data[i], config, false, false); + psString maskResolved = pmConfigConvertFilename(maskNames->data[i], config, false, false); + psString weightResolved = pmConfigConvertFilename(weightNames->data[i], config, false, false); + stack->imageFits->data[i] = psFitsOpen(imageResolved, "r"); + stack->maskFits->data[i] = psFitsOpen(maskResolved, "r"); + stack->weightFits->data[i] = psFitsOpen(weightResolved, "r"); + psFree(imageResolved); + psFree(maskResolved); + psFree(weightResolved); + if (!stack->imageFits->data[i] || !stack->maskFits->data[i] || !stack->weightFits->data[i]) { + psError(PS_ERR_UNKNOWN, false, "Unable to open convolved files %s, %s, %s", + (char*)imageNames->data[i], (char*)maskNames->data[i], (char*)weightNames->data[i]); + return false; + } + } + + int numThreads = psMetadataLookupS32(NULL, config->arguments, "-threads"); // Number of threads + + // Generate readouts for each input file in each file group + psArray *threads = psArrayAlloc(numThreads + 1); // A stack for each thread + for (int i = 0; i < threads->n; i++) { + psArray *readouts = psArrayAlloc(numInputs); // Input readouts + for (int j = 0; j < numInputs; j++) { + pmCell *cell = cells->data[j]; // Cell with data + if (!cell) { + continue; + } + readouts->data[j] = pmReadoutAlloc(cell); + } + threads->data[i] = ppStackThreadAlloc(readouts); + psFree(readouts); // Drop reference + } + stack->threads = threads; + + return stack; +} + + +ppStackThread *ppStackThreadRead(bool *status, ppStackThreadData *stack, pmConfig *config, int numChunk, + int overlap) +{ + assert(status); + PS_ASSERT_PTR_NON_NULL(stack, NULL); + PS_ASSERT_PTR_NON_NULL(config, NULL); + + *status = true; + + psArray *threads = stack->threads; // Threads for reading + psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, PPSTACK_RECIPE); // Recipe + if (!recipe) { + psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to find recipe %s", PPSTACK_RECIPE); + return NULL; + } + int rows = psMetadataLookupS32(NULL, recipe, "ROWS"); // Number of rows to read per chunk + if (rows <= 0) { + psError(PS_ERR_BAD_PARAMETER_VALUE, false, "ROWS is not set in the recipe."); + return NULL; + } + + // Select an available group + while (true) { + // check for any groups which can read data + for (int j = 0; j < threads->n; j++) { + ppStackThread *thread = threads->data[j]; + if (thread->read) { + continue; + } + + thread->firstScan = stack->lastScan; // Overlap is taken care of in pmReadoutReadChunk + thread->lastScan = stack->lastScan + rows; + stack->lastScan = thread->lastScan; + + psArray *readouts = thread->readouts; + + psTimerStart ("ppStackReadChunk"); + + psTrace("ppStack", 2, "Reading data for chunk %d into group %d....\n", numChunk, j); + for (int i = 0; i < readouts->n; i++) { + pmReadout *ro = readouts->data[i]; // Input readout + if (!ro) { + continue; + } + + // override the recorded last scan + ro->thisImageScan = thread->firstScan; + ro->thisWeightScan = thread->firstScan; + ro->thisMaskScan = thread->firstScan; + ro->lastImageScan = thread->lastScan; + ro->lastMaskScan = thread->lastScan; + ro->lastWeightScan = thread->lastScan; + ro->forceScan = true; + + psFits *imageFits = stack->imageFits->data[i]; // FITS file for image + psFits *maskFits = stack->maskFits->data[i]; // FITS file for mask + psFits *weightFits = stack->weightFits->data[i]; // FITS file for weight + + + bool keepReading = false; + if (pmReadoutMore(ro, imageFits, 0, rows, config)) { + keepReading = true; + if (!pmReadoutReadChunk(ro, imageFits, 0, rows, overlap, config)) { + psError(PS_ERR_IO, false, "Unable to read chunk %d for file PPSTACK.INPUT %d", + numChunk, i); + *status = false; + return NULL; + } + } + + if (pmReadoutMoreMask(ro, maskFits, 0, rows, config)) { + keepReading = true; + if (!pmReadoutReadChunkMask(ro, maskFits, 0, rows, overlap, config)) { + psError(PS_ERR_IO, false, "Unable to read chunk %d for file PPSTACK.INPUT.MASK %d", + numChunk, i); + *status = false; + return NULL; + } + } + + if (pmReadoutMoreWeight(ro, weightFits, 0, rows, config)) { + keepReading = true; + if (!pmReadoutReadChunkWeight(ro, weightFits, 0, rows, overlap, config)) { + psError(PS_ERR_IO, false, "Unable to read chunk %d for file PPSTACK.INPUT.WEIGHT %d", + numChunk, i); + *status = false; + return NULL; + } + } + if (!keepReading) { + return NULL; + } + } + + thread->read = thread->busy = true; + return thread; + } + + // Check for threads that are ready to read + bool wait = true; + for (int j = 0; j < threads->n; j++) { + ppStackThread *thread = threads->data[j]; + if (thread->busy) { + continue; + } + thread->read = false; + wait = false; + } + if (wait) { + // No threads currently available + usleep(THREAD_WAIT); + } + } + return NULL; +} + + +void ppStackThreadInit(void) +{ + static bool threaded = false; // Are we running threaded? + if (threaded) { + psAbort("Already running threaded."); + } + threaded = true; + + { + psThreadTask *task = psThreadTaskAlloc("PPSTACK_INITIAL_COMBINE", 5); + task->function = &ppStackReadoutInitialThread; + psThreadTaskAdd(task); + psFree(task); + } + + { + psThreadTask *task = psThreadTaskAlloc("PPSTACK_FINAL_COMBINE", 4); + task->function = &ppStackReadoutFinalThread; + psThreadTaskAdd(task); + psFree(task); + } + + return; +}