Index: trunk/ppSub/src/Makefile.am =================================================================== --- trunk/ppSub/src/Makefile.am (revision 21183) +++ trunk/ppSub/src/Makefile.am (revision 21257) @@ -4,12 +4,22 @@ ppSub_LDFLAGS = $(PSLIB_LIBS) $(PSMODULE_LIBS) $(PPSTATS_LIBS) $(PSPHOT_LIBS) $(PPSUB_LIBS) -ppSub_SOURCES = \ - ppSub.c \ - ppSubArguments.c \ - ppSubBackground.c \ - ppSubCamera.c \ - ppSubLoop.c \ - ppSubReadout.c \ - ppSubVersion.c +ppSub_SOURCES = \ + ppSub.c \ + ppSubArguments.c \ + ppSubVersion.c \ + ppSubBackground.c \ + ppSubCamera.c \ + ppSubLoop.c \ + ppSubReadout.c \ + ppSubDefineOutput.c \ + ppSubExtras.c \ + ppSubMakePSF.c \ + ppSubMatchPSFs.c \ + ppSubReadoutPhotometry.c \ + ppSubReadoutSubtract.c \ + ppSubReadoutUpdate.c \ + ppSubSetMasks.c \ + ppSubReadoutRenorm.c \ + ppSubVarianceFactors.c ppSubKernel_CPPFLAGS = $(PSLIB_CFLAGS) $(PSMODULE_CFLAGS) $(PPSUB_CFLAGS) Index: trunk/ppSub/src/ppSub.c =================================================================== --- trunk/ppSub/src/ppSub.c (revision 21183) +++ trunk/ppSub/src/ppSub.c (revision 21257) @@ -1,11 +1,2 @@ -#ifdef HAVE_CONFIG_H -#include -#endif - -#include -#include -#include -#include - #include "ppSub.h" Index: trunk/ppSub/src/ppSub.h =================================================================== --- trunk/ppSub/src/ppSub.h (revision 21183) +++ trunk/ppSub/src/ppSub.h (revision 21257) @@ -1,4 +1,15 @@ #ifndef PP_SUB_H #define PP_SUB_H + +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include +#include +#include +#include +#include #define PPSUB_RECIPE "PPSUB" /// Name of the recipe to use @@ -6,5 +17,5 @@ /// Setup the arguments parsing bool ppSubArgumentsSetup(int argc, char *argv[], ///< Command-line arguments - pmConfig *config ///< Configuration + pmConfig *config ///< Configuration ); @@ -27,7 +38,51 @@ ); +/// generate (if needed) and set or update the masks for input and reference images +bool ppSubSetMasks (pmConfig *config, ///< Configuration + const pmFPAview *view ///< View of active readout + ); + +/// Generate the PSF-matching kernel and convolve the images as needed. Most of this function +/// involves looking up the parameters in the recipe and supplying them to the function +/// pmSubtractionMatch() +bool ppSubMatchPSFs (pmConfig *config, ///< Configuration + const pmFPAview *view ///< View of active readout + ); + +/// generate the output readout and pass the kernel info to the header +bool ppSubDefineOutput(pmConfig *config, ///< Configuration + const pmFPAview *view ///< View of active readout + ); + +/// Calculate the variance factor for the output image based on the input images +bool ppSubVarianceFactors(pmConfig *config, ///< Configuration + psMetadata *stats, ///< Statistics, for output + const pmFPAview *view ///< View of active readout + ); + +/// Photometry stage 1: measure the PSF from the minuend image +bool ppSubMakePSF(pmConfig *config, ///< Configuration + const pmFPAview *view ///< View of active readout + ); + +/// Perform the actual image subtraction, update output concepts +bool ppSubReadoutSubtract(pmConfig *config, ///< Configuration + const pmFPAview *view ///< View of active readout + ); + + +/// Photometry stage 2: find and measure sources on the subtracted image +bool ppSubReadoutPhotometry(pmConfig *config, ///< Configuration + psMetadata *stats, ///< Statistics, for output + const pmFPAview *view ///< View of active readout + ); + +/// Renormalize, update headers and generate JPEGs +bool ppSubReadoutUpdate(pmConfig *config, ///< Configuration + const pmFPAview *view ///< View of active readout + ); + /// Higher-order background subtraction -bool ppSubBackground(pmReadout *ro, ///< Readout of interest - pmConfig *config, ///< Configuration +bool ppSubBackground(pmConfig *config, ///< Configuration const pmFPAview *view ///< View to readout ); @@ -37,4 +92,26 @@ ); +/// Generate and Set the masks if needed +bool ppSubSetMasks ( + pmConfig *config, ///< Configuration + const pmFPAview *view ///< view to readout + ); + +/// Renormalize readout for peak pixels +bool ppSubReadoutRenormPixels ( + pmConfig *config, ///< Configuration + psMetadata *recipe, ///< Recipe + pmReadout *readout ///< Readout + ); + +/// Renormalize readout for photometry analysis +bool ppSubReadoutRenormPhot ( + pmConfig *config, ///< Configuration + psMetadata *recipe, ///< Recipe + pmReadout *readout ///< Readout + ); + +// Copy every instance of a single keyword from one metadata to another +bool psMetadataCopySingle(psMetadata *target, psMetadata *source, const char *name); #endif Index: trunk/ppSub/src/ppSubArguments.c =================================================================== --- trunk/ppSub/src/ppSubArguments.c (revision 21183) +++ trunk/ppSub/src/ppSubArguments.c (revision 21257) @@ -1,12 +1,2 @@ -#ifdef HAVE_CONFIG_H -#include -#endif - -#include -#include -#include -#include -#include - #include "ppSub.h" @@ -172,4 +162,5 @@ bool ppSubArgumentsSetup(int argc, char *argv[], pmConfig *config) { + int argnum; // Argument Number assert(config); @@ -180,10 +171,7 @@ } - { - int arg; // Argument Number - if ((arg = psArgumentGet(argc, argv, "-psphot-visual"))) { - psArgumentRemove(arg, &argc, argv); - psMetadataAddBool(recipe, PS_LIST_TAIL, "PSPHOT.VISUAL", 0, "Visual guide to psphot?", true); - } + if ((argnum = psArgumentGet(argc, argv, "-psphot-visual"))) { + psArgumentRemove(argnum, &argc, argv); + psphotSetVisual(true); } @@ -381,4 +369,13 @@ } + // XXX move this to ppSubArguments + 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; + } + } + return true; Index: trunk/ppSub/src/ppSubBackground.c =================================================================== --- trunk/ppSub/src/ppSubBackground.c (revision 21183) +++ trunk/ppSub/src/ppSubBackground.c (revision 21257) @@ -1,41 +1,27 @@ -#ifdef HAVE_CONFIG_H -#include -#endif - -#include -#include -#include -#include #include "ppSub.h" -// Copied from ppImageSubtractBackground() -bool ppSubBackground(pmReadout *ro, pmConfig *config, const pmFPAview *view) +// Based on ppImageSubtractBackground() +bool ppSubBackground(pmConfig *config, const pmFPAview *view) { psAssert(config, "Need configuration"); psAssert(view, "Need view to chip"); - // The background model file may be defined, though the model hasn't been built - // If so, we will build it below. - bool status; // Status of lookup - pmFPAfile *modelFile = psMetadataLookupPtr(&status, config->files, "PSPHOT.BACKMDL"); // Background model - if (!status) { - psError(PS_ERR_UNEXPECTED_NULL, true, "PSPHOT.BACKMDL file is not defined"); - return false; - } + bool status; // Status of metadata lookups psMetadata *ppSubRecipe = psMetadataLookupPtr(NULL, config->recipes, PPSUB_RECIPE); psAssert(ppSubRecipe, "Need PPSUB recipe"); + psMetadata *psphotRecipe = psMetadataLookupPtr(NULL, config->recipes, PSPHOT_RECIPE); - psAssert(psphotRecipe, "Need PSPHOT recipe"); + psAssert(psphotRecipe, "Need PSPHOT recipe for binning"); - psString maskBadStr = psMetadataLookupStr(NULL, ppSubRecipe, "MASK.BAD"); // Name of bits to mask for bad - psImageMaskType maskBad = pmConfigMaskGet(maskBadStr, config); // Bits to mask for bad pixels + psImageMaskType maskBad = pmConfigMaskGet("BLANK", config); - // user-defined masks to test for good/bad pixels (build from recipe list if not yet set) - psMetadataAddImageMask(psphotRecipe, PS_LIST_TAIL, "MASK.PSPHOT", PS_META_REPLACE, "user-defined mask", maskBad); + // select the output readout + pmReadout *outRO = pmFPAfileThisReadout(config->files, view, "PPSUB.OUTPUT"); - psImage *image = ro->image, *mask = ro->mask; // Image and mask of interest + // select the model readout, if it exist already; if not, generate it + pmReadout *modelRO = pmFPAfileThisReadout(config->files, view, "PSPHOT.BACKMDL"); - pmReadout *modelRO = pmFPAviewThisReadout(view, modelFile->fpa); // Background model + // if necessary, generate the background model if (!modelRO) { // Create the background model @@ -44,6 +30,6 @@ return false; } - modelRO = (modelFile->mode == PM_FPA_MODE_INTERNAL) ? modelFile->readout : - pmFPAviewThisReadout(view, modelFile->fpa); + // select the model readout (should now exist) + modelRO = pmFPAfileThisReadout(config->files, view, "PSPHOT.BACKMDL"); if (!modelRO) { psError(PS_ERR_UNEXPECTED_NULL, true, "Unable to find background model"); @@ -54,4 +40,7 @@ "PSPHOT.BACKGROUND.BINNING"); // Binning for model psImage *modelImage = modelRO->image; // Background model + + psImage *image = outRO->image; // Image of interest + psImage *mask = outRO->mask; // Mask of interest // Do the background subtraction Index: trunk/ppSub/src/ppSubCamera.c =================================================================== --- trunk/ppSub/src/ppSubCamera.c (revision 21183) +++ trunk/ppSub/src/ppSubCamera.c (revision 21257) @@ -1,16 +1,7 @@ -#ifdef HAVE_CONFIG_H -#include -#endif - -#include -#include -#include -#include - #include "ppSub.h" bool ppSubCamera(pmConfig *config) { - bool status = false; // Status of definition + bool status = false; // result from pmFPAfileDefine operations // Input image @@ -119,4 +110,80 @@ } + // Convolved input image + pmFPAfile *inConv = pmFPAfileDefineFromFile(config, input, 1, 1, "PPSUB.INPUT.CONV"); + if (!inConv) { + psError(PS_ERR_IO, false, _("Unable to generate output file for PPSUB.INPUT.CONV")); + return false; + } + if (output->type != PM_FPA_FILE_IMAGE) { + psError(PS_ERR_IO, true, "PPSUB.INPUT.CONV is not of type IMAGE"); + return false; + } + // XXX should be based on recipe : inConv->save = true; + + // Convolved input mask + pmFPAfile *inConvMask = pmFPAfileDefineOutput(config, inConv->fpa, "PPSUB.INPUT.CONV.MASK"); + if (!inConvMask) { + psError(PS_ERR_IO, false, _("Unable to generate output file from PPSUB.INPUT.CONV.MASK")); + return false; + } + if (inConvMask->type != PM_FPA_FILE_MASK) { + psError(PS_ERR_IO, true, "PPSUB.INPUT.CONV.MASK is not of type MASK"); + return false; + } + // XXX should be based on recipe : inConvMask->save = true; + + // Convolved input weight + if (inputWeight) { + pmFPAfile *inConvWeight = pmFPAfileDefineOutput(config, inConv->fpa, "PPSUB.INPUT.CONV.WEIGHT"); + if (!inConvWeight) { + psError(PS_ERR_IO, false, _("Unable to generate output file from PPSUB.OUTPUT.WEIGHT")); + return false; + } + if (inConvWeight->type != PM_FPA_FILE_WEIGHT) { + psError(PS_ERR_IO, true, "PPSUB.INPUT.CONV.WEIGHT is not of type WEIGHT"); + return false; + } + // XXX should be based on recipe : inConvWeight->save = true; + } + + // Convolved ref image + pmFPAfile *refConv = pmFPAfileDefineFromFile(config, ref, 1, 1, "PPSUB.REF.CONV"); + if (!refConv) { + psError(PS_ERR_IO, false, _("Unable to generate output file for PPSUB.REF.CONV")); + return false; + } + if (output->type != PM_FPA_FILE_IMAGE) { + psError(PS_ERR_IO, true, "PPSUB.REF.CONV is not of type IMAGE"); + return false; + } + // XXX should be based on recipe : refConv->save = true; + + // Convolved ref mask + pmFPAfile *refConvMask = pmFPAfileDefineOutput(config, refConv->fpa, "PPSUB.REF.CONV.MASK"); + if (!refConvMask) { + psError(PS_ERR_IO, false, _("Unable to generate output file from PPSUB.REF.CONV.MASK")); + return false; + } + if (refConvMask->type != PM_FPA_FILE_MASK) { + psError(PS_ERR_IO, true, "PPSUB.REF.CONV.MASK is not of type MASK"); + return false; + } + // XXX should be based on recipe : refConvMask->save = true; + + // Convolved ref weight + if (refWeight) { + pmFPAfile *refConvWeight = pmFPAfileDefineOutput(config, refConv->fpa, "PPSUB.REF.CONV.WEIGHT"); + if (!refConvWeight) { + psError(PS_ERR_IO, false, _("Unable to generate output file from PPSUB.OUTPUT.WEIGHT")); + return false; + } + if (refConvWeight->type != PM_FPA_FILE_WEIGHT) { + psError(PS_ERR_IO, true, "PPSUB.REF.CONV.WEIGHT is not of type WEIGHT"); + return false; + } + // XXX should be based on recipe : refConvWeight->save = true; + } + // Output JPEGs pmFPAfile *jpeg1 = pmFPAfileDefineOutput(config, NULL, "PPSUB.OUTPUT.JPEG1"); @@ -192,8 +259,16 @@ } pmFPAfileActivate(config->files, false, "PSPHOT.PSF.LOAD"); + + pmFPAfile *psphotResid = pmFPAfileDefineOutputForFormat (config, NULL, "PSPHOT.RESID", output->cameraName, output->formatName); + if (!psphotResid) { + psError(PSPHOT_ERR_CONFIG, false, "Cannot find a rule for PSPHOT.RESID"); + return false; + } + // make this optional: psphotResid->save = true; } // Always do this, since we're using psphot for the background model. - { + // XXX Not certain that I need to generate a separate pmFPAfile for PSPHOT.INPUT + if (0) { pmFPAfile *psphot = pmFPAfileDefineFromFPA(config, output->fpa, 1, 1, "PSPHOT.INPUT"); if (!psphot) { Index: trunk/ppSub/src/ppSubDefineOutput.c =================================================================== --- trunk/ppSub/src/ppSubDefineOutput.c (revision 21257) +++ trunk/ppSub/src/ppSubDefineOutput.c (revision 21257) @@ -0,0 +1,58 @@ +#include "ppSub.h" + +bool ppSubDefineOutput (pmConfig *config, const pmFPAview *view) { + + bool mdok; + + // generate an output readout + 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(); + } + + // convolved input images + pmReadout *inConv = pmFPAfileThisReadout(config->files, view, "PPSUB.INPUT.CONV"); // Input readout + pmReadout *refConv = pmFPAfileThisReadout(config->files, view, "PPSUB.REF.CONV"); // Reference readout + + // Add kernel descrption to header. We don't know which readout gets the kernels because + // it depends on which is larger (the choice is set in ppSubMatchPSFs) + bool inputHasKernel = true; + pmSubtractionKernels *kernels = psMetadataLookupPtr(&mdok, inConv->analysis, PM_SUBTRACTION_ANALYSIS_KERNEL); // Subtraction kernel + if (!kernels) { + kernels = psMetadataLookupPtr(&mdok, refConv->analysis, PM_SUBTRACTION_ANALYSIS_KERNEL); + inputHasKernel = false; + } + 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 (inputHasKernel) { + outRO->analysis = psMetadataCopy(outRO->analysis, inConv->analysis); + } else { + outRO->analysis = psMetadataCopy(outRO->analysis, refConv->analysis); + } + +#ifdef TESTING + psImage *kernelImage = NULL; + + if (inputHasKernel) { + kernelImage = psMetadataLookupPtr(&mdok, inConv->analysis, "SUBTRACTION.KERNEL.IMAGE"); // Image of the kernels + } else { + kernelImage = psMetadataLookupPtr(&mdok, refConv->analysis, "SUBTRACTION.KERNEL.IMAGE"); + } + + psFits *fits = psFitsOpen("kernel.fits", "w"); + psFitsWriteImage(fits, NULL, kernelImage, 0, NULL); + psFitsClose(fits); +#endif + + return true; +} Index: trunk/ppSub/src/ppSubExtras.c =================================================================== --- trunk/ppSub/src/ppSubExtras.c (revision 21257) +++ trunk/ppSub/src/ppSubExtras.c (revision 21257) @@ -0,0 +1,24 @@ +#include "ppSub.h" + +// This file contains minor functions used in ppSub. some of these should be pushed into +// psModules or psLib. + +// Copy every instance of a single keyword from one metadata to another +bool psMetadataCopySingle(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; +} Index: trunk/ppSub/src/ppSubLoop.c =================================================================== --- trunk/ppSub/src/ppSubLoop.c (revision 21183) +++ trunk/ppSub/src/ppSubLoop.c (revision 21257) @@ -1,12 +1,2 @@ -#ifdef HAVE_CONFIG_H -#include -#endif - -#include -#include -#include -#include -#include - #include "ppSub.h" @@ -28,9 +18,4 @@ } psFree(resolved); - } - - if (!pmConfigMaskSetBits(NULL, NULL, config)) { - psError(PS_ERR_UNKNOWN, false, "Unable to determine mask value."); - return false; } @@ -141,5 +126,6 @@ return false; } - ppStatsFPA(stats, output->fpa, view, pmConfigMaskGet("MASK.VALUE", config), config); + psImageMaskType maskValue = pmConfigMaskGet("MASK.VALUE", config); + ppStatsFPA(stats, output->fpa, view, maskValue, config); } Index: trunk/ppSub/src/ppSubMakePSF.c =================================================================== --- trunk/ppSub/src/ppSubMakePSF.c (revision 21257) +++ trunk/ppSub/src/ppSubMakePSF.c (revision 21257) @@ -0,0 +1,79 @@ +#include "ppSub.h" + +// Photometry stage 1: measure the PSF from the minuend image +bool ppSubMakePSF (pmConfig *config, const pmFPAview *view) { + + bool mdok = false; + + // 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"); + + // Look up recipe values + psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, PPSUB_RECIPE); // Recipe for ppSub + psAssert(recipe, "We checked this earlier, so it should be here."); + + psMetadata *psphotRecipe = psMetadataLookupMetadata(NULL, config->recipes, PSPHOT_RECIPE); // Recipe for psphot + psAssert(recipe, "We checked this earlier, so it should be here."); + + if (!psMetadataLookupBool(NULL, recipe, "PHOTOMETRY")) return true; + + bool reverse = psMetadataLookupBool(NULL, config->arguments, "REVERSE"); // Reverse sense of subtraction? + + pmReadout *minuend = NULL; + pmFPAfile *minuendFile = NULL; + if (reverse) { + minuend = pmFPAfileThisReadout(config->files, view, "PPSUB.REF.CONV"); + minuendFile = psMetadataLookupPtr(&mdok, config->files, "PPSUB.REF.CONV"); + } else { + minuend = pmFPAfileThisReadout(config->files, view, "PPSUB.INPUT.CONV"); + minuendFile = psMetadataLookupPtr(&mdok, config->files, "PPSUB.INPUT.CONV"); + } + + // supply the minuend pmFPAfile to psphot as PSPHOT.INPUT: + psMetadataAddPtr (config->files, PS_LIST_TAIL, "PSPHOT.INPUT", PS_DATA_UNKNOWN | PS_META_REPLACE, "psphot input : view on another pmFPAfile", minuendFile); + + // old-style variance renormalization + if (!ppSubReadoutRenormPhot (config, recipe, minuend)) { + psError(PS_ERR_UNKNOWN, false, "failure in renormalization"); + return false; + } + + // extract the loaded sources from the associated readout + pmReadout *sourcesRO = pmFPAfileThisReadout(config->files, view, "PPSUB.SOURCES"); + psArray *sources = psMetadataLookupPtr(&mdok, sourcesRO->analysis, "PSPHOT.SOURCES"); + + // generate PSF from the supplied sources (assumes image is background-subtracted) + if (!psphotReadoutFindPSF(config, view, sources)) { + psError(PS_ERR_UNKNOWN, false, "Unable to perform photometry on subtracted image."); + return false; + } + + // Record the FWHM in the output header + pmReadout *outRO = pmFPAfileThisReadout(config->files, view, "PPSUB.OUTPUT"); // Output readout + pmHDU *hdu = pmHDUFromCell(outRO->parent); // HDU with header + psMetadataItemSupplement(hdu->header, psphotRecipe, "FWHM_MAJ"); + psMetadataItemSupplement(hdu->header, psphotRecipe, "FWHM_MIN"); + + return true; +} + +// XXX we used to need this, is it still needed? + +// pmCell *photCell = pmFPAfileThisCell(config->files, view, "PSPHOT.INPUT"); +// pmCellFreeReadouts(photCell); + +// 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; +// } + +// Blow away the sources psphot found --- they're irrelevant for the subtraction +// XXX is this still needed? These are now probably not being set +// pmReadout *photRO = pmFPAviewThisReadout(view, photFile->fpa); // Readout with sources +// psMetadataRemoveKey(photRO->analysis, "PSPHOT.SOURCES"); +// psMetadataRemoveKey(photRO->analysis, "PSPHOT.HEADER"); + Index: trunk/ppSub/src/ppSubMatchPSFs.c =================================================================== --- trunk/ppSub/src/ppSubMatchPSFs.c (revision 21257) +++ trunk/ppSub/src/ppSubMatchPSFs.c (revision 21257) @@ -0,0 +1,121 @@ +#include "ppSub.h" + +// Generate the PSF-matching kernel and convolve the images as needed. Most of this function +// involves looking up the parameters in the recipe and supplying them to the function +// pmSubtractionMatch(); + +bool ppSubMatchPSFs (pmConfig *config, const pmFPAview *view) { + + bool mdok; // Status of MD lookup + + // 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."); + + // input images + pmReadout *inRO = pmFPAfileThisReadout(config->files, view, "PPSUB.INPUT"); // Input readout + pmReadout *refRO = pmFPAfileThisReadout(config->files, view, "PPSUB.REF"); // Reference readout + + // get or generate output image holders + pmReadout *inConv = pmFPAfileThisReadout(config->files, view, "PPSUB.INPUT.CONV"); // Input convolved readout + if (!inConv) { + pmCell *cell = pmFPAfileThisCell(config->files, view, "PPSUB.INPUT.CONV"); // Input convolved readout + inConv = pmReadoutAlloc(cell); // Convolved version of input + } + pmReadout *refConv = pmFPAfileThisReadout(config->files, view, "PPSUB.REF.CONV"); // Reference convolved readout + if (!refConv) { + pmCell *cell = pmFPAfileThisCell(config->files, view, "PPSUB.REF.CONV"); // Input convolved readout + refConv = pmReadoutAlloc(cell); // Convolved version of input + } + + // options which control the generation of optimal parameters + 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 + + // Vector with FWHMs for optimum search + psVector *optWidths = NULL; + if (optimum) { + optWidths = psVectorCreate(optWidths, optMin, optMax, optStep, PS_TYPE_F32); + } + + // Sources in image : these must be loaded from previous analysis stages + pmReadout *sourcesRO = pmFPAfileThisReadout(config->files, view, "PPSUB.SOURCES"); // Readout with sources + psArray *sources = NULL; + if (sourcesRO) { + sources = psMetadataLookupPtr(&mdok, sourcesRO->analysis, "PSPHOT.SOURCES"); + } + + int footprint = psMetadataLookupS32(NULL, recipe, "STAMP.FOOTPRINT"); // Stamp half-size + int stride = psMetadataLookupS32(NULL, recipe, "STRIDE"); // Size of convolution patches + float regionSize = psMetadataLookupF32(NULL, recipe, "REGION.SIZE"); // Size of iso-kernel regs + float spacing = psMetadataLookupF32(NULL, recipe, "STAMP.SPACING"); // Typical stamp spacing + float threshold = psMetadataLookupF32(NULL, recipe, "STAMP.THRESHOLD"); // Threshold for stmps + + const char *stampsName = psMetadataLookupStr(&mdok, config->arguments, "STAMPS"); // Filename for stamps + + 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); + return false; + } + + int size = psMetadataLookupS32(NULL, recipe, "KERNEL.SIZE"); // Kernel half-size + int order = psMetadataLookupS32(NULL, recipe, "SPATIAL.ORDER"); // Spatial polynomial order + 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 + + int optOrder = psMetadataLookupS32(&mdok, recipe, "OPTIMUM.ORDER"); // Order for search + float optThresh = psMetadataLookupF32(&mdok, recipe, "OPTIMUM.TOL"); // Tolerance for search + + 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 badFrac = psMetadataLookupF32(NULL, recipe, "BADFRAC"); // Maximum bad fraction + float poorFrac = psMetadataLookupF32(&mdok, recipe, "POOR.FRACTION"); // Fraction for "poor" + + // XXX EAM : do we need to / want to define different values for BAD and POOR subtraction vs BAD and POOR warp? + psImageMaskType maskVal = pmConfigMaskGet("MASK.VALUE", config); // Bits to mask going in to pmSubtractionMatch + psImageMaskType maskPoor = pmConfigMaskGet("POOR.WARP", config); // Bits to mask for poor pixels + psImageMaskType maskBad = pmConfigMaskGet("BLANK", config); // Bits to mask for bad pixels + + bool dual = psMetadataLookupBool(&mdok, recipe, "DUAL"); // Dual convolution? + pmSubtractionMode subMode = dual ? PM_SUBTRACTION_MODE_DUAL : PM_SUBTRACTION_MODE_UNSURE; // Subtracn mode + + int threads = psMetadataLookupS32(NULL, config->arguments, "-threads"); // Number of threads + if (threads > 0) { + pmSubtractionThreadsInit(inRO, refRO); + } + + // 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."); + return false; + } + + psFree(optWidths); + + pmSubtractionThreadsFinalize(inRO, refRO); + + // XXX drop the pixels associated with inRO and refRO (now that we have inConv and refConf) + + psphotSaveImage (NULL, inRO->image, "inRO.fits"); + psphotSaveImage (NULL, refRO->image, "refRO.fits"); + psphotSaveImage (NULL, inConv->image, "inConv.fits"); + psphotSaveImage (NULL, refConv->image, "refConv.fits"); + + return true; +} 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; Index: trunk/ppSub/src/ppSubReadoutPhotometry.c =================================================================== --- trunk/ppSub/src/ppSubReadoutPhotometry.c (revision 21257) +++ trunk/ppSub/src/ppSubReadoutPhotometry.c (revision 21257) @@ -0,0 +1,145 @@ +#include "ppSub.h" + +// Photometry stage 2: find and measure sources on the subtracted image +bool ppSubReadoutPhotometry (pmConfig *config, psMetadata *stats, const pmFPAview *view) { + + bool mdok = false; + + // 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, "PHOTOMETRY")) { + pmFPAfileDropInternal(config->files, "PSPHOT.BACKMDL"); + pmFPAfileDropInternal(config->files, "PSPHOT.BACKMDL.STDEV"); + return true; + } + + // The PSF (measured in ppSubMakePSF) is stored on the chip->analysis of PSPHOT.INPUT + // In order to use an incoming PSF, it must be stored on the chip->analysis of PSPHOT.PSF.LOAD + pmChip *psfInputChip = pmFPAfileThisChip(config->files, view, "PSPHOT.INPUT"); + psAssert (psfInputChip, "should have been generated for ppSubMakePSF"); + + pmChip *psfLoadChip = pmFPAfileThisChip(config->files, view, "PSPHOT.PSF.LOAD"); + psAssert (psfLoadChip, "PSPHOT.PSF.LOAD should have been defined in ppSubCamera"); + + pmPSF *psf = psMetadataLookupPtr(NULL, psfInputChip->analysis, "PSPHOT.PSF"); + if (!psf) { + psError(PS_ERR_UNEXPECTED_NULL, true, "Unable to find PSF from psphot"); + return false; + } + psMetadataAddPtr(psfLoadChip->analysis, PS_LIST_TAIL, "PSPHOT.PSF", PS_DATA_UNKNOWN | PS_META_REPLACE, "PSF from matched addition", psf); + + // psphotReadoutMinimal performs the photometry analysis on PSPHOT.INPUT; we need to move + // around the pointers so PSPHOT.INPUT corresponds to the output image; previously, it was + // equivalent to the minuend image. + pmFPAfile *outputFile = psMetadataLookupPtr(&mdok, config->files, "PPSUB.OUTPUT"); + pmReadout *outRO = pmFPAviewThisReadout(view, outputFile->fpa); // Readout with the sources + + // XXX possibly rename this to PPSUB.RESID? + pmFPAfile *photFile = psMetadataLookupPtr(&mdok, config->files, "PSPHOT.RESID"); + pmReadout *photRO = pmFPAviewThisReadout(view, photFile->fpa); // Readout with the sources + if (!photRO) { + pmCell *cell = pmFPAfileThisCell(config->files, view, "PSPHOT.RESID"); // Output cell + photRO = pmReadoutAlloc(cell); // Output readout: subtraction + } + photRO->image = psImageCopy(photRO->image, outRO->image, PS_TYPE_F32); + if (outRO->weight) { + photRO->weight = psImageCopy(photRO->weight, outRO->weight, PS_TYPE_F32); + } else { + psFree (photRO->weight); + photRO->weight = NULL; + } + if (outRO->mask) { + photRO->mask = psImageCopy(photRO->mask, outRO->mask, PS_TYPE_IMAGE_MASK); + } else { + psFree (photRO->mask); + photRO->mask = NULL; + } + + // pmFPAfile *photFile = outputFile; + psMetadataAddPtr (config->files, PS_LIST_TAIL, "PSPHOT.INPUT", PS_DATA_UNKNOWN | PS_META_REPLACE, "psphot input : view on another pmFPAfile", photFile); + + // old-style variance renormalization + if (!ppSubReadoutRenormPhot (config, recipe, photRO)) { + psError(PS_ERR_UNKNOWN, false, "failure in renormalization"); + return false; + } + + if (!psphotReadoutMinimal(config, view)) { + psWarning("Unable to perform photometry on subtracted image."); + psErrorStackPrint(stderr, "Error stack from photometry:"); + psErrorClear(); + } + + if (stats) { + 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")); + } + + pmFPAfileDropInternal(config->files, "PSPHOT.BACKMDL"); + pmFPAfileDropInternal(config->files, "PSPHOT.BACKMDL.STDEV"); + + return true; +} + + // XXX not sure that this is still needed (only if psphotReadoutMinimal measures the background) + // 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"); + + +#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 Index: trunk/ppSub/src/ppSubReadoutRenorm.c =================================================================== --- trunk/ppSub/src/ppSubReadoutRenorm.c (revision 21257) +++ trunk/ppSub/src/ppSubReadoutRenorm.c (revision 21257) @@ -0,0 +1,51 @@ +#include "ppSub.h" + +bool ppSubReadoutRenormPixels (pmConfig *config, psMetadata *recipe, pmReadout *readout) { + + bool mdok = false; + + if (!psMetadataLookupBool(&mdok, recipe, "RENORM")) return true; + + // Statistics for renormalisation + psStatsOptions renormMean = psStatsOptionFromString(psMetadataLookupStr(&mdok, recipe, "RENORM.MEAN")); + psStatsOptions renormStdev = psStatsOptionFromString(psMetadataLookupStr(&mdok, 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; + } + psImageMaskType maskValue = pmConfigMaskGet("BLANK", config); // Bits to mask + if (!pmReadoutWeightRenormPixels(readout, maskValue, renormMean, renormStdev, NULL)) { + psError(PS_ERR_UNKNOWN, false, "Unable to renormalise variances."); + return false; + } + return true; +} + +// old-style variance renormalization +bool ppSubReadoutRenormPhot (pmConfig *config, psMetadata *recipe, pmReadout *readout) { + + bool mdok = false; + + if (!psMetadataLookupBool(&mdok, recipe, "RENORM")) return true; + + // Statistics for renormalisation + psStatsOptions renormMean = psStatsOptionFromString(psMetadataLookupStr(&mdok, recipe, "RENORM.MEAN")); + psStatsOptions renormStdev = psStatsOptionFromString(psMetadataLookupStr(&mdok, 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; + } + + // other renorm parameters + int renormNum = psMetadataLookupS32(&mdok, recipe, "RENORM.NUM"); // Number of samples for renormalisation + float renormWidth = psMetadataLookupS32(&mdok, recipe, "RENORM.WIDTH"); // Width of Gaussian phot + + psImageMaskType maskValue = pmConfigMaskGet("BLANK", config); // Bits to mask + if (!pmReadoutWeightRenormPhot(readout, maskValue, renormNum, renormWidth, renormMean, renormStdev, NULL)) { + psError(PS_ERR_UNKNOWN, false, "Unable to renormalise variances."); + return false; + } + return true; +} Index: trunk/ppSub/src/ppSubReadoutSubtract.c =================================================================== --- trunk/ppSub/src/ppSubReadoutSubtract.c (revision 21257) +++ trunk/ppSub/src/ppSubReadoutSubtract.c (revision 21257) @@ -0,0 +1,104 @@ +#include "ppSub.h" +#define WCS_TOLERANCE 0.001 // Tolerance for WCS + +bool ppSubReadoutSubtract (pmConfig *config, const pmFPAview *view) { + + bool mdok = false; + + // Look up recipe values + psMetadata *recipe = psMetadataLookupMetadata(&mdok, config->recipes, PPSUB_RECIPE); // Recipe for ppSim + psAssert(recipe, "We checked this earlier, so it should be here."); + + bool reverse = psMetadataLookupBool(&mdok, config->arguments, "REVERSE"); // Reverse sense of subtraction? + + // Subtraction is: minuend - subtrahend + pmReadout *minuend = NULL; + pmReadout *subtrahend = NULL; + if (reverse) { + minuend = pmFPAfileThisReadout(config->files, view, "PPSUB.REF.CONV"); + subtrahend = pmFPAfileThisReadout(config->files, view, "PPSUB.INPUT.CONV"); + } else { + minuend = pmFPAfileThisReadout(config->files, view, "PPSUB.INPUT.CONV"); + subtrahend = pmFPAfileThisReadout(config->files, view, "PPSUB.REF.CONV"); + } + +// XXX this should be replaced by SAVE options to the inConv and refConf pmFPAfiles +#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 + + // Do the actual subtraction + pmReadout *outRO = pmFPAfileThisReadout(config->files, view, "PPSUB.OUTPUT"); + outRO->image = (psImage *) psBinaryOp(outRO->image, minuend->image, "-", subtrahend->image); + if (minuend->weight && subtrahend->weight) { + outRO->weight = (psImage *) psBinaryOp(outRO->weight, minuend->weight, "+", subtrahend->weight); + } + outRO->mask = (psImage *) psBinaryOp(outRO->mask, minuend->mask, "|", subtrahend->mask); + + outRO->data_exists = true; + outRO->parent->data_exists = true; + outRO->parent->parent->data_exists = true; + + // copy concepts from the input to the output (XXX should this always use minuend?) + pmFPAfile *inFile = psMetadataLookupPtr (&mdok, config->files, "PPSUB.INPUT"); + pmFPA *inFPA = inFile->fpa; + + pmFPAfile *outFile = psMetadataLookupPtr (&mdok, config->files, "PPSUB.OUTPUT"); + pmFPA *outFPA = outFile->fpa; + + if (!pmFPACopyConcepts(outFPA, inFPA)) { + psError(PS_ERR_UNKNOWN, false, "Unable to copy concepts from input to output."); + psFree(outRO); + return false; + } + + // get the HDUs and output Chip to copy the astrometry + pmHDU *inHDU = inFPA->hdu; // input HDU + pmHDU *outHDU = outFPA->hdu; + pmChip *outChip = pmFPAfileThisChip(config->files, view, "PPSUB.OUTPUT"); + + if (!outHDU || !inHDU) { + psError(PS_ERR_UNKNOWN, false, "Unable to find HDU at FPA level to copy astrometry."); + return false; + } + + // Copy astrometry over + // It should get into the output images and photometry + if (!pmAstromReadWCS(outFPA, outChip, inHDU->header, 1.0)) { + psError(PS_ERR_UNKNOWN, false, "Unable to read WCS astrometry from input FPA."); + return false; + } + + if (!pmAstromWriteWCS(outHDU->header, outFPA, outChip, WCS_TOLERANCE)) { + psError(PS_ERR_UNKNOWN, false, "Unable to write WCS astrometry to output FPA."); + return false; + } + + return true; +} + +// XXX this test code was in place to check for and squash unexpected NANs + +#if 0 + pmReadoutMaskApply(outRO, maskBad); + + 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 + Index: trunk/ppSub/src/ppSubReadoutUpdate.c =================================================================== --- trunk/ppSub/src/ppSubReadoutUpdate.c (revision 21257) +++ trunk/ppSub/src/ppSubReadoutUpdate.c (revision 21257) @@ -0,0 +1,76 @@ +#include "ppSub.h" + +// Renormalize, update headers and generate JPEGs +bool ppSubReadoutUpdate (pmConfig *config, const pmFPAview *view) { + + bool mdok = false; + + // Look up recipe values + psMetadata *recipe = psMetadataLookupMetadata(&mdok, config->recipes, PPSUB_RECIPE); // Recipe for ppSim + psAssert(recipe, "We checked this earlier, so it should be here."); + + // select the output readout + pmReadout *outRO = pmFPAfileThisReadout(config->files, view, "PPSUB.OUTPUT"); + + // Renormalising for pixels, because that's what magic desires + if (!ppSubReadoutRenormPixels (config, recipe, outRO)) { + psError(PS_ERR_BAD_PARAMETER_VALUE, false, "failure in renormalization"); + return false; + } + + // select the output FPA and HDU to get the output header + pmFPAfile *outFile = psMetadataLookupPtr (&mdok, config->files, "PPSUB.OUTPUT"); + pmFPA *outFPA = outFile->fpa; + pmHDU *outHDU = outFPA->hdu; + + // 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 + { + psImageMaskType maskBad = pmConfigMaskGet("BLANK", config); // Bits to mask for bad pixels + + 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); + pmReadout *ro2 = pmReadoutAlloc(cell2); // Binned readouts + if (!pmReadoutRebin(ro1, outRO, maskBad, bin1, bin1)) { + psError(PS_ERR_UNKNOWN, false, "Unable to bin output (1st binning)"); + psFree(ro1); + psFree(ro2); + return false; + } + if (!pmReadoutRebin(ro2, ro1, 0, bin2, bin2)) { + psError(PS_ERR_UNKNOWN, false, "Unable to bin output (2nd binning)"); + psFree(ro1); + psFree(ro2); + 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 + + return true; +} Index: trunk/ppSub/src/ppSubSetMasks.c =================================================================== --- trunk/ppSub/src/ppSubSetMasks.c (revision 21257) +++ trunk/ppSub/src/ppSubSetMasks.c (revision 21257) @@ -0,0 +1,102 @@ +#include "ppSub.h" + +// this function generates (if needed) and sets or updates the masks for both the input and +// reference images. this function also has the code for interpolation over bad pixels, but it +// is currently if-def-ed out + +bool ppSubSetMasks (pmConfig *config, const pmFPAview *view) { + + psImageMaskType maskValue; + psImageMaskType markValue; + bool mdok = false; + + if (!pmConfigMaskSetBits(&maskValue, &markValue, config)) { + psError(PS_ERR_UNKNOWN, false, "Unable to determine mask value."); + return false; + } + + // set the mask bits needed by psphot (in psphot recipe) + psphotSetMaskRecipe (config, maskValue, markValue); + + // Look up recipe values + psMetadata *recipe = psMetadataLookupMetadata(&mdok, config->recipes, PPSUB_RECIPE); // Recipe for ppSim + psAssert(recipe, "We checked this earlier, so it should be here."); + + psImageMaskType satValue = pmConfigMaskGet("SAT", config); + psAssert (satValue, "SAT must be non-zero"); + + psImageMaskType badValue = pmConfigMaskGet("BAD", config); + psAssert (badValue, "BAD must be non-zero"); + + // input images + pmReadout *inRO = pmFPAfileThisReadout(config->files, view, "PPSUB.INPUT"); // Input readout + pmReadout *refRO = pmFPAfileThisReadout(config->files, view, "PPSUB.REF"); // Reference readout + + psImage *input = inRO->image; // Input image + psImage *reference = refRO->image; // Reference image + PS_ASSERT_IMAGES_SIZE_EQUAL(input, reference, false); + + // XXX use psImageCopy below to avoid caring about image dimensions + int numCols = input->numCols; + int numRows = input->numRows; + + // Generate masks if they don't exist + if (!inRO->mask) { + if (psMetadataLookupBool(NULL, recipe, "MASK.GENERATE")) { + pmReadoutSetMask(inRO, satValue, badValue); + } 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, satValue, badValue); + } else { + refRO->mask = psImageAlloc(numCols, numRows, PS_TYPE_IMAGE_MASK); + psImageInit(refRO->mask, 0); + } + } + + // Mask the NAN values + if (!pmReadoutMaskNonfinite(inRO, satValue)) { + psError(PS_ERR_UNKNOWN, false, "Unable to mask non-finite pixels in input."); + return false; + } + if (!pmReadoutMaskNonfinite(refRO, satValue)) { + psError(PS_ERR_UNKNOWN, false, "Unable to mask non-finite pixels in reference."); + return false; + } + +#if (0) + // 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."); + + bool mdok = false; + + 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" + + psImageMaskType maskPoor = pmConfigMaskGet("POOR.WARP", config); // Bits to mask for poor pixels + psImageMaskType maskBad = pmConfigMaskGet("BLANK", config); // Bits to mask for bad pixels + + // 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 + + return true; +} Index: trunk/ppSub/src/ppSubVarianceFactors.c =================================================================== --- trunk/ppSub/src/ppSubVarianceFactors.c (revision 21257) +++ trunk/ppSub/src/ppSubVarianceFactors.c (revision 21257) @@ -0,0 +1,103 @@ +#include "ppSub.h" + +// Calculate the variance factor for the output image based on the input images +bool ppSubVarianceFactors (pmConfig *config, psMetadata *stats, const pmFPAview *view) { + + // 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."); + + // convolved input images + pmCell *inCell = pmFPAfileThisCell(config->files, view, "PPSUB.INPUT.CONV"); // Input cell + pmCell *refCell = pmFPAfileThisCell(config->files, view, "PPSUB.REF.CONV"); // Reference cell + + pmReadout *inConv = pmFPAfileThisReadout(config->files, view, "PPSUB.INPUT.CONV"); // Input cell + pmReadout *refConv = pmFPAfileThisReadout(config->files, view, "PPSUB.REF.CONV"); // Reference cell + + pmReadout *outRO = pmFPAfileThisReadout(config->files, view, "PPSUB.OUTPUT"); // Output readout + + // 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."); + return false; + } + + float vfIn = psMetadataLookupF32(NULL, inCell->concepts, "CELL.VARFACTOR"); // Variance factor for input + if (!isfinite(vfIn)) { + vfIn = 1.0; + } + float vfRef = psMetadataLookupF32(NULL, refCell->concepts, "CELL.VARFACTOR"); // Variance factor for ref + if (!isfinite(vfRef)) { + vfRef = 1.0; + } + + if (kernels->mode == PM_SUBTRACTION_MODE_1) { + vfIn *= pmSubtractionVarianceFactor(kernels, 0.0, 0.0, false); + } + if (kernels->mode == PM_SUBTRACTION_MODE_2) { + vfRef *= pmSubtractionVarianceFactor(kernels, 0.0, 0.0, false); + } + if (kernels->mode == PM_SUBTRACTION_MODE_DUAL) { + vfIn *= pmSubtractionVarianceFactor(kernels, 0.0, 0.0, false); + 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 + + // mask these values when examining the background + psImageMaskType maskVal = pmConfigMaskGet("MASK.VALUE", config); // Bits to mask going in to pmSubtractionMatch + psImageMaskType maskBad = pmConfigMaskGet("BLANK", config); // Bits to mask for bad pixels + + // measure the mean of the convolved input variance image + 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(vfStats); + psFree(rng); + return false; + } + float inMeanVar = vfStats->robustMedian; // Mean variance of input + + // measure the mean of the convolved reference variance image + 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(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) { + psMetadataCopySingle(stats, outRO->analysis, PM_SUBTRACTION_ANALYSIS_MODE); + psMetadataCopySingle(stats, outRO->analysis, PM_SUBTRACTION_ANALYSIS_STAMPS); + psMetadataCopySingle(stats, outRO->analysis, PM_SUBTRACTION_ANALYSIS_DEV_MEAN); + psMetadataCopySingle(stats, outRO->analysis, PM_SUBTRACTION_ANALYSIS_DEV_RMS); + psMetadataCopySingle(stats, outRO->analysis, PM_SUBTRACTION_ANALYSIS_NORM); + psMetadataCopySingle(stats, outRO->analysis, PM_SUBTRACTION_ANALYSIS_BGDIFF); + psMetadataCopySingle(stats, outRO->analysis, PM_SUBTRACTION_ANALYSIS_MX); + psMetadataCopySingle(stats, outRO->analysis, PM_SUBTRACTION_ANALYSIS_MY); + psMetadataCopySingle(stats, outRO->analysis, PM_SUBTRACTION_ANALYSIS_MXX); + psMetadataCopySingle(stats, outRO->analysis, PM_SUBTRACTION_ANALYSIS_MXY); + psMetadataCopySingle(stats, outRO->analysis, PM_SUBTRACTION_ANALYSIS_MYY); + + psMetadataAddF32(stats, PS_LIST_TAIL, "TIME_MATCH", 0, "Time to match PSFs", + psTimerClear("PPSUB_MATCH")); + } + + return true; +} Index: trunk/ppSub/src/ppSubVersion.c =================================================================== --- trunk/ppSub/src/ppSubVersion.c (revision 21183) +++ trunk/ppSub/src/ppSubVersion.c (revision 21257) @@ -1,11 +1,2 @@ -#ifdef HAVE_CONFIG_H -#include -#endif - -#include -#include -#include -#include - #include "ppSub.h"