Index: trunk/psphot/src/psphotStackMatchPSFsUtils.c =================================================================== --- trunk/psphot/src/psphotStackMatchPSFsUtils.c (revision 27850) +++ trunk/psphot/src/psphotStackMatchPSFsUtils.c (revision 28013) @@ -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; @@ -51,6 +45,7 @@ } +# define SN_MIN 50.0 psArray *stackSourcesFilter(psArray *sources, // Source list to filter - int exclusion // Exclusion zone, pixels + int exclusion // Exclusion zone, pixels ) { @@ -68,4 +63,6 @@ continue; } + if (!source->peak) continue; + if (source->peak->SN < SN_MIN) continue; coordsFromSource(&x->data.F32[numGood], &y->data.F32[numGood], source); numGood++; @@ -83,4 +80,6 @@ continue; } + if (!source->peak) continue; + if (source->peak->SN < SN_MIN) continue; float xSource, ySource; // Coordinates of source coordsFromSource(&xSource, &ySource, source); @@ -90,5 +89,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 +103,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 +110,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 +120,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 +137,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 +152,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 +161,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 +179,4 @@ psImageCovarianceTransfer(readout->variance, readout->covariance); return pmReadoutVarianceRenormalise(readout, maskBad, num, minValid, maxValid); -#else - return true; -#endif } @@ -190,367 +185,424 @@ // 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() +# 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) { + + 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 *readoutOut, pmReadout *readoutSrc, 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(&mdok, 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, readoutSrc, stampSources, options->psf, maskVal | maskBad, footprint + size); + if (!fake) goto escape; + + dumpImage(fake, readoutSrc, index, "fake"); + dumpImage(readoutSrc, readoutSrc, index, "real"); + + if (threads) pmSubtractionThreadsInit(); + + // Do the image matching + pmSubtractionKernels *kernel = psMetadataLookupPtr(&mdok, readoutSrc->analysis, PM_SUBTRACTION_ANALYSIS_KERNEL); // Conv kernel + if (kernel) { + if (!pmSubtractionMatchPrecalc(NULL, readoutOut, fake, readoutSrc, readoutSrc->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, readoutOut, fake, readoutSrc, 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, readoutOut->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(readoutOut, readoutSrc, index, "conv"); + dumpImageDiff(readoutOut, fake, readoutSrc, 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 *readoutSrc, 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, readoutSrc->image, readoutSrc->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, readoutSrc->image->numCols, readoutSrc->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, readoutSrc->mask, PS_TYPE_IMAGE_MASK); + + // Add the background into the target image + psImage *bgImage = stackBackgroundModel(readoutSrc, 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; +}