Index: /trunk/psphot/src/psphotStackMatchPSFs.c =================================================================== --- /trunk/psphot/src/psphotStackMatchPSFs.c (revision 27849) +++ /trunk/psphot/src/psphotStackMatchPSFs.c (revision 27850) @@ -5,8 +5,4 @@ bool status = true; - // select the appropriate recipe information - psMetadata *recipe = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE); - psAssert (recipe, "missing recipe?"); - int num = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.INPUT.NUM"); psAssert (status, "programming error: must define PSPHOT.INPUT.NUM"); @@ -14,5 +10,5 @@ // loop over the available readouts for (int i = 0; i < num; i++) { - if (!psphotStackMatchPSFsReadout (config, view, i, recipe)) { + if (!psphotStackMatchPSFsReadout (config, view, i)) { psError (PSPHOT_ERR_CONFIG, false, "failed to find initial detections for PSPHOT.INPUT entry %d", i); return false; @@ -23,5 +19,5 @@ // convolve the image to match desired PSF -bool psphotStackMatchPSFsReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe) { +bool psphotStackMatchPSFsReadout (pmConfig *config, const pmFPAview *view, int index) { bool status; @@ -29,7 +25,4 @@ float NSIGMA_PEAK = 25.0; int NMAX = 0; - - Raw; - Cnv; // find the currently selected readout @@ -50,223 +43,26 @@ psAssert (maskVal, "missing mask value?"); - psphotStackMatch(); + /***** set up recipe options *****/ - return true; -} - -#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 - -// #define TESTING // Enable debugging output - -#ifdef TESTING -// Read a FITS image -static 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(PPSTACK_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); - psFitsClose(fits); - return false; - } - psFitsClose(fits); - - psFree(*target); - *target = image; - - return true; -} -#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 - ) -{ - psAssert(sources && sources->n > 0, "Require array of sources"); - if (exclusion <= 0) { - return psMemIncrRefCounter(sources); - } - - int num = sources->n; // Number of sources - psVector *x = psVectorAlloc(num, PS_TYPE_F32), *y = psVectorAlloc(num, PS_TYPE_F32); // Coordinates - int numGood = 0; // Number of good sources - for (int i = 0; i < num; i++) { - pmSource *source = sources->data[i]; // Source of interest - if (!source) { - continue; - } - coordsFromSource(&x->data.F32[numGood], &y->data.F32[numGood], source); - numGood++; - } - x->n = y->n = numGood; - - psTree *tree = psTreePlant(2, 2, PS_TREE_EUCLIDEAN, x, y); // kd tree - - psArray *filtered = psArrayAllocEmpty(numGood); // Filtered list of sources - psVector *coords = psVectorAlloc(2, PS_TYPE_F64); // Coordinates of source - int numFiltered = 0; // Number of filtered sources - for (int i = 0; i < num; i++) { - pmSource *source = sources->data[i]; // Source of interest - if (!source) { - continue; - } - 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 - psTrace("ppStack", 9, "Source at %.0lf,%.0lf has %ld sources in exclusion zone", - coords->data.F64[0], coords->data.F64[1], numWithin); - if (numWithin == 1) { - // Only itself inside the exclusion zone - filtered = psArrayAdd(filtered, filtered->n, source); - } else { - numFiltered++; - } - } - psFree(coords); - psFree(tree); - psFree(x); - psFree(y); - - psLogMsg("ppStack", PS_LOG_INFO, "Filtered out %d of %d sources", numFiltered, numGood); - - return filtered; -} - -// 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 - ) -{ - 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 - - 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 = psphotModelBackgroundReadoutNoFile(ro, config); // Binned background model - psImageBinning *binning = psMetadataLookupPtr(NULL, ro->analysis, - "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(PPSTACK_ERR_DATA, false, "Unable to unbin background model"); - psFree(binned); - psFree(unbinned); - return NULL; - } - psFree(binned); - - return unbinned; -} - -// Renormalise a readout's variance map -bool stackRenormaliseReadout(const pmConfig *config, // Configuration - pmReadout *readout // Readout to renormalise - ) -{ -#if 1 - bool mdok; // Status of metadata lookups - - psMetadata *recipe = psMetadataLookupPtr(NULL, config->recipes, PPSTACK_RECIPE); // Recipe for ppStack - psAssert(recipe, "Need PPSTACK recipe"); - - if (!psMetadataLookupBool(&mdok, recipe, "RENORM")) return true; - - int num = psMetadataLookupS32(&mdok, recipe, "RENORM.NUM"); - if (!mdok) { - psError(PPSTACK_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"); - return false; - } - float maxValid = psMetadataLookupF32(&mdok, recipe, "RENORM.MAX"); - if (!mdok) { - psError(PPSTACK_ERR_CONFIG, true, "RENORM.MAX is not set in the recipe"); - return false; - } - - psImageMaskType maskBad = pmConfigMaskGet("BLANK", config); // Bits to mask - - psImageCovarianceTransfer(readout->variance, readout->covariance); - return pmReadoutVarianceRenormalise(readout, maskBad, num, minValid, maxValid); -#else - return true; -#endif -} - -bool ppStackMatch(pmReadout *readout, ppStackOptions *options, int index, const pmConfig *config) -{ - assert(readout); - assert(options); - assert(config); - - 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 + 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 *ppsub = psMetadataLookupMetadata(NULL, config->recipes, "PPSUB"); // PPSUB recipe - int size = psMetadataLookupS32(NULL, ppsub, "KERNEL.SIZE"); // Kernel half-size + psMetadata *subRecipe = psMetadataLookupMetadata(NULL, config->recipes, "PPSUB"); // PPSUB recipe + psAssert(subRecipe, "recipe missing"); - psString maskValStr = psMetadataLookupStr(NULL, ppsub, "MASK.VAL"); // Name of bits to mask going in + 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, recipe, "MASK.POOR"); // Name of bits to mask for poor + 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, recipe, "MASK.BAD"); // Name of bits to mask for bad + 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, ppsub, "PENALTY"); // Penalty for wideness + float penalty = psMetadataLookupF32(NULL, subRecipe, "PENALTY"); // Penalty for wideness int threads = psMetadataLookupS32(NULL, config->arguments, "-threads"); // Number of threads @@ -276,363 +72,26 @@ } - // Match the PSF + // Image Matching (PSFs or just flux) if (options->convolve) { + // Full match of PSFs pmReadout *conv = pmReadoutAlloc(NULL); // Conv readout, for holding results temporarily -#ifdef TESTING - // 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. // Read previously produced kernel if (psMetadataLookupBool(NULL, config->arguments, "PPSTACK.DEBUG.STACK")) { - pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PPSTACK.CONV.KERNEL", index); - psAssert(file, "Require file"); + loadKernel(); + } else { + matchKernel(); + } // !DEBUG.STACK - 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 + saveMatchData(); - // 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); + saveChiSquare(); - 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); - } else { -#endif - - // Normal operations here - psAssert(options->psf, "Require target PSF"); - psAssert(options->sourceLists && options->sourceLists->data[index], "Require source list"); - - int order = psMetadataLookupS32(NULL, ppsub, "SPATIAL.ORDER"); // Spatial polynomial order - float regionSize = psMetadataLookupF32(NULL, ppsub, "REGION.SIZE"); // Size of iso-kernel regs - float spacing = psMetadataLookupF32(NULL, ppsub, "STAMP.SPACING"); // Typical stamp spacing - int footprint = psMetadataLookupS32(NULL, ppsub, "STAMP.FOOTPRINT"); // Stamp half-size - float threshold = psMetadataLookupF32(NULL, ppsub, "STAMP.THRESHOLD"); // Threshold for stmps - int stride = psMetadataLookupS32(NULL, ppsub, "STRIDE"); // Size of convolution patches - int iter = psMetadataLookupS32(NULL, ppsub, "ITER"); // Rejection iterations - float rej = psMetadataLookupF32(NULL, ppsub, "REJ"); // Rejection threshold - float kernelError = psMetadataLookupF32(NULL, ppsub, "KERNEL.ERR"); // Relative systematic error in kernel - float normFrac = psMetadataLookupF32(NULL, ppsub, "NORM.FRAC"); // Fraction of window for normalisn windw - float sysError = psMetadataLookupF32(NULL, ppsub, "SYS.ERR"); // Relative systematic error in images - float skyErr = psMetadataLookupF32(NULL, ppsub, "SKY.ERR"); // Additional error in sky - float covarFrac = psMetadataLookupF32(NULL, ppsub, "COVAR.FRAC"); // Fraction for covariance calculation - - const char *typeStr = psMetadataLookupStr(NULL, ppsub, "KERNEL.TYPE"); // Kernel type - pmSubtractionKernelsType type = pmSubtractionKernelsTypeFromString(typeStr); // Kernel type - psVector *widths = psMetadataLookupPtr(NULL, ppsub, "ISIS.WIDTHS"); // ISIS Gaussian widths - psVector *orders = psMetadataLookupPtr(NULL, ppsub, "ISIS.ORDERS"); // ISIS Polynomial orders - int inner = psMetadataLookupS32(NULL, ppsub, "INNER"); // Inner radius - int ringsOrder = psMetadataLookupS32(NULL, ppsub, "RINGS.ORDER"); // RINGS polynomial order - int binning = psMetadataLookupS32(NULL, ppsub, "SPAM.BINNING"); // Binning for SPAM kernel - float badFrac = psMetadataLookupF32(NULL, ppsub, "BADFRAC"); // Maximum bad fraction - bool optimum = psMetadataLookupBool(&mdok, ppsub, "OPTIMUM"); // Derive optimum parameters? - float optMin = psMetadataLookupF32(&mdok, ppsub, "OPTIMUM.MIN"); // Minimum width for search - float optMax = psMetadataLookupF32(&mdok, ppsub, "OPTIMUM.MAX"); // Maximum width for search - float optStep = psMetadataLookupF32(&mdok, ppsub, "OPTIMUM.STEP"); // Step for search - float optThresh = psMetadataLookupF32(&mdok, ppsub, "OPTIMUM.TOL"); // Tolerance for search - int optOrder = psMetadataLookupS32(&mdok, ppsub, "OPTIMUM.ORDER"); // Order for search - float poorFrac = psMetadataLookupF32(&mdok, ppsub, "POOR.FRACTION"); // Fraction for "poor" - - bool scale = psMetadataLookupBool(NULL, ppsub, "SCALE"); // Scale kernel parameters? - float scaleRef = psMetadataLookupF32(NULL, ppsub, "SCALE.REF"); // Reference for scaling - float scaleMin = psMetadataLookupF32(NULL, ppsub, "SCALE.MIN"); // Minimum for scaling - float scaleMax = psMetadataLookupF32(NULL, ppsub, "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 - 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); - -#if 1 - // Add the background into the target image - psImage *bgImage = stackBackgroundModel(readout, config); // Image of background - psBinaryOp(fake->image, fake->image, "+", bgImage); - psFree(bgImage); -#endif - -#ifdef TESTING - { - 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); - } -#endif - - 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); - } - - -#ifdef TESTING - { - 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); - } -#endif - - 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); -#ifdef TESTING - } -#endif - - // 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); - } - 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"); - - // 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); - } - - // 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)); - } + 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 + 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); + psWarning("Maximum deconvolution fraction (%f) exceeds limit (%f) --- rejecting image %d\n", deconv, deconvLimit, index); psFree(conv); return NULL; @@ -643,5 +102,5 @@ psFree(conv); } else { - // Match the normalisation + // 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)); @@ -671,5 +130,5 @@ // Measure the variance level for the weighting - if (psMetadataLookupBool(NULL, recipe, "WEIGHTS")) { + 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."); @@ -678,6 +137,5 @@ return false; } - options->weightings->data.F32[index] = 1.0 / (psStatsGetValue(bg, PS_STAT_ROBUST_MEDIAN) * - psImageCovarianceFactor(readout->covariance)); + options->weightings->data.F32[index] = 1.0 / (psStatsGetValue(bg, PS_STAT_ROBUST_MEDIAN) * psImageCovarianceFactor(readout->covariance)); } else { options->weightings->data.F32[index] = 1.0; @@ -689,16 +147,5 @@ psFree(bg); -#ifdef TESTING - { - 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); - } -#endif + dumpImage3(); return true; Index: /trunk/psphot/src/psphotStackMatchPSFsUtils.c =================================================================== --- /trunk/psphot/src/psphotStackMatchPSFsUtils.c (revision 27850) +++ /trunk/psphot/src/psphotStackMatchPSFsUtils.c (revision 27850) @@ -0,0 +1,556 @@ +/***** 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 + +// 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(PPSTACK_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); + psFitsClose(fits); + return false; + } + psFitsClose(fits); + + psFree(*target); + *target = image; + + return true; +} + +// Get coordinates from a source +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; +} + +psArray *stackSourcesFilter(psArray *sources, // Source list to filter + int exclusion // Exclusion zone, pixels + ) +{ + psAssert(sources && sources->n > 0, "Require array of sources"); + if (exclusion <= 0) { + return psMemIncrRefCounter(sources); + } + + int num = sources->n; // Number of sources + psVector *x = psVectorAlloc(num, PS_TYPE_F32), *y = psVectorAlloc(num, PS_TYPE_F32); // Coordinates + int numGood = 0; // Number of good sources + for (int i = 0; i < num; i++) { + pmSource *source = sources->data[i]; // Source of interest + if (!source) { + continue; + } + coordsFromSource(&x->data.F32[numGood], &y->data.F32[numGood], source); + numGood++; + } + x->n = y->n = numGood; + + psTree *tree = psTreePlant(2, 2, PS_TREE_EUCLIDEAN, x, y); // kd tree + + psArray *filtered = psArrayAllocEmpty(numGood); // Filtered list of sources + psVector *coords = psVectorAlloc(2, PS_TYPE_F64); // Coordinates of source + int numFiltered = 0; // Number of filtered sources + for (int i = 0; i < num; i++) { + pmSource *source = sources->data[i]; // Source of interest + if (!source) { + continue; + } + 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 + psTrace("ppStack", 9, "Source at %.0lf,%.0lf has %ld sources in exclusion zone", + coords->data.F64[0], coords->data.F64[1], numWithin); + if (numWithin == 1) { + // Only itself inside the exclusion zone + filtered = psArrayAdd(filtered, filtered->n, source); + } else { + numFiltered++; + } + } + psFree(coords); + psFree(tree); + psFree(x); + psFree(y); + + psLogMsg("ppStack", PS_LOG_INFO, "Filtered out %d of %d sources", numFiltered, numGood); + + return filtered; +} + +// 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 + ) +{ + 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 + + 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 = psphotModelBackgroundReadoutNoFile(ro, config); // Binned background model + psImageBinning *binning = psMetadataLookupPtr(NULL, ro->analysis, + "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(PPSTACK_ERR_DATA, false, "Unable to unbin background model"); + psFree(binned); + psFree(unbinned); + return NULL; + } + psFree(binned); + + return unbinned; +} + +// Renormalise a readout's variance map +bool stackRenormaliseReadout(const pmConfig *config, // Configuration + pmReadout *readout // Readout to renormalise + ) +{ +#if 1 + bool mdok; // Status of metadata lookups + + psMetadata *recipe = psMetadataLookupPtr(NULL, config->recipes, PPSTACK_RECIPE); // Recipe for ppStack + psAssert(recipe, "Need PPSTACK recipe"); + + if (!psMetadataLookupBool(&mdok, recipe, "RENORM")) return true; + + int num = psMetadataLookupS32(&mdok, recipe, "RENORM.NUM"); + if (!mdok) { + psError(PPSTACK_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"); + return false; + } + float maxValid = psMetadataLookupF32(&mdok, recipe, "RENORM.MAX"); + if (!mdok) { + psError(PPSTACK_ERR_CONFIG, true, "RENORM.MAX is not set in the recipe"); + return false; + } + + psImageMaskType maskBad = pmConfigMaskGet("BLANK", config); // Bits to mask + + psImageCovarianceTransfer(readout->variance, readout->covariance); + return pmReadoutVarianceRenormalise(readout, maskBad, num, minValid, maxValid); +#else + return true; +#endif +} + +// 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. + +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() + { + 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 + 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); + } + 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)); + } + +}