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