Index: trunk/ppStack/src/ppStackReadout.c =================================================================== --- trunk/ppStack/src/ppStackReadout.c (revision 15457) +++ trunk/ppStack/src/ppStackReadout.c (revision 15844) @@ -28,4 +28,8 @@ 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.MODE"); // Model for PSF + int psfOrder = psMetadataLookupS32(&mdok, config->arguments, "PSF.ORDER"); // Spatial order for PSF // Get the output target @@ -43,10 +47,51 @@ "^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 fileNum = 0; // Number of file 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); @@ -85,90 +130,89 @@ // Generate convolved version of input - pmReadout *convolved = pmReadoutAlloc(NULL); // Convolved version of input readout -#ifndef NO_CONVOLUTION - if (!ppStackMatch(convolved, ro, sources, 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; - } + 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); - } + 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 - convolved = psMemIncrRefCounter(ro); - sources = NULL; -#endif // NO_CONVOLVUTION - - -#if 0 - // Background subtraction, scaling and normalisation is performed automatically by the image matching - - // Brain-dead background subtraction - if (!psImageBackground(stats, 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)); - } -#endif + } 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) { @@ -330,11 +374,12 @@ } -#if 0 - // 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)); - } -#endif + 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)); + } + } + psMetadataAddF32(outCell->concepts, PS_LIST_TAIL, "CELL.EXPOSURE", PS_META_REPLACE, "Summed exposure time (sec)", totExposure);