Index: /branches/pap_branch_080207/ppStack/src/Makefile.am =================================================================== --- /branches/pap_branch_080207/ppStack/src/Makefile.am (revision 16381) +++ /branches/pap_branch_080207/ppStack/src/Makefile.am (revision 16382) @@ -9,5 +9,7 @@ ppStackCamera.c \ ppStackLoop.c \ + ppStackPSF.c \ ppStackReadout.c \ + ppStackPhotometry.c \ ppStackVersion.c \ ppStackMatch.c Index: /branches/pap_branch_080207/ppStack/src/ppStack.h =================================================================== --- /branches/pap_branch_080207/ppStack/src/ppStack.h (revision 16381) +++ /branches/pap_branch_080207/ppStack/src/ppStack.h (revision 16382) @@ -19,7 +19,20 @@ ); +// 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 + ); + +// Perform photometry on stack +bool ppStackPhotometry(pmConfig *config, // Configuration + const pmReadout *readout, // Readout to be photometered + const pmFPAview *view // View to readout ); @@ -35,6 +48,5 @@ /// 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 const pmReadout *sourcesRO, ///< Readout with sources const pmPSF *psf, ///< Target PSF Index: /branches/pap_branch_080207/ppStack/src/ppStackLoop.c =================================================================== --- /branches/pap_branch_080207/ppStack/src/ppStackLoop.c (revision 16381) +++ /branches/pap_branch_080207/ppStack/src/ppStackLoop.c (revision 16382) @@ -13,4 +13,5 @@ // 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", @@ -23,22 +24,18 @@ "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 for convolution; these should be turned on one by one -static char *convolutionFiles[] = { "PPSTACK.INPUT", "PPSTACK.INPUT.MASK", "PPSTACK.INPUT.WEIGHT", - "PPSTACK.OUTCONV", "PPSTACK.OUTCONV.MASK", "PPSTACK.OUTCONV.WEIGHT", - 0 }; - // Input files for the combination -static char *inCombineFiles = { "PPSTACK.INCONV", "PPSTACK.INCONV.MASK", "PPSTACK.INCONV.WEIGHT", 0 }; +static char *inCombineFiles[] = { "PPSTACK.INCONV", "PPSTACK.INCONV.MASK", "PPSTACK.INCONV.WEIGHT", 0 }; // Output files for the combination and photometry -static char *outCombineFiles = { "PPSTACK.OUTPUT", "PPSTACK.OUTPUT.MASK", "PPSTACK.OUTPUT.WEIGHT", - "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 }; +static char *outCombineFiles[] = { "PPSTACK.OUTPUT", "PPSTACK.OUTPUT.MASK", "PPSTACK.OUTPUT.WEIGHT", + "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 }; @@ -103,5 +100,4 @@ - // Iterate down the hierarchy, loading files; we can get away with this because we're working on skycells static pmFPAview *filesIterateDown(pmConfig *config // Configuration @@ -158,4 +154,6 @@ bool ppStackLoop(pmConfig *config) { + assert(config); + psMaskType maskBlank = psMetadataLookupU8(NULL, config->arguments, "MASK.BLANK"); // Mask for blank reg. @@ -182,7 +180,17 @@ 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 + + + + // XXX Will need to adjust this when we're doing the convolution + int overlap = 0; // Overlap between consecutive scans + + + // Preparation iteration: Load the sources, and get a target PSF model - psArray *sources = NULL; // Sources to use for PSF matching + pmReadout *sources = NULL; // Readout with sources to use for PSF matching pmPSF *targetPSF = NULL; // Target PSF { @@ -194,7 +202,6 @@ } - pmReadout *sourcesRO = pmFPAfileThisReadout(config->files, view, "PPSTACK.SOURCES"); // Readout // We want to hang on to the 'sources' even when its host FPA is blown away - sources = psMemIncrRefCounter(psMetadataLookupPtr(&mdok, sourcesRO->analysis, "PSPHOT.SOURCES")); + sources = psMemIncrRefCounter(pmFPAfileThisReadout(config->files, view, "PPSTACK.SOURCES")); if (!sources) { psError(PS_ERR_UNKNOWN, true, "Unable to find sources."); @@ -207,27 +214,44 @@ "^PPSTACK.INPUT$"); // Iterator psMetadataItem *fileItem; // Item from iteration - int fileNum = 0; // Number of file 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, ro->parent->parent->analysis, "PSPHOT.PSF"); // 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; } psListAdd(psfList, PS_LIST_TAIL, psf); + + pmCell *cell = pmFPAviewThisCell(view, inputFile->fpa); // Cell of interest + pmHDU *hdu = pmHDUFromCell(cell); + 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); - psArray *psfArray = psListToArray(psfList); // Array of PSFs + targetPSF = ppStackPSF(config, numCols, numRows, psfList); psFree(psfList); - targetPSF = pmPSFEnvelope(numCols, numRows, psfArray, psfInstances, psfRadius, psfModel, - psfOrder, psfOrder); - psFree(psfArray); if (!targetPSF) { psError(PS_ERR_UNKNOWN, false, "Unable to determine output PSF."); @@ -241,10 +265,12 @@ // Generate convolutions and write them to disk - psArray *cells = psArrayAlloc(numFiles); // Cells for convolved images --- a handle for reading again - for (int i = 0; i < numFiles; i++) { + psArray *cells = psArrayAlloc(num); // Cells for convolved images --- a handle for reading again + for (int i = 0; i < num; i++) { pmFPAfileActivate(config->files, false, NULL); psArray *files = fileActivationSingle(config, prepareFiles, true, i); pmFPAview *view = filesIterateDown(config); if (!view) { + psFree(sources); + psFree(targetPSF); return false; } @@ -254,7 +280,8 @@ // Background subtraction, scaling and normalisation is performed automatically by the image matching - if (!ppStackMatch(ro, sources, outPSF, config)) { + if (!ppStackMatch(readout, sources, targetPSF, config)) { psError(PS_ERR_UNKNOWN, false, "Unable to match image %d --- ignoring.", i); - psFree(convolved); + psFree(sources); + psFree(targetPSF); return false; } @@ -290,4 +317,6 @@ filesIterateUp(config); } + psFree(sources); + psFree(targetPSF); // Set files to read at the readout level --- then when we iterate down, the pmFPAfileIO stuff should take @@ -295,6 +324,6 @@ fileSetDataLevel(config, inCombineFiles, PM_FPA_LEVEL_CHUNK); pmFPAfileActivate(config->files, false, NULL); - fileActivation(inCombineFiles, true); - fileActivation(outCombineFiles, true); + fileActivation(config, inCombineFiles, true); + fileActivation(config, outCombineFiles, true); { @@ -304,18 +333,17 @@ return false; } - pmReadout *outRO = pmFPAfileThisReadout(config->files, view, "PPSTACK.OUTPUT"); // Output readout - psFree(view); - - psArray *readouts = psArrayAlloc(numFiles); // Readouts for convolved images - for (int i = 0; i < numFiles; i++) { + 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); // Read convolutions by chunks - int numRead = 0; // Number of inputs read - int numChunk = 0; // Chunk number bool more = true; // More to read? for (int numChunk = 0; more; numChunk++) { - for (int i = 0; i < numFiles; i++) { + for (int i = 0; i < num; i++) { pmReadout *readout = readouts->data[i]; assert(readout); @@ -336,14 +364,23 @@ psError(PS_ERR_IO, false, "Unable to read chunk %d for file %d", numChunk, i); psFree(readouts); - psFree(cells); + psFree(outRO); + psFree(view); + return false; } } +#if 0 if (!ppStackReadout(config, outRO, readouts)) { psError(PS_ERR_UNKNOWN, false, "Unable to stack images.\n"); + psFree(readouts); + psFree(outRO); + psFree(view); return false; } - - for (int i = 0; i < numFiles && more; i++) { +#else + printf("Stack...\n"); +#endif + + for (int i = 0; i < num && more; i++) { pmReadout *readout = readouts->data[i]; assert(readout); @@ -364,15 +401,27 @@ } } + psFree(readouts); // Get rid of the input files - fileActivation(outCombineFiles, false); + fileActivation(config, outCombineFiles, false); filesIterateUp(config); - ppStackPhotometry(config, outRO); + 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, true, "Unable to find HDU for output."); + psError(PS_ERR_UNKNOWN, false, "Unable to find HDU for output."); return false; } @@ -382,15 +431,12 @@ ppStackVersionMetadata(hdu->header); - // Statistics on output - if (stats) { - ppStatsFPA(stats, output->fpa, view, maskBlank, config); - } - // Write out the output files - filesActivation(outCombineFiles, true); + fileActivation(config, outCombineFiles, true); filesIterateUp(config); - } - - psFree(view); + + psFree(outRO); + psFree(view); + } + // Write out summary statistics Index: /branches/pap_branch_080207/ppStack/src/ppStackPhotometry.c =================================================================== --- /branches/pap_branch_080207/ppStack/src/ppStackPhotometry.c (revision 16381) +++ /branches/pap_branch_080207/ppStack/src/ppStackPhotometry.c (revision 16382) @@ -10,8 +10,8 @@ #include "ppStack.h" -bool ppStackPhotometry(const pmConfig *config, const pmReadout *readout, const pmFPAview *view) +bool ppStackPhotometry(pmConfig *config, const pmReadout *readout, const pmFPAview *view) { pmFPAfile *photFile = psMetadataLookupPtr(NULL, config->files, "PSPHOT.INPUT"); - pmFPACopy(photFile->fpa, outRO->parent->parent->parent); + pmFPACopy(photFile->fpa, readout->parent->parent->parent); if (!psphotReadout(config, view)) { Index: /branches/pap_branch_080207/ppStack/src/ppStackReadout.c =================================================================== --- /branches/pap_branch_080207/ppStack/src/ppStackReadout.c (revision 16381) +++ /branches/pap_branch_080207/ppStack/src/ppStackReadout.c (revision 16382) @@ -10,213 +10,66 @@ #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) +bool ppStackReadout(const pmConfig *config, pmReadout *outRO, const psArray *readouts) { + assert(config); + assert(outRO); + assert(readouts); + // 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 psMaskType maskBad = psMetadataLookupU8(NULL, config->arguments, "MASK.BAD"); // Value to mask 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 +// float threshold = psMetadataLookupF32(NULL, config->arguments, "THRESHOLD.MASK"); // Threshold for mask deconvolution + + 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 0 + // I don't think 'scale' is needed, since the convolution matches the scales + float scale = psMetadataLookupF32(&mdok, fpa->analysis, "PPSTACK.SCALE"); // Rel. scale if (!mdok || !isfinite(scale)) { - psWarning("No SCALE supplied for image %d --- set to unity.", fileNum); + psWarning("No SCALE supplied for image %d --- set to unity.", i); scale = 1.0; } +#endif 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 (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 +77,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 +84,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."); @@ -239,35 +92,14 @@ } - // Don't need the original any more! - psFree(ro); - // 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 +109,4 @@ psFree(cellList); psFree(stack); - psFree(outRO); return false; } @@ -304,7 +135,15 @@ #endif - for (int i = 0; i < stack->n; i++) { + + + // XXX THIS NEEDS WORK +#if 0 + for (int i = 0; i < num; i++) { pmStackData *data = stack->data[i]; // Data for this image pmReadout *readout = data->readout; // Readout for this image + + + // XXX Need some other mechanism to get the subtraction kernel + // Extract the regions and solutions used in the image matching @@ -344,4 +183,6 @@ psFree(regions); } +#endif + #ifdef REJECTION_FILES @@ -370,14 +211,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 +229,5 @@ 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 return true;