Index: /branches/eam_branches/gDiff.20200202/Makefile.am =================================================================== --- /branches/eam_branches/gDiff.20200202/Makefile.am (revision 41260) +++ /branches/eam_branches/gDiff.20200202/Makefile.am (revision 41260) @@ -0,0 +1,3 @@ +SUBDIRS = src + +CLEANFILES = *.pyc *~ core core.* Index: /branches/eam_branches/gDiff.20200202/autogen.sh =================================================================== --- /branches/eam_branches/gDiff.20200202/autogen.sh (revision 41260) +++ /branches/eam_branches/gDiff.20200202/autogen.sh (revision 41260) @@ -0,0 +1,103 @@ +#!/bin/sh +# Run this to generate all the initial makefiles, etc. + +srcdir=`dirname $0` +test -z "$srcdir" && srcdir=. + +ORIGDIR=`pwd` +cd $srcdir + +PROJECT=gDiff +TEST_TYPE=-f +# change this to be a unique filename in the top level dir +FILE=autogen.sh + +DIE=0 + +LIBTOOLIZE=libtoolize +ACLOCAL="aclocal $ACLOCAL_FLAGS" +AUTOHEADER=autoheader +AUTOMAKE=automake +AUTOCONF=autoconf + +($LIBTOOLIZE --version) < /dev/null > /dev/null 2>&1 || { + echo + echo "You must have $LIBTOOLIZE installed to compile $PROJECT." + echo "Download the appropriate package for your distribution," + echo "or get the source tarball at http://ftp.gnu.org/gnu/libtool/" + DIE=1 +} + +($ACLOCAL --version) < /dev/null > /dev/null 2>&1 || { + echo + echo "You must have $ACLOCAL installed to compile $PROJECT." + echo "Download the appropriate package for your distribution," + echo "or get the source tarball at http://ftp.gnu.org/gnu/automake/" + DIE=1 +} + +($AUTOHEADER --version) < /dev/null > /dev/null 2>&1 || { + echo + echo "You must have $AUTOHEADER installed to compile $PROJECT." + echo "Download the appropriate package for your distribution," + echo "or get the source tarball at http://ftp.gnu.org/gnu/autoconf/" + DIE=1 +} + +($AUTOMAKE --version) < /dev/null > /dev/null 2>&1 || { + echo + echo "You must have $AUTOMAKE installed to compile $PROJECT." + echo "Download the appropriate package for your distribution," + echo "or get the source tarball at http://ftp.gnu.org/gnu/automake/" + DIE=1 +} + +($AUTOCONF --version) < /dev/null > /dev/null 2>&1 || { + echo + echo "You must have $AUTOCONF installed to compile $PROJECT." + echo "Download the appropriate package for your distribution," + echo "or get the source tarball at http://ftp.gnu.org/gnu/autoconf/" + DIE=1 +} + +if test "$DIE" -eq 1; then + exit 1 +fi + +test $TEST_TYPE $FILE || { + echo "You must run this script in the top-level $PROJECT directory" + exit 1 +} + +if test -z "$*"; then + echo "I am going to run ./configure with no arguments - if you wish " + echo "to pass any to it, please specify them on the $0 command line." +fi + +$LIBTOOLIZE --copy --force || echo "$LIBTOOLIZE failed" +$ACLOCAL || echo "$ACLOCAL failed" +$AUTOHEADER || echo "$AUTOHEADER failed" +$AUTOMAKE --add-missing --force-missing --copy || echo "$AUTOMAKE failed" +$AUTOCONF || echo "$AUTOCONF failed" + +cd $ORIGDIR + +run_configure=true +for arg in $*; do + case $arg in + --no-configure) + run_configure=false + ;; + *) + ;; + esac +done + +if $run_configure; then + $srcdir/configure --enable-maintainer-mode "$@" + echo + echo "Now type 'make' to compile $PROJECT." +else + echo + echo "Now run 'configure' and 'make' to compile $PROJECT." +fi Index: /branches/eam_branches/gDiff.20200202/configure.ac =================================================================== --- /branches/eam_branches/gDiff.20200202/configure.ac (revision 41260) +++ /branches/eam_branches/gDiff.20200202/configure.ac (revision 41260) @@ -0,0 +1,47 @@ +dnl Process this file with autoconf to produce a configure script. +AC_PREREQ(2.61) + +AC_INIT([gDiff], [0.1.1], [ipp-support@ifa.hawaii.edu]) +AC_CONFIG_SRCDIR([src]) + +AM_INIT_AUTOMAKE([1.6 foreign dist-bzip2]) +AM_CONFIG_HEADER([src/config.h]) +AM_MAINTAINER_MODE + +IPP_STDLDFLAGS + +AC_LANG(C) +AC_GNU_SOURCE +AC_PROG_CC_C99 +AC_PROG_INSTALL +AC_PROG_LIBTOOL +AC_SYS_LARGEFILE + +PKG_CHECK_MODULES([PSLIB], [pslib >= 1.0.0]) +PKG_CHECK_MODULES([PSMODULE], [psmodules >= 1.0.0]) +PKG_CHECK_MODULES([PPSTATS], [ppStats >= 1.0.0]) +PKG_CHECK_MODULES([PSPHOT], [psphot >= 0.9.0]) + +AC_PATH_PROG([ERRORCODES], [psParseErrorCodes], [missing]) +if test "$ERRORCODES" = "missing" ; then + AC_MSG_ERROR([psParseErrorCodes is required]) +fi + +dnl Set CFLAGS for build +IPP_STDOPTS +IPP_STDCFLAGS + +IPP_VERSION + +AC_SUBST([GDIFF_CFLAGS]) +AC_SUBST([GDIFF_LIBS]) + +dnl add to AC_CONFIG_FILES below +dnl Doxyfile + +AC_CONFIG_FILES([ + Makefile + src/Makefile +]) + +AC_OUTPUT Index: /branches/eam_branches/gDiff.20200202/src/Makefile.am =================================================================== --- /branches/eam_branches/gDiff.20200202/src/Makefile.am (revision 41260) +++ /branches/eam_branches/gDiff.20200202/src/Makefile.am (revision 41260) @@ -0,0 +1,64 @@ +bin_PROGRAMS = gDiff + +# Force recompilation of gDiffVersion.c, since it gets the version information + +# replace cp with below +# pslib-setsvnversion.pl GDIFF gDiffVersionDefinitions.h.in gDiffVersionDefinitions.h +# cp gDiffVersionDefinitions.h.in gDiffVersionDefinitions.h +gDiffVersion.c: gDiffVersionDefinitions.h +gDiffVersionDefinitions.h: gDiffVersionDefinitions.h.in FORCE + ../pslib-setsvnversion.pl GDIFF gDiffVersionDefinitions.h.in gDiffVersionDefinitions.h +FORCE: ; + +gDiff_CPPFLAGS = $(PSLIB_CFLAGS) $(PSMODULE_CFLAGS) $(PPSTATS_CFLAGS) $(PSPHOT_CFLAGS) $(GDIFF_CFLAGS) +gDiff_LDFLAGS = $(PSLIB_LIBS) $(PSMODULE_LIBS) $(PPSTATS_LIBS) $(PSPHOT_LIBS) $(GDIFF_LIBS) + +gDiff_SOURCES = \ + gDiff.c \ + gDiffArguments.c \ + gDiffVersion.c \ + gDiffCamera.c \ + gDiffData.c \ + gDiffErrorCodes.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 + +noinst_HEADERS = \ + gDiff.h + + +### Error codes. +BUILT_SOURCES = gDiffErrorCodes.h gDiffErrorCodes.c +CLEANFILES = gDiffErrorCodes.h gDiffErrorCodes.c + +gDiffErrorCodes.h : gDiffErrorCodes.dat gDiffErrorCodes.h.in + $(ERRORCODES) --data=gDiffErrorCodes.dat --outdir=. gDiffErrorCodes.h + +gDiffErrorCodes.c : gDiffErrorCodes.dat gDiffErrorCodes.c.in gDiffErrorCodes.h + $(ERRORCODES) --data=gDiffErrorCodes.dat --outdir=. gDiffErrorCodes.c + +clean-local: + -rm -f TAGS + +# Tags for emacs +tags: + etags `find . -name \*.[ch] -print` Index: /branches/eam_branches/gDiff.20200202/src/config.h.in =================================================================== --- /branches/eam_branches/gDiff.20200202/src/config.h.in (revision 41260) +++ /branches/eam_branches/gDiff.20200202/src/config.h.in (revision 41260) @@ -0,0 +1,105 @@ +/* src/config.h.in. Generated from configure.ac by autoheader. */ + +/* Define to 1 if you have the header file. */ +#undef HAVE_DLFCN_H + +/* Define to 1 if you have the header file. */ +#undef HAVE_INTTYPES_H + +/* Define to 1 if you have the header file. */ +#undef HAVE_MEMORY_H + +/* Define to 1 if you have the header file. */ +#undef HAVE_STDINT_H + +/* Define to 1 if you have the header file. */ +#undef HAVE_STDLIB_H + +/* Define to 1 if you have the header file. */ +#undef HAVE_STRINGS_H + +/* Define to 1 if you have the header file. */ +#undef HAVE_STRING_H + +/* Define to 1 if you have the header file. */ +#undef HAVE_SYS_STAT_H + +/* Define to 1 if you have the header file. */ +#undef HAVE_SYS_TYPES_H + +/* Define to 1 if you have the header file. */ +#undef HAVE_UNISTD_H + +/* Define to the sub-directory in which libtool stores uninstalled libraries. + */ +#undef LT_OBJDIR + +/* Name of package */ +#undef PACKAGE + +/* Define to the address where bug reports for this package should be sent. */ +#undef PACKAGE_BUGREPORT + +/* Define to the full name of this package. */ +#undef PACKAGE_NAME + +/* Define to the full name and version of this package. */ +#undef PACKAGE_STRING + +/* Define to the one symbol short name of this package. */ +#undef PACKAGE_TARNAME + +/* Define to the home page for this package. */ +#undef PACKAGE_URL + +/* Define to the version of this package. */ +#undef PACKAGE_VERSION + +/* Define to 1 if you have the ANSI C header files. */ +#undef STDC_HEADERS + +/* Enable extensions on AIX 3, Interix. */ +#ifndef _ALL_SOURCE +# undef _ALL_SOURCE +#endif +/* Enable GNU extensions on systems that have them. */ +#ifndef _GNU_SOURCE +# undef _GNU_SOURCE +#endif +/* Enable threading extensions on Solaris. */ +#ifndef _POSIX_PTHREAD_SEMANTICS +# undef _POSIX_PTHREAD_SEMANTICS +#endif +/* Enable extensions on HP NonStop. */ +#ifndef _TANDEM_SOURCE +# undef _TANDEM_SOURCE +#endif +/* Enable general extensions on Solaris. */ +#ifndef __EXTENSIONS__ +# undef __EXTENSIONS__ +#endif + + +/* Version number of package */ +#undef VERSION + +/* Enable large inode numbers on Mac OS X 10.5. */ +#ifndef _DARWIN_USE_64_BIT_INODE +# define _DARWIN_USE_64_BIT_INODE 1 +#endif + +/* Number of bits in a file offset, on hosts where this is settable. */ +#undef _FILE_OFFSET_BITS + +/* Define for large files, on AIX-style hosts. */ +#undef _LARGE_FILES + +/* Define to 1 if on MINIX. */ +#undef _MINIX + +/* Define to 2 if the system does not provide POSIX.1 features except with + this defined. */ +#undef _POSIX_1_SOURCE + +/* Define to 1 if you need to in order for `stat' and other things to work. */ +#undef _POSIX_SOURCE Index: /branches/eam_branches/gDiff.20200202/src/gDiff.c =================================================================== --- /branches/eam_branches/gDiff.20200202/src/gDiff.c (revision 41260) +++ /branches/eam_branches/gDiff.20200202/src/gDiff.c (revision 41260) @@ -0,0 +1,93 @@ +/** @file gDiff.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" + +int main(int argc, char *argv[]) +{ + psLibInit(NULL); + psTimerStart("gDiff"); + + pmErrorRegister(); + gDiffErrorRegister(); + psphotErrorRegister(); + + gDiffData *data = NULL; // Processing data + pmConfig *config = pmConfigRead(&argc, argv, GDIFF_RECIPE); // Configuration + if (!config) goto die; + + gDiffVersionPrint(); // Print the version information + pmModelClassInit(); // returns false if already called, cannot otherwise fail + psphotInit(); // cannot return false + + data = gDiffDataAlloc(config); + + if (!gDiffArguments(argc, argv, data)) { + psError(psErrorCodeLast(), false, "Error reading arguments."); + goto die; + } + + // determine the camera and setup the I/O file handles + if (!gDiffCamera(data)) goto die; + + // actually calculate the difference: + if (!gDiffLoop(data)) goto die; + + die: + { + psExit exitValue = gDiffExitCode(PS_EXIT_SUCCESS); // 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 (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); + } + + 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(); + + 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 41260) +++ /branches/eam_branches/gDiff.20200202/src/gDiff.h (revision 41260) @@ -0,0 +1,198 @@ +/** @file gDiff.h + * + * @brief + * + * @ingroup gDiff + * + * @author IfA + * @version $Revision: 1.7 $ $Name: not supported by cvs2svn $ + * @date $Date: 2009-02-18 00:31:20 $ + * Copyright 2009 Institute for Astronomy, University of Hawaii + */ + +#ifndef PP_SUB_H +#define PP_SUB_H + +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include +#include +#include + +#include "gDiffErrorCodes.h" + +/// @addtogroup gDiff +/// @{ + +#define GDIFF_RECIPE "GDIFF" /// Name of the recipe to use +#define WCS_TOLERANCE 0.001 // Tolerance for WCS + +// Output files, for activation/deactivation +typedef enum { + GDIFF_FILES_INPUT = 0x01, // Input files + GDIFF_FILES_CONV = 0x02, // Convolved files (output) + GDIFF_FILES_SUB = 0x04, // Subtracted files (output) + GDIFF_FILES_INV = 0x08, // Inverse subtracted files (output) + GDIFF_FILES_PSF = 0x10, // PSF files (output) + GDIFF_FILES_PHOT_SUB = 0x20, // Subtraction photometry files (output) + GDIFF_FILES_PHOT_INV = 0x40, // Inverse subtraction photometry files (output) + GDIFF_FILES_PHOT = 0x80, // General photometry files (internal) + GDIFF_FILES_ALL = 0xFF, // All files +} gDiffFiles; + +/// Data for processing +typedef struct { + pmConfig *config; // Configuration + psErrorCode quality; // Quality code; 0 for no problem + bool photometry; // Perform photometry? + bool inverse; // Output inverse subtraction as well? + bool forcedPhot1; // perform forced photometry? + bool forcedPhot2; // perform forced photometry? + bool saveInConv; // Save convolved input? + bool saveRefConv; // Save convolved reference? + psString stamps; // Stamps file + pmPSF *psf; // Point Spread Function + psString statsName; // Name of statistics file + FILE *statsFile; // Statistics file + psMetadata *stats; // Statistics +} gDiffData; + +/// Constructor +gDiffData *gDiffDataAlloc(pmConfig *config ///< Configuration + ); + +/// Setup the arguments parsing +bool gDiffArguments(int argc, char *argv[], ///< Command-line arguments + gDiffData *data ///< Processing data + ); + +/// Parse the camera input +bool gDiffCamera(gDiffData *data ///< Processing data + ); + +/// Loop over the FPA hierarchy +bool gDiffLoop(gDiffData *data ///< Processing data + ); + +/// Generate (if needed) and set or update the masks for input and reference images +bool gDiffSetMasks(pmConfig *config ///< Configuration + ); + +// Calculate mask statistics +bool gDiffMaskStats(pmConfig *config, pmFPAview *view, psMetadata *stats); + +/// Generate the PSF-matching kernel and convolve the images as needed. Most of this function involves +/// looking up the parameters in the recipe and supplying them to the function pmSubtractionMatch() +bool gDiffMatchPSFs(gDiffData *data ///< Processing data + ); + +/// Threshold low pixels in image +bool gDiffLowThreshold(gDiffData *data ///< Processing data + ); + +/// Generate the output readout and pass the kernel info to the header +bool gDiffDefineOutput(const char *name,///< Name of output to define + pmConfig *config ///< Configuration + ); + +/// Photometry stage 1: measure the PSF from the minuend image +bool gDiffMakePSF(gDiffData *data ///< Processing data + ); + +/// Perform the actual image subtraction, update output concepts +bool gDiffReadoutSubtract(pmConfig *config ///< Configuration + ); + + +/// Photometry stage 2: find and measure sources on the subtracted image +bool gDiffReadoutPhotometry(const char *name, ///< Name of file to photometer + gDiffData *data ///< Processing data + ); + +bool gDiffInputDetections (bool *foundDetections, const char *sourcesName, const char *imageName, gDiffData *data); +bool gDiffReadoutForcedPhot(const char *outputName, const char *targetName, const char *sourceName, gDiffData *data); +bool psphotCopyResults (bool *foundDetections, pmFPAfile *target, pmFPAfile *source, pmFPAview *view); + +/// Higher-order background subtraction +bool gDiffBackground(pmConfig *config ///< Configuration + ); + +/// Perform Variance correction (rescale within a modest range) +bool gDiffVarianceRescale(pmConfig *config, ///< Configuration + gDiffData *data ///< Processing data + ); + +/// Put the program version information into a header +bool gDiffVersionHeader(psMetadata *header ///< Header to populate + ); + +/// Print version information +void gDiffVersionPrint(void); + +/// write the version info to a string +psString gDiffVersionLong(void); + +/// Mark the data quality as bad and prepare to suspend processing +void gDiffDataQuality(gDiffData *data, ///< Processing data + psErrorCode error,///< Error code + gDiffFiles files ///< Files to deactivate + ); + + +/// Activate or deactivate files +void gDiffFilesActivate(pmConfig *config, // Configuration + gDiffFiles files, // File to activate/deactivate + bool state // Activation state + ); + +/// Generate a view suitable for a readout +/// +/// Assumes we're working with skycells +pmFPAview *gDiffViewReadout(void); + +/// Iterate down the FPA hierarchy, opening files +bool gDiffFilesIterateDown(pmConfig *config, // Configuration + gDiffFiles files // Files to open + ); + +/// Iterate up the FPA hierarchy, closing files +bool gDiffFilesIterateUp(pmConfig *config, // Configuration + gDiffFiles files // Files to open + ); + +/// Collect statistics +bool gDiffReadoutStats(gDiffData *data // Processing data + ); + +/// Generate JPEG images +bool gDiffReadoutJpeg(pmConfig *config // Configuration + ); + +/// Generate JPEG images +bool gDiffResidualSampleJpeg(pmConfig *config); + +/// Generate inverse subtraction +bool gDiffReadoutInverse(pmConfig *config // Configuration + ); + + +// Copy every instance of a single keyword from one metadata to another +bool psMetadataCopySingle(psMetadata *target, psMetadata *source, const char *name); + +bool gDiffCopyPSF (pmFPAfile *output, pmFPAfile *input, pmFPAview *view); + +/// Return appropriate exit code +psExit gDiffExitCode(psExit exitValue // Current exit value + ); + +bool gDiffFlagNeighbors(pmConfig *config, pmFPAview *view, psArray *sources, bool matchRef); +bool gDiffMatchSources (psArray *objects, psArray *sources, float RADIUS, float MIN_SN); +bool gDiffSetSourceImageIDs (psArray *sources, int imageID); + +void gDiffSetThreads (void); + +///@} +#endif Index: /branches/eam_branches/gDiff.20200202/src/gDiffArguments.c =================================================================== --- /branches/eam_branches/gDiff.20200202/src/gDiffArguments.c (revision 41260) +++ /branches/eam_branches/gDiff.20200202/src/gDiffArguments.c (revision 41260) @@ -0,0 +1,159 @@ +#include "gDiff.h" + +// Print usage information and die +static void usage(const char *program, // Name of the program + psMetadata *arguments, // Command-line arguments + pmConfig *config // Configuration + ) +{ + fprintf(stderr, "\nPan-STARRS PSF-matched image subtraction\n\n"); + fprintf(stderr, "Usage: %s OUTPUT_ROOT \n" + "\t[-psf REFERENCE.psf.fits]\n\n" + "This subtracts the convolved REFERENCE from the INPUT, by default.\n\n", + program); + psArgumentHelp(arguments); + psFree(config); + pmConfigDone(); + psLibFinalize(); + exit(PS_EXIT_CONFIG_ERROR); +} + +/** + * Add a single filename to the arguments as an array, so that it can be used with pmFPAfileBindFromArgs, etc + */ +static void fileList(const char *file, // The symbolic name for the file + const char *name, // The name of the file + const char *comment, // Description of the file + pmConfig *config // Configuration + ) +{ + psArray *files = psArrayAlloc(1); // Array with file names + files->data[0] = psStringCopy(name); + psMetadataAddArray(config->arguments, PS_LIST_TAIL, file, 0, comment, files); + psFree(files); + return; +} + +void gDiffSetThreads (void) { + // gDiff does not have any of its own thread handlers + return; +} + +bool gDiffArguments(int argc, char *argv[], gDiffData *data) +{ + assert(data); + pmConfig *config = data->config; + assert(config); + + // generic arguments (version, dumpconfig) + PS_ARGUMENTS_GENERIC( gDiff, config, argc, argv ); + + // thread arguments + PS_ARGUMENTS_THREADS( gDiff, config, argc, argv ) + + int argNum = psArgumentGet(argc, argv, "-debug"); // Debugging argument number + if (argNum) { + psArgumentRemove(argNum, &argc, argv); + pmSubtractionRegions(true); + } + + + 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); + 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); + psMetadataAddBool(arguments, PS_LIST_TAIL, "-inverse", 0, "Generate inverse subtractions?", false); + psMetadataAddBool(arguments, PS_LIST_TAIL, "-save-inconv", 0, "Save input convolved images?", false); + psMetadataAddBool(arguments, PS_LIST_TAIL, "-save-refconv", 0, "Save reference convolved images?", false); + psMetadataAddBool(arguments, PS_LIST_TAIL, "-visual", 0, "Show diagnostic plots", NULL); + psMetadataAddBool(arguments, PS_LIST_TAIL, "-updatemode", 0, "update mode?", false); + + if (argc == 1 || !psArgumentParse(arguments, &argc, argv) || argc != 2) { + usage(argv[0], arguments, config); + } + + psMetadataAddStr(config->arguments, PS_LIST_TAIL, "OUTPUT", 0, "Name of the output image", argv[1]); + + + const char *inImage = psMetadataLookupStr(NULL, arguments, "-inimage"); // Name of input image + if (inImage && strlen(inImage) > 0) { + fileList("INPUT", inImage, "Name of the input image", config); + } + const char *inMask = psMetadataLookupStr(NULL, arguments, "-inmask"); // Name of input mask + if (inMask && strlen(inMask) > 0) { + fileList("INPUT.MASK", inMask, "Name of the input mask image", config); + } + const char *inVariance = psMetadataLookupStr(NULL, arguments, "-invariance"); // Name of input variance + if (inVariance && strlen(inVariance) > 0) { + fileList("INPUT.VARIANCE", inVariance, "Name of the input variance image", config); + } + const char *inSources = psMetadataLookupStr(NULL, arguments, "-insources"); // Name of input source list + if (inSources && strlen(inSources) > 0) { + fileList("INPUT.SOURCES", inSources, "Name of the input source list", config); + } + + const char *refImage = psMetadataLookupStr(NULL, arguments, "-refimage"); // Name of reference image + if (refImage && strlen(refImage) > 0) { + fileList("REF", refImage, "Name of the reference image", config); + } + const char *refMask = psMetadataLookupStr(NULL, arguments, "-refmask"); // Name of reference mask + if (refMask && strlen(refMask) > 0) { + fileList("REF.MASK", refMask, "Name of the reference mask image", config); + } + const char *refVariance = psMetadataLookupStr(NULL, arguments, "-refvariance"); // Name of ref variance + if (refVariance && strlen(refVariance) > 0) { + fileList("REF.VARIANCE", refVariance, "Name of the reference variance image", config); + } + const char *refSources = psMetadataLookupStr(NULL, arguments, "-refsources"); // Name of ref source list + if (refSources && strlen(refSources) > 0) { + fileList("REF.SOURCES", refSources, "Name of the reference source list", config); + } + + const char *kernel = psMetadataLookupStr(NULL, arguments, "-kernel"); // Name of kernel + if (kernel && strlen(kernel) > 0) { + fileList("KERNEL", kernel, "Name of the kernel to apply", config); + } + + data->stamps = psMemIncrRefCounter(psMetadataLookupStr(NULL, arguments, "-stamps")); + + data->statsName = psMemIncrRefCounter(psMetadataLookupStr(NULL, arguments, "-stats")); + if (data->statsName && strlen(data->statsName) > 0) { + psString resolved = pmConfigConvertFilename(data->statsName, config, true, true); // Resolved filename + if (!resolved) { + psError(psErrorCodeLast(), false, "Unable to resolve statistics file: %s", data->statsName); + return false; + } + data->statsFile = fopen(resolved, "w"); + if (!data->statsFile) { + psError(PS_ERR_IO, true, "Unable to open statistics file %s for writing.\n", resolved); + psFree(resolved); + return false; + } + psFree(resolved); + } + + data->saveInConv = psMetadataLookupBool(NULL, arguments, "-save-inconv"); + data->saveRefConv = psMetadataLookupBool(NULL, arguments, "-save-refconv"); + + if (psMetadataLookupBool(NULL, arguments, "-visual")) { + pmVisualSetVisual(true); + } + + psTrace("gDiff", 1, "Done reading command-line arguments\n"); + + return true; +} Index: /branches/eam_branches/gDiff.20200202/src/gDiffCamera.c =================================================================== --- /branches/eam_branches/gDiff.20200202/src/gDiffCamera.c (revision 41260) +++ /branches/eam_branches/gDiff.20200202/src/gDiffCamera.c (revision 41260) @@ -0,0 +1,477 @@ +#include "gDiff.h" + +bool checkFileruleFileSave (pmFPAfile *file, pmConfig *config) { + + bool status; + + psMetadata *filerule = pmConfigFileRule(config, file->camera, file->name); // File rule + if (!filerule) return false; + + char *myString = psMetadataLookupStr(&status, filerule, "FILE.SAVE"); + if (!myString) return false; + + // do not change the value from the default unless TRUE or FALSE are found + if (!strcasecmp(myString, "TRUE")) { + file->save = true; + return true; + } + if (!strcasecmp(myString, "FALSE")) { + file->save = false; + return true; + } + + return false; +} + +// Define an input file +static pmFPAfile *defineInputFile(bool *success, + pmConfig *config,// Configuration + pmFPAfile *bind, // File to which to bind, or NULL + char *filerule, // Name of file rule + char *argname, // Argument name + pmFPAfileType fileType // Type of file + ) +{ + bool status; + + *success = false; + + pmFPAfile *file = NULL; + + // look for the file on the argument list + if (bind) { + file = pmFPAfileBindFromArgs(&status, bind, config, filerule, argname); + } else { + file = pmFPAfileDefineFromArgs(&status, config, filerule, argname); + } + + if (!status) { + psError(psErrorCodeLast(), false, "Failed to load file definition for %s", filerule); + return false; + } + if (!file) { + // look for the file on the RUN metadata + file = pmFPAfileDefineFromRun(&status, bind, config, filerule); // File to return + if (!status) { + psError(psErrorCodeLast(), false, "Failed to load file definition for %s", filerule); + return NULL; + } + } + + if (!file) { + // no file defined + *success = true; + return NULL; + } + + if (file->type != fileType) { + psError(GDIFF_ERR_CONFIG, true, "%s is not of type %s", filerule, pmFPAfileStringFromType(fileType)); + return NULL; + } + + *success = true; + return file; +} + +// Define an output file +static pmFPAfile *defineOutputFile(pmConfig *config, // Configuration + pmFPAfile *template, // File to use as basis for definition + bool source, // Is template a source (T), or for binding (F)? + char *filerule, // Name of file rule + pmFPAfileType fileType // Type of file + ) +{ + + pmFPAfile *file = source ? pmFPAfileDefineFromFile(config, template, 1, 1, filerule) : + pmFPAfileDefineOutput(config, template->fpa, filerule); + if (!file) { + psError(psErrorCodeLast(), false, _("Unable to generate output file from %s"), filerule); + return NULL; + } + if (file->type != fileType) { + psError(GDIFF_ERR_CONFIG, true, "%s is not of type %s", filerule, pmFPAfileStringFromType(fileType)); + return NULL; + } + + return file; +} + +// Define an output file that will be used in a calculation +// This means it might already be available in the RUN metadata +static pmFPAfile *defineCalcFile(pmConfig *config, // Configuration + pmFPAfile *bind, // File to which to bind, or NULL + char *filerule, // Name of file rule + const char *argname, // Argument name for file + pmFPAfileType fileType // Type of file + ) +{ + bool status; + + // look for the file on the RUN metadata + pmFPAfile *file = pmFPAfileDefineFromRun(&status, NULL, config, filerule); // File to return + if (!status) { + psError(psErrorCodeLast(), false, "Failed to load file definition for %s", filerule); + return NULL; + } + file = pmFPAfileBindFromArgs(&status, bind, config, filerule, argname); + if (!status) { + psError(psErrorCodeLast(), false, "Failed to load file definition for %s", filerule); + return NULL; + } + if (file) { + // It's an input + file->save = false; + } + + // define new version of file + if (!file) { + file = pmFPAfileDefineOutput(config, bind ? bind->fpa : NULL, filerule); + if (!status) { + psError(psErrorCodeLast(), false, "Failed to load file definition for %s", filerule); + return false; + } + if (file) { + // It's an output (respect filerules) + checkFileruleFileSave(file, config); + } + } + + if (!file) { + return NULL; + } + + if (file->type != fileType) { + psError(GDIFF_ERR_CONFIG, true, "%s is not of type %s", filerule, pmFPAfileStringFromType(fileType)); + return NULL; + } + + return file; +} + +bool gDiffCamera(gDiffData *data) +{ + bool success = true; + + psAssert(data, "Require processing data"); + pmConfig *config = data->config; + psAssert(config, "Require configuration"); + + // Input image + pmFPAfile *input = defineInputFile(&success, config, NULL, "GDIFF.INPUT", "INPUT", PM_FPA_FILE_IMAGE); + if (!success) { + psError(psErrorCodeLast(), false, "Failed to build FPA from GDIFF.INPUT"); + return false; + } + + defineInputFile(&success, config, input, "GDIFF.INPUT.MASK", "INPUT.MASK", PM_FPA_FILE_MASK); + if (!success) { + psError(psErrorCodeLast(), false, "Failed to build FPA from GDIFF.INPUT.MASK"); + return false; + } + + pmFPAfile *inVar = defineInputFile(&success, config, input, "GDIFF.INPUT.VARIANCE", "INPUT.VARIANCE", PM_FPA_FILE_VARIANCE); + if (!success) { + psError(psErrorCodeLast(), false, "Failed to build FPA from GDIFF.INPUT.VARIANCE"); + return false; + } + + defineInputFile(&success, config, NULL, "GDIFF.INPUT.SOURCES", "INPUT.SOURCES", PM_FPA_FILE_CMF); + if (!success) { + psError(psErrorCodeLast(), false, "Failed to build FPA from GDIFF.INPUT.SOURCES"); + return false; + } + + // Reference image + pmFPAfile *ref = defineInputFile(&success, config, NULL, "GDIFF.REF", "REF", PM_FPA_FILE_IMAGE); + if (!success) { + psError(psErrorCodeLast(), false, "Failed to build FPA from GDIFF.REF"); + return false; + } + + defineInputFile(&success, config, ref, "GDIFF.REF.MASK", "REF.MASK", PM_FPA_FILE_MASK); + if (!success) { + psError(psErrorCodeLast(), false, "Failed to build FPA from GDIFF.REF.MASK"); + return false; + } + + pmFPAfile *refVar = defineInputFile(&success, config, ref, "GDIFF.REF.VARIANCE", "REF.VARIANCE", PM_FPA_FILE_VARIANCE); + if (!success) { + psError(psErrorCodeLast(), false, "Failed to build FPA from GDIFF.REF.VARIANCE"); + return false; + } + + defineInputFile(&success, config, NULL, "GDIFF.REF.SOURCES", "REF.SOURCES", PM_FPA_FILE_CMF); + if (!success) { + psError(psErrorCodeLast(), false, "Failed to build FPA from GDIFF.REF.SOURCES"); + return false; + } + + // Now that the camera has been determined, we can read the recipe + psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, GDIFF_RECIPE); // Recipe for ppSim + if (!recipe) { + psError(GDIFF_ERR_CONFIG, false, "Unable to find recipe %s", GDIFF_RECIPE); + return false; + } + if (psMetadataLookupBool(NULL, config->arguments, "-photometry")) { + psMetadataAddBool(recipe, PS_LIST_TAIL, "PHOTOMETRY", PS_META_REPLACE, + "Perform photometry?", true); + } + if (psMetadataLookupBool(NULL, config->arguments, "-inverse")) { + psMetadataAddBool(recipe, PS_LIST_TAIL, "INVERSE", PS_META_REPLACE, + "Generate inverse subtractions?", true); + } + if (psMetadataLookupBool(NULL, config->arguments, "-forced-phot")) { + psMetadataAddBool(recipe, PS_LIST_TAIL, "FORCED.PHOTOMETRY.BOTH", PS_META_REPLACE, "Perform forced photometry?", true); + } + if (psMetadataLookupBool(NULL, config->arguments, "-forced-input1")) { + psMetadataAddBool(recipe, PS_LIST_TAIL, "FORCED.PHOTOMETRY.INPUT1", PS_META_REPLACE, "Perform forced photometry?", true); + } + if (psMetadataLookupBool(NULL, config->arguments, "-forced-input2")) { + psMetadataAddBool(recipe, PS_LIST_TAIL, "FORCED.PHOTOMETRY.INPUT2", PS_META_REPLACE, "Perform forced photometry?", true); + } + + data->inverse = psMetadataLookupBool(NULL, recipe, "INVERSE"); + data->photometry = psMetadataLookupBool(NULL, recipe, "PHOTOMETRY"); + data->forcedPhot1 = psMetadataLookupBool(NULL, recipe, "FORCED.PHOTOMETRY.BOTH") || psMetadataLookupBool(NULL, recipe, "FORCED.PHOTOMETRY.INPUT1"); + data->forcedPhot2 = psMetadataLookupBool(NULL, recipe, "FORCED.PHOTOMETRY.BOTH") || psMetadataLookupBool(NULL, recipe, "FORCED.PHOTOMETRY.INPUT2"); + + // Convolved input image + pmFPAfile *inConvImage = defineOutputFile(config, input, true, "GDIFF.INPUT.CONV", PM_FPA_FILE_IMAGE); + pmFPAfile *inConvMask = defineOutputFile(config, inConvImage, false, "GDIFF.INPUT.CONV.MASK", + PM_FPA_FILE_MASK); + if (!inConvImage || !inConvMask) { + psError(psErrorCodeLast(), false, "Unable to define output files"); + return false; + } + // these ->save values below are set by command-line arguments (in gDiffArguments.c) : -save-inconv, -save-refconv + checkFileruleFileSave (inConvImage, config); + inConvImage->save = inConvImage->save || data->saveInConv; + + checkFileruleFileSave (inConvMask, config); + inConvMask->save = inConvMask->save || data->saveInConv; + if (inVar) { + pmFPAfile *inConvVar = defineOutputFile(config, inConvImage, false, "GDIFF.INPUT.CONV.VARIANCE", + PM_FPA_FILE_VARIANCE); + if (!inConvVar) { + psError(psErrorCodeLast(), false, "Unable to define output files"); + return false; + } + checkFileruleFileSave (inConvVar, config); + inConvVar->save = inConvVar->save || data->saveInConv; + } + + // Convolved ref image + pmFPAfile *refConvImage = defineOutputFile(config, input, true, "GDIFF.REF.CONV", PM_FPA_FILE_IMAGE); + pmFPAfile *refConvMask = defineOutputFile(config, refConvImage, false, "GDIFF.REF.CONV.MASK", + PM_FPA_FILE_MASK); + if (!refConvImage || !refConvMask) { + psError(psErrorCodeLast(), false, "Unable to define output files"); + return false; + } + checkFileruleFileSave (refConvImage, config); + refConvImage->save = refConvImage->save || data->saveRefConv; + + checkFileruleFileSave (refConvMask, config); + refConvMask->save = refConvMask->save || data->saveRefConv; + if (refVar) { + pmFPAfile *refConvVar = defineOutputFile(config, refConvImage, false, "GDIFF.REF.CONV.VARIANCE", + PM_FPA_FILE_VARIANCE); + if (!refConvVar) { + psError(psErrorCodeLast(), false, "Unable to define output files"); + return false; + } + checkFileruleFileSave (refConvVar, config); + refConvVar->save = refConvVar->save || data->saveRefConv; + } + + // Output image + pmFPAfile *output = defineOutputFile(config, inConvImage, true, "GDIFF.OUTPUT", PM_FPA_FILE_IMAGE); + pmFPAfile *outMask = defineOutputFile(config, output, false, "GDIFF.OUTPUT.MASK", PM_FPA_FILE_MASK); + if (!output || !outMask) { + psError(psErrorCodeLast(), false, "Unable to define output files"); + return false; + } + checkFileruleFileSave (output, config); + checkFileruleFileSave (outMask, config); + pmFPAfile *outVar = NULL; + if (inVar && refVar) { + outVar = defineOutputFile(config, output, false, "GDIFF.OUTPUT.VARIANCE", + PM_FPA_FILE_VARIANCE); + if (!outVar) { + psError(psErrorCodeLast(), false, "Unable to define output files"); + return false; + } + checkFileruleFileSave (outVar, config); + } + // If we are in update mode unconditionally save the output files + bool updateMode = psMetadataLookupBool(NULL, config->arguments, "-updatemode"); + if (updateMode) { + output->save = true; + outMask->save = true; + if (outVar) { + outVar->save = true; + } + } + + pmFPAfile *inverse = NULL; // Inverse output image + if (data->inverse) { + // Inverse output image + inverse = defineOutputFile(config, output, true, "GDIFF.INVERSE", PM_FPA_FILE_IMAGE); + pmFPAfile *invMask = defineOutputFile(config, inverse, false, "GDIFF.INVERSE.MASK", + PM_FPA_FILE_MASK); + if (!inverse || !invMask) { + psError(psErrorCodeLast(), false, "Unable to define output files"); + return false; + } + checkFileruleFileSave (inverse, config); + checkFileruleFileSave (invMask, config); + if (inVar && refVar) { + pmFPAfile *invVar = defineOutputFile(config, inverse, false, "GDIFF.INVERSE.VARIANCE", + PM_FPA_FILE_VARIANCE); + if (!invVar) { + psError(psErrorCodeLast(), false, "Unable to define output files"); + return false; + } + checkFileruleFileSave(invVar, config); + } + } + + + // Output JPEGs + pmFPAfile *jpeg1 = pmFPAfileDefineOutput(config, NULL, "GDIFF.OUTPUT.JPEG1"); + if (!jpeg1) { + psError(psErrorCodeLast(), false, _("Unable to generate output file from GDIFF.OUTPUT.JPEG1")); + return false; + } + if (jpeg1->type != PM_FPA_FILE_JPEG) { + psError(psErrorCodeLast(), true, "GDIFF.OUTPUT.JPEG1 is not of type JPEG"); + return false; + } + checkFileruleFileSave(jpeg1, config); + pmFPAfile *jpeg2 = pmFPAfileDefineOutput(config, NULL, "GDIFF.OUTPUT.JPEG2"); + if (!jpeg2) { + psError(psErrorCodeLast(), false, _("Unable to generate output file from GDIFF.OUTPUT.JPEG2")); + return false; + } + if (jpeg2->type != PM_FPA_FILE_JPEG) { + psError(psErrorCodeLast(), true, "GDIFF.OUTPUT.JPEG2 is not of type JPEG"); + return false; + } + checkFileruleFileSave(jpeg2, config); + + // Output residual JPEG + pmFPAfile *jpeg3 = pmFPAfileDefineOutput(config, NULL, "GDIFF.OUTPUT.RESID.JPEG"); + if (!jpeg3) { + psError(psErrorCodeLast(), false, _("Unable to generate output file from GDIFF.OUTPUT.RESID.JPEG")); + return false; + } + if (jpeg3->type != PM_FPA_FILE_JPEG) { + psError(psErrorCodeLast(), true, "GDIFF.OUTPUT.RESID.JPEG is not of type JPEG"); + return false; + } + checkFileruleFileSave(jpeg3, config); + + // Output subtraction kernel + pmFPAfile *kernel = defineCalcFile(config, output, "GDIFF.OUTPUT.KERNELS", "KERNEL", + PM_FPA_FILE_SUBKERNEL); + if (!kernel) { + psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to define file GDIFF.OUTPUT.KERNELS"); + return false; + } + + // 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); + } + + 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); + } + } + + return true; +} + + Index: /branches/eam_branches/gDiff.20200202/src/gDiffData.c =================================================================== --- /branches/eam_branches/gDiff.20200202/src/gDiffData.c (revision 41260) +++ /branches/eam_branches/gDiff.20200202/src/gDiffData.c (revision 41260) @@ -0,0 +1,52 @@ +#include "gDiff.h" + + +static void gDiffDataFree(gDiffData *data) +{ + psAssert(!data->statsFile, "Statistics file still open."); + psFree(data->statsName); + psFree(data->stamps); + psFree(data->psf); + psFree(data->stats); + psFree(data->config); + + return; +} + +gDiffData *gDiffDataAlloc(pmConfig *config) +{ + gDiffData *data = psAlloc(sizeof(gDiffData)); // Processing data, to return + psMemSetDeallocator(data, (psFreeFunc)gDiffDataFree); + + data->config = config; + data->quality = 0; + data->photometry = false; + data->inverse = false; + data->saveInConv = false; + data->saveRefConv = false; + data->stamps = NULL; + data->psf = NULL; + data->statsName = NULL; + data->statsFile = NULL; + data->stats = psMetadataAlloc(); + psMetadataAddS32(data->stats, PS_LIST_TAIL, "QUALITY", 0, "Data quality", 0); + + return data; +} + + +void gDiffDataQuality(gDiffData *data, psErrorCode error, gDiffFiles files) +{ + psAssert(data, "Require processing data"); + + if (psMetadataLookupS32(NULL, data->stats, "QUALITY") == 0) { + data->quality = error; + psMetadataAddS32(data->stats, PS_LIST_TAIL, "QUALITY", PS_META_REPLACE, "Data quality", error); + } + + gDiffFilesActivate(data->config, files, false); + + psErrorClear(); + + return; +} Index: /branches/eam_branches/gDiff.20200202/src/gDiffDefineOutput.c =================================================================== --- /branches/eam_branches/gDiff.20200202/src/gDiffDefineOutput.c (revision 41260) +++ /branches/eam_branches/gDiff.20200202/src/gDiffDefineOutput.c (revision 41260) @@ -0,0 +1,283 @@ +/** @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 + */ + +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include +#include +#include + +#include "gDiff.h" + +bool gDiffCopyWarpToChip (psMetadata *tgtHeader, psMetadata *srcHeader); + +bool gDiffDefineOutput(const char *name, pmConfig *config) +{ + psAssert(name, "Require name"); + psAssert(config, "Require configuration"); + + bool status = false; + + pmFPAview *view = gDiffViewReadout(); // View to readout + pmCell *outCell = pmFPAfileThisCell(config->files, view, name); // Output cell + pmFPA *outFPA = outCell->parent->parent; // Output FPA + pmHDU *outHDU = outFPA->hdu; // Output HDU + if (!outHDU->header) { + outHDU->header = psMetadataAlloc(); + } + + // The output readout may already be present if we read in the convolution kernel + pmReadout *outRO = NULL; // Output readout + if (outCell->readouts && outCell->readouts->n > 0 && outCell->readouts->data[0]) { + outRO = psMemIncrRefCounter(outCell->readouts->data[0]); + } else { + outRO = pmReadoutAlloc(outCell); + } + + // Convolved input images + psMetadata *recipe = psMetadataLookupPtr(&status, config->recipes, GDIFF_RECIPE); + bool noConvolve = psMetadataLookupBool(&status, recipe, "NOCONVOLVE"); // Do not use convolved images. + bool reverse = psMetadataLookupBool(&status, config->arguments, "REVERSE"); // Reverse sense of subtraction? + bool mapFromPositive = psMetadataLookupBool(&status, config->arguments, "CHIP.MAP.FROM.POSITIVE"); + // mapFromPositive = true means "always grab the warp->chip map from the positive image (eg, A for A - B, B for B - A) + // mapFromPositive = false means "always grab the warp->chip map from the first image (eg, GDIFF.INPUT if not 'reverse') + + pmReadout *inConv; + if (noConvolve) { + inConv = pmFPAfileThisReadout(config->files, view, "GDIFF.INPUT"); // Input readout + } else { + inConv = pmFPAfileThisReadout(config->files, view, "GDIFF.INPUT.CONV"); // Input readout + } + pmReadout *refConv; + if (noConvolve) { + refConv = pmFPAfileThisReadout(config->files, view, "GDIFF.REF"); // Reference readout + } else { + refConv = pmFPAfileThisReadout(config->files, view, "GDIFF.REF.CONV"); // Reference readout + } + + // Add kernel descrption to header. + // We don't know which readout has the kernels because it depends on which PSF is larger + bool mdok; // Status of MD lookup + psMetadata *analysis = inConv->analysis; // Analysis metadata with kernel information + pmHDU *hdu = pmHDUFromCell(inConv->parent); + pmSubtractionKernels *kernels = psMetadataLookupPtr(&mdok, analysis, PM_SUBTRACTION_ANALYSIS_KERNEL); // Subtraction kernel + + if (!kernels) { + hdu = pmHDUFromCell(refConv->parent); + analysis = refConv->analysis; + kernels = psMetadataLookupPtr(&mdok, analysis, PM_SUBTRACTION_ANALYSIS_KERNEL); + } + if (!kernels && !noConvolve) { + psError(GDIFF_ERR_UNKNOWN, true, "Unable to find SUBTRACTION.KERNEL"); + psFree(outRO); + psFree(view); + return false; + } + psAssert (hdu, "unable to find HDU"); + if (!noConvolve) { + psMetadataAddStr(outHDU->header, PS_LIST_TAIL, "GDIFF.KERNEL", 0, "Subtraction kernel", kernels->description); + } + outHDU->header = psMetadataCopy(outHDU->header, hdu->header); + + // in warp-stack mode, we should always use the warp + // in warp-warp mode, we should use the positive warp for the positive subtraction + + bool normalDiff = true; + if (!strcmp (name, "GDIFF.INVERSE")) { + normalDiff = false; + } + + if (reverse) { + // normal = GDIFF.REF - GDIFF.INPUT + // inverse = GDIFF.INPUT - GDIFF.REF + if (mapFromPositive) { + pmCell *cell_Pos = normalDiff ? pmFPAfileThisCell(config->files, view, "GDIFF.REF") : pmFPAfileThisCell(config->files, view, "GDIFF.INPUT"); + pmHDU *hdu_Pos = pmHDUFromCell(cell_Pos); + gDiffCopyWarpToChip (outHDU->header, hdu_Pos->header); + } else { + pmCell *cell_Pos = pmFPAfileThisCell(config->files, view, "GDIFF.REF"); + pmHDU *hdu_Pos = pmHDUFromCell(cell_Pos); + gDiffCopyWarpToChip (outHDU->header, hdu_Pos->header); + } + } else { + // normal = GDIFF.INPUT - GDIFF.REF + // inverse = GDIFF.REF - GDIFF.INPUT + if (mapFromPositive) { + pmCell *cell_Pos = normalDiff ? pmFPAfileThisCell(config->files, view, "GDIFF.INPUT") : pmFPAfileThisCell(config->files, view, "GDIFF.REF"); + pmHDU *hdu_Pos = pmHDUFromCell(cell_Pos); + gDiffCopyWarpToChip (outHDU->header, hdu_Pos->header); + } else { + pmCell *cell_Pos = pmFPAfileThisCell(config->files, view, "GDIFF.INPUT"); + pmHDU *hdu_Pos = pmHDUFromCell(cell_Pos); + gDiffCopyWarpToChip (outHDU->header, hdu_Pos->header); + } + } + + // Add additional data to the header + pmFPAfile *refFile = psMetadataLookupPtr(NULL, config->files, "GDIFF.REF"); // Reference file + pmFPAfile *inFile = psMetadataLookupPtr(NULL, config->files, "GDIFF.INPUT"); // Input file + + // save the names of the input and reference image in the header + psString refBase = psStringFileBasename(refFile->origname); // Basename of reference + psMetadataAddStr(outHDU->header, PS_LIST_TAIL, "GDIFF.REFERENCE", 0, "Subtraction reference", refBase); + psFree(refBase); + + psString inBase = psStringFileBasename(inFile->origname); // Basename of input + psMetadataAddStr(outHDU->header, PS_LIST_TAIL, "GDIFF.INPUT", 0, "Subtraction input", inBase); + psFree(inBase); + + gDiffVersionHeader(outHDU->header); + + outRO->analysis = psMetadataCopy(outRO->analysis, analysis); + + psFree(outRO); + psFree(view); + + return true; +} + +// we have 4 sets of header keywords to copy: +// SRC_nnnn, SEC_nnnn, MPX_nnnn, MPY_nnnn + +bool gDiffCopyWarpToChip (psMetadata *tgtHeader, psMetadata *srcHeader) { + + char keyword[80]; + + bool status = false; + + int Nchip = 0; + while (true) { + snprintf (keyword, 80, "SRC_%04d", Nchip); + char *string = psMetadataLookupStr (&status, srcHeader, keyword); + if (!status) { + break; + } + psMetadataAddStr(tgtHeader, PS_LIST_TAIL, keyword, PS_META_REPLACE, "input image", string); + Nchip ++; + } + + for (int i = 0; i < Nchip; i++) { + snprintf (keyword, 80, "SEC_%04d", i); + char *string = psMetadataLookupStr (&status, srcHeader, keyword); + if (!status) { + psWarning ("cannot find keyword %s\n", keyword); + continue; + } + psMetadataAddStr(tgtHeader, PS_LIST_TAIL, keyword, PS_META_REPLACE, "input image", string); + } + + for (int i = 0; i < Nchip; i++) { + snprintf (keyword, 80, "MPX_%04d", i); + char *string = psMetadataLookupStr (&status, srcHeader, keyword); + if (!status) { + psWarning ("cannot find keyword %s\n", keyword); + continue; + } + psMetadataAddStr(tgtHeader, PS_LIST_TAIL, keyword, PS_META_REPLACE, "input image", string); + } + + for (int i = 0; i < Nchip; i++) { + snprintf (keyword, 80, "MPY_%04d", i); + char *string = psMetadataLookupStr (&status, srcHeader, keyword); + if (!status) { + psWarning ("cannot find keyword %s\n", keyword); + continue; + } + psMetadataAddStr(tgtHeader, PS_LIST_TAIL, keyword, PS_META_REPLACE, "input image", string); + } + + // Copy PSREFCAT from the source to the target as well + snprintf(keyword, 80, "PSREFCAT"); + char *string = psMetadataLookupStr(&status, srcHeader, keyword); + if (status) { + psMetadataAddStr(tgtHeader, PS_LIST_TAIL, keyword, PS_META_REPLACE, "input image reference catalog", string); + } + + return true; +} + +bool gDiffCopyWarpToChip_Alt (psMetadata *tgtHeader, psMetadata *srcHeader, bool isPositive) { + + char srcKeyword[80], tgtKeyword[80]; + + bool status = false; + + int Nchip = 0; + while (true) { + snprintf (srcKeyword, 80, "SRC_%04d", Nchip); + if (isPositive) { + snprintf (tgtKeyword, 80, "SRCP_%03d", Nchip); + } else { + snprintf (tgtKeyword, 80, "SRCM_%03d", Nchip); + } + + char *string = psMetadataLookupStr (&status, srcHeader, srcKeyword); + if (!status) { + break; + } + psMetadataAddStr(tgtHeader, PS_LIST_TAIL, srcKeyword, PS_META_REPLACE, "input image", string); + Nchip ++; + } + + for (int i = 0; i < Nchip; i++) { + snprintf (srcKeyword, 80, "SEC_%04d", i); + if (isPositive) { + snprintf (tgtKeyword, 80, "SECP_%03d", Nchip); + } else { + snprintf (tgtKeyword, 80, "SECM_%03d", Nchip); + } + + char *string = psMetadataLookupStr (&status, srcHeader, srcKeyword); + if (!status) { + psWarning ("cannot find keyword %s\n", srcKeyword); + continue; + } + psMetadataAddStr(tgtHeader, PS_LIST_TAIL, srcKeyword, PS_META_REPLACE, "input image", string); + } + + for (int i = 0; i < Nchip; i++) { + snprintf (srcKeyword, 80, "MPX_%04d", i); + if (isPositive) { + snprintf (tgtKeyword, 80, "MPXP_%03d", Nchip); + } else { + snprintf (tgtKeyword, 80, "MPXM_%03d", Nchip); + } + + char *string = psMetadataLookupStr (&status, srcHeader, srcKeyword); + if (!status) { + psWarning ("cannot find keyword %s\n", srcKeyword); + continue; + } + psMetadataAddStr(tgtHeader, PS_LIST_TAIL, srcKeyword, PS_META_REPLACE, "input image", string); + } + + for (int i = 0; i < Nchip; i++) { + snprintf (srcKeyword, 80, "MPY_%04d", i); + if (isPositive) { + snprintf (tgtKeyword, 80, "MPYP_%03d", Nchip); + } else { + snprintf (tgtKeyword, 80, "MPYM_%03d", Nchip); + } + + char *string = psMetadataLookupStr (&status, srcHeader, srcKeyword); + if (!status) { + psWarning ("cannot find keyword %s\n", srcKeyword); + continue; + } + psMetadataAddStr(tgtHeader, PS_LIST_TAIL, srcKeyword, PS_META_REPLACE, "input image", string); + } + return true; +} + Index: /branches/eam_branches/gDiff.20200202/src/gDiffErrorCodes.c.in =================================================================== --- /branches/eam_branches/gDiff.20200202/src/gDiffErrorCodes.c.in (revision 41260) +++ /branches/eam_branches/gDiff.20200202/src/gDiffErrorCodes.c.in (revision 41260) @@ -0,0 +1,35 @@ +/** @file gDiffErrorCodes.c.in + * @brief + * @ingroup gDiff + * + * @author IfA + * @version $Revision: 1.3 $ $Name: not supported by cvs2svn $ + * @date $Date: 2009-02-05 20:44:04 $ + * Copyright 2009 Institute for Astronomy, University of Hawaii + */ + +#include "pslib.h" +#include "gDiffErrorCodes.h" + +/* + * The line + { GDIFF_ERR_$X{ErrorCode}, "$X{ErrorDescription}"}, + * (without the Xs) + * will be replaced by values from errorCodes.dat + */ +void gDiffErrorRegister(void) +{ + static psErrorDescription errors[] = { + { GDIFF_ERR_BASE, "First value we use; lower values belong to psLib" }, + { GDIFF_ERR_${ErrorCode}, "${ErrorDescription}"}, + }; + static int nerror = GDIFF_ERR_NERROR - GDIFF_ERR_BASE; ///< number of values in enum + + for (int i = 0; i < nerror; i++) { + psErrorDescription *tmp = psAlloc(sizeof(psErrorDescription)); + *tmp = errors[i]; + psErrorRegister(tmp, 1); + psFree(tmp); /* it's on the internal list */ + } + nerror = 0; // don't register more than once +} Index: /branches/eam_branches/gDiff.20200202/src/gDiffErrorCodes.h.in =================================================================== --- /branches/eam_branches/gDiff.20200202/src/gDiffErrorCodes.h.in (revision 41260) +++ /branches/eam_branches/gDiff.20200202/src/gDiffErrorCodes.h.in (revision 41260) @@ -0,0 +1,28 @@ +/** @file gDiffErrorCodes.h.in + * @brief + * @ingroup gDiff + * + * @author IfA + * @version $Revision: 1.2 $ $Name: not supported by cvs2svn $ + * @date $Date: 2009-02-05 20:44:04 $ + * Copyright 2009 Institute for Astronomy, University of Hawaii + */ + +#if !defined(GDIFF_ERROR_CODES_H) +#define GDIFF_ERROR_CODES_H +/* + * The line + * GDIFF_ERR_$X{ErrorCode}, + * (without the X) + * + * will be replaced by values from errorCodes.dat + */ +typedef enum { + GDIFF_ERR_BASE = 14000, + GDIFF_ERR_${ErrorCode}, + GDIFF_ERR_NERROR +} gDiffErrorCode; + +void gDiffErrorRegister(void); + +#endif Index: /branches/eam_branches/gDiff.20200202/src/gDiffExit.c =================================================================== --- /branches/eam_branches/gDiff.20200202/src/gDiffExit.c (revision 41260) +++ /branches/eam_branches/gDiff.20200202/src/gDiffExit.c (revision 41260) @@ -0,0 +1,72 @@ +#include +#include +#include + +#include "gDiff.h" + +psExit gDiffExitCode(psExit exitValue) +{ + if (exitValue != PS_EXIT_SUCCESS) { + return exitValue; + } + + // gcc -Wswitch complains here if err is declared as type psErrorCode + // the collection of ps*ErrorCode values are enums defined separately for + // each module (psphot, pswarp, etc). the lowest type, psErrorCode is only the base set and does + // not include the possible psphot values + + // for now, to get around this, we just use an int for the switch + + // psErrorCode errorCode = psErrorCodeLast(); // Error code + int errorCode = psErrorCodeLast(); // Error code + if (errorCode != PS_ERR_NONE) { + pmFPAfileFreeSetStrict(false); + psErrorStackPrint(stderr, "Error in subtraction:"); + switch (errorCode) { + case GDIFF_ERR_UNKNOWN: + case PS_ERR_UNKNOWN: + psLogMsg("gDiff", PS_LOG_WARN, "Unknown error code: %x", errorCode); + exitValue = PS_EXIT_UNKNOWN_ERROR; + break; + case PS_ERR_IO: + case PS_ERR_DB_CLIENT: + case PS_ERR_DB_SERVER: + case PS_ERR_BAD_FITS: + case PS_ERR_OS_CALL_FAILED: + case PM_ERR_SYS: + case GDIFF_ERR_IO: + psLogMsg("gDiff", PS_LOG_WARN, "I/O error code: %x", errorCode); + exitValue = PS_EXIT_SYS_ERROR; + break; + case PS_ERR_BAD_PARAMETER_VALUE: + case PS_ERR_BAD_PARAMETER_TYPE: + case PS_ERR_BAD_PARAMETER_NULL: + case PS_ERR_BAD_PARAMETER_SIZE: + case GDIFF_ERR_ARGUMENTS: + case GDIFF_ERR_CONFIG: + psLogMsg("gDiff", PS_LOG_WARN, "Configuration error code: %x", errorCode); + exitValue = PS_EXIT_CONFIG_ERROR; + break; + case PSPHOT_ERR_PSF: + case GDIFF_ERR_DATA: + case GDIFF_ERR_NO_OVERLAP: + psLogMsg("gDiff", PS_LOG_WARN, "Data error code: %x", errorCode); + exitValue = PS_EXIT_DATA_ERROR; + break; + case PS_ERR_UNEXPECTED_NULL: + case PS_ERR_PROGRAMMING: + case GDIFF_ERR_NOT_IMPLEMENTED: + case GDIFF_ERR_PROG: + psLogMsg("gDiff", PS_LOG_WARN, "Programming error code: %x", errorCode); + exitValue = PS_EXIT_PROG_ERROR; + break; + default: + // It's a programming error if we're not dealing with the error correctly + psLogMsg("gDiff", PS_LOG_WARN, "Unrecognised error code: %x", errorCode); + exitValue = PS_EXIT_PROG_ERROR; + break; + } + } + + return exitValue; +} Index: /branches/eam_branches/gDiff.20200202/src/gDiffFiles.c =================================================================== --- /branches/eam_branches/gDiff.20200202/src/gDiffFiles.c (revision 41260) +++ /branches/eam_branches/gDiff.20200202/src/gDiffFiles.c (revision 41260) @@ -0,0 +1,235 @@ +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include +#include +#include + +#include "gDiff.h" + + +// 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 }; + +// 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", + NULL }; + +// Subtraction files +static const char *subFiles[] = { "GDIFF.OUTPUT", "GDIFF.OUTPUT.MASK", "GDIFF.OUTPUT.VARIANCE", + "GDIFF.OUTPUT.JPEG1", "GDIFF.OUTPUT.JPEG2", "GDIFF.OUTPUT.RESID.JPEG", + NULL }; + +// Subtraction photometry +static const char *subPhotFiles[] = { "GDIFF.OUTPUT.SOURCES", NULL }; + +// Inverse subtraction files +static const char *invFiles[] = { "GDIFF.INVERSE", "GDIFF.INVERSE.MASK", "GDIFF.INVERSE.VARIANCE", NULL }; + +// Inverse subtraction photometry +static const char *invPhotFiles[] = { "GDIFF.INVERSE.SOURCES", NULL }; + +// PSF files +static const char *psfFiles[] = { "PSPHOT.PSF.SAVE", NULL }; + +// Calculation (may be either input or output) files +static const char *calcFiles[] = { "GDIFF.OUTPUT.KERNELS", NULL }; + + +// Activate/deactivate a list of files +static void filesActivate(pmConfig *config, // Configuration + const char **files, // List of files + bool state // Activation status to set + ) +{ + for (int i = 0; files[i]; i++) { + pmFPAfileActivate(config->files, state, files[i]); + } + return; +} + +// Activate/deactivate a list of files depending on their 'save' boolean. +// This is so we can activate/deactivate the 'calculation' files, which may be either input or output, which +// is indicated by their 'save' boolean. +static void filesActivateSave(pmConfig *config, // Configuration + const char **files, // List of files + bool save, // Activate when this save state is set + bool state // Activation status to set + ) +{ + for (int i = 0; files[i]; i++) { + pmFPAfile *file = pmFPAfileSelectSingle(config->files, files[i], 0); + if (file && file->save == save) { + pmFPAfileActivate(config->files, state, files[i]); + } + } + return; +} + +void gDiffFilesActivate(pmConfig *config, gDiffFiles files, bool state) +{ + psAssert(config, "Require configuration"); + + if (files & GDIFF_FILES_INPUT) { + filesActivate(config, inputFiles, state); + filesActivateSave(config, calcFiles, false, state); + } + if (files & GDIFF_FILES_CONV) { + filesActivate(config, convFiles, state); + } + if (files & GDIFF_FILES_SUB) { + filesActivate(config, subFiles, state); + filesActivateSave(config, calcFiles, true, state); + } + if (files & GDIFF_FILES_INV) { + filesActivate(config, invFiles, state); + } + if (files & GDIFF_FILES_PHOT_SUB) { + filesActivate(config, subPhotFiles, state); + } + if (files & GDIFF_FILES_PHOT_INV) { + filesActivate(config, invPhotFiles, state); + } + if (files & GDIFF_FILES_PSF) { + filesActivate(config, psfFiles, state); + } + if (files & GDIFF_FILES_PHOT) { + psphotFilesActivate(config, state); + } + + return; +} + + +pmFPAview *gDiffViewReadout(void) +{ + pmFPAview *view = pmFPAviewAlloc(0); + view->chip = view->cell = view->readout = 0; + return view; +} + +bool gDiffFilesIterateDown(pmConfig *config, gDiffFiles files) +{ + psAssert(config, "Require configuration"); + + gDiffFilesActivate(config, GDIFF_FILES_ALL, false); + gDiffFilesActivate(config, files, true); + + pmFPAview *view = pmFPAviewAlloc(0);// View to FPA top + + // FPA + if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) { + psFree (view); + return false; + } + + // Chip + view->chip = 0; + if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) { + psFree (view); + return false; + } + + // Cell + view->cell = 0; + if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) { + psFree (view); + return false; + } + + // Readout + view->readout = 0; + if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) { + psFree (view); + return false; + } + + gDiffFilesActivate(config, GDIFF_FILES_ALL, false); + + psFree (view); + return true; +} + +bool gDiffFilesIterateUp(pmConfig *config, gDiffFiles files) +{ + psAssert(config, "Require configuration"); + + gDiffFilesActivate(config, GDIFF_FILES_ALL, false); + gDiffFilesActivate(config, files, true); + + pmFPAview *view = gDiffViewReadout(); // View to readout + + // Readout + if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) { + psFree (view); + return false; + } + + // Cell + view->readout = -1; + if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) { + psFree (view); + return false; + } + + // Chip + view->cell = -1; + if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) { + psFree (view); + return false; + } + + // FPA + view->chip = -1; + if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) { + psFree (view); + return false; + } + + psFree(view); + + // Disable writing of data now that we've written things out + if (files & GDIFF_FILES_CONV) { + pmFPAview *view = gDiffViewReadout(); // View to readout + { + pmReadout *ro = pmFPAfileThisReadout(config->files, view, "GDIFF.INPUT.CONV"); + if (ro) { + ro->data_exists = ro->parent->data_exists = ro->parent->parent->data_exists = false; + } + } + { + pmReadout *ro = pmFPAfileThisReadout(config->files, view, "GDIFF.REF.CONV"); + if (ro) { + ro->data_exists = ro->parent->data_exists = ro->parent->parent->data_exists = false; + } + } + psFree(view); + } + if (files & GDIFF_FILES_SUB) { + pmFPAview *view = gDiffViewReadout(); // View to readout + pmReadout *ro = pmFPAfileThisReadout(config->files, view, "GDIFF.OUTPUT"); + if (ro) { + ro->data_exists = ro->parent->data_exists = ro->parent->parent->data_exists = false; + } + psFree(view); + } + // Disable writing of data now that we've written things out + if (files & GDIFF_FILES_INV) { + pmFPAview *view = gDiffViewReadout(); // View to readout + pmReadout *ro = pmFPAfileThisReadout(config->files, view, "GDIFF.INVERSE"); + if (ro) { + ro->data_exists = ro->parent->data_exists = ro->parent->parent->data_exists = false; + } + psFree(view); + } + + gDiffFilesActivate(config, GDIFF_FILES_ALL, false); + + psFree (view); + return true; +} Index: /branches/eam_branches/gDiff.20200202/src/gDiffLoop.c =================================================================== --- /branches/eam_branches/gDiff.20200202/src/gDiffLoop.c (revision 41260) +++ /branches/eam_branches/gDiff.20200202/src/gDiffLoop.c (revision 41260) @@ -0,0 +1,336 @@ +#include "gDiff.h" + +bool dumpout(pmConfig *config, char *name) +{ + pmFPAview *view = gDiffViewReadout(); // View to readout + pmReadout *out = pmFPAfileThisReadout(config->files, view, "GDIFF.OUTPUT"); + psphotSaveImage (NULL, out->image, name); + psFree(view); + return true; +} + + +bool gDiffLoop(gDiffData *data) +{ + 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; +} Index: /branches/eam_branches/gDiff.20200202/src/gDiffVersion.c =================================================================== --- /branches/eam_branches/gDiff.20200202/src/gDiffVersion.c (revision 41260) +++ /branches/eam_branches/gDiff.20200202/src/gDiffVersion.c (revision 41260) @@ -0,0 +1,102 @@ +#include "gDiff.h" +#include +#include "gDiffVersionDefinitions.h" + +#ifndef GDIFF_VERSION +#error "GDIFF_VERSION is not set" +#endif +#ifndef GDIFF_BRANCH +#error "GDIFF_BRANCH is not set" +#endif +#ifndef GDIFF_SOURCE +#error "GDIFF_SOURCE is not set" +#endif + +psString gDiffVersion(void) { + char *value = NULL; + psStringAppend(&value, "%s@%s", GDIFF_BRANCH, GDIFF_VERSION); + return value; +} + +psString gDiffSource(void) { + return psStringCopy(GDIFF_SOURCE); +} + +psString gDiffVersionLong(void) { + psString version = gDiffVersion(); // Version, to return + psString source = gDiffSource(); // Source + + psStringPrepend(&version, "gDiff "); + psStringAppend(&version, " from %s, built %s, %s", source, __DATE__, __TIME__); + psFree(source); + +#ifdef __OPTIMIZE__ + psStringAppend(&version, " optimised"); +#else + psStringAppend(&version, " unoptimised"); +#endif + + return version; +}; + +bool gDiffVersionHeader(psMetadata *header) { + PS_ASSERT_METADATA_NON_NULL(header, false); + + psTime *time = psTimeGetNow(PS_TIME_TAI); // The time now + psString timeString = psTimeToISO(time); // The time in an ISO string + psFree(time); + psString history = NULL; // History string + psStringAppend(&history, "gDiff at %s", timeString); + psFree(timeString); + psMetadataAddStr(header, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK, NULL, history); + psFree(history); + + psLibVersionHeader(header); + psModulesVersionHeader(header); + psphotVersionHeader(header); + ppStatsVersionHeader(header); + + psString version = gDiffVersion(); // Software version + psString source = gDiffSource(); // Software source + + psMetadataAddStr(header, PS_LIST_TAIL, "GDIFF_V", PS_META_REPLACE, NULL, GDIFF_VERSION); + + psStringPrepend(&version, "gDiff version: "); + psStringPrepend(&source, "gDiff source: "); + + psMetadataAddStr(header, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK, NULL, version); + psMetadataAddStr(header, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK, NULL, source); + + psFree(version); + psFree(source); + + return true; +} + +void gDiffVersionPrint(void) { + psTime *time = psTimeGetNow(PS_TIME_TAI); // The time now + psString timeString = psTimeToISO(time); // The time in an ISO string + psFree(time); + psLogMsg("gDiff", PS_LOG_INFO, "gDiff at %s", timeString); + psFree(timeString); + + psString pslib = psLibVersionLong();// psLib version + psString psmodules = psModulesVersionLong(); // psModules version + psString psphot = psphotVersionLong(); // psphot version + psString ppStats = ppStatsVersionLong(); // ppStats version + psString gDiff = gDiffVersionLong(); // gDiff version + + psLogMsg("gDiff", PS_LOG_INFO, "%s", pslib); + psLogMsg("gDiff", PS_LOG_INFO, "%s", psmodules); + psLogMsg("gDiff", PS_LOG_INFO, "%s", psphot); + psLogMsg("gDiff", PS_LOG_INFO, "%s", ppStats); + psLogMsg("gDiff", PS_LOG_INFO, "%s", gDiff); + + psFree(pslib); + psFree(psmodules); + psFree(psphot); + psFree(ppStats); + psFree(gDiff); + + return; +} Index: /branches/eam_branches/gDiff.20200202/src/gDiffVersionDefinitions.h.in =================================================================== --- /branches/eam_branches/gDiff.20200202/src/gDiffVersionDefinitions.h.in (revision 41260) +++ /branches/eam_branches/gDiff.20200202/src/gDiffVersionDefinitions.h.in (revision 41260) @@ -0,0 +1,8 @@ +#ifndef GDIFF_VERSION_DEFINITIONS_H +#define GDIFF_VERSION_DEFINITIONS_H + +#define GDIFF_VERSION @GDIFF_VERSION@ // SVN version +#define GDIFF_BRANCH @GDIFF_BRANCH@ // SVN branch +#define GDIFF_SOURCE @GDIFF_SOURCE@ // SVN source + +#endif