Index: trunk/ppSub/doc/outline.txt =================================================================== --- trunk/ppSub/doc/outline.txt (revision 26899) +++ trunk/ppSub/doc/outline.txt (revision 26899) @@ -0,0 +1,26 @@ + +bool ppSubMatchPSFs(ppSubData *data) + - setup huge number of recipe values + - load or find sources + - pmSubtractionMatch + - data normalization + +bool pmSubtractionMatch() + - setup + - loop over iso-kernel regions + - identify stamp sources + - extract stamp pixels + - generate kernel basis functions + - iterate to reject stamps + - pmSubtractionCalculateEquation + - pmSubtractionSolveEquation + - pmSubtractionCalculateDeviations + - pmSubtractionRejectStamps + +bool pmSubtractionCalculateEquationStamp() + - convolve stamp with each kernel basis function + - generate the polynomial terms for the basis functions + - sum the elements of the chisq equation + +bool pmSubtractionSolveEquation() + Index: trunk/ppSub/src/ppSub.c =================================================================== --- trunk/ppSub/src/ppSub.c (revision 26571) +++ trunk/ppSub/src/ppSub.c (revision 26899) @@ -24,4 +24,15 @@ int main(int argc, char *argv[]) { + +# if 0 + psLibInit(NULL); + pmVisualSetVisual(true); + for (int order = 2; order < 11; order ++) { + pmSubtractionDeconvolutionTest (order); + } + psLibFinalize(); + exit (1); +# endif + psExit exitValue = PS_EXIT_SUCCESS; // Exit value psTimerStart("ppSub"); Index: trunk/ppSub/src/ppSub.h =================================================================== --- trunk/ppSub/src/ppSub.h (revision 26571) +++ trunk/ppSub/src/ppSub.h (revision 26899) @@ -152,4 +152,7 @@ ); +/// Generate JPEG images +bool ppSubResidualSampleJpeg(pmConfig *config); + /// Generate inverse subtraction bool ppSubReadoutInverse(pmConfig *config // Configuration @@ -160,4 +163,6 @@ bool psMetadataCopySingle(psMetadata *target, psMetadata *source, const char *name); +bool ppSubCopyPSF (pmFPAfile *output, pmFPAfile *input, pmFPAview *view); + ///@} #endif Index: trunk/ppSub/src/ppSubBackground.c =================================================================== --- trunk/ppSub/src/ppSubBackground.c (revision 26571) +++ trunk/ppSub/src/ppSubBackground.c (revision 26899) @@ -42,5 +42,5 @@ if (!modelRO) { // Create the background model - if (!psphotModelBackground(config, view, "PPSUB.OUTPUT")) { + if (!psphotModelBackgroundReadoutFileIndex(config, view, "PPSUB.OUTPUT", 0)) { psError(PS_ERR_UNKNOWN, false, "Unable to model background"); psFree(view); Index: trunk/ppSub/src/ppSubCamera.c =================================================================== --- trunk/ppSub/src/ppSubCamera.c (revision 26571) +++ trunk/ppSub/src/ppSubCamera.c (revision 26899) @@ -23,5 +23,6 @@ // Define an input file -static pmFPAfile *defineInputFile(pmConfig *config,// Configuration +static pmFPAfile *defineInputFile(bool *success, + pmConfig *config,// Configuration pmFPAfile *bind, // File to which to bind, or NULL char *filerule, // Name of file rule @@ -32,5 +33,8 @@ bool status; + *success = false; + pmFPAfile *file = NULL; + // look for the file on the argument list if (bind) { @@ -39,4 +43,5 @@ file = pmFPAfileDefineFromArgs(&status, config, filerule, argname); } + if (!status) { psError(PS_ERR_UNKNOWN, false, "Failed to load file definition for %s", filerule); @@ -53,4 +58,6 @@ if (!file) { + // no file defined + *success = true; return NULL; } @@ -61,4 +68,5 @@ } + *success = true; return file; } @@ -137,4 +145,6 @@ bool ppSubCamera(ppSubData *data) { + bool success = true; + psAssert(data, "Require processing data"); pmConfig *config = data->config; @@ -142,25 +152,52 @@ // Input image - pmFPAfile *input = defineInputFile(config, NULL, "PPSUB.INPUT", "INPUT", PM_FPA_FILE_IMAGE); - if (!input) { + pmFPAfile *input = defineInputFile(&success, config, NULL, "PPSUB.INPUT", "INPUT", PM_FPA_FILE_IMAGE); + if (!success) { psError(PS_ERR_IO, false, "Failed to build FPA from PPSUB.INPUT"); return false; } - defineInputFile(config, input, "PPSUB.INPUT.MASK", "INPUT.MASK", PM_FPA_FILE_MASK); - pmFPAfile *inVar = defineInputFile(config, input, "PPSUB.INPUT.VARIANCE", "INPUT.VARIANCE", - PM_FPA_FILE_VARIANCE); - defineInputFile(config, NULL, "PPSUB.INPUT.SOURCES", "INPUT.SOURCES", PM_FPA_FILE_CMF); + + defineInputFile(&success, config, input, "PPSUB.INPUT.MASK", "INPUT.MASK", PM_FPA_FILE_MASK); + if (!success) { + psError(PS_ERR_IO, false, "Failed to build FPA from PPSUB.INPUT.MASK"); + return false; + } + + pmFPAfile *inVar = defineInputFile(&success, config, input, "PPSUB.INPUT.VARIANCE", "INPUT.VARIANCE", PM_FPA_FILE_VARIANCE); + if (!success) { + psError(PS_ERR_IO, false, "Failed to build FPA from PPSUB.INPUT.VARIANCE"); + return false; + } + + defineInputFile(&success, config, NULL, "PPSUB.INPUT.SOURCES", "INPUT.SOURCES", PM_FPA_FILE_CMF); + if (!success) { + psError(PS_ERR_IO, false, "Failed to build FPA from PPSUB.INPUT.SOURCES"); + return false; + } // Reference image - pmFPAfile *ref = defineInputFile(config, NULL, "PPSUB.REF", "REF", PM_FPA_FILE_IMAGE); - if (!ref) { + pmFPAfile *ref = defineInputFile(&success, config, NULL, "PPSUB.REF", "REF", PM_FPA_FILE_IMAGE); + if (!success) { psError(PS_ERR_IO, false, "Failed to build FPA from PPSUB.REF"); return false; } - defineInputFile(config, ref, "PPSUB.REF.MASK", "REF.MASK", PM_FPA_FILE_MASK); - pmFPAfile *refVar = defineInputFile(config, ref, "PPSUB.REF.VARIANCE", "REF.VARIANCE", - PM_FPA_FILE_VARIANCE); - defineInputFile(config, NULL, "PPSUB.REF.SOURCES", "REF.SOURCES", PM_FPA_FILE_CMF); - + + defineInputFile(&success, config, ref, "PPSUB.REF.MASK", "REF.MASK", PM_FPA_FILE_MASK); + if (!success) { + psError(PS_ERR_IO, false, "Failed to build FPA from PPSUB.REF.MASK"); + return false; + } + + pmFPAfile *refVar = defineInputFile(&success, config, ref, "PPSUB.REF.VARIANCE", "REF.VARIANCE", PM_FPA_FILE_VARIANCE); + if (!success) { + psError(PS_ERR_IO, false, "Failed to build FPA from PPSUB.REF.VARIANCE"); + return false; + } + + defineInputFile(&success, config, NULL, "PPSUB.REF.SOURCES", "REF.SOURCES", PM_FPA_FILE_CMF); + if (!success) { + psError(PS_ERR_IO, false, "Failed to build FPA from PPSUB.REF.SOURCES"); + return false; + } // Now that the camera has been determined, we can read the recipe @@ -298,4 +335,16 @@ jpeg2->save = true; + // Output residual JPEG + pmFPAfile *jpeg3 = pmFPAfileDefineOutput(config, NULL, "PPSUB.OUTPUT.RESID.JPEG"); + if (!jpeg3) { + psError(PS_ERR_IO, false, _("Unable to generate output file from PPSUB.OUTPUT.RESID.JPEG")); + return false; + } + if (jpeg3->type != PM_FPA_FILE_JPEG) { + psError(PS_ERR_IO, true, "PPSUB.OUTPUT.RESID.JPEG is not of type JPEG"); + return false; + } + jpeg3->save = true; + // Output subtraction kernel pmFPAfile *kernel = defineCalcFile(config, NULL, "PPSUB.OUTPUT.KERNELS", PM_FPA_FILE_SUBKERNEL); @@ -306,5 +355,5 @@ // psPhot input - if (data->photometry) { + if (data->photometry || 1) { psphotModelClassInit(); // load implementation-specific models @@ -318,8 +367,10 @@ return false; } + // specify the number of psphot input images + psMetadataAddS32 (config->arguments, PS_LIST_TAIL, "PSPHOT.INPUT.NUM", PS_META_REPLACE, "number of inputs", 1); pmFPAfileActivate(config->files, false, "PSPHOT.INPUT"); - // Internal-ish file for getting the PSF from the minuend - pmFPAfile *psf = pmFPAfileDefineOutputFromFile(config, psphot, "PSPHOT.PSF.LOAD"); + // Internal file for getting the PSF from the minuend + pmFPAfile *psf = pmFPAfileDefineFromFPA(config, output->fpa, 1, 1, "PSPHOT.PSF.LOAD"); if (!psf) { psError(PS_ERR_IO, false, "Failed to build FPA from PSPHOT.PSF.LOAD"); Index: trunk/ppSub/src/ppSubDefineOutput.c =================================================================== --- trunk/ppSub/src/ppSubDefineOutput.c (revision 26571) +++ trunk/ppSub/src/ppSubDefineOutput.c (revision 26899) @@ -51,7 +51,9 @@ bool mdok; // Status of MD lookup psMetadata *analysis = inConv->analysis; // Analysis metadata with kernel information - pmSubtractionKernels *kernels = psMetadataLookupPtr(&mdok, analysis, - PM_SUBTRACTION_ANALYSIS_KERNEL); // Subtraction kernel + pmHDU *hdu = pmHDUFromCell(inConv->parent); + pmSubtractionKernels *kernels = psMetadataLookupPtr(&mdok, analysis, PM_SUBTRACTION_ANALYSIS_KERNEL); // Subtraction kernel + if (!kernels) { + hdu = pmHDUFromCell(refConv->parent); analysis = refConv->analysis; kernels = psMetadataLookupPtr(&mdok, analysis, PM_SUBTRACTION_ANALYSIS_KERNEL); @@ -62,6 +64,7 @@ return false; } - psMetadataAddStr(outHDU->header, PS_LIST_TAIL, "PPSUB.KERNEL", 0, "Subtraction kernel", - kernels->description); + psAssert (hdu, "unable to find HDU"); + psMetadataAddStr(outHDU->header, PS_LIST_TAIL, "PPSUB.KERNEL", 0, "Subtraction kernel", kernels->description); + outHDU->header = psMetadataCopy(outHDU->header, hdu->header); // Add additional data to the header Index: trunk/ppSub/src/ppSubFiles.c =================================================================== --- trunk/ppSub/src/ppSubFiles.c (revision 26571) +++ trunk/ppSub/src/ppSubFiles.c (revision 26899) @@ -23,5 +23,5 @@ // Subtraction files static const char *subFiles[] = { "PPSUB.OUTPUT", "PPSUB.OUTPUT.MASK", "PPSUB.OUTPUT.VARIANCE", - "PPSUB.OUTPUT.JPEG1", "PPSUB.OUTPUT.JPEG2", + "PPSUB.OUTPUT.JPEG1", "PPSUB.OUTPUT.JPEG2", "PPSUB.OUTPUT.RESID.JPEG", NULL }; Index: trunk/ppSub/src/ppSubLoop.c =================================================================== --- trunk/ppSub/src/ppSubLoop.c (revision 26571) +++ trunk/ppSub/src/ppSubLoop.c (revision 26899) @@ -27,5 +27,5 @@ pmConfigCamerasCull(config, NULL); - pmConfigRecipesCull(config, "PPSUB,PPSTATS,PSPHOT,MASKS,JPEG"); + pmConfigRecipesCull(config, "PPSUB,PPSTATS,PSPHOT,PSASTRO,MASKS,JPEG"); @@ -54,4 +54,9 @@ // Can't do anything at all return true; + } + // generate the residual stamp grid for visualization + if (!ppSubResidualSampleJpeg(config)) { + psError(PS_ERR_UNKNOWN, false, "Unable to update."); + return false; } @@ -130,4 +135,5 @@ } + // generate the binned image used to write the jpeg if (!ppSubReadoutJpeg(config)) { psError(PS_ERR_UNKNOWN, false, "Unable to update."); Index: trunk/ppSub/src/ppSubMakePSF.c =================================================================== --- trunk/ppSub/src/ppSubMakePSF.c (revision 26571) +++ trunk/ppSub/src/ppSubMakePSF.c (revision 26899) @@ -22,4 +22,6 @@ #include "ppSub.h" +psArray *ppSubSelectPSFSources(psArray *sources); + bool ppSubMakePSF(ppSubData *data) { @@ -61,20 +63,33 @@ return false; } + pmReadout *photRO = pmFPAviewThisReadout(view, photFile->fpa); // Readout to photometer - if (psMetadataLookup(photRO->analysis, "PSPHOT.SOURCES")) { - psMetadataRemoveKey(photRO->analysis, "PSPHOT.SOURCES"); - } - - // psphotSaveImage (photFile->fpa->hdu->header, photRO->image, "findpsf.im.fits"); - // psphotSaveImage (photFile->fpa->hdu->header, photRO->variance, "findpsf.wt.fits"); - // psphotSaveImage (photFile->fpa->hdu->header, photRO->mask, "findpsf.mk.fits"); - - // XXX for testing, can i dump a complete copy of the psphot sources here? - // I actually only need the X,Y coordinates to test this. + + // we need to remove any existing PSPHOT.DETECTIONS (why not do this in psphot?) + if (psMetadataLookup(photRO->analysis, "PSPHOT.DETECTIONS")) { + psMetadataRemoveKey(photRO->analysis, "PSPHOT.DETECTIONS"); + } + if (psMetadataLookup(photRO->parent->parent->analysis, "PSPHOT.PSF")) { + psMetadataRemoveKey(photRO->parent->parent->analysis, "PSPHOT.PSF"); + } + +# ifdef TESTING + // XXX for testing, dump these images: + psphotSaveImage (NULL, photRO->image, "findpsf.im.fits"); + psphotSaveImage (NULL, photRO->variance, "findpsf.wt.fits"); + psphotSaveImage (NULL, photRO->mask, "findpsf.mk.fits"); +# endif // Extract the loaded sources from the associated readout, and generate PSF // Here, we assume the image is background-subtracted - psArray *sources = psMetadataLookupPtr(&mdok, minuend->analysis, "PSPHOT.SOURCES"); - if (!psphotReadoutFindPSF(config, view, sources)) { + pmDetections *detections = psMetadataLookupPtr(&mdok, minuend->analysis, "PSPHOT.DETECTIONS"); + psAssert(detections, "missing detections?"); + psArray *sources = detections->allSources; + psAssert(sources, "missing sources?"); + + // XXX filter sources? limit the total number and return only brighter objects? + // use flags to toss totally bogus entries? + psArray *goodSources = ppSubSelectPSFSources (sources); + if (!psphotReadoutFindPSF(config, view, goodSources)) { // This is likely a data quality issue // XXX Split into multiple cases using error codes? @@ -83,5 +98,14 @@ ppSubDataQuality(data, PSPHOT_ERR_PSF, PPSUB_FILES_PHOT_SUB | PPSUB_FILES_PHOT_INV); psFree(view); + psFree(goodSources); return true; + } + + // save the resulting PSF information on the pmFPAfile PSPHOT.PSF.LOAD + pmFPAfile *psfFile = psMetadataLookupPtr(&mdok, config->files, "PSPHOT.PSF.LOAD"); // PSF file + if (!ppSubCopyPSF (psfFile, photFile, view)) { + psErrorStackPrint(stderr, "PSF was not generated"); + psWarning("PSF was not generated --- suspect bad data quality."); + ppSubDataQuality(data, psErrorCodeLast(), PPSUB_FILES_PHOT_SUB | PPSUB_FILES_PHOT_INV); } @@ -89,11 +113,88 @@ psMetadataRemoveKey(photRO->analysis, "PSPHOT.HEADER"); - if (!data->quality) { - data->psf = psMemIncrRefCounter(psMetadataLookupPtr(NULL, photRO->parent->parent->analysis, - "PSPHOT.PSF")); - } - + psFree(goodSources); psFree(view); return true; } + +bool ppSubCopyPSF (pmFPAfile *output, pmFPAfile *input, pmFPAview *view) { + + bool mdok = false; + + pmChip *inputChip = pmFPAviewThisChip(view, input->fpa); // Chip with PSF info + pmChip *outputChip = pmFPAviewThisChip(view, output->fpa); // Chip to store PSF info + + pmReadout *inputRO = pmFPAviewThisReadout(view, input->fpa); // Readout with PSF info + pmReadout *outputRO = pmFPAviewThisReadout(view, output->fpa); // Readout to store PSF info + + if (!outputRO) { + pmCell *outputCell = pmFPAviewThisCell(view, output->fpa); + outputRO = pmReadoutAlloc(outputCell); + outputRO->image = psMemIncrRefCounter (inputRO->image); + } + + // copy the PSF-related data to PSPHOT.PSF.LOAD for safe-keeping + psMetadata *psfRegions = psMetadataAlloc(); + + int nRegions = psMetadataLookupS32 (&mdok, inputRO->analysis, "PSF.CLUMP.NREGIONS"); + if (!nRegions) { + psErrorStackPrint(stderr, "No PSF available"); + return false; + } + psMetadataAddS32 (psfRegions, PS_LIST_TAIL, "PSF.CLUMP.NREGIONS", PS_META_REPLACE, "psf clump regions", nRegions); + + for (int i = 0; i < nRegions; i++) { + char fieldName[80]; + snprintf (fieldName, 80, "PSF.CLUMP.REGION.%03d", i); + psMetadata *regionMD = psMetadataLookupPtr (&mdok, inputRO->analysis, fieldName); + if (!regionMD) { + psWarning ("missing psf clump region metadata for entry %d", i); + continue; + } + psMetadataAddMetadata (psfRegions, PS_LIST_TAIL, fieldName, PS_META_REPLACE, "psf clump region", regionMD); + } + + // XXX why am I replacing the entire analysis MD? + psFree(outputRO->analysis); + outputRO->analysis = psfRegions; + + // copy the PSF model data + pmPSF *psf = psMetadataLookupPtr(NULL, inputChip->analysis, "PSPHOT.PSF"); // PSF for photometry + if (!psf) { + psErrorStackPrint(stderr, "No PSF available"); + return false; + } + + psMetadataAddPtr(outputChip->analysis, PS_LIST_TAIL, "PSPHOT.PSF", PS_DATA_UNKNOWN | PS_META_REPLACE, "PSF from ppSubMakePSF", psf); + + return true; +} + + +// XXX hardwired MAX for now +# define MAX_NPSF 500 + +psArray *ppSubSelectPSFSources(psArray *sources){ + + sources = psArraySort (sources, pmSourceSortBySN); + + psArray *subset = psArrayAllocEmpty(MAX_NPSF); + + int nPSF = 0; + for (int i = 0; (nPSF < MAX_NPSF) && (i < sources->n); i++) { + + pmSource *source = sources->data[i]; + if (!source) continue; + + // skip non-astronomical objects (very likely defects) + if (source->type == PM_SOURCE_TYPE_DEFECT) continue; + if (source->type == PM_SOURCE_TYPE_SATURATED) continue; + if (source->mode & PM_SOURCE_MODE_SATSTAR) continue; + + psArrayAdd (subset, 100, source); + nPSF++; + } + + return subset; +} Index: trunk/ppSub/src/ppSubMatchPSFs.c =================================================================== --- trunk/ppSub/src/ppSubMatchPSFs.c (revision 26571) +++ trunk/ppSub/src/ppSubMatchPSFs.c (revision 26899) @@ -18,6 +18,9 @@ #include #include +#include #include "ppSub.h" + +#define COVAR_FRAC 0.01 // Truncation fraction for covariance matrix // Normalise a region on an image @@ -37,4 +40,133 @@ } +// Measure the PSF for an image +static float subImagePSF(ppSubData *data, // Processing data + const pmReadout *ro, // Readout for which to measure PSF + psArray *sources // Sources with positions at which to measure PSF + ) +{ + psAssert(data, "Require processing data"); + pmConfig *config = data->config; // Configuration + psAssert(config, "Require configuration"); + + psAssert(ro, "Need readout"); + psAssert(sources, "Need sources."); + + pmFPAfile *photFile = psMetadataLookupPtr(NULL, config->files, "PSPHOT.INPUT"); // Photometry file + psAssert(photFile, "Need photometry file."); + if (!pmFPACopy(photFile->fpa, ro->parent->parent->parent)) { + psError(PS_ERR_UNKNOWN, false, "Unable to copy FPA for photometry"); + return false; + } + + pmFPAview *view = ppSubViewReadout(); // View to readout + pmReadout *photRO = pmFPAviewThisReadout(view, photFile->fpa); // Readout to photometer + + if (psMetadataLookup(photRO->analysis, "PSPHOT.DETECTIONS")) { + psMetadataRemoveKey(photRO->analysis, "PSPHOT.DETECTIONS"); + } + if (psMetadataLookup(photRO->parent->parent->analysis, "PSPHOT.PSF")) { + psMetadataRemoveKey(photRO->parent->parent->analysis, "PSPHOT.PSF"); + } + + // Extract the loaded sources from the associated readout, and generate PSF + // Here, we assume the image is background-subtracted + if (!psphotReadoutFindPSF(config, view, sources)) { + psErrorStackPrint(stderr, "Unable to determine PSF"); + psWarning("Unable to determine PSF."); + psFree(view); + return true; + } + + psFree(view); + + psMetadata *header = psMetadataLookupMetadata(NULL, photRO->analysis, "PSPHOT.HEADER"); + psAssert(header, "Require header."); + float fwhm = psMetadataLookupF32(NULL, header, "FWHM_MAJ"); + + return fwhm; +} + +// Scale the kernel parameters according to the PSFs +static bool subScaleKernel(ppSubData *data, // Processing data + psVector *kernelWidths, // Widths for kernel + int *kernelSize, // Size of kernel + int *stampSize // Size of stamps (footprint) + ) +{ + psAssert(data, "Require processing data"); + pmConfig *config = data->config; // Configuration + psAssert(config, "Require configuration"); + + // Nothing to do if pre-calculated kernel exists + pmFPAview *view = ppSubViewReadout(); // View to readout + pmReadout *kernelRO = pmFPAfileThisReadout(config->files, view, "PPSUB.OUTPUT.KERNELS"); // RO with kernel + if (kernelRO) { + psFree(view); + return true; + } + + // Look up recipe values + psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, PPSUB_RECIPE); // Recipe for ppSim + psAssert(recipe, "We checked this earlier, so it should be here."); + if (!psMetadataLookupBool(NULL, recipe, "SCALE")) { + // No scaling requested + psFree(view); + return true; + } + + // Input images + pmReadout *inRO = pmFPAfileThisReadout(config->files, view, "PPSUB.INPUT"); // Input readout + pmReadout *refRO = pmFPAfileThisReadout(config->files, view, "PPSUB.REF"); // Reference readout + + // Input sources + pmReadout *inSourceRO = pmFPAfileThisReadout(config->files, view, "PPSUB.INPUT.SOURCES"); + pmReadout *refSourceRO = pmFPAfileThisReadout(config->files, view, "PPSUB.REF.SOURCES"); + // XXX assert on inSourcesRO and refSourcesRO? + + pmDetections *inDetections = psMetadataLookupPtr(NULL, inSourceRO->analysis, "PSPHOT.DETECTIONS"); // Input sources + pmDetections *refDetections = psMetadataLookupPtr(NULL, refSourceRO->analysis, "PSPHOT.DETECTIONS"); // Ref sources + + psFree(view); + + if (!inDetections || !refDetections) { + psWarning("Unable to scale kernel, since no sources were provided."); + return true; + } + + psArray *inSources = inDetections->allSources; + psAssert (inSources, "missing inSources?"); + + psArray *refSources = refDetections->allSources; + psAssert (refSources, "missing refSources?"); + + float inFWHM = subImagePSF(data, inRO, inSources); // FWHM for input + float refFWHM = subImagePSF(data, refRO, refSources); // FWHM for reference + + psLogMsg("ppSub", PS_LOG_INFO, "Input FWHM: %f\nReference FWHM: %f\n", inFWHM, refFWHM); + if (!isfinite(inFWHM) || !isfinite(refFWHM)) { + psWarning("Unable to scale kernel, since unable to measure PSFs."); + return true; + } + + float scaleRef = psMetadataLookupF32(NULL, recipe, "SCALE.REF"); // Reference for scaling + float scaleMin = psMetadataLookupF32(NULL, recipe, "SCALE.MIN"); // Minimum for scaling + float scaleMax = psMetadataLookupF32(NULL, recipe, "SCALE.MAX"); // Maximum for scaling + if (!isfinite(scaleRef) || !isfinite(scaleMin) || !isfinite(scaleMax)) { + psError(PPSUB_ERR_ARGUMENTS, false, + "Scale parameters (SCALE.REF=%f, SCALE.MIN=%f, SCALE.MAX=%f) not set in recipe.", + scaleRef, scaleMin, scaleMax); + return false; + } + + if (!pmSubtractionParamsScale(kernelSize, stampSize, kernelWidths, inFWHM, refFWHM, + scaleRef, scaleMin, scaleMax)) { + psError(PS_ERR_UNKNOWN, false, "Unable to scale parameters."); + return false; + } + + return true; +} + bool ppSubMatchPSFs(ppSubData *data) @@ -71,37 +203,44 @@ pmReadout *kernelRO = pmFPAfileThisReadout(config->files, view, "PPSUB.OUTPUT.KERNELS"); // RO with kernel - psFree(view); - // Sources in image, used for stamps: these must be loaded from previous analysis stages pmReadout *inSourceRO = pmFPAfileThisReadout(config->files, view, "PPSUB.INPUT.SOURCES"); pmReadout *refSourceRO = pmFPAfileThisReadout(config->files, view, "PPSUB.REF.SOURCES"); - psArray *inSources = inSourceRO ? psMetadataLookupPtr(&mdok, inSourceRO->analysis, "PSPHOT.SOURCES") : - NULL; // Source list from input image - psArray *refSources = refSourceRO ? psMetadataLookupPtr(&mdok, refSourceRO->analysis, "PSPHOT.SOURCES") : - NULL ; // Source list from reference image - - psArray *sources = NULL; // Merged list of sources - if (inSources && refSources) { + + pmDetections *inDetections = inSourceRO ? psMetadataLookupPtr(&mdok, inSourceRO->analysis, "PSPHOT.DETECTIONS") : NULL; // Input sources + pmDetections *refDetections = refSourceRO ? psMetadataLookupPtr(&mdok, refSourceRO->analysis, "PSPHOT.DETECTIONS") : NULL; // Ref sources + + psFree(view); + + pmDetections *detections = NULL; // Merged detection set + if (inDetections && refDetections) { + psArray *inSources = inDetections->allSources; + psArray *refSources = refDetections->allSources; + + psAssert (inSources, "missing in sources?"); + psAssert (refSources, "missing ref sources?"); + + detections = pmDetectionsAlloc(); float radius = psMetadataLookupF32(NULL, recipe, "SOURCE.RADIUS"); // Matching radius psArray *lists = psArrayAlloc(2); // Source lists lists->data[0] = psMemIncrRefCounter(inSources); lists->data[1] = psMemIncrRefCounter(refSources); - sources = pmSourceMatchMerge(lists, radius); + detections->allSources = pmSourceMatchMerge(lists, radius); psFree(lists); - if (!sources) { + if (!detections->allSources) { psError(PS_ERR_UNKNOWN, false, "Unable to merge source lists"); + psFree(detections); return false; } - } else if (inSources) { - sources = psMemIncrRefCounter(inSources); - } else if (refSources) { - sources = psMemIncrRefCounter(refSources); - } - - psMetadataAddArray(inConv->analysis, PS_LIST_TAIL, "PSPHOT.SOURCES", PS_META_REPLACE, - "Merged source list", sources); - psMetadataAddArray(refConv->analysis, PS_LIST_TAIL, "PSPHOT.SOURCES", PS_META_REPLACE, - "Merged source list", sources); - psFree(sources); // Drop reference + } + if (!detections && inDetections) { + detections = psMemIncrRefCounter(inDetections); + } + if (!detections && refDetections) { + detections = psMemIncrRefCounter(refDetections); + } + + psMetadataAddPtr(inConv->analysis, PS_LIST_TAIL, "PSPHOT.DETECTIONS", PS_META_REPLACE | PS_DATA_UNKNOWN, "Merged source list", detections); + psMetadataAddPtr(refConv->analysis, PS_LIST_TAIL, "PSPHOT.DETECTIONS", PS_META_REPLACE | PS_DATA_UNKNOWN, "Merged source list", detections); + psFree(detections); // Drop reference int footprint = psMetadataLookupS32(NULL, recipe, "STAMP.FOOTPRINT"); // Stamp half-size @@ -131,5 +270,8 @@ float rej = psMetadataLookupF32(NULL, recipe, "REJ"); // Rejection threshold float kernelErr = psMetadataLookupF32(NULL, recipe, "KERNEL.ERR"); // Relative systematic error in kernel + float normFrac = psMetadataLookupF32(NULL, recipe, "NORM.FRAC"); // Fraction of window for normalisn windw float sysErr = psMetadataLookupF32(NULL, recipe, "SYS.ERR"); // Relative systematic error in images + float skyErr = psMetadataLookupF32(NULL, recipe, "SKY.ERR"); // Additional error in sky + float covarFrac = psMetadataLookupF32(NULL, recipe, "COVAR.FRAC"); // Fraction for covariance calculation float badFrac = psMetadataLookupF32(NULL, recipe, "BADFRAC"); // Maximum bad fraction @@ -169,7 +311,12 @@ break; default: - psErrorStackPrint(stderr, "Invalid value for -convolve"); + psError(PS_ERR_BAD_PARAMETER_VALUE, false, "Invalid value for -convolve"); return false; } + } + + if (!subScaleKernel(data, widths, &size, &footprint)) { + psError(PS_ERR_UNKNOWN, false, "Unable to scale kernel parameters"); + return false; } @@ -177,4 +324,11 @@ if (threads > 0) { pmSubtractionThreadsInit(inRO, refRO); + } + + if (inRO->covariance) { + psImageCovarianceTruncate(inRO->covariance, COVAR_FRAC); + } + if (refRO->covariance) { + psImageCovarianceTruncate(refRO->covariance, COVAR_FRAC); } @@ -183,12 +337,33 @@ if (kernelRO) { success = pmSubtractionMatchPrecalc(inConv, refConv, inRO, refRO, kernelRO->analysis, - stride, kernelErr, maskVal, maskBad, maskPoor, poorFrac, badFrac); + stride, kernelErr, covarFrac, maskVal, maskBad, maskPoor, + poorFrac, badFrac); } else { success = pmSubtractionMatch(inConv, refConv, inRO, refRO, footprint, stride, regionSize, - spacing, threshold, sources, data->stamps, type, size, order, + spacing, threshold, detections->allSources, data->stamps, type, size, order, widths, orders, inner, ringsOrder, binning, penalty, optimum, - optWidths, optOrder, optThresh, iter, rej, sysErr, kernelErr, maskVal, - maskBad, maskPoor, poorFrac, badFrac, subMode); - } + optWidths, optOrder, optThresh, iter, rej, normFrac, + sysErr, skyErr, kernelErr, covarFrac, maskVal, maskBad, maskPoor, + poorFrac, badFrac, subMode); + } + +# ifdef TESTING + // XXX for testing + psphotSaveImage (NULL, refRO->image, "refRO.im.fits"); + psphotSaveImage (NULL, refRO->variance, "refRO.wt.fits"); + psphotSaveImage (NULL, refRO->mask, "refRO.mk.fits"); + + psphotSaveImage (NULL, inRO->image, "inRO.im.fits"); + psphotSaveImage (NULL, inRO->variance, "inRO.wt.fits"); + psphotSaveImage (NULL, inRO->mask, "inRO.mk.fits"); + + psphotSaveImage (NULL, inConv->image, "inConv.im.fits"); + psphotSaveImage (NULL, inConv->variance, "inConv.wt.fits"); + psphotSaveImage (NULL, inConv->mask, "inConv.mk.fits"); + + psphotSaveImage (NULL, refConv->image, "refConv.im.fits"); + psphotSaveImage (NULL, refConv->variance, "refConv.wt.fits"); + psphotSaveImage (NULL, refConv->mask, "refConv.mk.fits"); +# endif psFree(optWidths); @@ -261,4 +436,11 @@ pmConceptsCopyFPA(refConv->parent->parent->parent, refRO->parent->parent->parent, true, true); + if (inConv->covariance) { + psImageCovarianceTruncate(inConv->covariance, COVAR_FRAC); + } + if (refConv->covariance) { + psImageCovarianceTruncate(refConv->covariance, COVAR_FRAC); + } + if (inConv->variance) { psImageCovarianceTransfer(inConv->variance, inConv->covariance); Index: trunk/ppSub/src/ppSubReadoutJpeg.c =================================================================== --- trunk/ppSub/src/ppSubReadoutJpeg.c (revision 26571) +++ trunk/ppSub/src/ppSubReadoutJpeg.c (revision 26899) @@ -49,2 +49,94 @@ return true; } + +bool ppSubResidualSampleJpeg(pmConfig *config) +{ + return true; + + pmFPAview *view = ppSubViewReadout(); // View to readout + + // we save sample difference stamps on the convolved image readout->analysis metadata + pmReadout *inConv = pmFPAfileThisReadout(config->files, view, "PPSUB.INPUT.CONV"); // Input convolved + + psArray *samples = psArrayAllocEmpty(16); // Array of sample stamp images + + // we may have multiple entries with the same name; pull them off into a separate array: + psString regex = NULL; // Regular expression + psStringAppend(®ex, "^%s$", "SUBTRACTION.SAMPLE.STAMP.SET"); + psMetadataIterator *iter = psMetadataIteratorAlloc(inConv->analysis, PS_LIST_HEAD, regex); // Iterator + psFree(regex); + + psMetadataItem *item = NULL;// Item from iteration + while ((item = psMetadataGetAndIncrement(iter))) { + assert(item->type == PS_DATA_ARRAY); + psArray *sampleStamps = item->data.V; + samples = psArrayAdd(samples, 16, sampleStamps); + } + psFree(iter); + psAssert (samples, "no sample stamps?"); + psAssert (samples->n, "no sample stamps?"); + + // get the kernel sizes + psArray *kernels = samples->data[0]; + psAssert (kernels, "no valid kernel?"); + psAssert (kernels->n, "no valid kernel?"); + + psImage *kernel = kernels->data[0]; + psAssert (kernel, "missing valid kernel?"); + + int DX = kernel->numCols; + int DY = kernel->numRows; + + // each array contains up to 9 sample stamps. generate an image mosaicking these together in 3x3 blocks. + int innerBorder = 1; + int outerBorder = 2; + int NXblock = sqrt(samples->n); + int NYblock = samples->n / NXblock; + if (samples->n % NXblock) NYblock ++; + + int NXpix = NXblock * (3 * (DX + innerBorder) + outerBorder); + int NYpix = NYblock * (3 * (DY + innerBorder) + outerBorder); + + // output cell + pmCell *cell = pmFPAfileThisCell(config->files, view, "PPSUB.OUTPUT.RESID.JPEG"); + pmReadout *readout = pmReadoutAlloc(cell); + readout->image = psImageAlloc(NXpix, NYpix, PS_TYPE_F32); + + cell->data_exists = true; + cell->parent->data_exists = true; + + psImageInit (readout->image, 0.0); + + for (int i = 0; i < samples->n; i++) { + + int xBlock = i % NXblock; + int yBlock = i / NYblock; + + psArray *kernels = samples->data[i]; + + for (int j = 0; j < kernels->n; j++) { + + psImage *kernel = kernels->data[j]; + + int xSub = j % 3; + int ySub = j / 3; + + int xPix = xBlock * (3 * (DX + innerBorder) + outerBorder) + xSub * (DX + innerBorder); + int yPix = yBlock * (3 * (DX + innerBorder) + outerBorder) + ySub * (DY + innerBorder); + + for (int y = 0; y < kernel->numRows; y++) { + for (int x = 0; x < kernel->numCols; x++) { + readout->image->data.F32[y + yPix][x + xPix] = kernel->data.F32[y][x]; + } + } + } + } + + { + psFits *fits = psFitsOpen ("resid.stamps.fits", "w"); + psFitsWriteImage (fits, NULL, readout->image, 0, NULL); + psFitsClose (fits); + } + + return true; +} Index: trunk/ppSub/src/ppSubReadoutPhotometry.c =================================================================== --- trunk/ppSub/src/ppSubReadoutPhotometry.c (revision 26571) +++ trunk/ppSub/src/ppSubReadoutPhotometry.c (revision 26899) @@ -24,4 +24,6 @@ bool ppSubReadoutPhotometry(const char *name, ppSubData *data) { + bool mdok = false; + psAssert(data, "Require processing data"); pmConfig *config = data->config; // Configuration @@ -42,31 +44,10 @@ } - // 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 + // select the view of interest pmFPAview *view = ppSubViewReadout(); // View to readout - pmChip *psfInputChip = pmFPAfileThisChip(config->files, view, "PSPHOT.INPUT"); // Chip with PSF - psAssert (psfInputChip, "should have been generated for ppSubMakePSF"); - pmChip *psfLoadChip = pmFPAfileThisChip(config->files, view, "PSPHOT.PSF.LOAD"); // Chip to have PSF - psAssert (psfLoadChip, "PSPHOT.PSF.LOAD should have been defined in ppSubCamera"); - pmPSF *psf = psMetadataLookupPtr(NULL, psfInputChip->analysis, "PSPHOT.PSF"); // PSF for photometry - if (!psf) { - psErrorStackPrint(stderr, "No PSF available"); - psWarning("No PSF available --- suspect bad data quality."); - ppSubDataQuality(data, psErrorCodeLast(), PPSUB_FILES_PHOT_SUB | PPSUB_FILES_PHOT_INV); - return true; - } - psMetadataAddPtr(psfLoadChip->analysis, PS_LIST_TAIL, "PSPHOT.PSF", PS_DATA_UNKNOWN | PS_META_REPLACE, - "PSF from ppSubMakePSF", psf); - - bool mdok = false; // 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. - pmReadout *inRO = pmFPAfileThisReadout(config->files, view, name); // Readout with image and sources - if (psMetadataLookup(inRO->analysis, "PSPHOT.SOURCES")) { - psMetadataRemoveKey(inRO->analysis, "PSPHOT.SOURCES"); - } - + // around the pointers so PSPHOT.INPUT corresponds to the output image of interest (on one + // pass this is the subtraction image, in another it is negative of the subtraction pmFPAfile *photFile = psMetadataLookupPtr(&mdok, config->files, "PSPHOT.INPUT"); // Photometry file pmFPAfile *inFile = psMetadataLookupPtr(&mdok, config->files, name); // Input file @@ -76,18 +57,27 @@ return false; } + + // drop references to PSPHOT.DETECTIONS on both of these (why is this needed for both??) pmReadout *photRO = pmFPAviewThisReadout(view, photFile->fpa); // Readout to photometer - if (psMetadataLookup(photRO->analysis, "PSPHOT.SOURCES")) { - psMetadataRemoveKey(photRO->analysis, "PSPHOT.SOURCES"); + if (psMetadataLookup(photRO->analysis, "PSPHOT.DETECTIONS")) { + psMetadataRemoveKey(photRO->analysis, "PSPHOT.DETECTIONS"); + } + pmReadout *inRO = pmFPAfileThisReadout(config->files, view, name); // Readout with image and sources + if (psMetadataLookup(inRO->analysis, "PSPHOT.DETECTIONS")) { + psMetadataRemoveKey(inRO->analysis, "PSPHOT.DETECTIONS"); } - psMetadataAddPtr(photRO->parent->parent->analysis, PS_LIST_TAIL, "PSPHOT.PSF", - PS_META_REPLACE | PS_DATA_UNKNOWN, "Point-spread function", data->psf); - - psFree(photRO->analysis); - photRO->analysis = psMetadataAlloc(); + // grab the PSF information from the pmFPAfile PSPHOT.PSF.LOAD + pmFPAfile *psfFile = psMetadataLookupPtr(&mdok, config->files, "PSPHOT.PSF.LOAD"); // PSF file + ppSubCopyPSF (photFile, psfFile, view); // psphotSaveImage (photFile->fpa->hdu->header, photRO->image, "findsrc.im.fits"); // psphotSaveImage (photFile->fpa->hdu->header, photRO->variance, "findsrc.wt.fits"); // psphotSaveImage (photFile->fpa->hdu->header, photRO->mask, "findsrc.mk.fits"); + + // erase the overlays from a previous psphot-related step + if (pmVisualIsVisual()) { + // psphotVisualEraseOverlays (1, "all"); + } if (!psphotReadoutMinimal(config, view)) { @@ -100,8 +90,10 @@ // If no sources were found, there's no error, but we want to trigger 'bad quality' - psArray *sources = psMetadataLookupPtr(NULL, photRO->analysis, "PSPHOT.SOURCES"); // Sources - if (!sources) { + pmDetections *detections = psMetadataLookupPtr(NULL, photRO->analysis, "PSPHOT.DETECTIONS"); // Sources + if (!detections) { ppSubDataQuality(data, PSPHOT_ERR_DATA, PPSUB_FILES_PHOT_SUB | PPSUB_FILES_PHOT_INV); } + psArray *sources = detections->allSources; + psAssert (sources, "missing sources?"); if (data->stats) { @@ -118,6 +110,6 @@ if (!data->quality) { - if (!psMetadataCopySingle(inRO->analysis, photRO->analysis, "PSPHOT.SOURCES")) { - psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to copy PSPHOT.SOURCES"); + if (!psMetadataCopySingle(inRO->analysis, photRO->analysis, "PSPHOT.DETECTIONS")) { + psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to copy PSPHOT.DETECTIONS"); return false; } @@ -125,4 +117,15 @@ psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to copy PSPHOT.HEADER"); return false; + } + if (!psMetadataCopySingle(inRO->analysis, photRO->analysis, PM_DETEFF_ANALYSIS)) { + psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to copy Detection Efficiency"); + return false; + } + + // Ensure photometry information is put in the header + pmHDU *hdu = pmHDUFromReadout(inRO); // HDU for readout + if (hdu) { + psMetadata *photHeader = psMetadataLookupMetadata(NULL, inRO->analysis, "PSPHOT.HEADER"); // Header + hdu->header = psMetadataCopy(hdu->header, photHeader); } } @@ -134,13 +137,10 @@ } - - - - #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 + pmDetections *detections = psMetadataLookupPtr(NULL, photRO->analysis, "PSPHOT.DETECTIONS"); // Sources + psArray *sources = detections->allSources; FILE *sourceFile = fopen("sources.dat", "w"); // File for sources fprintf(sourceFile, Index: trunk/ppSub/src/ppSubThreshold.c =================================================================== --- trunk/ppSub/src/ppSubThreshold.c (revision 26571) +++ trunk/ppSub/src/ppSubThreshold.c (revision 26899) @@ -27,4 +27,5 @@ psImageMaskType maskIgnore, // Ignore pixels with this mask psImageMaskType maskThresh, // Give pixels this mask if below threshold + psRegion *region, // Region of interest const char *description // Description of image ) @@ -37,9 +38,11 @@ psImage *image = ro->image; // Image psImage *mask = ro->mask; // Mask - int numCols = ro->image->numCols, numRows = ro->image->numRows; // Size of image + + psImage *subImage = psImageSubset(image, *region); // Image with region of interest + psImage *subMask = psImageSubset(mask, *region); // Maks with region of interest psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV); // Statistics psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS); // Random number generator - if (!psImageBackground(stats, NULL, image, mask, maskIgnore, rng)) { + if (!psImageBackground(stats, NULL, subImage, subMask, maskIgnore, rng)) { psError(PS_ERR_UNKNOWN, false, "Unable to determine threshold."); psFree(rng); @@ -49,10 +52,13 @@ psFree(rng); + psFree(subImage); + psFree(subMask); + float threshold = stats->robustMedian - thresh * stats->robustStdev; // Threshold below which to clip psFree(stats); psLogMsg("ppSub", PS_LOG_INFO, "Masking pixels below %f in %s", threshold, description); - for (int y = 0; y < numRows; y++) { - for (int x = 0; x < numCols; x++) { + for (int y = region->y0; y < region->y1; y++) { + for (int x = region->x0; x < region->x1; x++) { if (mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & maskIgnore) { continue; @@ -94,8 +100,4 @@ return false; } - if (!lowThreshold(in, thresh, maskVal, maskThresh, "input convolved image")) { - psError(PS_ERR_UNKNOWN, false, "Unable to threshold input image."); - return false; - } pmReadout *ref = pmFPAfileThisReadout(config->files, view, "PPSUB.REF.CONV"); // Reference image @@ -104,8 +106,21 @@ return false; } - if (!lowThreshold(ref, thresh, maskVal, maskThresh, "reference convolved image")) { - psError(PS_ERR_UNKNOWN, false, "Unable to threshold input image."); - return false; + + psMetadataIterator *regIter = psMetadataIteratorAlloc(in->analysis, PS_LIST_HEAD, + "^" PM_SUBTRACTION_ANALYSIS_REGION "$"); + psMetadataItem *regItem; // Item with region + while ((regItem = psMetadataGetAndIncrement(regIter))) { + psAssert(regItem->type == PS_DATA_REGION && regItem->data.V, "Expect region type"); + psRegion *region = regItem->data.V; // Region of interest + if (!lowThreshold(in, thresh, maskVal, maskThresh, region, "input convolved image")) { + psError(PS_ERR_UNKNOWN, false, "Unable to threshold input image."); + return false; + } + if (!lowThreshold(ref, thresh, maskVal, maskThresh, region, "reference convolved image")) { + psError(PS_ERR_UNKNOWN, false, "Unable to threshold input image."); + return false; + } } + psFree(regIter); psFree(view);