Index: trunk/ppSub/src/ppSubReadout.c =================================================================== --- trunk/ppSub/src/ppSubReadout.c (revision 21183) +++ trunk/ppSub/src/ppSubReadout.c (revision 21257) @@ -1,723 +1,53 @@ -#ifdef HAVE_CONFIG_H -#include -#endif - -#include -#include -#include -#include -#include - #include "ppSub.h" - -#define WCS_TOLERANCE 0.001 // Tolerance for WCS -//#define TESTING // For test output - -#define SOURCE_MASK (PM_SOURCE_MODE_FAIL | PM_SOURCE_MODE_SATSTAR | PM_SOURCE_MODE_BLEND | \ - PM_SOURCE_MODE_BADPSF | PM_SOURCE_MODE_DEFECT | PM_SOURCE_MODE_SATURATED | \ - PM_SOURCE_MODE_CR_LIMIT | PM_SOURCE_MODE_EXT_LIMIT) // Mask to apply to input sources - - -// Copy every instance of a single keyword from one metadata to another -static bool metadataCopySingle(psMetadata *target, psMetadata *source, const char *name) -{ - PS_ASSERT_METADATA_NON_NULL(target, false); - PS_ASSERT_METADATA_NON_NULL(source, false); - PS_ASSERT_STRING_NON_EMPTY(name, false); - - psString regex = NULL; // Regular expression - psStringAppend(®ex, "^%s$", name); - psMetadataIterator *iter = psMetadataIteratorAlloc(source, PS_LIST_HEAD, regex); // Iterator - psFree(regex); - psMetadataItem *item; // Item from iteration - while ((item = psMetadataGetAndIncrement(iter))) { - psMetadataAddItem(target, item, PS_LIST_TAIL, PS_META_DUPLICATE_OK); - } - psFree(iter); - - return true; -} - bool ppSubReadout(pmConfig *config, psMetadata *stats, const pmFPAview *view) { psTimerStart("PPSUB_MATCH"); - pmReadout *inRO = pmFPAfileThisReadout(config->files, view, "PPSUB.INPUT"); // Input readout - pmReadout *refRO = pmFPAfileThisReadout(config->files, view, "PPSUB.REF"); // Reference readout - pmReadout *sourcesRO = pmFPAfileThisReadout(config->files, view, "PPSUB.SOURCES"); // Readout with sources - pmReadout *inConv = pmReadoutAlloc(NULL); // Convolved version of input - pmReadout *refConv = pmReadoutAlloc(NULL); // Convolved version of reference - pmCell *outCell = pmFPAfileThisCell(config->files, view, "PPSUB.OUTPUT"); // Output cell - pmReadout *outRO = pmReadoutAlloc(outCell); // Output readout: subtraction - pmFPA *outFPA = outCell->parent->parent; // Output FPA - pmHDU *outHDU = outFPA->hdu; // Output HDU - if (!outHDU->header) { - outHDU->header = psMetadataAlloc(); - } - psImage *input = inRO->image; // Input image - psImage *reference = refRO->image; // Reference image - PS_ASSERT_IMAGES_SIZE_EQUAL(input, reference, false); - - int threads = psMetadataLookupS32(NULL, config->arguments, "-threads"); // Number of threads - if (threads > 0) { - if (!psThreadPoolInit(threads)) { - psError(PS_ERR_UNKNOWN, false, "Unable to setup %d threads", threads); - return false; - } - pmSubtractionThreadsInit(inRO, refRO); - } - - // 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."); - - const char *typeStr = psMetadataLookupStr(NULL, recipe, "KERNEL.TYPE"); // Kernel type - psAssert(typeStr, "We put it here in ppSubArguments.c"); - 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); - psFree(inConv); - psFree(refConv); - psFree(outRO); + if (!ppSubSetMasks (config, view)) { + psError(PS_ERR_UNKNOWN, false, "Unable to subtract background."); return false; } - bool mdok; // Status of MD lookup - int size = psMetadataLookupS32(NULL, recipe, "KERNEL.SIZE"); // Kernel half-size - int order = psMetadataLookupS32(NULL, recipe, "SPATIAL.ORDER"); // Spatial polynomial order - float regionSize = psMetadataLookupF32(NULL, recipe, "REGION.SIZE"); // Size of iso-kernel regs - float spacing = psMetadataLookupF32(NULL, recipe, "STAMP.SPACING"); // Typical stamp spacing - int footprint = psMetadataLookupS32(NULL, recipe, "STAMP.FOOTPRINT"); // Stamp half-size - float threshold = psMetadataLookupF32(NULL, recipe, "STAMP.THRESHOLD"); // Threshold for stmps - int stride = psMetadataLookupS32(NULL, recipe, "STRIDE"); // Size of convolution patches - 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 - bool reverse = psMetadataLookupBool(NULL, config->arguments, "REVERSE"); // Reverse sense of subtraction? - psVector *widths = psMetadataLookupPtr(NULL, recipe, "ISIS.WIDTHS"); // ISIS Gaussian widths - psVector *orders = psMetadataLookupPtr(NULL, recipe, "ISIS.ORDERS"); // ISIS Polynomial orders - int inner = psMetadataLookupS32(NULL, recipe, "INNER"); // Inner radius - int ringsOrder = psMetadataLookupS32(NULL, recipe, "RINGS.ORDER"); // RINGS polynomial order - int binning = psMetadataLookupS32(NULL, recipe, "SPAM.BINNING"); // Binning for SPAM kernel - float penalty = psMetadataLookupF32(NULL, recipe, "PENALTY"); // Penalty for wideness - psString maskValStr = psMetadataLookupStr(NULL, recipe, "MASK.IN"); // 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 - 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 - psImageMaskType maskBad = pmConfigMaskGet(maskBadStr, config); // Bits to mask for bad pixels - float badFrac = psMetadataLookupF32(NULL, recipe, "BADFRAC"); // Maximum bad fraction - const char *stampsName = psMetadataLookupStr(&mdok, config->arguments, "STAMPS"); // Filename for stamps - - bool optimum = psMetadataLookupBool(&mdok, recipe, "OPTIMUM"); // Derive optimum parameters? - float optMin = psMetadataLookupF32(&mdok, recipe, "OPTIMUM.MIN"); // Minimum width for search - float optMax = psMetadataLookupF32(&mdok, recipe, "OPTIMUM.MAX"); // Maximum width for search - float optStep = psMetadataLookupF32(&mdok, recipe, "OPTIMUM.STEP"); // Step for search - float optThresh = psMetadataLookupF32(&mdok, recipe, "OPTIMUM.TOL"); // Tolerance for search - int optOrder = psMetadataLookupS32(&mdok, recipe, "OPTIMUM.ORDER"); // Order for search - bool dual = psMetadataLookupBool(&mdok, recipe, "DUAL"); // Dual convolution? - - // Statistics for renormalisation - psStatsOptions renormMean = psStatsOptionFromString(psMetadataLookupStr(NULL, recipe, "RENORM.MEAN")); - psStatsOptions renormStdev = psStatsOptionFromString(psMetadataLookupStr(NULL, recipe, "RENORM.STDEV")); - if (renormMean == PS_STAT_NONE || renormStdev == PS_STAT_NONE) { - psError(PS_ERR_BAD_PARAMETER_VALUE, false, - "Unable to parse renormalisation statistics from recipe."); - return false; - } - int renormNum = psMetadataLookupS32(&mdok, recipe, "RENORM.NUM"); // Number of samples for renormalisation - float renormWidth = psMetadataLookupS32(&mdok, recipe, "RENORM.WIDTH"); // Width of Gaussian phot - - psString interpModeStr = psMetadataLookupStr(&mdok, recipe, "INTERPOLATION"); // Interpolation mode - psImageInterpolateMode interpMode = psImageInterpolateModeFromString(interpModeStr); // Interp - if (interpMode == PS_INTERPOLATE_NONE) { - psError(PS_ERR_BAD_PARAMETER_VALUE, false, "Unknown interpolation mode: %s", interpModeStr); - return false; - } - float poorFrac = psMetadataLookupF32(&mdok, recipe, "POOR.FRACTION"); // Fraction for "poor" - - pmSubtractionMode subMode = dual ? PM_SUBTRACTION_MODE_DUAL : PM_SUBTRACTION_MODE_UNSURE; // Subtracn mode - - // Generate masks if they don't exist - int numCols = input->numCols, numRows = input->numRows; // Image dimensions - if (!inRO->mask) { - if (psMetadataLookupBool(NULL, recipe, "MASK.GENERATE")) { - pmReadoutSetMask(inRO, pmConfigMaskGet("SAT", config), pmConfigMaskGet("BAD", config)); - } else { - inRO->mask = psImageAlloc(numCols, numRows, PS_TYPE_IMAGE_MASK); - psImageInit(inRO->mask, 0); - } - } - if (!refRO->mask) { - if (psMetadataLookupBool(NULL, recipe, "MASK.GENERATE")) { - pmReadoutSetMask(refRO, pmConfigMaskGet("SAT", config), pmConfigMaskGet("BAD", config)); - } else { - refRO->mask = psImageAlloc(numCols, numRows, PS_TYPE_IMAGE_MASK); - psImageInit(refRO->mask, 0); - } - } - - if (!pmReadoutMaskNonfinite(inRO, pmConfigMaskGet("SAT", config))) { - psError(PS_ERR_UNKNOWN, false, "Unable to mask non-finite pixels in input."); - return false; - } - if (!pmReadoutMaskNonfinite(refRO, pmConfigMaskGet("SAT", config))) { - psError(PS_ERR_UNKNOWN, false, "Unable to mask non-finite pixels in reference."); + if (!ppSubMatchPSFs (config, view)) { + psError(PS_ERR_UNKNOWN, false, "Unable to subtract background."); return false; } - psVector *optWidths = NULL; // Vector with FWHMs for optimum search - if (optimum) { - optWidths = psVectorCreate(optWidths, optMin, optMax, optStep, PS_TYPE_F32); - } - - psArray *sources = NULL; // Sources in image - if (sourcesRO) { - sources = psMetadataLookupPtr(&mdok, sourcesRO->analysis, "PSPHOT.SOURCES"); - } - -#if 0 - // Interpolate over bad pixels, so the bad pixels don't explode - if (!pmReadoutInterpolateBadPixels(inRO, maskVal, interpMode, poorFrac, maskPoor, maskBad)) { - psError(PS_ERR_UNKNOWN, false, "Unable to interpolate bad pixels for input image."); - return false; - } - if (!pmReadoutInterpolateBadPixels(refRO, maskVal, interpMode, poorFrac, maskPoor, maskBad)) { - psError(PS_ERR_UNKNOWN, false, "Unable to interpolate bad pixels for reference image."); - return false; - } - maskVal |= maskBad; -#endif - - // Match the PSFs - if (!pmSubtractionMatch(inConv, refConv, inRO, refRO, footprint, stride, regionSize, spacing, threshold, - sources, stampsName, type, size, order, widths, orders, inner, ringsOrder, - binning, penalty, optimum, optWidths, optOrder, optThresh, iter, rej, sys, - maskVal, maskBad, maskPoor, poorFrac, badFrac, subMode)) { - psError(PS_ERR_UNKNOWN, false, "Unable to match images."); - psFree(inConv); - psFree(refConv); - psFree(outRO); - return false; - } - psFree(optWidths); - - pmSubtractionThreadsFinalize(inRO, refRO); - - // Add kernel descrption to header - pmSubtractionKernels *kernels = psMetadataLookupPtr(&mdok, inConv->analysis, - PM_SUBTRACTION_ANALYSIS_KERNEL); // Subtraction kernel - if (!kernels) { - kernels = psMetadataLookupPtr(&mdok, refConv->analysis, PM_SUBTRACTION_ANALYSIS_KERNEL); - } - if (!kernels) { - psError(PS_ERR_UNEXPECTED_NULL, true, "Unable to find SUBTRACTION.KERNEL"); - psFree(inConv); - psFree(refConv); - psFree(outRO); - return false; - } - psMetadataAddStr(outHDU->header, PS_LIST_TAIL, "PPSUB.KERNEL", 0, - "Subtraction kernel", kernels->description); - - { - if (psMetadataLookup(inConv->analysis, PM_SUBTRACTION_ANALYSIS_KERNEL)) { - outRO->analysis = psMetadataCopy(outRO->analysis, inConv->analysis); - } else if (psMetadataLookup(refConv->analysis, PM_SUBTRACTION_ANALYSIS_KERNEL)) { - outRO->analysis = psMetadataCopy(outRO->analysis, refConv->analysis); - } else { - psWarning("Unable to find subtraction kernel in analysis metadata."); - } - - // Update variance factors - // It's not possible to do this perfectly, since we're combining different images: - // vf_sum sigma_sum^2 = vf_1 * sigma_1^2 + vf_2 * sigma_2^2 - // It's not possible to write vf_sum as a function of vf_1 and vf_2 with no reference to the sigmas. - // Instead, we're going to cheat. - bool mdok; // Status of MD lookup - pmSubtractionKernels *kernels = psMetadataLookupPtr(&mdok, outRO->analysis, - PM_SUBTRACTION_ANALYSIS_KERNEL); // Kernels - if (!mdok) { - psError(PS_ERR_UNEXPECTED_NULL, true, "Unable to find subtraction kernels."); - psFree(inConv); - psFree(refConv); - psFree(outRO); - return false; - } - float vfIn = psMetadataLookupF32(NULL, inRO->parent->concepts, - "CELL.VARFACTOR"); // Variance factor for input - if (!isfinite(vfIn)) { - vfIn = 1.0; - } - float vfRef = psMetadataLookupF32(NULL, refRO->parent->concepts, - "CELL.VARFACTOR"); // Variance factor for ref - if (!isfinite(vfRef)) { - vfRef = 1.0; - } - if (kernels->mode == PM_SUBTRACTION_MODE_1 || kernels->mode == PM_SUBTRACTION_MODE_DUAL) { - vfIn *= pmSubtractionVarianceFactor(kernels, 0.0, 0.0, false); - } - if (kernels->mode == PM_SUBTRACTION_MODE_2) { - vfRef *= pmSubtractionVarianceFactor(kernels, 0.0, 0.0, false); - } else if (kernels->mode == PM_SUBTRACTION_MODE_DUAL) { - vfRef *= pmSubtractionVarianceFactor(kernels, 0.0, 0.0, true); - } - - psStats *vfStats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN); // Statistics - psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS, 0); // Random number generator - if (!psImageBackground(vfStats, NULL, inConv->weight, inConv->mask, maskVal | maskBad, rng)) { - psError(PS_ERR_UNKNOWN, false, "Unable to measure mean variance for convolved input"); - psFree(inConv); - psFree(refConv); - psFree(outRO); - psFree(vfStats); - psFree(rng); - return false; - } - float inMeanVar = vfStats->robustMedian; // Mean variance of input - if (!psImageBackground(vfStats, NULL, refConv->weight, refConv->mask, maskVal | maskBad, rng)) { - psError(PS_ERR_UNKNOWN, false, "Unable to measure mean variance for convolved reference"); - psFree(inConv); - psFree(refConv); - psFree(outRO); - psFree(vfStats); - psFree(rng); - return false; - } - float refMeanVar = vfStats->robustMedian; // Mean variance of reference - psFree(rng); - psFree(vfStats); - - float vf = (vfIn * inMeanVar + vfRef * refMeanVar) / (inMeanVar + refMeanVar); // Resulting var factor - psMetadataItem *vfItem = psMetadataLookup(outRO->parent->concepts, "CELL.VARFACTOR"); - vfItem->data.F32 = vf; - - // Statistics on the matching - if (stats) { - metadataCopySingle(stats, outRO->analysis, PM_SUBTRACTION_ANALYSIS_MODE); - metadataCopySingle(stats, outRO->analysis, PM_SUBTRACTION_ANALYSIS_STAMPS); - metadataCopySingle(stats, outRO->analysis, PM_SUBTRACTION_ANALYSIS_DEV_MEAN); - metadataCopySingle(stats, outRO->analysis, PM_SUBTRACTION_ANALYSIS_DEV_RMS); - metadataCopySingle(stats, outRO->analysis, PM_SUBTRACTION_ANALYSIS_NORM); - metadataCopySingle(stats, outRO->analysis, PM_SUBTRACTION_ANALYSIS_BGDIFF); - metadataCopySingle(stats, outRO->analysis, PM_SUBTRACTION_ANALYSIS_MX); - metadataCopySingle(stats, outRO->analysis, PM_SUBTRACTION_ANALYSIS_MY); - metadataCopySingle(stats, outRO->analysis, PM_SUBTRACTION_ANALYSIS_MXX); - metadataCopySingle(stats, outRO->analysis, PM_SUBTRACTION_ANALYSIS_MXY); - metadataCopySingle(stats, outRO->analysis, PM_SUBTRACTION_ANALYSIS_MYY); - - psMetadataAddF32(stats, PS_LIST_TAIL, "TIME_MATCH", 0, "Time to match PSFs", - psTimerClear("PPSUB_MATCH")); - } - } - - -#ifdef TESTING - psImage *kernelImage = psMetadataLookupPtr(&mdok, inConv->analysis, - "SUBTRACTION.KERNEL.IMAGE"); // Image of the kernels - if (!kernelImage) { - kernelImage = psMetadataLookupPtr(&mdok, refConv->analysis, "SUBTRACTION.KERNEL.IMAGE"); - } - psFits *fits = psFitsOpen("kernel.fits", "w"); - psFitsWriteImage(fits, NULL, kernelImage, 0, NULL); - psFitsClose(fits); -#endif - - // Subtraction is: minuend - subtrahend - pmReadout *minuend = inConv; - pmReadout *subtrahend = refConv; - - if (reverse) { - pmReadout *temp = subtrahend; - subtrahend = minuend; - minuend = temp; - } - -#ifdef TESTING - { - pmReadoutMaskApply(minuend, maskVal); - psFits *fits = psFitsOpen("minuend.fits", "w"); - psFitsWriteImage(fits, NULL, minuend->image, 0, NULL); - psFitsClose(fits); - } - { - pmReadoutMaskApply(subtrahend, maskVal); - psFits *fits = psFitsOpen("subtrahend.fits", "w"); - psFitsWriteImage(fits, NULL, subtrahend->image, 0, NULL); - psFitsClose(fits); - } -#endif - - // Photometry is to be performed in two stages: - // 1. Measure the PSF using the PSF-matched images - // 2. Find and measure sources on the subtracted image - psTimerStart("PPSUB_PHOT"); - - // Photometry stage 1: measure the PSF - pmPSF *psf = NULL; // PSF for photometry - if (psMetadataLookupBool(NULL, recipe, "PHOTOMETRY")) { - outRO->image = psImageCopy(outRO->image, minuend->image, PS_TYPE_F32); - if (minuend->weight) { - outRO->weight = psImageCopy(outRO->weight, minuend->weight, PS_TYPE_F32); - } - if (minuend->mask) { - outRO->mask = psImageCopy(outRO->mask, minuend->mask, PS_TYPE_IMAGE_MASK); - } - outRO->data_exists = outCell->data_exists = outCell->parent->data_exists = true; - - if (psMetadataLookupBool(&mdok, recipe, "RENORM")) { - psImageMaskType maskValue = pmConfigMaskGet("BLANK", config); // Bits to mask - if (!pmReadoutWeightRenormPhot(outRO, maskValue, renormNum, renormWidth, - renormMean, renormStdev, NULL)) { - psError(PS_ERR_UNKNOWN, false, "Unable to renormalise variances."); - psFree(outRO); - return false; - } - } - - pmFPAfile *photFile = psMetadataLookupPtr(NULL, config->files, "PSPHOT.INPUT"); - pmFPACopy(photFile->fpa, outRO->parent->parent->parent); - - // We have a list of sources, so we could use those to redetermine the PSF model. - // Until Gene makes the necessary adaptations to psphot, we will simply redetermine the PSF model from - // scratch. - - if (psMetadataLookupBool(&mdok, recipe, "PSPHOT.VISUAL")) { - psphotSetVisual(true); - } - - // Need to ensure aperture residual is not calculated - psMetadata *psphotRecipe = psMetadataLookupMetadata(NULL, config->recipes, PSPHOT_RECIPE); // Recipe - if (!psphotRecipe) { - psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to find %s recipe.", PSPHOT_RECIPE); - return false; - } - const char *breakpoint = psMetadataLookupStr(&mdok, psphotRecipe, "BREAK_POINT"); // Current break point - psMetadataAddStr(psphotRecipe, PS_LIST_TAIL, "BREAK_POINT", PS_META_REPLACE, - "ALTERED break point for psphot operations", "PSFMODEL"); - - // set maskValue and markValue in the psphot recipe - psImageMaskType maskValue = maskVal; - psImageMaskType markValue = pmConfigMaskGet("MARK.VALUE", config); // Bits to use for marking - psMetadataAddImageMask(psphotRecipe, PS_LIST_TAIL, "MASK.PSPHOT", PS_META_REPLACE, "Bits to mask", maskValue); - psMetadataAddImageMask(psphotRecipe, PS_LIST_TAIL, "MARK.PSPHOT", PS_META_REPLACE, "Bits to use for marking", - markValue); - - if (!psphotReadout(config, view)) { - psWarning("Unable to perform photometry on subtracted image."); - psErrorStackPrint(stderr, "Error stack from photometry:"); - psErrorClear(); - } - - if (!pmFPAfileDropInternal(config->files, "PSPHOT.BACKMDL") || - !pmFPAfileDropInternal (config->files, "PSPHOT.BACKMDL.STDEV") || - !pmFPAfileDropInternal (config->files, "PSPHOT.BACKGND")) { - psError(PS_ERR_UNKNOWN, false, "Unable to drop PSPHOT internal files."); - return false; - } - - if (breakpoint) { - psMetadataAddStr(psphotRecipe, PS_LIST_TAIL, "BREAK_POINT", PS_META_REPLACE, - "RESTORED break point for psphot operations", breakpoint); - } - - // Blow away the sources psphot found --- they're irrelevant for the subtraction - pmReadout *photRO = pmFPAviewThisReadout(view, photFile->fpa); // Readout with sources - psMetadataRemoveKey(photRO->analysis, "PSPHOT.SOURCES"); - psMetadataRemoveKey(photRO->analysis, "PSPHOT.HEADER"); - - psf = psMemIncrRefCounter(psMetadataLookupPtr(NULL, photRO->parent->parent->analysis, "PSPHOT.PSF")); - if (!psf) { - psError(PS_ERR_UNEXPECTED_NULL, true, "Unable to find PSF from psphot"); - return false; - } - - // Record the FWHM - pmHDU *hdu = pmHDUFromCell(outRO->parent); // HDU with header - psMetadataItemSupplement(hdu->header, psphotRecipe, "FWHM_MAJ"); - psMetadataItemSupplement(hdu->header, psphotRecipe, "FWHM_MIN"); - - - pmCell *photCell = pmFPAfileThisCell(config->files, view, "PSPHOT.INPUT"); - pmCellFreeReadouts(photCell); - } - - // Do the actual subtraction - outRO->image = (psImage*)psBinaryOp(outRO->image, minuend->image, "-", subtrahend->image); - if (inConv->weight && refConv->weight) { - outRO->weight = (psImage*)psBinaryOp(outRO->weight, inConv->weight, "+", refConv->weight); - } - outRO->mask = (psImage*)psBinaryOp(outRO->mask, inConv->mask, "|", refConv->mask); - - psFree(inConv); - psFree(refConv); - - outRO->data_exists = outCell->data_exists = outCell->parent->data_exists = true; - - // Copy astrometry over - // It should get into the output images and photometry - pmFPA *refFPA = refRO->parent->parent->parent; // Reference FPA - pmHDU *refHDU = refFPA->hdu; // Reference HDU - if (!outHDU || !refHDU) { - psWarning("Unable to find HDU at FPA level to copy astrometry."); - } else if (!pmAstromReadWCS(outFPA, outCell->parent, refHDU->header, 1.0)) { - psWarning("Unable to read WCS astrometry from reference FPA."); - psErrorClear(); - } else if (!pmAstromWriteWCS(outHDU->header, outFPA, outCell->parent, WCS_TOLERANCE)) { - psWarning("Unable to write WCS astrometry to output FPA."); - psErrorClear(); - } - -#if 0 - pmReadoutMaskApply(outRO, maskBad); -#endif - -#ifdef TESTING - for (int y = 0; y < outRO->image->numRows; y++) { - for (int x = 0; x < outRO->image->numCols; x++) { - if (isnan(outRO->image->data.F32[y][x]) && !(outRO->mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & maskVal)) { - printf("Unmasked NAN at %d %d --> %d\n", x, y, outRO->mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x]); - } - } - } -#endif - - if (!pmFPACopyConcepts(outCell->parent->parent, inRO->parent->parent->parent)) { - psError(PS_ERR_UNKNOWN, false, "Unable to copy concepts from input to output."); - psFree(outRO); + if (!ppSubDefineOutput (config, view)) { + psError(PS_ERR_UNKNOWN, false, "Unable to subtract background."); return false; } -#if 0 - // XXX Under development + if (!ppSubVarianceFactors (config, stats, view)) { + psError(PS_ERR_UNKNOWN, false, "Unable to subtract background."); + return false; + } - // QA: photometry of known sources with measured PSF - if (sourcesRO && psMetadataLookupBool(NULL, recipe, "PHOTOMETRY")) { - sources = psMetadataLookupPtr(&mdok, sourcesRO->analysis, "PSPHOT.SOURCES"); + if (!ppSubMakePSF(config, view)) { + psError(PS_ERR_UNKNOWN, false, "Unable to subtract background."); + return false; + } - - - // For each source, need to: - // * generate appropriate model - // * select pixels - // * define fit radius - // Don't need to measure sky (psphotFitSourcesLinear does that) - - psArray *dummy = psArrayAlloc(1); // Dummy array of sources, for psphotFitSourcesLinear with 1 source - for (long i = 0; i < sources->n; i++) { - pmSource *source = sources->data[i]; - if (!source || (source->mode & SOURCE_MASK)) { - continue; - } - - float origMag = source->psfMag; // Original magnitude - - // Coordinates of source - float x = source->modelPSF->params->data.F32[PM_PAR_XPOS]; - float y = source->modelPSF->params->data.F32[PM_PAR_YPOS]; - - model = pmModelFromPSFforXY(psf, x, y, NAN); - if (!fakeModel) { - psWarning("Can't generate model"); - continue; - } - if (!pmSourceDefinePixels(fakeSource, readout, x, y, fakeRadius)) { - psErrorClear(); - continue; - } -pmModel *pmModelFromPSFforXY ( - const pmPSF *psf, - float Xo, - float Yo, - float Io - ); - -bool pmSourceDefinePixels( - pmSource *mySource, ///< source to be re-defined - const pmReadout *readout, ///< base the source on this readout - psF32 x, ///< center coords of source - psF32 y, ///< center coords of source - psF32 Radius ///< size of box on source -); - - -bool psphotFitSourcesLinear (pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf, bool final) { - - } -#endif - - // Photometry stage 2: find and measure sources on the subtracted image - if (psMetadataLookupBool(NULL, recipe, "PHOTOMETRY")) { - // The PSF should already be stored for the readout - pmFPAfile *photFile = psMetadataLookupPtr(NULL, config->files, "PSPHOT.INPUT"); - pmFPACopy(photFile->fpa, outFPA); - - pmReadout *psfRO = pmFPAfileThisReadout(config->files, view, "PSPHOT.PSF.LOAD"); - if (!psfRO) { - psError(PS_ERR_UNEXPECTED_NULL, true, "Unable to find readout for file PSPHOT.PSF.LOAD"); - return false; - } - psMetadataAddPtr(psfRO->parent->parent->analysis, PS_LIST_TAIL, "PSPHOT.PSF", PS_DATA_UNKNOWN, - "PSF from matched addition", psf); - psFree(psf); - - if (psMetadataLookupBool(&mdok, recipe, "PSPHOT.VISUAL")) { - psphotSetVisual(true); - } - - if (psMetadataLookupBool(&mdok, recipe, "RENORM")) { - psImageMaskType maskValue = pmConfigMaskGet("BLANK", config); // Bits to mask - if (!pmReadoutWeightRenormPhot(outRO, maskValue, renormNum, renormWidth, - renormMean, renormStdev, NULL)) { - psError(PS_ERR_UNKNOWN, false, "Unable to renormalise variances."); - psFree(outRO); - return false; - } - } - - // Need to ensure aperture residual is not calculated - psMetadata *psphotRecipe = psMetadataLookupMetadata(NULL, config->recipes, PSPHOT_RECIPE); // Recipe - psMetadataItem *item = psMetadataLookup(psphotRecipe, "MEASURE.APTREND"); // Item determining aptrend - if (!item) { - psWarning("Unable to find MEASURE.APTREND in psphot recipe"); - psErrorClear(); - } else { - item->data.B = false; - } - - if (!psphotReadout(config, view)) { - psWarning("Unable to perform photometry on subtracted image."); - psErrorStackPrint(stderr, "Error stack from photometry:"); - psErrorClear(); - } - - if (!pmFPAfileDropInternal(config->files, "PSPHOT.BACKMDL") || - !pmFPAfileDropInternal (config->files, "PSPHOT.BACKMDL.STDEV") || - !pmFPAfileDropInternal (config->files, "PSPHOT.BACKGND")) { - psError(PS_ERR_UNKNOWN, false, "Unable to drop PSPHOT internal files."); - return false; - } - - pmFPAfileActivate(config->files, false, "PSPHOT.INPUT"); - pmFPAfileActivate(config->files, false, "PSPHOT.LOAD.PSF"); - - if (stats) { - pmReadout *photRO = pmFPAviewThisReadout(view, photFile->fpa); // Readout with the sources - psArray *sources = psMetadataLookupPtr(NULL, photRO->analysis, "PSPHOT.SOURCES"); // Sources - psMetadataAddS32(stats, PS_LIST_TAIL, "NUM_SOURCES", 0, "Number of sources detected", - sources ? sources->n : 0); - psMetadataAddF32(stats, PS_LIST_TAIL, "TIME_PHOT", 0, "Time to do photometry", - psTimerClear("PPSUB_PHOT")); - } - -#ifdef TESTING - // Record data about sources: not everything gets into the output CMF files - { - pmReadout *photRO = pmFPAviewThisReadout(view, photFile->fpa); // Readout with the sources - psArray *sources = psMetadataLookupPtr(NULL, photRO->analysis, "PSPHOT.SOURCES"); // Sources - FILE *sourceFile = fopen("sources.dat", "w"); // File for sources - fprintf(sourceFile, - "# x y mag mag_err psf_chisq cr_nsigma ext_nsigma psf_qf flags m_x m_y m_xx m_xy m_yy\n"); - for (int i = 0; i < sources->n; i++) { - pmSource *source = sources->data[i]; - if (!source) { - continue; - } - - float x, y; // Position of source - float chi2; // chi^2 for source - if (source->modelPSF) { - x = source->modelPSF->params->data.F32[PM_PAR_XPOS]; - y = source->modelPSF->params->data.F32[PM_PAR_YPOS]; - chi2 = source->modelPSF->chisq; - } else if (source->peak) { - x = source->peak->xf; - y = source->peak->yf; - chi2 = NAN; - } else { - psWarning("No position available for source."); - continue; - } - - float xMoment = NAN, yMoment = NAN, xxMoment = NAN, xyMoment = NAN, yyMoment = NAN; - if (source->moments) { - xMoment = source->moments->Mx; - yMoment = source->moments->My; - xxMoment = source->moments->Mxx; - xyMoment = source->moments->Mxy; - yyMoment = source->moments->Myy; - } - - fprintf(sourceFile, "%f %f %f %f %f %f %f %f %d %f %f %f %f %f\n", - x, y, source->psfMag, source->errMag, chi2, source->crNsigma, source->extNsigma, - source->pixWeight, source->mode, xMoment, yMoment, xxMoment, xyMoment, yyMoment); - } - fclose(sourceFile); - } -#endif - + if (!ppSubReadoutSubtract (config, view)) { + psError(PS_ERR_UNKNOWN, false, "Unable to subtract background."); + return false; } // Higher order background subtraction using psphot - if (!ppSubBackground(outRO, config, view)) { + if (!ppSubBackground(config, view)) { psError(PS_ERR_UNKNOWN, false, "Unable to subtract background."); - psFree(outRO); + return false; + } + + if (!ppSubReadoutPhotometry (config, stats, view)) { + psError(PS_ERR_UNKNOWN, false, "Unable to subtract background."); return false; } - // Renormalising for pixels, because that's what magic desires - if (psMetadataLookupBool(&mdok, recipe, "RENORM")) { - psImageMaskType maskValue = pmConfigMaskGet("BLANK", config); // Bits to mask - if (!pmReadoutWeightRenormPixels(outRO, maskValue, renormMean, renormStdev, NULL)) { - psError(PS_ERR_UNKNOWN, false, "Unable to renormalise variances."); - psFree(outRO); - return false; - } + if (!ppSubReadoutUpdate (config, view)) { + psError(PS_ERR_UNKNOWN, false, "Unable to subtract background."); + return false; } - - // Add additional data to the header - pmFPAfile *refFile = psMetadataLookupPtr(NULL, config->files, "PPSUB.REF"); // Reference file - pmFPAfile *inFile = psMetadataLookupPtr(NULL, config->files, "PPSUB.INPUT"); // Input file - psMetadataAddStr(outHDU->header, PS_LIST_TAIL, "PPSUB.REFERENCE", 0, - "Subtraction reference", refFile->filename); - psMetadataAddStr(outHDU->header, PS_LIST_TAIL, "PPSUB.INPUT", 0, - "Subtraction input", inFile->filename); - - - // Generate binned JPEGs - { - int bin1 = psMetadataLookupS32(NULL, recipe, "BIN1"); // First binning level - int bin2 = psMetadataLookupS32(NULL, recipe, "BIN2"); // Second binning level - - // Target cells - pmCell *cell1 = pmFPAfileThisCell(config->files, view, "PPSUB.OUTPUT.JPEG1"); - pmCell *cell2 = pmFPAfileThisCell(config->files, view, "PPSUB.OUTPUT.JPEG2"); - - pmReadout *ro1 = pmReadoutAlloc(cell1), *ro2 = pmReadoutAlloc(cell2); // Binned readouts - if (!pmReadoutRebin(ro1, outRO, maskBad, bin1, bin1) || !pmReadoutRebin(ro2, ro1, 0, bin2, bin2)) { - psError(PS_ERR_UNKNOWN, false, "Unable to bin output."); - psFree(ro1); - psFree(ro2); - psFree(outRO); - return false; - } - psFree(ro1); - psFree(ro2); - } - -#ifdef TESTING - // Significance image - { - psImage *sig = (psImage*)psBinaryOp(NULL, outRO->image, "*", outRO->image); - psBinaryOp(sig, sig, "/", outRO->weight); - psFits *fits = psFitsOpen("significance.fits", "w"); - psFitsWriteImage(fits, NULL, sig, 0, NULL); - psFitsClose(fits); - psFree(sig); - } -#endif - - psFree(outRO); return true;