Index: branches/tap_branches/ppSub/src/ppSubMatchPSFs.c =================================================================== --- branches/tap_branches/ppSub/src/ppSubMatchPSFs.c (revision 25900) +++ branches/tap_branches/ppSub/src/ppSubMatchPSFs.c (revision 27838) @@ -18,6 +18,9 @@ #include #include +#include #include "ppSub.h" + +#define COVAR_FRAC 0.01 // Truncation fraction for covariance matrix // Normalise a region on an image @@ -37,4 +40,148 @@ } +// Measure the PSF for an image +static float subImagePSF(ppSubData *data, // Processing data + const pmReadout *ro, // Readout for which to measure PSF + psArray *sources // Sources with positions at which to measure PSF + ) +{ + psAssert(data, "Require processing data"); + pmConfig *config = data->config; // Configuration + psAssert(config, "Require configuration"); + + psAssert(ro, "Need readout"); + psAssert(sources, "Need sources."); + + pmFPAfile *photFile = psMetadataLookupPtr(NULL, config->files, "PSPHOT.INPUT"); // Photometry file + psAssert(photFile, "Need photometry file."); + if (!pmFPACopy(photFile->fpa, ro->parent->parent->parent)) { + psError(PPSUB_ERR_CONFIG, false, "Unable to copy FPA for photometry"); + return false; + } + + pmFPAview *view = ppSubViewReadout(); // View to readout + pmReadout *photRO = pmFPAviewThisReadout(view, photFile->fpa); // Readout to photometer + + if (psMetadataLookup(photRO->analysis, "PSPHOT.DETECTIONS")) { + psMetadataRemoveKey(photRO->analysis, "PSPHOT.DETECTIONS"); + } + if (psMetadataLookup(photRO->parent->parent->analysis, "PSPHOT.PSF")) { + psMetadataRemoveKey(photRO->parent->parent->analysis, "PSPHOT.PSF"); + } + + // Extract the loaded sources from the associated readout, and generate PSF + // Here, we assume the image is background-subtracted + if (!psphotReadoutFindPSF(config, view, sources)) { + psErrorStackPrint(stderr, "Unable to determine PSF"); + psWarning("Unable to determine PSF."); + psFree(view); + return true; + } + + psFree(view); + + psMetadata *header = psMetadataLookupMetadata(NULL, photRO->analysis, "PSPHOT.HEADER"); + psAssert(header, "Require header."); + float fwhm = psMetadataLookupF32(NULL, header, "FWHM_MAJ"); + + if (psMetadataLookup(photRO->analysis, "PSPHOT.DETECTIONS")) { + psMetadataRemoveKey(photRO->analysis, "PSPHOT.DETECTIONS"); + } + if (psMetadataLookup(photRO->parent->parent->analysis, "PSPHOT.PSF")) { + psMetadataRemoveKey(photRO->parent->parent->analysis, "PSPHOT.PSF"); + } + + return fwhm; +} + +// Scale the kernel parameters according to the PSFs +static bool subScaleKernel(ppSubData *data, // Processing data + psVector *kernelWidths, // Widths for kernel + int *kernelSize, // Size of kernel + int *stampSize // Size of stamps (footprint) + ) +{ + psAssert(data, "Require processing data"); + pmConfig *config = data->config; // Configuration + psAssert(config, "Require configuration"); + + // Nothing to do if pre-calculated kernel exists + pmFPAview *view = ppSubViewReadout(); // View to readout + pmReadout *kernelRO = pmFPAfileThisReadout(config->files, view, "PPSUB.OUTPUT.KERNELS"); // RO with kernel + if (kernelRO) { + psFree(view); + return true; + } + + // Look up recipe values + psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, PPSUB_RECIPE); // Recipe for ppSim + psAssert(recipe, "We checked this earlier, so it should be here."); + if (!psMetadataLookupBool(NULL, recipe, "SCALE")) { + // No scaling requested + psFree(view); + return true; + } + + // Input images + pmReadout *inRO = pmFPAfileThisReadout(config->files, view, "PPSUB.INPUT"); // Input readout + pmReadout *refRO = pmFPAfileThisReadout(config->files, view, "PPSUB.REF"); // Reference readout + + // Input sources + pmReadout *inSourceRO = pmFPAfileThisReadout(config->files, view, "PPSUB.INPUT.SOURCES"); + pmReadout *refSourceRO = pmFPAfileThisReadout(config->files, view, "PPSUB.REF.SOURCES"); + if (!inSourceRO || !refSourceRO) { + psWarning("Unable to scale kernel, since no sources were provided."); + return true; + } + + pmDetections *inDetections = psMetadataLookupPtr(NULL, inSourceRO->analysis, "PSPHOT.DETECTIONS"); // Input sources + pmDetections *refDetections = psMetadataLookupPtr(NULL, refSourceRO->analysis, "PSPHOT.DETECTIONS"); // Ref sources + + psFree(view); + + if (!inDetections || !refDetections) { + psWarning("Unable to scale kernel, since no sources were provided."); + return true; + } + + psArray *inSources = inDetections->allSources; + psAssert (inSources, "missing inSources?"); + + psArray *refSources = refDetections->allSources; + psAssert (refSources, "missing refSources?"); + + float inFWHM = psMetadataLookupF32(NULL, inRO->parent->parent->concepts, "CHIP.SEEING"); // FWHM for input + if (!isfinite(inFWHM) || inFWHM == 0.0) { + inFWHM = subImagePSF(data, inRO, inSources); + } + float refFWHM = psMetadataLookupF32(NULL, refRO->parent->parent->concepts, "CHIP.SEEING"); // FWHM for ref + if (!isfinite(refFWHM) || refFWHM == 0.0) { + refFWHM = subImagePSF(data, refRO, refSources); + } + psLogMsg("ppSub", PS_LOG_INFO, "Input FWHM: %f\nReference FWHM: %f\n", inFWHM, refFWHM); + if (!isfinite(inFWHM) || !isfinite(refFWHM)) { + psWarning("Unable to scale kernel, since unable to measure PSFs."); + return true; + } + + float scaleRef = psMetadataLookupF32(NULL, recipe, "SCALE.REF"); // Reference for scaling + float scaleMin = psMetadataLookupF32(NULL, recipe, "SCALE.MIN"); // Minimum for scaling + float scaleMax = psMetadataLookupF32(NULL, recipe, "SCALE.MAX"); // Maximum for scaling + if (!isfinite(scaleRef) || !isfinite(scaleMin) || !isfinite(scaleMax)) { + psError(PPSUB_ERR_ARGUMENTS, false, + "Scale parameters (SCALE.REF=%f, SCALE.MIN=%f, SCALE.MAX=%f) not set in recipe.", + scaleRef, scaleMin, scaleMax); + return false; + } + + if (!pmSubtractionParamsScale(kernelSize, stampSize, kernelWidths, inFWHM, refFWHM, + scaleRef, scaleMin, scaleMax)) { + psError(PPSUB_ERR_DATA, false, "Unable to scale parameters."); + return false; + } + + return true; +} + bool ppSubMatchPSFs(ppSubData *data) @@ -71,37 +218,44 @@ pmReadout *kernelRO = pmFPAfileThisReadout(config->files, view, "PPSUB.OUTPUT.KERNELS"); // RO with kernel - psFree(view); - // Sources in image, used for stamps: these must be loaded from previous analysis stages pmReadout *inSourceRO = pmFPAfileThisReadout(config->files, view, "PPSUB.INPUT.SOURCES"); pmReadout *refSourceRO = pmFPAfileThisReadout(config->files, view, "PPSUB.REF.SOURCES"); - psArray *inSources = inSourceRO ? psMetadataLookupPtr(&mdok, inSourceRO->analysis, "PSPHOT.SOURCES") : - NULL; // Source list from input image - psArray *refSources = refSourceRO ? psMetadataLookupPtr(&mdok, refSourceRO->analysis, "PSPHOT.SOURCES") : - NULL ; // Source list from reference image - - psArray *sources = NULL; // Merged list of sources - if (inSources && refSources) { + + pmDetections *inDetections = inSourceRO ? psMetadataLookupPtr(&mdok, inSourceRO->analysis, "PSPHOT.DETECTIONS") : NULL; // Input sources + pmDetections *refDetections = refSourceRO ? psMetadataLookupPtr(&mdok, refSourceRO->analysis, "PSPHOT.DETECTIONS") : NULL; // Ref sources + + psFree(view); + + pmDetections *detections = NULL; // Merged detection set + if (inDetections && refDetections) { + psArray *inSources = inDetections->allSources; + psArray *refSources = refDetections->allSources; + + psAssert (inSources, "missing in sources?"); + psAssert (refSources, "missing ref sources?"); + + detections = pmDetectionsAlloc(); float radius = psMetadataLookupF32(NULL, recipe, "SOURCE.RADIUS"); // Matching radius psArray *lists = psArrayAlloc(2); // Source lists lists->data[0] = psMemIncrRefCounter(inSources); lists->data[1] = psMemIncrRefCounter(refSources); - sources = pmSourceMatchMerge(lists, radius); + detections->allSources = pmSourceMatchMerge(lists, radius); psFree(lists); - if (!sources) { - psError(PS_ERR_UNKNOWN, false, "Unable to merge source lists"); + if (!detections->allSources) { + psError(PPSUB_ERR_DATA, false, "Unable to merge source lists"); + psFree(detections); return false; } - } else if (inSources) { - sources = psMemIncrRefCounter(inSources); - } else if (refSources) { - sources = psMemIncrRefCounter(refSources); - } - - psMetadataAddArray(inConv->analysis, PS_LIST_TAIL, "PSPHOT.SOURCES", PS_META_REPLACE, - "Merged source list", sources); - psMetadataAddArray(refConv->analysis, PS_LIST_TAIL, "PSPHOT.SOURCES", PS_META_REPLACE, - "Merged source list", sources); - psFree(sources); // Drop reference + } + if (!detections && inDetections) { + detections = psMemIncrRefCounter(inDetections); + } + if (!detections && refDetections) { + detections = psMemIncrRefCounter(refDetections); + } + + psMetadataAddPtr(inConv->analysis, PS_LIST_TAIL, "PSPHOT.DETECTIONS", PS_META_REPLACE | PS_DATA_UNKNOWN, "Merged source list", detections); + psMetadataAddPtr(refConv->analysis, PS_LIST_TAIL, "PSPHOT.DETECTIONS", PS_META_REPLACE | PS_DATA_UNKNOWN, "Merged source list", detections); + psFree(detections); // Drop reference int footprint = psMetadataLookupS32(NULL, recipe, "STAMP.FOOTPRINT"); // Stamp half-size @@ -115,5 +269,5 @@ pmSubtractionKernelsType type = pmSubtractionKernelsTypeFromString(typeStr); // Type of kernel if (type == PM_SUBTRACTION_KERNEL_NONE) { - psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Unrecognised kernel type: %s", typeStr); + psError(PPSUB_ERR_ARGUMENTS, true, "Unrecognised kernel type: %s", typeStr); return false; } @@ -130,5 +284,9 @@ int iter = psMetadataLookupS32(NULL, recipe, "ITER"); // Rejection iterations float rej = psMetadataLookupF32(NULL, recipe, "REJ"); // Rejection threshold - float sys = psMetadataLookupF32(NULL, recipe, "SYS"); // Relative systematic error in kernel + float kernelErr = psMetadataLookupF32(NULL, recipe, "KERNEL.ERR"); // Relative systematic error in kernel + float normFrac = psMetadataLookupF32(NULL, recipe, "NORM.FRAC"); // Fraction of window for normalisn windw + float sysErr = psMetadataLookupF32(NULL, recipe, "SYS.ERR"); // Relative systematic error in images + float skyErr = psMetadataLookupF32(NULL, recipe, "SKY.ERR"); // Additional error in sky + float covarFrac = psMetadataLookupF32(NULL, recipe, "COVAR.FRAC"); // Fraction for covariance calculation float badFrac = psMetadataLookupF32(NULL, recipe, "BADFRAC"); // Maximum bad fraction @@ -168,13 +326,29 @@ break; default: - psErrorStackPrint(stderr, "Invalid value for -convolve"); + psError(PPSUB_ERR_ARGUMENTS, false, "Invalid value for -convolve"); return false; } } + if (!subScaleKernel(data, widths, &size, &footprint)) { + psError(PPSUB_ERR_DATA, false, "Unable to scale kernel parameters"); + return false; + } + int threads = psMetadataLookupS32(NULL, config->arguments, "-threads"); // Number of threads if (threads > 0) { - pmSubtractionThreadsInit(inRO, refRO); - } + pmSubtractionThreadsInit(); + } + + if (inRO->covariance) { + psImageCovarianceTruncate(inRO->covariance, COVAR_FRAC); + } + if (refRO->covariance) { + psImageCovarianceTruncate(refRO->covariance, COVAR_FRAC); + } + + // Data doesn't exist until it's created in pmSubtraction... + inConv->data_exists = inConv->parent->data_exists = inConv->parent->parent->data_exists = false; + refConv->data_exists = refConv->parent->data_exists = refConv->parent->parent->data_exists = false; // Match the PSFs @@ -182,15 +356,36 @@ if (kernelRO) { success = pmSubtractionMatchPrecalc(inConv, refConv, inRO, refRO, kernelRO->analysis, - stride, sys, maskVal, maskBad, maskPoor, poorFrac, badFrac); + stride, kernelErr, covarFrac, maskVal, maskBad, maskPoor, + poorFrac, badFrac); } else { success = pmSubtractionMatch(inConv, refConv, inRO, refRO, footprint, stride, regionSize, - spacing, threshold, sources, data->stamps, type, size, order, - widths, orders, inner, ringsOrder, binning, penalty, optimum, - optWidths, optOrder, optThresh, iter, rej, sys, maskVal, - maskBad, maskPoor, poorFrac, badFrac, subMode); - } + spacing, threshold, detections ? detections->allSources : NULL, + data->stamps, type, size, order, widths, orders, inner, ringsOrder, + binning, penalty, optimum, optWidths, optOrder, optThresh, iter, rej, + normFrac, sysErr, skyErr, kernelErr, covarFrac, + maskVal, maskBad, maskPoor, poorFrac, badFrac, subMode); + } + +# ifdef TESTING + // XXX for testing + psphotSaveImage (NULL, refRO->image, "refRO.im.fits"); + psphotSaveImage (NULL, refRO->variance, "refRO.wt.fits"); + psphotSaveImage (NULL, refRO->mask, "refRO.mk.fits"); + + psphotSaveImage (NULL, inRO->image, "inRO.im.fits"); + psphotSaveImage (NULL, inRO->variance, "inRO.wt.fits"); + psphotSaveImage (NULL, inRO->mask, "inRO.mk.fits"); + + psphotSaveImage (NULL, inConv->image, "inConv.im.fits"); + psphotSaveImage (NULL, inConv->variance, "inConv.wt.fits"); + psphotSaveImage (NULL, inConv->mask, "inConv.mk.fits"); + + psphotSaveImage (NULL, refConv->image, "refConv.im.fits"); + psphotSaveImage (NULL, refConv->variance, "refConv.wt.fits"); + psphotSaveImage (NULL, refConv->mask, "refConv.mk.fits"); +# endif psFree(optWidths); - pmSubtractionThreadsFinalize(inRO, refRO); + pmSubtractionThreadsFinalize(); if (!success) { @@ -207,5 +402,5 @@ return true; } else { - psError(PS_ERR_UNKNOWN, false, "Unable to match images."); + psError(PPSUB_ERR_DATA, false, "Unable to match images."); return false; } @@ -260,4 +455,11 @@ pmConceptsCopyFPA(refConv->parent->parent->parent, refRO->parent->parent->parent, true, true); + if (inConv->covariance) { + psImageCovarianceTruncate(inConv->covariance, COVAR_FRAC); + } + if (refConv->covariance) { + psImageCovarianceTruncate(refConv->covariance, COVAR_FRAC); + } + if (inConv->variance) { psImageCovarianceTransfer(inConv->variance, inConv->covariance);