Index: /trunk/ppStack/src/Makefile.am =================================================================== --- /trunk/ppStack/src/Makefile.am (revision 16604) +++ /trunk/ppStack/src/Makefile.am (revision 16605) @@ -9,5 +9,7 @@ ppStackCamera.c \ ppStackLoop.c \ + ppStackPSF.c \ ppStackReadout.c \ + ppStackPhotometry.c \ ppStackVersion.c \ ppStackMatch.c Index: /trunk/ppStack/src/ppStack.h =================================================================== --- /trunk/ppStack/src/ppStack.h (revision 16604) +++ /trunk/ppStack/src/ppStack.h (revision 16605) @@ -19,7 +19,22 @@ ); +// Determine target PSF for input images +pmPSF *ppStackPSF(const pmConfig *config, // Configuration + int numCols, int numRows, // Size of image + const psList *list // List of input PSFs + ); + // Perform stacking on a readout -bool ppStackReadout(pmConfig *config, // Configuration - const pmFPAview *view // View for readout +bool ppStackReadout(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 photometry on stack +bool ppStackPhotometry(pmConfig *config, // Configuration + const pmReadout *readout, // Readout to be photometered + const pmFPAview *view // View to readout ); @@ -35,6 +50,7 @@ /// Convolve image to match specified seeing -bool ppStackMatch(pmReadout *output, ///< Convolved readout - const pmReadout *input, ///< Readout to be convolved +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 pmReadout *sourcesRO, ///< Readout with sources const pmPSF *psf, ///< Target PSF Index: /trunk/ppStack/src/ppStackArguments.c =================================================================== --- /trunk/ppStack/src/ppStackArguments.c (revision 16604) +++ /trunk/ppStack/src/ppStackArguments.c (revision 16605) @@ -116,4 +116,5 @@ psMetadataAddU8(arguments, PS_LIST_TAIL, "-mask-blank", 0, "Mask value for blank region", 0); psMetadataAddF32(arguments, PS_LIST_TAIL, "-threshold-mask", 0, "Threshold for mask deconvolution", NAN); + psMetadataAddS32(arguments, PS_LIST_TAIL, "-rows", 0, "Rows to read at once", 128); psMetadataAddBool(arguments, PS_LIST_TAIL, "-photometry", 0, "Do photometry on stacked image?", false); psMetadataAddS32(arguments, PS_LIST_TAIL, "-psf-instances", 0, "Number of instances for PSF generation", 5); @@ -160,4 +161,5 @@ VALUE_ARG_RECIPE_MASK("-mask-blank", "MASK.BLANK"); VALUE_ARG_RECIPE_FLOAT("-threshold-mask", "THRESHOLD.MASK", F32); + VALUE_ARG_RECIPE_INT("-rows", "ROWS", S32, 0); VALUE_ARG_RECIPE_INT("-psf-instances", "PSF.INSTANCES", S32, 0); Index: /trunk/ppStack/src/ppStackCamera.c =================================================================== --- /trunk/ppStack/src/ppStackCamera.c (revision 16604) +++ /trunk/ppStack/src/ppStackCamera.c (revision 16605) @@ -10,4 +10,62 @@ #include "ppStack.h" + + +#if 0 +// Define an output convolved image file +static pmFPAfile *defineOutputConvolved(const char *name, // FPA file name + pmFPA *fpa, // FPA to bind + const pmConfig *config, // Configuration + pmFPAfileType type // Expected type + ) +{ + pmFPAfile *file = pmFPAfileDefineOutput(config, fpa, name); + if (!file) { + psError(PS_ERR_UNKNOWN, false, "Unable to define output convolved file %s", name); + return NULL; + } + if (file->type != PM_FPA_FILE_IMAGE) { + psError(PS_ERR_IO, true, "PPSTACK.OUTCONV is not of type %s", pmFPAfileStringFromType(type)); + return NULL; + } + + return file; +} + +// Define an input convolved image file, using the output as a basis +static pmFPAfile *defineInputConvolved(const char *inputName, // Input FPA file name + pmFPAfile *outFile, // Corresponding output FPA file + pmConfig *config, // Configuration + pmFPAfileType type // Expected type + ) +{ + pmFPAview *view = pmFPAviewAlloc(0); // View into sky cells + view->chip = view->cell = view->readout = 0; + + psString imageName = pmFPANameFromRule(outFile->filerule, outFile->fpa, view); + psArray *imageNames = psArrayAlloc(1); + imageNames->data[0] = imageName; + psMetadataAddArray(config->arguments, PS_LIST_TAIL, "INCONV.FILENAMES", PS_META_REPLACE, + "Filenames of input convolved image files", imageNames); + psFree(imageNames); + bool found = false; // Found the file? + pmFPAfile *imageFile = pmFPAfileDefineFromArgs(&found, config, "PPSTACK.INCONV", + "INCONV.FILENAMES"); + psMetadataRemoveKey(config->arguments, "INCONV.FILENAMES"); + if (!imageFile || !found) { + psError(PS_ERR_UNKNOWN, false, "Unable to define %s file", inputName); + return NULL; + } + if (imageFile->type != type) { + psError(PS_ERR_IO, true, "PPSTACK.INCONV is not of type %s", + pmFPAfileStringFromType(type)); + return NULL; + } + + return imageFile; +} +#endif + + bool ppStackCamera(pmConfig *config) @@ -138,4 +196,28 @@ } } + +#if 0 + // Output convolved files + pmFPAfile *outconvImage = defineOutputConvolved("PPSTACK.OUTCONV", imageFile->fpa, config, + PM_FPA_FILE_IMAGE); + pmFPAfile *outconvMask = defineOutputConvolved("PPSTACK.OUTCONV.MASK", imageFile->fpa, config, + PM_FPA_FILE_MASK); + pmFPAfile *outconvWeight = defineOutputConvolved("PPSTACK.OUTCONV.WEIGHT", imageFile->fpa, config, + PM_FPA_FILE_WEIGHT); + if (!outconvImage || !outconvMask || !outconvWeight) { + return false; + } + + // Input convolved files + pmFPAfile *inconvImage = defineInputConvolved("PPSTACK.INCONV", outconvImage, config, + PM_FPA_FILE_IMAGE); + pmFPAfile *inconvMask = defineInputConvolved("PPSTACK.INCONV.MASK", outconvMask, config, + PM_FPA_FILE_MASK); + pmFPAfile *inconvWeight = defineInputConvolved("PPSTACK.INCONV.WEIGHT", outconvWeight, config, + PM_FPA_FILE_WEIGHT); + if (!inconvImage || !inconvMask || !inconvWeight) { + return false; + } +#endif psMetadataAddF32(imageFile->fpa->analysis, PS_LIST_TAIL, "PPSTACK.WEIGHTING", 0, @@ -241,4 +323,5 @@ } + // Output PSF if (havePSFs) { Index: /trunk/ppStack/src/ppStackLoop.c =================================================================== --- /trunk/ppStack/src/ppStackLoop.c (revision 16604) +++ /trunk/ppStack/src/ppStackLoop.c (revision 16605) @@ -11,6 +11,176 @@ #include "ppStack.h" +// Here follows lists of files for activation/deactivation at various stages. Each must be NULL-terminated. + +#if 0 +// All files in the system +static char *allFiles[] = { "PPSTACK.INPUT", "PPSTACK.INPUT.MASK", "PPSTACK.INPUT.WEIGHT", + "PPSTACK.OUTPUT", "PPSTACK.OUTPUT.MASK", "PPSTACK.OUTPUT.WEIGHT", + "PSPHOT.PSF.SAVE", "PSPHOT.PSF.LOAD", "PPSTACK.SOURCES", + "PSPHOT.OUTPUT", "PSPHOT.RESID", "PSPHOT.BACKMDL", "PSPHOT.BACKMDL.STDEV", + "PSPHOT.BACKGND", "PSPHOT.BACKSUB", "SOURCE.PLOT.MOMENTS", + "SOURCE.PLOT.PSFMODEL", "SOURCE.PLOT.APRESID", "PSPHOT.INPUT.CMF", + 0 }; +#endif + +// Files required in preparation for convolution +static char *prepareFiles[] = { "PSPHOT.PSF.LOAD", "PPSTACK.SOURCES", 0 }; + +// Files required for the convolution +static char *convolveFiles[] = { "PPSTACK.INPUT", "PPSTACK.INPUT.MASK", "PPSTACK.INPUT.WEIGHT", 0 }; + +// Output files for the combination +static char *combineFiles[] = { "PPSTACK.OUTPUT", "PPSTACK.OUTPUT.MASK", "PPSTACK.OUTPUT.WEIGHT", 0 }; + +// Files for photometry +static char *photFiles[] = { "PSPHOT.OUTPUT", "PSPHOT.RESID", "PSPHOT.BACKMDL", + "PSPHOT.BACKMDL.STDEV", "PSPHOT.BACKGND", "PSPHOT.BACKSUB", + "SOURCE.PLOT.MOMENTS", "SOURCE.PLOT.PSFMODEL", "SOURCE.PLOT.APRESID", + "PSPHOT.INPUT.CMF", 0 }; + +//#define CONVOLVED_ALREADY // Already have the convolution products --- testing + + +// Activate/deactivate a list of files +static void fileActivation(pmConfig *config, // Configuration + char **files, // Files to turn on/off + bool state // Activation state + ) +{ + assert(config); + assert(files); + + for (int i = 0; files[i] != NULL; i++) { + pmFPAfileActivate(config->files, state, files[i]); + } + return; +} + +// Activate/deactivate a single element for a list; return array of files +static psArray *fileActivationSingle(pmConfig *config, // Configuration + char **files, // Files to turn on/off + bool state, // Activation state + int num // Number of file in sequence + ) +{ + assert(config); + assert(files); + + psList *list = psListAlloc(NULL); // List of files + + for (int i = 0; files[i] != NULL; i++) { + pmFPAfile *file = pmFPAfileActivateSingle(config->files, state, files[i], num); // Activated file + psListAdd(list, PS_LIST_TAIL, file); + } + + psArray *array = psListToArray(list); + psFree(list); + + return array; +} + +#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 + ) +{ + assert(config); + + pmFPAview *view = pmFPAviewAlloc(0);// Pointer into FPA hierarchy + if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) { + return NULL; + } + view->chip = 0; + if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) { + return NULL; + } + view->cell = 0; + if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) { + return NULL; + } + view->readout = 0; + if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) { + return NULL; + } + return view; +} + +// Iterate up the hierarchy, writing files; we can get away with this because we're working on skycells +static bool filesIterateUp(pmConfig *config // Configuration + ) +{ + assert(config); + + pmFPAview *view = pmFPAviewAlloc(0);// Pointer into FPA hierarchy + view->chip = view->cell = view->readout = 0; + if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) { + return false; + } + view->readout = -1; + if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) { + return false; + } + view->cell = -1; + if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) { + return false; + } + view->chip = -1; + if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) { + return false; + } + return true; +} + +#ifndef CONVOLVED_ALREADY +// Write an image to a FITS file +static bool writeImage(const char *name, // Name of image + psMetadata *header, // Header + const psImage *image // Image + ) +{ + assert(name); + assert(image); + + psFits *fits = psFitsOpen(name, "w"); + if (!fits) { + psError(PS_ERR_IO, false, "Unable to open FITS file to write image."); + return false; + } + if (!psFitsWriteImage(fits, header, image, 0, NULL)) { + psError(PS_ERR_IO, false, "Unable to write FITS image."); + return false; + } + psFitsClose(fits); + return true; +} +#endif + + bool ppStackLoop(pmConfig *config) { + assert(config); + psMaskType maskBlank = psMetadataLookupU8(NULL, config->arguments, "MASK.BLANK"); // Mask for blank reg. @@ -32,82 +202,267 @@ } - pmFPAfile *input = psMetadataLookupPtr(NULL, config->files, "PPSTACK.INPUT"); // Token input file - if (!input) { - psError(PS_ERR_UNEXPECTED_NULL, false, "Can't find input data!\n"); - return false; - } pmFPAfile *output = psMetadataLookupPtr(NULL, config->files, "PPSTACK.OUTPUT"); // Output file if (!output) { - psError(PS_ERR_UNEXPECTED_NULL, false, "Can't find output data!\n"); - return false; - } - - pmFPAview *view = pmFPAviewAlloc(0); // Pointer into FPA hierarchy - pmHDU *lastHDU = NULL; // Last HDU that was updated - - // Iterate over the FPA hierarchy - if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) { - return false; - } - - pmChip *chip; // Chip of interest - while ((chip = pmFPAviewNextChip(view, input->fpa, 1)) != NULL) { - if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) { - return false; - } - - pmCell *cell; // Cell of interest - while ((cell = pmFPAviewNextCell(view, input->fpa, 1)) != NULL) { - if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) { + psError(PS_ERR_UNEXPECTED_NULL, false, "Can't find output data!"); + return false; + } + int num = psMetadataLookupS32(NULL, config->arguments, "INPUTS.NUM"); // Number of inputs + int numScans = psMetadataLookupS32(NULL, config->arguments, "ROWS"); // Number of scans to read at once + psMetadata *ppsub = psMetadataLookupMetadata(NULL, config->recipes, "PPSUB"); // PPSUB recipe + int overlap = 2 * psMetadataLookupS32(NULL, ppsub, "KERNEL.SIZE"); // Overlap by kernel size between consecutive scans + + // Preparation iteration: Load the sources, and get a target PSF model + pmReadout *sources = NULL; // Readout with sources to use for PSF matching + pmPSF *targetPSF = NULL; // Target PSF + { + pmFPAfileActivate(config->files, false, NULL); + fileActivation(config, prepareFiles, true); + pmFPAview *view = filesIterateDown(config); + if (!view) { + return false; + } + + // We want to hang on to the 'sources' even when its host FPA is blown away + sources = psMemIncrRefCounter(pmFPAfileThisReadout(config->files, view, "PPSTACK.SOURCES")); + if (!sources) { + psError(PS_ERR_UNKNOWN, true, "Unable to find sources."); + psFree(view); + return false; + } + + // Generate list of PSFs + psMetadataIterator *fileIter = psMetadataIteratorAlloc(config->files, PS_LIST_HEAD, + "^PPSTACK.INPUT$"); // Iterator + psMetadataItem *fileItem; // Item from iteration + psList *psfList = psListAlloc(NULL); // List of PSFs for PSF envelope + 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 + pmChip *chip = pmFPAviewThisChip(view, inputFile->fpa); // The chip: holds the PSF + pmPSF *psf = psMetadataLookupPtr(NULL, chip->analysis, "PSPHOT.PSF"); // PSF + if (!psf) { + psError(PS_ERR_UNKNOWN, false, "Unable to find PSF."); + psFree(view); + psFree(sources); + psFree(fileIter); + psFree(psfList); return false; } - - // Put version information into the header + psListAdd(psfList, PS_LIST_TAIL, psf); + + pmCell *cell = pmFPAviewThisCell(view, inputFile->fpa); // Cell of interest pmHDU *hdu = pmHDUFromCell(cell); - if (hdu && hdu != lastHDU) { - if (!hdu->header) { - hdu->header = psMetadataAlloc(); - } - ppStackVersionMetadata(hdu->header); - lastHDU = hdu; - } - - pmReadout *readout; // Readout of interest - while ((readout = pmFPAviewNextReadout(view, input->fpa, 1))) { - if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) { + assert(hdu && hdu->header); + int naxis1 = psMetadataLookupS32(NULL, hdu->header, "NAXIS1"); // Number of columns + int naxis2 = psMetadataLookupS32(NULL, hdu->header, "NAXIS2"); // Number of rows + if (naxis1 <= 0 || naxis2 <= 0) { + psError(PS_ERR_UNKNOWN, false, "Unable to determine size of image from PSF."); + psFree(view); + psFree(sources); + psFree(fileIter); + psFree(psfList); + return false; + } + if (numCols == 0 && numRows == 0) { + numCols = naxis1; + numRows = naxis2; + } + } + psFree(fileIter); + psFree(view); + + targetPSF = ppStackPSF(config, numCols, numRows, psfList); + psFree(psfList); + if (!targetPSF) { + psError(PS_ERR_UNKNOWN, false, "Unable to determine output PSF."); + psFree(sources); + return false; + } + + filesIterateUp(config); + } + + const char *suffix = "conv.fits"; // Suffix for convolved images; ultimately this will be from recipe + const char *outName = psMetadataLookupStr(NULL, config->arguments, "OUTPUT"); // Output root + assert(outName); + psArray *imageNames = psArrayAlloc(num); + psArray *maskNames = psArrayAlloc(num); + psArray *weightNames = psArrayAlloc(num); + for (int i = 0; i < num; i++) { + psString imageName = NULL, maskName = NULL, weightName = NULL; // Names for convolved images + psStringAppend(&imageName, "%s.im-%d.%s", outName, i, suffix); + psStringAppend(&maskName, "%s.mk-%d.%s", outName, i, suffix); + psStringAppend(&weightName, "%s.wt-%d.%s", outName, i, suffix); + imageNames->data[i] = imageName; + maskNames->data[i] = maskName; + weightNames->data[i] = weightName; + } + + // Generate convolutions and write them to disk + psArray *cells = psArrayAlloc(num); // Cells for convolved images --- a handle for reading again + psArray *subKernels = psArrayAlloc(num); // Subtraction kernels --- required in the stacking + psArray *subRegions = psArrayAlloc(num); // Subtraction regions --- required in the stacking + for (int i = 0; i < num; i++) { + pmFPAfileActivate(config->files, false, NULL); + psArray *files = fileActivationSingle(config, convolveFiles, true, i); + pmFPAview *view = filesIterateDown(config); + if (!view) { + psFree(sources); + psFree(targetPSF); + return false; + } + + pmReadout *readout = pmFPAviewThisReadout(view, ((pmFPAfile*)files->data[0])->fpa); // Input readout + psFree(view); + +#ifndef CONVOLVED_ALREADY + // Background subtraction, scaling and normalisation is performed automatically by the image matching + psArray *regions = NULL, *kernels = NULL; // Regions and kernels used in subtraction + if (!ppStackMatch(readout, ®ions, &kernels, sources, targetPSF, config)) { + psError(PS_ERR_UNKNOWN, false, "Unable to match image %d --- ignoring.", i); + psFree(sources); + psFree(targetPSF); + return false; + } + + subRegions->data[i] = regions; + subKernels->data[i] = kernels; + + // Write the temporary convolved files + pmHDU *hdu = readout->parent->parent->parent->hdu; // HDU for convolved image + assert(hdu); + writeImage(imageNames->data[i], hdu->header, readout->image); + writeImage(maskNames->data[i], hdu->header, readout->mask); + writeImage(weightNames->data[i], hdu->header, readout->weight); +#endif + + cells->data[i] = psMemIncrRefCounter(readout->parent); + filesIterateUp(config); + } + psFree(sources); + psFree(targetPSF); + + // Stack the convolved files + { + pmFPAfileActivate(config->files, false, NULL); + fileActivation(config, combineFiles, true); + if (psMetadataLookupBool(&mdok, config->arguments, "PHOTOMETRY")) { + fileActivation(config, photFiles, true); + } + pmFPAview *view = filesIterateDown(config); + if (!view) { + psFree(cells); + psFree(subKernels); + psFree(subRegions); + 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++) { + readouts->data[i] = pmReadoutAlloc(cells->data[i]); // 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++) { + imageFits->data[i] = psFitsOpen(imageNames->data[i], "r"); + maskFits->data[i] = psFitsOpen(maskNames->data[i], "r"); + weightFits->data[i] = psFitsOpen(weightNames->data[i], "r"); + } + + // Read convolutions by chunks + bool more = true; // More to read? + 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, overlap) || + !pmReadoutReadChunkMask(readout, maskFits->data[i], 0, numScans, overlap) || + !pmReadoutReadChunkWeight(readout, weightFits->data[i], 0, numScans, overlap)) { + 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; } - - // Perform the analysis - if (!ppStackReadout(config, view)) { - psError(PS_ERR_UNKNOWN, false, "Unable to stack images.\n"); - return false; - } - - if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) { - return false; - } - } - - // Perform statistics on the cell - if (stats) { - ppStatsFPA(stats, output->fpa, view, maskBlank, config); - } - - if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) { + } + +#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 (!ppStackReadout(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; } - } - - if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) { - return false; - } - } - - if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) { - return false; - } - - psFree(view); + + 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); + psFree(subKernels); + psFree(subRegions); + for (int i = 0; i < num; i++) { + psFitsClose(imageFits->data[i]); + psFitsClose(maskFits->data[i]); + psFitsClose(weightFits->data[i]); + } + + if (psMetadataLookupBool(&mdok, config->arguments, "PHOTOMETRY") && + !ppStackPhotometry(config, outRO, view)) { + psError(PS_ERR_UNKNOWN, false, "Unable to perform photometry on output."); + psFree(outRO); + psFree(view); + return false; + } + + // Statistics on output + if (stats) { + ppStatsFPA(stats, outCell->parent->parent, view, maskBlank, config); + } + + // 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 + fileActivation(config, combineFiles, true); + filesIterateUp(config); + + psFree(outRO); + psFree(view); + } + // Write out summary statistics Index: /trunk/ppStack/src/ppStackMatch.c =================================================================== --- /trunk/ppStack/src/ppStackMatch.c (revision 16604) +++ /trunk/ppStack/src/ppStackMatch.c (revision 16605) @@ -9,9 +9,19 @@ #include "ppStack.h" +#define ARRAY_BUFFER 16 // Number to add to array at a time + + //#define TESTING -bool ppStackMatch(pmReadout *output, const pmReadout *input, const pmReadout *sourcesRO, - const pmPSF *psf, const pmConfig *config) +bool ppStackMatch(pmReadout *readout, psArray **regions, psArray **kernels, + const pmReadout *sourcesRO, const pmPSF *psf, const pmConfig *config) { + assert(readout); + assert(regions && !*regions); + assert(kernels && !*kernels); + assert(sourcesRO); + assert(psf); + assert(config); + // Look up appropriate values from the ppSub recipe bool mdok; // Status of MD lookup @@ -68,5 +78,5 @@ pmReadout *fake = pmReadoutAlloc(NULL); // Fake readout with target PSF - if (!pmReadoutFakeFromSources(fake, input->image->numCols, input->image->numRows, sources, NULL, NULL, + if (!pmReadoutFakeFromSources(fake, readout->image->numCols, readout->image->numRows, sources, NULL, NULL, psf, powf(10.0, -0.4 * maxMag), 0, false)) { psError(PS_ERR_UNKNOWN, false, "Unable to generate fake image with target PSF."); @@ -85,5 +95,6 @@ // Do the image matching - if (!pmSubtractionMatch(output, NULL, input, fake, footprint, regionSize, spacing, threshold, + pmReadout *output = pmReadoutAlloc(NULL); // Output readout, for holding results temporarily + if (!pmSubtractionMatch(output, NULL, readout, fake, footprint, regionSize, spacing, threshold, sources, stampsName, type, size, order, widths, orders, inner, ringsOrder, binning, optimum, optWidths, optOrder, optThresh, iter, rej, maskBad, @@ -92,4 +103,5 @@ psFree(fake); psFree(optWidths); + psFree(output); return false; } @@ -97,4 +109,51 @@ psFree(optWidths); + // Replace original images with convolved + psFree(readout->image); + psFree(readout->mask); + psFree(readout->weight); + readout->image = psMemIncrRefCounter(output->image); + readout->mask = psMemIncrRefCounter(output->mask); + readout->weight = psMemIncrRefCounter(output->weight); + + // Extract the regions and solutions used in the image matching + // This stops them from being freed when we iterate back up the FPA + *regions = psArrayAllocEmpty(ARRAY_BUFFER); // Array of regions + { + psString regex = NULL; // Regular expression + psStringAppend(®ex, "^%s$", PM_SUBTRACTION_ANALYSIS_REGION); + psMetadataIterator *iter = psMetadataIteratorAlloc(output->analysis, PS_LIST_HEAD, regex); // Iterator + psFree(regex); + 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); + } + *kernels = psArrayAllocEmpty(ARRAY_BUFFER); // Array of kernels + { + psString regex = NULL; // Regular expression + psStringAppend(®ex, "^%s$", PM_SUBTRACTION_ANALYSIS_KERNEL); + psMetadataIterator *iter = psMetadataIteratorAlloc(output->analysis, PS_LIST_HEAD, regex); // Iterator + psFree(regex); + psMetadataItem *item = NULL;// Item from iteration + while ((item = psMetadataGetAndIncrement(iter))) { + assert(item->type == PS_DATA_UNKNOWN); + // Set the normalisation dimensions, since these will be otherwise unavailable when reading the + // images by scans. + pmSubtractionKernels *kernel = item->data.V; // Kernel used in subtraction + kernel->numCols = readout->image->numCols; + kernel->numRows = readout->image->numRows; + + *kernels = psArrayAdd(*kernels, ARRAY_BUFFER, kernel); + } + psFree(iter); + } + assert((*regions)->n == (*kernels)->n); + + + psFree(output); + return true; } Index: /trunk/ppStack/src/ppStackPSF.c =================================================================== --- /trunk/ppStack/src/ppStackPSF.c (revision 16605) +++ /trunk/ppStack/src/ppStackPSF.c (revision 16605) @@ -0,0 +1,30 @@ +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include +#include + +#include "ppStack.h" + +pmPSF *ppStackPSF(const pmConfig *config, int numCols, int numRows, const psList *list) +{ + // Get the recipe values + int psfInstances = psMetadataLookupS32(NULL, 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(NULL, config->arguments, "PSF.ORDER"); // Spatial order for PSF + + // Solve for the target PSF + psArray *array = psListToArray(list); // Array of PSFs + pmPSF *psf = pmPSFEnvelope(numCols, numRows, array, psfInstances, psfRadius, psfModel, + psfOrder, psfOrder); + psFree(array); + if (!psf) { + psError(PS_ERR_UNKNOWN, false, "Unable to determine output PSF."); + return NULL; + } + + return psf; +} Index: /trunk/ppStack/src/ppStackPhotometry.c =================================================================== --- /trunk/ppStack/src/ppStackPhotometry.c (revision 16605) +++ /trunk/ppStack/src/ppStackPhotometry.c (revision 16605) @@ -0,0 +1,27 @@ +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include +#include +#include + +#include "ppStack.h" + +bool ppStackPhotometry(pmConfig *config, const pmReadout *readout, const pmFPAview *view) +{ + pmFPAfile *photFile = psMetadataLookupPtr(NULL, config->files, "PSPHOT.INPUT"); + pmFPACopy(photFile->fpa, readout->parent->parent->parent); + + if (!psphotReadout(config, view)) { + // Clear the error, so that the output files are written. + psWarning("Unable to perform photometry on stacked image."); + psErrorStackPrint(stderr, "Error stack from photometry:"); + psErrorClear(); + } + + pmFPAfileActivate(config->files, false, "PSPHOT.INPUT"); + + return true; +} Index: /trunk/ppStack/src/ppStackReadout.c =================================================================== --- /trunk/ppStack/src/ppStackReadout.c (revision 16604) +++ /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;