Index: trunk/ppStack/src/ppStack.h =================================================================== --- trunk/ppStack/src/ppStack.h (revision 15816) +++ trunk/ppStack/src/ppStack.h (revision 15844) @@ -38,4 +38,5 @@ const pmReadout *input, ///< Readout to be convolved const pmReadout *sourcesRO, ///< Readout with sources + const pmPSF *psf, ///< Target PSF const pmConfig *config ///< Configuration ); Index: trunk/ppStack/src/ppStackArguments.c =================================================================== --- trunk/ppStack/src/ppStackArguments.c (revision 15816) +++ trunk/ppStack/src/ppStackArguments.c (revision 15844) @@ -25,4 +25,5 @@ "\tMASK(STR): Mask filename\n" "\tWEIGHT(STR) Weight map filename\n" + "\tPSF(STR) PSF filename\n" "\tWEIGHTING(F32): Relative weighting to be applied\n" "\tSCALE(F32): Relative scaling to be applied\n", @@ -117,4 +118,8 @@ psMetadataAddF32(arguments, PS_LIST_TAIL, "-threshold-mask", 0, "Threshold for mask deconvolution", NAN); 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); + psMetadataAddF32(arguments, PS_LIST_TAIL, "-psf-radius", 0, "Radius for PSF generation", 20.0); + psMetadataAddStr(arguments, PS_LIST_TAIL, "-psf-model", 0, "Model name for PSF generation", "PS_MODEL_RGAUSS"); + psMetadataAddS32(arguments, PS_LIST_TAIL, "-psf-order", 0, "Spatial order for PSF generation", 3); if (argc == 1 || !psArgumentParse(arguments, &argc, argv) || argc != 3) { @@ -165,4 +170,9 @@ VALUE_ARG_RECIPE_FLOAT("-threshold-mask", "THRESHOLD.MASK", F32); + VALUE_ARG_RECIPE_INT("-psf-instances", "PSF.INSTANCES", S32, 0); + VALUE_ARG_RECIPE_FLOAT("-psf-radius", "PSF.RADIUS", F32); + VALUE_ARG_RECIPE_INT("-psf-order", "PSF.ORDER", S32, 0); + valueArgStr(config, arguments, "-psf-model", "PSF.MODEL", config->arguments); + if (psMetadataLookupBool(NULL, arguments, "-photometry") || psMetadataLookupBool(NULL, recipe, "PHOTOMETRY")) { Index: trunk/ppStack/src/ppStackCamera.c =================================================================== --- trunk/ppStack/src/ppStackCamera.c (revision 15816) +++ trunk/ppStack/src/ppStackCamera.c (revision 15844) @@ -14,4 +14,5 @@ { bool haveWeights = false; // Do we have weight maps? + bool havePSFs = false; // Do we have PSFs? psMetadata *inputs = psMetadataLookupMetadata(NULL, config->arguments, "INPUTS"); // The inputs info @@ -39,4 +40,5 @@ psString mask = psMetadataLookupStr(&mdok, input, "MASK"); // Name of mask psString weight = psMetadataLookupStr(&mdok, input, "WEIGHT"); // Name of weight map + psString psf = psMetadataLookupStr(&mdok, input, "PSF"); // Name of PSF float weighting = psMetadataLookupF32(&mdok, input, "WEIGHTING"); // Relative weighting @@ -115,4 +117,26 @@ } + // Optionally add the psf file + if (psf && strlen(psf) > 0) { + havePSFs = true; + psArray *psfFiles = psArrayAlloc(1); // Array of filenames for this FPA + psfFiles->data[0] = psMemIncrRefCounter(psf); + psMetadataAddArray(config->arguments, PS_LIST_TAIL, "PSF.FILENAMES", PS_META_REPLACE, + "Filenames of PSF files", psfFiles); + psFree(psfFiles); + + bool status; + pmFPAfile *psfFile = pmFPAfileBindFromArgs(&status, imageFile, config, "PSPHOT.PSF.LOAD", + "PSF.FILENAMES"); + if (!status) { + psError(PS_ERR_UNKNOWN, false, "Unable to define file from psf %d (%s)", i, psf); + return false; + } + if (psfFile->type != PM_FPA_FILE_PSF) { + psError(PS_ERR_IO, true, "PSPHOT.PSF.LOAD is not of type PSF"); + return false; + } + } + psMetadataAddF32(imageFile->fpa->analysis, PS_LIST_TAIL, "PPSTACK.WEIGHTING", 0, "Relative weighting for image", weighting); @@ -129,4 +153,7 @@ if (psMetadataLookup(config->arguments, "WEIGHT.FILENAMES")) { psMetadataRemoveKey(config->arguments, "WEIGHT.FILENAMES"); + } + if (psMetadataLookup(config->arguments, "PSF.FILENAMES")) { + psMetadataRemoveKey(config->arguments, "PSF.FILENAMES"); } @@ -214,4 +241,18 @@ } + // Output PSF + if (havePSFs) { + pmFPAfile *outPSF = pmFPAfileDefineOutput(config, output->fpa, "PSPHOT.PSF.SAVE"); + if (!outPSF) { + psError(PS_ERR_IO, false, _("Unable to generate output file from PSPHOT.PSF.SAVE")); + return false; + } + if (outPSF->type != PM_FPA_FILE_PSF) { + psError(PS_ERR_IO, true, "PSPHOT.PSF.SAVE is not of type PSF"); + return false; + } + outPSF->save = true; + } + // Sources for use as stamps bool status = false; // Found the file? Index: trunk/ppStack/src/ppStackMatch.c =================================================================== --- trunk/ppStack/src/ppStackMatch.c (revision 15816) +++ trunk/ppStack/src/ppStackMatch.c (revision 15844) @@ -12,5 +12,5 @@ bool ppStackMatch(pmReadout *output, const pmReadout *input, const pmReadout *sourcesRO, - const pmConfig *config) + const pmPSF *psf, const pmConfig *config) { // Look up appropriate values from the ppSub recipe @@ -46,5 +46,4 @@ // These values are specified specifically for stacking const char *stampsName = psMetadataLookupStr(NULL, config->arguments, "STAMPS"); // Filename for stamps - float target = psMetadataLookupF32(NULL, config->arguments, "TARGET"); // Target PSF width psVector *optWidths = NULL; // Vector with FWHMs for optimum search @@ -66,6 +65,14 @@ } } - pmReadout *fake = pmReadoutFakeFromSources(input->image->numCols, input->image->numRows, - sources, target, powf(10.0, -0.4 * maxMag)); + + pmReadout *fake = pmReadoutAlloc(NULL); // Fake readout with target PSF + + if (!pmReadoutFakeFromSources(fake, input->image->numCols, input->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."); + psFree(fake); + psFree(optWidths); + return false; + } #ifdef TESTING Index: trunk/ppStack/src/ppStackReadout.c =================================================================== --- trunk/ppStack/src/ppStackReadout.c (revision 15816) +++ 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);