Index: /trunk/ppStack/src/ppStack.h =================================================================== --- /trunk/ppStack/src/ppStack.h (revision 16692) +++ /trunk/ppStack/src/ppStack.h (revision 16693) @@ -2,5 +2,6 @@ #define PP_STACK_H -#define PPSTACK_RECIPE "PPSTACK" +#define PPSTACK_RECIPE "PPSTACK" // Name of the recipe +#define PPSTACK_REJECTED_PIXELS "PPSTACK.PIXELS" // Name of rejected pixels metadata items #include @@ -26,9 +27,16 @@ // Perform stacking on a readout -bool ppStackReadout(const pmConfig *config, // Configuration +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 + ); + +// Perform stacking on a readout +bool ppStackReadoutFinal(const pmConfig *config, // Configuration + pmReadout *outRO, // Output readout + const psArray *readouts, // Input readouts + const psArray *rejected // Array with pixels rejected in each image ); Index: /trunk/ppStack/src/ppStackLoop.c =================================================================== --- /trunk/ppStack/src/ppStackLoop.c (revision 16692) +++ /trunk/ppStack/src/ppStackLoop.c (revision 16693) @@ -413,5 +413,5 @@ #endif - if (!ppStackReadout(config, outRO, readouts, subRegions, subKernels)) { + if (!ppStackReadoutInitial(config, outRO, readouts, subRegions, subKernels)) { psError(PS_ERR_UNKNOWN, false, "Unable to stack images.\n"); psFree(readouts); @@ -432,7 +432,76 @@ } - psFree(readouts); + // Reset for the second read + // Extract the rejection lists psFree(subKernels); psFree(subRegions); + psArray *rejected = psArrayAlloc(num); // Rejected pixels + for (int i = 0; i < num; i++) { + pmReadout *ro = readouts->data[i]; // Readout of interest + pmReadoutFreeData(ro); + + psPixels *rejects = NULL; // Rejection list for this readout + psMetadataIterator *iter = psMetadataIteratorAlloc(ro->analysis, PS_LIST_HEAD, + "^" PPSTACK_REJECTED_PIXELS "$"); // Iterator + psMetadataItem *item; + while ((item = psMetadataGetAndIncrement(iter))) { + psPixels *pixels = item->data.V; // Rejected pixels + psTrace("ppStack", 5, "Adding %ld rejected pixels to image %d", pixels->n, i); + rejects = psPixelsConcatenate(rejects, pixels); + } + psFree(iter); + psTrace("ppStack", 5, "%ld rejected pixels rejected from image %d", rejects->n, i); + psMetadataRemoveKey(ro->analysis, PPSTACK_REJECTED_PIXELS); + rejected->data[i] = rejects; + } + + + // Read convolutions by chunks + more = true; + for (int numChunk = 0; more; numChunk++) { + for (int i = 0; i < num; i++) { + pmReadout *readout = readouts->data[i]; + assert(readout); + + if (!pmReadoutReadChunk(readout, imageFits->data[i], 0, numScans, 0) || + !pmReadoutReadChunkMask(readout, maskFits->data[i], 0, numScans, 0) || + !pmReadoutReadChunkWeight(readout, weightFits->data[i], 0, numScans, 0)) { + 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; + } + } + +#ifdef TESTING + { + pmReadout *ro = readouts->data[0]; + psTrace("ppStack", 1, "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]; + assert(readout); + more &= pmReadoutMore(readout, imageFits->data[i], 0, numScans); + more &= pmReadoutMoreMask(readout, maskFits->data[i], 0, numScans); + more &= pmReadoutMoreWeight(readout, weightFits->data[i], 0, numScans); + } + } + + psFree(readouts); + for (int i = 0; i < num; i++) { psFitsClose(imageFits->data[i]); Index: /trunk/ppStack/src/ppStackReadout.c =================================================================== --- /trunk/ppStack/src/ppStackReadout.c (revision 16692) +++ /trunk/ppStack/src/ppStackReadout.c (revision 16693) @@ -16,9 +16,8 @@ #define COMBINED_FILES // Write combined images? -static int sectionNum = 0; // Section number; for debugging outputs - - -bool ppStackReadout(const pmConfig *config, pmReadout *outRO, const psArray *readouts, - const psArray *regions, const psArray *kernels) + + +bool ppStackReadoutInitial(const pmConfig *config, pmReadout *outRO, const psArray *readouts, + const psArray *regions, const psArray *kernels) { assert(config); @@ -29,4 +28,5 @@ assert(readouts->n == regions->n); assert(regions->n == kernels->n); + static int sectionNum = 0; // Section number; for debugging outputs @@ -43,4 +43,105 @@ int kernelSize = psMetadataLookupS32(NULL, ppsub, "KERNEL.SIZE"); // Kernel half-size + + int num = readouts->n; // Number of inputs + psArray *stack = psArrayAlloc(num); // Array for stacking + + for (int i = 0; i < num; i++) { + pmReadout *ro = readouts->data[i]; + assert(ro); + pmFPA *fpa = ro->parent->parent->parent; // Parent FPA + + float weighting = psMetadataLookupF32(&mdok, fpa->analysis, "PPSTACK.WEIGHTING"); // Relative weight + if (!mdok || !isfinite(weighting)) { + psWarning("No WEIGHTING supplied for image %d --- set to unity.", i); + weighting = 1.0; + } + + // Ensure there is a mask, or pmStackCombine will complain + if (!ro->mask) { + ro->mask = psImageAlloc(ro->image->numCols, ro->image->numRows, PS_TYPE_MASK); + psImageInit(ro->mask, 0); + } + + stack->data[i] = pmStackDataAlloc(ro, weighting); + } + + if (!pmStackCombine(outRO, stack, maskBad, maskBlank, kernelSize, iter, combineRej, useVariance, safe)) { + psError(PS_ERR_UNKNOWN, false, "Unable to combine input readouts with rejection."); + psFree(stack); + return false; + } + +#ifdef COMBINED_FILES + { + psString name = NULL; // Name of image + psStringAppend(&name, "combined_initial_%03d.fits", sectionNum); + psFits *fits = psFitsOpen(name, "w"); + psFree(name); + psFitsWriteImage(fits, NULL, outRO->image, 0, NULL); + psFitsClose(fits); + } +#endif +#ifdef INSPECTION_FILES + for (int i = 0; i < stack->n; i++) { + pmStackData *data = stack->data[i]; // Data for this image + psImage *inspected = psPixelsToMask(NULL, data->pixels, + psRegionSet(0, outRO->image->numCols - 1, + 0, outRO->image->numRows - 1), + maskBlank); + psString name = NULL; // Name of image + psStringAppend(&name, "inspect_%03d_%03d.fits", sectionNum, i); + psFits *fits = psFitsOpen(name, "w"); + psFree(name); + psFitsWriteImage(fits, NULL, inspected, 0, NULL); + psFree(inspected); + psFitsClose(fits); + } +#endif + + // Reject pixels + for (int i = 0; i < num; i++) { + pmStackData *data = stack->data[i]; // Data for this image + pmReadout *readout = data->readout; // Readout of interest + int col0 = readout->col0, row0 = readout->row0; // Offset for readout + int numCols = readout->image->numCols, numRows = readout->image->numRows; // Size of image + + psRegion *valid = psRegionAlloc(col0, col0 + numCols, row0, row0 + numRows); // Valid region for rej + psPixels *reject = pmStackReject(data->pixels, valid, threshold, regions->data[i], + kernels->data[i]); // Pixels to reject + psFree(valid); + + psMetadataAddPtr(readout->analysis, PS_LIST_TAIL, PPSTACK_REJECTED_PIXELS, + PS_DATA_PIXELS | PS_META_DUPLICATE_OK, "Rejected pixels from initial combination", + reject); + psFree(reject); // Drop reference + } + + psFree(stack); + + sectionNum++; + + return true; +} + + + +bool ppStackReadoutFinal(const pmConfig *config, pmReadout *outRO, const psArray *readouts, + const psArray *rejected) +{ + assert(config); + assert(outRO); + assert(readouts); + assert(rejected); + assert(readouts->n == rejected->n); + + static int sectionNum = 0; // Section number; for debugging outputs + + + // Get the recipe values + bool mdok; // Status of MD lookup + psMaskType maskBad = psMetadataLookupU8(NULL, config->arguments, "MASK.BAD"); // Value to mask + psMaskType maskBlank = psMetadataLookupU8(NULL, config->arguments, "MASK.BLANK"); // Mask for blank reg. + bool useVariance = psMetadataLookupBool(&mdok, config->arguments, "VARIANCE"); // Use variance for rejection? int num = readouts->n; // Number of inputs @@ -110,78 +211,32 @@ psListAdd(cellList, PS_LIST_TAIL, ro->parent); - stack->data[i] = pmStackDataAlloc(ro, weighting); - } - - if (!pmStackCombine(outRO, stack, maskBad, maskBlank, kernelSize, iter, combineRej, useVariance, safe)) { - psError(PS_ERR_UNKNOWN, false, "Unable to combine input readouts with rejection."); - psFree(fpaList); - psFree(cellList); - psFree(stack); - return false; - } - -#ifdef COMBINED_FILES - { - psString name = NULL; // Name of image - psStringAppend(&name, "combined1_%03d.fits", sectionNum); - psFits *fits = psFitsOpen(name, "w"); - psFree(name); - psFitsWriteImage(fits, NULL, outRO->image, 0, NULL); - psFitsClose(fits); - } -#endif -#ifdef INSPECTION_FILES - for (int i = 0; i < stack->n; i++) { - pmStackData *data = stack->data[i]; // Data for this image - psImage *inspected = psPixelsToMask(NULL, data->pixels, - psRegionSet(0, data->readout->image->numCols - 1, - 0, data->readout->image->numRows - 1), - maskBlank); - psString name = NULL; // Name of image - psStringAppend(&name, "inspect_%03d_%03d.fits", sectionNum, i); - psFits *fits = psFitsOpen(name, "w"); - psFree(name); - psFitsWriteImage(fits, NULL, inspected, 0, NULL); - psFree(inspected); - psFitsClose(fits); - } -#endif - - // Reject pixels - for (int i = 0; i < num; i++) { - pmStackData *data = stack->data[i]; // Data for this image - pmReadout *readout = data->readout; // Readout of interest - int col0 = readout->col0, row0 = readout->row0; // Offset for readout - int numCols = readout->image->numCols, numRows = readout->image->numRows; // Size of image - - psRegion *valid = psRegionAlloc(col0, col0 + numCols, row0, row0 + numRows); // Valid region for rej - psPixels *reject = pmStackReject(data->pixels, valid, threshold, regions->data[i], - kernels->data[i]); // Pixels to reject - psFree(valid); - psFree(data->pixels); - data->pixels = reject; + pmStackData *data = pmStackDataAlloc(ro, weighting); + data->pixels = psMemIncrRefCounter(rejected->data[i]); + stack->data[i] = data; } #ifdef REJECTION_FILES - for (int i = 0; i < stack->n; i++) { - pmStackData *data = stack->data[i]; // Data for this image - if (!data) { - continue; - } - psImage *rejected = psPixelsToMask(NULL, data->pixels, - psRegionSet(0, data->readout->image->numCols - 1, - 0, data->readout->image->numRows - 1), - maskBlank); - psString name = NULL; // Name of image - psStringAppend(&name, "reject_%03d_%03d.fits", sectionNum, i); - psFits *fits = psFitsOpen(name, "w"); - psFree(name); - psFitsWriteImage(fits, NULL, rejected, 0, NULL); - psFree(rejected); - psFitsClose(fits); - } -#endif - - if (!pmStackCombine(outRO, stack, maskBad, maskBlank, kernelSize, 0, combineRej, useVariance, false)) { + if (sectionNum == 0) { + for (int i = 0; i < stack->n; i++) { + pmStackData *data = stack->data[i]; // Data for this image + if (!data) { + continue; + } + psImage *reject = psPixelsToMask(NULL, data->pixels, + psRegionSet(0, outRO->image->numCols - 1, + 0, outRO->image->numRows - 1), + maskBlank); + psString name = NULL; // Name of image + psStringAppend(&name, "reject_%03d.fits", i); + psFits *fits = psFitsOpen(name, "w"); + psFree(name); + psFitsWriteImage(fits, NULL, reject, 0, NULL); + psFree(reject); + psFitsClose(fits); + } + } +#endif + + if (!pmStackCombine(outRO, stack, maskBad, maskBlank, 0, 0, NAN, useVariance, false)) { psError(PS_ERR_UNKNOWN, false, "Unable to combine input readouts."); psFree(fpaList); @@ -194,5 +249,5 @@ { psString name = NULL; // Name of image - psStringAppend(&name, "combined2_%03d.fits", sectionNum); + psStringAppend(&name, "combined_final_%03d.fits", sectionNum); psFits *fits = psFitsOpen(name, "w"); psFree(name); @@ -219,6 +274,4 @@ psFree(stack); - sectionNum++; - return true; }