Index: branches/tap_branches/ppStack/src/ppStackMatch.c =================================================================== --- branches/tap_branches/ppStack/src/ppStackMatch.c (revision 25900) +++ branches/tap_branches/ppStack/src/ppStackMatch.c (revision 27838) @@ -18,5 +18,5 @@ #define COVAR_FRAC 0.01 // Truncation fraction for covariance matrix -//#define TESTING // Enable debugging output +// #define TESTING // Enable debugging output #ifdef TESTING @@ -31,10 +31,10 @@ psFree(resolved); if (!fits) { - psError(PS_ERR_IO, false, "Unable to open previously produced image: %s", name); + 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(PS_ERR_IO, false, "Unable to read previously produced image: %s", name); + psError(PPSTACK_ERR_IO, false, "Unable to read previously produced image: %s", name); psFitsClose(fits); return false; @@ -115,4 +115,6 @@ psFree(coords); psFree(tree); + psFree(x); + psFree(y); psLogMsg("ppStack", PS_LOG_INFO, "Filtered out %d of %d sources", numFiltered, numGood); @@ -144,5 +146,5 @@ psMetadataAddU8(psphotRecipe, PS_LIST_TAIL, "MASK.PSPHOT", PS_META_REPLACE, "user-defined mask", maskBad); - psImage *binned = psphotBackgroundModel(ro, config); // Binned background model + psImage *binned = psphotModelBackgroundReadoutNoFile(ro, config); // Binned background model psImageBinning *binning = psMetadataLookupPtr(NULL, ro->analysis, "PSPHOT.BACKGROUND.BINNING"); // Binning for model @@ -150,14 +152,12 @@ 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"); + psError(PPSTACK_ERR_DATA, false, "Unable to unbin background model"); + psFree(binned); 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; + psFree(binned); + + return unbinned; } @@ -167,4 +167,5 @@ ) { +#if 1 bool mdok; // Status of metadata lookups @@ -176,15 +177,15 @@ int num = psMetadataLookupS32(&mdok, recipe, "RENORM.NUM"); if (!mdok) { - psError(PS_ERR_BAD_PARAMETER_VALUE, true, "RENORM.NUM is not set in the recipe"); + 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(PS_ERR_BAD_PARAMETER_VALUE, true, "RENORM.MIN is not set in the recipe"); + 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(PS_ERR_BAD_PARAMETER_VALUE, true, "RENORM.MAX is not set in the recipe"); + psError(PPSTACK_ERR_CONFIG, true, "RENORM.MAX is not set in the recipe"); return false; } @@ -192,5 +193,9 @@ psImageMaskType maskBad = pmConfigMaskGet("BLANK", config); // Bits to mask + psImageCovarianceTransfer(readout->variance, readout->covariance); return pmReadoutVarianceRenormalise(readout, maskBad, num, minValid, maxValid); +#else + return true; +#endif } @@ -212,5 +217,5 @@ int size = psMetadataLookupS32(NULL, ppsub, "KERNEL.SIZE"); // Kernel half-size - psString maskValStr = psMetadataLookupStr(NULL, recipe, "MASK.VAL"); // Name of bits to mask going in + psString maskValStr = psMetadataLookupStr(NULL, ppsub, "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 @@ -224,5 +229,5 @@ if (!pmReadoutMaskNonfinite(readout, maskVal)) { - psError(PS_ERR_UNKNOWN, false, "Unable to mask non-finite pixels in readout."); + psError(psErrorCodeLast(), false, "Unable to mask non-finite pixels in readout."); return false; } @@ -251,5 +256,5 @@ psFree(resolved); if (!fits || !pmReadoutReadSubtractionKernels(conv, fits)) { - psError(PS_ERR_IO, false, "Unable to read previously produced kernel"); + psError(PPSTACK_ERR_IO, false, "Unable to read previously produced kernel"); psFitsClose(fits); return false; @@ -257,8 +262,8 @@ psFitsClose(fits); - if (!readImage(&readout->image, options->imageNames->data[index], config) || - !readImage(&readout->mask, options->maskNames->data[index], config) || - !readImage(&readout->variance, options->varianceNames->data[index], config)) { - psError(PS_ERR_IO, false, "Unable to read previously produced image."); + 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; } @@ -269,9 +274,11 @@ PM_SUBTRACTION_ANALYSIS_KERNEL); - pmSubtractionAnalysis(readout->analysis, NULL, kernels, region, + 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; @@ -292,5 +299,10 @@ int iter = psMetadataLookupS32(NULL, ppsub, "ITER"); // Rejection iterations float rej = psMetadataLookupF32(NULL, ppsub, "REJ"); // Rejection threshold - float sysError = psMetadataLookupF32(NULL, ppsub, "SYS"); // Relative systematic error in kernel + 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 @@ -309,4 +321,16 @@ 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 @@ -322,5 +346,5 @@ psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS); // Random number generator if (!psImageBackground(bg, NULL, readout->image, readout->mask, maskVal | maskBad, rng)) { - psError(PS_ERR_UNKNOWN, false, "Can't measure background for image."); + psError(PPSTACK_ERR_DATA, false, "Can't measure background for image."); psFree(fake); psFree(optWidths); @@ -338,8 +362,9 @@ 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(PS_ERR_UNKNOWN, false, "Unable to generate fake image with target PSF."); + psError(PPSTACK_ERR_DATA, false, "Unable to generate fake image with target PSF."); psFree(fake); psFree(optWidths); @@ -347,11 +372,14 @@ 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 @@ -379,5 +407,5 @@ if (threads > 0) { - pmSubtractionThreadsInit(readout, fake); + pmSubtractionThreadsInit(); } @@ -387,7 +415,7 @@ if (kernel) { if (!pmSubtractionMatchPrecalc(NULL, conv, fake, readout, readout->analysis, - stride, sysError, maskVal, maskBad, maskPoor, + stride, kernelError, covarFrac, maskVal, maskBad, maskPoor, poorFrac, badFrac)) { - psError(PS_ERR_UNKNOWN, false, "Unable to convolve images."); + psError(psErrorCodeLast(), false, "Unable to convolve images."); psFree(fake); psFree(optWidths); @@ -395,26 +423,46 @@ psFree(conv); if (threads > 0) { - pmSubtractionThreadsFinalize(readout, fake); + pmSubtractionThreadsFinalize(); } return false; } } else { - if (!pmSubtractionMatch(NULL, conv, fake, readout, footprint, stride, regionSize, spacing, - threshold, stampSources, stampsName, type, size, order, widths, - orders, inner, ringsOrder, binning, penalty, - optimum, optWidths, optOrder, optThresh, iter, rej, sysError, - maskVal, maskBad, maskPoor, poorFrac, badFrac, - PM_SUBTRACTION_MODE_2)) { - psError(PS_ERR_UNKNOWN, false, "Unable to match images."); + // 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(readout, fake); + 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 @@ -447,5 +495,5 @@ if (threads > 0) { - pmSubtractionThreadsFinalize(readout, fake); + pmSubtractionThreadsFinalize(); } @@ -487,10 +535,5 @@ while ((item = psMetadataGetAndIncrement(iter))) { assert(item->type == PS_DATA_UNKNOWN); - // Set the normalisation dimensions, since these will be otherwise unavailable when reading - // the images by scans. pmSubtractionKernels *kernel = item->data.V; // Kernel used in subtraction - kernel->numCols = readout->image->numCols; - kernel->numRows = readout->image->numRows; - kernels = psArrayAdd(kernels, ARRAY_BUFFER, kernel); } @@ -518,10 +561,35 @@ } + // 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)); + } + // Reject image completely if the maximum deconvolution fraction exceeds the limit float deconv = psMetadataLookupF32(NULL, conv->analysis, PM_SUBTRACTION_ANALYSIS_DECONV_MAX); // Max deconvolution fraction if (deconv > deconvLimit) { - psWarning("Maximum deconvolution fraction (%f) exceeds limit (%f) --- rejecting\n", - deconv, deconvLimit); + psWarning("Maximum deconvolution fraction (%f) exceeds limit (%f) --- rejecting image %d\n", + deconv, deconvLimit, index); psFree(conv); return NULL; @@ -542,13 +610,13 @@ psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS); // Random number generator if (!psImageBackground(bg, NULL, readout->image, readout->mask, maskVal | maskBad, rng)) { - psWarning("Can't measure background for image."); - psErrorClear(); + psWarning("Can't measure background for image."); + psErrorClear(); } else { - if (!psMetadataLookupBool(NULL, config->arguments, "PPSTACK.SKIP.BG.SUB")) { - psLogMsg("ppStack", PS_LOG_INFO, "Correcting convolved image background by %lf (+/- %lf)", - psStatsGetValue(bg, PS_STAT_ROBUST_MEDIAN), psStatsGetValue(bg, PS_STAT_ROBUST_STDEV)); - (void)psBinaryOp(readout->image, readout->image, "-", - psScalarAlloc(psStatsGetValue(bg, PS_STAT_ROBUST_MEDIAN), PS_TYPE_F32)); - } + if (!psMetadataLookupBool(NULL, config->arguments, "PPSTACK.SKIP.BG.SUB")) { + psLogMsg("ppStack", PS_LOG_INFO, "Correcting convolved image background by %lf (+/- %lf)", + psStatsGetValue(bg, PS_STAT_ROBUST_MEDIAN), psStatsGetValue(bg, PS_STAT_ROBUST_STDEV)); + (void)psBinaryOp(readout->image, readout->image, "-", + psScalarAlloc(psStatsGetValue(bg, PS_STAT_ROBUST_MEDIAN), PS_TYPE_F32)); + } } @@ -560,14 +628,18 @@ // Measure the variance level for the weighting - if (!psImageBackground(bg, NULL, readout->variance, readout->mask, maskVal | maskBad, rng)) { - psError(PS_ERR_UNKNOWN, false, "Can't measure mean variance for image."); - psFree(rng); - psFree(bg); - return false; - } - options->weightings->data.F32[index] = 1.0 / (psStatsGetValue(bg, PS_STAT_ROBUST_MEDIAN) * - psImageCovarianceFactor(readout->covariance)); - psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "PPSTACK.WEIGHTING", 0, - "Weighting by 1/noise^2 for stack", options->weightings->data.F32[index]); + if (psMetadataLookupBool(NULL, recipe, "WEIGHTS")) { + if (!psImageBackground(bg, NULL, readout->variance, readout->mask, maskVal | maskBad, rng)) { + psError(PPSTACK_ERR_DATA, false, "Can't measure mean variance for image."); + psFree(rng); + psFree(bg); + return false; + } + options->weightings->data.F32[index] = 1.0 / (psStatsGetValue(bg, PS_STAT_ROBUST_MEDIAN) * + psImageCovarianceFactor(readout->covariance)); + } else { + options->weightings->data.F32[index] = 1.0; + } + psLogMsg("ppStack", PS_LOG_INFO, "Weighting for image %d is %f\n", + index, options->weightings->data.F32[index]); psFree(rng);