Index: trunk/ppStack/src/ppStackMatch.c =================================================================== --- trunk/ppStack/src/ppStackMatch.c (revision 21183) +++ trunk/ppStack/src/ppStackMatch.c (revision 21199) @@ -6,4 +6,5 @@ #include #include +#include #include "ppStack.h" @@ -47,4 +48,20 @@ #endif +// Get coordinates from a source +static void coordsFromSource(float *x, float *y, const pmSource *source) +{ + assert(x && y); + assert(source); + + if (source->modelPSF) { + *x = source->modelPSF->params->data.F32[PM_PAR_XPOS]; + *y = source->modelPSF->params->data.F32[PM_PAR_YPOS]; + } else { + *x = source->peak->xf; + *y = source->peak->yf; + } + return; +} + static psArray *stackSourcesFilter(psArray *sources, // Source list to filter int exclusion // Exclusion zone, pixels @@ -64,6 +81,5 @@ continue; } - x->data.F32[numGood] = source->peak->xf; - y->data.F32[numGood] = source->peak->yf; + coordsFromSource(&x->data.F32[numGood], &y->data.F32[numGood], source); numGood++; } @@ -80,6 +96,9 @@ continue; } - coords->data.F64[0] = source->peak->xf; - coords->data.F64[1] = source->peak->yf; + float xSource, ySource; // Coordinates of source + coordsFromSource(&xSource, &ySource, source); + + coords->data.F64[0] = xSource; + coords->data.F64[1] = ySource; long numWithin = psTreeWithin(tree, coords, exclusion); // Number within exclusion zone @@ -101,55 +120,45 @@ } -#define BG_SIZE 64 // Large box half-size in which to measure background - -// Generate a background model of the readout we're matching -psImage *stackBackgroundModel(pmReadout *ro, psImageMaskType maskVal, const psArray *sources, int size) +// Add background into the fake image +// Based on ppSubBackground() +static psImage *stackBackgroundModel(pmReadout *ro, // Readout for which to generate background model + const pmConfig *config // Configuration + ) { - psImage *image = ro->image, *mask = ro->mask; // Image and mask of readout + psAssert(ro && ro->image, "Need readout image"); + psAssert(config, "Need configuration"); + + psImage *image = ro->image; // Image of interest int numCols = image->numCols, numRows = image->numRows; // Size of image - psImage *bgImage = psImageAlloc(numCols, numRows, PS_TYPE_F32); // Background image - psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN); // Statistics for measuring background - psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS, 0); // Random number generator - - for (int i = 0; i < sources->n; i++) { - pmSource *source = sources->data[i]; // Source of interest - if (!source) { - continue; - } - - int x = source->peak->xf + 0.5, y = source->peak->yf + 0.5; // Coordinates of source - - int xMin = PS_MAX(x - BG_SIZE, 0), xMax = PS_MIN(x + BG_SIZE, numCols); - int yMin = PS_MAX(y - BG_SIZE, 0), yMax = PS_MIN(y + BG_SIZE, numCols); - - psRegion region = psRegionSet(xMin, xMax, yMin, yMax); // Region of interest - psImage *subImage = psImageSubset(image, region); // Subimage - psImage *subMask = psImageSubset(mask, region); // Subimage mask - if (!psImageBackground(stats, NULL, subImage, subMask, maskVal, rng)) { - // Can't do anything about it - psErrorClear(); - continue; - } - psFree(subImage); - psFree(subMask); - - float value = stats->robustMedian; // Value to set background to - - int uMin = PS_MAX(x - size, 0), uMax = PS_MIN(x + size, numCols); - int vMin = PS_MAX(y - size, 0), vMax = PS_MIN(y + size, numCols); - - for (int v = vMin; v < vMax; v++) { - for (int u = uMin; u < uMax; u++) { - bgImage->data.F32[v][u] = value; - } - } - } - - psFree(stats); - psFree(rng); - - return bgImage; + psMetadata *ppStackRecipe = psMetadataLookupPtr(NULL, config->recipes, PPSTACK_RECIPE); + psAssert(ppStackRecipe, "Need PPSTACK recipe"); + psMetadata *psphotRecipe = psMetadataLookupPtr(NULL, config->recipes, PSPHOT_RECIPE); + psAssert(psphotRecipe, "Need PSPHOT recipe"); + + psString maskBadStr = psMetadataLookupStr(NULL, ppStackRecipe, "MASK.BAD");// Name of bits to mask for bad + psMaskType maskBad = pmConfigMaskGet(maskBadStr, config); // Bits to mask for bad pixels + + // user-defined masks to test for good/bad pixels (build from recipe list if not yet set) + psMetadataAddU8(psphotRecipe, PS_LIST_TAIL, "MASK.PSPHOT", PS_META_REPLACE, "user-defined mask", maskBad); + + psImage *binned = psphotBackgroundModel(ro, config); // Binned background model + psImageBinning *binning = psMetadataLookupPtr(NULL, psphotRecipe, + "PSPHOT.BACKGROUND.BINNING"); // Binning for model + psAssert(binning, "Need binning parameters"); + psImage *unbinned = psImageAlloc(numCols, numRows, PS_TYPE_F32); // Unbinned background model + if (!psImageUnbin(unbinned, binned, binning)) { + psError(PS_ERR_UNKNOWN, false, "Unable to unbin background model"); + psFree(unbinned); + return NULL; + } + + // XXX should these really be here?? (probably not...) + // pmFPAfileDropInternal(config->files, "PSPHOT.BACKMDL"); + // pmFPAfileDropInternal(config->files, "PSPHOT.BACKMDL.STDEV"); + + return unbinned; } + bool ppStackMatch(pmReadout *readout, psArray **regions, psArray **kernels, float *chi2, float *weighting, @@ -163,4 +172,6 @@ psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, PPSTACK_RECIPE); // ppStack recipe psAssert(recipe, "We've thrown an error on this before."); + + float deconvLimit = psMetadataLookupF32(NULL, recipe, "DECONV.LIMIT"); // Limit on deconvolution fraction // Look up appropriate values from the ppSub recipe @@ -195,31 +206,29 @@ assert(outName); // Read convolution kernel - { - psString filename = NULL; // Output filename - psStringAppend(&filename, "%s.%d.kernel", outName, numInput); - 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(output, fits)) { - psError(PS_ERR_IO, false, "Unable to read previously produced kernel"); - psFitsClose(fits); - return false; - } + psString filename = NULL; // Output filename + psStringAppend(&filename, "%s.%d.kernel", outName, numInput); + 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(output, fits)) { + psError(PS_ERR_IO, false, "Unable to read previously produced kernel"); psFitsClose(fits); - - // Add in variance factor - pmSubtractionKernels *kernels = psMetadataLookupPtr(NULL, output->analysis, - PM_SUBTRACTION_ANALYSIS_KERNEL); // Kernels - float vf = pmSubtractionVarianceFactor(kernels, 0.0, 0.0, false); // Variance factor - psMetadataItem *vfItem = psMetadataLookup(readout->parent->concepts, "CELL.VARFACTOR"); - if (!isfinite(vf)) { - vf = 1.0; - } - if (isfinite(vfItem->data.F32)) { - vfItem->data.F32 *= vf; - } else { - vfItem->data.F32 = vf; - } + return false; + } + psFitsClose(fits); + + // Add in variance factor + pmSubtractionKernels *kernels = psMetadataLookupPtr(NULL, output->analysis, + PM_SUBTRACTION_ANALYSIS_KERNEL); // Kernels + float vf = pmSubtractionVarianceFactor(kernels, 0.0, 0.0, false); // Variance factor + psMetadataItem *vfItem = psMetadataLookup(readout->parent->concepts, "CELL.VARFACTOR"); + if (!isfinite(vf)) { + vf = 1.0; + } + if (isfinite(vfItem->data.F32)) { + vfItem->data.F32 *= vf; + } else { + vfItem->data.F32 = vf; } @@ -244,4 +253,10 @@ psFree(maskName); psFree(weightName); + + psRegion *region = psMetadataLookupPtr(NULL, output->analysis, + PM_SUBTRACTION_ANALYSIS_REGION); // Convolution region + + pmSubtractionAnalysis(readout->analysis, kernels, region, + readout->image->numCols, readout->image->numRows); } else { #endif @@ -285,30 +300,4 @@ } - // Generate a fake image to match to - float minFlux = INFINITY; // Minimum flux for fake image - { - psStats *bg = psStatsAlloc(PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV); // Statistics for bg - if (!psImageBackground(bg, NULL, readout->image, readout->mask, maskVal, rng)) { - psWarning("Can't measure background for image."); - psErrorClear(); - } else { - minFlux = bg->robustStdev; - } - psFree(bg); - } - -#if 0 - float maxMag = -INFINITY; // Maximum magnitude of sources - for (int i = 0; i < sources->n; i++) { - pmSource *source = sources->data[i]; // Source of interest - if (source->psfMag > maxMag && source->psfMag <= MAG_IGNORE && - !(source->mode & SOURCE_MASK)) { - maxMag = source->psfMag; - } - } - minFlux = PS_MIN(FAINT_SOURCE_FRAC * powf(10.0, -0.4 * maxMag), minFlux); -#endif - - #if 0 // Testing the normalisation of the fake image @@ -318,7 +307,7 @@ pmSource *source = pmSourceAlloc(); // Source sources->data[0] = source; - source->peak = pmPeakAlloc(50, 50, 10000, PM_PEAK_LONE); + source->peak = pmPeakAlloc(500, 500, 10000, PM_PEAK_LONE); source->psfMag = -13.0; - pmReadoutFakeFromSources(test, 100, 100, sources, NULL, NULL, psf, 0.1, 0, false, true); + pmReadoutFakeFromSources(test, 1000, 1000, sources, NULL, NULL, psf, 0.1, 0, false, true); float sum = 0.0; for (int y = 0; y < test->image->numRows; y++) { @@ -342,5 +331,6 @@ if (!pmReadoutFakeFromSources(fake, readout->image->numCols, readout->image->numRows, - stampSources, NULL, NULL, psf, minFlux, 0, false, true)) { + stampSources, NULL, NULL, psf, NAN, footprint + size, + false, true)) { psError(PS_ERR_UNKNOWN, false, "Unable to generate fake image with target PSF."); psFree(fake); @@ -351,24 +341,25 @@ // Add the background into the target image - psImage *bgImage = stackBackgroundModel(readout, maskVal, stampSources, footprint); // Image of background + psImage *bgImage = stackBackgroundModel(readout, config); // Image of background psBinaryOp(fake->image, fake->image, "+", bgImage); psFree(bgImage); - #ifdef TESTING { + pmHDU *hdu = pmHDUFromCell(readout->parent); psString name = NULL; psStringAppend(&name, "fake_%03d.fits", numInput); psFits *fits = psFitsOpen(name, "w"); psFree(name); - psFitsWriteImage(fits, NULL, fake->image, 0, NULL); + psFitsWriteImage(fits, hdu->header, fake->image, 0, NULL); psFitsClose(fits); } { + pmHDU *hdu = pmHDUFromCell(readout->parent); psString name = NULL; psStringAppend(&name, "real_%03d.fits", numInput); psFits *fits = psFitsOpen(name, "w"); psFree(name); - psFitsWriteImage(fits, NULL, readout->image, 0, NULL); + psFitsWriteImage(fits, hdu->header, readout->image, 0, NULL); psFitsClose(fits); } @@ -410,12 +401,14 @@ #ifdef TESTING { + pmHDU *hdu = pmHDUFromCell(readout->parent); psString name = NULL; psStringAppend(&name, "conv_%03d.fits", numInput); psFits *fits = psFitsOpen(name, "w"); psFree(name); - psFitsWriteImage(fits, NULL, output->image, 0, NULL); + psFitsWriteImage(fits, hdu->header, output->image, 0, NULL); psFitsClose(fits); } { + pmHDU *hdu = pmHDUFromCell(readout->parent); psString name = NULL; psStringAppend(&name, "diff_%03d.fits", numInput); @@ -423,5 +416,5 @@ psFree(name); psBinaryOp(fake->image, output->image, "-", fake->image); - psFitsWriteImage(fits, NULL, fake->image, 0, NULL); + psFitsWriteImage(fits, hdu->header, fake->image, 0, NULL); psFitsClose(fits); } @@ -548,4 +541,14 @@ } + // Reject image completely if the maximum deconvolution fraction exceeds the limit + float deconv = psMetadataLookupF32(NULL, output->analysis, + PM_SUBTRACTION_ANALYSIS_DECONV_MAX); // Maximum deconvolution fraction + if (deconv > deconvLimit) { + psWarning("Maximum deconvolution fraction (%f) exceeds limit (%f) --- rejecting\n", + deconv, deconvLimit); + psFree(output); + return NULL; + } + // Renormalise the variances if desired if (renorm) { @@ -580,4 +583,74 @@ psFree(bg); + +#if 0 +#define RADIUS 10 // Radius of photometry +#define MIN_ERR 0.05 // Minimum photometric error, mag +#define MAX_MAG -13 // Maximum magnitude for source + + // Ensure the normalisation is correct + // XXX Ideally, would like to do proper PSF-fit photometry, but will settle for aperture photometry + int numSources = sources->n; // Number of sources + psVector *ratio = psVectorAlloc(numSources, PS_TYPE_F32); // Flux ratios for sources + psVector *ratioMask = psVectorAlloc(numSources, PS_TYPE_MASK); // Mask for flux ratios + psVectorInit(ratioMask, 0xFF); + psImage *image = readout->image; // Image of interest + psImage *mask = readout->mask; // Mask of interest + int numCols = image->numCols, numRows = image->numRows; // Size of image + for (int i = 0; i < numSources; i++) { + pmSource *source = sources->data[i]; // Source of interest + if (!source || source->mode & SOURCE_MASK || !isfinite(source->psfMag) || !isfinite(source->errMag) || + source->errMag > MIN_ERR || source->psfMag > MAX_MAG) { + continue; + } + + float xSrc, ySrc; // Source coordinates + coordsFromSource(&xSrc, &ySrc, source); + int xMin = PS_MAX(0, xSrc - RADIUS), xMax = PS_MIN(numCols - 1, xSrc + RADIUS); // Bounds in x + int yMin = PS_MAX(0, ySrc - RADIUS), yMax = PS_MIN(numRows - 1, ySrc + RADIUS); // Bounds in y + int numPix = 0; // Number of pixels + float sum = 0.0; // Sum of pixels + for (int y = yMin; y <= yMax; y++) { + for (int x = xMin; x <= xMax; x++) { + if (mask->data.PS_TYPE_MASK_DATA[y][x] & maskBad) { + continue; + } + sum += image->data.F32[y][x]; + numPix++; + } + } + if (sum >= 0 && numPix > 0) { + float mag = -2.5 * log10(sum * M_PI * PS_SQR(RADIUS) / numPix); // Instrumental magnitude + ratio->data.F32[i] = mag - source->psfMag; + ratioMask->data.PS_TYPE_MASK_DATA[i] = 0; + } + } + + psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV); // Statistics + if (!psVectorStats(stats, ratio, NULL, ratioMask, 0xFF)) { + psWarning("Unable to measure normalisation --- assuming correct."); + } else { + psLogMsg("ppStack", PS_LOG_INFO, "Renormalising image by %f (+/- %f) mag\n", + stats->robustMedian, stats->robustStdev); + float norm = powf(10.0, -0.4 * stats->robustMedian); // Normalisation to apply + psBinaryOp(image, image, "*", psScalarAlloc(norm, PS_TYPE_F32)); + } + psFree(stats); + psFree(ratio); + psFree(ratioMask); +#endif + +#ifdef TESTING + { + pmHDU *hdu = pmHDUFromCell(readout->parent); + psString name = NULL; + psStringAppend(&name, "convolved_%03d.fits", numInput); + psFits *fits = psFitsOpen(name, "w"); + psFree(name); + psFitsWriteImage(fits, hdu->header, output->image, 0, NULL); + psFitsClose(fits); + } +#endif + psFree(output);