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