Index: /branches/eam_branches/gDiff.20200202/src/Makefile.am =================================================================== --- /branches/eam_branches/gDiff.20200202/src/Makefile.am (revision 41412) +++ /branches/eam_branches/gDiff.20200202/src/Makefile.am (revision 41413) @@ -8,5 +8,5 @@ gDiffVersion.c: gDiffVersionDefinitions.h gDiffVersionDefinitions.h: gDiffVersionDefinitions.h.in FORCE - ../pslib-setsvnversion.pl GDIFF gDiffVersionDefinitions.h.in gDiffVersionDefinitions.h + pslib-setsvnversion.pl GDIFF gDiffVersionDefinitions.h.in gDiffVersionDefinitions.h FORCE: ; @@ -17,29 +17,28 @@ gDiff.c \ gDiffArguments.c \ - gDiffVersion.c \ + gDiffBackground.c \ gDiffCamera.c \ gDiffData.c \ + gDiffDefineOutput.c \ gDiffErrorCodes.c \ + gDiffExit.c \ + gDiffExtras.c \ + gDiffFiles.c \ + gDiffFlagNeighbors.c \ + gDiffInputDetections.c \ gDiffLoop.c \ - gDiffExit.c \ - gDiffFiles.c \ - gDiffDefineOutput.c - -# gDiffBackground.c \ -# gDiffVarianceRescale.c \ -# gDiffInputDetections.c \ -# gDiffExtras.c \ -# gDiffFlagNeighbors.c \ -# gDiffMakePSF.c \ -# gDiffMatchPSFs.c \ -# gDiffMaskStats.c \ -# gDiffReadoutInverse.c \ -# gDiffReadoutJpeg.c \ -# gDiffReadoutPhotometry.c \ -# gDiffReadoutForcedPhot.c \ -# gDiffReadoutStats.c \ -# gDiffReadoutSubtract.c \ -# gDiffSetMasks.c \ -# gDiffThreshold.c + gDiffMakePSF.c \ + gDiffMaskStats.c \ + gDiffMatchPSFs.c \ + gDiffReadoutForcedPhot.c \ + gDiffReadoutInverse.c \ + gDiffReadoutJpeg.c \ + gDiffReadoutPhotometry.c \ + gDiffReadoutStats.c \ + gDiffReadoutSubtract.c \ + gDiffSetMasks.c \ + gDiffThreshold.c \ + gDiffVarianceRescale.c \ + gDiffVersion.c noinst_HEADERS = \ Index: /branches/eam_branches/gDiff.20200202/src/gDiff.c =================================================================== --- /branches/eam_branches/gDiff.20200202/src/gDiff.c (revision 41412) +++ /branches/eam_branches/gDiff.20200202/src/gDiff.c (revision 41413) @@ -18,4 +18,6 @@ psphotErrorRegister(); + psExit exitValue = PS_EXIT_SUCCESS; // default value + gDiffData *data = NULL; // Processing data pmConfig *config = pmConfigRead(&argc, argv, GDIFF_RECIPE); // Configuration @@ -26,6 +28,8 @@ psphotInit(); // cannot return false + // the 'data' container carries pointers to all data structures used by gDiff data = gDiffDataAlloc(config); + // parse the command-line options if (!gDiffArguments(argc, argv, data)) { psError(psErrorCodeLast(), false, "Error reading arguments."); @@ -33,5 +37,5 @@ } - // determine the camera and setup the I/O file handles + // determine the camera and setup the (many) I/O file handles if (!gDiffCamera(data)) goto die; @@ -40,54 +44,52 @@ die: - { - psExit exitValue = gDiffExitCode(PS_EXIT_SUCCESS); // Exit code + exitValue = gDiffExitCode(exitValue); // Exit code - if (data && data->stats && data->statsFile) { - psMetadataAddF32(data->stats, PS_LIST_TAIL, "TIME_DIFF", 0, "Total time (sec)", psTimerMark("gDiff")); - psString stats = psMetadataConfigFormat(data->stats); // Statistics to output - if (!stats || strlen(stats) == 0) { - psError(GDIFF_ERR_IO, false, "Unable to format statistics file"); - } else if (fprintf(data->statsFile, "%s", stats) != strlen(stats)) { - psError(GDIFF_ERR_IO, true, "Unable to write statistics file"); - } - psFree(stats); - if (fclose(data->statsFile) == EOF) { - psError(GDIFF_ERR_IO, true, "Unable to close statistics file"); - } - data->statsFile = NULL; - pmConfigRunFilenameAddWrite(data->config, "STATS", data->statsName); - exitValue = gDiffExitCode(exitValue); - } + if (data && data->stats && data->statsFile) { + psMetadataAddF32(data->stats, PS_LIST_TAIL, "TIME_DIFF", 0, "Total time (sec)", psTimerMark("gDiff")); + psString stats = psMetadataConfigFormat(data->stats); // Statistics to output + if (!stats || strlen(stats) == 0) { + psError(GDIFF_ERR_IO, false, "Unable to format statistics file"); + } else if (fprintf(data->statsFile, "%s", stats) != strlen(stats)) { + psError(GDIFF_ERR_IO, true, "Unable to write statistics file"); + } + psFree(stats); + if (fclose(data->statsFile) == EOF) { + psError(GDIFF_ERR_IO, true, "Unable to close statistics file"); + } + data->statsFile = NULL; + pmConfigRunFilenameAddWrite(data->config, "STATS", data->statsName); + exitValue = gDiffExitCode(exitValue); + } - if (config && !gDiffFilesIterateUp(config, GDIFF_FILES_ALL)) { - psError(psErrorCodeLast(), false, "Unable to close files."); - exitValue = gDiffExitCode(exitValue); - pmFPAfileFreeSetStrict(false); - } + if (config && !gDiffFilesIterateUp(config, GDIFF_FILES_ALL)) { + psError(psErrorCodeLast(), false, "Unable to close files."); + exitValue = gDiffExitCode(exitValue); + pmFPAfileFreeSetStrict(false); + } - if (data) { - psString dump_file = psMetadataLookupStr(NULL, data->config->arguments, "DUMP_CONFIG"); - if (dump_file) { - if (!pmConfigDump(data->config, dump_file)) { - psError(GDIFF_ERR_IO, false, "Unable to dump configuration."); - exitValue = gDiffExitCode(exitValue); - } - } - psFree(data); - } + if (data) { + psString dump_file = psMetadataLookupStr(NULL, data->config->arguments, "DUMP_CONFIG"); + if (dump_file) { + if (!pmConfigDump(data->config, dump_file)) { + psError(GDIFF_ERR_IO, false, "Unable to dump configuration."); + exitValue = gDiffExitCode(exitValue); + } + } + psFree(data); + } - psTrace("gDiff", 1, "Finished at %f sec\n", psTimerMark("gDiff")); - psLogMsg("gDiff", PS_LOG_INFO, "Complete gDiff run: %f sec\n", psTimerMark("gDiff")); - psTimerStop(); + psTrace("gDiff", 1, "Finished at %f sec\n", psTimerMark("gDiff")); + psLogMsg("gDiff", PS_LOG_INFO, "Complete gDiff run: %f sec\n", psTimerMark("gDiff")); + psTimerStop(); - pmVisualClose(); // close plot windows, if -visual is set - pmModelClassCleanup(); - pmConfigDone(); - pmVisualCleanup (); - psLibFinalize(); + pmVisualClose(); // close plot windows, if -visual is set + pmModelClassCleanup(); + pmConfigDone(); + pmVisualCleanup (); + psLibFinalize(); - fprintf (stderr, "found %d leaks at %s\n", psMemCheckLeaks (0, NULL, stdout, false), "gDiff"); - exitValue = gDiffExitCode(exitValue); - exit(exitValue); - } + fprintf (stderr, "found %d leaks at %s\n", psMemCheckLeaks (0, NULL, stdout, false), "gDiff"); + exitValue = gDiffExitCode(exitValue); + exit(exitValue); } Index: /branches/eam_branches/gDiff.20200202/src/gDiff.h =================================================================== --- /branches/eam_branches/gDiff.20200202/src/gDiff.h (revision 41412) +++ /branches/eam_branches/gDiff.20200202/src/gDiff.h (revision 41413) @@ -21,4 +21,5 @@ #include #include +#include #include Index: /branches/eam_branches/gDiff.20200202/src/gDiffArguments.c =================================================================== --- /branches/eam_branches/gDiff.20200202/src/gDiffArguments.c (revision 41412) +++ /branches/eam_branches/gDiff.20200202/src/gDiffArguments.c (revision 41413) @@ -1,2 +1,10 @@ +/** @file gDiffArguments.c + * @brief: image subtraction based on matched Gaussians profiles + * @author Eugene Magnier @ IfA + * @version $Revision: 1.12 $ + * @date $Date: 2020-01-25 $ + * Copyright 2020 Institute for Astronomy, University of Hawaii + */ + #include "gDiff.h" @@ -60,21 +68,21 @@ psMetadata *arguments = config->arguments; // Command-line arguments - psMetadataAddStr(arguments, PS_LIST_TAIL, "-inimage", 0, "Input image", NULL); - psMetadataAddStr(arguments, PS_LIST_TAIL, "-inmask", 0, "Input mask image", NULL); - psMetadataAddStr(arguments, PS_LIST_TAIL, "-invariance", 0, "Input variance image", NULL); - psMetadataAddStr(arguments, PS_LIST_TAIL, "-insources", 0, "Input source list", NULL); - psMetadataAddStr(arguments, PS_LIST_TAIL, "-refimage", 0, "Reference image", NULL); - psMetadataAddStr(arguments, PS_LIST_TAIL, "-refmask", 0, "Reference mask image", NULL); - psMetadataAddStr(arguments, PS_LIST_TAIL, "-refvariance", 0, "Reference variance image", NULL); - psMetadataAddStr(arguments, PS_LIST_TAIL, "-refsources", 0, "Reference source list", NULL); - psMetadataAddStr(arguments, PS_LIST_TAIL, "-kernel", 0, "Precalculated kernel to apply", NULL); - psMetadataAddStr(arguments, PS_LIST_TAIL, "-stats", 0, "Statistics file", NULL); - psMetadataAddStr(arguments, PS_LIST_TAIL, "-stamps", 0, "Stamps filename; x,y on each line", NULL); - psMetadataAddS32(arguments, PS_LIST_TAIL, "-convolve", 0, "Image to convolve [1 or 2]", 0); + psMetadataAddStr (arguments, PS_LIST_TAIL, "-inimage", 0, "Input image", NULL); + psMetadataAddStr (arguments, PS_LIST_TAIL, "-inmask", 0, "Input mask image", NULL); + psMetadataAddStr (arguments, PS_LIST_TAIL, "-invariance", 0, "Input variance image", NULL); + psMetadataAddStr (arguments, PS_LIST_TAIL, "-insources", 0, "Input source list", NULL); + psMetadataAddStr (arguments, PS_LIST_TAIL, "-refimage", 0, "Reference image", NULL); + psMetadataAddStr (arguments, PS_LIST_TAIL, "-refmask", 0, "Reference mask image", NULL); + psMetadataAddStr (arguments, PS_LIST_TAIL, "-refvariance", 0, "Reference variance image", NULL); + psMetadataAddStr (arguments, PS_LIST_TAIL, "-refsources", 0, "Reference source list", NULL); + psMetadataAddStr (arguments, PS_LIST_TAIL, "-kernel", 0, "Precalculated kernel to apply", NULL); + psMetadataAddStr (arguments, PS_LIST_TAIL, "-stats", 0, "Statistics file", NULL); + psMetadataAddStr (arguments, PS_LIST_TAIL, "-stamps", 0, "Stamps filename; x,y on each line", NULL); + psMetadataAddS32 (arguments, PS_LIST_TAIL, "-convolve", 0, "Image to convolve [1 or 2]", 0); psMetadataAddBool(arguments, PS_LIST_TAIL, "-photometry", 0, "Perform photometry?", NULL); psMetadataAddBool(arguments, PS_LIST_TAIL, "-forced-phot", 0, "Perform forced photometry?", NULL); psMetadataAddBool(arguments, PS_LIST_TAIL, "-forced-input1", 0, "Perform forced photometry?", NULL); psMetadataAddBool(arguments, PS_LIST_TAIL, "-forced-input2", 0, "Perform forced photometry?", NULL); - psMetadataAddF32(arguments, PS_LIST_TAIL, "-zp", 0, "Zero point for photometry", NAN); + psMetadataAddF32 (arguments, PS_LIST_TAIL, "-zp", 0, "Zero point for photometry", NAN); psMetadataAddBool(arguments, PS_LIST_TAIL, "-inverse", 0, "Generate inverse subtractions?", false); psMetadataAddBool(arguments, PS_LIST_TAIL, "-save-inconv", 0, "Save input convolved images?", false); Index: /branches/eam_branches/gDiff.20200202/src/gDiffBackground.c =================================================================== --- /branches/eam_branches/gDiff.20200202/src/gDiffBackground.c (revision 41413) +++ /branches/eam_branches/gDiff.20200202/src/gDiffBackground.c (revision 41413) @@ -0,0 +1,76 @@ +/** @file gDiffBackground.c + * @brief: image subtraction based on matched Gaussians profiles + * @author Eugene Magnier @ IfA + * @version $Revision: 1.12 $ + * @date $Date: 2020-01-25 $ + * Copyright 2020 Institute for Astronomy, University of Hawaii + */ + +#include "gDiff.h" + +bool gDiffBackground(pmConfig *config) +{ + psAssert(config, "Require configuration"); + + bool mdok; // Status of metadata lookups + + psMetadata *gDiffRecipe = psMetadataLookupPtr(NULL, config->recipes, GDIFF_RECIPE); // Recipe for gDiff + psAssert(gDiffRecipe, "Need GDIFF recipe"); + + psMetadata *psphotRecipe = psMetadataLookupPtr(NULL, config->recipes, PSPHOT_RECIPE); // Recipe for psphot + psAssert(psphotRecipe, "Need PSPHOT recipe for binning"); + + bool doApplyMaskNaN = psMetadataLookupBool(&mdok, gDiffRecipe, "APPLY.PIXELNAN"); // NaN the pixels underneath masks + + psImageMaskType maskBad = pmConfigMaskGet("BLANK", config); // Bits to mask + + pmFPAview *view = gDiffViewReadout(); // View to readout + pmReadout *outRO = pmFPAfileThisReadout(config->files, view, "GDIFF.OUTPUT"); // Output image + + // Generate the background model + if (!psphotModelBackground(config, view, "GDIFF.OUTPUT")) { + psError(psErrorCodeLast(), false, "Unable to model background"); + psFree(view); + return false; + } + + // select the model readout (should now exist) + pmReadout *modelRO = pmFPAfileThisReadout(config->files, view, "PSPHOT.BACKMDL"); + if (!modelRO) { + psError(psErrorCodeLast(), false, "Unable to find background model"); + psFree(view); + return false; + } + psFree(view); + + // get background model parameters + psImageBinning *binning = psMetadataLookupPtr(&mdok, modelRO->analysis, "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 + int numCols = image->numCols, numRows = image->numRows; // Size of image + for (int y = 0; y < numRows; y++) { + for (int x = 0; x < numCols; x++) { + // special case 1: NAN the masked pixels + if (doApplyMaskNaN && mask && (mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & maskBad)) { + image->data.F32[y][x] = NAN; + continue; + } + + float value = psImageUnbinPixel(x + 0.5, y + 0.5, modelImage, binning); // Background value + + // special case 2: NAN & mask pixels without a valid background model + if (!isfinite(value)) { + image->data.F32[y][x] = NAN; + mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] |= maskBad; + continue; + } + + // pixel & background are good: subtract + image->data.F32[y][x] -= value; + } + } + return true; +} Index: /branches/eam_branches/gDiff.20200202/src/gDiffCamera.c =================================================================== --- /branches/eam_branches/gDiff.20200202/src/gDiffCamera.c (revision 41412) +++ /branches/eam_branches/gDiff.20200202/src/gDiffCamera.c (revision 41413) @@ -1,2 +1,10 @@ +/** @file gDiffCamera.c + * @brief: image subtraction based on matched Gaussians profiles + * @author Eugene Magnier @ IfA + * @version $Revision: 1.12 $ + * @date $Date: 2020-01-25 $ + * Copyright 2020 Institute for Astronomy, University of Hawaii + */ + #include "gDiff.h" @@ -207,4 +215,19 @@ } + /* in ppSub, there is a recipe value 'PHOTOMETRY' which is supposed to enable / disable photometry analysis. + This is set to FALSE in the recipes by default. + + In the script diff_skycell.pl, there is a global variable 'do_photom' which is TRUE by default, + but is set to FALSE if a) the reduction class is NOCONVDIFF or b) the run_state is not new (i.e, 'update' mode) + (the script changes the call in an inconsistent way, if do_photom is true, it appends -photometry, but + if it is false, it appends -Db PHOTOMETRY FALSE. this could be made consistent. + (also note that it is bad practice to have hard-wired reduction classes) + + But both of these situations are ignored because the test in ppSubCamera for data->photometry + (set if PHOTOMETRY is true or -photometry is supplied) is ignored and the photometric analysis + is always performed. + + */ + // Now that the camera has been determined, we can read the recipe psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, GDIFF_RECIPE); // Recipe for ppSim @@ -380,94 +403,108 @@ } - // psPhot input - if (data->photometry || 1) { - psphotModelClassInit(); // load implementation-specific models - - pmFPAfile *psphot = pmFPAfileDefineFromFPA(config, output->fpa, 1, 1, "PSPHOT.INPUT"); - if (!psphot) { - psError(psErrorCodeLast(), false, "Failed to build FPA from PSPHOT.INPUT"); - return false; - } - if (psphot->type != PM_FPA_FILE_IMAGE) { - psError(psErrorCodeLast(), true, "PSPHOT.INPUT is not of type IMAGE"); - return false; - } - pmFPAfileActivate(config->files, false, "PSPHOT.INPUT"); - - // Internal file for getting the PSF from the minuend - pmFPAfile *psf = pmFPAfileDefineFromFPA(config, output->fpa, 1, 1, "PSPHOT.PSF.LOAD"); - if (!psf) { - psError(psErrorCodeLast(), false, "Failed to build FPA from PSPHOT.PSF.LOAD"); - return false; - } - if (psf->type != PM_FPA_FILE_PSF) { - psError(psErrorCodeLast(), true, "PSPHOT.PSF.LOAD is not of type PSF"); - return false; - } - pmFPAfileActivate(config->files, false, "PSPHOT.PSF.LOAD"); - - if (!psphotDefineFiles(config, psphot)) { - psError(psErrorCodeLast(), false, "Unable to set up psphot files."); - return false; - } - - // Deactivate psphot output sources --- we want to define output source files of our own - pmFPAfile *psphotOutput = pmFPAfileSelectSingle(config->files, "PSPHOT.OUTPUT", 0); - psphotOutput->save = false; // this one should NOT be set - - pmFPAfile *outSources = defineOutputFile(config, output, false, "GDIFF.OUTPUT.SOURCES", - PM_FPA_FILE_CMF); - if (!outSources) { - psError(psErrorCodeLast(), false, "Unable to set up output source file."); - return false; - } - checkFileruleFileSave(outSources, config); - - if (data->inverse) { - pmFPAfile *invSources = defineOutputFile(config, inverse, false, "GDIFF.INVERSE.SOURCES", - PM_FPA_FILE_CMF); - if (!invSources) { - psError(psErrorCodeLast(), false, "Unable to set up inverse source file."); - return false; - } - checkFileruleFileSave(invSources, config); - } - - // files need to do the forced photometry on the positions of sources in the positive images - if (data->forcedPhot1) { - // this pmFPAfile is used to carry sources detected in the positive image #1 - pmFPAfile *posSources1 = defineOutputFile(config, input, true, "GDIFF.POS1.SOURCES", PM_FPA_FILE_CMF); - if (!posSources1) { - psError(psErrorCodeLast(), false, "Unable to set up forced source file."); - return false; - } - checkFileruleFileSave(posSources1, config); - - // this pmFPAfile is used to carry sources detected in the diff image @ the positions from positive image #1 - pmFPAfile *frcSources1 = defineOutputFile(config, input, true, "GDIFF.FORCED1.SOURCES", PM_FPA_FILE_CMF); - if (!frcSources1) { - psError(psErrorCodeLast(), false, "Unable to set up forced source file."); - return false; - } - checkFileruleFileSave(frcSources1, config); + /* in ppSub, there is a recipe value 'PHOTOMETRY' which is supposed to enable / disable photometry analysis. + This is set to FALSE in the recipes by default. + + In the script diff_skycell.pl, there is a global variable 'do_photom' which is TRUE by default, + but is set to FALSE if a) the reduction class is NOCONVDIFF or b) the run_state is not new (i.e, 'update' mode) + (the script changes the call in an inconsistent way, if do_photom is true, it appends -photometry, but + if it is false, it appends -Db PHOTOMETRY FALSE. this could be made consistent. + (also note that it is bad practice to have hard-wired reduction classes) + + But both of these situations are ignored because the test in ppSubCamera for data->photometry + (set if PHOTOMETRY is true or -photometry is supplied) is ignored and the photometric analysis + is always performed. + + Should we rationalize this? + */ + + // define the psphot I/O files + psphotModelClassInit(); // load implementation-specific models + + pmFPAfile *psphot = pmFPAfileDefineFromFPA(config, output->fpa, 1, 1, "PSPHOT.INPUT"); + if (!psphot) { + psError(psErrorCodeLast(), false, "Failed to build FPA from PSPHOT.INPUT"); + return false; + } + if (psphot->type != PM_FPA_FILE_IMAGE) { + psError(psErrorCodeLast(), true, "PSPHOT.INPUT is not of type IMAGE"); + return false; + } + pmFPAfileActivate(config->files, false, "PSPHOT.INPUT"); + + // Internal file for getting the PSF from the minuend + pmFPAfile *psf = pmFPAfileDefineFromFPA(config, output->fpa, 1, 1, "PSPHOT.PSF.LOAD"); + if (!psf) { + psError(psErrorCodeLast(), false, "Failed to build FPA from PSPHOT.PSF.LOAD"); + return false; + } + if (psf->type != PM_FPA_FILE_PSF) { + psError(psErrorCodeLast(), true, "PSPHOT.PSF.LOAD is not of type PSF"); + return false; + } + pmFPAfileActivate(config->files, false, "PSPHOT.PSF.LOAD"); + + if (!psphotDefineFiles(config, psphot)) { + psError(psErrorCodeLast(), false, "Unable to set up psphot files."); + return false; + } + + // Deactivate psphot output sources --- we want to define output source files of our own + pmFPAfile *psphotOutput = pmFPAfileSelectSingle(config->files, "PSPHOT.OUTPUT", 0); + psphotOutput->save = false; // this one should NOT be set + + pmFPAfile *outSources = defineOutputFile(config, output, false, "GDIFF.OUTPUT.SOURCES", + PM_FPA_FILE_CMF); + if (!outSources) { + psError(psErrorCodeLast(), false, "Unable to set up output source file."); + return false; + } + checkFileruleFileSave(outSources, config); + + if (data->inverse) { + pmFPAfile *invSources = defineOutputFile(config, inverse, false, "GDIFF.INVERSE.SOURCES", + PM_FPA_FILE_CMF); + if (!invSources) { + psError(psErrorCodeLast(), false, "Unable to set up inverse source file."); + return false; } - - if (data->forcedPhot2) { - // this pmFPAfile is used to carry sources detected in the positive image #2 - pmFPAfile *posSources2 = defineOutputFile(config, ref, true, "GDIFF.POS2.SOURCES", PM_FPA_FILE_CMF); - if (!posSources2) { - psError(psErrorCodeLast(), false, "Unable to set up forced source file."); - return false; - } - checkFileruleFileSave(posSources2, config); - - // this pmFPAfile is used to carry sources detected in the diff image @ the positions from positive image #2 - pmFPAfile *frcSources2 = defineOutputFile(config, ref, true, "GDIFF.FORCED2.SOURCES", PM_FPA_FILE_CMF); - if (!frcSources2) { - psError(psErrorCodeLast(), false, "Unable to set up forced source file."); - return false; - } - checkFileruleFileSave(frcSources2, config); - } + checkFileruleFileSave(invSources, config); + } + + // files need to do the forced photometry on the positions of sources in the positive images + if (data->forcedPhot1) { + // this pmFPAfile is used to carry sources detected in the positive image #1 + pmFPAfile *posSources1 = defineOutputFile(config, input, true, "GDIFF.POS1.SOURCES", PM_FPA_FILE_CMF); + if (!posSources1) { + psError(psErrorCodeLast(), false, "Unable to set up forced source file."); + return false; + } + checkFileruleFileSave(posSources1, config); + + // this pmFPAfile is used to carry sources detected in the diff image @ the positions from positive image #1 + pmFPAfile *frcSources1 = defineOutputFile(config, input, true, "GDIFF.FORCED1.SOURCES", PM_FPA_FILE_CMF); + if (!frcSources1) { + psError(psErrorCodeLast(), false, "Unable to set up forced source file."); + return false; + } + checkFileruleFileSave(frcSources1, config); + } + + if (data->forcedPhot2) { + // this pmFPAfile is used to carry sources detected in the positive image #2 + pmFPAfile *posSources2 = defineOutputFile(config, ref, true, "GDIFF.POS2.SOURCES", PM_FPA_FILE_CMF); + if (!posSources2) { + psError(psErrorCodeLast(), false, "Unable to set up forced source file."); + return false; + } + checkFileruleFileSave(posSources2, config); + + // this pmFPAfile is used to carry sources detected in the diff image @ the positions from positive image #2 + pmFPAfile *frcSources2 = defineOutputFile(config, ref, true, "GDIFF.FORCED2.SOURCES", PM_FPA_FILE_CMF); + if (!frcSources2) { + psError(psErrorCodeLast(), false, "Unable to set up forced source file."); + return false; + } + checkFileruleFileSave(frcSources2, config); } Index: /branches/eam_branches/gDiff.20200202/src/gDiffData.c =================================================================== --- /branches/eam_branches/gDiff.20200202/src/gDiffData.c (revision 41412) +++ /branches/eam_branches/gDiff.20200202/src/gDiffData.c (revision 41413) @@ -1,2 +1,10 @@ +/** @file gDiffData.c + * @brief: image subtraction based on matched Gaussians profiles + * @author Eugene Magnier @ IfA + * @version $Revision: 1.12 $ + * @date $Date: 2020-01-25 $ + * Copyright 2020 Institute for Astronomy, University of Hawaii + */ + #include "gDiff.h" Index: /branches/eam_branches/gDiff.20200202/src/gDiffDefineOutput.c =================================================================== --- /branches/eam_branches/gDiff.20200202/src/gDiffDefineOutput.c (revision 41412) +++ /branches/eam_branches/gDiff.20200202/src/gDiffDefineOutput.c (revision 41413) @@ -1,22 +1,9 @@ /** @file gDiffDefineOutput.c - * - * @brief - * - * @ingroup gDiff - * - * @author IfA - * @version $Revision: 1.3 $ $Name: not supported by cvs2svn $ - * @date $Date: 2009-02-18 00:31:20 $ - * Copyright 2009 Institute for Astronomy, University of Hawaii + * @brief: image subtraction based on matched Gaussians profiles + * @author Eugene Magnier @ IfA + * @version $Revision: 1.12 $ + * @date $Date: 2020-01-25 $ + * Copyright 2020 Institute for Astronomy, University of Hawaii */ - -#ifdef HAVE_CONFIG_H -#include -#endif - -#include -#include -#include -#include #include "gDiff.h" Index: /branches/eam_branches/gDiff.20200202/src/gDiffErrorCodes.dat =================================================================== --- /branches/eam_branches/gDiff.20200202/src/gDiffErrorCodes.dat (revision 41412) +++ /branches/eam_branches/gDiff.20200202/src/gDiffErrorCodes.dat (revision 41413) @@ -3,5 +3,5 @@ # BASE = 14000 First value we use; lower values belong to psLib -UNKNOWN Unknown ppSub error code +UNKNOWN Unknown gDiff error code NOT_IMPLEMENTED Desired feature is not yet implemented ARGUMENTS Incorrect arguments Index: /branches/eam_branches/gDiff.20200202/src/gDiffExit.c =================================================================== --- /branches/eam_branches/gDiff.20200202/src/gDiffExit.c (revision 41412) +++ /branches/eam_branches/gDiff.20200202/src/gDiffExit.c (revision 41413) @@ -1,5 +1,9 @@ -#include -#include -#include +/** @file gDiffExit.c + * @brief: image subtraction based on matched Gaussians profiles + * @author Eugene Magnier @ IfA + * @version $Revision: 1.12 $ + * @date $Date: 2020-01-25 $ + * Copyright 2020 Institute for Astronomy, University of Hawaii + */ #include "gDiff.h" @@ -16,8 +20,10 @@ // not include the possible psphot values - // for now, to get around this, we just use an int for the switch + // thus, we cannot do this: + // psErrorCode errorCode = psErrorCodeLast(); // Error code - // psErrorCode errorCode = psErrorCodeLast(); // Error code + // for now, to get around this, we just use an int for the switch: int errorCode = psErrorCodeLast(); // Error code + if (errorCode != PS_ERR_NONE) { pmFPAfileFreeSetStrict(false); Index: /branches/eam_branches/gDiff.20200202/src/gDiffExtras.c =================================================================== --- /branches/eam_branches/gDiff.20200202/src/gDiffExtras.c (revision 41413) +++ /branches/eam_branches/gDiff.20200202/src/gDiffExtras.c (revision 41413) @@ -0,0 +1,88 @@ +/** @file gDiffExtras.c + * @brief: image subtraction based on matched Gaussians profiles + * @author Eugene Magnier @ IfA + * @version $Revision: 1.12 $ + * @date $Date: 2020-01-25 $ + * Copyright 2020 Institute for Astronomy, University of Hawaii + */ + +#include "gDiff.h" + +// This file contains minor functions used in gDiff. +// Some of these should be pushed into psModules or psLib. + +// copy psphot result (detections, header, deteff data) from one pmFPAfile to another +bool psphotCopyResults (bool *foundDetections, pmFPAfile *target, pmFPAfile *source, pmFPAview *view) { + + pmReadout *sourceRO = pmFPAviewThisReadout(view, source->fpa); + psAssert(sourceRO, "programming error: source readout not defined"); + + pmReadout *targetRO = pmFPAviewThisReadout(view, target->fpa); + if (!targetRO) { + pmCell *targetCell = pmFPAviewThisCell(view, target->fpa); + targetRO = pmReadoutAlloc(targetCell); + psAssert(targetRO, "programming error: could not make target readout"); + } + + // If no sources were found, there's no error, but we need to inform the calling function + pmDetections *detections = psMetadataLookupPtr(NULL, sourceRO->analysis, "PSPHOT.DETECTIONS"); // Sources + if (!detections) { + *foundDetections = false; + return true; + } else { + *foundDetections = true; + } + + psArray *sources = detections->allSources; + psAssert (sources, "missing sources?"); + + if (!psMetadataCopySingle(targetRO->analysis, sourceRO->analysis, "PSPHOT.DETECTIONS")) { + psError(GDIFF_ERR_PROG, false, "Unable to copy PSPHOT.DETECTIONS"); + return false; + } + if (!psMetadataCopySingle(targetRO->analysis, sourceRO->analysis, "PSPHOT.HEADER")) { + psError(GDIFF_ERR_PROG, false, "Unable to copy PSPHOT.HEADER"); + return false; + } + if (!psMetadataCopySingle(targetRO->analysis, sourceRO->analysis, PM_DETEFF_ANALYSIS)) { + psError(GDIFF_ERR_PROG, false, "Unable to copy Detection Efficiency"); + return false; + } + + // Ensure photometry information is put in the header + // XXX create one if it does not exist? + pmHDU *hdu = pmHDUFromReadout(targetRO); // HDU for readout + if (hdu) { + psMetadata *header = psMetadataLookupMetadata(NULL, targetRO->analysis, "PSPHOT.HEADER"); // Header + hdu->header = psMetadataCopy(hdu->header, header); + } + + targetRO->data_exists = true; + targetRO->parent->data_exists = true; + targetRO->parent->parent->data_exists = true; + + return true; +} + +/** + * 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: /branches/eam_branches/gDiff.20200202/src/gDiffFiles.c =================================================================== --- /branches/eam_branches/gDiff.20200202/src/gDiffFiles.c (revision 41412) +++ /branches/eam_branches/gDiff.20200202/src/gDiffFiles.c (revision 41413) @@ -1,10 +1,9 @@ -#ifdef HAVE_CONFIG_H -#include -#endif - -#include -#include -#include -#include +/** @file gDiffFiles.c + * @brief: image subtraction based on matched Gaussians profiles + * @author Eugene Magnier @ IfA + * @version $Revision: 1.12 $ + * @date $Date: 2020-01-25 $ + * Copyright 2020 Institute for Astronomy, University of Hawaii + */ #include "gDiff.h" @@ -12,11 +11,16 @@ // Input files -static const char *inputFiles[] = { "GDIFF.INPUT", "GDIFF.INPUT.MASK", "GDIFF.INPUT.VARIANCE", - "GDIFF.INPUT.SOURCES", "GDIFF.REF", "GDIFF.REF.MASK", - "GDIFF.REF.VARIANCE", "GDIFF.REF.SOURCES", NULL }; +static const char *inputFiles[] = { "GDIFF.INPUT", + "GDIFF.INPUT.MASK", + "GDIFF.INPUT.VARIANCE", + "GDIFF.INPUT.SOURCES", + "GDIFF.REF", + "GDIFF.REF.MASK", + "GDIFF.REF.VARIANCE", + "GDIFF.REF.SOURCES", NULL }; // Convolved files static const char *convFiles[] = { "GDIFF.INPUT.CONV", "GDIFF.INPUT.CONV.MASK", "GDIFF.INPUT.CONV.VARIANCE", - "GDIFF.REF.CONV", "GDIFF.REF.CONV.MASK", "GDIFF.REF.CONV.VARIANCE", + "GDIFF.REF.CONV", "GDIFF.REF.CONV.MASK", "GDIFF.REF.CONV.VARIANCE", NULL }; @@ -234,2 +238,3 @@ return true; } + Index: /branches/eam_branches/gDiff.20200202/src/gDiffFlagNeighbors.c =================================================================== --- /branches/eam_branches/gDiff.20200202/src/gDiffFlagNeighbors.c (revision 41413) +++ /branches/eam_branches/gDiff.20200202/src/gDiffFlagNeighbors.c (revision 41413) @@ -0,0 +1,283 @@ +/** @file gDiffFlagNeighbors.c + * @brief: image subtraction based on matched Gaussians profiles + * @author Eugene Magnier @ IfA + * @version $Revision: 1.12 $ + * @date $Date: 2020-01-25 $ + * Copyright 2020 Institute for Astronomy, University of Hawaii + */ + +#include "gDiff.h" + +# define MIN_ERROR 0.02 /* padding for self-match flux test */ +# define MAX_OFFSET 5.0 /* delta mag greater than this is not a self match */ + +// match the diff sources to the detections from the positive or negative input images +// matchRef : false if this is the standard diff image, true if this is the inverse (in-ref vs ref-in) +// this is needed to choose which value is positive and which is negative relative to the difference +bool gDiffFlagNeighbors(pmConfig *config, pmFPAview *view, psArray *sources, bool subInverse) { + + bool status; + + psAssert(config, "Require config"); + psAssert(view, "Require config"); + psAssert (sources, "missing sources?"); + + // Look up recipe values + psMetadata *recipe = psMetadataLookupMetadata(&status, config->recipes, GDIFF_RECIPE); // should this be psphot? + psAssert(recipe, "We checked this earlier, so it should be here."); + + float radius = psMetadataLookupF32(&status, recipe, "INPUT.MATCH.RADIUS"); + psAssert (status, "cannot find INPUT.MATCH.RADIUS"); + + float minSN = psMetadataLookupF32(&status, recipe, "INPUT.MATCH.MIN.SN"); + psAssert (status, "cannot find INPUT.MATCH.MIN.SN"); + + // Input sources + pmReadout *inSourceRO = pmFPAfileThisReadout(config->files, view, "GDIFF.INPUT.SOURCES"); + pmDetections *inDetections = inSourceRO ? psMetadataLookupPtr(NULL, inSourceRO->analysis, "PSPHOT.DETECTIONS") : NULL; // Input sources + if (!inDetections) { + psWarning("Unable to filter detections based on image A."); + } + + pmReadout *refSourceRO = pmFPAfileThisReadout(config->files, view, "GDIFF.REF.SOURCES"); + pmDetections *refDetections = refSourceRO ? psMetadataLookupPtr(NULL, refSourceRO->analysis, "PSPHOT.DETECTIONS") : NULL; // Ref sources + if (!refDetections) { + psWarning("Unable to filter detections based on image B."); + } + if (!inDetections && !refDetections) { + return true; + } + + psArray *inSources = inDetections->allSources; + psAssert (inSources, "missing inSources?"); + + psArray *refSources = refDetections->allSources; + psAssert (refSources, "missing refSources?"); + + gDiffSetSourceImageIDs(sources, 0); + gDiffSetSourceImageIDs(inSources, 1); + gDiffSetSourceImageIDs(refSources, 2); + + // define the objects from the diff sources (1-to-1) + psArray *objects = psArrayAllocEmpty(sources->n); + for (int i = 0; i < sources->n; i++) { + pmSource *src = sources->data[i]; + pmPhotObj *obj = pmPhotObjAlloc(); + pmPhotObjAddSource(obj, src); + psArrayAdd (objects, 100, obj); + psFree (obj); + } + + // Match the sources to both inSources and refSources. We want to find the closest source + // in the positive source lists to each source in the diff source list. + // XXX we need to either (a) apply S/N and other filters before we pass the in/ref sources + // to this function, or (b) add those filters in the function + gDiffMatchSources (objects, inSources, radius, minSN); + gDiffMatchSources (objects, refSources, radius, minSN); + + // now mark or flag the sources with matches that meet some criterion + for (int i = 0; i < objects->n; i++) { + pmPhotObj *obj = objects->data[i]; + + // the first entry should be from image 0 (the diff image). + // we should have only 1 of 4 possible results: + // a) only one entry -- the diff source (no positive matches) + // b) 2 matches -- diff and a source from image 1 + // c) 2 matches -- diff and a source from image 2 + // d) 3 matches -- diff and a source from image 1 & 2 + + switch (obj->sources->n) { + case 1: { + pmSource *source = obj->sources->data[0]; + psAssert (source->imageID == 0, "error in pmPhotObj construction: first entry is not image ID 0"); + break; + } + case 2: { + pmSource *source = obj->sources->data[0]; + psAssert (source->imageID == 0, "error in pmPhotObj construction: first entry not 0 (case 2)"); + + pmSource *sourceM1 = obj->sources->data[1]; + psAssert ((sourceM1->imageID == 1) || (sourceM1->imageID == 2), "error in pmPhotObj construction (second entry is not from 1 or 2"); + + // which image gives the match, the positive or negative one? + bool positive = false; + if (subInverse) { + positive = (sourceM1->imageID == 2); + } else { + positive = (sourceM1->imageID == 1); + } + + // check if source and sourceM1 are likely to be the same flux (ie, isolated detection in one image) + if (positive && isfinite(sourceM1->psfMag) && isfinite(source->psfMag)) { + float dMag = sourceM1->psfMag - source->psfMag; + float nSig = fabs(dMag) / hypot(MIN_ERROR, source->psfMagErr); + if (nSig < MAX_OFFSET) { + source->mode2 |= PM_SOURCE_MODE2_DIFF_SELF_MATCH; + } + } + + // only one match. set a flag? + source->mode2 |= PM_SOURCE_MODE2_DIFF_WITH_SINGLE; + if (source->diffStats) { + float SN1 = isfinite(sourceM1->psfMagErr) ? 1.0 / sourceM1->psfMagErr : NAN; + if (positive) { + source->diffStats->SNp = SN1; + source->diffStats->Rp = hypot(source->peak->xf - sourceM1->peak->xf, source->peak->yf - sourceM1->peak->yf); + } else { + source->diffStats->SNm = SN1; + source->diffStats->Rm = hypot(source->peak->xf - sourceM1->peak->xf, source->peak->yf - sourceM1->peak->yf); + } + } + break; + } + case 3: { + pmSource *source = obj->sources->data[0]; + psAssert (source->imageID == 0, "error in pmPhotObj construction: first entry not 0 (case 3)"); + + pmSource *sourceM1 = obj->sources->data[1]; + pmSource *sourceM2 = obj->sources->data[2]; + psAssert ((sourceM1->imageID == 1) || (sourceM1->imageID == 2), "error in pmPhotObj construction (case 3.1)"); + psAssert ((sourceM2->imageID == 1) || (sourceM2->imageID == 2), "error in pmPhotObj construction (case 3.2)"); + psAssert (sourceM1->imageID != sourceM2->imageID, "error in pmPhotObj construction (case 3.3)"); + + // we don't know which matched source (sourceM1 or sourceM2) is the positive detection, + // and the answer depends on both the ID and the subtraction order (in-ref vs ref-in) + // positive is true if sourceM1 is the positive detection + bool positive = false; + if (subInverse) { + positive = (sourceM1->imageID == 2); + } else { + positive = (sourceM1->imageID == 1); + } + + // check for self-detection + float magPSF1 = positive ? sourceM1->psfMag : sourceM2->psfMag; + if (isfinite(magPSF1) && isfinite(source->psfMag)) { + float dMag = magPSF1 - source->psfMag; + float nSig = fabs(dMag) / hypot(MIN_ERROR, source->psfMagErr); + if (nSig < MAX_OFFSET) { + source->mode2 |= PM_SOURCE_MODE2_DIFF_SELF_MATCH; + } + } + + source->mode2 |= PM_SOURCE_MODE2_DIFF_WITH_DOUBLE; + if (source->diffStats) { + float SN1 = isfinite(sourceM1->psfMagErr) ? 1.0 / sourceM1->psfMagErr : NAN; + float SN2 = isfinite(sourceM2->psfMagErr) ? 1.0 / sourceM2->psfMagErr : NAN; + if (positive) { + source->diffStats->SNp = SN1; + source->diffStats->Rp = hypot(source->peak->xf - sourceM1->peak->xf, source->peak->yf - sourceM1->peak->yf); + source->diffStats->SNm = SN2; + source->diffStats->Rm = hypot(source->peak->xf - sourceM2->peak->xf, source->peak->yf - sourceM2->peak->yf); + } else { + source->diffStats->SNp = SN2; + source->diffStats->Rp = hypot(source->peak->xf - sourceM2->peak->xf, source->peak->yf - sourceM2->peak->yf); + source->diffStats->SNm = SN1; + source->diffStats->Rm = hypot(source->peak->xf - sourceM1->peak->xf, source->peak->yf - sourceM1->peak->yf); + } + } + break; + } + default: + psAbort ("error in pmPhotObj construction (unexpected matches)"); + } + } + psFree (objects); + return true; +} + +// XXX based on psphotMatchSourcesToObjects +# define NEXT1 { i++; continue; } +# define NEXT2 { j++; continue; } +bool gDiffMatchSources (psArray *objects, psArray *sources, float RADIUS, float MIN_SN) { + + float dx, dy; + + float RADIUS2 = RADIUS*RADIUS; + + // sort the source list by X + sources = psArraySort (sources, pmSourceSortByX); + objects = psArraySort (objects, pmPhotObjSortByX); + + // match sources to existing objects + + psLogMsg ("psphot", PS_LOG_DETAIL, "attempt to match sources (%ld vs %ld) (SN > %f)", sources->n, objects->n, MIN_SN); + + int nMatch = 0; + int i = 0; + int j = 0; + for (i = j = 0; (i < sources->n) && (j < objects->n); ) { + + pmSource *src = sources->data[i]; + pmPhotObj *obj = objects->data[j]; + + // fprintf (stderr, "%f, %f vs %f, %f\n", src->peak->xf, src->peak->yf, obj->x, obj->y); + + if (!src) NEXT1; + if (!src->peak) NEXT1; + if (!isfinite(src->peak->xf)) NEXT1; + if (!isfinite(src->peak->yf)) NEXT1; + if (!isfinite(src->peak->detValue)) NEXT1; + if (sqrt(src->peak->detValue) < MIN_SN) NEXT1; + + if (!obj) NEXT2; + if (!isfinite(obj->x)) NEXT2; + if (!isfinite(obj->y)) NEXT2; + + dx = src->peak->xf - obj->x; + if (dx < -1.02*RADIUS) NEXT1; + if (dx > +1.02*RADIUS) NEXT2; + + // we are within match range, look for closest match to object j + int Imin = -1; + float Rmin = RADIUS2; + for (int I = i; (dx < +1.02*RADIUS) && (I < sources->n); I++) { + + src = sources->data[I]; + + if (!src) continue; + if (!src->peak) continue; + if (!isfinite(src->peak->xf)) continue; + if (!isfinite(src->peak->yf)) continue; + if (!isfinite(src->peak->detValue)) continue; + if (sqrt(src->peak->detValue) < MIN_SN) continue; + + dx = src->peak->xf - obj->x; + dy = src->peak->yf - obj->y; + + float dr = dx*dx + dy*dy; + if (dr > RADIUS2) continue; + if (dr > Rmin) continue; + Rmin = dr; + Imin = I; + // fprintf (stderr, "j: %d, I: %d, Imin: %d, dr: %f, dx: %f, dy: %f\n", j, I, Imin, dr, dx, dy); + } + + // no match, try next object + if (Imin == -1) { + // fprintf (stderr, "*** missed j: %d, Imin: %d, obj x,y: %f, %f\n", j, Imin, obj->x, obj->y); + j++; + continue; + } + src = sources->data[Imin]; + + // fprintf (stderr, "j: %d, Imin: %d, obj x,y: %f, %f src x,y: %f, %f, SN: %f, ID: %d\n", j, Imin, obj->x, obj->y, src->peak->xf, src->peak->yf, sqrt(src->peak->detValue), src->id); + + // add source to object + pmPhotObjAddSource (obj, src); + nMatch ++; + j++; + } + + psLogMsg ("psphot", PS_LOG_DETAIL, "matched sources (%d matches)", nMatch); + return true; +} + +bool gDiffSetSourceImageIDs (psArray *sources, int imageID) { + + for (int i = 0; i < sources->n; i++) { + pmSource *source = sources->data[i]; + source->imageID = imageID; + } + return true; +} Index: /branches/eam_branches/gDiff.20200202/src/gDiffInputDetections.c =================================================================== --- /branches/eam_branches/gDiff.20200202/src/gDiffInputDetections.c (revision 41413) +++ /branches/eam_branches/gDiff.20200202/src/gDiffInputDetections.c (revision 41413) @@ -0,0 +1,137 @@ +/** @file gDiffInputDetections.c + * @brief: image subtraction based on matched Gaussians profiles + * @author Eugene Magnier @ IfA + * @version $Revision: 1.12 $ + * @date $Date: 2020-01-25 $ + * Copyright 2020 Institute for Astronomy, University of Hawaii + */ + +#include "gDiff.h" + +// internal functions: +pmHDU *getHDUandLevel (pmFPALevel *level, pmFPA *fpa, pmFPAview *view); +pmHDU *setHDUatLevel (pmFPALevel level, pmFPA *fpa, pmFPAview *view, char *extname); + +bool gDiffInputDetections (bool *foundDetections, const char *sourcesName, const char *imageName, gDiffData *data) { + + bool mdok = false; + + psAssert(data, "Require processing data"); + pmConfig *config = data->config; // Configuration + psAssert(config, "Require configuration"); + + // select the view of interest + pmFPAview *view = gDiffViewReadout(); // View to readout + + // point PSPHOT.INPUT at the positive image + pmFPAfile *imageFile = psMetadataLookupPtr(&mdok, config->files, imageName); // Image to photometer + pmFPAfile *sourcesFile = psMetadataLookupPtr(&mdok, config->files, sourcesName); // Place results her + pmFPAfile *photFile = psMetadataLookupPtr(&mdok, config->files, "PSPHOT.INPUT"); // holder for the working image + + // create a psphot working file, so we don't damage the image or affect existing sources + if (!pmFPACopy(photFile->fpa, imageFile->fpa)) { + psError(GDIFF_ERR_CONFIG, false, "Unable to copy FPA for photometry"); + psFree(view); + return false; + } + + pmFPALevel hduLevel = PM_FPA_LEVEL_NONE; + pmHDU *imageHDU = getHDUandLevel (&hduLevel, imageFile->fpa, view); + + if (imageHDU) { + pmFPALevel outLevel = PM_FPA_LEVEL_NONE; + pmHDU *sourcesHDU = getHDUandLevel (&outLevel, sourcesFile->fpa, view); + if (!sourcesHDU) { + sourcesHDU = setHDUatLevel (hduLevel, sourcesFile->fpa, view, imageHDU->extname); + } + + if (!sourcesHDU->header) { + sourcesHDU->header = psMetadataAlloc(); + } + if (!psMetadataCopy(sourcesHDU->header, imageHDU->header)) { + psError(GDIFF_ERR_PROG, false, "Unable to copy header"); + psFree(view); + return false; + } + } + + if (!psphotReadout(config, view, "PSPHOT.INPUT")) { + psErrorStackPrint(stderr, "Unable to perform photometry on image"); + psWarning("Unable to perform photometry on image --- suspect bad data quality."); + gDiffDataQuality(data, psErrorCodeLast(), GDIFF_FILES_PHOT_SUB | GDIFF_FILES_PHOT_INV); + } + // save the outputs on GDIFF.POS1.SOURCES + if (!psphotCopyResults (foundDetections, sourcesFile, photFile, view)) { + psError(GDIFF_ERR_PROG, false, "Unable to copy psphot outputs"); + psFree(view); + return false; + } + // if no sources were found here, we report that back and let them handle it + + psFree(view); + return true; +} + +/**** internal functions ****/ + +pmHDU *getHDUandLevel (pmFPALevel *level, pmFPA *fpa, pmFPAview *view) { + + pmHDU *hdu = NULL; + *level = PM_FPA_LEVEL_NONE; + + pmCell *cell = pmFPAviewThisCell(view, fpa); + if (cell) { + hdu = cell->hdu; // cell + if (hdu) { + *level = PM_FPA_LEVEL_CELL; + return hdu; + } + } + + pmChip *chip = pmFPAviewThisChip(view, fpa); + if (chip) { + hdu = chip->hdu; // chip + if (hdu) { + *level = PM_FPA_LEVEL_CHIP; + return hdu; + } + } + + hdu = fpa->hdu; // fpa + if (hdu) { + *level = PM_FPA_LEVEL_FPA; + return hdu; + } + return NULL; +} + +pmHDU *setHDUatLevel (pmFPALevel level, pmFPA *fpa, pmFPAview *view, char *extname) { + + switch (level) { + case (PM_FPA_LEVEL_FPA): { + if (!fpa->hdu) { + fpa->hdu = pmHDUAlloc(extname); + } + return fpa->hdu; + } + case (PM_FPA_LEVEL_CHIP): { + pmChip *chip = pmFPAviewThisChip(view, fpa); + if (!chip->hdu) { + chip->hdu = pmHDUAlloc(extname); + } + return chip->hdu; + } + case (PM_FPA_LEVEL_CELL): { + pmCell *cell = pmFPAviewThisCell(view, fpa); + if (!cell->hdu) { + cell->hdu = pmHDUAlloc(extname); + } + return cell->hdu; + } + default: + return NULL; + } + return NULL; +} + + Index: /branches/eam_branches/gDiff.20200202/src/gDiffLoop.c =================================================================== --- /branches/eam_branches/gDiff.20200202/src/gDiffLoop.c (revision 41412) +++ /branches/eam_branches/gDiff.20200202/src/gDiffLoop.c (revision 41413) @@ -1,3 +1,376 @@ +/** @file gDiffLoop.c + * @brief: image subtraction based on matched Gaussians profiles + * @author Eugene Magnier @ IfA + * @version $Revision: 1.12 $ + * @date $Date: 2020-01-25 $ + * Copyright 2020 Institute for Astronomy, University of Hawaii + */ + #include "gDiff.h" + +# define HDUCHECK(A) \ + fprintf (stderr, "%s\n", A); \ + listhdu (config, "GDIFF.OUTPUT"); listhdu (config, "PSPHOT.INPUT"); + +// internal functions: +bool dumpout(pmConfig *config, char *name); +bool listhdu(pmConfig *config, char *name); + +bool gDiffLoop(gDiffData *data) +{ + bool mdok = false; + bool success = true; + + psAssert(data, "Require processing data"); + pmConfig *config = data->config; // Configuration + psAssert(config, "Require configuration."); + + HDUCHECK ("0a: start gDiffLoop"); + + psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, GDIFF_RECIPE); // Recipe for ppSim + psAssert(recipe, "We checked this earlier, so it should be here."); + bool noConvolve = psMetadataLookupBool(&mdok, recipe, "NOCONVOLVE"); // Do not use convolved images. + + pmConfigCamerasCull(config, NULL); + pmConfigRecipesCull(config, "GDIFF,PPSTATS,PSPHOT,PSASTRO,MASKS,JPEG"); + + pmFPAfile *input = psMetadataLookupPtr(NULL, config->files, "GDIFF.INPUT"); + pmFPAfile *reference = psMetadataLookupPtr(NULL, config->files, "GDIFF.REF"); + pmFPAfile *output = psMetadataLookupPtr(NULL, config->files, "GDIFF.OUTPUT"); + psAssert(input && reference && output, "Require files"); + + // load input files : we are required to load a complete set of files that describe the input. + // we are not able to process subdivisions of an image (e.g., cells vs chips) + if (!gDiffFilesIterateDown(config, GDIFF_FILES_INPUT | GDIFF_FILES_CONV)) { + psError(GDIFF_ERR_IO, false, "Unable to load files."); + return false; + } + + psTimerStart("GDIFF_MATCH"); + + HDUCHECK ("0b: gDiffSetMasks"); + + if (!gDiffSetMasks(config)) { + psError(psErrorCodeLast(), false, "Unable to set masks."); + return false; + } + + HDUCHECK ("0b.1: gDiffInputDetections"); + if (data->forcedPhot1) { + bool foundDetections = false; + if (!gDiffInputDetections(&foundDetections, "GDIFF.POS1.SOURCES", "GDIFF.INPUT", data)) { + psError(psErrorCodeLast(), false, "Unable to measure positive detections (1)"); + success = false; + goto ESCAPE; + } + // if nothing was found, don't bother doing the forced photometry below + if (!foundDetections) { + psWarning ("no sources found in positive image 1, skipping forced photometry"); + data->forcedPhot1 = false; + } + } + HDUCHECK ("0b.2: gDiffInputDetections"); + if (data->forcedPhot2) { + // XXX why is this only done for Phot2 and not Phot1 ?? + + // Change the recipe to use a higher nsigma limit and quit after pass1 + psMetadata *psphotRecipe = psMetadataLookupPtr (NULL, config->recipes, PSPHOT_RECIPE); + + psF32 nsigma_peak_save = psMetadataLookupF32 (NULL, psphotRecipe, "PEAKS_NSIGMA_LIMIT"); + char *breakPt_save = psMetadataLookupStr (NULL, psphotRecipe, "BREAK_POINT"); + + psF32 pos2_nsigma_peak = psMetadataLookupF32 (&mdok, psphotRecipe, "PEAKS_POS2_NSIGMA_LIMIT"); + if (!mdok) { + psWarning("PEAKS_POS2_NSIGMA_LIMIT not found in psphotRecipe. Will use 25.\n"); + pos2_nsigma_peak = 25.; + } + psMetadataAddF32(psphotRecipe, PS_LIST_TAIL, "PEAKS_NSIGMA_LIMIT", PS_META_REPLACE, "", pos2_nsigma_peak); + psMetadataAddStr(psphotRecipe, PS_LIST_TAIL, "BREAK_POINT", PS_META_REPLACE, "", "PASS1"); + + bool foundDetections = false; + if (!gDiffInputDetections(&foundDetections, "GDIFF.POS2.SOURCES", "GDIFF.REF", data)) { + psError(psErrorCodeLast(), false, "Unable to measure positive detections (2)"); + psMetadataAddF32(psphotRecipe, PS_LIST_TAIL, "PEAKS_NSIGMA_LIMIT", PS_META_REPLACE, "", nsigma_peak_save); + psMetadataAddStr(psphotRecipe, PS_LIST_TAIL, "BREAK_POINT", PS_META_REPLACE, "", breakPt_save); + success = false; + goto ESCAPE; + } + psMetadataAddF32(psphotRecipe, PS_LIST_TAIL, "PEAKS_NSIGMA_LIMIT", PS_META_REPLACE, "", nsigma_peak_save); + psMetadataAddStr(psphotRecipe, PS_LIST_TAIL, "BREAK_POINT", PS_META_REPLACE, "", breakPt_save); + // if nothing was found, don't bother doing the forced photometry below + if (!foundDetections) { + psWarning ("no sources found in positive image 2, skipping forced photometry"); + data->forcedPhot2 = false; + } + } + + // XXX if it exists, use the POS1, POS2 successs for the FWHMs + + /**** The PSF matching process is performed in this function below. **** + **** This is the core of the PSF-matched image subtraction code. ****/ + + HDUCHECK ("1: before gDiffMatchPSFs"); + + if (!gDiffMatchPSFs(data)) { + psError(psErrorCodeLast(), false, "Unable to match PSFs."); + success = false; + goto ESCAPE; + } + if (data->quality) { + // Can't do anything at all + success = false; + goto ESCAPE; + } + + // the functions in this block are not critical for initial testing, so block them out for now + + // generate the residual stamp grid for visualization + if (!gDiffResidualSampleJpeg(config)) { + psError(psErrorCodeLast(), false, "Unable to update."); + success = false; + goto ESCAPE; + } + + // XXX add in a positive image detection step here (if needed) + // XXX EAM what does this mean??? + + psMetadataAddF32(data->stats, PS_LIST_TAIL, "TIME_MATCH", 0, "Time to match PSFs", psTimerClear("GDIFF_MATCH")); + + // Close input files (freeing up space) : for the noConvolve case, we cannot close the + // inputs since we will use them for subtraction below. wait until later to do the work + if (!noConvolve) { + if (!gDiffFilesIterateUp(config, GDIFF_FILES_INPUT)) { + psError(GDIFF_ERR_IO, false, "Unable to close input files."); + success = false; + goto ESCAPE; + } + } + if (!gDiffLowThreshold(data)) { + psError(psErrorCodeLast(), false, "Unable to threshold images."); + success = false; + goto ESCAPE; + } + + HDUCHECK ("2: after gDiffLowThreshold"); + + // Set up subtraction files + if (!gDiffFilesIterateDown(config, GDIFF_FILES_SUB)) { + psError(GDIFF_ERR_IO, false, "Unable to set up subtraction files."); + success = false; + goto ESCAPE; + } + + if (!gDiffDefineOutput("GDIFF.OUTPUT", config)) { + psError(psErrorCodeLast(), false, "Unable to define output."); + success = false; + goto ESCAPE; + } + + HDUCHECK ("2a: gDiffMakePSF"); + + // the functions in this block are not critical for initial testing, so block them out for now + if (!data->quality && !gDiffMakePSF(data)) { + psError(psErrorCodeLast(), false, "Unable to generate PSF."); + success = false; + goto ESCAPE; + } + + // Now we've got a PSF, blow away detections in case they're confused with real output detections + { + pmFPAview *view = gDiffViewReadout(); // View to readout + pmReadout *out = pmFPAfileThisReadout(config->files, view, "GDIFF.OUTPUT"); + if (psMetadataLookup(out->analysis, "PSPHOT.DETECTIONS")) { + psMetadataRemoveKey(out->analysis, "PSPHOT.DETECTIONS"); + } + psFree(view); + } + + HDUCHECK ("3: before gDiffReadoutSubtract"); + + if (!gDiffReadoutSubtract(config)) { + psError(psErrorCodeLast(), false, "Unable to subtract images."); + success = false; + goto ESCAPE; + } + + // Close convolved files + if (!gDiffFilesIterateUp(config, GDIFF_FILES_PSF | GDIFF_FILES_CONV)) { + psError(GDIFF_ERR_IO, false, "Unable to close input files."); + success = false; + goto ESCAPE; + } + // Close unconvolved input files if we did the noConvolve subtraction + // (otherwise they are closed above) + if (noConvolve && !gDiffFilesIterateUp(config, GDIFF_FILES_INPUT)) { + psError(GDIFF_ERR_IO, false, "Unable to close input files."); + success = false; + goto ESCAPE; + } + + HDUCHECK ("4: before gDiffBackground"); + + // Higher order background subtraction using psphot + if (!gDiffBackground(config)) { + psError(psErrorCodeLast(), false, "Unable to subtract background."); + success = false; + goto ESCAPE; + } + // dumpout(config, "diff.2b.fits"); + + HDUCHECK ("5: before gDiffVarianceRescale"); + + // Perform Variance correction (rescale within a modest range) + if (!gDiffVarianceRescale(config, data)) { + psError(psErrorCodeLast(), false, "Unable to rescale variance."); + success = false; + goto ESCAPE; + } + // dumpout(config, "diff.2c.fits"); + + if (data->quality) { + // Done all we can do up to this point + success = false; + goto ESCAPE; + } + + if (!gDiffFilesIterateDown(config, GDIFF_FILES_PHOT_SUB)) { + psError(GDIFF_ERR_IO, false, "Unable to set up photometry files."); + success = false; + goto ESCAPE; + } + + HDUCHECK ("6: before gDiffReadoutPhotometry"); + + if (!data->quality && !gDiffReadoutPhotometry("GDIFF.OUTPUT", data)) { + psError(psErrorCodeLast(), false, "Unable to perform photometry."); + success = false; + goto ESCAPE; + } + // dumpout(config, "diff.3.fits"); + + HDUCHECK ("7: before gDiffReadoutForcedPhot (1)"); + + // forced photometry for positive image 1 + if (data->forcedPhot1 && !data->quality) { + if (!gDiffReadoutForcedPhot("GDIFF.FORCED1.SOURCES", "GDIFF.OUTPUT", "GDIFF.POS1.SOURCES", data)) { + psError(psErrorCodeLast(), false, "Unable to perform photometry."); + success = false; + goto ESCAPE; + } + } + + // forced photometry for positive image 2 + if (data->forcedPhot2 && !data->quality) { + if (!gDiffReadoutForcedPhot("GDIFF.FORCED2.SOURCES", "GDIFF.OUTPUT", "GDIFF.POS2.SOURCES", data)) { + psError(psErrorCodeLast(), false, "Unable to perform photometry."); + success = false; + goto ESCAPE; + } + } + + if (!gDiffFilesIterateUp(config, GDIFF_FILES_PHOT_SUB)) { + psError(GDIFF_ERR_IO, false, "Unable to set up photometry files."); + success = false; + goto ESCAPE; + } + + // Perform statistics on the cell + if (!gDiffReadoutStats(data)) { + psError(psErrorCodeLast(), false, "Unable to collect statistics"); + success = false; + goto ESCAPE; + } + // Do Mask Stats + { + pmFPAview *view = gDiffViewReadout(); // View to readout + if (!gDiffMaskStats(config, view,data->stats)) { + psError(psErrorCodeLast(), false, "Unable to generate mask statistics"); + success = false; + psFree(view); + goto ESCAPE; + } + psFree(view); + } + // dumpout(config, "diff.4.fits"); + + // generate the binned image used to write the jpeg + if (!gDiffReadoutJpeg(config)) { + psError(psErrorCodeLast(), false, "Unable to update."); + success = false; + goto ESCAPE; + } + + if (data->inverse) { + // Set up inverse subtraction files + if (!gDiffFilesIterateDown(config, GDIFF_FILES_INV)) { + psError(GDIFF_ERR_IO, false, "Unable to set up inverse files."); + success = false; + goto ESCAPE; + } + + if (data->inverse && !gDiffDefineOutput("GDIFF.INVERSE", config)) { + psError(psErrorCodeLast(), false, "Unable to define inverse."); + success = false; + goto ESCAPE; + } + + if (!gDiffReadoutInverse(config)) { + psError(psErrorCodeLast(), false, "Unable to invert images."); + success = false; + goto ESCAPE; + } + + //MEH -- need to also clear out detections or inv.cmf corrupted on update + pmFPAview *viewinv = gDiffViewReadout(); // View to readout + pmReadout *inv = pmFPAfileThisReadout(config->files, viewinv, "GDIFF.INVERSE"); + if (psMetadataLookup(inv->analysis, "PSPHOT.DETECTIONS")) { + psMetadataRemoveKey(inv->analysis, "PSPHOT.DETECTIONS"); + } + psFree(viewinv); + + // Close subtraction files and open inverse photometry files + if (!gDiffFilesIterateUp(config, GDIFF_FILES_SUB)) { + psError(GDIFF_ERR_IO, false, "Unable to close subtraction files."); + success = false; + goto ESCAPE; + } + + if (!gDiffFilesIterateDown(config, GDIFF_FILES_PHOT_INV)) { + psError(GDIFF_ERR_IO, false, "Unable to set up inverse files."); + success = false; + goto ESCAPE; + } + + if (!data->quality && !gDiffReadoutPhotometry("GDIFF.INVERSE", data)) { + psError(psErrorCodeLast(), false, "Unable to perform photometry."); + success = false; + goto ESCAPE; + } + + // Close inverse subtraction files + if (!gDiffFilesIterateUp(config, GDIFF_FILES_INV | GDIFF_FILES_PHOT_INV)) { + psError(GDIFF_ERR_IO, false, "Unable to close subtraction files."); + success = false; + goto ESCAPE; + } + } else { + // Close subtraction files + // dumpout(config, "diff.5.fits"); + if (!gDiffFilesIterateUp(config, GDIFF_FILES_SUB)) { + psError(GDIFF_ERR_IO, false, "Unable to close subtraction files."); + success = false; + goto ESCAPE; + } + } + + ESCAPE: + pmFPAfileDropInternal(config->files, "PSPHOT.BACKGND"); + pmFPAfileDropInternal(config->files, "PSPHOT.BACKMDL"); + pmFPAfileDropInternal(config->files, "PSPHOT.BACKMDL.STDEV"); + + return success; +} + +/**** internal functions ****/ bool dumpout(pmConfig *config, char *name) @@ -10,327 +383,8 @@ } - -bool gDiffLoop(gDiffData *data) +bool listhdu(pmConfig *config, char *name) { - bool mdok = false; - bool success = true; - - psAssert(data, "Require processing data"); - pmConfig *config = data->config; // Configuration - psAssert(config, "Require configuration."); - - psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, GDIFF_RECIPE); // Recipe for ppSim - psAssert(recipe, "We checked this earlier, so it should be here."); - bool noConvolve = psMetadataLookupBool(&mdok, recipe, "NOCONVOLVE"); // Do not use convolved images. - - pmConfigCamerasCull(config, NULL); - pmConfigRecipesCull(config, "GDIFF,PPSTATS,PSPHOT,PSASTRO,MASKS,JPEG"); - - pmFPAfile *input = psMetadataLookupPtr(NULL, config->files, "GDIFF.INPUT"); - pmFPAfile *reference = psMetadataLookupPtr(NULL, config->files, "GDIFF.REF"); - pmFPAfile *output = psMetadataLookupPtr(NULL, config->files, "GDIFF.OUTPUT"); - psAssert(input && reference && output, "Require files"); - - if (!gDiffFilesIterateDown(config, GDIFF_FILES_INPUT | GDIFF_FILES_CONV)) { - psError(GDIFF_ERR_IO, false, "Unable to load files."); - return false; - } - - psTimerStart("GDIFF_MATCH"); - - if (!gDiffSetMasks(config)) { - psError(psErrorCodeLast(), false, "Unable to set masks."); - return false; - } - - if (data->forcedPhot1) { - bool foundDetections = false; - if (!gDiffInputDetections(&foundDetections, "GDIFF.POS1.SOURCES", "GDIFF.INPUT", data)) { - psError(psErrorCodeLast(), false, "Unable to measure positive detections (1)"); - success = false; - goto ESCAPE; - } - // if nothing was found, don't bother doing the forced photometry below - if (!foundDetections) { - psWarning ("no sources found in positive image 1, skipping forced photometry"); - data->forcedPhot1 = false; - } - } - if (data->forcedPhot2) { - // Change the recipe to use a higher nsigma limit and quit after pass1 - psMetadata *psphotRecipe = psMetadataLookupPtr (NULL, config->recipes, PSPHOT_RECIPE); - - psF32 nsigma_peak_save = psMetadataLookupF32 (NULL, psphotRecipe, "PEAKS_NSIGMA_LIMIT"); - char *breakPt_save = psMetadataLookupStr (NULL, psphotRecipe, "BREAK_POINT"); - - psF32 pos2_nsigma_peak = psMetadataLookupF32 (&mdok, psphotRecipe, "PEAKS_POS2_NSIGMA_LIMIT"); - if (!mdok) { - psWarning("PEAKS_POS2_NSIGMA_LIMIT not found in psphotRecipe. Will use 25.\n"); - pos2_nsigma_peak = 25.; - } - psMetadataAddF32(psphotRecipe, PS_LIST_TAIL, "PEAKS_NSIGMA_LIMIT", PS_META_REPLACE, "", pos2_nsigma_peak); - psMetadataAddStr(psphotRecipe, PS_LIST_TAIL, "BREAK_POINT", PS_META_REPLACE, "", "PASS1"); - - bool foundDetections = false; - if (!gDiffInputDetections(&foundDetections, "GDIFF.POS2.SOURCES", "GDIFF.REF", data)) { - psError(psErrorCodeLast(), false, "Unable to measure positive detections (2)"); - psMetadataAddF32(psphotRecipe, PS_LIST_TAIL, "PEAKS_NSIGMA_LIMIT", PS_META_REPLACE, "", nsigma_peak_save); - psMetadataAddStr(psphotRecipe, PS_LIST_TAIL, "BREAK_POINT", PS_META_REPLACE, "", breakPt_save); - success = false; - goto ESCAPE; - } - psMetadataAddF32(psphotRecipe, PS_LIST_TAIL, "PEAKS_NSIGMA_LIMIT", PS_META_REPLACE, "", nsigma_peak_save); - psMetadataAddStr(psphotRecipe, PS_LIST_TAIL, "BREAK_POINT", PS_META_REPLACE, "", breakPt_save); - // if nothing was found, don't bother doing the forced photometry below - if (!foundDetections) { - psWarning ("no sources found in positive image 2, skipping forced photometry"); - data->forcedPhot2 = false; - } - } - - // XXX if it exists, use the POS1, POS2 successs for the FWHMs - if (!gDiffMatchPSFs(data)) { - psError(psErrorCodeLast(), false, "Unable to match PSFs."); - success = false; - goto ESCAPE; - } - - if (data->quality) { - // Can't do anything at all - success = false; - goto ESCAPE; - } - // generate the residual stamp grid for visualization - if (!gDiffResidualSampleJpeg(config)) { - psError(psErrorCodeLast(), false, "Unable to update."); - success = false; - goto ESCAPE; - } - - // XXX add in a positive image detection step here (if needed) - - - psMetadataAddF32(data->stats, PS_LIST_TAIL, "TIME_MATCH", 0, "Time to match PSFs", - psTimerClear("GDIFF_MATCH")); - - // Close input files (freeing up space) : for the noConvolve case, we cannot close the - // inputs since we will use them for subtraction below. wait until later to do the work - if (!noConvolve) { - if (!gDiffFilesIterateUp(config, GDIFF_FILES_INPUT)) { - psError(GDIFF_ERR_IO, false, "Unable to close input files."); - success = false; - goto ESCAPE; - } - } - if (!gDiffLowThreshold(data)) { - psError(psErrorCodeLast(), false, "Unable to threshold images."); - success = false; - goto ESCAPE; - } - - // Set up subtraction files - if (!gDiffFilesIterateDown(config, GDIFF_FILES_SUB)) { - psError(GDIFF_ERR_IO, false, "Unable to set up subtraction files."); - success = false; - goto ESCAPE; - } - - if (!gDiffDefineOutput("GDIFF.OUTPUT", config)) { - psError(psErrorCodeLast(), false, "Unable to define output."); - success = false; - goto ESCAPE; - } - - if (!data->quality && !gDiffMakePSF(data)) { - psError(psErrorCodeLast(), false, "Unable to generate PSF."); - success = false; - goto ESCAPE; - } - - // Now we've got a PSF, blow away detections in case they're confused with real output detections - { - pmFPAview *view = gDiffViewReadout(); // View to readout - pmReadout *out = pmFPAfileThisReadout(config->files, view, "GDIFF.OUTPUT"); - if (psMetadataLookup(out->analysis, "PSPHOT.DETECTIONS")) { - psMetadataRemoveKey(out->analysis, "PSPHOT.DETECTIONS"); - } - psFree(view); - } - - if (!gDiffReadoutSubtract(config)) { - psError(psErrorCodeLast(), false, "Unable to subtract images."); - success = false; - goto ESCAPE; - } - // dumpout(config, "diff.1.fits"); - - // Close convolved files - if (!gDiffFilesIterateUp(config, GDIFF_FILES_PSF | GDIFF_FILES_CONV)) { - psError(GDIFF_ERR_IO, false, "Unable to close input files."); - success = false; - goto ESCAPE; - } - // dumpout(config, "diff.2a.fits"); - if (noConvolve) { - if (!gDiffFilesIterateUp(config, GDIFF_FILES_INPUT)) { - psError(GDIFF_ERR_IO, false, "Unable to close input files."); - success = false; - goto ESCAPE; - } - } - - // Higher order background subtraction using psphot - if (!gDiffBackground(config)) { - psError(psErrorCodeLast(), false, "Unable to subtract background."); - success = false; - goto ESCAPE; - } - // dumpout(config, "diff.2b.fits"); - - // Perform Variance correction (rescale within a modest range) - if (!gDiffVarianceRescale(config, data)) { - psError(psErrorCodeLast(), false, "Unable to rescale variance."); - success = false; - goto ESCAPE; - } - // dumpout(config, "diff.2c.fits"); - - if (data->quality) { - // Done all we can do up to this point - success = false; - goto ESCAPE; - } - - if (!gDiffFilesIterateDown(config, GDIFF_FILES_PHOT_SUB)) { - psError(GDIFF_ERR_IO, false, "Unable to set up photometry files."); - success = false; - goto ESCAPE; - } - - if (!data->quality && !gDiffReadoutPhotometry("GDIFF.OUTPUT", data)) { - psError(psErrorCodeLast(), false, "Unable to perform photometry."); - success = false; - goto ESCAPE; - } - // dumpout(config, "diff.3.fits"); - - // forced photometry for positive image 1 - if (data->forcedPhot1 && !data->quality) { - if (!gDiffReadoutForcedPhot("GDIFF.FORCED1.SOURCES", "GDIFF.OUTPUT", "GDIFF.POS1.SOURCES", data)) { - psError(psErrorCodeLast(), false, "Unable to perform photometry."); - success = false; - goto ESCAPE; - } - } - - // forced photometry for positive image 2 - if (data->forcedPhot2 && !data->quality) { - if (!gDiffReadoutForcedPhot("GDIFF.FORCED2.SOURCES", "GDIFF.OUTPUT", "GDIFF.POS2.SOURCES", data)) { - psError(psErrorCodeLast(), false, "Unable to perform photometry."); - success = false; - goto ESCAPE; - } - } - - if (!gDiffFilesIterateUp(config, GDIFF_FILES_PHOT_SUB)) { - psError(GDIFF_ERR_IO, false, "Unable to set up photometry files."); - success = false; - goto ESCAPE; - } - - // Perform statistics on the cell - if (!gDiffReadoutStats(data)) { - psError(psErrorCodeLast(), false, "Unable to collect statistics"); - success = false; - goto ESCAPE; - } - // Do Mask Stats - { - pmFPAview *view = gDiffViewReadout(); // View to readout - if (!gDiffMaskStats(config, view,data->stats)) { - psError(psErrorCodeLast(), false, "Unable to generate mask statistics"); - success = false; - psFree(view); - goto ESCAPE; - } - psFree(view); - } - // dumpout(config, "diff.4.fits"); - - // generate the binned image used to write the jpeg - if (!gDiffReadoutJpeg(config)) { - psError(psErrorCodeLast(), false, "Unable to update."); - success = false; - goto ESCAPE; - } - - if (data->inverse) { - // Set up inverse subtraction files - if (!gDiffFilesIterateDown(config, GDIFF_FILES_INV)) { - psError(GDIFF_ERR_IO, false, "Unable to set up inverse files."); - success = false; - goto ESCAPE; - } - - if (data->inverse && !gDiffDefineOutput("GDIFF.INVERSE", config)) { - psError(psErrorCodeLast(), false, "Unable to define inverse."); - success = false; - goto ESCAPE; - } - - if (!gDiffReadoutInverse(config)) { - psError(psErrorCodeLast(), false, "Unable to invert images."); - success = false; - goto ESCAPE; - } - - //MEH -- need to also clear out detections or inv.cmf corrupted on update - pmFPAview *viewinv = gDiffViewReadout(); // View to readout - pmReadout *inv = pmFPAfileThisReadout(config->files, viewinv, "GDIFF.INVERSE"); - if (psMetadataLookup(inv->analysis, "PSPHOT.DETECTIONS")) { - psMetadataRemoveKey(inv->analysis, "PSPHOT.DETECTIONS"); - } - psFree(viewinv); - - // Close subtraction files and open inverse photometry files - if (!gDiffFilesIterateUp(config, GDIFF_FILES_SUB)) { - psError(GDIFF_ERR_IO, false, "Unable to close subtraction files."); - success = false; - goto ESCAPE; - } - - if (!gDiffFilesIterateDown(config, GDIFF_FILES_PHOT_INV)) { - psError(GDIFF_ERR_IO, false, "Unable to set up inverse files."); - success = false; - goto ESCAPE; - } - - if (!data->quality && !gDiffReadoutPhotometry("GDIFF.INVERSE", data)) { - psError(psErrorCodeLast(), false, "Unable to perform photometry."); - success = false; - goto ESCAPE; - } - - // Close inverse subtraction files - if (!gDiffFilesIterateUp(config, GDIFF_FILES_INV | GDIFF_FILES_PHOT_INV)) { - psError(GDIFF_ERR_IO, false, "Unable to close subtraction files."); - success = false; - goto ESCAPE; - } - } else { - // Close subtraction files - // dumpout(config, "diff.5.fits"); - if (!gDiffFilesIterateUp(config, GDIFF_FILES_SUB)) { - psError(GDIFF_ERR_IO, false, "Unable to close subtraction files."); - success = false; - goto ESCAPE; - } - } - - ESCAPE: - pmFPAfileDropInternal(config->files, "PSPHOT.BACKGND"); - pmFPAfileDropInternal(config->files, "PSPHOT.BACKMDL"); - pmFPAfileDropInternal(config->files, "PSPHOT.BACKMDL.STDEV"); - - return success; + pmFPAfile *file = psMetadataLookupPtr(NULL, config->files, name); + fprintf (stderr, "HDU: %s, %llx\n", name, (long long unsigned int) file->fpa->hdu); + return true; } Index: /branches/eam_branches/gDiff.20200202/src/gDiffMakePSF.c =================================================================== --- /branches/eam_branches/gDiff.20200202/src/gDiffMakePSF.c (revision 41413) +++ /branches/eam_branches/gDiff.20200202/src/gDiffMakePSF.c (revision 41413) @@ -0,0 +1,190 @@ +/** @file gDiffMakePSF.c + * @brief: image subtraction based on matched Gaussians profiles + * @author Eugene Magnier @ IfA + * @version $Revision: 1.12 $ + * @date $Date: 2020-01-25 $ + * Copyright 2020 Institute for Astronomy, University of Hawaii + */ + +#include "gDiff.h" + +psArray *gDiffSelectPSFSources(psArray *sources); + +bool gDiffMakePSF(gDiffData *data) +{ + psAssert(data, "Require processing data"); + pmConfig *config = data->config; // Configuration + psAssert(config, "Require configuration"); + + if (!data->photometry) { + return true; + } + + psTimerStart("GDIFF_PHOT"); + + bool reverse = psMetadataLookupBool(NULL, config->arguments, "REVERSE"); // Reverse sense of subtraction? + + bool mdok = false; // Status of MD lookup + pmReadout *minuend = NULL; // Image that will be positive following subtraction + pmFPAfile *minuendFile = NULL; // File for minuend image + pmFPAview *view = gDiffViewReadout(); // View to readout + + psMetadata *recipe = psMetadataLookupMetadata(&mdok, config->recipes, GDIFF_RECIPE); // Recipe for gDiff + bool noConvolve = psMetadataLookupBool(&mdok, recipe, "NOCONVOLVE"); // Do not use convolved images. + + if (noConvolve) { + // if we do not convolve, we need to copy the detections to the image for analysis + pmReadout *psfSourcesRO = NULL; // readout containing loaded sources for psf model + psWarning("Not using Convolved images because NOCONVOLVE is TRUE\n"); + if (reverse) { + minuend = pmFPAfileThisReadout(config->files, view, "GDIFF.REF"); + minuendFile = psMetadataLookupPtr(&mdok, config->files, "GDIFF.REF"); + psfSourcesRO = pmFPAfileThisReadout(config->files, view, "GDIFF.REF.SOURCES"); + } else { + minuend = pmFPAfileThisReadout(config->files, view, "GDIFF.INPUT"); + minuendFile = psMetadataLookupPtr(&mdok, config->files, "GDIFF.INPUT"); + psfSourcesRO = pmFPAfileThisReadout(config->files, view, "GDIFF.INPUT.SOURCES"); + } + psAssert (psfSourcesRO, "missing readout with sources"); + pmDetections *psfDetections = psMetadataLookupPtr(&mdok, psfSourcesRO->analysis, "PSPHOT.DETECTIONS"); + psMetadataAddPtr(minuend->analysis, PS_LIST_TAIL, "PSPHOT.DETECTIONS", PS_META_REPLACE | PS_DATA_UNKNOWN, "Merged source list", psfDetections); + } else { + printf("Using Convolved images because NOCONVOLVE is FALSE\n"); + if (reverse) { + minuend = pmFPAfileThisReadout(config->files, view, "GDIFF.REF.CONV"); + minuendFile = psMetadataLookupPtr(&mdok, config->files, "GDIFF.REF.CONV"); + } else { + minuend = pmFPAfileThisReadout(config->files, view, "GDIFF.INPUT.CONV"); + minuendFile = psMetadataLookupPtr(&mdok, config->files, "GDIFF.INPUT.CONV"); + } + } + + pmFPAfile *photFile = psMetadataLookupPtr(&mdok, config->files, "PSPHOT.INPUT"); // Photometry file + if (!pmFPACopy(photFile->fpa, minuendFile->fpa)) { + psError(GDIFF_ERR_CONFIG, false, "Unable to copy FPA for photometry"); + psFree(view); + return false; + } + + pmReadout *photRO = pmFPAviewThisReadout(view, photFile->fpa); // Readout to photometer + + // 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 + pmDetections *detections = psMetadataLookupPtr(&mdok, minuend->analysis, "PSPHOT.DETECTIONS"); + if (!detections || !detections->allSources) { + psError(GDIFF_ERR_CONFIG, true, "No sources from which to determine PSF."); + psFree(view); + return false; + } + psArray *sources = detections->allSources; + + // XXX filter sources? limit the total number and return only brighter objects? + // use flags to toss totally bogus entries? + psArray *goodSources = gDiffSelectPSFSources (sources); + + if (!psphotReadoutFindPSF(config, view, "PSPHOT.INPUT", goodSources)) { + // This is likely a data quality issue + // XXX Split into multiple cases using error codes? + psErrorStackPrint(stderr, "Unable to determine PSF"); + psWarning("Unable to determine PSF --- suspect bad data quality."); + gDiffDataQuality(data, PSPHOT_ERR_PSF, GDIFF_FILES_PHOT_SUB | GDIFF_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 (!gDiffCopyPSF(psfFile, photFile, view)) { + psErrorStackPrint(stderr, "PSF was not generated"); + psWarning("PSF was not generated --- suspect bad data quality."); + gDiffDataQuality(data, psErrorCodeLast(), GDIFF_FILES_PHOT_SUB | GDIFF_FILES_PHOT_INV); + } + + // Get rid of the generated header; it will be regenerated by the real photometry run + psMetadataRemoveKey(photRO->analysis, "PSPHOT.HEADER"); + + psFree(goodSources); + psFree(view); + + return true; +} + +bool gDiffCopyPSF(pmFPAfile *output, pmFPAfile *input, pmFPAview *view) +{ + 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); + psFree(outputRO); // I have a copy on the outputCell + } + + // Copy the PSF-related data + psMetadataIterator *iter = psMetadataIteratorAlloc(inputRO->analysis, PS_LIST_HEAD, "^PSF.*"); + psMetadataItem *item; // Item from iteration + while ((item = psMetadataGetAndIncrement(iter))) { + psMetadataAddItem(outputRO->analysis, item, PS_LIST_TAIL, PS_META_REPLACE); + } + psFree(iter); + + // 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 gDiffMakePSF", psf); + + return true; +} + + +// XXX hardwired MAX for now +# define MAX_NPSF 500 + +psArray *gDiffSelectPSFSources(psArray *sources){ + + sources = psArraySort (sources, pmSourceSortByFlux); + + 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: /branches/eam_branches/gDiff.20200202/src/gDiffMaskStats.c =================================================================== --- /branches/eam_branches/gDiff.20200202/src/gDiffMaskStats.c (revision 41413) +++ /branches/eam_branches/gDiff.20200202/src/gDiffMaskStats.c (revision 41413) @@ -0,0 +1,57 @@ +/** @file gDiffReadoutForcedPhot.c + * @brief: run forced photometry on the specified image (i.e., the input images) + * @author Eugene Magnier @ IfA + * @version $Revision: 1.12 $ + * @date $Date: 2020-01-25 $ + * Copyright 2020 Institute for Astronomy, University of Hawaii + */ + +#include "gDiff.h" + +// XXX EAM : should't this function be called as part of gDiffReadoutStats?? + +bool gDiffMaskStats(pmConfig *config, pmFPAview *view, psMetadata *stats) +{ + PS_ASSERT_PTR_NON_NULL(view, false); + PS_ASSERT_PTR_NON_NULL(config, false); + + pmReadout *readout = pmFPAfileThisReadout(config->files, view, "GDIFF.OUTPUT"); + if (!stats || !readout || !readout->data_exists) { + // Nothing to process + return(true); + } + + bool status; + + psMetadata *recipe = psMetadataLookupMetadata(&status, config->recipes, GDIFF_RECIPE); + psImageMaskType staticMaskVal = psMetadataLookupImageMask(&status, recipe, "MASKSTAT.STATIC"); + psImageMaskType magicMaskVal = psMetadataLookupImageMask(&status, recipe, "MASKSTAT.MAGIC"); + psImageMaskType dynamicMaskVal = psMetadataLookupImageMask(&status, recipe, "MASKSTAT.DYNAMIC"); + psImageMaskType advisoryMaskVal = psMetadataLookupImageMask(&status, recipe, "MASKSTAT.ADVISORY"); + + psS32 Npix_valid = 0; + psS32 Npix_static = 0; + psS32 Npix_magic = 0; + psS32 Npix_dynamic = 0; + psS32 Npix_advisory = 0; + + psImage *mask = readout->mask; // Mask of interest; + if (!pmSingleImageMaskStats(mask,&Npix_valid,&Npix_static,&Npix_magic, + &Npix_dynamic,&Npix_advisory, + staticMaskVal,magicMaskVal, + dynamicMaskVal,advisoryMaskVal)) { + psError(PS_ERR_UNKNOWN, false, "Unable to calculate masks for readout."); + return(false); + } + psMetadataAddS32(stats, PS_LIST_TAIL,"MASKFRAC_NPIX", 0, + "Number of valid pixels", Npix_valid); + psMetadataAddF32(stats,PS_LIST_TAIL, "MASKFRAC_STATIC", 0, + "Fraction of pixels statically masked", (float) Npix_static / Npix_valid); + psMetadataAddF32(stats,PS_LIST_TAIL, "MASKFRAC_DYNAMIC", 0, + "Fraction of pixels dynamically masked", (float) Npix_dynamic / Npix_valid); + psMetadataAddF32(stats,PS_LIST_TAIL, "MASKFRAC_MAGIC", 0, + "Fraction of pixels magically masked", (float) Npix_magic / Npix_valid); + psMetadataAddF32(stats,PS_LIST_TAIL, "MASKFRAC_ADVISORY", 0, + "Fraction of pixels masked as an advisory", (float) Npix_advisory / Npix_valid); + return(true); +} Index: /branches/eam_branches/gDiff.20200202/src/gDiffMatchPSFs.c =================================================================== --- /branches/eam_branches/gDiff.20200202/src/gDiffMatchPSFs.c (revision 41413) +++ /branches/eam_branches/gDiff.20200202/src/gDiffMatchPSFs.c (revision 41413) @@ -0,0 +1,552 @@ +/** @file gDiffMatchPSFs.c + * @brief: image subtraction based on matched Gaussians profiles + * @author Eugene Magnier @ IfA + * @version $Revision: 1.12 $ + * @date $Date: 2020-01-25 $ + * Copyright 2020 Institute for Astronomy, University of Hawaii + */ + +/***** this is the core of the PSF matching process + ***** this file needs to be carefully vetted + *****/ + +#include "gDiff.h" + +#define COVAR_FRAC 0.01 // Truncation fraction for covariance matrix + +// internal functions: +static void normaliseRegion(psImage *image, const psRegion *region, float norm); +static float subImagePSF(gDiffData *data, const pmReadout *ro, psArray *sources); +static bool subScaleKernel(gDiffData *data, psVector *kernelWidths, int *kernelSize, int *stampSize); +pmSubtractionMode subModeFromString (char *string); + +// generate the matched PSFs +bool gDiffMatchPSFs(gDiffData *data) +{ + bool mdok = false; + + psAssert(data, "Require processing data"); + pmConfig *config = data->config; // Configuration + psAssert(config, "Require configuration"); + + // Look up recipe values + psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, GDIFF_RECIPE); // Recipe for ppSim + psAssert(recipe, "We checked this earlier, so it should be here."); + + // XXX EAM : for better code consistency downstream, the NOCONVOLVE case should copy the + // unconvolved images to their convolved counterparts + + bool noConvolve = psMetadataLookupBool(&mdok, recipe, "NOCONVOLVE"); // Do not use convolved images. + if (noConvolve) { + psWarning("not matching PSFs because NOCONVOLVE is TRUE\n"); + return true; + } + + pmFPAview *view = gDiffViewReadout(); // View to readout + + // Input images + pmReadout *inRO = pmFPAfileThisReadout(config->files, view, "GDIFF.INPUT"); // Input readout + pmReadout *refRO = pmFPAfileThisReadout(config->files, view, "GDIFF.REF"); // Reference readout + + // Output image holders + pmReadout *inConv = pmFPAfileThisReadout(config->files, view, "GDIFF.INPUT.CONV"); // Input convolved + if (!inConv) { + pmCell *cell = pmFPAfileThisCell(config->files, view, "GDIFF.INPUT.CONV"); // Cell for convolved input + inConv = pmReadoutAlloc(cell); + psFree(inConv); + } + pmReadout *refConv = pmFPAfileThisReadout(config->files, view, "GDIFF.REF.CONV"); // Reference convolved + if (!refConv) { + pmCell *cell = pmFPAfileThisCell(config->files, view, "GDIFF.REF.CONV"); // Cell for convolved ref. + refConv = pmReadoutAlloc(cell); + psFree(refConv); + } + + // Load pre-calculated kernel, if available + pmReadout *kernelRO = pmFPAfileThisReadout(config->files, view, "GDIFF.OUTPUT.KERNELS"); // RO with kernel + + // Sources in image, used for stamps: these must be loaded from previous analysis stages + pmReadout *inSourceRO = pmFPAfileThisReadout(config->files, view, "GDIFF.INPUT.SOURCES"); + pmReadout *refSourceRO = pmFPAfileThisReadout(config->files, view, "GDIFF.REF.SOURCES"); + + // get the detection tables from the readouts above + 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); + + // XXX MEH changed to get only match (true to cull single), no apparent need of unmatched sources and can cause trouble + detections->allSources = pmSourceMatchMerge(lists, radius, true); + psFree(lists); + if (!detections->allSources) { + psFree(detections); + int error = psErrorCodeLast(); // Error code + if (error == PM_ERR_OBJECTS) { + psErrorStackPrint(stderr, "Unable to match source lists"); + psWarning("Unable to match source lists --- suspect bad data quality."); + gDiffDataQuality(data, error, GDIFF_FILES_ALL); + return true; + } else { + // XXX is this a data quality error as well? (how else can pmSourceMatchMerge fail?) + psError(error, false, "Unable to merge source lists"); + return false; + } + } + } + + /** XXX EAM: the original design allowed detections to be supplied from either image. + ** however, as noted above, unmatched detections can cause poor PSF matching. + ** Should the analysis *require* valid matched detections from both images? + ** if so, !detections should raise a 'quality' error + **/ + + 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 + + /** XXX EAM: this is a huge number of parameters which we grab all at once + ** rationalize this somehow? + **/ + + 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 *typeStr = psMetadataLookupStr(NULL, recipe, "KERNEL.TYPE"); // Kernel type + psAssert(typeStr, "We put it here in gDiffArguments.c"); + pmSubtractionKernelsType type = pmSubtractionKernelsTypeFromString(typeStr); // Type of kernel + if (type == PM_SUBTRACTION_KERNEL_NONE) { + psError(GDIFF_ERR_ARGUMENTS, 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 iter = psMetadataLookupS32(NULL, recipe, "ITER"); // Rejection iterations + 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 + float poorFrac = psMetadataLookupF32(&mdok, recipe, "POOR.FRACTION"); // Fraction for "poor" + + // 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 + psVector *optWidths = NULL; // Vector with FWHMs for optimum search + if (optimum) { + optWidths = psVectorCreate(optWidths, optMin, optMax, optStep, PS_TYPE_F32); + } + int optOrder = psMetadataLookupS32(&mdok, recipe, "OPTIMUM.ORDER"); // Order for search + float optThresh = psMetadataLookupF32(&mdok, recipe, "OPTIMUM.TOL"); // Tolerance for search + + psImageMaskType maskVal = pmConfigMaskGet("MASK.VALUE", config); // Bits to mask in inputs + psImageMaskType maskPoor = pmConfigMaskGet("CONV.POOR", config); // Bits to mask for poor pixels + psImageMaskType maskBad = pmConfigMaskGet("BLANK", config); // Bits to mask for bad pixels + + // we have three things to control the convolution choice: + + // 1) in the recipe, the keyword DUAL can be TRUE or FALSE. if TRUE, DUAL is used. If not + // DUAL, then we are doing SINGLE convolution (but see CONVOLVE.TARGET below). + + // 2) if the -convolve option is given on the command line, the specified input (1 or 2) is + // the one convolved. Allowed values are 1 or 2. (XXX EAM: what happens if DUAL is TRUE and -convolve has + // a value of 1 or 2? Is -convolve equal to CONVOLVE.TARGET? + + // 3) if the -convolve option is NOT given, then the value of the recipe keyword + // CONVOLVE.TARGET is used. this may have the value DUAL, AUTO, SINGLE1, SINGLE2, 1, 2. + // thus, DUAL convolution can be turned on with either DUAL = T or CONVOLVE.TARGET = DUAL + + bool dual = psMetadataLookupBool(&mdok, recipe, "DUAL"); // Dual convolution? + pmSubtractionMode subMode; // Subtraction mode + if (dual) { + subMode = PM_SUBTRACTION_MODE_DUAL; + } else { + char *convolveName = psMetadataLookupStr(&mdok, recipe, "CONVOLVE.TARGET"); // recipe value for target + int convolve = psMetadataLookupS32(&mdok, config->arguments, "-convolve"); // override with command-line option + switch (convolve) { + case 0: + // convolve is 0 if it is not supplied on the command line + subMode = subModeFromString(convolveName); + break; + case 1: + subMode = PM_SUBTRACTION_MODE_1; + break; + case 2: + subMode = PM_SUBTRACTION_MODE_2; + break; + default: + psError(GDIFF_ERR_ARGUMENTS, false, "Invalid value for -convolve"); + return false; + } + } + + if (!subScaleKernel(data, widths, &size, &footprint)) { + psError(GDIFF_ERR_DATA, false, "Unable to scale kernel parameters"); + return false; + } + + int threads = psMetadataLookupS32(NULL, config->arguments, "NTHREADS"); // Number of threads + if (threads > 0) { + pmSubtractionThreadsInit(); + } + + if (inRO->covariance) { + psKernel *truncated = psImageCovarianceTruncate(inRO->covariance, COVAR_FRAC); + psFree(inRO->covariance); + inRO->covariance = truncated; + } + if (refRO->covariance) { + psKernel *truncated = psImageCovarianceTruncate(refRO->covariance, COVAR_FRAC); + psFree(refRO->covariance); + refRO->covariance = truncated; + } + + // Data doesn't exist until it's created in pmSubtraction... + inConv->data_exists = inConv->parent->data_exists = inConv->parent->parent->data_exists = false; + refConv->data_exists = refConv->parent->data_exists = refConv->parent->parent->data_exists = false; + + /***** this function call below is the core of gDiffMatchPSFs *****/ + + // Match the PSFs + bool success = false; // Operation was successful? + if (kernelRO) { + success = pmSubtractionMatchPrecalc(inConv, refConv, inRO, refRO, kernelRO->analysis, + stride, kernelErr, covarFrac, maskVal, maskBad, maskPoor, + poorFrac, badFrac); + } else { + success = pmSubtractionMatch(inConv, refConv, inRO, refRO, footprint, stride, regionSize, + spacing, threshold, detections ? detections->allSources : NULL, + data->stamps, type, size, order, widths, orders, inner, ringsOrder, + binning, penalty, optimum, optWidths, optOrder, optThresh, iter, rej, + normFrac, sysErr, skyErr, kernelErr, covarFrac, + maskVal, maskBad, maskPoor, poorFrac, badFrac, subMode); + } + +# ifdef TESTING + // XXX for testing + psphotSaveImage (NULL, refRO->image, "refRO.im.fits"); + psphotSaveImage (NULL, refRO->variance, "refRO.wt.fits"); + psphotSaveImage (NULL, refRO->mask, "refRO.mk.fits"); + + psphotSaveImage (NULL, inRO->image, "inRO.im.fits"); + psphotSaveImage (NULL, inRO->variance, "inRO.wt.fits"); + psphotSaveImage (NULL, inRO->mask, "inRO.mk.fits"); + + psphotSaveImage (NULL, inConv->image, "inConv.im.fits"); + psphotSaveImage (NULL, inConv->variance, "inConv.wt.fits"); + psphotSaveImage (NULL, inConv->mask, "inConv.mk.fits"); + + psphotSaveImage (NULL, refConv->image, "refConv.im.fits"); + psphotSaveImage (NULL, refConv->variance, "refConv.wt.fits"); + psphotSaveImage (NULL, refConv->mask, "refConv.mk.fits"); +# endif + + psFree(optWidths); + pmSubtractionThreadsFinalize(); + + if (!success) { + int error = psErrorCodeLast(); // Error code + switch (error) { + case PM_ERR_STAMPS: + psErrorStackPrint(stderr, "Unable to find stamps"); + psWarning("Unable to find stamps --- suspect bad data quality."); + gDiffDataQuality(data, error, GDIFF_FILES_ALL); + return true; + case PM_ERR_SMALL_AREA: + psErrorStackPrint(stderr, "Insufficient area for PSF matching"); + psWarning("Insufficient area for PSF matching --- suspect bad data quality."); + gDiffDataQuality(data, error, GDIFF_FILES_ALL); + return true; + case PM_ERR_DATA: + // XXX EAM : why was this block dropped? + // XX psErrorStackPrint(stderr, "Unable to solve for the kernel to match images"); + // XX psWarning("Failed to find PSF match kernel --- suspect bad data quality."); + // XX gDiffDataQuality(data, error, GDIFF_FILES_ALL); + // XX return true; + default: + psError(GDIFF_ERR_DATA, false, "Unable to match images."); + return false; + } + } + + // Need to be careful with the normalisation + // We will normalise everything to the normalisation of the *input* image + { + // Since the entries are MULTI, we have to retrieve them differently + psMetadataIterator *regIter = psMetadataIteratorAlloc(inConv->analysis, PS_LIST_HEAD, "^" PM_SUBTRACTION_ANALYSIS_REGION "$"); + psMetadataIterator *modeIter = psMetadataIteratorAlloc(inConv->analysis, PS_LIST_HEAD, "^" PM_SUBTRACTION_ANALYSIS_MODE "$"); + psMetadataIterator *normIter = psMetadataIteratorAlloc(inConv->analysis, PS_LIST_HEAD, "^" PM_SUBTRACTION_ANALYSIS_NORM "$"); + 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 + psMetadataItem *modeItem = psMetadataGetAndIncrement(modeIter); // Item with mode + psAssert(modeItem && modeItem->type == PS_DATA_S32, "Expect subtraction mode"); + pmSubtractionMode mode = modeItem->data.S32; // Subtraction mode + psMetadataItem *normItem = psMetadataGetAndIncrement(normIter); // Item with normalisation + psAssert(normItem && normItem->type == PS_DATA_F32 && isfinite(normItem->data.F32), "Expect normalisation"); + float norm = normItem->data.F32; // Normalisation + + switch (mode) { + case PM_SUBTRACTION_MODE_1: // Convolved the input to match template + case PM_SUBTRACTION_MODE_DUAL: // Convolved both; template should have flux conserved + psLogMsg("gDiff", PS_LOG_INFO, "Correcting image for normalisation of %f\n", norm); + normaliseRegion(inConv->image, region, 1.0 / norm); + normaliseRegion(refConv->image, region, 1.0 / norm); + normaliseRegion(inConv->variance, region, 1.0 / PS_SQR(norm)); + normaliseRegion(refConv->variance, region, 1.0 / PS_SQR(norm)); + break; + case PM_SUBTRACTION_MODE_2: // Convolved the template to match input + // We're already happy! + psLogMsg("gDiff", PS_LOG_INFO, "Image normalisation is correct\n"); + break; + default: + psAbort("Invalid subtraction mode: %x", mode); + } + } + psFree(regIter); + psFree(modeIter); + psFree(normIter); + } + + pmConceptsCopyFPA(inConv->parent->parent->parent, inRO->parent->parent->parent, true, true); + pmConceptsCopyFPA(refConv->parent->parent->parent, refRO->parent->parent->parent, true, true); + + if (inConv->covariance) { + psKernel *truncated = psImageCovarianceTruncate(inConv->covariance, COVAR_FRAC); + psFree(inConv->covariance); + inConv->covariance = truncated; + } + if (refConv->covariance) { + psKernel *truncated = psImageCovarianceTruncate(refConv->covariance, COVAR_FRAC); + psFree(refConv->covariance); + refConv->covariance = truncated; + } + + if (inConv->variance) { + psImageCovarianceTransfer(inConv->variance, inConv->covariance); + } + if (refConv->variance) { + psImageCovarianceTransfer(refConv->variance, refConv->covariance); + } + + return true; +} + +/***** internal functions *****/ + +// Normalise a region on an image +static void normaliseRegion(psImage *image, // Image to normalise + const psRegion *region, // Region of image to normalise + float norm // Normalisation + ) +{ + if (!image) { + return; + } + psAssert(region, "Expect region"); + psImage *subImage = psImageSubset(image, *region); // Sub-image + psBinaryOp(subImage, subImage, "*", psScalarAlloc(norm, PS_TYPE_F32)); + psFree(subImage); + return; +} + +// Measure the PSF for an image +static float subImagePSF(gDiffData *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(GDIFF_ERR_CONFIG, false, "Unable to copy FPA for photometry"); + return NAN; + } + + pmFPAview *view = gDiffViewReadout(); // 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, "PSPHOT.INPUT", sources)) { + psErrorStackPrint(stderr, "Unable to determine PSF"); + psWarning("Unable to determine PSF."); + psFree(view); + return NAN; + } + psFree(view); + + psMetadata *header = psMetadataLookupMetadata(NULL, photRO->analysis, "PSPHOT.HEADER"); + psAssert(header, "Require header."); + float fwhm = psMetadataLookupF32(NULL, header, "FWHM_MAJ"); + + if (psMetadataLookup(photRO->analysis, "PSPHOT.DETECTIONS")) { + psMetadataRemoveKey(photRO->analysis, "PSPHOT.DETECTIONS"); + } + if (psMetadataLookup(photRO->parent->parent->analysis, "PSPHOT.PSF")) { + psMetadataRemoveKey(photRO->parent->parent->analysis, "PSPHOT.PSF"); + } + + return fwhm; +} + +// Scale the kernel parameters according to the PSFs +static bool subScaleKernel(gDiffData *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 = gDiffViewReadout(); // View to readout + pmReadout *kernelRO = pmFPAfileThisReadout(config->files, view, "GDIFF.OUTPUT.KERNELS"); // RO with kernel + if (kernelRO) { + psFree(view); + return true; + } + + // Look up recipe values + psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, GDIFF_RECIPE); // Recipe for ppSim + psAssert(recipe, "We checked this earlier, so it should be here."); + + // Input images + pmReadout *inRO = pmFPAfileThisReadout(config->files, view, "GDIFF.INPUT"); // Input readout + pmReadout *refRO = pmFPAfileThisReadout(config->files, view, "GDIFF.REF"); // Reference readout + + // Input sources + pmReadout *inSourceRO = pmFPAfileThisReadout(config->files, view, "GDIFF.INPUT.SOURCES"); + pmReadout *refSourceRO = pmFPAfileThisReadout(config->files, view, "GDIFF.REF.SOURCES"); + if (!inSourceRO || !refSourceRO) { + psError(GDIFF_ERR_DATA, false, "Unable to scale kernel, since no sources were provided."); + psFree(view); + return false; + } + + 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) { + psError(GDIFF_ERR_DATA, false, "Unable to set FWHM or scale kernel, since no sources were provided."); + return false; + } + + psArray *inSources = inDetections->allSources; + psAssert (inSources, "missing inSources?"); + + psArray *refSources = refDetections->allSources; + psAssert (refSources, "missing refSources?"); + + float inFWHM = psMetadataLookupF32(NULL, inRO->parent->parent->concepts, "CHIP.SEEING"); // FWHM for input + if (!isfinite(inFWHM) || inFWHM == 0.0) { + inFWHM = subImagePSF(data, inRO, inSources); + } + float refFWHM = psMetadataLookupF32(NULL, refRO->parent->parent->concepts, "CHIP.SEEING"); // FWHM for ref + if (!isfinite(refFWHM) || refFWHM == 0.0) { + refFWHM = subImagePSF(data, refRO, refSources); + } + psLogMsg("gDiff", PS_LOG_INFO, "Input FWHM: %f\nReference FWHM: %f\n", inFWHM, refFWHM); + if (!isfinite(inFWHM) || !isfinite(refFWHM)) { + psErrorStackPrint(stderr, "Cannot determine FHWM for images, giving up."); + int error = psErrorCodeLast(); // Error code + gDiffDataQuality(data, error, GDIFF_FILES_ALL); + return true; + } + + // we need to register the FWHM values for use downstream + pmSubtractionSetFWHMs(inFWHM, refFWHM); + + // is auto-scaling needed? + bool scale = psMetadataLookupBool(NULL, recipe, "SCALE"); + 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 + + // XXX if SCALE is FALSE, do we just skip this function? are the FHWMs used if not auto-scaled? + if (scale && (!isfinite(scaleRef) || !isfinite(scaleMin) || !isfinite(scaleMax))) { + psError(GDIFF_ERR_ARGUMENTS, false, + "auto-scale selected but scale parameters (SCALE.REF=%f, SCALE.MIN=%f, SCALE.MAX=%f) not set in recipe.", + scaleRef, scaleMin, scaleMax); + return false; + } + + pmSubtractionParamScaleOptions(scale, scaleRef, scaleMin, scaleMax); + + // if (!pmSubtractionParamsScale(kernelSize, stampSize, kernelWidths)) { + // psError(GDIFF_ERR_DATA, false, "Unable to scale parameters."); + // return false; + // } + + return true; +} + +pmSubtractionMode subModeFromString (char *string) { + + if (!strcasecmp(string, "AUTO")) return PM_SUBTRACTION_MODE_UNSURE; + if (!strcasecmp(string, "DUAL")) return PM_SUBTRACTION_MODE_DUAL; + if (!strcasecmp(string, "SINGLE_AUTO")) return PM_SUBTRACTION_MODE_SINGLE_AUTO; + if (!strcasecmp(string, "SINGLE1")) return PM_SUBTRACTION_MODE_1; + if (!strcasecmp(string, "SINGLE2")) return PM_SUBTRACTION_MODE_2; + if (!strcasecmp(string, "1")) return PM_SUBTRACTION_MODE_1; + if (!strcasecmp(string, "2")) return PM_SUBTRACTION_MODE_2; + + return PM_SUBTRACTION_MODE_UNSURE; +} + Index: /branches/eam_branches/gDiff.20200202/src/gDiffReadoutForcedPhot.c =================================================================== --- /branches/eam_branches/gDiff.20200202/src/gDiffReadoutForcedPhot.c (revision 41413) +++ /branches/eam_branches/gDiff.20200202/src/gDiffReadoutForcedPhot.c (revision 41413) @@ -0,0 +1,93 @@ +/** @file gDiffReadoutForcedPhot.c + * @brief: run forced photometry on the specified image (i.e., the input images) + * @author Eugene Magnier @ IfA + * @version $Revision: 1.12 $ + * @date $Date: 2020-01-25 $ + * Copyright 2020 Institute for Astronomy, University of Hawaii + */ + +#include "gDiff.h" + +// XXX EAM : does this function use the correct PSF for the given image?? +// XXX EAM : is this function applied to the convolved input images or the unconvolved ones? + +// run forced photometry on the image in 'targetName' at the positions of the sources in 'sourceName' +// this function is only called if (data->forcedPhot1 || data->forcedPhot2) (gDiffLoop.c:186,194) +bool gDiffReadoutForcedPhot(const char *outputName, const char *targetName, const char *sourceName, gDiffData *data) +{ + bool foundDetections = false; + bool mdok = false; + + psAssert(data, "Require processing data"); + pmConfig *config = data->config; // Configuration + psAssert(config, "Require configuration"); + + // select the view of interest + pmFPAview *view = gDiffViewReadout(); // View to readout + + pmFPAfile *sourceFile = psMetadataLookupPtr(&mdok, config->files, sourceName); // positive image sources + pmFPAfile *targetFile = psMetadataLookupPtr(&mdok, config->files, targetName); // diff image pixels + pmFPAfile *outputFile = psMetadataLookupPtr(&mdok, config->files, outputName); // result file + + psAssert (sourceFile, "failed to construct sourceName %s", sourceName); + psAssert (targetFile, "failed to construct targetName %s", targetName); + psAssert (outputFile, "failed to construct outputName %s", outputName); + + // copy the image data to PSPHOT.INPUT so that psphotReadoutForcedKnownSources does + // not affect the prior results on targetName (ie, diff detections) + pmFPAfile *photFile = psMetadataLookupPtr(&mdok, config->files, "PSPHOT.INPUT"); // Photometry file + if (!pmFPACopy(photFile->fpa, targetFile->fpa)) { + psError(GDIFF_ERR_CONFIG, false, "Unable to copy FPA for photometry"); + psFree(view); + return false; + } + + // select the reference position sources from sourceName + pmReadout *sourceRO = pmFPAviewThisReadout(view, sourceFile->fpa); + psAssert(sourceRO, "programming error: source readout not defined"); + + pmDetections *sourceDet = psMetadataLookupPtr(NULL, sourceRO->analysis, "PSPHOT.DETECTIONS"); // Sources + if (!sourceDet) { + // XXX remove the pixels from photFile? + // XXX other cleanup operations? + psFree(view); + return true; + } + psArray *sources = sourceDet->allSources; + psAssert (sources, "missing sources?"); + + // grab the PSF information from the pmFPAfile PSPHOT.PSF.LOAD + pmFPAfile *psfFile = psMetadataLookupPtr(&mdok, config->files, "PSPHOT.PSF.LOAD"); // PSF file + gDiffCopyPSF (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 (!psphotReadoutForcedKnownSources(config, view, "PSPHOT.INPUT", sources)) { + // This is likely a data quality issue + // XXX Split into multiple cases using error codes? + psErrorStackPrint(stderr, "Unable to perform photometry on image"); + psWarning("Unable to perform photometry on image --- suspect bad data quality."); + gDiffDataQuality(data, psErrorCodeLast(), GDIFF_FILES_PHOT_SUB | GDIFF_FILES_PHOT_INV); + } + // save the outputs on GDIFF.POS1.SOURCES + if (!psphotCopyResults (&foundDetections, outputFile, photFile, view)) { + psError(GDIFF_ERR_PROG, false, "Unable to copy psphot outputs"); + psFree(view); + return false; + } + + float newTime = psTimerClear("GDIFF_PHOT"); // Time for photometry + float oldTime = psMetadataLookupF32(&mdok, data->stats, "TIME_PHOT"); // Previous time for photometry + float elapsed = isfinite(oldTime) ? oldTime + newTime : newTime; + psMetadataAddF32(data->stats, PS_LIST_TAIL, "TIME_PHOT", PS_META_REPLACE, "Time to do photometry", elapsed); + + psFree(view); + return true; +} Index: /branches/eam_branches/gDiff.20200202/src/gDiffReadoutInverse.c =================================================================== --- /branches/eam_branches/gDiff.20200202/src/gDiffReadoutInverse.c (revision 41413) +++ /branches/eam_branches/gDiff.20200202/src/gDiffReadoutInverse.c (revision 41413) @@ -0,0 +1,60 @@ +/** @file gDiffReadoutInverse.c + * @brief: generate the inverse image (-1 * output) and output structures + * @author Eugene Magnier @ IfA + * @version $Revision: 1.12 $ + * @date $Date: 2020-01-25 $ + * Copyright 2020 Institute for Astronomy, University of Hawaii + */ + +#include "gDiff.h" + +bool gDiffReadoutInverse(pmConfig *config) +{ + psAssert(config, "Require configuration"); + // Check configuration for convolve option + psMetadata *recipe = psMetadataLookupPtr(NULL, config->recipes, GDIFF_RECIPE); + bool noConvolve = psMetadataLookupBool(NULL, recipe, "NOCONVOLVE"); // Do not use convolved images. + + pmFPAview *view = gDiffViewReadout(); // View to readout + pmReadout *outRO = pmFPAfileThisReadout(config->files, view, "GDIFF.OUTPUT"); + pmReadout *invRO = pmFPAfileThisReadout(config->files, view, "GDIFF.INVERSE"); + + invRO->image = (psImage*)psBinaryOp(invRO->image, outRO->image, "*", psScalarAlloc(-1.0, PS_TYPE_F32)); + invRO->mask = psMemIncrRefCounter(outRO->mask); + invRO->variance = psMemIncrRefCounter(outRO->variance); + invRO->covariance = psMemIncrRefCounter(outRO->covariance); + invRO->analysis = psMetadataCopy(invRO->analysis, outRO->analysis); + + invRO->data_exists = invRO->parent->data_exists = invRO->parent->parent->data_exists = true; + + // Get concepts from reference + pmFPAfile *refFile; + if (noConvolve) { + refFile = psMetadataLookupPtr(NULL, config->files, "GDIFF.REF"); // File with concepts + } + else { + refFile = psMetadataLookupPtr(NULL, config->files, "GDIFF.REF.CONV"); // File with concepts + } + pmFPA *invFPA = invRO->parent->parent->parent; // Inverse FPA + pmConceptsCopyFPA(invFPA, refFile->fpa, true, true); + + // Get astrometry from (forward) subtraction + pmChip *outChip = outRO->parent->parent; // Output chip + pmFPA *outFPA = outChip->parent; // Output FPA + pmChip *invChip = invRO->parent->parent; // Inverse chip + pmHDU *invHDU = invFPA->hdu; // Inverse HDU + if (!pmAstromWriteWCS(invHDU->header, outFPA, outChip, WCS_TOLERANCE)) { + psError(psErrorCodeLast(), false, "Unable to write WCS astrometry to GDIFF.INVERSE."); + psFree(view); + return false; + } + // Read from newly written astrometry so that it exists in the "inverse" FPA (for sources) + if (!pmAstromReadWCS(invFPA, invChip, invHDU->header, 1.0)) { + psError(psErrorCodeLast(), false, "Unable to read WCS astrometry."); + psFree(view); + return false; + } + + psFree(view); + return true; +} Index: /branches/eam_branches/gDiff.20200202/src/gDiffReadoutJpeg.c =================================================================== --- /branches/eam_branches/gDiff.20200202/src/gDiffReadoutJpeg.c (revision 41413) +++ /branches/eam_branches/gDiff.20200202/src/gDiffReadoutJpeg.c (revision 41413) @@ -0,0 +1,145 @@ +/** @file gDiffReadoutJpeg.c + * @brief generate JPEG images + * @author Eugene Magnier @ IfA + * @version $Revision: 1.12 $ + * @date $Date: 2020-01-25 $ + * Copyright 2020 Institute for Astronomy, University of Hawaii + */ + +#include "gDiff.h" + +bool gDiffReadoutJpeg(pmConfig *config) +{ + psAssert(config, "Require configuration"); + + pmFPAview *view = gDiffViewReadout(); // View to readout + pmReadout *outRO = pmFPAfileThisReadout(config->files, view, "GDIFF.OUTPUT"); // Readout to jpeg + + // Look up recipe values + psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, GDIFF_RECIPE); // Recipe for ppSim + psAssert(recipe, "We checked this earlier, so it should be here."); + + // 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, "GDIFF.OUTPUT.JPEG1"); // Rebinned cell once + pmCell *cell2 = pmFPAfileThisCell(config->files, view, "GDIFF.OUTPUT.JPEG2"); // Rebinned cell twice + psFree(view); + + pmReadout *ro1 = pmReadoutAlloc(cell1), *ro2 = pmReadoutAlloc(cell2); // Binned readouts + if (!pmReadoutRebin(ro1, outRO, maskBad, bin1, bin1)) { + psError(GDIFF_ERR_DATA, false, "Unable to bin output (1st binning)"); + psFree(ro1); + psFree(ro2); + return false; + } + if (!pmReadoutRebin(ro2, ro1, 0, bin2, bin2)) { + psError(GDIFF_ERR_DATA, false, "Unable to bin output (2nd binning)"); + psFree(ro1); + psFree(ro2); + return false; + } + psFree(ro1); + psFree(ro2); + + return true; +} + +// XXX EAM : this function is mis-named: it generates a FITS image +bool gDiffResidualSampleJpeg(pmConfig *config) +{ + // XXX EAM : the option to run this should be in the recipe or command-line options + return true; + + pmFPAview *view = gDiffViewReadout(); // View to readout + + // we save sample difference stamps on the convolved image readout->analysis metadata + pmReadout *inConv = pmFPAfileThisReadout(config->files, view, "GDIFF.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, "GDIFF.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); + } + + psFree(view); + return true; +} Index: /branches/eam_branches/gDiff.20200202/src/gDiffReadoutPhotometry.c =================================================================== --- /branches/eam_branches/gDiff.20200202/src/gDiffReadoutPhotometry.c (revision 41413) +++ /branches/eam_branches/gDiff.20200202/src/gDiffReadoutPhotometry.c (revision 41413) @@ -0,0 +1,133 @@ +/** @file gDiffReadoutPhotometry.c + * @brief: runs psphotReadoutMinimal on the specified image + * @author Eugene Magnier @ IfA + * @version $Revision: 1.12 $ + * @date $Date: 2020-01-25 $ + * Copyright 2020 Institute for Astronomy, University of Hawaii + */ + +#include "gDiff.h" + +/* This function runs psphotReadoutMinimal on the image specified by 'name' (may be the subtracted image or the inverse) + * The psphot call uses the supplied PSF which is generated in gDiffMakePSF using the positive image (minuend) + * which is the convolved image if convolution was performed. + */ + +bool gDiffReadoutPhotometry(const char *name, gDiffData *data) +{ + bool mdok = false; + + psAssert(data, "Require processing data"); + pmConfig *config = data->config; // Configuration + psAssert(config, "Require configuration"); + + if (!data->photometry) { + return true; + } + + // select the view of interest + pmFPAview *view = gDiffViewReadout(); // View to readout + + // psphotReadoutMinimal performs the photometry analysis on PSPHOT.INPUT; we need to move + // 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 + if (!pmFPACopy(photFile->fpa, inFile->fpa)) { + psError(GDIFF_ERR_CONFIG, false, "Unable to copy FPA for photometry"); + psFree(view); + return false; + } + + // psphot may need the PHU header + if (!photFile->fpa->hdu) { + photFile->fpa->hdu = psMemIncrRefCounter (inFile->fpa->hdu); + } + + // 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.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"); + } + + // grab the PSF information from the pmFPAfile PSPHOT.PSF.LOAD + pmFPAfile *psfFile = psMetadataLookupPtr(&mdok, config->files, "PSPHOT.PSF.LOAD"); // PSF file + + gDiffCopyPSF (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, "PSPHOT.INPUT")) { + // This is likely a data quality issue + // XXX Split into multiple cases using error codes? + psErrorStackPrint(stderr, "Unable to perform photometry on image"); + psWarning("Unable to perform photometry on image --- suspect bad data quality."); + gDiffDataQuality(data, psErrorCodeLast(), GDIFF_FILES_PHOT_SUB | GDIFF_FILES_PHOT_INV); + } + + // If no sources were found, there's no error, but we want to trigger 'bad quality' + psWarning("no sources found: why is this being set to bad quality??"); + pmDetections *detections = psMetadataLookupPtr(NULL, photRO->analysis, "PSPHOT.DETECTIONS"); // Sources + if (!detections) { + gDiffDataQuality(data, PSPHOT_ERR_DATA, GDIFF_FILES_PHOT_SUB | GDIFF_FILES_PHOT_INV); + } + psArray *sources = detections->allSources; + psAssert (sources, "missing sources?"); + + // a likely source of false positives is poor subtractions. this results in + // detections in the wings (or cores) of bright(er) stars found in both images. + // flag detections based on their distance from the bright(er) input sources. + bool subInverse = !strcasecmp(name, "GDIFF.INVERSE"); + gDiffFlagNeighbors (config, view, sources, subInverse); + + if (data->stats) { + bool mdok; + int numSources = psMetadataLookupS32(&mdok, data->stats, "NUM_SOURCES"); // Number of sources + numSources += sources ? sources->n : 0; + psMetadataAddS32(data->stats, PS_LIST_TAIL, "NUM_SOURCES", PS_META_REPLACE, + "Total number of sources detected", numSources); + float newTime = psTimerClear("GDIFF_PHOT"); // Time for photometry + float oldTime = psMetadataLookupF32(&mdok, data->stats, "TIME_PHOT"); // Previous time for photometry + psMetadataAddF32(data->stats, PS_LIST_TAIL, "TIME_PHOT", PS_META_REPLACE, "Time to do photometry", + isfinite(oldTime) ? oldTime + newTime : newTime); + } + + if (!data->quality) { + if (!psMetadataCopySingle(inRO->analysis, photRO->analysis, "PSPHOT.DETECTIONS")) { + psError(GDIFF_ERR_PROG, false, "Unable to copy PSPHOT.DETECTIONS"); + psFree(view); + return false; + } + if (!psMetadataCopySingle(inRO->analysis, photRO->analysis, "PSPHOT.HEADER")) { + psError(GDIFF_ERR_PROG, false, "Unable to copy PSPHOT.HEADER"); + psFree(view); + return false; + } + if (!psMetadataCopySingle(inRO->analysis, photRO->analysis, PM_DETEFF_ANALYSIS)) { + psError(GDIFF_ERR_PROG, false, "Unable to copy Detection Efficiency"); + psFree(view); + 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); + } + } + + psFree(view); + return true; +} Index: /branches/eam_branches/gDiff.20200202/src/gDiffReadoutStats.c =================================================================== --- /branches/eam_branches/gDiff.20200202/src/gDiffReadoutStats.c (revision 41413) +++ /branches/eam_branches/gDiff.20200202/src/gDiffReadoutStats.c (revision 41413) @@ -0,0 +1,50 @@ +/** @file gDiffReadoutStats.c + * @brief: calculate and save stats on the subtraction image + * @author Eugene Magnier @ IfA + * @version $Revision: 1.12 $ + * @date $Date: 2020-01-25 $ + * Copyright 2020 Institute for Astronomy, University of Hawaii + */ + +#include "gDiff.h" + +bool gDiffReadoutStats(gDiffData *data) +{ + psAssert(data, "Require processing data"); + pmConfig *config = data->config; // Configuration + psAssert(config, "Require configuration"); + + if (!data->statsFile) { + // Nothing to do + return true; + } + + pmFPAfile *output = psMetadataLookupPtr(NULL, config->files, "GDIFF.OUTPUT"); // Output file + if (!output) { + psError(GDIFF_ERR_PROG, true, "Unable to find file GDIFF.OUTPUT.\n"); + return false; + } + psImageMaskType maskValue = pmConfigMaskGet("MASK.VALUE", config); + pmFPAview *view = gDiffViewReadout(); // View to readout + ppStatsFPA(data->stats, output->fpa, view, maskValue, config); + + pmReadout *outRO = pmFPAviewThisReadout(view, output->fpa); // Readout of interest + psFree(view); + + // Statistics on the matching + psMetadataCopySingle(data->stats, outRO->analysis, PM_SUBTRACTION_ANALYSIS_MODE); + psMetadataCopySingle(data->stats, outRO->analysis, PM_SUBTRACTION_ANALYSIS_STAMPS); + psMetadataCopySingle(data->stats, outRO->analysis, PM_SUBTRACTION_ANALYSIS_DEV_MEAN); + psMetadataCopySingle(data->stats, outRO->analysis, PM_SUBTRACTION_ANALYSIS_DEV_RMS); + psMetadataCopySingle(data->stats, outRO->analysis, PM_SUBTRACTION_ANALYSIS_NORM); + psMetadataCopySingle(data->stats, outRO->analysis, PM_SUBTRACTION_ANALYSIS_BGDIFF); + psMetadataCopySingle(data->stats, outRO->analysis, PM_SUBTRACTION_ANALYSIS_MX); + psMetadataCopySingle(data->stats, outRO->analysis, PM_SUBTRACTION_ANALYSIS_MY); + psMetadataCopySingle(data->stats, outRO->analysis, PM_SUBTRACTION_ANALYSIS_MXX); + psMetadataCopySingle(data->stats, outRO->analysis, PM_SUBTRACTION_ANALYSIS_MXY); + psMetadataCopySingle(data->stats, outRO->analysis, PM_SUBTRACTION_ANALYSIS_MYY); + psMetadataCopySingle(data->stats, outRO->analysis, PM_SUBTRACTION_ANALYSIS_DECONV_MAX); + psMetadataCopySingle(data->stats, outRO->analysis, PM_SUBTRACTION_ANALYSIS_CONVOL_MAX); + + return true; +} Index: /branches/eam_branches/gDiff.20200202/src/gDiffReadoutSubtract.c =================================================================== --- /branches/eam_branches/gDiff.20200202/src/gDiffReadoutSubtract.c (revision 41413) +++ /branches/eam_branches/gDiff.20200202/src/gDiffReadoutSubtract.c (revision 41413) @@ -0,0 +1,149 @@ +/** @file gDiffReadoutSubtract.c + * @brief: image subtraction based on matched Gaussians profiles + * @author Eugene Magnier @ IfA + * @version $Revision: 1.12 $ + * @date $Date: 2020-01-25 $ + * Copyright 2020 Institute for Astronomy, University of Hawaii + */ + +#include "gDiff.h" + +bool gDiffReadoutSubtract(pmConfig *config) +{ + bool mdok = false; // Status of MD lookup + + psAssert(config, "Require configuration"); + + // Look up recipe values + psMetadata *recipe = psMetadataLookupMetadata(&mdok, config->recipes, GDIFF_RECIPE); // Recipe for gDiff + psAssert(recipe, "We checked this earlier, so it should be here."); + + bool reverse = psMetadataLookupBool(&mdok, config->arguments, "REVERSE"); // Reverse sense of subtraction? + bool noConvolve = psMetadataLookupBool(&mdok, recipe, "NOCONVOLVE"); // Do not use convolved images. + bool addPair = psMetadataLookupBool(&mdok, recipe, "ADD.NOT.SUBTRACT"); // add instead of subtracting + bool doApplyMaskNaN = psMetadataLookupBool(&mdok, recipe, "APPLY.PIXELNAN"); // NaN the pixels underneath masks + + pmFPAview *view = gDiffViewReadout(); // View to readout + + // XXX EAM : noConvolve should generate a copy of the input images in the .CONV files, so we do not need to + // make these ad-hoc tests + + // Subtraction is: minuend - subtrahend + pmReadout *minuend = NULL; // Positive image + pmReadout *subtrahend = NULL; // Negative image + if (!noConvolve) { + printf("Using Convolved images because NOCONVOLVE is FALSE\n"); + if (reverse) { + minuend = pmFPAfileThisReadout(config->files, view, "GDIFF.REF.CONV"); + subtrahend = pmFPAfileThisReadout(config->files, view, "GDIFF.INPUT.CONV"); + } else { + minuend = pmFPAfileThisReadout(config->files, view, "GDIFF.INPUT.CONV"); + subtrahend = pmFPAfileThisReadout(config->files, view, "GDIFF.REF.CONV"); + } + } else { + psWarning("Not using Convolved images because NOCONVOLVE is TRUE\n"); + if (reverse) { + minuend = pmFPAfileThisReadout(config->files, view, "GDIFF.REF"); + subtrahend = pmFPAfileThisReadout(config->files, view, "GDIFF.INPUT"); + } else { + minuend = pmFPAfileThisReadout(config->files, view, "GDIFF.INPUT"); + subtrahend = pmFPAfileThisReadout(config->files, view, "GDIFF.REF"); + } + } + + // Do the actual subtraction + pmReadout *outRO = pmFPAfileThisReadout(config->files, view, "GDIFF.OUTPUT"); + + if (addPair) { + outRO->image = (psImage*)psBinaryOp(outRO->image, minuend->image, "+", subtrahend->image); + } else { + outRO->image = (psImage*)psBinaryOp(outRO->image, minuend->image, "-", subtrahend->image); + } + outRO->mask = (psImage*)psBinaryOp(outRO->mask, minuend->mask, "|", subtrahend->mask); + outRO->variance = (psImage*)psBinaryOp(outRO->variance, minuend->variance, "+", subtrahend->variance); + + // NAN the masked pixels in the diff image (pixels masked in A are not yet NAN'ed in B) + psImageMaskType maskVal = pmConfigMaskGet("MASK.VALUE", config) | pmConfigMaskGet("BLANK", config); // Bits to mask in inputs + if(doApplyMaskNaN) { + for (int iy = 0; iy < outRO->image->numRows; iy++) { + for (int ix = 0; ix < outRO->image->numCols; ix++) { + if ((outRO->mask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] & maskVal) == 0) continue; + outRO->image->data.F32[iy][ix] = NAN; + } + } + } + + // Measure the variance scales + psStats *varStats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN); // Statistics for variance images + psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS); // Random number generator + psImageBackground(varStats, NULL, minuend->variance, minuend->mask, maskVal, rng); + float minuendVar = varStats->robustMedian; // Mean variance for minuend + psImageBackground(varStats, NULL, subtrahend->variance, subtrahend->mask, maskVal, rng); + float subtrahendVar = varStats->robustMedian; // Mean variance for subtrahend + psFree(varStats); + psFree(rng); + + // XXX EAM : I suspect that the weighted averaging is giving the wrong answer for + // single-direction PSF matching (and maybe DUAL as well). I am testing the difference by + // generating A+B images instead of A-B images and then checking the significance of the + // sources in the image + + // Combine the covariances + // These are weighted by the appropriate mean variance. This is probably not perfectly correct, but it + // does seem to reproduce the correct magnitude limit in psphot. + psArray *covars = psArrayAlloc(2); // Covariance pseudo-matrices + psVector *covarWeights = psVectorAlloc(2, PS_TYPE_F32); // Weights for covariances + covars->data[0] = psMemIncrRefCounter(minuend->covariance); + covars->data[1] = psMemIncrRefCounter(subtrahend->covariance); + covarWeights->data.F32[0] = minuendVar; + covarWeights->data.F32[1] = subtrahendVar; +# if (0) + outRO->covariance = psImageCovarianceAverageWeighted(covars, covarWeights); +# else + outRO->covariance = psImageCovarianceAverage(covars); +# endif + psFree(covars); + psFree(covarWeights); + + psImageCovarianceTransfer(outRO->variance, outRO->covariance); + + outRO->data_exists = true; + outRO->parent->data_exists = true; + outRO->parent->parent->data_exists = true; + + pmSubtractionVisualShowSubtraction(minuend->image, subtrahend->image, outRO->image); + + // Copy concepts from the input to the output + pmFPAfile *inFile = psMetadataLookupPtr(&mdok, config->files, "GDIFF.INPUT"); // Input file + pmFPA *inFPA = inFile->fpa; // Input FPA + pmFPAfile *outFile = psMetadataLookupPtr(&mdok, config->files, "GDIFF.OUTPUT"); // Output file + pmFPA *outFPA = outFile->fpa; // Output FPA + if (!pmConceptsCopyFPA(outFPA, inFPA, true, true)) { + psError(GDIFF_ERR_CONFIG, false, "Unable to copy concepts from input to output."); + psFree(outRO); + psFree(view); + return false; + } + + // Copy astrometry over + // It should find its way into the output images and photometry + pmHDU *inHDU = inFPA->hdu; // Input HDU + pmHDU *outHDU = outFPA->hdu; // Output HDU + pmChip *outChip = pmFPAfileThisChip(config->files, view, "GDIFF.OUTPUT"); // Output chip + psFree(view); + + if (!outHDU || !inHDU) { + psError(GDIFF_ERR_PROG, true, "Unable to find HDU at FPA level to copy astrometry."); + return false; + } + if (!pmAstromReadWCS(outFPA, outChip, inHDU->header, 1.0)) { + psError(psErrorCodeLast(), false, "Unable to read WCS astrometry from input FPA."); + return false; + } + if (!pmAstromWriteWCS(outHDU->header, outFPA, outChip, WCS_TOLERANCE)) { + psError(psErrorCodeLast(), false, "Unable to write WCS astrometry to output FPA."); + return false; + } + + return true; +} Index: /branches/eam_branches/gDiff.20200202/src/gDiffSetMasks.c =================================================================== --- /branches/eam_branches/gDiff.20200202/src/gDiffSetMasks.c (revision 41413) +++ /branches/eam_branches/gDiff.20200202/src/gDiffSetMasks.c (revision 41413) @@ -0,0 +1,229 @@ +/** @file gDiffSetMasks.c + * @brief: image subtraction based on matched Gaussians profiles + * @author Eugene Magnier @ IfA + * @version $Revision: 1.12 $ + * @date $Date: 2020-01-25 $ + * Copyright 2020 Institute for Astronomy, University of Hawaii + */ + +#include "gDiff.h" + +// Structure to hold the properties of a mask value +typedef struct { + char *badMaskName; // name for "bad" (i.e., mask me please) pixels + char *fallbackName; // Fallback name in case a bad mask name is not defined + psImageMaskType defaultMaskValue; // Default value in case a bad mask name and its fallback are not defined + bool isBad; // include this value as part of the MASK.VALUE entry (generically bad) +} pmConfigMaskInfo; + +static pmConfigMaskInfo skycellmasks[] = { + // Features of the detector + { "DETECTOR", NULL, 0x01, false }, // Something is wrong with the detector + { "FLAT", "DETECTOR", 0x01, false }, // Pixel doesn't flat-field properly + { "DARK", "DETECTOR", 0x01, false }, // Pixel doesn't dark-subtract properly + { "BLANK", "DETECTOR", 0x01, true }, // Pixel doesn't contain valid data + { "CTE", "DETECTOR", 0x01, false }, // Pixel has poor CTE + { "BURNTOOL", NULL, 0x04, false }, // Pixel has been touched by burntool + // Invalid signal ranges + { "SAT", NULL, 0x02, true }, // Pixel is saturated or non-linear + { "LOW", "SAT", 0x02, true }, // Pixel is low + { "SUSPECT", NULL, 0x04, false }, // Pixel is suspected of being bad + // Non-astronomical structures + { "CR", NULL, 0x08, true }, // Pixel contains a cosmic ray + { "SPIKE", NULL, 0x08, false }, // Pixel contains a diffraction spike + { "GHOST", NULL, 0x08, false }, // Pixel contains an optical ghost + { "STREAK", NULL, 0x08, false }, // Pixel contains streak data + { "STARCORE", NULL, 0x08, false }, // Pixel contains a bright star core + // Effects of convolution and interpolation + { "CONV.BAD", NULL, 0x02, true }, // Pixel is bad after convolution with a bad pixel + { "CONV.POOR", NULL, 0x04, false }, // Pixel is poor after convolution with a bad pixel +}; + +// internal function: +bool gDiffMaskSetInMetadata(psImageMaskType *outMaskValue, psImageMaskType *outMarkValue, psMetadata *source); + +bool gDiffSetMasks(pmConfig *config) +{ + psAssert(config, "Require configuration"); + + // Look up recipe values + psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, GDIFF_RECIPE); // Recipe for ppSim + psAssert(recipe, "We checked this earlier, so it should be here."); + bool doApplyMaskNaN = psMetadataLookupBool(NULL, recipe, "APPLY.PIXELNAN"); // NaN the pixels underneath masks + + psImageMaskType maskValue, markValue; // Mask values + if(doApplyMaskNaN) { + if (!pmConfigMaskSetBits(&maskValue, &markValue, config)) { + psError(GDIFF_ERR_CONFIG, false, "Unable to determine mask value."); + return false; + } + } else { + psMetadata *maskrecipe = psMetadataLookupMetadata(NULL, config->recipes, "MASKS"); // The recipe + if (!maskrecipe) { + psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to find MASKS recipe."); + return false; + } + if (!gDiffMaskSetInMetadata(&maskValue, &markValue, maskrecipe)) { + psError(GDIFF_ERR_CONFIG, false, "Unable to determine mask value."); + return false; + } + } + + // Set the mask bits needed by psphot (in psphot recipe) + psphotSetMaskRecipe(config, maskValue, markValue); + + + psImageMaskType satValue = pmConfigMaskGet("SAT", config); + psAssert(satValue, "SAT must be non-zero"); + + psImageMaskType lowValue = pmConfigMaskGet("LOW", config); + if (!lowValue) { + // Look up old name for backward compatability + lowValue = pmConfigMaskGet("BAD", config); + } + psAssert(lowValue, "LOW or BAD must be non-zero"); + + // Input images + pmFPAview *view = gDiffViewReadout(); // View to readout + pmReadout *inRO = pmFPAfileThisReadout(config->files, view, "GDIFF.INPUT"); // Input readout + pmReadout *refRO = pmFPAfileThisReadout(config->files, view, "GDIFF.REF"); // Reference readout + + psImage *input = inRO->image; // Input image + psImage *reference = refRO->image; // Reference image + PS_ASSERT_IMAGES_SIZE_EQUAL(input, reference, false); + int numCols = input->numCols, numRows = input->numRows; // Size of image + + // Generate masks if they don't exist + if (!inRO->mask) { + if (psMetadataLookupBool(NULL, recipe, "MASK.GENERATE")) { + pmReadoutSetMask(inRO, satValue, lowValue); + } 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, lowValue); + } else { + refRO->mask = psImageAlloc(numCols, numRows, PS_TYPE_IMAGE_MASK); + psImageInit(refRO->mask, 0); + } + } + + // Mask the NAN values (USE BLANK instead of SAT?) + if (!pmReadoutMaskInvalid(inRO, maskValue, satValue)) { + psError(GDIFF_ERR_DATA, false, "Unable to mask non-finite pixels in input."); + psFree(view); + return false; + } + if (!pmReadoutMaskInvalid(refRO, maskValue, satValue)) { + psError(GDIFF_ERR_DATA, false, "Unable to mask non-finite pixels in reference."); + psFree(view); + return false; + } + +#if 0 + // Interpolation over bad pixels: this takes a while! + bool mdok = false; + psString interpModeStr = psMetadataLookupStr(&mdok, recipe, "INTERPOLATION"); // Interpolation mode + psImageInterpolateMode interpMode = psImageInterpolateModeFromString(interpModeStr); // Interp + if (interpMode == PS_INTERPOLATE_NONE) { + psError(GDIFF_ERR_CONFIG, false, "Unknown interpolation mode: %s", interpModeStr); + psFree(view); + return false; + } + float poorFrac = psMetadataLookupF32(&mdok, recipe, "POOR.FRACTION"); // Fraction for "poor" + + psImageMaskType maskPoor = pmConfigMaskGet("CONV.POOR", 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(GDIFF_ERR_DATA, false, "Unable to interpolate bad pixels for input image."); + psFree(view); + return false; + } + if (!pmReadoutInterpolateBadPixels(refRO, maskVal, interpMode, poorFrac, maskPoor, maskBad)) { + psError(GDIFF_ERR_DATA, false, "Unable to interpolate bad pixels for reference image."); + psFree(view); + return false; + } + maskVal |= maskBad; +#endif + + psFree(view); + return true; +} + +// XXX EAM : this function seems to be identical to pmConfigMask.c: pmConfigMaskSetInMetadata +// why is a different function needed here? +// XXX EAM : answer: because a few of the "isBad" values are defined differently for ppSub vs other analysis +// this seems like a clear case where a recipe should have been used. + +bool gDiffMaskSetInMetadata(psImageMaskType *outMaskValue, // Value of MASK.VALUE, returned + psImageMaskType *outMarkValue, // Value of MARK.VALUE, returned + psMetadata *source // Source of mask bits + ) +{ + PS_ASSERT_METADATA_NON_NULL(source, false); + + // Ensure all the bad mask names exist, and set the value to catch all bad pixels + psImageMaskType maskValue = 0; // Value to mask to catch all the bad pixels + psImageMaskType allMasks = 0; // Value to mask to catch all masked bits (to set MARK) + + int nMasks = sizeof (skycellmasks) / sizeof (pmConfigMaskInfo); + + for (int i = 0; i < nMasks; i++) { + bool mdok; // Status of MD lookup + psImageMaskType value = psMetadataLookupImageMaskFromGeneric(&mdok, source, skycellmasks[i].badMaskName); // Value of mask + if (!mdok) { + psWarning ("problem with mask value %s\n", skycellmasks[i].badMaskName); + } + + if (!value) { + if (skycellmasks[i].fallbackName) { + value = psMetadataLookupImageMaskFromGeneric(&mdok, source, skycellmasks[i].fallbackName); + } + if (!value) { + value = skycellmasks[i].defaultMaskValue; + } + psMetadataAddImageMask(source, PS_LIST_TAIL, skycellmasks[i].badMaskName, PS_META_REPLACE, NULL, value); + } + if (skycellmasks[i].isBad) { + maskValue |= value; + } + allMasks |= value; + } + + // search for an unset bit to use for MARK: + psImageMaskType markValue = 0x00; + psImageMaskType markTrial = 0x01; + + int nBits = sizeof(psImageMaskType) * 8; + for (int i = 0; !markValue && (i < nBits); i++) { + if (allMasks & markTrial) { + markTrial <<= 1; + } else { + markValue = markTrial; + } + } + if (!markValue) { + psError (PS_ERR_UNKNOWN, true, "Unable to define the MARK bit mask: all bits taken!"); + return false; + } + + // update the list with the results + psMetadataAddImageMask(source, PS_LIST_TAIL, "MASK.VALUE", PS_META_REPLACE, NULL, maskValue); + psMetadataAddImageMask(source, PS_LIST_TAIL, "MARK.VALUE", PS_META_REPLACE, NULL, markValue); + + if (outMaskValue) { + *outMaskValue = maskValue; + } + if (outMarkValue) { + *outMarkValue = markValue; + } + + return true; +} + Index: /branches/eam_branches/gDiff.20200202/src/gDiffThreshold.c =================================================================== --- /branches/eam_branches/gDiff.20200202/src/gDiffThreshold.c (revision 41413) +++ /branches/eam_branches/gDiff.20200202/src/gDiffThreshold.c (revision 41413) @@ -0,0 +1,138 @@ +/** @file gDiffThreshold.c + * @brief Mask pixels below threshold + * @author Eugene Magnier @ IfA + * @version $Revision: 1.12 $ + * @date $Date: 2020-01-25 $ + * Copyright 2020 Institute for Astronomy, University of Hawaii + */ + +#include "gDiff.h" + +// internal function: +static bool lowThreshold(pmReadout *ro, float thresh, psImageMaskType maskIgnore, psImageMaskType maskThresh, psRegion *region, const char *description); + +bool gDiffLowThreshold(gDiffData *data) +{ + psAssert(data, "Require processing data"); + pmConfig *config = data->config; // Configuration + psAssert(config, "Require configuration"); + + // Look up recipe values + psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, GDIFF_RECIPE); // Recipe for ppSim + psAssert(recipe, "We checked this earlier, so it should be here."); + + psImageMaskType maskVal = pmConfigMaskGet("MASK.VALUE", config); // Bits to mask in inputs + psImageMaskType maskThresh = pmConfigMaskGet("LOW", config); // Bits to mask for low pixels + + float thresh = psMetadataLookupF32(NULL, recipe, "LOW.THRESHOLD"); // Threshold for masking + if (!isfinite(thresh)) { + return true; + } + + pmFPAview *view = gDiffViewReadout(); // View to readout + + // Input images + bool noConvolve = psMetadataLookupBool(NULL, recipe, "NOCONVOLVE"); + + pmReadout *in; + if (noConvolve) { + in = pmFPAfileThisReadout(config->files, view, "GDIFF.INPUT"); // Input image + } + else { + in = pmFPAfileThisReadout(config->files, view, "GDIFF.INPUT.CONV"); // Input image + } + if (!in) { + psError(GDIFF_ERR_UNKNOWN, false, "Unable to find readout."); + psFree(view); + return false; + } + + pmReadout *ref; + if (noConvolve) { + ref = pmFPAfileThisReadout(config->files, view, "GDIFF.REF"); // Reference image + } + else { + ref = pmFPAfileThisReadout(config->files, view, "GDIFF.REF.CONV"); // Reference image + } + if (!ref) { + psError(GDIFF_ERR_UNKNOWN, false, "Unable to find readout."); + psFree(view); + 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(psErrorCodeLast(), false, "Unable to threshold input image."); + psFree(view); + return false; + } + if (!lowThreshold(ref, thresh, maskVal, maskThresh, region, "reference convolved image")) { + psError(psErrorCodeLast(), false, "Unable to threshold input image."); + psFree(view); + return false; + } + } + psFree(regIter); + + psFree(view); + + return true; +} + +/*** internal functions ***/ +// Apply low threshold to image +static bool lowThreshold( + pmReadout *ro, // Readout to threshold + float thresh, // Threshold to apply + 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 + ) +{ + PM_ASSERT_READOUT_NON_NULL(ro, false); + PM_ASSERT_READOUT_IMAGE(ro, false); + PM_ASSERT_READOUT_MASK(ro, false); + PS_ASSERT_FLOAT_LARGER_THAN(thresh, 0.0, false); + + psImage *image = ro->image; // Image + psImage *mask = ro->mask; // Mask + + 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, subImage, subMask, maskIgnore, rng)) { + psError(GDIFF_ERR_DATA, false, "Unable to determine threshold."); + psFree(rng); + psFree(stats); + return false; + } + psFree(rng); + + psFree(subImage); + psFree(subMask); + + float threshold = stats->robustMedian - thresh * stats->robustStdev; // Threshold below which to clip + psFree(stats); + psLogMsg("gDiff", PS_LOG_INFO, "Masking pixels below %f in %s", threshold, description); + + 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; + } + if (image->data.F32[y][x] < threshold) { + mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] |= maskThresh; + } + } + } + + return true; +} + Index: /branches/eam_branches/gDiff.20200202/src/gDiffVarianceRescale.c =================================================================== --- /branches/eam_branches/gDiff.20200202/src/gDiffVarianceRescale.c (revision 41413) +++ /branches/eam_branches/gDiff.20200202/src/gDiffVarianceRescale.c (revision 41413) @@ -0,0 +1,61 @@ +/** @file gDiffVarianceRescale.c + * @brief: image subtraction based on matched Gaussians profiles + * @author Eugene Magnier @ IfA + * @version $Revision: 1.12 $ + * @date $Date: 2020-01-25 $ + * Copyright 2020 Institute for Astronomy, University of Hawaii + */ + +#include "gDiff.h" + +bool gDiffVarianceRescale(pmConfig *config, gDiffData *data) +{ + psAssert(config, "Require configuration"); + + bool mdok; // Status of metadata lookups + + psMetadata *recipe = psMetadataLookupPtr(NULL, config->recipes, GDIFF_RECIPE); // Recipe for gDiff + psAssert(recipe, "Need GDIFF recipe"); + + if (!psMetadataLookupBool(&mdok, recipe, "RENORM")) return true; + + int num = psMetadataLookupS32(&mdok, recipe, "RENORM.NUM"); + if (!mdok) { + psError(GDIFF_ERR_ARGUMENTS, true, "RENORM.NUM is not set in the recipe"); + return false; + } + float minValid = psMetadataLookupF32(&mdok, recipe, "RENORM.MIN"); + if (!mdok) { + psError(GDIFF_ERR_ARGUMENTS, true, "RENORM.MIN is not set in the recipe"); + return false; + } + float maxValid = psMetadataLookupF32(&mdok, recipe, "RENORM.MAX"); + if (!mdok) { + psError(GDIFF_ERR_ARGUMENTS, true, "RENORM.MAX is not set in the recipe"); + return false; + } + + psImageMaskType maskBad = pmConfigMaskGet("BLANK", config); // Bits to mask + + pmFPAview *view = gDiffViewReadout(); // View to readout + pmReadout *readout = pmFPAfileThisReadout(config->files, view, "GDIFF.OUTPUT"); // Output image + + if (!readout->variance) { + // Nothing to renormalise + psWarning("Renormalisation of the variance requested, but no variance provided."); + psFree(view); + return true; + } + + if (!pmReadoutVarianceRenormalise(readout, maskBad, num, minValid, maxValid)) { + psErrorStackPrint(stderr, "Unable to renormalise variances"); + psWarning("Unable to renormalise variances --- suspect bad data quality."); + // Allow the convolved and subtracted images to be written + gDiffDataQuality(data, GDIFF_ERR_VARIANCE, + GDIFF_FILES_INPUT | GDIFF_FILES_PHOT_SUB | GDIFF_FILES_PHOT_INV | + GDIFF_FILES_PSF | GDIFF_FILES_PHOT); + } + + psFree(view); + return true; +} Index: /branches/eam_branches/gDiff.20200202/src/gDiffVersion.c =================================================================== --- /branches/eam_branches/gDiff.20200202/src/gDiffVersion.c (revision 41412) +++ /branches/eam_branches/gDiff.20200202/src/gDiffVersion.c (revision 41413) @@ -1,4 +1,11 @@ +/** @file gDiffVersion.c + * @brief: image subtraction based on matched Gaussians profiles + * @author Eugene Magnier @ IfA + * @version $Revision: 1.12 $ + * @date $Date: 2020-01-25 $ + * Copyright 2020 Institute for Astronomy, University of Hawaii + */ + #include "gDiff.h" -#include #include "gDiffVersionDefinitions.h"