Index: trunk/ppStack/src/ppStackReadout.c =================================================================== --- trunk/ppStack/src/ppStackReadout.c (revision 15849) +++ trunk/ppStack/src/ppStackReadout.c (revision 16605) @@ -10,16 +10,25 @@ #include "ppStack.h" -#define ARRAY_BUFFER 16 // Number to add to array at a time #define WCS_TOLERANCE 0.001 // Tolerance for WCS -//#define NO_CONVOLUTION // Don't perform convolutions? -//#define CONVOLUTION_FILES // Write convolutions? //#define REJECTION_FILES // Write rejection mask? //#define INSPECTION_FILES // Write inspection mask? -bool ppStackReadout(pmConfig *config, const pmFPAview *view) +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) { + assert(config); + assert(outRO); + assert(readouts); + assert(regions); + assert(kernels); + assert(readouts->n == regions->n); + assert(regions->n == kernels->n); + + // Get the recipe values - bool mdok; // Status of MD lookup int iter = psMetadataLookupS32(NULL, config->arguments, "ITER"); // Rejection iterations float combineRej = psMetadataLookupF32(NULL, config->arguments, "COMBINE.REJ"); // Combination threshold @@ -27,196 +36,41 @@ psMaskType maskBlank = psMetadataLookupU8(NULL, config->arguments, "MASK.BLANK"); // Mask for blank reg. float threshold = psMetadataLookupF32(NULL, config->arguments, "THRESHOLD.MASK"); // Threshold for mask deconvolution - bool photometry = psMetadataLookupBool(&mdok, config->arguments, "PHOTOMETRY"); // Perform photometry? - int psfInstances = psMetadataLookupS32(&mdok, config->arguments, "PSF.INSTANCES"); // Number of instances for PSF - float psfRadius = psMetadataLookupF32(NULL, config->arguments, "PSF.RADIUS"); // Radius for PSF - const char *psfModel = psMetadataLookupStr(NULL, config->arguments, "PSF.MODEL"); // Model for PSF - int psfOrder = psMetadataLookupS32(&mdok, config->arguments, "PSF.ORDER"); // Spatial order for PSF - - // Get the output target - pmCell *outCell = pmFPAfileThisCell(config->files, view, "PPSTACK.OUTPUT"); // Output cell - pmReadout *outRO = pmReadoutAlloc(outCell); // Output readout - pmFPA *outFPA = outCell->parent->parent; // Output FPA - - int num = psMetadataLookupS32(NULL, config->arguments, "INPUTS.NUM"); // Number of inputs - assert(num > 0); - - // Get the input sources - psArray *stack = psArrayAllocEmpty(ARRAY_BUFFER); // The stack of inputs - pmReadout *sources = pmFPAfileThisReadout(config->files, view, "PPSTACK.SOURCES"); // Sources for stamps - psMetadataIterator *fileIter = psMetadataIteratorAlloc(config->files, PS_LIST_HEAD, - "^PPSTACK.INPUT$"); // Iterator over input files - psMetadataItem *fileItem; // Item from iteration - int fileNum = 0; // Number of file - psList *psfList = psListAlloc(NULL); // List of PSFs for PSF envelope - pmPSF *outPSF = NULL; // Ouptut PSF - int numCols = 0, numRows = 0; // Size of image - while ((fileItem = psMetadataGetAndIncrement(fileIter))) { - assert(fileItem->type == PS_DATA_UNKNOWN); - pmFPAfile *inputFile = fileItem->data.V; // An input file - pmReadout *ro = pmFPAviewThisReadout(view, inputFile->fpa); // Corresponding readout - pmCell *cell = ro->parent; // The cell - pmChip *chip = cell->parent; // The chip: holds the PSF - - bool mdok; // Status of MD lookup - pmPSF *psf = psMetadataLookupPtr(&mdok, chip->analysis, "PSPHOT.PSF"); - if (mdok && psf) { - psListAdd(psfList, PS_LIST_TAIL, psf); - if (numCols == 0 && numRows == 0) { - numCols = ro->image->numCols; - numRows = ro->image->numRows; - } - } - } - - bool convolve = false; // Convolve the input images? - if (psfList->n > 0) { - psArray *psfArray = psListToArray(psfList); // Array of PSFs - outPSF = pmPSFEnvelope(numCols, numRows, psfArray, psfInstances, psfRadius, psfModel, - psfOrder, psfOrder); - psFree(psfArray); - if (!outPSF) { - psError(PS_ERR_UNKNOWN, false, "Unable to determine output PSF."); - // XXX Free stuff - return false; - } - convolve = true; - PS_ASSERT_PTR_NON_NULL(sources, false); - } - - // Iterate through again to get the convolved images (or not) - - psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN); // Statistics - psRandom *rng = psRandomAlloc(0, PS_RANDOM_TAUS); // Random number generator + + int num = readouts->n; // Number of inputs + psArray *stack = psArrayAlloc(num); // Array for stacking + + pmCell *outCell = outRO->parent; // Output cell + pmChip *outChip = outCell->parent; // Output chip + pmFPA *outFPA = outChip->parent; // Output FPA + float totExposure = 0.0; // Total exposure time psList *fpaList = psListAlloc(NULL); // List of input FPAs, for concept averaging psList *cellList = psListAlloc(NULL); // List of input cells, for concept averaging - - psMetadataIteratorSet(fileIter, PS_LIST_HEAD); - fileNum = 0; - while ((fileItem = psMetadataGetAndIncrement(fileIter))) { - assert(fileItem->type == PS_DATA_UNKNOWN); - pmFPAfile *inputFile = fileItem->data.V; // An input file - pmReadout *ro = pmFPAviewThisReadout(view, inputFile->fpa); // Corresponding readout + for (int i = 0; i < num; i++) { + pmReadout *ro = readouts->data[i]; + assert(ro); + pmFPA *fpa = ro->parent->parent->parent; // Parent FPA bool mdok; // Status of MD lookup - float weighting = psMetadataLookupF32(&mdok, inputFile->fpa->analysis, - "PPSTACK.WEIGHTING"); // Relative weighting + float weighting = psMetadataLookupF32(&mdok, fpa->analysis, "PPSTACK.WEIGHTING"); // Relative weight if (!mdok || !isfinite(weighting)) { - psWarning("No WEIGHTING supplied for image %d --- set to unity.", fileNum); + psWarning("No WEIGHTING supplied for image %d --- set to unity.", i); weighting = 1.0; } - float scale = psMetadataLookupF32(&mdok, inputFile->fpa->analysis, "PPSTACK.SCALE"); // Rel. scale - if (!mdok || !isfinite(scale)) { - psWarning("No SCALE supplied for image %d --- set to unity.", fileNum); - scale = 1.0; - } float exposure = psMetadataLookupF32(NULL, ro->parent->concepts, "CELL.EXPOSURE"); // Exposure time -#if 0 if (!isfinite(exposure)) { psError(PS_ERR_BAD_PARAMETER_VALUE, true, "CELL.EXPOSURE is not set for input file %ld", stack->n); - psFree(stats); - psFree(rng); - psFree(fileIter); psFree(fpaList); psFree(cellList); - psFree(outRO); psFree(stack); return false; } -#endif totExposure += exposure; // Total exposure time - // Generate convolved version of input - pmReadout *convolved; - if (convolve) { - convolved = pmReadoutAlloc(NULL); // Convolved version of input readout - // Background subtraction, scaling and normalisation is performed automatically by the image - // matching - if (!ppStackMatch(convolved, ro, sources, outPSF, config)) { - psWarning("Unable to match image %d --- ignoring.", fileNum); - psErrorClear(); - psFree(convolved); - // XXX Free the bad image so it's not taking up good memory! - continue; - } - -#ifdef CONVOLUTION_FILES - if (convolved->image) { - psString name = NULL; // Name of image - psStringAppend(&name, "convolved%03d_image.fits", fileNum); - psFits *fits = psFitsOpen(name, "w"); - psFree(name); - psFitsWriteImage(fits, NULL, convolved->image, 0, NULL); - psFitsClose(fits); - } - if (convolved->mask) { - psString name = NULL; // Name of image - psStringAppend(&name, "convolved%03d_mask.fits", fileNum); - psFits *fits = psFitsOpen(name, "w"); - psFree(name); - psFitsWriteImage(fits, NULL, convolved->mask, 0, NULL); - psFitsClose(fits); - } - if (convolved->weight) { - psString name = NULL; // Name of image - psStringAppend(&name, "convolved%03d_weight.fits", fileNum); - psFits *fits = psFitsOpen(name, "w"); - psFree(name); - psFitsWriteImage(fits, NULL, convolved->weight, 0, NULL); - psFitsClose(fits); - } -#endif // CONVOLUTION_FILES - - } else { - // No convolution --- just use the unconvolved images as the "convolved" images, with some tweaks - convolved = psMemIncrRefCounter(ro); - sources = NULL; - - // Brain-dead background subtraction - if (!psImageBackground(stats, NULL, ro->image, ro->mask, maskBad, rng)) { - psError(PS_ERR_UNKNOWN, false, "Unable to get image background on image %ld", stack->n); - psFree(stats); - psFree(rng); - psFree(fileIter); - psFree(fpaList); - psFree(cellList); - psFree(stack); - psFree(outRO); - psFree(convolved); - return false; - } - psTrace("ppStack", 3, "Background for image %d is %f\n", fileNum, stats->robustMedian); - (void)psBinaryOp(ro->image, ro->image, "-", psScalarAlloc(stats->robustMedian, PS_TYPE_F32)); - - // Apply scaling - (void)psBinaryOp(ro->image, ro->image, "*", psScalarAlloc(1.0 / scale, PS_TYPE_F32)); - - // Normalise each input by the exposure time - float exposure = psMetadataLookupF32(NULL, ro->parent->concepts, "CELL.EXPOSURE");// Exposure time - if (!isfinite(exposure)) { - psError(PS_ERR_BAD_PARAMETER_VALUE, true, - "CELL.EXPOSURE is not set for input file %ld", stack->n); - psFree(stats); - psFree(rng); - psFree(fileIter); - psFree(fpaList); - psFree(cellList); - psFree(outRO); - psFree(convolved); - psFree(stack); - return false; - } - - (void)psBinaryOp(ro->image, ro->image, "/", psScalarAlloc(exposure, PS_TYPE_F32)); - if (ro->weight) { - (void)psBinaryOp(ro->weight, ro->weight, "/", psScalarAlloc(exposure, PS_TYPE_F32)); - } - } - - if (fileNum == 0) { +#if 0 + if (i == 0) { // Copy astrometry over - pmFPA *fpa = ro->parent->parent->parent; // Template FPA pmHDU *hdu = fpa->hdu; // Template HDU pmHDU *outHDU = outFPA->hdu; // Output HDU @@ -224,5 +78,5 @@ psWarning("Unable to find HDU at FPA level to copy astrometry."); } else { - if (!pmAstromReadWCS(outFPA, outCell->parent, hdu->header, 1.0)) { + if (!pmAstromReadWCS(outFPA, outChip, hdu->header, 1.0)) { psErrorClear(); psWarning("Unable to read WCS astrometry from input FPA."); @@ -231,5 +85,5 @@ outHDU->header = psMetadataAlloc(); } - if (!pmAstromWriteWCS(outHDU->header, outFPA, outCell->parent, WCS_TOLERANCE)) { + if (!pmAstromWriteWCS(outHDU->header, outFPA, outChip, WCS_TOLERANCE)) { psErrorClear(); psWarning("Unable to write WCS astrometry to output FPA."); @@ -238,36 +92,16 @@ } } - - // Don't need the original any more! - psFree(ro); +#endif // Ensure there is a mask, or pmStackCombine will complain - if (!convolved->mask) { - convolved->mask = psImageAlloc(convolved->image->numCols, convolved->image->numRows, - PS_TYPE_MASK); - psImageInit(convolved->mask, 0); - } - - psListAdd(fpaList, PS_LIST_TAIL, inputFile->fpa); + if (!ro->mask) { + ro->mask = psImageAlloc(ro->image->numCols, ro->image->numRows, PS_TYPE_MASK); + psImageInit(ro->mask, 0); + } + + psListAdd(fpaList, PS_LIST_TAIL, fpa); psListAdd(cellList, PS_LIST_TAIL, ro->parent); - pmStackData *data = pmStackDataAlloc(convolved, weighting); // Data to stack - psFree(convolved); - psArrayAdd(stack, ARRAY_BUFFER, data); - psFree(data); // Drop reference - - fileNum++; - } - psFree(fileIter); - psFree(stats); - psFree(rng); - - if (fileNum == 0) { - psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Not enough good files to combine."); - psFree(fpaList); - psFree(cellList); - psFree(stack); - psFree(outRO); - return false; + stack->data[i] = pmStackDataAlloc(ro, weighting); } @@ -277,5 +111,4 @@ psFree(cellList); psFree(stack); - psFree(outRO); return false; } @@ -283,5 +116,8 @@ #ifdef INSPECTION_FILES { - psFits *fits = psFitsOpen("combined.fits", "w"); + psString name = NULL; // Name of image + psStringAppend(&name, "combined_%03d.fits", sectionNum); + psFits *fits = psFitsOpen(name, "w"); + psFree(name); psFitsWriteImage(fits, NULL, outRO->image, 0, NULL); psFitsClose(fits); @@ -291,9 +127,9 @@ pmStackData *data = stack->data[i]; // Data for this image psImage *inspected = psPixelsToMask(NULL, data->pixels, - psRegionSet(0, data->readout->image->numCols, - 0, data->readout->image->numRows), + psRegionSet(0, data->readout->image->numCols - 1, + 0, data->readout->image->numRows - 1), maskBlank); psString name = NULL; // Name of image - psStringAppend(&name, "inspect%03d.fits", i); + psStringAppend(&name, "inspect_%03d_%03d.fits", sectionNum, i); psFits *fits = psFitsOpen(name, "w"); psFree(name); @@ -304,43 +140,17 @@ #endif - for (int i = 0; i < stack->n; i++) { + // Reject pixels + for (int i = 0; i < num; i++) { pmStackData *data = stack->data[i]; // Data for this image - pmReadout *readout = data->readout; // Readout for this image - - // Extract the regions and solutions used in the image matching - psArray *regions = psArrayAllocEmpty(ARRAY_BUFFER); // Array of regions - { - psMetadataIterator *iter = psMetadataIteratorAlloc(readout->analysis, PS_LIST_HEAD, - "^SUBTRACTION.REGION$"); // Iterator - psMetadataItem *item = NULL;// Item from iteration - while ((item = psMetadataGetAndIncrement(iter))) { - assert(item->type == PS_DATA_REGION); - regions = psArrayAdd(regions, ARRAY_BUFFER, item->data.V); - } - psFree(iter); - } - psArray *solutions = psArrayAllocEmpty(ARRAY_BUFFER); // Array of solutions - { - psMetadataIterator *iter = psMetadataIteratorAlloc(readout->analysis, PS_LIST_HEAD, - "^SUBTRACTION.SOLUTION$"); // Iterator - psMetadataItem *item = NULL;// Item from iteration - while ((item = psMetadataGetAndIncrement(iter))) { - assert(item->type == PS_DATA_VECTOR); - solutions = psArrayAdd(solutions, ARRAY_BUFFER, item->data.V); - } - psFree(iter); - } - assert(regions->n == solutions->n); - - pmSubtractionKernels *kernels = psMetadataLookupPtr(NULL, readout->analysis, - "SUBTRACTION.KERNEL"); // Kernels - - psPixels *reject = pmStackReject(data->pixels, threshold, regions, - solutions, kernels); // List of pixels to reject + 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; - - psFree(solutions); - psFree(regions); } @@ -352,9 +162,9 @@ } psImage *rejected = psPixelsToMask(NULL, data->pixels, - psRegionSet(0, data->readout->image->numCols, - 0, data->readout->image->numRows), + psRegionSet(0, data->readout->image->numCols - 1, + 0, data->readout->image->numRows - 1), maskBlank); psString name = NULL; // Name of image - psStringAppend(&name, "reject%03d.fits", i); + psStringAppend(&name, "reject_%03d_%03d.fits", sectionNum, i); psFits *fits = psFitsOpen(name, "w"); psFree(name); @@ -370,14 +180,5 @@ psFree(cellList); psFree(stack); - psFree(outRO); return false; - } - - if (!convolve) { - // Restore image to counts using the total exposure time - (void)psBinaryOp(outRO->image, outRO->image, "*", psScalarAlloc(totExposure, PS_TYPE_F32)); - if (outRO->weight) { - (void)psBinaryOp(outRO->weight, outRO->weight, "*", psScalarAlloc(totExposure, PS_TYPE_F32)); - } } @@ -397,19 +198,7 @@ psFree(cellList); - if (photometry) { - pmFPAfile *photFile = psMetadataLookupPtr(NULL, config->files, "PSPHOT.INPUT"); - pmFPACopy(photFile->fpa, outRO->parent->parent->parent); - - if (!psphotReadout(config, view)) { - psWarning("Unable to perform photometry on stacked image."); - psErrorStackPrint(stderr, "Error stack from photometry:"); - psErrorClear(); - } - - pmFPAfileActivate(config->files, false, "PSPHOT.INPUT"); - } - psFree(stack); - psFree(outRO); // Drop reference + + sectionNum++; return true;