Index: /branches/eam_branches/psphot.20100506/doc/stack.txt =================================================================== --- /branches/eam_branches/psphot.20100506/doc/stack.txt (revision 27875) +++ /branches/eam_branches/psphot.20100506/doc/stack.txt (revision 27876) @@ -1,2 +1,21 @@ + +20100506: + + we may load RAW and/or CNV images. + + * options->convolve: + * if only one of RAW or CNV exists, convolve it to match target PSF + * if both exist, convolve desired one + + * if matching PSF exists, use it (unless told to re-measure) + + * output goes to OUT (which is then used by psphot analysis routines) + + * detect + + * if (RAW) ? RAW : OUT + + * + 20100503: Index: /branches/eam_branches/psphot.20100506/src/Makefile.am =================================================================== --- /branches/eam_branches/psphot.20100506/src/Makefile.am (revision 27875) +++ /branches/eam_branches/psphot.20100506/src/Makefile.am (revision 27876) @@ -95,4 +95,6 @@ psphotFitSourcesLinearStack.c \ psphotSourceMatch.c \ + psphotStackMatchPSFs.c \ + psphotStackMatchPSFsUtils.c \ psphotCleanup.c Index: /branches/eam_branches/psphot.20100506/src/psphot.h =================================================================== --- /branches/eam_branches/psphot.20100506/src/psphot.h (revision 27875) +++ /branches/eam_branches/psphot.20100506/src/psphot.h (revision 27876) @@ -367,3 +367,73 @@ int pmPhotObjSortByX (const void **a, const void **b); +/// Options for stacking process +typedef struct { + // Setup + + int num; // Number of inputs + bool convolve; // Convolve images? + // psArray *convImages, *convMasks, *convVariances; // Filenames for the temporary convolved images + + // bool matchZPs; // Adjust relative fluxes based on transparency analysis? + // bool photometry; // Perform photometry? + // psMetadata *stats; // Statistics for output + // FILE *statsFile; // File to which to write statistics + // psArray *origImages, *origMasks, *origVariances; // Filenames of the original images + // psArray *origCovars; // Original covariances matrices + // int quality; // Bad data quality flag + + // Prepare + pmPSF *psf; // Target PSF + psVector *inputSeeing; // Input seeing FWHMs + psVector *inputMask; // Mask for inputs + + float targetSeeing; // Target seeing FWHM + psArray *sourceLists; // Individual lists of sources for matching + psVector *norm; // Normalisation for each image + psArray *psfs; + + // psVector *exposures; // Exposure times + // float sumExposure; // Sum of exposure times + // float zp; // Zero point for output + // psVector *inputMask; // Mask for inputs + // psArray *sources; // Matched sources + + // Convolve + psArray *kernels; // PSF-matching kernels --- required in the stacking + psArray *regions; // PSF-matching regions --- required in the stacking + psVector *matchChi2; // chi^2 for stamps from matching + psVector *weightings; // Combination weightings for images (1/noise^2) + // psArray *cells; // Cells for convolved images --- a handle for reading again + // int numCols, numRows; // Size of image + // psArray *convCovars; // Convolved covariance matrices + + // Combine initial + // pmReadout *outRO; // Output readout + // pmReadout *expRO; // Exposure readout + // psArray *inspect; // Array of arrays of pixels to inspect + + // Rejection + // psArray *rejected; // Rejected pixels +} psphotStackOptions; + +/*** psphotStackMatchPSF prototypes ***/ +bool psphotStackMatchPSFs (pmConfig *config, const pmFPAview *view); +bool psphotStackMatchPSFsReadout (pmConfig *config, const pmFPAview *view, psphotStackOptions *options, int index); + +// psphotStackMatchPSFsUtils +psVector *SetOptWidths (bool *optimum, psMetadata *recipe); +pmReadout *makeFakeReadout(pmConfig *config, pmReadout *raw, psArray *sources, pmPSF *psf, psImageMaskType maskVal, int fullSize); +bool rescaleData(pmReadout *readout, pmConfig *config, psphotStackOptions *options, int index); +bool renormKernel(pmReadout *readout, psphotStackOptions *options, int index); +bool saveMatchData (pmReadout *readout, psphotStackOptions *options, int index); +bool matchKernel(pmConfig *config, pmReadout *cnv, pmReadout *raw, psphotStackOptions *options, int index); +bool dumpImageDiff(pmReadout *readoutConv, pmReadout *readoutFake, pmReadout *readoutRef, int index, char *rootname); +bool dumpImage(pmReadout *readoutOut, pmReadout *readoutRef, int index, char *rootname); +bool loadKernel (pmConfig *config, pmReadout *readoutCnv, psphotStackOptions *options, int index); +bool stackRenormaliseReadout(const pmConfig *config, pmReadout *readout); +psImage *stackBackgroundModel(pmReadout *ro, const pmConfig *config); +psArray *stackSourcesFilter(psArray *sources, int exclusion); +void coordsFromSource(float *x, float *y, const pmSource *source); +bool readImage(psImage **target, const char *name, const pmConfig *config); + #endif Index: /branches/eam_branches/psphot.20100506/src/psphotErrorCodes.dat =================================================================== --- /branches/eam_branches/psphot.20100506/src/psphotErrorCodes.dat (revision 27875) +++ /branches/eam_branches/psphot.20100506/src/psphotErrorCodes.dat (revision 27876) @@ -11,4 +11,5 @@ APERTURE Problem with aperture photometry SKY Problem in sky determination +IO Problem in data I/O # these errors correspond to standard exit conditions ARGUMENTS Incorrect arguments Index: /branches/eam_branches/psphot.20100506/src/psphotSetMaskBits.c =================================================================== --- /branches/eam_branches/psphot.20100506/src/psphotSetMaskBits.c (revision 27875) +++ /branches/eam_branches/psphot.20100506/src/psphotSetMaskBits.c (revision 27876) @@ -37,2 +37,4 @@ return true; } + +// XXX should these be in config->analysis or somewhere else besides 'recipe'? Index: /branches/eam_branches/psphot.20100506/src/psphotStackArguments.c =================================================================== --- /branches/eam_branches/psphot.20100506/src/psphotStackArguments.c (revision 27875) +++ /branches/eam_branches/psphot.20100506/src/psphotStackArguments.c (revision 27876) @@ -11,4 +11,6 @@ if (psArgumentGet(argc, argv, "-help")) writeHelpInfo(stdout); if (psArgumentGet(argc, argv, "-h")) writeHelpInfo(stdout); + + if (psArgumentGet(argc, argv, "-template")) dumpTemplate(); // load config data from default locations @@ -84,7 +86,8 @@ "\n" "where INPUTS.mdc contains various METADATAs, each with:\n" - "\tIMAGE(STR): Image filename\n" - "\tMASK(STR): Mask filename\n" - "\tVARIANCE(STR): Variance map filename\n" + "\tIMAGE : Image filename\n" + "\tMASK : Mask filename\n" + "\tVARIANCE : Variance map filename\n" + "(use -template to generate a sample input.mdc file)\n" "OUTROOT is the 'root name' for output files\n" "\n" @@ -144,2 +147,30 @@ } +static void dumpTemplate(void) { + + fprintf (stdout, "# this line is required for multiple INPUT blocks to be accepted\n"); + fprintf (stdout, "INPUT MULTI\n\n"); + + fprintf (stdout, "# copy and repeat the following block as needed (one per input image set)\n"); + fprintf (stdout, "# either RAW (unconvolved) or CNV (convolved) input images are required\n"); + fprintf (stdout, "# if both are supplied, by default RAW is used for detection, CNV is convolved (further) to match target PSF\n"); + fprintf (stdout, "# if MASK or VARIANCE images are not supplied, they will be generated\n"); + fprintf (stdout, "# if MASK or VARIANCE images are not supplied, they will be generated\n"); + fprintf (stdout, "# PSF may be supplied for the convolution target\n"); + fprintf (stdout, "INPUT METADATA\n"); + fprintf (stdout, " RAW:IMAGE STR file.im.fits # signal image filename\n"); + fprintf (stdout, " RAW:MASK STR file.mk.fits # mask image filename\n"); + fprintf (stdout, " RAW:VARIANCE STR file.wt.fits # variance image filename\n"); + fprintf (stdout, " RAW:PSF STR file.psf.fits # psf from input unconvolved image\n"); + + fprintf (stdout, " CNV:IMAGE STR file.im.fits # signal image filename\n"); + fprintf (stdout, " CNV:MASK STR file.mk.fits # mask image filename\n"); + fprintf (stdout, " CNV:VARIANCE STR file.wt.fits # variance image filename\n"); + fprintf (stdout, " CNV:PSF STR file.psf.fits # psf from input convolved image\n"); + + fprintf (stdout, " SOURCES STR file.cmf # measured source positions\n"); + fprintf (stdout, "END\n"); + + psLibFinalize(); + exit(PS_EXIT_SUCCESS); +} Index: /branches/eam_branches/psphot.20100506/src/psphotStackImageLoop.c =================================================================== --- /branches/eam_branches/psphot.20100506/src/psphotStackImageLoop.c (revision 27875) +++ /branches/eam_branches/psphot.20100506/src/psphotStackImageLoop.c (revision 27876) @@ -14,6 +14,10 @@ pmReadout *readout; - pmFPAfile *input = psMetadataLookupPtr (&status, config->files, "PSPHOT.INPUT"); - if (!status) { + pmFPAfile *inputRaw = psMetadataLookupPtr (&status, config->files, "PSPHOT.STACK.INPUT.RAW"); + pmFPAfile *inputCnv = psMetadataLookupPtr (&status, config->files, "PSPHOT.STACK.INPUT.CNV"); + + pmFPAfile *input = inputRaw ? inputRaw : inputCnv; + + if (!input) { psError(PSPHOT_ERR_PROG, false, "Can't find input data!"); return false; @@ -29,10 +33,10 @@ psLogMsg ("psphot", 4, "Chip %d: %x %x\n", view->chip, chip->file_exists, chip->process); if (! chip->process || ! chip->file_exists) { continue; } - if (!pmFPAfileIOChecks (config, view, PM_FPA_BEFORE)) ESCAPE ("failed input for Chip in psphotStack."); + // if (!pmFPAfileIOChecks (config, view, PM_FPA_BEFORE)) ESCAPE ("failed input for Chip in psphotStack."); // there is now only a single chip (multiple readouts?). loop over it and process while ((cell = pmFPAviewNextCell (view, input->fpa, 1)) != NULL) { psLogMsg ("psphot", 5, "Cell %d: %x %x\n", view->cell, cell->file_exists, cell->process); - if (!pmFPAfileIOChecks (config, view, PM_FPA_BEFORE)) ESCAPE ("failed input for Cell in psphotStack."); + // if (!pmFPAfileIOChecks (config, view, PM_FPA_BEFORE)) ESCAPE ("failed input for Cell in psphotStack."); // process each of the readouts @@ -41,6 +45,5 @@ if (! readout->data_exists) { continue; } -# if (0) - // uncomment to generate matched psfs + // PSF matching if (!psphotStackMatchPSFs (config, view)) { psError(psErrorCodeLast(), false, "failure in psphotStackMatchPSFs for chip %d, cell %d, readout %d\n", view->chip, view->cell, view->readout); @@ -48,5 +51,4 @@ return false; } -# endif // XXX for now, we assume there is only a single chip in the PHU: Index: /branches/eam_branches/psphot.20100506/src/psphotStackMatchPSFs.c =================================================================== --- /branches/eam_branches/psphot.20100506/src/psphotStackMatchPSFs.c (revision 27875) +++ /branches/eam_branches/psphot.20100506/src/psphotStackMatchPSFs.c (revision 27876) @@ -1,14 +1,57 @@ # include "psphotInternal.h" -bool psphotStackMatchPSFs (pmConfig *config, const pmFPAview *view, bool firstPass) +bool psphotStackMatchPSFs (pmConfig *config, const pmFPAview *view) { bool status = true; + + psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, "PSPHOT"); + psAssert(recipe, "We've thrown an error on this before."); int num = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.INPUT.NUM"); psAssert (status, "programming error: must define PSPHOT.INPUT.NUM"); - // loop over the available readouts + // 'options' carries info needed to perform the stack matching + psphotStackOptions *options = psphotStackOptionsAlloc(num); + + options->convolve = psMetadataLookupBool (&status, recipe, "PSPHOT.STACK.MATCH.PSF"); + psAssert (status, "PSPHOT.STACK.MATCH.PSF not in recipe"); + + if (options->convolve) { + char *convolveSource = psMetadataLookupStr (&status, recipe, "PSPHOT.STACK.MATCH.PSF.SOURCE"); + options->convolveSource = psphotStackConvolveSourceFromString (convolveSource); + if (options->convolveSource == PSPHOT_CNV_SRC_NONE) { + psError (PSPHOT_ERR_CONFIG, true, "stack convolution source not defined in recipe"); + return false; + } + } + + // loop over the available readouts (ignore chisq image) for (int i = 0; i < num; i++) { - if (!psphotStackMatchPSFsReadout (config, view, i)) { + if (!psphotStackMatchPSFsPrepare (config, view, options, i)) { + psError (PSPHOT_ERR_CONFIG, false, "failed to set PSF matching options for entry %d", i); + return false; + } + } + + // Generate target PSF + if (options->convolve) { + options->psf = ppStackPSF(config, options->numCols, options->numRows, options->psfs, options->inputMask); + if (!options->psf) { + psError(psErrorCodeLast(), false, "Unable to determine output PSF."); + return false; + } + psMetadataAddPtr(config->arguments, PS_LIST_TAIL, "PSF.TARGET", PS_DATA_UNKNOWN, "Target PSF for stack", options->psf); + options->targetSeeing = pmPSFtoFWHM(options->psf, 0.5 * options->numCols, 0.5 * options->numRows); // FWHM for target + psLogMsg("psphotStack", PS_LOG_INFO, "Target seeing FWHM: %f\n", options->targetSeeing); + + // XXX is this needed to supply the psf to psphot?? + // pmChip *outChip = pmFPAfileThisChip(config->files, view, "PPSTACK.TARGET.PSF"); // Output chip + // psMetadataAddPtr(outChip->analysis, PS_LIST_TAIL, "PSPHOT.PSF", PS_DATA_UNKNOWN, "Target PSF", options->psf); + // outChip->data_exists = true; + } + + // loop over the available readouts (ignore chisq image) + for (int i = 0; i < num; i++) { + if (!psphotStackMatchPSFsReadout (config, view, options, i)) { psError (PSPHOT_ERR_CONFIG, false, "failed to find initial detections for PSPHOT.INPUT entry %d", i); return false; @@ -19,16 +62,11 @@ // convolve the image to match desired PSF -bool psphotStackMatchPSFsReadout (pmConfig *config, const pmFPAview *view, int index) { - - bool status; - int pass; - float NSIGMA_PEAK = 25.0; - int NMAX = 0; +bool psphotStackMatchPSFsReadout (pmConfig *config, const pmFPAview *view, psphotStackOptions *options, int index) { // find the currently selected readout - pmFPAfile *fileRaw = pmFPAfileSelectSingle(config->files, "PSPHOT.INPUT", index); // File of interest + pmFPAfile *fileRaw = pmFPAfileSelectSingle(config->files, "PSPHOT.STACK.INPUT.RAW", index); // File of interest psAssert (fileRaw, "missing file?"); - pmFPAfile *fileCnv = pmFPAfileSelectSingle(config->files, "PSPHOT.INPUT.CONV", index); // File of interest + pmFPAfile *fileCnv = pmFPAfileSelectSingle(config->files, "PSPHOT.STACK.INPUT.CNV", index); // File of interest psAssert (fileCnv, "missing file?"); @@ -39,33 +77,7 @@ psAssert (readoutCnv, "missing readout?"); - // user-defined masks to test for good/bad pixels (build from recipe list if not yet set) - psImageMaskType maskVal = psMetadataLookupImageMask(&status, recipe, "MASK.PSPHOT"); // Mask value for bad pixels - psAssert (maskVal, "missing mask value?"); - - /***** set up recipe options *****/ - - psMetadata *stackRecipe = psMetadataLookupMetadata(NULL, config->recipes, PPSTACK_RECIPE); // ppStack recipe - psAssert(stackRecipe, "We've thrown an error on this before."); - - // Look up appropriate values from the ppSub recipe - psMetadata *subRecipe = psMetadataLookupMetadata(NULL, config->recipes, "PPSUB"); // PPSUB recipe - psAssert(subRecipe, "recipe missing"); - - int size = psMetadataLookupS32(NULL, subRecipe, "KERNEL.SIZE"); // Kernel half-size - - float deconvLimit = psMetadataLookupF32(NULL, stackRecipe, "DECONV.LIMIT"); // Limit on deconvolution fraction - - psString maskValStr = psMetadataLookupStr(NULL, subRecipe, "MASK.VAL"); // Name of bits to mask going in - psImageMaskType maskVal = pmConfigMaskGet(maskValStr, config); // Bits to mask going in to pmSubtractionMatch - psString maskPoorStr = psMetadataLookupStr(NULL, stackRecipe, "MASK.POOR"); // Name of bits to mask for poor - psImageMaskType maskPoor = pmConfigMaskGet(maskPoorStr, config); // Bits to mask for poor pixels - psString maskBadStr = psMetadataLookupStr(NULL, stackRecipe, "MASK.BAD"); // Name of bits to mask for bad - psImageMaskType maskBad = pmConfigMaskGet(maskBadStr, config); // Bits to mask for bad pixels - - bool mdok; // Status of MD lookup - float penalty = psMetadataLookupF32(NULL, subRecipe, "PENALTY"); // Penalty for wideness - int threads = psMetadataLookupS32(NULL, config->arguments, "-threads"); // Number of threads - - if (!pmReadoutMaskNonfinite(readout, maskVal)) { + // set NAN pixels to 'SAT' + psImageMaskType maskVal = pmConfigMaskGet("SAT", config); + if (!pmReadoutMaskNonfinite(readoutRaw, maskVal)) { psError(psErrorCodeLast(), false, "Unable to mask non-finite pixels in readout."); return false; @@ -74,79 +86,29 @@ // Image Matching (PSFs or just flux) if (options->convolve) { - // Full match of PSFs - pmReadout *conv = pmReadoutAlloc(NULL); // Conv readout, for holding results temporarily - - // Read previously produced kernel - if (psMetadataLookupBool(NULL, config->arguments, "PPSTACK.DEBUG.STACK")) { - loadKernel(); - } else { - matchKernel(); - } // !DEBUG.STACK - - saveMatchData(); - - saveChiSquare(); - - renormKernel(); - - // Reject image completely if the maximum deconvolution fraction exceeds the limit - float deconv = psMetadataLookupF32(NULL, conv->analysis, PM_SUBTRACTION_ANALYSIS_DECONV_MAX); // Max deconvolution fraction - if (deconv > deconvLimit) { - psWarning("Maximum deconvolution fraction (%f) exceeds limit (%f) --- rejecting image %d\n", deconv, deconvLimit, index); - psFree(conv); - return NULL; - } - - readout->analysis = psMetadataCopy(readout->analysis, conv->analysis); - - psFree(conv); + matchKernel(config, readoutCnv, readoutRaw, options, index); + saveMatchData(readoutCnv, options, index); + // renormKernel(readoutCnv, options, index); } else { - // only match the flux - float norm = powf(10.0, -0.4 * options->norm->data.F32[index]); // Normalisation - psBinaryOp(readout->image, readout->image, "*", psScalarAlloc(norm, PS_TYPE_F32)); - psBinaryOp(readout->variance, readout->variance, "*", psScalarAlloc(PS_SQR(norm), PS_TYPE_F32)); + // only match the flux (NO! not for multi-filter, at least!) + // XXX do not generate readoutCnv in this case? + // float norm = powf(10.0, -0.4 * options->norm->data.F32[index]); // Normalisation + // psBinaryOp(readoutRaw->image, readoutRaw->image, "*", psScalarAlloc(norm, PS_TYPE_F32)); + // psBinaryOp(readoutRaw->variance, readoutRaw->variance, "*", psScalarAlloc(PS_SQR(norm), PS_TYPE_F32)); } - // Ensure the background value is zero - psStats *bg = psStatsAlloc(PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV); // Statistics for background - psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS); // Random number generator - if (!psImageBackground(bg, NULL, readout->image, readout->mask, maskVal | maskBad, rng)) { - psWarning("Can't measure background for image."); - psErrorClear(); - } else { - if (!psMetadataLookupBool(NULL, config->arguments, "PPSTACK.SKIP.BG.SUB")) { - psLogMsg("ppStack", PS_LOG_INFO, "Correcting convolved image background by %lf (+/- %lf)", - psStatsGetValue(bg, PS_STAT_ROBUST_MEDIAN), psStatsGetValue(bg, PS_STAT_ROBUST_STDEV)); - (void)psBinaryOp(readout->image, readout->image, "-", - psScalarAlloc(psStatsGetValue(bg, PS_STAT_ROBUST_MEDIAN), PS_TYPE_F32)); - } - } + rescaleData(readoutCnv, config, options, index); - if (!stackRenormaliseReadout(config, readout)) { - psFree(rng); - psFree(bg); - return false; - } - - // Measure the variance level for the weighting - if (psMetadataLookupBool(NULL, stackRecipe, "WEIGHTS")) { - if (!psImageBackground(bg, NULL, readout->variance, readout->mask, maskVal | maskBad, rng)) { - psError(PPSTACK_ERR_DATA, false, "Can't measure mean variance for image."); - psFree(rng); - psFree(bg); - return false; - } - options->weightings->data.F32[index] = 1.0 / (psStatsGetValue(bg, PS_STAT_ROBUST_MEDIAN) * psImageCovarianceFactor(readout->covariance)); - } else { - options->weightings->data.F32[index] = 1.0; - } - psLogMsg("ppStack", PS_LOG_INFO, "Weighting for image %d is %f\n", - index, options->weightings->data.F32[index]); - - psFree(rng); - psFree(bg); - - dumpImage3(); + dumpImage(readoutCnv, readoutRaw, index, "convolved"); return true; } + + +# if (0) +// Read previously produced kernel +if (psMetadataLookupBool(NULL, config->arguments, "PPSTACK.DEBUG.STACK")) { + loadKernel(config, readoutCnv, options, index); +} else { + matchKernel(config, readoutCnv, readoutRaw, options, index); +} +# endif Index: /branches/eam_branches/psphot.20100506/src/psphotStackMatchPSFsPrepare.c =================================================================== --- /branches/eam_branches/psphot.20100506/src/psphotStackMatchPSFsPrepare.c (revision 27876) +++ /branches/eam_branches/psphot.20100506/src/psphotStackMatchPSFsPrepare.c (revision 27876) @@ -0,0 +1,119 @@ +# include "psphotInternal.h" + +// set up the stacking parameters +bool psphotStackMatchPSFsPrepare (pmConfig *config, const pmFPAview *view, psphotStackOptions *options, int index) { + + if (!options->convolve) { + psLogMsg ("psphotStack", PS_LOG_INFO, "psf-matching NOT selected"); + return true; + } + + // which image do we want to convolve? RAW, CNV, AUTO? + // find the currently selected readout + pmFPAfile *fileRaw = pmFPAfileSelectSingle(config->files, "PSPHOT.STACK.INPUT.RAW", index); // File of interest + pmFPAfile *fileCnv = pmFPAfileSelectSingle(config->files, "PSPHOT.STACK.INPUT.CNV", index); // File of interest + + pmFPAfile *fileSrc = NULL; + pmPSF *psf = NULL; + + switch (options->convolveSource) { + case AUTO: + fileSrc = fileCnv ? fileCnv : fileRaw; + break; + + case RAW: + fileSrc = fileRaw; + break; + + case CNV: + fileSrc = fileCnv; + break; + + default: + psAbort("impossible case"); + } + if (!fileSrc) { + psError(PSPHOT_ERR_CONFIG, "desired convolution source is missing (cnv : %llx, raw : %llx)", (long long) fileCnv, (long long) fileRaw); + return false; + } + + // select the corresponding input psf + { + pmChip *chip = pmFPAviewThisChip(view, fileSrc->fpa); // The chip holds the PSF + pmPSF *psf = psMetadataLookupPtr(NULL, chip->analysis, "PSPHOT.PSF"); // PSF + if (!psf) { + // XXX if we were not supplied a PSF, we should be able to generate one by calling psphot + psError(PPSTACK_ERR_PROG, false, "Unable to find PSF."); + return false; + } + options->psfs->data[index] = psMemIncrRefCounter(psf); + + // find the image size + pmCell *cell = pmFPAviewThisCell(view, fileSrc->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 + psAssert ((naxis1 > 0) && (naxis2 > 0), "Unable to determine size of image from PSF."); + if (!options->numCols) { + options->numCols = naxis1; + } else { + if (options->numCols != naxis1) { + psError (PSPHOT_ERR_CONFIG, true, "mismatched sizes (NAXIS1) for input images"); + return false; + } + } + if (!options->numRows) { + options->numRows = naxis2; + } else { + if (options->numRows != naxis2) { + psError (PSPHOT_ERR_CONFIG, true, "mismatched sizes (NAXIS2) for input images"); + return false; + } + } + } + + // load the sources (used to find reference sources for the kernel stamps) + { + pmFPAfile *outputImage = pmFPAfileDefineOutput(config, NULL, "PSPHOT.STACK.OUTPUT.IMAGE"); + pmReadout *readout = pmFPAviewThisReadout(view, outputImage->fpa); // Readout with sources + detections = psMetadataLookupPtr(NULL, readout->analysis, "PSPHOT.DETECTIONS"); // Sources + if (!detections || !detections->allSources) { + psWarning("No detections found for image %d --- rejecting.", i); + options->inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PPSTACK_MASK_CAL; + continue; + } + psAssert (detections->allSources, "missing sources?"); + options->sourceLists->data[i] = psMemIncrRefCounter(detections->allSources); + } + + // determine the input seeing + { + // XXX set this based on the mode (+1 for poly, +0 for map) + float xNum = PS_MAX(psf->trendNx, 1); + float yNum = PS_MAX(psf->trendNy, 1); // Number of realisations + float sumFWHM = 0.0; // FWHM for image + int numFWHM = 0; // Number of FWHM measurements + for (float y = 0; y < yNum; y += 1.0) { + float yPos = options->numRows * ((y + 0.5) / yNum); + for (float x = 0; x < xNum; x++) { + float xPos = options->numCols * ((x + 0.5) / xNum); + float fwhm = pmPSFtoFWHM(psf, xPos, yPos); // FWHM for image + if (isfinite(fwhm)) { + sumFWHM += fwhm; + numFWHM++; + } + } + } + if (numFWHM == 0) { + options->inputSeeing->data.F32[index] = NAN; + options->inputMask->data.PS_TYPE_VECTOR_MASK_DATA[index] = PPSTACK_MASK_PSF; + psLogMsg("ppStack", PS_LOG_INFO, "Unable to measure PSF FWHM for image %d --- rejected.", index); + } else { + options->inputSeeing->data.F32[index] = sumFWHM / (float)numFWHM; + } + psLogMsg ("psphotStack", PS_LOG_INFO, "Input Seeing for %d: %f\n", index, options->inputSeeing->data.F32[index]); + } + + return true; +} Index: /branches/eam_branches/psphot.20100506/src/psphotStackMatchPSFsUtils.c =================================================================== --- /branches/eam_branches/psphot.20100506/src/psphotStackMatchPSFsUtils.c (revision 27875) +++ /branches/eam_branches/psphot.20100506/src/psphotStackMatchPSFsUtils.c (revision 27876) @@ -1,10 +1,4 @@ -/***** defines *****/ - -#define ARRAY_BUFFER 16 // Number to add to array at a time -#define MAG_IGNORE 50 // Ignore magnitudes fainter than this --- they're not real! -#define FAKE_SIZE 1 // Size of fake convolution kernel -#define SOURCE_MASK (PM_SOURCE_MODE_FAIL | PM_SOURCE_MODE_DEFECT | PM_SOURCE_MODE_SATURATED | PM_SOURCE_MODE_CR_LIMIT | PM_SOURCE_MODE_EXT_LIMIT) // Mask to apply to input sources -#define NOISE_FRACTION 0.01 // Set minimum flux to this fraction of noise -#define COVAR_FRAC 0.01 // Truncation fraction for covariance matrix +# include "psphotInternal.h" +# define ARRAY_BUFFER 16 // Number to add to array at a time // XXX better name @@ -18,10 +12,10 @@ psFree(resolved); if (!fits) { - psError(PPSTACK_ERR_IO, false, "Unable to open previously produced image: %s", name); + psError(PSPHOT_ERR_IO, false, "Unable to open previously produced image: %s", name); return false; } psImage *image = psFitsReadImage(fits, psRegionSet(0,0,0,0), 0); // Image of interest if (!image) { - psError(PPSTACK_ERR_IO, false, "Unable to read previously produced image: %s", name); + psError(PSPHOT_ERR_IO, false, "Unable to read previously produced image: %s", name); psFitsClose(fits); return false; @@ -52,5 +46,5 @@ psArray *stackSourcesFilter(psArray *sources, // Source list to filter - int exclusion // Exclusion zone, pixels + int exclusion // Exclusion zone, pixels ) { @@ -90,5 +84,5 @@ long numWithin = psTreeWithin(tree, coords, exclusion); // Number within exclusion zone - psTrace("ppStack", 9, "Source at %.0lf,%.0lf has %ld sources in exclusion zone", + psTrace("psphotStack", 9, "Source at %.0lf,%.0lf has %ld sources in exclusion zone", coords->data.F64[0], coords->data.F64[1], numWithin); if (numWithin == 1) { @@ -104,5 +98,5 @@ psFree(y); - psLogMsg("ppStack", PS_LOG_INFO, "Filtered out %d of %d sources", numFiltered, numGood); + psLogMsg("psphotStack", PS_LOG_INFO, "Filtered out %d of %d sources", numFiltered, numGood); return filtered; @@ -111,5 +105,5 @@ // Add background into the fake image // Based on ppSubBackground() -static psImage *stackBackgroundModel(pmReadout *ro, // Readout for which to generate background model +psImage *stackBackgroundModel(pmReadout *ro, // Readout for which to generate background model const pmConfig *config // Configuration ) @@ -121,7 +115,7 @@ int numCols = image->numCols, numRows = image->numRows; // Size of image - psMetadata *ppStackRecipe = psMetadataLookupPtr(NULL, config->recipes, PPSTACK_RECIPE); + psMetadata *ppStackRecipe = psMetadataLookupPtr(NULL, config->recipes, "PPSTACK"); psAssert(ppStackRecipe, "Need PPSTACK recipe"); - psMetadata *psphotRecipe = psMetadataLookupPtr(NULL, config->recipes, PSPHOT_RECIPE); + psMetadata *psphotRecipe = psMetadataLookupPtr(NULL, config->recipes, "PSPHOT"); psAssert(psphotRecipe, "Need PSPHOT recipe"); @@ -138,5 +132,5 @@ psImage *unbinned = psImageAlloc(numCols, numRows, PS_TYPE_F32); // Unbinned background model if (!psImageUnbin(unbinned, binned, binning)) { - psError(PPSTACK_ERR_DATA, false, "Unable to unbin background model"); + psError(PSPHOT_ERR_DATA, false, "Unable to unbin background model"); psFree(binned); psFree(unbinned); @@ -153,8 +147,7 @@ ) { -#if 1 bool mdok; // Status of metadata lookups - psMetadata *recipe = psMetadataLookupPtr(NULL, config->recipes, PPSTACK_RECIPE); // Recipe for ppStack + psMetadata *recipe = psMetadataLookupPtr(NULL, config->recipes, "PPSTACK"); // Recipe for ppStack psAssert(recipe, "Need PPSTACK recipe"); @@ -163,15 +156,15 @@ int num = psMetadataLookupS32(&mdok, recipe, "RENORM.NUM"); if (!mdok) { - psError(PPSTACK_ERR_CONFIG, true, "RENORM.NUM is not set in the recipe"); + psError(PSPHOT_ERR_CONFIG, true, "RENORM.NUM is not set in the recipe"); return false; } float minValid = psMetadataLookupF32(&mdok, recipe, "RENORM.MIN"); if (!mdok) { - psError(PPSTACK_ERR_CONFIG, true, "RENORM.MIN is not set in the recipe"); + psError(PSPHOT_ERR_CONFIG, true, "RENORM.MIN is not set in the recipe"); return false; } float maxValid = psMetadataLookupF32(&mdok, recipe, "RENORM.MAX"); if (!mdok) { - psError(PPSTACK_ERR_CONFIG, true, "RENORM.MAX is not set in the recipe"); + psError(PSPHOT_ERR_CONFIG, true, "RENORM.MAX is not set in the recipe"); return false; } @@ -181,7 +174,4 @@ psImageCovarianceTransfer(readout->variance, readout->covariance); return pmReadoutVarianceRenormalise(readout, maskBad, num, minValid, maxValid); -#else - return true; -#endif } @@ -190,367 +180,422 @@ // It implicitly assumes the output root name is the same between invocations. -bool loadKernel () { - pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PPSTACK.CONV.KERNEL", index); - psAssert(file, "Require file"); - - pmFPAview *view = pmFPAviewAlloc(0); // View to readout of interest - view->chip = view->cell = view->readout = 0; - psString filename = pmFPAfileNameFromRule(file->filerule, file, view); // Filename of interest - - // Read convolution kernel - psString resolved = pmConfigConvertFilename(filename, config, false, false); // Resolved filename - psFree(filename); - psFits *fits = psFitsOpen(resolved, "r"); // FITS file for subtraction kernel - psFree(resolved); - if (!fits || !pmReadoutReadSubtractionKernels(conv, fits)) { - psError(PPSTACK_ERR_IO, false, "Unable to read previously produced kernel"); - psFitsClose(fits); - return false; - } - psFitsClose(fits); - - if (!readImage(&readout->image, options->convImages->data[index], config) || - !readImage(&readout->mask, options->convMasks->data[index], config) || - !readImage(&readout->variance, options->convVariances->data[index], config)) { - psError(PPSTACK_ERR_IO, false, "Unable to read previously produced image."); - return false; - } - - psRegion *region = psMetadataLookupPtr(NULL, conv->analysis, - PM_SUBTRACTION_ANALYSIS_REGION); // Convolution region - pmSubtractionKernels *kernels = psMetadataLookupPtr(NULL, conv->analysis, - PM_SUBTRACTION_ANALYSIS_KERNEL); - - pmSubtractionAnalysis(conv->analysis, NULL, kernels, region, - readout->image->numCols, readout->image->numRows); - - psKernel *kernel = pmSubtractionKernel(kernels, 0.0, 0.0, false); // Convolution kernel - bool oldThreads = psImageCovarianceSetThreads(true); // Old thread setting - psKernel *covar = psImageCovarianceCalculate(kernel, readout->covariance); // Covariance matrix - psImageCovarianceSetThreads(oldThreads); - psFree(readout->covariance); - readout->covariance = covar; - psFree(kernel); -} - -bool dumpImage() { - // XXX should be optional - { - pmHDU *hdu = pmHDUFromCell(readout->parent); - psString name = NULL; - psStringAppend(&name, "fake_%03d.fits", index); - pmStackVisualPlotTestImage(fake->image, name); - psFits *fits = psFitsOpen(name, "w"); - psFree(name); - psFitsWriteImage(fits, hdu->header, fake->image, 0, NULL); - psFitsClose(fits); - } - { - pmHDU *hdu = pmHDUFromCell(readout->parent); - psString name = NULL; - psStringAppend(&name, "real_%03d.fits", index); - pmStackVisualPlotTestImage(readout->image, name); - psFits *fits = psFitsOpen(name, "w"); - psFree(name); - psFitsWriteImage(fits, hdu->header, readout->image, 0, NULL); - psFitsClose(fits); - } -} - -bool dumpImage2() { - // XXX should be optional - - { - pmHDU *hdu = pmHDUFromCell(readout->parent); - psString name = NULL; - psStringAppend(&name, "conv_%03d.fits", index); - pmStackVisualPlotTestImage(conv->image, name); - psFits *fits = psFitsOpen(name, "w"); - psFree(name); - psFitsWriteImage(fits, hdu->header, conv->image, 0, NULL); - psFitsClose(fits); - } - { - pmHDU *hdu = pmHDUFromCell(readout->parent); - psString name = NULL; - psStringAppend(&name, "diff_%03d.fits", index); - pmStackVisualPlotTestImage(fake->image, name); - psFits *fits = psFitsOpen(name, "w"); - psFree(name); - psBinaryOp(fake->image, conv->image, "-", fake->image); - psFitsWriteImage(fits, hdu->header, fake->image, 0, NULL); - psFitsClose(fits); - } -} - -bool dumpImage3() +bool loadKernel (pmConfig *config, pmReadout *readoutCnv, psphotStackOptions *options, int index) { + + // Read the convolution kernel from the saved file + pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PPSTACK.CONV.KERNEL", index); + psAssert(file, "Require file"); + + pmFPAview *view = pmFPAviewAlloc(0); // View to readout of interest + view->chip = view->cell = view->readout = 0; + psString filename = pmFPAfileNameFromRule(file->filerule, file, view); // Filename of interest + + // Read convolution kernel data + psString resolved = pmConfigConvertFilename(filename, config, false, false); // Resolved filename + psFree(filename); + psFits *fits = psFitsOpen(resolved, "r"); // FITS file for subtraction kernel + psFree(resolved); + if (!fits || !pmReadoutReadSubtractionKernels(readoutCnv, fits)) { + psError(PSPHOT_ERR_IO, false, "Unable to read previously produced kernel"); + psFitsClose(fits); + return false; + } + psFitsClose(fits); + + // read the convolved pixels (image, mask, variance) -- names are pre-defined + if (!readImage(&readoutCnv->image, options->convImages->data[index], config) || + !readImage(&readoutCnv->mask, options->convMasks->data[index], config) || + !readImage(&readoutCnv->variance, options->convVariances->data[index], config)) { + psError(PSPHOT_ERR_IO, false, "Unable to read previously produced image."); + return false; + } + + // XXX ??? not sure what is happening here -- consult Paul + psRegion *region = psMetadataLookupPtr(NULL, readoutCnv->analysis, PM_SUBTRACTION_ANALYSIS_REGION); // Convolution region + pmSubtractionKernels *kernels = psMetadataLookupPtr(NULL, readoutCnv->analysis, PM_SUBTRACTION_ANALYSIS_KERNEL); + + pmSubtractionAnalysis(readoutCnv->analysis, NULL, kernels, region, readoutCnv->image->numCols, readoutCnv->image->numRows); + + psKernel *kernel = pmSubtractionKernel(kernels, 0.0, 0.0, false); // Convolution kernel + + // update the covariance matrix + // XXX why is this needed if we have correctly read the saved data? + bool oldThreads = psImageCovarianceSetThreads(true); // Old thread setting + psKernel *covar = psImageCovarianceCalculate(kernel, readoutCnv->covariance); // Covariance matrix + psImageCovarianceSetThreads(oldThreads); + psFree(readoutCnv->covariance); + readoutCnv->covariance = covar; + psFree(kernel); + return true; +} + +bool dumpImage(pmReadout *readoutOut, pmReadout *readoutRef, int index, char *rootname) { + + pmHDU *hdu = pmHDUFromCell(readoutRef->parent); + psString name = NULL; + psStringAppend(&name, "%s_%03d.fits", rootname, index); + pmStackVisualPlotTestImage(readoutOut->image, name); + psFits *fits = psFitsOpen(name, "w"); + psFree(name); + psFitsWriteImage(fits, hdu->header, readoutOut->image, 0, NULL); + psFitsClose(fits); + return true; +} + +bool dumpImageDiff(pmReadout *readoutConv, pmReadout *readoutFake, pmReadout *readoutRef, int index, char *rootname) { + + pmHDU *hdu = pmHDUFromCell(readoutRef->parent); + psString name = NULL; + psStringAppend(&name, "%s_%03d.fits", rootname, index); + pmStackVisualPlotTestImage(readoutFake->image, name); + psFits *fits = psFitsOpen(name, "w"); + psFree(name); + psBinaryOp(readoutFake->image, readoutConv->image, "-", readoutFake->image); + psFitsWriteImage(fits, hdu->header, readoutFake->image, 0, NULL); + psFitsClose(fits); + return true; +} + +// perform the bulk of the PSF-matching +bool matchKernel(pmConfig *config, pmReadout *readoutCnv, pmReadout *readoutRaw, psphotStackOptions *options, int index) { + + bool mdok; + + psAssert(options->psf, "Require target PSF"); + psAssert(options->sourceLists && options->sourceLists->data[index], "Require source list"); + + psMetadata *stackRecipe = psMetadataLookupMetadata(NULL, config->recipes, "PPSTACK"); // ppStack recipe + psAssert(stackRecipe, "We've thrown an error on this before."); + + // Look up appropriate values from the ppSub recipe + psMetadata *subRecipe = psMetadataLookupMetadata(NULL, config->recipes, "PPSUB"); // PPSUB recipe + psAssert(subRecipe, "recipe missing"); + + psString maskValStr = psMetadataLookupStr(NULL, subRecipe, "MASK.VAL"); // Name of bits to mask going in + psString maskPoorStr = psMetadataLookupStr(NULL, stackRecipe, "MASK.POOR"); // Name of bits to mask for poor + psString maskBadStr = psMetadataLookupStr(NULL, stackRecipe, "MASK.BAD"); // Name of bits to mask for bad + + psImageMaskType maskVal = pmConfigMaskGet(maskValStr, config); // Bits to mask going in to pmSubtractionMatch + psImageMaskType maskPoor = pmConfigMaskGet(maskPoorStr, config); // Bits to mask for poor pixels + psImageMaskType maskBad = pmConfigMaskGet(maskBadStr, config); // Bits to mask for bad pixels + + float penalty = psMetadataLookupF32(NULL, subRecipe, "PENALTY"); // Penalty for wideness + int threads = psMetadataLookupS32(NULL, config->arguments, "-threads"); // Number of threads + + int order = psMetadataLookupS32(NULL, subRecipe, "SPATIAL.ORDER"); // Spatial polynomial order + float regionSize = psMetadataLookupF32(NULL, subRecipe, "REGION.SIZE"); // Size of iso-kernel regs + float spacing = psMetadataLookupF32(NULL, subRecipe, "STAMP.SPACING"); // Typical stamp spacing + + int footprint = psMetadataLookupS32(NULL, subRecipe, "STAMP.FOOTPRINT"); // Stamp half-size + int size = psMetadataLookupS32(NULL, subRecipe, "KERNEL.SIZE"); // Kernel half-size + + float threshold = psMetadataLookupF32(NULL, subRecipe, "STAMP.THRESHOLD"); // Threshold for stmps + int stride = psMetadataLookupS32(NULL, subRecipe, "STRIDE"); // Size of convolution patches + int iter = psMetadataLookupS32(NULL, subRecipe, "ITER"); // Rejection iterations + float rej = psMetadataLookupF32(NULL, subRecipe, "REJ"); // Rejection threshold + float kernelError = psMetadataLookupF32(NULL, subRecipe, "KERNEL.ERR"); // Relative systematic error in kernel + float normFrac = psMetadataLookupF32(NULL, subRecipe, "NORM.FRAC"); // Fraction of window for normalisn windw + float sysError = psMetadataLookupF32(NULL, subRecipe, "SYS.ERR"); // Relative systematic error in images + float skyErr = psMetadataLookupF32(NULL, subRecipe, "SKY.ERR"); // Additional error in sky + float covarFrac = psMetadataLookupF32(NULL, subRecipe, "COVAR.FRAC"); // Fraction for covariance calculation + + const char *typeStr = psMetadataLookupStr(NULL, subRecipe, "KERNEL.TYPE"); // Kernel type + pmSubtractionKernelsType type = pmSubtractionKernelsTypeFromString(typeStr); // Kernel type + psVector *widths = psMetadataLookupPtr(NULL, subRecipe, "ISIS.WIDTHS"); // ISIS Gaussian widths + psVector *orders = psMetadataLookupPtr(NULL, subRecipe, "ISIS.ORDERS"); // ISIS Polynomial orders + int inner = psMetadataLookupS32(NULL, subRecipe, "INNER"); // Inner radius + int ringsOrder = psMetadataLookupS32(NULL, subRecipe, "RINGS.ORDER"); // RINGS polynomial order + int binning = psMetadataLookupS32(NULL, subRecipe, "SPAM.BINNING"); // Binning for SPAM kernel + float badFrac = psMetadataLookupF32(NULL, subRecipe, "BADFRAC"); // Maximum bad fraction + float optThresh = psMetadataLookupF32(&mdok, subRecipe, "OPTIMUM.TOL"); // Tolerance for search + int optOrder = psMetadataLookupS32(&mdok, subRecipe, "OPTIMUM.ORDER"); // Order for search + float poorFrac = psMetadataLookupF32(&mdok, subRecipe, "POOR.FRACTION"); // Fraction for "poor" + + bool scale = psMetadataLookupBool(NULL, subRecipe, "SCALE"); // Scale kernel parameters? + float scaleRef = psMetadataLookupF32(NULL, subRecipe, "SCALE.REF"); // Reference for scaling + float scaleMin = psMetadataLookupF32(NULL, subRecipe, "SCALE.MIN"); // Minimum for scaling + float scaleMax = psMetadataLookupF32(NULL, subRecipe, "SCALE.MAX"); // Maximum for scaling + if (!isfinite(scaleRef) || !isfinite(scaleMin) || !isfinite(scaleMax)) { + psError(PSPHOT_ERR_CONFIG, false, + "Scale parameters (SCALE.REF=%f, SCALE.MIN=%f, SCALE.MAX=%f) not set in PPSUB recipe.", + scaleRef, scaleMin, scaleMax); + return false; + } + + // These values are specified specifically for stacking + const char *stampsName = psMetadataLookupStr(NULL, config->arguments, "STAMPS");// Stamps filename + + psVector *widthsCopy = NULL; + psVector *optWidths = NULL; + pmReadout *fake = NULL; + psArray *stampSources = NULL; + + bool optimum = false; + optWidths = SetOptWidths(&optimum, subRecipe); // Vector with FWHMs for optimum search + + // For the sake of stamps, remove nearby sources + stampSources = stackSourcesFilter(options->sourceLists->data[index], footprint); // Filtered list of sources + + fake = makeFakeReadout(config, readoutRaw, stampSources, options->psf, maskVal | maskBad, footprint + size); + if (!fake) goto escape; + + dumpImage(fake, readoutRaw, index, "fake"); + dumpImage(readoutRaw, readoutRaw, index, "real"); + + if (threads) pmSubtractionThreadsInit(); + + // Do the image matching + pmSubtractionKernels *kernel = psMetadataLookupPtr(&mdok, readoutRaw->analysis, PM_SUBTRACTION_ANALYSIS_KERNEL); // Conv kernel + if (kernel) { + if (!pmSubtractionMatchPrecalc(NULL, readoutCnv, fake, readoutRaw, readoutRaw->analysis, stride, kernelError, covarFrac, maskVal, maskBad, maskPoor, poorFrac, badFrac)) { + psError(psErrorCodeLast(), false, "Unable to convolve images."); + goto escape; + } + } else { + // Scale the input parameters + widthsCopy = psVectorCopy(NULL, widths, PS_TYPE_F32); // Copy of kernel widths + if (scale && !pmSubtractionParamsScale(&size, &footprint, widthsCopy, options->inputSeeing->data.F32[index], options->targetSeeing, scaleRef, scaleMin, scaleMax)) { + psError(psErrorCodeLast(), false, "Unable to scale kernel parameters"); + goto escape; + } + + if (!pmSubtractionMatch(NULL, readoutCnv, fake, readoutRaw, footprint, stride, regionSize, spacing, threshold, stampSources, stampsName, type, size, order, widthsCopy, orders, inner, ringsOrder, binning, penalty, optimum, optWidths, optOrder, optThresh, iter, rej, normFrac, sysError, skyErr, kernelError, covarFrac, maskVal, maskBad, maskPoor, poorFrac, badFrac, PM_SUBTRACTION_MODE_2)) { + psError(psErrorCodeLast(), false, "Unable to match images."); + goto escape; + } + } + + // Reject image completely if the maximum deconvolution fraction exceeds the limit + float deconvLimit = psMetadataLookupF32(NULL, stackRecipe, "DECONV.LIMIT"); // Limit on deconvolution fraction + float deconv = psMetadataLookupF32(NULL, readoutCnv->analysis, PM_SUBTRACTION_ANALYSIS_DECONV_MAX); // Max deconvolution fraction + if (deconv > deconvLimit) { + psWarning("Maximum deconvolution fraction (%f) exceeds limit (%f) --- rejecting image %d\n", deconv, deconvLimit, index); + goto escape; + } + + dumpImage(readoutCnv, readoutRaw, index, "conv"); + dumpImageDiff(readoutCnv, fake, readoutRaw, index, "diff"); + + psFree(fake); + psFree(optWidths); + psFree(stampSources); + psFree(widthsCopy); + pmSubtractionThreadsFinalize(); + return true; + +escape: + psFree(fake); + psFree(optWidths); + psFree(stampSources); + psFree(widthsCopy); + pmSubtractionThreadsFinalize(); + return false; +} + +// Extract the regions and solutions used in the image matching +// This stops them from being freed when we iterate back up the FPA +// Record the chi-square value +// XXX this function may not be needed for psphotStack +bool saveMatchData (pmReadout *readout, psphotStackOptions *options, int index) { + + psArray *regions = options->regions->data[index] = psArrayAllocEmpty(ARRAY_BUFFER); // Match regions { - pmHDU *hdu = pmHDUFromCell(readout->parent); - psString name = NULL; - psStringAppend(&name, "convolved_%03d.fits", index); - pmStackVisualPlotTestImage(readout->image, name); - psFits *fits = psFitsOpen(name, "w"); - psFree(name); - psFitsWriteImage(fits, hdu->header, readout->image, 0, NULL); - psFitsClose(fits); - } - -bool matchKernel() { - // Normal operations here - psAssert(options->psf, "Require target PSF"); - psAssert(options->sourceLists && options->sourceLists->data[index], "Require source list"); - - int order = psMetadataLookupS32(NULL, subRecipe, "SPATIAL.ORDER"); // Spatial polynomial order - float regionSize = psMetadataLookupF32(NULL, subRecipe, "REGION.SIZE"); // Size of iso-kernel regs - float spacing = psMetadataLookupF32(NULL, subRecipe, "STAMP.SPACING"); // Typical stamp spacing - int footprint = psMetadataLookupS32(NULL, subRecipe, "STAMP.FOOTPRINT"); // Stamp half-size - float threshold = psMetadataLookupF32(NULL, subRecipe, "STAMP.THRESHOLD"); // Threshold for stmps - int stride = psMetadataLookupS32(NULL, subRecipe, "STRIDE"); // Size of convolution patches - int iter = psMetadataLookupS32(NULL, subRecipe, "ITER"); // Rejection iterations - float rej = psMetadataLookupF32(NULL, subRecipe, "REJ"); // Rejection threshold - float kernelError = psMetadataLookupF32(NULL, subRecipe, "KERNEL.ERR"); // Relative systematic error in kernel - float normFrac = psMetadataLookupF32(NULL, subRecipe, "NORM.FRAC"); // Fraction of window for normalisn windw - float sysError = psMetadataLookupF32(NULL, subRecipe, "SYS.ERR"); // Relative systematic error in images - float skyErr = psMetadataLookupF32(NULL, subRecipe, "SKY.ERR"); // Additional error in sky - float covarFrac = psMetadataLookupF32(NULL, subRecipe, "COVAR.FRAC"); // Fraction for covariance calculation - - const char *typeStr = psMetadataLookupStr(NULL, subRecipe, "KERNEL.TYPE"); // Kernel type - pmSubtractionKernelsType type = pmSubtractionKernelsTypeFromString(typeStr); // Kernel type - psVector *widths = psMetadataLookupPtr(NULL, subRecipe, "ISIS.WIDTHS"); // ISIS Gaussian widths - psVector *orders = psMetadataLookupPtr(NULL, subRecipe, "ISIS.ORDERS"); // ISIS Polynomial orders - int inner = psMetadataLookupS32(NULL, subRecipe, "INNER"); // Inner radius - int ringsOrder = psMetadataLookupS32(NULL, subRecipe, "RINGS.ORDER"); // RINGS polynomial order - int binning = psMetadataLookupS32(NULL, subRecipe, "SPAM.BINNING"); // Binning for SPAM kernel - float badFrac = psMetadataLookupF32(NULL, subRecipe, "BADFRAC"); // Maximum bad fraction - bool optimum = psMetadataLookupBool(&mdok, subRecipe, "OPTIMUM"); // Derive optimum parameters? - float optMin = psMetadataLookupF32(&mdok, subRecipe, "OPTIMUM.MIN"); // Minimum width for search - float optMax = psMetadataLookupF32(&mdok, subRecipe, "OPTIMUM.MAX"); // Maximum width for search - float optStep = psMetadataLookupF32(&mdok, subRecipe, "OPTIMUM.STEP"); // Step for search - float optThresh = psMetadataLookupF32(&mdok, subRecipe, "OPTIMUM.TOL"); // Tolerance for search - int optOrder = psMetadataLookupS32(&mdok, subRecipe, "OPTIMUM.ORDER"); // Order for search - float poorFrac = psMetadataLookupF32(&mdok, subRecipe, "POOR.FRACTION"); // Fraction for "poor" - - bool scale = psMetadataLookupBool(NULL, subRecipe, "SCALE"); // Scale kernel parameters? - float scaleRef = psMetadataLookupF32(NULL, subRecipe, "SCALE.REF"); // Reference for scaling - float scaleMin = psMetadataLookupF32(NULL, subRecipe, "SCALE.MIN"); // Minimum for scaling - float scaleMax = psMetadataLookupF32(NULL, subRecipe, "SCALE.MAX"); // Maximum for scaling - if (!isfinite(scaleRef) || !isfinite(scaleMin) || !isfinite(scaleMax)) { - psError(PPSTACK_ERR_CONFIG, false, - "Scale parameters (SCALE.REF=%f, SCALE.MIN=%f, SCALE.MAX=%f) not set in PPSUB recipe.", - scaleRef, scaleMin, scaleMax); - return false; - } - - - // These values are specified specifically for stacking - const char *stampsName = psMetadataLookupStr(NULL, config->arguments, "STAMPS");// Stamps filename - - psVector *optWidths = NULL; // Vector with FWHMs for optimum search - if (optimum) { - optWidths = psVectorCreate(optWidths, optMin, optMax, optStep, PS_TYPE_F32); - } - - pmReadout *fake = pmReadoutAlloc(NULL); // Fake readout with target PSF - - psStats *bg = psStatsAlloc(PS_STAT_ROBUST_STDEV); // Statistics for background - psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS); // Random number generator - if (!psImageBackground(bg, NULL, readout->image, readout->mask, maskVal | maskBad, rng)) { - psError(PPSTACK_ERR_DATA, false, "Can't measure background for image."); - psFree(fake); - psFree(optWidths); - psFree(conv); - psFree(bg); - psFree(rng); - return false; - } - float minFlux = NOISE_FRACTION * bg->robustStdev; // Minimum flux level for fake image + psString regex = NULL; // Regular expression + psStringAppend(®ex, "^%s$", PM_SUBTRACTION_ANALYSIS_REGION); + psMetadataIterator *iter = psMetadataIteratorAlloc(readout->analysis, PS_LIST_HEAD, regex); + 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); + } + + psArray *kernels = options->kernels->data[index] = psArrayAllocEmpty(ARRAY_BUFFER); // Match kernels + { + psString regex = NULL; // Regular expression + psStringAppend(®ex, "^%s$", PM_SUBTRACTION_ANALYSIS_KERNEL); + psMetadataIterator *iter = psMetadataIteratorAlloc(readout->analysis, PS_LIST_HEAD, regex); + psFree(regex); + psMetadataItem *item = NULL;// Item from iteration + while ((item = psMetadataGetAndIncrement(iter))) { + assert(item->type == PS_DATA_UNKNOWN); + pmSubtractionKernels *kernel = item->data.V; // Kernel used in subtraction + kernels = psArrayAdd(kernels, ARRAY_BUFFER, kernel); + } + psFree(iter); + } + psAssert((regions)->n == (kernels)->n, "Number of match regions and kernels should match"); + + // Record chi^2 + { + double sum = 0.0; // Sum of chi^2 + int num = 0; // Number of measurements of chi^2 + psString regex = NULL; // Regular expression + psStringAppend(®ex, "^%s$", PM_SUBTRACTION_ANALYSIS_KERNEL); + psMetadataIterator *iter = psMetadataIteratorAlloc(readout->analysis, PS_LIST_HEAD, regex); + psFree(regex); + psMetadataItem *item = NULL;// Item from iteration + while ((item = psMetadataGetAndIncrement(iter))) { + assert(item->type == PS_DATA_UNKNOWN); + pmSubtractionKernels *kernels = item->data.V; // Convolution kernels + sum += kernels->mean; + num++; + } + psFree(iter); + options->matchChi2->data.F32[index] = sum / (psImageCovarianceFactor(readout->covariance) * num); + } + + return true; +} + +// Kernel normalisation for convolved readout +bool renormKernel(pmReadout *readout, psphotStackOptions *options, int index) { + + double sum = 0.0; // Sum of chi^2 + int num = 0; // Number of measurements of chi^2 + psString regex = NULL; // Regular expression + psStringAppend(®ex, "^%s$", PM_SUBTRACTION_ANALYSIS_NORM); + psMetadataIterator *iter = psMetadataIteratorAlloc(readout->analysis, PS_LIST_HEAD, regex); + psFree(regex); + psMetadataItem *item = NULL;// Item from iteration + while ((item = psMetadataGetAndIncrement(iter))) { + assert(item->type == PS_TYPE_F32); + float norm = item->data.F32; // Normalisation + sum += norm; + num++; + } + psFree(iter); + float conv = sum/num; // Mean normalisation from convolution + float stars = powf(10.0, -0.4 * options->norm->data.F32[index]); // Normalisation from stars + float renorm = stars / conv; // Renormalisation to apply + psLogMsg("psphotStack", PS_LOG_INFO, "Renormalising image %d by %f (kernel: %f, stars: %f)\n", index, renorm, conv, stars); + + psBinaryOp(readout->image, readout->image, "*", psScalarAlloc(renorm, PS_TYPE_F32)); + psBinaryOp(readout->variance, readout->variance, "*", psScalarAlloc(PS_SQR(renorm), PS_TYPE_F32)); + return true; +} + +// adjust scaling for readout (remove background, ..., determine weighting) +bool rescaleData(pmReadout *readout, pmConfig *config, psphotStackOptions *options, int index) { + + psMetadata *stackRecipe = psMetadataLookupMetadata(NULL, config->recipes, "PPSTACK"); // ppStack recipe + psAssert(stackRecipe, "We've thrown an error on this before."); + + // Look up appropriate values from the ppSub recipe + psMetadata *subRecipe = psMetadataLookupMetadata(NULL, config->recipes, "PPSUB"); // PPSUB recipe + psAssert(subRecipe, "recipe missing"); + + psString maskValStr = psMetadataLookupStr(NULL, subRecipe, "MASK.VAL"); // Name of bits to mask going in + psString maskBadStr = psMetadataLookupStr(NULL, stackRecipe, "MASK.BAD"); // Name of bits to mask for bad + + psImageMaskType maskVal = pmConfigMaskGet(maskValStr, config); // Bits to mask going in to pmSubtractionMatch + psImageMaskType maskBad = pmConfigMaskGet(maskBadStr, config); // Bits to mask for bad pixels + + // Ensure the background value is zero + psStats *bg = psStatsAlloc(PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV); // Statistics for background + psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS); // Random number generator + + // XXX why is this in config->arguments and not recipe? + if (!psMetadataLookupBool(NULL, config->arguments, "PPSTACK.SKIP.BG.SUB")) { + if (!psImageBackground(bg, NULL, readout->image, readout->mask, maskVal | maskBad, rng)) { + psAbort("Can't measure background for image."); + // XXX we used to clear error: why is this acceptable? psErrorClear(); + } + + float value = psStatsGetValue(bg, PS_STAT_ROBUST_MEDIAN); + float stdev = psStatsGetValue(bg, PS_STAT_ROBUST_STDEV); + + psLogMsg("psphotStack", PS_LOG_INFO, "Correcting convolved image background by %lf (+/- %lf)", value, stdev); + psBinaryOp(readout->image, readout->image, "-", psScalarAlloc(value, PS_TYPE_F32)); + } + + if (!stackRenormaliseReadout(config, readout)) { + psFree(rng); + psFree(bg); + return false; + } + + // Measure the variance level for the weighting + if (psMetadataLookupBool(NULL, stackRecipe, "WEIGHTS")) { + if (!psImageBackground(bg, NULL, readout->variance, readout->mask, maskVal | maskBad, rng)) { + psError(PSPHOT_ERR_DATA, false, "Can't measure mean variance for image."); psFree(rng); psFree(bg); - - // For the sake of stamps, remove nearby sources - psArray *stampSources = stackSourcesFilter(options->sourceLists->data[index], - footprint); // Filtered list of sources - - bool oldThreads = pmReadoutFakeThreads(true); // Old threading state - if (!pmReadoutFakeFromSources(fake, readout->image->numCols, readout->image->numRows, - stampSources, SOURCE_MASK, NULL, NULL, options->psf, - minFlux, footprint + size, false, true)) { - psError(PPSTACK_ERR_DATA, false, "Unable to generate fake image with target PSF."); - psFree(fake); - psFree(optWidths); - psFree(conv); - return false; - } - pmReadoutFakeThreads(oldThreads); - - fake->mask = psImageCopy(NULL, readout->mask, PS_TYPE_IMAGE_MASK); - - // Add the background into the target image - psImage *bgImage = stackBackgroundModel(readout, config); // Image of background - psBinaryOp(fake->image, fake->image, "+", bgImage); - psFree(bgImage); - - dumpImage(); - - if (threads > 0) { - pmSubtractionThreadsInit(); - } - - // Do the image matching - pmSubtractionKernels *kernel = psMetadataLookupPtr(&mdok, readout->analysis, PM_SUBTRACTION_ANALYSIS_KERNEL); // Conv kernel - if (kernel) { - if (!pmSubtractionMatchPrecalc(NULL, conv, fake, readout, readout->analysis, - stride, kernelError, covarFrac, maskVal, maskBad, maskPoor, - poorFrac, badFrac)) { - psError(psErrorCodeLast(), false, "Unable to convolve images."); - psFree(fake); - psFree(optWidths); - psFree(stampSources); - psFree(conv); - if (threads > 0) { - pmSubtractionThreadsFinalize(); - } - return false; - } - } else { - // Scale the input parameters - psVector *widthsCopy = psVectorCopy(NULL, widths, PS_TYPE_F32); // Copy of kernel widths - if (scale && !pmSubtractionParamsScale(&size, &footprint, widthsCopy, - options->inputSeeing->data.F32[index], - options->targetSeeing, scaleRef, scaleMin, scaleMax)) { - psError(psErrorCodeLast(), false, "Unable to scale kernel parameters"); - psFree(fake); - psFree(optWidths); - psFree(stampSources); - psFree(conv); - psFree(widthsCopy); - if (threads > 0) { - pmSubtractionThreadsFinalize(); - } - return false; - } - - if (!pmSubtractionMatch(NULL, conv, fake, readout, footprint, stride, regionSize, spacing, - threshold, stampSources, stampsName, type, size, order, widthsCopy, - orders, inner, ringsOrder, binning, penalty, - optimum, optWidths, optOrder, optThresh, iter, rej, normFrac, - sysError, skyErr, kernelError, covarFrac, maskVal, maskBad, maskPoor, - poorFrac, badFrac, PM_SUBTRACTION_MODE_2)) { - psError(psErrorCodeLast(), false, "Unable to match images."); - psFree(fake); - psFree(optWidths); - psFree(stampSources); - psFree(conv); - psFree(widthsCopy); - if (threads > 0) { - pmSubtractionThreadsFinalize(); - } - return false; - } - psFree(widthsCopy); - } - - dumpImage2(); - - psFree(fake); - psFree(optWidths); - psFree(stampSources); - - if (threads > 0) { - pmSubtractionThreadsFinalize(); - } - - // Replace original images with convolved - psFree(readout->image); - psFree(readout->mask); - psFree(readout->variance); - psFree(readout->covariance); - readout->image = psMemIncrRefCounter(conv->image); - readout->mask = psMemIncrRefCounter(conv->mask); - readout->variance = psMemIncrRefCounter(conv->variance); - readout->covariance = psImageCovarianceTruncate(conv->covariance, COVAR_FRAC); - -} - -bool saveMatchData () { - // Extract the regions and solutions used in the image matching - // This stops them from being freed when we iterate back up the FPA - psArray *regions = options->regions->data[index] = psArrayAllocEmpty(ARRAY_BUFFER); // Match regions - { - psString regex = NULL; // Regular expression - psStringAppend(®ex, "^%s$", PM_SUBTRACTION_ANALYSIS_REGION); - psMetadataIterator *iter = psMetadataIteratorAlloc(conv->analysis, PS_LIST_HEAD, regex); - 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); + return false; } - psArray *kernels = options->kernels->data[index] = psArrayAllocEmpty(ARRAY_BUFFER); // Match kernels - { - psString regex = NULL; // Regular expression - psStringAppend(®ex, "^%s$", PM_SUBTRACTION_ANALYSIS_KERNEL); - psMetadataIterator *iter = psMetadataIteratorAlloc(conv->analysis, PS_LIST_HEAD, regex); - psFree(regex); - psMetadataItem *item = NULL;// Item from iteration - while ((item = psMetadataGetAndIncrement(iter))) { - assert(item->type == PS_DATA_UNKNOWN); - pmSubtractionKernels *kernel = item->data.V; // Kernel used in subtraction - kernels = psArrayAdd(kernels, ARRAY_BUFFER, kernel); - } - psFree(iter); - } - psAssert((regions)->n == (kernels)->n, "Number of match regions and kernels should match"); -} - -bool saveChiSquare() { - // Record chi^2 - { - double sum = 0.0; // Sum of chi^2 - int num = 0; // Number of measurements of chi^2 - psString regex = NULL; // Regular expression - psStringAppend(®ex, "^%s$", PM_SUBTRACTION_ANALYSIS_KERNEL); - psMetadataIterator *iter = psMetadataIteratorAlloc(conv->analysis, PS_LIST_HEAD, regex); - psFree(regex); - psMetadataItem *item = NULL;// Item from iteration - while ((item = psMetadataGetAndIncrement(iter))) { - assert(item->type == PS_DATA_UNKNOWN); - pmSubtractionKernels *kernels = item->data.V; // Convolution kernels - sum += kernels->mean; - num++; - } - psFree(iter); - options->matchChi2->data.F32[index] = sum / (psImageCovarianceFactor(readout->covariance) * num); - } - -} - -bool renormKernel() { - // Kernel normalisation - { - double sum = 0.0; // Sum of chi^2 - int num = 0; // Number of measurements of chi^2 - psString regex = NULL; // Regular expression - psStringAppend(®ex, "^%s$", PM_SUBTRACTION_ANALYSIS_NORM); - psMetadataIterator *iter = psMetadataIteratorAlloc(conv->analysis, PS_LIST_HEAD, regex); - psFree(regex); - psMetadataItem *item = NULL;// Item from iteration - while ((item = psMetadataGetAndIncrement(iter))) { - assert(item->type == PS_TYPE_F32); - float norm = item->data.F32; // Normalisation - sum += norm; - num++; - } - psFree(iter); - float conv = sum/num; // Mean normalisation from convolution - float stars = powf(10.0, -0.4 * options->norm->data.F32[index]); // Normalisation from stars - float renorm = stars / conv; // Renormalisation to apply - psLogMsg("ppStack", PS_LOG_INFO, "Renormalising image %d by %f (kernel: %f, stars: %f)\n", - index, renorm, conv, stars); - psBinaryOp(readout->image, readout->image, "*", psScalarAlloc(renorm, PS_TYPE_F32)); - psBinaryOp(readout->variance, readout->variance, "*", psScalarAlloc(PS_SQR(renorm), PS_TYPE_F32)); - } - -} + options->weightings->data.F32[index] = 1.0 / (psStatsGetValue(bg, PS_STAT_ROBUST_MEDIAN) * psImageCovarianceFactor(readout->covariance)); + } else { + options->weightings->data.F32[index] = 1.0; + } + psLogMsg("psphotStack", PS_LOG_INFO, "Weighting for image %d is %f\n", index, options->weightings->data.F32[index]); + + psFree(rng); + psFree(bg); + return true; +} + +# define NOISE_FRACTION 0.01 // Set minimum flux to this fraction of noise +# define SOURCE_MASK (PM_SOURCE_MODE_FAIL | PM_SOURCE_MODE_DEFECT | PM_SOURCE_MODE_SATURATED | PM_SOURCE_MODE_CR_LIMIT | PM_SOURCE_MODE_EXT_LIMIT) // Mask to apply to input sources + +// generate a fake readout against which to PSF match +pmReadout *makeFakeReadout(pmConfig *config, pmReadout *readoutRaw, psArray *sources, pmPSF *psf, psImageMaskType maskVal, int fullSize) { + + pmReadout *fake = pmReadoutAlloc(NULL); // Fake readout with target PSF + + psStats *bg = psStatsAlloc(PS_STAT_ROBUST_STDEV); // Statistics for background + psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS); // Random number generator + if (!psImageBackground(bg, NULL, readoutRaw->image, readoutRaw->mask, maskVal, rng)) { + psError(PSPHOT_ERR_DATA, false, "Can't measure background for image."); + psFree(fake); + psFree(bg); + psFree(rng); + return NULL; + } + float minFlux = NOISE_FRACTION * bg->robustStdev; // Minimum flux level for fake image + psFree(rng); + psFree(bg); + + bool oldThreads = pmReadoutFakeThreads(true); // Old threading state + if (!pmReadoutFakeFromSources(fake, readoutRaw->image->numCols, readoutRaw->image->numRows, sources, SOURCE_MASK, NULL, NULL, psf, minFlux, fullSize, false, true)) { + psError(PSPHOT_ERR_DATA, false, "Unable to generate fake image with target PSF."); + psFree(fake); + return NULL; + } + pmReadoutFakeThreads(oldThreads); + + fake->mask = psImageCopy(NULL, readoutRaw->mask, PS_TYPE_IMAGE_MASK); + + // Add the background into the target image + psImage *bgImage = stackBackgroundModel(readoutRaw, config); // Image of background + psBinaryOp(fake->image, fake->image, "+", bgImage); + psFree(bgImage); + + return fake; +} + +// set the widths +psVector *SetOptWidths (bool *optimum, psMetadata *recipe) { + + bool status; + + *optimum = psMetadataLookupBool(&status, recipe, "OPTIMUM"); // Derive optimum parameters? + psAssert (status, "missing recipe value %s", "OPTIMUM"); + + psVector *optWidths = NULL; // Vector with FWHMs for optimum search + + if (*optimum) { + float optMin = psMetadataLookupF32(&status, recipe, "OPTIMUM.MIN"); // Minimum width for search + psAssert (status, "missing recipe value %s", "OPTIMUM.MIN"); + + float optMax = psMetadataLookupF32(&status, recipe, "OPTIMUM.MAX"); // Maximum width for search + psAssert (status, "missing recipe value %s", "OPTIMUM.MAX"); + + float optStep = psMetadataLookupF32(&status, recipe, "OPTIMUM.STEP"); // Step for search + psAssert (status, "missing recipe value %s", "OPTIMUM.STEP"); + + optWidths = psVectorCreate(optWidths, optMin, optMax, optStep, PS_TYPE_F32); + } + + return optWidths; +} Index: /branches/eam_branches/psphot.20100506/src/psphotStackOptions.c =================================================================== --- /branches/eam_branches/psphot.20100506/src/psphotStackOptions.c (revision 27876) +++ /branches/eam_branches/psphot.20100506/src/psphotStackOptions.c (revision 27876) @@ -0,0 +1,69 @@ +# include "psphotInternal.h" + +typedef enum { + PSPHOT_CNV_SRC_NONE, + PSPHOT_CNV_SRC_AUTO, + PSPHOT_CNV_SRC_CNV, + PSPHOT_CNV_SRC_RAW, +} psphotStackConvolveSource; + +static void psphotStackOptionsFree (psphotStackOptions *options) { + + if (options == NULL) return; + + psFree (options->psf); + psFree (options->inputSeeing); + psFree (options->sourceLists); + psFree (options->norm); + psFree (options->kernels); + psFree (options->regions); + psFree (options->matchChi2); + psFree (options->weightings); + + return; +} + +psphotStackOptions *psphotStackOptionsAlloc (int num) { + + psphotStackOptions *options = (psphotStackOptions *) psAlloc(sizeof(psphotStackOptions)); + psMemSetDeallocator(options, (psFreeFunc) psphotStackOptionsFree); + + options->numCols = 0; + options->numRows = 0; + + options->num = num; + options->psf = NULL; + options->convolve = false; + + options->psfs = psArrayAlloc(num); + + options->inputSeeing = psVectorAlloc(num, PS_TYPE_F32); + psVectorInit(options->inputSeeing, NAN); + + options->matchChi2 = psVectorAlloc(num, PS_TYPE_F32); + psVectorInit(options->matchChi2, NAN); + + options->inputMask = psVectorAlloc(num, PS_TYPE_VECTOR_MASK); // Mask for inputs + psVectorInit(options->inputMask, 0); + + options->targetSeeing = NAN; + + options->sourceLists = psArrayAlloc(num); // Individual lists of sources for matching + options->norm = psVectorAlloc(num, PS_TYPE_F32); + psVectorInit(options->norm, NAN); + + options->kernels = NULL; + options->regions = NULL; + options->matchChi2 = NULL; + options->weightings = NULL; + + return options; +} + +psphotStackConvolveSource psphotStackConvolveSourceFromString (const char *string) { + + if (!strcasecmp(string, "AUTO")) return PSPHOT_CNV_SRC_AUTO; + if (!strcasecmp(string, "CNV")) return PSPHOT_CNV_SRC_CNV; + if (!strcasecmp(string, "RAW")) return PSPHOT_CNV_SRC_RAW; + return PSPHOT_CNV_SRC_NONE; +} Index: /branches/eam_branches/psphot.20100506/src/psphotStackParseCamera.c =================================================================== --- /branches/eam_branches/psphot.20100506/src/psphotStackParseCamera.c (revision 27875) +++ /branches/eam_branches/psphot.20100506/src/psphotStackParseCamera.c (revision 27876) @@ -25,40 +25,125 @@ psMetadata *input = item->data.md; // The input metadata of interest - // look for 'IMAGE', 'MASK', 'VARIANCE' in folder (only 'IMAGE' is required) - - psString image = psMetadataLookupStr(NULL, input, "IMAGE"); // Name of image - if (!image || strlen(image) == 0) { - psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Component %s lacks IMAGE of type STR", item->name); + pmFPAfile *rawInputFile = NULL; + pmFPAfile *cnvInputFile = NULL; + + // RAW (unconvolved) input data (RAW:IMAGE, RAW:MASK, RAW:VARIANCE, RAW:PSF) + psString rawImage = psMetadataLookupStr(NULL, input, "RAW:IMAGE"); // Name of image + if (rawImage && strlen(rawImage) > 0) { + pmFPAfile *rawInputFile = defineFile(config, NULL, "PSPHOT.STACK.INPUT.RAW", rawImage, PM_FPA_FILE_IMAGE); // File for image + if (!rawInputFile) { + psError(PS_ERR_UNKNOWN, false, "Unable to define file from image %d (%s)", i, rawImage); + return false; + } + psString mask = psMetadataLookupStr(&status, input, "RAW:MASK"); // Name of mask + if (mask && strlen(mask) > 0) { + if (!defineFile(config, rawInputFile, "PSPHOT.STACK.MASK.RAW", mask, PM_FPA_FILE_MASK)) { + psError(PS_ERR_UNKNOWN, false, "Unable to define file from mask %d (%s)", i, mask); + return false; + } + } + psString variance = psMetadataLookupStr(&status, input, "RAW:VARIANCE"); // Name of variance map + if (variance && strlen(variance) > 0) { + if (!defineFile(config, rawInputFile, "PSPHOT.STACK.VARIANCE.RAW", variance, PM_FPA_FILE_VARIANCE)) { + psError(PS_ERR_UNKNOWN, false, "Unable to define file from variance %d (%s)", i, variance); + return false; + } + } + psString psf = psMetadataLookupStr(&status, input, "RAW:PSF"); // Name of mask + if (psf && strlen(psf) > 0) { + if (!defineFile(config, rawInputFile, "PSPHOT.STACK.PSF.RAW", psf, PM_FPA_FILE_PSF)) { + psError(PS_ERR_UNKNOWN, false, "Unable to define file from psf %d (%s)", i, psf); + return false; + } + } + } + + // CNV (convolved) input data (CNV:IMAGE, CNV:MASK, CNV:VARIANCE, CNV:PSF) + psString cnvImage = psMetadataLookupStr(NULL, input, "CNV:IMAGE"); // Name of image + if (cnvImage && strlen(cnvImage) > 0) { + pmFPAfile *cnvInputFile = defineFile(config, NULL, "PSPHOT.STACK.INPUT.CNV", cnvImage, PM_FPA_FILE_IMAGE); // File for image + if (!cnvInputFile) { + psError(PS_ERR_UNKNOWN, false, "Unable to define file from image %d (%s)", i, cnvImage); + return false; + } + psString mask = psMetadataLookupStr(&status, input, "CNV:MASK"); // Name of mask + if (mask && strlen(mask) > 0) { + if (!defineFile(config, cnvInputFile, "PSPHOT.STACK.MASK.CNV", mask, PM_FPA_FILE_MASK)) { + psError(PS_ERR_UNKNOWN, false, "Unable to define file from mask %d (%s)", i, mask); + return false; + } + } + psString variance = psMetadataLookupStr(&status, input, "CNV:VARIANCE"); // Name of variance map + if (variance && strlen(variance) > 0) { + if (!defineFile(config, cnvInputFile, "PSPHOT.STACK.VARIANCE.CNV", variance, PM_FPA_FILE_VARIANCE)) { + psError(PS_ERR_UNKNOWN, false, "Unable to define file from variance %d (%s)", i, variance); + return false; + } + } + psString psf = psMetadataLookupStr(&status, input, "CNV:PSF"); // Name of mask + if (psf && strlen(psf) > 0) { + if (!defineFile(config, cnvInputFile, "PSPHOT.STACK.PSF.CNV", psf, PM_FPA_FILE_PSF)) { + psError(PS_ERR_UNKNOWN, false, "Unable to define file from psf %d (%s)", i, psf); + return false; + } + } + } + + if (!rawInputFile && !cnvInputFile) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Component %s (%d) lacks IMAGE of type STR", item->name, i); return false; } - pmFPAfile *imageFile = defineFile(config, NULL, "PSPHOT.INPUT", image, PM_FPA_FILE_IMAGE); // File for image - if (!imageFile) { - psError(PS_ERR_UNKNOWN, false, "Unable to define file from image %d (%s)", i, image); - return false; - } - - psString mask = psMetadataLookupStr(&status, input, "MASK"); // Name of mask - if (mask && strlen(mask) > 0) { - if (!defineFile(config, imageFile, "PSPHOT.MASK", mask, PM_FPA_FILE_MASK)) { - psError(PS_ERR_UNKNOWN, false, "Unable to define file from mask %d (%s)", i, mask); - return false; - } - } - - psString variance = psMetadataLookupStr(&status, input, "VARIANCE"); // Name of variance map - if (variance && strlen(variance) > 0) { - if (!defineFile(config, imageFile, "PSPHOT.VARIANCE", variance, PM_FPA_FILE_VARIANCE)) { - psError(PS_ERR_UNKNOWN, false, "Unable to define file from variance %d (%s)", i, variance); - return false; - } - } - // the output sources are carried on the input->fpa structures - pmFPAfile *outsources = pmFPAfileDefineOutputFromFile (config, imageFile, "PSPHOT.STACK.OUTPUT"); - if (!outsources) { - psError(PSPHOT_ERR_CONFIG, false, "Cannot find a rule for PSPHOT.STACK.OUTPUT"); - return false; - } - outsources->save = true; - outsources->fileID = i; // this is used to generate output names + + psString sources = psMetadataLookupStr(&status, input, "SOURCES"); // Name of mask + pmFPAfile *srcInputFile = rawInputFile ? rawInputFile : cnvInputFile; + if (sources && strlen(sources) > 0) { + if (!defineFile(config, srcInputFile, "PSPHOT.STACK.SOURCES.CNV", sources, PM_FPA_FILE_CMF)) { + psError(PS_ERR_UNKNOWN, false, "Unable to define file from sources %d (%s)", i, sources); + return false; + } + } + + // generate an pmFPAimage for the output convolved image + // XXX output of these files should be optional + { + pmFPAfile *outputImage = pmFPAfileDefineOutput(config, NULL, "PSPHOT.STACK.OUTPUT.IMAGE"); + if (!outputImage) { + psError(PSPHOT_ERR_CONFIG, false, "Trouble defining PSPHOT.STACK.OUTPUT.IMAGE"); + return false; + } + outputImage->save = true; + outsources->fileID = i; // this is used to generate output names + + pmFPAfile *outputMask = pmFPAfileDefineOutput(config, outputImage->fpa, "PSPHOT.STACK.OUTPUT.MASK"); + if (!outputMask) { + psError(PS_ERR_IO, false, _("Unable to generate output file from PSPHOT.STACK.OUTPUT.MASK")); + return NULL; + } + if (outputMask->type != PM_FPA_FILE_MASK) { + psError(PS_ERR_IO, true, "PSPHOT.STACK.OUTPUT.MASK is not of type MASK"); + return NULL; + } + outputMask->save = true; + + pmFPAfile *outputVariance = pmFPAfileDefineOutput(config, outputImage->fpa, "PSPHOT.STACK.OUTPUT.VARIANCE"); + if (!outputVariance) { + psError(PS_ERR_IO, false, _("Unable to generate output file from PSPHOT.STACK.OUTPUT.VARIANCE")); + return NULL; + } + if (outputVariance->type != PM_FPA_FILE_VARIANCE) { + psError(PS_ERR_IO, true, "PSPHOT.STACK.OUTPUT.VARIANCE is not of type VARIANCE"); + return NULL; + } + outputVariance->save = true; + + // the output sources are carried on the outputImage->fpa structures + pmFPAfile *outsources = pmFPAfileDefineOutputFromFile (config, outputImage, "PSPHOT.STACK.OUTPUT"); + if (!outsources) { + psError(PSPHOT_ERR_CONFIG, false, "Cannot find a rule for PSPHOT.STACK.OUTPUT"); + return false; + } + outsources->save = true; + outsources->fileID = i; // this is used to generate output names + } } psMetadataRemoveKey(config->arguments, "FILENAMES"); @@ -71,41 +156,35 @@ // generate an pmFPAimage for the chisqImage - pmFPAfile *chisqImage = pmFPAfileDefineOutput(config, NULL, "PSPHOT.CHISQ.IMAGE"); - if (!chisqImage) { - psError(PSPHOT_ERR_CONFIG, false, "Trouble defining PSPHOT.CHISQ.IMAGE"); - return false; - } - chisqImage->save = true; - - pmFPAfile *chisqMask = pmFPAfileDefineOutput(config, chisqImage->fpa, "PSPHOT.CHISQ.MASK"); - if (!chisqMask) { - psError(PS_ERR_IO, false, _("Unable to generate output file from PSPHOT.CHISQ.MASK")); - return NULL; - } - if (chisqMask->type != PM_FPA_FILE_MASK) { - psError(PS_ERR_IO, true, "PSPHOT.CHISQ.MASK is not of type MASK"); - return NULL; - } - chisqMask->save = true; - - pmFPAfile *chisqVariance = pmFPAfileDefineOutput(config, chisqImage->fpa, "PSPHOT.CHISQ.VARIANCE"); - if (!chisqVariance) { - psError(PS_ERR_IO, false, _("Unable to generate output file from PSPHOT.CHISQ.VARIANCE")); - return NULL; - } - if (chisqVariance->type != PM_FPA_FILE_VARIANCE) { - psError(PS_ERR_IO, true, "PSPHOT.CHISQ.VARIANCE is not of type VARIANCE"); - return NULL; - } - chisqVariance->save = true; - -# if (0) - // define the additional input/output files associated with psphot - // XXX figure out which files are needed by psphotStack - if (false && !psphotDefineFiles (config, input)) { - psError(PSPHOT_ERR_CONFIG, false, "Trouble defining the additional input/output files"); - return false; - } -# endif + // XXX output of these files should be optional + { + pmFPAfile *chisqImage = pmFPAfileDefineOutput(config, NULL, "PSPHOT.CHISQ.IMAGE"); + if (!chisqImage) { + psError(PSPHOT_ERR_CONFIG, false, "Trouble defining PSPHOT.CHISQ.IMAGE"); + return false; + } + chisqImage->save = true; + + pmFPAfile *chisqMask = pmFPAfileDefineOutput(config, chisqImage->fpa, "PSPHOT.CHISQ.MASK"); + if (!chisqMask) { + psError(PS_ERR_IO, false, _("Unable to generate output file from PSPHOT.CHISQ.MASK")); + return NULL; + } + if (chisqMask->type != PM_FPA_FILE_MASK) { + psError(PS_ERR_IO, true, "PSPHOT.CHISQ.MASK is not of type MASK"); + return NULL; + } + chisqMask->save = true; + + pmFPAfile *chisqVariance = pmFPAfileDefineOutput(config, chisqImage->fpa, "PSPHOT.CHISQ.VARIANCE"); + if (!chisqVariance) { + psError(PS_ERR_IO, false, _("Unable to generate output file from PSPHOT.CHISQ.VARIANCE")); + return NULL; + } + if (chisqVariance->type != PM_FPA_FILE_VARIANCE) { + psError(PS_ERR_IO, true, "PSPHOT.CHISQ.VARIANCE is not of type VARIANCE"); + return NULL; + } + chisqVariance->save = true; + } psTrace("psphot", 1, "Done with psphotStackParseCamera...\n"); @@ -145,2 +224,29 @@ } + +/*** + * + * psphotStack : + + * * inputs: + * * unconvolved images + * * raw convolved images + * * psfs (unconvolved or convolved?) + * * sources + + * optionally convolve the unconvolved or the raw inputs + * optionally perform no convolutions + * optionally save the psf-matched images + + */ + + +# if (0) + // define the additional input/output files associated with psphot + // XXX figure out which files are needed by psphotStack + if (false && !psphotDefineFiles (config, input)) { + psError(PSPHOT_ERR_CONFIG, false, "Trouble defining the additional input/output files"); + return false; + } +# endif +