Index: trunk/psphot/src/psphotStackMatchPSFsUtils.c =================================================================== --- trunk/psphot/src/psphotStackMatchPSFsUtils.c (revision 31154) +++ trunk/psphot/src/psphotStackMatchPSFsUtils.c (revision 32348) @@ -2,30 +2,10 @@ # define ARRAY_BUFFER 16 // Number to add to array at a time -// XXX better name -bool readImage(psImage **target, // Target for image - const char *name, // Name of FITS file - const pmConfig *config // Configuration - ) -{ - psString resolved = pmConfigConvertFilename(name, config, false, false); // Resolved filename - psFits *fits = psFitsOpen(resolved, "r"); - psFree(resolved); - if (!fits) { - 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(PSPHOT_ERR_IO, false, "Unable to read previously produced image: %s", name); - psFitsClose(fits); - return false; - } - psFitsClose(fits); - - psFree(*target); - *target = image; - - return true; -} +psVector *SetOptWidths (bool *optimum, psMetadata *recipe); +pmReadout *makeFakeReadout(pmConfig *config, pmReadout *raw, psArray *sources, pmPSF *psf, psImageMaskType maskVal, int fullSize); +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); // Get coordinates from a source @@ -148,6 +128,6 @@ // Renormalise a readout's variance map -bool stackRenormaliseReadout(const pmConfig *config, // Configuration - pmReadout *readout // Readout to renormalise +bool psphotStackRenormaliseVariance(const pmConfig *config, // Configuration + pmReadout *readout // Readout to renormalise ) { @@ -180,59 +160,4 @@ return pmReadoutVarianceRenormalise(readout, maskBad, num, minValid, maxValid); } - -// This is a hack to use the temporary convolved images and kernel generated previously. -// This makes the 'matching' operation much faster, allowing debugging of the stack process easier. -// It implicitly assumes the output root name is the same between invocations. - -# if (0) -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; -} -# endif bool dumpImage(pmReadout *readoutOut, pmReadout *readoutRef, int index, char *rootname) { @@ -362,6 +287,6 @@ widthsCopy = psVectorCopy(NULL, widths, PS_TYPE_F32); // Copy of kernel widths - // we need to register the FWHM values for use downstream - pmSubtractionSetFWHMs(options->targetSeeing, options->inputSeeing->data.F32[index]); + // we need to register the FWHM values for use by pmSubtraction code + pmSubtractionSetFWHMs(options->targetSeeing->data.F32[0], options->inputSeeing->data.F32[index]); pmSubtractionParamScaleOptions(scale, scaleRef, scaleMin, scaleMax); @@ -469,90 +394,4 @@ } -// 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); - 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("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