Changeset 35563

Timestamp:

May 9, 2013, 12:26:48 PM (13 years ago)

Author:

eugene

Message:

major upgrades to pswarp to enable skycell warp -> chip

Location:

trunk/pswarp

Files:

• doc/notes.20130406.txt (copied) (copied from branches/eam_branches/ipp-20130419/pswarp/doc/notes.20130406.txt )
• src/Makefile.am (modified) (2 diffs)
• src/pswarp.c (modified) (3 diffs)
• src/pswarp.h (modified) (4 diffs)
• src/pswarpArguments.c (modified) (5 diffs)
• src/pswarpCleanup.c (copied) (copied from branches/eam_branches/ipp-20130419/pswarp/src/pswarpCleanup.c )
• src/pswarpDefine.c (modified) (4 diffs)
• src/pswarpDefineBackground.c (copied) (copied from branches/eam_branches/ipp-20130419/pswarp/src/pswarpDefineBackground.c )
• src/pswarpDefineLayout.c (copied) (copied from branches/eam_branches/ipp-20130419/pswarp/src/pswarpDefineLayout.c )
• src/pswarpDefineSkycell.c (modified) (2 diffs)
• src/pswarpExit.c (deleted)
• src/pswarpExitCode.c (copied) (copied from branches/eam_branches/ipp-20130419/pswarp/src/pswarpExitCode.c )
• src/pswarpLoadAstrometry.c (copied) (copied from branches/eam_branches/ipp-20130419/pswarp/src/pswarpLoadAstrometry.c )
• src/pswarpLoop.c (modified) (3 diffs)
• src/pswarpLoopBackground.c (copied) (copied from branches/eam_branches/ipp-20130419/pswarp/src/pswarpLoopBackground.c )
• src/pswarpLoopSkycell.c (copied) (copied from branches/eam_branches/ipp-20130419/pswarp/src/pswarpLoopSkycell.c )
• src/pswarpMakePSF.c (copied) (copied from branches/eam_branches/ipp-20130419/pswarp/src/pswarpMakePSF.c )
• src/pswarpMaskStats.c (modified) (2 diffs)
• src/pswarpMatchRange.c (modified) (1 diff)
• src/pswarpOptions.c (copied) (copied from branches/eam_branches/ipp-20130419/pswarp/src/pswarpOptions.c )
• src/pswarpOverlaps.c (copied) (copied from branches/eam_branches/ipp-20130419/pswarp/src/pswarpOverlaps.c )
• src/pswarpParseCamera.c (modified) (9 diffs)
• src/pswarpPixelsLit.c (modified) (1 diff)
• src/pswarpStatsFile.c (copied) (copied from branches/eam_branches/ipp-20130419/pswarp/src/pswarpStatsFile.c )
• src/pswarpTransformReadout.c (modified) (6 diffs)
• src/pswarpTransformTile.c (modified) (3 diffs)
• src/pswarpUpdateMetadata.c (copied) (copied from branches/eam_branches/ipp-20130419/pswarp/src/pswarpUpdateMetadata.c )
• src/pswarpUpdateStatistics.c (copied) (copied from branches/eam_branches/ipp-20130419/pswarp/src/pswarpUpdateStatistics.c )
• src/pswarpVersion.c (modified) (1 diff)

trunk/pswarp/src/Makefile.am

@@ -28 +28 @@
         pswarp.c                        \
         pswarpArguments.c               \
-        pswarpExit.c                    \
+        pswarpExitCode.c                \
         pswarpLoop.c                    \
+        pswarpLoopBackground.c          \
+        pswarpCleanup.c                 \
+        pswarpStatsFile.c               \
         pswarpDefine.c                  \
-        pswarpDefineSkycell.c           \
+        pswarpDefineBackground.c        \
+        pswarpDefineSkycell.c           \
+        pswarpDefineLayout.c            \
         pswarpErrorCodes.c              \
+        pswarpLoadAstrometry.c          \
+        pswarpMakePSF.c                 \
         pswarpMapGrid.c                 \
         pswarpMaskStats.c               \
         pswarpMatchRange.c              \
+        pswarpOverlaps.c                \
+        pswarpOptions.c                 \
         pswarpParseCamera.c             \
         pswarpPixelsLit.c               \
@@ -43 +52 @@
         pswarpTransformSources.c        \
         pswarpTransformTile.c           \
+        pswarpUpdateStatistics.c        \
+        pswarpUpdateMetadata.c          \
         pswarpVersion.c                 \
         pswarpFiles.c

trunk/pswarp/src/pswarp.c

@@ -12 +12 @@
 
 #include "pswarp.h"
-#include "pswarpFileNames.h"
 
 int main (int argc, char **argv)
@@ -22 +21 @@
     psphotInit();
 
-    const char *statsName = NULL;       // Filename for statistics
-    psMetadata *stats = NULL;           // Container for statistics
-    FILE *statsFile = NULL;             // File stream for statistics
-
     pmConfig *config = pswarpArguments(argc, argv);
     if (!config) {
-        goto DIE;
+        pswarpCleanup(config, NULL);
     }
+
+    pswarpStatsFile *statsFile = pswarpStatsFileOpen (config);
 
     pswarpVersionPrint();
@@ -35 +32 @@
     // load identify the data sources
     if (!pswarpParseCamera(config)) {
-        goto DIE;
+        pswarpCleanup(config, statsFile);
     }
 
     if (!pswarpOptions(config)) {
-        goto DIE;
+        pswarpCleanup(config, statsFile);
     }
 
-    // load the skycell layout information
-    if (!pswarpDefine(config)) {
-        goto DIE;
-    }
-    if (psMetadataLookupBool(NULL,config->arguments,"BACKGROUND.MODEL")) {
-      if (!pswarpDefineBackground(config)) {
-        goto DIE;
-      }
+    // load the input & output astrometry, find the output overlaps, generate the output pixels
+    if (!pswarpDefineLayout(config)) {
+        pswarpCleanup(config, statsFile);
     }
 
-    // Open the statistics file
-    bool mdok;                          ///< Status of MD lookup
-    statsName = psMetadataLookupStr(&mdok, config->arguments, "STATS"); ///< Filename for statistics
-    if (mdok && statsName && strlen(statsName) > 0) {
-        psString resolved = pmConfigConvertFilename(statsName, config, true, true);
-        statsFile = fopen(resolved, "w");
-        if (!statsFile) {
-            psError(PS_ERR_IO, true, "Unable to open statistics file %s for writing.\n", resolved);
-            psFree(resolved);
-            goto DIE;
-        }
-        psFree(resolved);
-        stats = psMetadataAlloc();
-        psMetadataAddS32(stats, PS_LIST_TAIL, "QUALITY", 0, "No problems", 0);
+    // load input pixels & and warp
+    if (!pswarpLoop(config, statsFile->md)) {
+        pswarpCleanup(config, statsFile);
     }
 
-    // load and warp
-    if (!pswarpLoop(config, stats)) {
-        goto DIE;
-    }
-    if (psMetadataLookupBool(NULL,config->arguments,"BACKGROUND.MODEL")) {
-      if (!pswarpLoopBackground(config, stats)) {
-        fprintf(stderr,"Dying!\n");
-        goto DIE;
-      }
+    // load input pixels & and warp
+    if (!pswarpLoopBackground(config, statsFile->md)) {
+        pswarpCleanup(config, statsFile);
     }
 
+    // output and free 
+    // NOTE: pswarpCleanup calls exit
     psLogMsg("pswarp", PS_LOG_INFO, "complete pswarp run: %f sec\n", psTimerMark("pswarp"));
-
-DIE:
-    {
-        psExit exitValue = pswarpExitCode(PS_EXIT_SUCCESS); // Exit code
-
-        // Ensure everything is written out, at every level
-        pswarpFileActivation(config, detectorFiles, true);
-        pswarpFileActivation(config, skycellFiles, true);
-        pswarpFileActivation(config, photFiles, true);
-        pswarpFileActivation(config, independentFiles, true);
-        if (!pswarpIOChecksAfter(config)) {
-            psError(psErrorCodeLast(), false, "Unable to write files.");
-            exitValue = pswarpExitCode(exitValue);
-            pmFPAfileFreeSetStrict(false);
-        }
-
-        // Write out summary statistics
-        if (stats && statsFile) {
-            psMetadataAddF32(stats, PS_LIST_TAIL, "DT_WARP", 0, "Time for warp completion",
-                             psTimerMark("pswarp"));
-
-            const char *statsMDC = psMetadataConfigFormat(stats);
-            if (!statsMDC) {
-                psError(psErrorCodeLast(), false, "Unable to get statistics file.");
-                exitValue = pswarpExitCode(exitValue);
-            }
-            psFree(stats);
-            if (fprintf(statsFile, "%s", statsMDC) != strlen(statsMDC)) {
-                psError(PSWARP_ERR_IO, true, "Unable to write statistics file.");
-                exitValue = pswarpExitCode(exitValue);
-            }
-            psFree(statsMDC);
-            if (fclose(statsFile) == EOF) {
-                psError(PSWARP_ERR_IO, true, "Unable to close statistics file.");
-                exitValue = pswarpExitCode(exitValue);
-            }
-            pmConfigRunFilenameAddWrite(config, "STATS", statsName);
-        }
-
-        // Dump configuration
-        psString dump_file = psMetadataLookupStr(&mdok, config->arguments, "DUMP_CONFIG");
-        if (dump_file) {
-            if (!pmConfigDump(config, dump_file)) {
-                psError(psErrorCodeLast(), false, "Unable to dump configuration");
-                exitValue = pswarpExitCode(exitValue);
-            }
-        }
-
-        psThreadPoolFinalize();
-        psMemCheckCorruption(stderr, true);
-
-        psFree(config);
-
-        psTimerStop();
-        pmVisualClose();
-        pmModelClassCleanup();
-        pmConceptsDone();
-        pmConfigDone();
-        psLibFinalize();
-
-        exitValue = pswarpExitCode(exitValue);
-        return exitValue;
-    }
+    pswarpCleanup(config, statsFile);
 }

trunk/pswarp/src/pswarp.h

@@ -22 +22 @@
 #include <psmodules.h>
 #include <psphot.h>
+#include <ppStats.h>
 
 #define THREADED 1
@@ -78 +79 @@
 } pswarpTransformTileArgs;
 
+typedef struct {
+    FILE *f;                            // File stream for statistics
+    char *name;                         // Filename for statistics
+    psMetadata *md;                     // Container for statistics
+} pswarpStatsFile;
+
 pswarpTransformTileArgs *pswarpTransformTileArgsAlloc(void);
 bool pswarpTransformTile (pswarpTransformTileArgs *args);
@@ -88 +95 @@
 bool pswarpDefineBackground (pmConfig *config);
 bool pswarpLoop (pmConfig *config, psMetadata *stats);
+bool pswarpLoopSkycell (pmConfig *config, psMetadata *stats);
 bool pswarpLoopBackground (pmConfig *config, psMetadata *stats);
 psExit pswarpExitCode(psExit exitValue);
-bool pswarpTransformReadout (pmReadout *output, pmReadout *input, pmConfig *config);
+bool pswarpTransformReadout (pmReadout *output, pmReadout *input, pmConfig *config, bool backgroundWarp);
 bool pswarpTransformSources(pmReadout *output, pmReadout *input, pmConfig *config);
 
@@ -160 +168 @@
     );
 
+pswarpStatsFile *pswarpStatsFileAlloc ();
+pswarpStatsFile *pswarpStatsFileOpen (pmConfig *config);
+bool pswarpStatsFileSave (pmConfig *config, pswarpStatsFile *statsFile);
+
+// cleanup memory and exit
+void pswarpCleanup (pmConfig *config, pswarpStatsFile *statsFile) PS_ATTR_NORETURN;
+
+// load the astrometry header info and generate the tranformations 
+bool pswarpDefineLayout (pmConfig *config);
+
+// XXX function for testing
+bool pswarpDumpOutput (pmConfig *config);
+
+// structure to describe approximate bounds for each fpa element (eg, chip)
+// P,Q are a locally linear projection
+typedef struct {
+    psVector *Pmin; 
+    psVector *Pmax;
+    psVector *Qmin;
+    psVector *Qmax;
+} pswarpBounds;
+
+bool pswarpFindOverlap (pmFPA *input, pmFPA *output, pswarpBounds *src, pswarpBounds *tgt);
+psProjection *pswarpLocalFrame (pmFPA *fpa);
+pswarpBounds *pswarpMakeBounds (pmFPA *fpa, psProjection *frame);
+bool pswarpBoundsAppend(pswarpBounds *bounds, float Pmin, float Pmax, float Qmin, float Qmax);
+pswarpBounds *pswarpBoundsAlloc();
+
+bool pswarpLoadAstrometry (pmFPAfile *target, pmFPAfile *astrom, pmConfig *config);
+
+bool pswarpTransformToTarget (pmFPA *output, pmReadout *input, pmConfig *config, bool backgroundWarp);
+bool pswarpMakePSF (pmConfig *config, pmFPAfile *output, psMetadata *stats);
+bool pswarpUpdateStatistics (pmFPA *output, psMetadata *stats, pmFPA *input, pmFPA *astrom, pmConfig *config);
+bool pswarpUpdateMetadata (pmFPA *output, pmFPA *skycell, pmFPA *input, pmFPA *astrom, pmConfig *config, bool fullImage);
+
+// XXX functions in pswarpParseCamera
+
+bool AddStringAsArray (psMetadata *md, char *string, char *name);
+
+pmFPAfile *pswarpDefineInputFile(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 pswarpParseSingleInput (pmConfig *config);
+bool pswarpParseMultiInput (pmConfig *config, psMetadata *fileListMD);
+
+bool pswarpModifyChipAstrom (pmConfig *config, pmFPAview *view, pmChip *chip, pmFPAfile *astrom, bool bilevelAstrometry, double xBin, double yBin);
+bool pswarpGetInputScales (double *xBin, double *yBin, pmConfig *config, pmFPAview *view, pmChip *chip);

trunk/pswarp/src/pswarpArguments.c

@@ -12 +12 @@
 
 # include "pswarp.h"
-# include <glob.h>
-
 
 static void usage (void) {
-    fprintf(stderr, "USAGE: pswarp [-file image(s)] [-list imagelist] [options] (output) (skycell)\n");
+    fprintf(stderr, "USAGE: pswarp [-input input.mdc] [-file image(s)] [-list imagelist] [options] (output) (skycell)\n");
     fprintf(stderr, "  options:\n");
+    fprintf(stderr, "    [-input input.mdc] : input image information in a metadata file\n");
+    fprintf(stderr, "    [-file input.fits[,input.fits]] : input image to be warped\n");
     fprintf(stderr, "    [-astrom astrom.cmp] : provide an alternative astrometry calibration\n");
     fprintf(stderr, "    [-mask mask.fits] : provide a corresponding mask image\n");
@@ -24 +24 @@
     exit(PS_EXIT_CONFIG_ERROR);
 }
-
 
 pmConfig *pswarpArguments (int argc, char **argv) {
@@ -55 +54 @@
     }
 
+    // XXX move to the single group below?
     pmConfigFileSetsMD(config->arguments, &argc, argv, "ASTROM",   "-astrom", "-astromlist");
 
-    {
-        int arg;                        ///< Argument Number
-        if ((arg = psArgumentGet(argc, argv, "-psphot-visual"))) {
-            psArgumentRemove(arg, &argc, argv);
-            pmVisualSetVisual(true);
-        }
+    // turn on psphot visualization
+    if ((N = psArgumentGet(argc, argv, "-psphot-visual"))) {
+        psArgumentRemove(N, &argc, argv);
+        pmVisualSetVisual(true);
     }
 
@@ -83 +81 @@
     pswarpSetThreads();
 
-    pmConfigFileSetsMD (config->arguments, &argc, argv, "INPUT", "-file", "-list");
-    pmConfigFileSetsMD (config->arguments, &argc, argv, "MASK", "-mask", "-masklist");
-    pmConfigFileSetsMD (config->arguments, &argc, argv, "VARIANCE", "-variance", "-variancelist");
-    pmConfigFileSetsMD (config->arguments, &argc, argv, "BACKGROUND", "-background", "-bkglist");
+    // there are three mutually exclusive ways of providing the input
+    // 1) supply -file (filename) [-mask .. -variance ..] on the command line
+    // 2) supply -input (input.mdc) on the command line
+    // 3) load inputs from RUN config info
+
+    // below, we check first for -file then for -input.  failure to find either implies use of
+    // the configuration metadata file.
+
+    bool singleInput = pmConfigFileSetsMD (config->arguments, &argc, argv, "INPUT", "-file", "-list");
+    if (singleInput) {
+        if (psArgumentGet(argc, argv, "-input")) {
+            psErrorStackPrint(stderr, "error in arguments : -input and -file / -list are mutually exclusive");
+            exit(PS_EXIT_CONFIG_ERROR);
+        }       
+        pmConfigFileSetsMD (config->arguments, &argc, argv, "MASK", "-mask", "-masklist");
+        pmConfigFileSetsMD (config->arguments, &argc, argv, "VARIANCE", "-variance", "-variancelist");
+        pmConfigFileSetsMD (config->arguments, &argc, argv, "BACKGROUND", "-background", "-bkglist");
+    } else {
+        // find the input data file (an mdc file)
+        if ((N = psArgumentGet(argc, argv, "-input"))) {
+            if (argc <= N+1) {
+                psErrorStackPrint(stderr, "Expected to see 1 more argument; saw %d", argc - 1);
+                exit(PS_EXIT_CONFIG_ERROR);
+            }
+            psArgumentRemove(N, &argc, argv);
+
+            unsigned int numBad = 0;                     // Number of bad lines
+            psMetadata *inputs = psMetadataConfigRead(NULL, &numBad, argv[N], false); // Input file info
+            if (!inputs || numBad > 0) {
+                psErrorStackPrint(stderr, "Unable to cleanly read MDC file with inputs.");
+                exit(PS_EXIT_CONFIG_ERROR);
+            }
+            psMetadataAddMetadata(config->arguments, PS_LIST_TAIL, "INPUTS", 0, "Metadata with input details", inputs);
+            psFree(inputs);
+
+            psArgumentRemove(N, &argc, argv);
+        }
+    }
+    if (argc != 3) {
+        usage();
+    }
+    if (psErrorCodeLast() != PS_ERR_NONE) {
+        psErrorStackPrint(stderr, "error in arguments");
+        exit(PS_EXIT_CONFIG_ERROR);
+    }
     
-    if (argc != 3) {
-        usage();
-    }
-
-
     psArray *array;
 
@@ -107 +141 @@
     return config;
 }
-
-/**
- * Parse the recipe and format into the arguments
- */
-bool pswarpOptions(pmConfig *config)
-{
-    // Select the appropriate recipe
-    psMetadata *recipe  = psMetadataLookupPtr(NULL, config->recipes, PSWARP_RECIPE);
-    if (!recipe) {
-        psError(PSWARP_ERR_CONFIG, true, "Can't find %s recipe!\n", PSWARP_RECIPE);
-        return false;
-    }
-
-    // Get grid size
-    bool status;                        ///< Status of MD lookup
-    int nGridX = psMetadataLookupS32(&status, recipe, "GRID.NX");
-    if (!status || nGridX <= 0) {
-        nGridX = 128;
-        psWarning("GRID.NX is not set in the recipe --- defaulting to %d", nGridX);
-    }
-    int nGridY = psMetadataLookupS32(&status, recipe, "GRID.NY");
-    if (!status) {
-        nGridY = 128;
-        psWarning("GRID.NY is not set in the recipe --- defaulting to %d", nGridY);
-    }
-
-    // Get interpolation mode
-    const char *name = psMetadataLookupStr (&status, recipe, "INTERPOLATION.MODE"); ///< Name of interp mode
-    if (!name) {
-        name = "BILINEAR";
-        psLogMsg("pswarp", 3, "defaulting to %s interpolation", name);
-    }
-    psImageInterpolateMode interpolationMode = psImageInterpolateModeFromString(name); ///< Mode for interp.
-    if (interpolationMode == PS_INTERPOLATE_NONE) {
-        interpolationMode = PS_INTERPOLATE_BILINEAR;
-        psLogMsg ("pswarp", 3,
-                  "Unknown interpolation mode %s, defaulting to bilinear interpolation\n", name);
-        name = "BILINEAR";
-    }
-
-    int numKernels = psMetadataLookupS32(&status, recipe, "INTERPOLATION.NUM");
-    if (!status) {
-        numKernels = 0;
-        psWarning("INTERPOLATION.NUM is not set in the recipe --- defaulting to %d", numKernels);
-    }
-
-    float poorFrac = psMetadataLookupF32(&status, recipe, "POOR.FRAC"); ///< Frac of bad flux for a "poor"
-    if (!status) {
-        poorFrac = 0.0;
-        psWarning("POOR.FRAC is not set in the %s recipe --- defaulting to %f.", PSWARP_RECIPE, poorFrac);
-    }
-
-    bool PSF = psMetadataLookupBool(&status, recipe, "PSF"); ///< Generate a PSF model?
-    if (!status) {
-        PSF = true;
-        psWarning("PSF is not set in the %s recipe --- defaulting to TRUE.", PSWARP_RECIPE);
-    }
-
-    bool doBKG = psMetadataLookupBool(&status,recipe, "BACKGROUND.MODEL"); ///< Generate the warped background model?
-    if (!status) {
-      doBKG = false;
-      psWarning("BACKGROUND.MODEL is not set in the %s recipe -- defaulting to FALSE.", PSWARP_RECIPE);
-    }
-    int bkgXgrid = psMetadataLookupS32(&status,recipe, "BKG.XGRID"); ///< Xsize of background model
-    if (!status) {
-      bkgXgrid = 10;
-      psWarning("BKG.XGRID is not set in the %s recipe -- defaultint to %d.",PSWARP_RECIPE,bkgXgrid);
-    }
-    int bkgYgrid = psMetadataLookupS32(&status,recipe, "BKG.YGRID"); ///< Xsize of background model
-    if (!status) {
-      bkgYgrid = 10;
-      psWarning("BKG.YGRID is not set in the %s recipe -- defaultint to %d.",PSWARP_RECIPE,bkgYgrid);
-    }
-
-    
-    // Set recipe values in the recipe (since we've possibly altered some)
-    psMetadataAddS32(recipe, PS_LIST_TAIL, "GRID.NX", PS_META_REPLACE,
-                     "Iso-astrom grid spacing in x", nGridX);
-    psMetadataAddS32(recipe, PS_LIST_TAIL, "GRID.NY", PS_META_REPLACE,
-                     "Iso-astrom grid spacing in y", nGridY);
-    psMetadataAddStr(recipe, PS_LIST_TAIL, "INTERPOLATION.MODE", PS_META_REPLACE,
-                     "Interpolation mode", name);
-    psMetadataAddS32(recipe, PS_LIST_TAIL, "INTERPOLATION.NUM", PS_META_REPLACE,
-                     "Interpolation pre-calculated kernels", numKernels);
-    psMetadataAddF32(recipe, PS_LIST_TAIL, "POOR.FRAC", PS_META_REPLACE,
-                     "Fraction of bad flux for a pixel to be marked as poor", poorFrac);
-    psMetadataAddBool(recipe, PS_LIST_TAIL, "PSF", PS_META_REPLACE, "Generate a PSF Model?", PSF);
-    psMetadataAddBool(recipe, PS_LIST_TAIL, "BACKGROUND.MODEL", PS_META_REPLACE, "Generate the warped background model?", doBKG);
-    psMetadataAddS32(recipe, PS_LIST_TAIL, "BKG.XGRID", PS_META_REPLACE, "Xsize of background model", bkgXgrid);
-    psMetadataAddS32(recipe, PS_LIST_TAIL, "BKG.YGRID", PS_META_REPLACE, "Ysize of background model", bkgYgrid);
-    
-    
-    // Set recipe values in the arguments
-    psMetadataAddS32(config->arguments, PS_LIST_TAIL, "GRID.NX", 0,
-                     "Iso-astrom grid spacing in x", nGridX);
-    psMetadataAddS32(config->arguments, PS_LIST_TAIL, "GRID.NY", 0,
-                     "Iso-astrom grid spacing in y", nGridY);
-    psMetadataAddS32(config->arguments, PS_LIST_TAIL, "INTERPOLATION.MODE", 0,
-                     "Interpolation mode", interpolationMode);
-    psMetadataAddS32(config->arguments, PS_LIST_TAIL, "INTERPOLATION.NUM", 0,
-                     "Interpolation pre-calculated kernels", numKernels);
-    psMetadataAddF32(config->arguments, PS_LIST_TAIL, "POOR.FRAC", 0,
-                     "Fraction of bad flux for a pixel to be marked as poor", poorFrac);
-    psMetadataAddBool(config->arguments, PS_LIST_TAIL, "PSF", PS_META_REPLACE, "Generate a PSF Model?", PSF);
-    psMetadataAddBool(config->arguments, PS_LIST_TAIL, "BACKGROUND.MODEL", PS_META_REPLACE, "Generate the warped background model?", doBKG);
-    psMetadataAddS32(config->arguments, PS_LIST_TAIL, "BKG.XGRID", PS_META_REPLACE, "Xsize of background model", bkgXgrid);
-    psMetadataAddS32(config->arguments, PS_LIST_TAIL, "BKG.YGRID", PS_META_REPLACE, "Ysize of background model", bkgYgrid);
-
-    psTrace("pswarp", 1, "Done with pswarpArguments...\n");
-
-    return (config);
-}

trunk/pswarp/src/pswarpDefine.c

@@ -1 +1 @@
 /** @file pswarpDefine.c
  *
- *  @brief
- *
+ *  @brief generates the output images
  *  @ingroup pswarp
  *
@@ -13 +12 @@
 # include "pswarp.h"
 
-/**
- * loads the skycell metadata
- */
 bool pswarpDefine (pmConfig *config) {
 
@@ -42 +38 @@
     }
     
-    // open the full skycell file; no need to defer different depths. only load the header data
-    pmFPAview *view = pmFPAviewAlloc (0);
-    pmFPAfileOpen (skycell, view, config);
+    // Read the output astrometry
+    // XXX this section loops over the input skycell to load the header info
+    {
+        pmFPAview *view = pmFPAviewAlloc(0);
 
-    // Read header and create target
-    {
-        if (!pmFPAReadHeaderSet(skycell->fpa, skycell->fits, config)) {
-            psError(psErrorCodeLast(), false, "Unable to read headers for skycell.");
-            psFree(view);
-            return false;
-        }
-        view->chip = 0;
-        view->cell = 0;
-        view->readout = 0;
-        pmCell *source = pmFPAfileThisCell(config->files, view, "PSWARP.SKYCELL"); ///< Source cell
-        pmHDU *hdu = pmHDUFromCell(source); ///< HDU for source
-        if (!hdu || !hdu->header) {
-            psError(PM_ERR_PROG, false, "Unable to find header for sky cell.");
-            psFree(view);
-            return false;
-        }
-        int numCols = psMetadataLookupS32(NULL, hdu->header, "NAXIS1"); ///< Number of columns
-        int numRows = psMetadataLookupS32(NULL, hdu->header, "NAXIS2"); ///< Number of rows
-        if ((numCols == 0) || (numRows == 0)) {
-            psError(PSWARP_ERR_DATA, false, "skycell has invalid dimensions %d x %d", numCols, numRows);
-            psFree(view);
-            return false;
-        }
+        pmChip *chip;
+        if (!pmFPAfileRead (skycell, view, config)) {
+            psError(PS_ERR_IO, false, "failed READ at FPA %s", skycell->name);
+            psFree(view);
+            return false;
+        }
+        while ((chip = pmFPAviewNextChip (view, skycell->fpa, 1)) != NULL) {
+            psTrace ("pswarp", 4, "Chip %d: %x %x\n", view->chip, chip->file_exists, chip->process);
+            if (!chip->process || !chip->file_exists) { continue; }
+            if (!pmFPAfileRead (skycell, view, config)) {
+                psError(psErrorCodeLast(), false, "failed READ at CHIP %s", skycell->name);
+                psFree(view);
+                return false;
+            }
+            pmCell *cell;
+            while ((cell = pmFPAviewNextCell (view, skycell->fpa, 1)) != NULL) {
+                psTrace ("pswarp", 4, "Cell %d: %x %x\n", view->cell, cell->file_exists, cell->process);
+                if (!cell->process || !cell->file_exists) { continue; }
+                if (!pmFPAfileRead (skycell, view, config)) {
+                    psError(psErrorCodeLast(), false, "failed READ at CELL %s", skycell->name);
+                    psFree(view);
+                    return false;
+                }
+                if (!pmFPAfileClose(skycell, view)) {
+                    psError(psErrorCodeLast(), false, "failed CLOSE at CELL %s", skycell->name);
+                    psFree (view);
+                    return false;
+                }
 
-        pmCell *target = pmFPAviewThisCell(view, output->fpa); ///< Target cell
-        pmReadout *readout = pmReadoutAlloc(target); ///< Target readout
-        readout->image = psImageAlloc(numCols, numRows, PS_TYPE_F32);
-        psImageInit(readout->image, NAN);
-        psFree(readout);                // Drop reference
-
-        bool status = false;
-        psMetadataItemSupplement(&status, target->concepts, source->concepts, "CELL.XBIN");
-        psMetadataItemSupplement(&status, target->concepts, source->concepts, "CELL.YBIN");
-        psMetadataItemSupplement(&status, target->concepts, source->concepts, "CELL.XSIZE");
-        psMetadataItemSupplement(&status, target->concepts, source->concepts, "CELL.YSIZE");
-        psMetadataItemSupplement(&status, target->concepts, source->concepts, "CELL.XPARITY");
-        psMetadataItemSupplement(&status, target->concepts, source->concepts, "CELL.YPARITY");
-        psMetadataItemSupplement(&status, target->concepts, source->concepts, "CELL.X0");
-        psMetadataItemSupplement(&status, target->concepts, source->concepts, "CELL.Y0");
-
-        
+                // we've got the output astrom header
+                pmHDU *hdu = pmHDUFromCell(cell); ///< HDU for source
+                if (!hdu || !hdu->header) {
+                    psError(PM_ERR_PROG, false, "Unable to find header for sky cell.");
+                    psFree(view);
+                    return false;
+                }
+                int numCols = psMetadataLookupS32(NULL, hdu->header, "NAXIS1"); ///< Number of columns
+                int numRows = psMetadataLookupS32(NULL, hdu->header, "NAXIS2"); ///< Number of rows
+                if ((numCols == 0) || (numRows == 0)) {
+                    psError(PSWARP_ERR_DATA, false, "skycell has invalid dimensions %d x %d", numCols, numRows);
+                    psFree(view);
+                    return false;
+                }
+                
+                // XXX generate the output pixels
+                pmCell *target = pmFPAviewThisCell(view, output->fpa); ///< Target cell
+                pmReadout *readout = pmReadoutAlloc(target); ///< Target readout
+                readout->image = psImageAlloc(numCols, numRows, PS_TYPE_F32);
+                psImageInit(readout->image, NAN);
+                psFree(readout);                // Drop reference
+                
+                bool status = false;
+                psMetadataItemSupplement(&status, target->concepts, cell->concepts, "CELL.XBIN");
+                psMetadataItemSupplement(&status, target->concepts, cell->concepts, "CELL.YBIN");
+                psMetadataItemSupplement(&status, target->concepts, cell->concepts, "CELL.XSIZE");
+                psMetadataItemSupplement(&status, target->concepts, cell->concepts, "CELL.YSIZE");
+                psMetadataItemSupplement(&status, target->concepts, cell->concepts, "CELL.XPARITY");
+                psMetadataItemSupplement(&status, target->concepts, cell->concepts, "CELL.YPARITY");
+                psMetadataItemSupplement(&status, target->concepts, cell->concepts, "CELL.X0");
+                psMetadataItemSupplement(&status, target->concepts, cell->concepts, "CELL.Y0");
+            }
+            if (!pmFPAfileClose(skycell, view)) {
+                psError(psErrorCodeLast(), false, "failed CLOSE at CHIP %s", skycell->name);
+                psFree (view);
+                return false;
+            }
+        }
+        if (!pmFPAfileClose(skycell, view)) {
+            psError(psErrorCodeLast(), false, "failed CLOSE at FPA %s", skycell->name);
+            psFree (view);
+            return false;
+        }
+        psFree(view);
     }
 
-    // XXX this is not a sufficient test
+    // the section below converts the header astrometry info (in skycell) to the pmFPA
+    // astrometry structures, saving them on the output fpa structure.
+
+    // XXX this block should be modified to loop over all chips
+
+    pmFPAview *view = pmFPAviewAlloc (0);
     view->chip = 0;
     view->cell = 0;
@@ -134 +167 @@
 }
 
-bool pswarpDefineBackground (pmConfig *config) {
-
-    // load the PSWARP recipe
-    psMetadata *recipe = psMetadataLookupPtr (NULL, config->recipes, PSWARP_RECIPE);
-    if (!recipe) {
-        psError(PSWARP_ERR_CONFIG, true, "Can't find PSWARP recipe!\n");
-        return false;
-    }
-
-    // select the input data sources
-    pmFPAfile *skycell = psMetadataLookupPtr (NULL, config->files, "PSWARP.SKYCELL");
-    if (!skycell) {
-        psError(PSWARP_ERR_CONFIG, false, "Can't find skycell data!\n");
-        return false;
-    }
-    pmFPAfile *output = psMetadataLookupPtr(NULL, config->files, "PSWARP.OUTPUT.BKGMODEL");
-    if (!output) {
-        psError(PSWARP_ERR_CONFIG, false, "Can't find output data!\n");
-        return false;
-    }
-    pmFPAfile *input = psMetadataLookupPtr(NULL, config->files, "PSWARP.INPUT");
-    if (!input) {
-        psError(PSWARP_ERR_CONFIG, false, "Can't find input data!\n");
-        return false;
-    }
-    
-    // open the full skycell file; no need to defer different depths. only load the header data
-    pmFPAview *view = pmFPAviewAlloc (0);
-    pmFPAfileOpen (skycell, view, config);
-
-    // Read header and create target
-    {
-        if (!pmFPAReadHeaderSet(skycell->fpa, skycell->fits, config)) {
-            psError(psErrorCodeLast(), false, "Unable to read headers for skycell.");
-            psFree(view);
-            return false;
-        }
-        view->chip = 0;
-        view->cell = 0;
-        view->readout = 0;
-        pmCell *source = pmFPAfileThisCell(config->files, view, "PSWARP.SKYCELL"); ///< Source cell
-        pmHDU *hdu = pmHDUFromCell(source); ///< HDU for source
-        if (!hdu || !hdu->header) {
-            psError(PM_ERR_PROG, false, "Unable to find header for sky cell.");
-            psFree(view);
-            return false;
-        }
-        int numCols = psMetadataLookupS32(NULL, hdu->header, "NAXIS1"); ///< Number of columns
-        int numRows = psMetadataLookupS32(NULL, hdu->header, "NAXIS2"); ///< Number of rows
-        if ((numCols == 0) || (numRows == 0)) {
-            psError(PSWARP_ERR_DATA, false, "skycell has invalid dimensions %d x %d", numCols, numRows);
-            psFree(view);
-            return false;
-        }
-
-        pmCell *target = pmFPAviewThisCell(view, output->fpa); ///< Target cell
-        pmReadout *readout = pmReadoutAlloc(target); ///< Target readout
-        readout->image = psImageAlloc(numCols / output->xBin, numRows / output->yBin, PS_TYPE_F32);
-        psImageInit(readout->image, NAN);
-        psFree(readout);                // Drop reference
-
-        bool status = false;
-        psMetadataItemSupplement(&status, target->concepts, source->concepts, "CELL.XBIN");
-        psMetadataItemSupplement(&status, target->concepts, source->concepts, "CELL.YBIN");
-/*      psMetadataAddS32(target->concepts,PS_LIST_TAIL,"CELL.XBIN",PS_META_REPLACE,"",output->xBin); */
-/*      psMetadataAddS32(target->concepts,PS_LIST_TAIL,"CELL.YBIN",PS_META_REPLACE,"",output->yBin); */
-        psMetadataItemSupplement(&status, target->concepts, source->concepts, "CELL.XSIZE");
-        psMetadataItemSupplement(&status, target->concepts, source->concepts, "CELL.YSIZE");
-/*      psMetadataAddS32(target->concepts,PS_LIST_TAIL,"CELL.XSIZE",PS_META_REPLACE,"",numCols / output->xBin); */
-/*      psMetadataAddS32(target->concepts,PS_LIST_TAIL,"CELL.YSIZE",PS_META_REPLACE,"",numRows / output->yBin); */
-        psMetadataItemSupplement(&status, target->concepts, source->concepts, "CELL.XPARITY");
-        psMetadataItemSupplement(&status, target->concepts, source->concepts, "CELL.YPARITY");
-        psMetadataItemSupplement(&status, target->concepts, source->concepts, "CELL.X0");
-        psMetadataItemSupplement(&status, target->concepts, source->concepts, "CELL.Y0");
-
-        
-    }
-
-    // XXX this is not a sufficient test
-    view->chip = 0;
-    view->cell = 0;
-    view->readout = -1;
-    pmHDU *phu = pmFPAviewThisPHU(view, skycell->fpa); ///< Skycell PHU
-    pmHDU *hdu = pmFPAviewThisHDU(view, skycell->fpa); ///< Skycell header
-
-    pmAstromWCS *WCS = pmAstromWCSfromHeader(hdu->header);
-
-    double cd1f = 1.0 * output->xBin;
-    double cd2f = 1.0 * output->yBin;
-    
-    WCS->crpix1 = WCS->crpix1 / cd1f;
-    WCS->crpix2 = WCS->crpix2 / cd2f;
-    
-    WCS->cdelt1 *= cd1f;
-    WCS->cdelt2 *= cd2f;
-
-    WCS->trans->x->coeff[1][0] *= cd1f;
-    WCS->trans->x->coeff[0][1] *= cd2f;
-    WCS->trans->y->coeff[1][0] *= cd1f;
-    WCS->trans->y->coeff[0][1] *= cd2f;
-    
-    
-    pmAstromWCStoHeader (hdu->header,WCS);
-
-    bool bilevelAstrometry = false;
-    if (phu) {
-        char *ctype = psMetadataLookupStr(NULL, phu->header, "CTYPE1");
-        if (ctype) {
-            bilevelAstrometry = !strcmp(&ctype[4], "-DIS");
-        }
-    }
-
-    // We read from the skycell into the output.  i.e., the output receives the desired astrometry.
-    pmChip *outputChip = pmFPAviewThisChip(view, output->fpa); ///< Chip in the output
-    if (bilevelAstrometry) {
-        if (!pmAstromReadBilevelMosaic(output->fpa, phu->header)) {
-            psError(psErrorCodeLast(), false, "Unable to read bilevel mosaic astrometry for skycell.");
-            psFree(view);
-            return false;
-        }
-        if (!pmAstromReadBilevelChip(outputChip, hdu->header)) {
-            psError(psErrorCodeLast(), false, "Unable to read bilevel chip astrometry for skycell.");
-            psFree(view);
-            return false;
-        }
-    } else {
-        // we use a default FPA pixel scale of 1.0
-      if (!pmAstromReadWCS(output->fpa, outputChip, hdu->header, 1.0)) {
-            psError(psErrorCodeLast(), false, "Unable to read WCS astrometry for skycell.");
-            psFree(view);
-            return false;
-        }
-    }
-
-    
-    view->chip = view->cell = view->readout = -1;
-    pmFPAAddSourceFromView(output->fpa, view, output->format);
-
-
-    psFree (view);
-    return true;
-}

trunk/pswarp/src/pswarpDefineSkycell.c

@@ -73 +73 @@
     psFitsClose(fits);
 
-    // We need to force the format for the skycell to be equivalent to SIMPLE.  Determine
-    // the current format from the header; Determine camera if not specified already
-    // XXX EAM : this operation should be defined as a pmConfig function (pmConfigCopy?)
-    skyConfig = pmConfigAlloc();
-    skyConfig->user = psMemIncrRefCounter(config->user);
-    skyConfig->system = psMemIncrRefCounter(config->system);
+    { // this bit is unique to pswarpDefineSkycell:
 
-    psFree (skyConfig->files);
-    skyConfig->files = psMemIncrRefCounter(config->files);
-    psFree (skyConfig->arguments);
-    skyConfig->arguments = psMemIncrRefCounter(config->arguments);
+      // We need to force the format for the skycell to be equivalent to SIMPLE.  Determine
+      // the current format from the header; Determine camera if not specified already
+      // XXX EAM : this operation should be defined as a pmConfig function (pmConfigCopy?)
+      skyConfig = pmConfigAlloc();
+      skyConfig->user = psMemIncrRefCounter(config->user);
+      skyConfig->system = psMemIncrRefCounter(config->system);
+
+      psFree (skyConfig->files);
+      skyConfig->files = psMemIncrRefCounter(config->files);
+      psFree (skyConfig->arguments);
+      skyConfig->arguments = psMemIncrRefCounter(config->arguments);
+    }
 
     format = pmConfigCameraFormatFromHeader (NULL, NULL, NULL, skyConfig, phu, false);
@@ -94 +97 @@
 
     // Record the camera name of the skycell, so we can save its configuration
-    psMetadataAddStr(config->arguments, PS_LIST_TAIL, "SKYCELL.CAMERA", 0, "Name of camera for skycell",
-                     skyConfig->cameraName);
+    // XXX this can be put outside of the f
+    psMetadataAddStr(config->arguments, PS_LIST_TAIL, "SKYCELL.CAMERA", 0, "Name of camera for skycell", skyConfig->cameraName);
 
     // build the template fpa, set up the basic view

trunk/pswarp/src/pswarpLoop.c

@@ -1 +1 @@
 /** @file pswarpLoop.c
  *
- *  @brief
- *
+ *  @brief main processing loop for pswarp
  *  @ingroup pswarp
  *
@@ -12 +11 @@
 
 #include "pswarp.h"
-#include <ppStats.h>
-#include "pswarpFileNames.h"            // Lists of file rules used at different stages
 
-#define WCS_NONLIN_TOL 0.001            // Non-linear tolerance for header WCS
-#define TESTING 0                       // Testing output?
-
-
-
-// Loop over the inputs, warp them to the output skycell and then write out the output.
+// Loop over the inputs, warp them to the output skycell and then update metadata
 bool pswarpLoop(pmConfig *config, psMetadata *stats)
 {
-    bool status;
-    bool mdok;                          // Status of MD lookup
-
-    const char *skyCamera = psMetadataLookupStr(NULL, config->arguments,
-                                                "SKYCELL.CAMERA");  // Name of camera for skycell
-    pmConfigCamerasCull(config, skyCamera);
-    pmConfigRecipesCull(config, "PSWARP,PPSTATS,PSPHOT,PSASTRO,MASKS,JPEG");
 
     // load the recipe
+    bool status = false;
     psMetadata *recipe = psMetadataLookupPtr (&status, config->recipes, PSWARP_RECIPE);
     if (!recipe) {
@@ -38 +24 @@
     }
 
-    if (!pswarpSetMaskBits(config)) {
-        psError(psErrorCodeLast(), false, "failed to set mask bits");
-        return NULL;
-    }
-
-    // output mask bits
-    psImageMaskType maskValue = psMetadataLookupImageMask(&status, recipe, "MASK.OUTPUT");
-    psAssert (status, "MASK.OUTPUT was not defined");
-
     // select the input data sources
-    pmFPAfile *input = psMetadataLookupPtr(NULL, config->files, "PSWARP.INPUT");
-    if (!input) {
-        psError(PSWARP_ERR_CONFIG, true, "Can't find input data!\n");
+    pmFPAfile *output = psMetadataLookupPtr(&status, config->files, "PSWARP.OUTPUT");
+    if (!output) {
+        psError(PSWARP_ERR_CONFIG, true, "Can't find output data!\n");
         return false;
     }
 
     // use the external astrometry source if supplied
-    pmFPAfile *astrom = psMetadataLookupPtr(NULL, config->files, "PSWARP.ASTROM");
-    if (!astrom) {
-        astrom = input;
-    }
-
-    if (astrom->camera != input->camera) {
-        psError(PSWARP_ERR_DATA, true, "Input camera and astrometry camera do not match.");
+    pmFPAfile *skycell = psMetadataLookupPtr(&status, config->files, "PSWARP.SKYCELL");
+    if (!skycell) {
+        psError(PSWARP_ERR_DATA, true, "Cannot find output astrometry.");
         return false;
     }
 
-    // select the output readout
     pmFPAview *view = pmFPAviewAlloc(0);
-    view->chip = 0;
-    view->cell = 0;
-    view->readout = 0;
-    pmReadout *output = pmFPAfileThisReadout(config->files, view, "PSWARP.OUTPUT");
-    if (!output) {
-        psError(PSWARP_ERR_CONFIG, true, "Can't find output data!\n");
+
+    int nInputs = psMetadataLookupS32(&status, config->arguments, "NUM_INPUTS");
+    if (!status) {
+        psError(PSWARP_ERR_DATA, true, "number of inputs is not defined (programming error)");
         return false;
     }
-    psFree (view);
 
+    // load in the input pixel data (ex. background model)
+    pmFPAfileActivate(config->files, false, NULL);
+    pmFPAfileActivate(config->files, true, "PSWARP.INPUT");
+    pmFPAfileActivate(config->files, true, "PSWARP.MASK");
+    pmFPAfileActivate(config->files, true, "PSWARP.VARIANCE");
 
-    // Turn all skycell files on to generate them, and then turn them off for the loop over the input images
-    // the input, which is in a different format.
-    {
-        pswarpFileActivation(config, detectorFiles, false);
-        pswarpFileActivation(config, photFiles, false);
-        pswarpFileActivation(config, independentFiles, false);
-        pswarpFileActivation(config, skycellFiles, true);
-        if (!pswarpIOChecksBefore(config)) {
-            psError(psErrorCodeLast(), false, "Unable to read files.");
-            goto DONE;
-        }
-        pswarpFileActivation(config, skycellFiles, false);
+    // We re-activate the CMF load so we can transform the sources as well as the pixels.
+    // We only need to read in these if the astrometry source is CMF.
+    pmFPAfileActivate(config->files, true, "PSWARP.ASTROM");
+
+    // loop over this section once per input group
+    for (int i = 0; i < nInputs; i++) {
+        // select the input data sources
+        pmFPAfile *input = pmFPAfileSelectSingle(config->files, "PSWARP.INPUT", i);
+        if (!input) {
+            psError(PSWARP_ERR_CONFIG, true, "Can't find input data!\n");
+            return false;
+        }
+
+        // select the input data sources
+        pmFPAfile *astrom = pmFPAfileSelectSingle(config->files, "PSWARP.ASTROM", i);
+        if (!astrom) {
+            astrom = input;
+        }
+
+        pmFPAviewReset (view);
+
+        // files associated with the science image
+        if (!pmFPAfileIOChecks (config, view, PM_FPA_BEFORE)) {
+            psError(psErrorCodeLast(), false, "Unable to read files.");
+            goto FAIL;
+        }
+
+        // *** main transformation block 
+        // *** this section loops over the input chips/cells and reads them one at a time
+        // *** the output chips/cells are filled where appropriate, but not yet written to disk 
+        pmChip *chip;
+        while ((chip = pmFPAviewNextChip (view, input->fpa, 1)) != NULL) {
+            psTrace ("pswarp", 4, "Chip %d: %x %x\n", view->chip, chip->file_exists, chip->process);
+            if (!chip->process || !chip->file_exists) { continue; }
+            if (!pmFPAfileIOChecks (config, view, PM_FPA_BEFORE)) {
+                psError(psErrorCodeLast(), false, "Unable to read files.");
+                goto FAIL;
+            }
+
+            pmCell *cell;
+            while ((cell = pmFPAviewNextCell (view, input->fpa, 1)) != NULL) {
+                psTrace ("pswarp", 4, "Cell %d: %x %x\n", view->cell, cell->file_exists, cell->process);
+                if (!cell->process || !cell->file_exists) { continue; }
+                if (!pmFPAfileIOChecks (config, view, PM_FPA_BEFORE)) {
+                    psError(psErrorCodeLast(), false, "Unable to read files.");
+                    goto FAIL;
+                }
+
+                // process each of the readouts
+                pmReadout *readout;
+                while ((readout = pmFPAviewNextReadout(view, input->fpa, 1)) != NULL) {
+                    if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) {
+                        psError(psErrorCodeLast(), false, "Unable to read files.");
+                        goto FAIL;
+                    }
+                    if (!readout->data_exists) {
+                        continue;
+                    }
+
+                    // Copy the detections from the astrometry carrier to the input, so they can be accessed by
+                    // pswarpTransformReadout
+                    if (astrom != input) {
+                        pmReadout *astromRO = pmFPAviewThisReadout(view, astrom->fpa); // Readout for astrometry
+                        pmDetections *detections = psMetadataLookupPtr(&status, astromRO->analysis, "PSPHOT.DETECTIONS"); // Sources from astrometry
+                        if (detections) {
+                            psMetadataAddPtr(readout->analysis, PS_LIST_TAIL, "PSPHOT.DETECTIONS", PS_DATA_ARRAY, "Sources from input astrometry", detections);
+                        }
+                    }
+
+                    pswarpTransformToTarget(output->fpa, readout, config, false);
+                
+                    if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) {
+                        psError(psErrorCodeLast(), false, "Unable to write files.");
+                        goto FAIL;
+                    }
+                }
+                if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) {
+                    psError(psErrorCodeLast(), false, "Unable to write files.");
+                    goto FAIL;
+                }
+            }
+            if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) {
+                psError(psErrorCodeLast(), false, "Unable to write files.");
+                goto FAIL;
+            }
+        }
+        if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) {
+            psError(psErrorCodeLast(), false, "Unable to write files.");
+            goto FAIL;
+        }
+
+        if (!pswarpUpdateStatistics (output->fpa, stats, input->fpa, astrom->fpa, config)) {
+            psError(psErrorCodeLast(), false, "problem generating statistics.");
+            goto FAIL;
+        }
+        if (!pswarpUpdateMetadata (output->fpa, skycell->fpa, input->fpa, astrom->fpa, config, true)) {
+            psError(psErrorCodeLast(), false, "problem generating statistics.");
+            goto FAIL;
+        }
     }
 
-    // Read the input astrometry
-    // XXX rather than use the activations here, this should just explicitly loop over the desired filerule
-    {
-        pmFPAfileActivate(config->files, true, "PSWARP.ASTROM");
-
-        pmChip *chip;
-        pmFPAview *view = pmFPAviewAlloc(0);
-        if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) {
-            psError(psErrorCodeLast(), false, "Unable to read files.");
-            goto DONE;
-        }
-        while ((chip = pmFPAviewNextChip (view, input->fpa, 1)) != NULL) {
-            psTrace ("pswarp", 4, "Chip %d: %x %x\n", view->chip, chip->file_exists, chip->process);
-            if (!chip->process || !chip->file_exists) { continue; }
-            if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) {
-                psError(psErrorCodeLast(), false, "Unable to read files.");
-                goto DONE;
-            }
-            pmCell *cell;
-            while ((cell = pmFPAviewNextCell (view, input->fpa, 1)) != NULL) {
-                psTrace ("pswarp", 4, "Cell %d: %x %x\n", view->cell, cell->file_exists, cell->process);
-                if (!cell->process || !cell->file_exists) { continue; }
-                if (!pmFPAfileIOChecks (config, view, PM_FPA_BEFORE) ||
-                    !pmFPAfileIOChecks (config, view, PM_FPA_AFTER)) {
-                    psError(psErrorCodeLast(), false, "Unable to read files.");
-                    goto DONE;
-                }
-            }
-            if (!pmFPAfileIOChecks (config, view, PM_FPA_AFTER)) {
-                psError(psErrorCodeLast(), false, "Unable to write files.");
-                goto DONE;
-            }
-        }
-        if (!pmFPAfileIOChecks (config, view, PM_FPA_AFTER)) {
-            psError(psErrorCodeLast(), false, "Unable to write files.");
-            goto DONE;
-        }
-        psFree(view);
-
-        pswarpFileActivation(config, detectorFiles, true);
-        pmFPAfileActivate(config->files, false, "PSWARP.ASTROM");
+    if (!pswarpMakePSF (config, output, stats)) {
+      psError(psErrorCodeLast(), false, "problem generating PSF.");
+      goto FAIL;
     }
 
-    // Turn on the source output --- we need to get rid of these so that we can measure the PSF
-    pmFPAfileActivate(config->files, true, "PSWARP.OUTPUT.SOURCES");
+    psFree(view);
+    return true;
 
-    // Don't care about the skycell anymore --- we've read it, and that's all we need to do.
-    pmFPAfileActivate(config->files, false, "PSWARP.SKYCELL");
-    view = pmFPAviewAlloc(0);
+ FAIL:
+    psFree (view);
+    return false;
+}
 
-    // find the FPA phu
-    bool bilevelAstrometry = false;
-    pmHDU *phu = pmFPAviewThisPHU(view, astrom->fpa);
-    if (phu) {
-        char *ctype = psMetadataLookupStr(NULL, phu->header, "CTYPE1");
-        if (ctype) {
-            bilevelAstrometry = !strcmp (&ctype[4], "-DIS");
-        }
-    }
-    if (bilevelAstrometry) {
-        if (!pmAstromReadBilevelMosaic(input->fpa, phu->header)) {
-            psError(psErrorCodeLast(), false, "Unable to read bilevel mosaic astrometry for input FPA.");
-            psFree(view);
-            psFree(stats);
-            goto DONE;
-        }
-    }
+// once the output fpa elements have been built, loop over the fpa and generate stats
+// for each readout
+bool pswarpTransformToTarget (pmFPA *output, pmReadout *input, pmConfig *config, bool backgroundWarp)  {
 
-    psList *cells = psListAlloc(NULL);  // List of cells, for concepts averaging
-
-    // files associated with the science image
-    if (!pmFPAfileIOChecks (config, view, PM_FPA_BEFORE)) {
-        psError(psErrorCodeLast(), false, "Unable to read files.");
-        goto DONE;
-    }
-
+    pmFPAview *view = pmFPAviewAlloc(0);
+    
     pmChip *chip;
-    while ((chip = pmFPAviewNextChip (view, input->fpa, 1)) != NULL) {
+    while ((chip = pmFPAviewNextChip (view, output, 1)) != NULL) {
         psTrace ("pswarp", 4, "Chip %d: %x %x\n", view->chip, chip->file_exists, chip->process);
         if (!chip->process || !chip->file_exists) { continue; }
-        if (!pmFPAfileIOChecks (config, view, PM_FPA_BEFORE)) {
-            psError(psErrorCodeLast(), false, "Unable to read files.");
-            goto DONE;
-        }
 
-        // read WCS data from the corresponding header
-        pmHDU *hdu = pmFPAviewThisHDU (view, astrom->fpa);
-
-        
-        if (bilevelAstrometry) {
-            if (!pmAstromReadBilevelChip (chip, hdu->header)) {
-                psError(psErrorCodeLast(), false, "Unable to read bilevel chip astrometry for input FPA.");
-                psFree(view);
-                psFree(stats);
-                goto DONE;
-            }
-        } else {
-            // we use a default FPA pixel scale of 1.0
-            if (!pmAstromReadWCS (input->fpa, chip, hdu->header, 1.0)) {
-                psError(psErrorCodeLast(), false, "Unable to read WCS astrometry for input FPA.");
-                psFree(view);
-                psFree(stats);
-                goto DONE;
-            }
-        }
-        
         pmCell *cell;
-        while ((cell = pmFPAviewNextCell (view, input->fpa, 1)) != NULL) {
+        while ((cell = pmFPAviewNextCell (view, output, 1)) != NULL) {
             psTrace ("pswarp", 4, "Cell %d: %x %x\n", view->cell, cell->file_exists, cell->process);
             if (!cell->process || !cell->file_exists) { continue; }
-            if (!pmFPAfileIOChecks (config, view, PM_FPA_BEFORE)) {
-                psError(psErrorCodeLast(), false, "Unable to read files.");
-                goto DONE;
-            }
-
-            psListAdd(cells, PS_LIST_TAIL, cell);
 
             // process each of the readouts
             pmReadout *readout;
-            while ((readout = pmFPAviewNextReadout(view, input->fpa, 1)) != NULL) {
-                if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) {
-                    psError(psErrorCodeLast(), false, "Unable to read files.");
-                    goto DONE;
-                }
-                if (!readout->data_exists) {
-                    continue;
-                }
-
-                // Copy the detections from the astrometry carrier to the input, so they can be accessed by
-                // pswarpTransformReadout
-                pmReadout *astromRO = pmFPAviewThisReadout(view, astrom->fpa); // Readout for astrometry
-                pmDetections *detections = psMetadataLookupPtr(&mdok, astromRO->analysis, "PSPHOT.DETECTIONS"); // Sources from astrometry
-                if (detections) {
-                    psMetadataAddPtr(readout->analysis, PS_LIST_TAIL, "PSPHOT.DETECTIONS", PS_DATA_ARRAY, "Sources from input astrometry", detections);
-                }
-
-                pswarpTransformReadout(output, readout, config);
-                
-                if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) {
-                    psError(psErrorCodeLast(), false, "Unable to write files.");
-                    goto DONE;
-                }
-            }
-            if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) {
-                psError(psErrorCodeLast(), false, "Unable to write files.");
-                goto DONE;
-            }
-        }
-        if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) {
-            psError(psErrorCodeLast(), false, "Unable to write files.");
-            goto DONE;
-        }
+            while ((readout = pmFPAviewNextReadout(view, output, 1)) != NULL) {
+                pswarpTransformReadout (readout, input, config, backgroundWarp);
+            }
+        }
     }
-
-    if (!output->data_exists) {
-        psWarning("No overlap between input and skycell.");
-        if (stats) {
-            psMetadataAddS32(stats, PS_LIST_TAIL, "QUALITY", PS_META_REPLACE,
-                             "No overlap between input and skycell", PSWARP_ERR_NO_OVERLAP);
-        }
-        psphotFilesActivate(config, false);
-        psFree(cells);
-        psFree(view);
-        goto DONE;
-    }
-    
-    pmCell *outCell = output->parent;   ///< Output cell
-    pmChip *outChip = outCell->parent;  ///< Output chip
-    pmFPA *outFPA = outChip->parent;    ///< Output FP
-
-    if (!pswarpPixelsLit(output, stats, config)) {
-        psError(psErrorCodeLast(), false, "Unable to calculate pixel regions.");
-        psFree(cells);
-        psFree(view);
-        goto DONE;
-    }
-    bool doStats = psMetadataLookupBool(&mdok,recipe,"MASK.STATS");
-    if (doStats) {
-      if (!pswarpMaskStats(output, stats, config)) {
-        psError(psErrorCodeLast(), false, "Unable to calculate mask stats.");
-        psFree(cells);
-        psFree(view);
-        goto DONE;
-      }
-    }
-    // Set covariance matrix for output
-    {
-        psList *covariances = psMetadataLookupPtr(&mdok, output->analysis,
-                                                  PSWARP_ANALYSIS_COVARIANCES); // Covariance matrices
-        psAssert(covariances, "Should be there");
-        psArray *covars = psListToArray(covariances); // Array of covariance matrices
-        psKernel *covar = psImageCovarianceAverage(covars);
-        psFree(covars);
-        psMetadataRemoveKey(output->analysis, PSWARP_ANALYSIS_COVARIANCES);
-
-        // Correct covariance matrix scale for the mean (square root of the) Jacobian
-        double jacobian = psMetadataLookupF64(NULL, output->analysis, PSWARP_ANALYSIS_JACOBIAN); // Jacobian
-        int goodPixels = psMetadataLookupS32(NULL, output->analysis, PSWARP_ANALYSIS_GOODPIX);   // Good pixels
-        jacobian /= goodPixels;
-        output->covariance = psImageCovarianceScale(covar, jacobian);
-        psFree(covar);
-
-        if (output->variance) {
-            psImageCovarianceTransfer(output->variance, output->covariance);
-        }
-    }
-
-    if (!pmConceptsAverageCells(outCell, cells, NULL, NULL, false)) {
-        psError(psErrorCodeLast(), false, "Unable to average cell concepts.");
-        psFree(stats);
-        psFree(cells);
-        psFree(view);
-        goto DONE;
-    }
-    psFree(cells);
-
-    psRegion *trimsec = psMetadataLookupPtr(NULL, outCell->concepts, "CELL.TRIMSEC"); ///< Trim section
-    trimsec->x0 = trimsec->x1 = trimsec->y0 = trimsec->y1 = 0; ///< All pixels
-
-    if (!psMetadataCopy(outFPA->concepts, input->fpa->concepts)) {
-        psError(psErrorCodeLast(), false, "Unable to copy FPA concepts from input to output.");
-        psFree(stats);
-        psFree(view);
-        goto DONE;
-    }
-
-    // Update ZP from the astrometry
-    {
-        psMetadataItem *item = psMetadataLookup(outFPA->concepts, "FPA.ZP");
-        item->data.F32 = psMetadataLookupF32(NULL, astrom->fpa->concepts, "FPA.ZP");
-    }
-
-    pmHDU *hdu = outFPA->hdu;           ///< HDU for the output warped image
-
-    // Copy header from target
-    {
-        pmFPAview *skyView = pmFPAviewAlloc(0); ///< View into skycell
-        skyView->chip = skyView->cell = 0;
-        pmCell *cell = pmFPAfileThisCell(config->files, skyView, "PSWARP.SKYCELL"); // Skycell cell
-        psFree(skyView);
-        pmHDU *skyHDU = pmHDUFromCell(cell); ///< HDU
-        if (!skyHDU) {
-            psError(PSWARP_ERR_DATA, false, "Unable to find skycell HDU.");
-            psFree(view);
-            goto DONE;
-        }
-        hdu->header = psMetadataCopy(hdu->header, skyHDU->header);
-    }
-
-    pswarpVersionHeader(hdu->header);
-    
-    if (!pmAstromWriteWCS(hdu->header, outFPA, outChip, WCS_NONLIN_TOL)) {
-        psError(psErrorCodeLast(), false, "Unable to generate WCS header.");
-        psFree(stats);
-        goto DONE;
-    }
-
-    if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) {
-        psError(psErrorCodeLast(), false, "Unable to write files.");
-        goto DONE;
-    }
-
-    // Done with the detector side of things
-    pswarpFileActivation(config, detectorFiles, false);
-    pswarpFileActivation(config, independentFiles, false);
-
-
-    // We need a new PSF model for the warped frame.  It would be good to generate this analytically, but
-    // that's going to be tricky.  We have a list of sources, so we use those to redetermine the PSF model.
-
-    if (psMetadataLookupBool(&mdok, recipe, "PSF")) {
-        pswarpFileActivation(config, photFiles, true);
-        if (!pswarpIOChecksBefore(config)) {
-            psError(psErrorCodeLast(), false, "Unable to read files.");
-            goto DONE;
-        }
-
-        // supply the readout and fpa of interest to psphot
-        pmFPAfile *photFile = psMetadataLookupPtr(NULL, config->files, "PSPHOT.INPUT");
-        pmFPACopy(photFile->fpa, outFPA);
-
-        pmFPAview *view = pmFPAviewAlloc(0); ///< View into skycell
-        view->chip = view->cell = view->readout = 0;
-
-        // grab the sources of interest from the storage location (pmFPAfile PSPHOT.INPUT.CMF)
-        psArray *sources = psphotLoadPSFSources (config, view);
-        if (!sources) {
-            psError(psErrorCodeLast(), false, "No sources supplied to measure PSF");
-            goto DONE;
-        }
-
-        pmModelClassSetLimits(PM_MODEL_LIMITS_STRICT);
-
-        // measure the PSF using these sources
-        if (!psphotReadoutFindPSF(config, view, "PSPHOT.INPUT", sources)) {
-            // 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.");
-            if (stats && psMetadataLookupS32(NULL, stats, "QUALITY") == 0) {
-                psMetadataAddS32(stats, PS_LIST_TAIL, "QUALITY", PS_META_REPLACE,
-                                 "Unable to determine PSF", psErrorCodeLast());
-            }
-            psErrorClear();
-            psphotFilesActivate(config, false);
-        }
-
-        // Ensure seeing is carried over
-        pmChip *photChip = pmFPAviewThisChip(view, photFile->fpa);                 // Chip with seeing
-        psMetadataItem *item = psMetadataLookup(outChip->concepts, "CHIP.SEEING"); // Concept with seeing
-        item->data.F32 = psMetadataLookupF32(NULL, photChip->concepts, "CHIP.SEEING");
-
-// XXX EAM : put this in a visualization function
-#if (TESTING)
-        {
-            #define PSF_SIZE 20         ///< Half-size of PSF
-            #define PSF_FLUX 10000      ///< Central flux for PSF
-            pmChip *photChip = pmFPAviewThisChip(view, photFile->fpa);
-            pmPSF *psf = psMetadataLookupPtr(NULL, photChip->analysis, "PSPHOT.PSF");
-            psImage *image = psImageAlloc(2 * PSF_SIZE + 1, 2 * PSF_SIZE + 1, PS_TYPE_F32);
-            psImageInit(image, 0);
-            pmModel *model = pmModelFromPSFforXY(psf, PSF_SIZE, PSF_SIZE, PSF_FLUX);
-            pmModelAdd(image, NULL, model, PM_MODEL_OP_FULL, 0);
-            psFree(model);
-            psFits *fits = psFitsOpen("psf.fits", "w");
-            psFitsWriteImage(fits, NULL, image, 0, NULL);
-            psFitsClose(fits);
-            psFree(image);
-        }
-#endif
-
-        psFree(view);
-    }
-
-    // Perform statistics on the output image
-    if (stats) {
-        if (!ppStatsFPA(stats, output->parent->parent->parent, view, maskValue, config)) {
-            psWarning("Unable to perform statistics on warped image.");
-        }
-    }
-    
-
-    // Add MD5 information for readout
-    const char *chipName = psMetadataLookupStr(NULL, output->parent->parent->concepts, "CHIP.NAME");
-    const char *cellName = psMetadataLookupStr(NULL, output->parent->concepts, "CELL.NAME");
-    psString headerName = NULL; ///< Header name for MD5
-    psStringAppend(&headerName, "MD5_%s_%s_%d", chipName, cellName, view->readout);
-    psVector *md5 = psImageMD5(output->image); ///< md5 hash
-    psString md5string = psMD5toString(md5); ///< String
-    psFree(md5);
-    psMetadataAddStr(hdu->header, PS_LIST_TAIL, headerName, PS_META_REPLACE,
-                     "Image MD5", md5string);
-    psFree(md5string);
-    psFree(headerName);
-    psFree(view);
-
- DONE:
-
     return true;
 }
 
-// Loop over the inputs, warp them to the output skycell and then write out the output.
-bool pswarpLoopBackground(pmConfig *config, psMetadata *stats)
-{
-    bool status;
-    bool mdok;                          // Status of MD lookup
-    const char *skyCamera = psMetadataLookupStr(NULL, config->arguments,
-                                                "SKYCELL.CAMERA");  // Name of camera for skycell
-    pmConfigCamerasCull(config, skyCamera);
-    pmConfigRecipesCull(config, "PSWARP,PPSTATS,PSPHOT,PSASTRO,MASKS,JPEG");
-
-    // load the recipe
-    psMetadata *recipe = psMetadataLookupPtr (&status, config->recipes, PSWARP_RECIPE);
-    if (!recipe) {
-        psError(PSWARP_ERR_CONFIG, false, "missing recipe %s", PSWARP_RECIPE);
-        return false;
-    }
-
-    if (!pswarpSetMaskBits(config)) {
-        psError(psErrorCodeLast(), false, "failed to set mask bits");
-        return NULL;
-    }
-
-    // select the input data sources
-    pmFPAfile *input = psMetadataLookupPtr(NULL, config->files, "PSWARP.BKGMODEL");
-    if (!input) {
-        psError(PSWARP_ERR_CONFIG, true, "Can't find input data!\n");
-        return false;
-    }
-
-    // use the external astrometry source if supplied
-    pmFPAfile *astrom = psMetadataLookupPtr(NULL, config->files, "PSWARP.ASTROM");
-    if (!astrom) {
-        astrom = input;
-    }
-
-    if (astrom->camera != input->camera) {
-        psError(PSWARP_ERR_DATA, true, "Input camera and astrometry camera do not match.");
-        return false;
-    }
-
-    // select the output readout
-    pmFPAview *view = pmFPAviewAlloc(0);
-    view->chip = 0;
-    view->cell = 0;
-    view->readout = 0;
-    pmReadout *output = pmFPAfileThisReadout(config->files, view, "PSWARP.OUTPUT.BKGMODEL");
-    if (!output) {
-        psError(PSWARP_ERR_CONFIG, true, "Can't find output background data!\n");
-        return false;
-    }
-    psFree (view);
-    // Turn all skycell files on to generate them, and then turn them off for the loop over the input images
-    // the input, which is in a different format.
-    {
-        pswarpFileActivation(config, detectorFiles, false);
-        pswarpFileActivation(config, photFiles, false);
-        pswarpFileActivation(config, independentFiles, false);
-        pswarpFileActivation(config, skycellFiles, true);
-        if (!pswarpIOChecksBefore(config)) {
-            psError(psErrorCodeLast(), false, "Unable to read files.");
-            goto DONE;
-        }
-        pswarpFileActivation(config, skycellFiles, false);
-    }
-    // Read the input astrometry
-    // XXX rather than use the activations here, this should just explicitly loop over the desired filerule
-    {
-
-      pmFPAfileActivate(config->files, true, "PSWARP.ASTROM");
-
-        pmChip *chip;
-        pmFPAview *view = pmFPAviewAlloc(0);
-        if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) {
-            psError(psErrorCodeLast(), false, "Unable to read files.");
-            goto DONE;
-        }
-
-        while ((chip = pmFPAviewNextChip (view, input->fpa, 1)) != NULL) {
-#if 0
-          // This needs to be removed because otherwise it throws an error of duplicate PSPHOT.DETECTIONS.
-            if (!chip->process || !chip->file_exists) { continue; }
-            if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) {
-                psError(psErrorCodeLast(), false, "Unable to read files.");
-                goto DONE;
-            }
-#endif
-            pmCell *cell;
-            while ((cell = pmFPAviewNextCell (view, input->fpa, 1)) != NULL) {
-              psTrace ("pswarp", 4, "ACell %d: %x %x %d\n", view->cell, cell->file_exists, cell->process,psErrorCodeLast());
-                if (!cell->process || !cell->file_exists) { continue; }
-                if (!pmFPAfileIOChecks (config, view, PM_FPA_BEFORE) ||
-                    !pmFPAfileIOChecks (config, view, PM_FPA_AFTER)) {
-                    psError(psErrorCodeLast(), false, "Unable to read files.");
-                    goto DONE;
-                }
-            }
-            if (!pmFPAfileIOChecks (config, view, PM_FPA_AFTER)) {
-                psError(psErrorCodeLast(), false, "Unable to write files.");
-                goto DONE;
-            }
-        }
-        if (!pmFPAfileIOChecks (config, view, PM_FPA_AFTER)) {
-            psError(psErrorCodeLast(), false, "Unable to write files.");
-            goto DONE;
-        }
-        psFree(view);
-        pswarpFileActivation(config, detectorFiles, true);
-        pmFPAfileActivate(config->files, false, "PSWARP.ASTROM");
-    }
-
-    // Don't care about the skycell anymore --- we've read it, and that's all we need to do.
-    pmFPAfileActivate(config->files, false, "PSWARP.SKYCELL");
-    view = pmFPAviewAlloc(0);
-
-    // find the FPA phu
-    bool bilevelAstrometry = false;
-    pmHDU *phu = pmFPAviewThisPHU(view, astrom->fpa);
-
-    //    pmAstromWCS *WCSF = pmAstromWCSfromHeader(phu->header);
-    
-    if (phu) {
-        char *ctype = psMetadataLookupStr(NULL, phu->header, "CTYPE1");
-        if (ctype) {
-            bilevelAstrometry = !strcmp (&ctype[4], "-DIS");
-        }
-    }
-    if (bilevelAstrometry) {
-        if (!pmAstromReadBilevelMosaic(input->fpa, phu->header)) {
-            psError(psErrorCodeLast(), false, "Unable to read bilevel mosaic astrometry for input FPA.");
-            psFree(view);
-            psFree(stats);
-            goto DONE;
-        }
-    }
-
-    psList *cells = psListAlloc(NULL);  // List of cells, for concepts averaging
-
-    // files associated with the science image
-    if (!pmFPAfileIOChecks (config, view, PM_FPA_BEFORE)) {
-        psError(psErrorCodeLast(), false, "Unable to read files.");
-        goto DONE;
-    }
-    pmChip *chip;
-    while ((chip = pmFPAviewNextChip (view, input->fpa, 1)) != NULL) {
-        psTrace ("pswarp", 4, "DChip %d: %x %x\n", view->chip, chip->file_exists, chip->process);
-        if (!chip->process || !chip->file_exists) { continue; }
-        if (!pmFPAfileIOChecks (config, view, PM_FPA_BEFORE)) {
-            psError(psErrorCodeLast(), false, "Unable to read files.");
-            goto DONE;
-        }
-
-        // Pull information from the header of the background files so we can use it to set values.
-        pmHDU *hdu = pmFPAviewThisHDU(view,input->fpa);
-        psMetadata *header = hdu->header;
-        
-        int IMAXIS1 = psMetadataLookupS32(NULL,header,"IMNAXIS1");
-        int IMAXIS2 = psMetadataLookupS32(NULL,header,"IMNAXIS2");
-        int NAXIS1 = psMetadataLookupS32(NULL,header,"NAXIS1");
-        int NAXIS2 = psMetadataLookupS32(NULL,header,"NAXIS2");
-        char *CCDSUM = psMetadataLookupStr(NULL,header,"CCDSUM");
-        int CCDSUM1 = atoi(strtok(CCDSUM," "));
-        int CCDSUM2 = atoi(strtok(NULL," "));
-        
-        psMetadataAddF32(config->arguments,PS_LIST_TAIL,"BKG_WARP_XOFFSET", PS_META_REPLACE,
-                         "xoffset for background model data", (NAXIS1 * CCDSUM1 - IMAXIS1) / (2.0 * CCDSUM1));
-        psMetadataAddF32(config->arguments,PS_LIST_TAIL,"BKG_WARP_YOFFSET", PS_META_REPLACE,
-                         "yoffset for background model data", (NAXIS2 * CCDSUM2 - IMAXIS2) / (2.0 * CCDSUM2));
-        psTrace("pswarp",5,"%d %d %d %d %d %d %g %g %d %d",
-                psMetadataLookupS32(NULL,header,"IMNAXIS1"),
-                psMetadataLookupS32(NULL,header,"IMNAXIS2"),
-                psMetadataLookupS32(NULL,header,"NAXIS1"),
-                psMetadataLookupS32(NULL,header,"NAXIS2"),
-                CCDSUM1,CCDSUM2,
-                psMetadataLookupF32(NULL,config->arguments,"BKG_WARP_XOFFSET"),
-                psMetadataLookupF32(NULL,config->arguments,"BKG_WARP_YOFFSET"),
-                psMetadataLookupS32(NULL,config->arguments,"BKG.XGRID"),
-                psMetadataLookupS32(NULL,config->arguments,"BKG.YGRID"));
-        
-        
-        // read WCS data from the corresponding header
-        hdu = pmFPAviewThisHDU (view, astrom->fpa);
-
-        pmAstromWCS *WCS = pmAstromWCSfromHeader(hdu->header);
-
-        double cd1f = (1.0 * CCDSUM1);// * (1.0 * psMetadataLookupS32(NULL,config->arguments,"BKG.XGRID"));
-        double cd2f = (1.0 * CCDSUM2);// * (1.0 * psMetadataLookupS32(NULL,config->arguments,"BKG.YGRID"));
-
-        WCS->cdelt1 *= cd1f;
-        WCS->cdelt2 *= cd2f;
-        WCS->crpix1 = WCS->crpix1 / cd1f;
-        WCS->crpix2 = WCS->crpix2 / cd2f;
-
-        // WCS->trans->x->nX/nY
-        for (int q = 0; q <= WCS->trans->x->nX; q++) {
-          for (int r = 0; r <= WCS->trans->x->nY; r++) {
-            WCS->trans->x->coeff[q][r] *= pow(cd1f,q) * pow(cd2f,r);
-          }
-        }
-        for (int q = 0; q <= WCS->trans->y->nX; q++) {
-          for (int r = 0; r <= WCS->trans->y->nY; r++) {
-            WCS->trans->y->coeff[q][r] *= pow(cd1f,q) * pow(cd2f,r);
-          }
-        }
-        pmAstromWCStoHeader (hdu->header,WCS);
-        
-        if (bilevelAstrometry) {
-            if (!pmAstromReadBilevelChip (chip, hdu->header)) {
-                psError(psErrorCodeLast(), false, "Unable to read bilevel chip astrometry for input FPA.");
-                psFree(view);
-                psFree(stats);
-                goto DONE;
-            }
-        } else {
-            // we use a default FPA pixel scale of 1.0
-            if (!pmAstromReadWCS (astrom->fpa, chip, hdu->header, 1.0)) {
-                psError(psErrorCodeLast(), false, "Unable to read WCS astrometry for input FPA.");
-                psFree(view);
-                psFree(stats);
-                goto DONE;
-            }
-        }
-        // Modify structure here.
-
-        
-        pmCell *cell;
-        while ((cell = pmFPAviewNextCell (view, input->fpa, 1)) != NULL) {
-            psTrace ("pswarp", 4, "DCell %d: %x %x\n", view->cell, cell->file_exists, cell->process);
-            if (!cell->process || !cell->file_exists) { continue; }
-            if (!pmFPAfileIOChecks (config, view, PM_FPA_BEFORE)) {
-                psError(psErrorCodeLast(), false, "Unable to read files.");
-                goto DONE;
-            }
-
-            psListAdd(cells, PS_LIST_TAIL, cell);
-
-            // process each of the readouts
-            pmReadout *readout;
-            while ((readout = pmFPAviewNextReadout(view, input->fpa, 1)) != NULL) {
-                if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) {
-                    psError(psErrorCodeLast(), false, "Unable to read files.");
-                    goto DONE;
-                }
-                if (!readout->data_exists) {
-                    continue;
-                }
-
-                for (int x = 0; x < readout->image->numCols; x++) {
-                  for (int y = 0; y < readout->image->numRows; y++) {
-                    readout->image->data.F32[y][x] = readout->image->data.F32[y][x] * (cd1f * cd2f) /
-                      (psMetadataLookupS32(&mdok,config->arguments,"BKG.XGRID") *
-                       psMetadataLookupS32(&mdok,config->arguments,"BKG.YGRID"));
-                  }
-                }
-                
-                psMetadataAddS32(config->arguments,PS_LIST_TAIL, "INTERPOLATION.MODE", PS_META_REPLACE, "", 8);
-
-                psMetadataAddBool(config->arguments,PS_LIST_TAIL, "BACKGROUND_WARPING", PS_META_REPLACE, "", true);
-                pswarpTransformReadout(output, readout, config);
-                psMetadataAddBool(config->arguments,PS_LIST_TAIL, "BACKGROUND_WARPING", PS_META_REPLACE, "", false);
-
-                if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) {
-                    psError(psErrorCodeLast(), false, "Unable to write files.");
-                    goto DONE;
-                }
-            }
-            if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) {
-                psError(psErrorCodeLast(), false, "Unable to write files.");
-                goto DONE;
-            }
-        }
-        if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) {
-            psError(psErrorCodeLast(), false, "Unable to write files.");
-            goto DONE;
-        }
-    }
-    if (!output->data_exists) {
-        psWarning("No overlap between input and skycell.");
-        psphotFilesActivate(config, false);
-        psFree(cells);
-        psFree(view);
-        goto DONE;
-    }
-    pmCell *outCell = output->parent;   ///< Output cell
-    pmChip *outChip = outCell->parent;  ///< Output chip
-    pmFPA *outFPA = outChip->parent;    ///< Output FP
-
-/*     if (!pswarpPixelsLit(output, stats, config)) { */
-/*         psError(psErrorCodeLast(), false, "Unable to calculate pixel regions."); */
-/*         psFree(cells); */
-/*         psFree(view); */
-/*         goto DONE; */
-/*     } */
-    psRegion *trimsec = psMetadataLookupPtr(NULL, outCell->concepts, "CELL.TRIMSEC"); ///< Trim section
-    trimsec->x0 = trimsec->x1 = trimsec->y0 = trimsec->y1 = 0; ///< All pixels
-
-    if (!psMetadataCopy(outFPA->concepts, input->fpa->concepts)) {
-        psError(psErrorCodeLast(), false, "Unable to copy FPA concepts from input to output.");
-        psFree(stats);
-        psFree(view);
-        goto DONE;
-    }
-
-    pmHDU *hdu = outFPA->hdu;           ///< HDU for the output warped image
-
-    // Copy header from target
-    {
-        pmFPAview *skyView = pmFPAviewAlloc(0); ///< View into skycell
-        skyView->chip = skyView->cell = 0;
-        pmCell *cell = pmFPAfileThisCell(config->files, skyView, "PSWARP.SKYCELL"); // Skycell cell
-        psFree(skyView);
-        pmHDU *skyHDU = pmHDUFromCell(cell); ///< HDU
-        if (!skyHDU) {
-            psError(PSWARP_ERR_DATA, false, "Unable to find skycell HDU.");
-            psFree(view);
-            goto DONE;
-        }
-        hdu->header = psMetadataCopy(hdu->header, skyHDU->header);
-    }
-    pswarpVersionHeader(hdu->header);
-    
-    if (!pmAstromWriteWCS(hdu->header, outFPA, outChip, WCS_NONLIN_TOL)) {
-        psError(psErrorCodeLast(), false, "Unable to generate WCS header.");
-        psFree(stats);
-        goto DONE;
-    }
-    if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) {
-        psError(psErrorCodeLast(), false, "Unable to write files.");
-        goto DONE;
-    }
-    // Done with the detector side of things
-    pswarpFileActivation(config, detectorFiles, false);
-    pswarpFileActivation(config, independentFiles, false);
-
-
-    // Add MD5 information for readout
-    const char *chipName = psMetadataLookupStr(NULL, output->parent->parent->concepts, "CHIP.NAME");
-    const char *cellName = psMetadataLookupStr(NULL, output->parent->concepts, "CELL.NAME");
-    psString headerName = NULL; ///< Header name for MD5
-    psStringAppend(&headerName, "MD5_%s_%s_%d", chipName, cellName, view->readout);
-    psVector *md5 = psImageMD5(output->image); ///< md5 hash
-    psString md5string = psMD5toString(md5); ///< String
-    psFree(md5);
-    psMetadataAddStr(hdu->header, PS_LIST_TAIL, headerName, PS_META_REPLACE,
-                     "Image MD5", md5string);
-    psFree(md5string);
-    psFree(headerName);
-    psFree(view);
- DONE:
-
-    return true;
-}

trunk/pswarp/src/pswarpMaskStats.c

@@ -1 +1 @@
 #include "pswarp.h"
-#include <ppStats.h>
 
 bool pswarpMaskStats(const pmReadout *readout, psMetadata *stats, const pmConfig *config)
@@ -31 +30 @@
                               dynamicMaskVal,advisoryMaskVal)) {
     psError(PS_ERR_UNKNOWN, false, "Unable to calculate masks for readout.");
-    return(false);
+    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);
+
+  // XXX with multiple inputs (eg, output stacks -> exposure), these only represent the last input
+  psMetadataAddS32(stats, PS_LIST_TAIL,"MASKFRAC_NPIX",     PS_META_REPLACE, "Number of valid pixels", Npix_valid);
+  psMetadataAddF32(stats,PS_LIST_TAIL, "MASKFRAC_STATIC",   PS_META_REPLACE, "Fraction of pixels statically masked", (float) Npix_static / Npix_valid);
+  psMetadataAddF32(stats,PS_LIST_TAIL, "MASKFRAC_DYNAMIC",  PS_META_REPLACE, "Fraction of pixels dynamically masked", (float) Npix_dynamic / Npix_valid);
+  psMetadataAddF32(stats,PS_LIST_TAIL, "MASKFRAC_MAGIC",    PS_META_REPLACE, "Fraction of pixels magically masked", (float) Npix_magic / Npix_valid);
+  psMetadataAddF32(stats,PS_LIST_TAIL, "MASKFRAC_ADVISORY", PS_META_REPLACE, "Fraction of pixels masked as an advisory", (float) Npix_advisory / Npix_valid);
+  return true;
 }

trunk/pswarp/src/pswarpMatchRange.c

@@ -115 +115 @@
     psPlaneTransformApply (srcFP, fpaSrc->fromTPA, srcTP); 
     psPlaneTransformApply (srcPix, chipSrc->fromFPA, srcFP); 
-    fprintf (stderr, "-> %f,%f ", srcPix->x, srcPix->y);
+    fprintf (stderr, "-> %f,%f\n", srcPix->x, srcPix->y);
 # endif
 

trunk/pswarp/src/pswarpParseCamera.c

@@ -13 +13 @@
 #include "pswarp.h"
 
-// Define an input file
-static pmFPAfile *defineInputFile(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;
-
-    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) {
-        return NULL;
-    }
-
-    if (file->type != fileType) {
-        psError(PSWARP_ERR_CONFIG, true, "%s is not of type %s", filerule, pmFPAfileStringFromType(fileType));
-        return NULL;
-    }
-
-    return file;
-}
-
-
-
+static bool foundAstrom = false;
+static bool foundVariance = false;
+static bool foundBackground = false;
 
 bool pswarpParseCamera(pmConfig *config)
 {
     psAssert(config, "Require configuration");
-    bool mdok;                          // Status of MD lookup
-
-    // The input image(s) is required: it defines the camera
-    pmFPAfile *input = defineInputFile(config, NULL, "PSWARP.INPUT", "INPUT", PM_FPA_FILE_IMAGE);
-    if (!input) {
-        psError(psErrorCodeLast(), false, "Failed to build FPA from PSWARP.INPUT");
-        return false;
-    }
-
-    pmFPAfile *astrom = defineInputFile(config, NULL, "PSWARP.ASTROM", "ASTROM", PM_FPA_FILE_CMF);
-    psLogMsg("pswarp", PS_LOG_INFO, "Astrometry source: %s", astrom ? "supplied" : "header");
-
-    pmFPAfile *inMask = defineInputFile(config, input, "PSWARP.MASK", "MASK", PM_FPA_FILE_MASK);
-    if (!inMask) {
-        psLogMsg("pswarp", PS_LOG_INFO, "No mask supplied");
-    }
-
-    pmFPAfile *inVariance = defineInputFile(config, input, "PSWARP.VARIANCE", "VARIANCE",
-                                            PM_FPA_FILE_VARIANCE);
-    if (!inVariance) {
-        psLogMsg("pswarp", PS_LOG_INFO, "No variance supplied");
-    }
-
-    pmFPAfile *inBackground = defineInputFile(config, NULL, "PSWARP.BKGMODEL", "BACKGROUND",
-                                              PM_FPA_FILE_IMAGE);
-    if (!inBackground) {
-      psLogMsg("pswarp", PS_LOG_INFO, "No background models supplied");
-    }
-    
+    bool status;                          // Status of MD lookup
+
+    // *** parse the input information (either from -file or from -input)
+
+    // if INPUTS exists, we have a metadata file to provide input, variance, mask, etc
+    psMetadata *inputFile = psMetadataLookupPtr(&status, config->arguments, "INPUTS");
+    if (inputFile) {
+        if (!pswarpParseMultiInput (config, inputFile)) {
+            psError(PSWARP_ERR_CONFIG, false, "failed to parse multi-input file");
+            return false;
+        }
+    } else {
+        if (!pswarpParseSingleInput (config)) {
+            psError(PSWARP_ERR_CONFIG, false, "failed to parse multi-input file");
+            return false;
+        }
+    }
+
+    // once we have read the input mask headers (if any), we can use that to set the mask bits
+    if (!pswarpSetMaskBits(config)) {
+        psError(psErrorCodeLast(), false, "failed to set mask bits");
+        return false;
+    }
+
+    // *** parse the output information (output target, output astrometry [skycell])
+
     // The input skycell is a required argument: it defines the output image
-    // XXX we may need a different skycell structure here
+    pmConfig *skyConfig = NULL;
     pmFPAfile *skycell = NULL;
-    pmConfig *skyConfig = NULL;
-    bool status = pswarpDefineSkycell(&skycell, &skyConfig, config, "PSWARP.SKYCELL", "SKYCELL");
-    if (!status) {
-        psError(psErrorCodeLast(), false, "Failed to build FPA from PSWARP.SKYCELL");
-        return false;
-    }
+
+    skycell = pmFPAfileDefineNewConfig(&status, &skyConfig, config, "PSWARP.SKYCELL", "SKYCELL");
+    if (!status || !skycell) {
+      psError(psErrorCodeLast(), false, "Failed to build FPA from PSWARP.SKYCELL");
+      return false;
+      }
+    psMetadataAddStr(config->arguments, PS_LIST_TAIL, "SKYCELL.CAMERA", 0, "Name of camera for skycell", skyConfig->cameraName);
     psFree(skyConfig);
 
@@ -107 +64 @@
     }
 
+    // use the skycell camera or the input camera?
+    bool useInputCamera = !strcmp (skycell->cameraName, "SIMPLE");
+
     // The output skycell
-    pmFPAfile *output = pmFPAfileDefineSkycell(config, NULL, "PSWARP.OUTPUT");
+    pmFPAfile *output = useInputCamera ? 
+        pmFPAfileDefineSkycell(config, NULL, "PSWARP.OUTPUT") :
+        pmFPAfileDefineOutputForFormat(config, NULL, "PSWARP.OUTPUT", skycell->cameraName, skycell->formatName);
     if (!output) {
         psError(psErrorCodeLast(), false, "Failed to build FPA from PSWARP.OUTPUT");
@@ -114 +76 @@
     }
     output->save = true;
-
-    pmFPAfile *outMask = pmFPAfileDefineSkycell(config, output->fpa, "PSWARP.OUTPUT.MASK");
+    pmFPAAddSourceFromFormat(output->fpa, output->format); // ** builds the HDUs, is this OK?
+
+    pmFPAfile *outMask = useInputCamera ?
+        pmFPAfileDefineSkycell(config, output->fpa, "PSWARP.OUTPUT.MASK"):
+        pmFPAfileDefineOutputForFormat(config, output->fpa, "PSWARP.OUTPUT.MASK", skycell->cameraName, skycell->formatName);
     if (!outMask) {
         psError(psErrorCodeLast(), false, "Failed to build FPA from PSWARP.OUTPUT.MASK");
@@ -122 +87 @@
     outMask->save = true;
 
-    if (inVariance) {
-        pmFPAfile *outVariance = pmFPAfileDefineSkycell(config, output->fpa, "PSWARP.OUTPUT.VARIANCE");
+    // only create an output variance in we supply an input variance
+    if (foundVariance) {
+        pmFPAfile *outVariance = useInputCamera ? 
+            pmFPAfileDefineSkycell(config, output->fpa, "PSWARP.OUTPUT.VARIANCE"):
+            pmFPAfileDefineOutputForFormat(config, output->fpa, "PSWARP.OUTPUT.VARIANCE", skycell->cameraName, skycell->formatName);
         if (!outVariance) {
             psError(psErrorCodeLast(), false, "Failed to build FPA from PSWARP.OUTPUT.VARIANCE");
@@ -131 +99 @@
     }
 
-    if (astrom && psMetadataLookupBool(&mdok, recipe, "SOURCES")) {
-        pmFPAfile *outSources = pmFPAfileDefineSkycell(config, output->fpa, "PSWARP.OUTPUT.SOURCES");
+    // XXX we assume input sources come from the input astrom description, but this need not be true (we could define an input sources file)
+    if (foundAstrom && psMetadataLookupBool(&status, recipe, "SOURCES")) {
+        pmFPAfile *outSources = useInputCamera ? 
+            pmFPAfileDefineSkycell(config, output->fpa, "PSWARP.OUTPUT.SOURCES"):
+            pmFPAfileDefineOutputForFormat(config, output->fpa, "PSWARP.OUTPUT.SOURCES", skycell->cameraName, skycell->formatName);
         if (!outSources) {
             psError(psErrorCodeLast(), false, "Failed to build FPA from PSWARP.OUTPUT.SOURCES");
@@ -140 +111 @@
     }
 
-    if (inBackground && psMetadataLookupBool(&mdok, recipe, "BACKGROUND.MODEL")) {
-      pmFPAfile *outBackground = pmFPAfileDefineSkycell(config, NULL, "PSWARP.OUTPUT.BKGMODEL");
-/*       pmFPAfile *outBackground = pmFPAfileDefineFromFPA(config,output->fpa, */
-/*                                                      psMetadataLookupS32(&mdok,recipe,"BKG.XGRID"), */
-/*                                                      psMetadataLookupS32(&mdok,recipe,"BKG.YGRID"), */
-/*                                                      "PSWARP.OUTPUT.BKGMODEL"); */
-      outBackground->xBin = psMetadataLookupS32(&mdok, recipe, "BKG.XGRID");
-      outBackground->yBin = psMetadataLookupS32(&mdok, recipe, "BKG.YGRID");
-      
-      if (!outBackground) {
-        psError(psErrorCodeLast(), false, "Failed to build FPA from PSWARP.OUTPUT.BKGMODEL");
-        return false;
-      }
-      outBackground->save = true;
-    }
-    
-    if (psMetadataLookupBool(&mdok, recipe, "PSF")) {
-        // This file, PSPHOT.INPUT, is just used as a carrier; output files (eg, PSPHOT.RESID) are defined by
-        // psphotDefineFiles
-        pmFPAfile *psphotInput = pmFPAfileDefineSkycell(config, NULL, "PSPHOT.INPUT");
+    // only create an output background if an input background is supplied
+    if (foundBackground && psMetadataLookupBool(&status, recipe, "BACKGROUND.MODEL")) {
+        // pmFPAfile *outBackground = pmFPAfileDefineSkycell(config, NULL, "PSWARP.OUTPUT.BKGMODEL");
+        pmFPAfile *outBackground = useInputCamera ? 
+            pmFPAfileDefineSkycell(config, NULL, "PSWARP.OUTPUT.BKGMODEL"):
+            pmFPAfileDefineOutputForFormat(config, NULL, "PSWARP.OUTPUT.BKGMODEL", skycell->cameraName, skycell->formatName);
+        if (!outBackground) {
+            psError(psErrorCodeLast(), false, "Failed to build FPA from PSWARP.OUTPUT.BKGMODEL");
+            return false;
+        }
+        outBackground->xBin = psMetadataLookupS32(&status, recipe, "BKG.XGRID");
+        outBackground->yBin = psMetadataLookupS32(&status, recipe, "BKG.YGRID");
+        outBackground->save = true;
+        pmFPAAddSourceFromFormat(outBackground->fpa, outBackground->format); // ** builds the HDUs, is this OK?
+    }
+
+    if (psMetadataLookupBool(&status, recipe, "PSF")) {
+        // This file, PSPHOT.INPUT, is just used as a carrier; output files (eg, PSPHOT.RESID) are defined by psphotDefineFiles
+        pmFPAfile *psphotInput = useInputCamera ? 
+            pmFPAfileDefineSkycell(config, NULL, "PSPHOT.INPUT"):
+            pmFPAfileDefineOutputForFormat(config, NULL, "PSPHOT.INPUT", skycell->cameraName, skycell->formatName);
         if (!psphotInput) {
             psError(psErrorCodeLast(), false, _("Unable to generate output file from PSPHOT.INPUT"));
@@ -165 +137 @@
         }
         psphotInput->src = psMemIncrRefCounter(output->fpa);
-        // specify the number of psphot input images
+        pmFPAAddSourceFromFormat(psphotInput->fpa, psphotInput->format); // ** builds the HDUs, is this OK?
+
+        // specify the number of psphot input images (psphotReadout loops over all input images)
         psMetadataAddS32 (config->arguments, PS_LIST_TAIL, "PSPHOT.INPUT.NUM", PS_META_REPLACE, "number of inputs", 1);
 
-        pmFPAfile *psphotInSources = pmFPAfileDefineSkycell(config, output->fpa, "PSPHOT.INPUT.CMF");
+        // the input sources (read from the input astrometry file) are transformed (in pswarpLoop) to the readout->analysis
+        // entries of the output file PSWARP.OUTPUT.SOURCES, associated with the main output pmFPAfile PSWARP.OUTPUT
+
+        // the PSPHOT.INPUT.CMF is used to supply those sources to psphot (specifically psphotLoadPSFSources).  
+
+        // pmFPAfile *psphotInSources = pmFPAfileDefineSkycell(config, output->fpa, "PSPHOT.INPUT.CMF");
+        pmFPAfile *psphotInSources = useInputCamera ? 
+            pmFPAfileDefineSkycell(config, output->fpa, "PSPHOT.INPUT.CMF"):
+            pmFPAfileDefineOutputForFormat(config, output->fpa, "PSPHOT.INPUT.CMF", skycell->cameraName, skycell->formatName);
         if (!psphotInSources) {
             psError(psErrorCodeLast(), false, _("Unable to generate output file from PSPHOT.INPUT.CMF"));
@@ -195 +177 @@
         psphotSources->save = false;
     }
+ 
+    // only keep the necessary (output) camera configs and relevant recipes
+    const char *skyCamera = psMetadataLookupStr(NULL, config->arguments, "SKYCELL.CAMERA");
+    pmConfigCamerasCull(config, skyCamera);
+    pmConfigRecipesCull(config, "PSWARP,PPSTATS,PSPHOT,PSASTRO,MASKS,JPEG");
+
+    psTrace("pswarp", 1, "Done with pswarpParseCamera...\n");
+    return true;
+}
+
+bool pswarpParseSingleInput (pmConfig *config) {
+
+    // The input image(s) is required: it defines the camera
+    pmFPAfile *input = pswarpDefineInputFile(config, NULL, "PSWARP.INPUT", "INPUT", PM_FPA_FILE_IMAGE);
+    if (!input) {
+        psError(psErrorCodeLast(), false, "Failed to build FPA from PSWARP.INPUT");
+        return false;
+    }
+
+    // Define the input astrometry file(s) (may be either WCS or CMF)
+    bool status = false;
+    pmFPAfile *astrom = pmFPAfileDefineFromArgs(&status, config, "PSWARP.ASTROM", "ASTROM");
+    if (!astrom) {
+        // look for the file on the RUN metadata
+        astrom = pmFPAfileDefineFromRun(&status, NULL, config, "PSWARP.ASTROM"); // File to return
+        if (!status) {
+            psLogMsg("pswarp", PS_LOG_INFO, "Failed to load file definition for %s", "PSWARP.ASTROM");
+        }
+    }
+    if (astrom) {
+        if ((astrom->type != PM_FPA_FILE_CMF) && (astrom->type != PM_FPA_FILE_WCS)) {
+            psLogMsg("pswarp", PS_LOG_INFO, "%s is neither CMF nor WCS", "PSWARP.ASTROM");
+            // XXX raise an error here??
+        } else {
+            foundAstrom = true;
+        }
+    }
+    psLogMsg("pswarp", PS_LOG_INFO, "Astrometry source: %s", astrom ? "supplied" : "header");
+
+    // Define the input mask file(s)
+    pmFPAfile *inMask = pswarpDefineInputFile(config, input, "PSWARP.MASK", "MASK", PM_FPA_FILE_MASK);
+    if (!inMask) {
+        psLogMsg("pswarp", PS_LOG_INFO, "No mask supplied");
+    }
+
+    // Define the input variance file(s)
+    pmFPAfile *inVariance = pswarpDefineInputFile(config, input, "PSWARP.VARIANCE", "VARIANCE", PM_FPA_FILE_VARIANCE);
+    if (inVariance) {
+        foundVariance = true;
+    } else {
+        psLogMsg("pswarp", PS_LOG_INFO, "No variance supplied");
+    }
+
+    pmFPAfile *inBackground = pswarpDefineInputFile(config, NULL, "PSWARP.BKGMODEL", "BACKGROUND", PM_FPA_FILE_IMAGE);
+    if (inBackground) {
+        foundBackground = true;
+    } else {
+      // cannot do the background model if an input is not supplied.
+      psMetadataAddBool (config->arguments, PS_LIST_TAIL, "BACKGROUND.MODEL", PS_META_REPLACE, "no input background supplied", false);
+      psLogMsg("pswarp", PS_LOG_INFO, "No background models supplied");
+    }
 
     // Chip selection: turn on only the chips specified
-    char *chipLine = psMetadataLookupStr(&mdok, config->arguments, "CHIP_SELECTIONS");
-    if (mdok) {
+    char *chipLine = psMetadataLookupStr(&status, config->arguments, "CHIP_SELECTIONS");
+    if (status) {
         psArray *chips = psStringSplitArray (chipLine, ",", false);
         if (chips->n > 0) {
@@ -213 +256 @@
     }
 
-    psTrace("pswarp", 1, "Done with pswarpParseCamera...\n");
+    psMetadataAddS32(config->arguments, PS_LIST_TAIL, "NUM_INPUTS", PS_META_REPLACE, "input file sets", 1);
     return true;
 }
+
+// read input file information from a metadata file
+// XXX for now, in the multi-input format, we disable PSF, BACKGROUND, and CHIP_SELECTIONS
+bool pswarpParseMultiInput (pmConfig *config, psMetadata *fileListMD) {
+
+    // the multi-input file consists of a set of metadata blocks, each with entries 
+    // for INPUT and (possibly) ASTROM, MASK, VARIANCE, ???
+
+    bool status = false;
+
+    for (int i = 0; i < fileListMD->list->n; i++) {
+        psMetadataItem *item = psMetadataGet(fileListMD, i);
+        if (item->type != PS_DATA_METADATA) {
+            psError(PS_ERR_BAD_PARAMETER_TYPE, true, "Component %s of the input metadata is not of type METADATA", item->name);
+            return false;
+        }
+
+        // next input metadata block of interest
+        psMetadata *fileBlockMD = item->data.md;
+
+        psString inputName = psMetadataLookupStr(&status, fileBlockMD, "INPUT"); // Name of image
+        if (!inputName || !strlen(inputName)) {
+            psError(PS_ERR_UNKNOWN, false, "input file not found in metadata block %d (%s)", i, item->name);
+            return false;
+        }
+        AddStringAsArray (config->arguments, inputName, "FILENAMES");
+
+        pmFPAfile *input = pswarpDefineInputFile (config, NULL, "PSWARP.INPUT", "FILENAMES", PM_FPA_FILE_IMAGE);
+        if (!input) {
+            psError(psErrorCodeLast(), false, "Failed to build FPA from PSWARP.INPUT");
+            return false;
+        }
+
+        psString astromName = psMetadataLookupStr(&status, fileBlockMD, "ASTROM"); // Name of astrom file
+        if (astromName && strlen(astromName)) {
+            AddStringAsArray (config->arguments, astromName, "FILENAMES");
+
+            pmFPAfile *astrom = pmFPAfileDefineFromArgs (&status, config, "PSWARP.ASTROM", "FILENAMES");
+            if (!status) {
+                psLogMsg("pswarp", PS_LOG_INFO, "Failed to load file definition for %s", "PSWARP.ASTROM");
+            }
+            if (astrom) {
+                if ((astrom->type != PM_FPA_FILE_CMF) && (astrom->type != PM_FPA_FILE_WCS)) {
+                    psLogMsg("pswarp", PS_LOG_INFO, "%s is neither CMF nor WCS", "PSWARP.ASTROM");
+                } else {
+                    foundAstrom = true;
+                }
+            }
+            psLogMsg("pswarp", PS_LOG_INFO, "Astrometry source: %s", astrom ? "supplied" : "header");
+        }
+
+        pmFPAfile *mask = NULL;
+        psString maskName = psMetadataLookupStr(&status, fileBlockMD, "MASK"); // Name of mask
+        if (maskName && strlen(maskName)) {
+            AddStringAsArray (config->arguments, maskName, "FILENAMES");
+
+            mask = pswarpDefineInputFile (config, input, "PSWARP.MASK", "FILENAMES", PM_FPA_FILE_MASK);
+        }
+        if (!mask) {
+            psLogMsg("pswarp", PS_LOG_INFO, "No mask supplied");
+        }
+
+        pmFPAfile *variance = NULL;
+        psString varianceName = psMetadataLookupStr(&status, fileBlockMD, "VARIANCE"); // Name of variance
+        if (varianceName && strlen(varianceName)) {
+            AddStringAsArray (config->arguments, varianceName, "FILENAMES");
+
+            variance = pswarpDefineInputFile (config, input, "PSWARP.VARIANCE", "FILENAMES", PM_FPA_FILE_VARIANCE);
+        }
+        if (variance) {
+            foundVariance = true;
+        } else {
+            psLogMsg("pswarp", PS_LOG_INFO, "No variance supplied");
+        }
+
+        // XXX for now we do not have an option to supply a background model set in multi input mode
+        psMetadataAddBool (config->arguments, PS_LIST_TAIL, "BACKGROUND.MODEL", PS_META_REPLACE, "no input background supplied", false);
+    }
+
+    psMetadataAddS32(config->arguments, PS_LIST_TAIL, "NUM_INPUTS", PS_META_REPLACE, "input file sets", fileListMD->list->n);
+    return true;
+}
+
+// Define an input file based on name in config->arguments or config RUN block
+pmFPAfile *pswarpDefineInputFile(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;
+
+    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) {
+        return NULL;
+    }
+
+    if (file->type != fileType) {
+        psError(PSWARP_ERR_CONFIG, true, "%s is not of type %s", filerule, pmFPAfileStringFromType(fileType));
+        return NULL;
+    }
+
+    return file;
+}
+
+bool AddStringAsArray (psMetadata *md, char *string, char *name) {
+
+    psArray *dummy = psArrayAlloc(1);   // Dummy array of filenames for this FPA
+    dummy->data[0] = psStringCopy(string);
+
+    psMetadataAddArray(md, PS_LIST_TAIL, name, PS_META_REPLACE, "string added as array", dummy);
+    psFree(dummy);
+    return true;
+}
+

trunk/pswarp/src/pswarpPixelsLit.c

@@ -64 +64 @@
     for (int y = 0; y < numRows; y++) {
         for (int x = 0; x < numCols; x++) {
-            if (!(mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & maskValue)) {
-                if (y > yMax) {
-                    yMax = y;
-                }
-                if (y < yMin) {
-                    yMin = y;
-                }
-                if (x > xMax) {
-                    xMax = x;
-                }
-                if (x < xMin) {
-                    xMin = x;
-                }
-            }
+            if (mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & maskValue) { continue; }
+            xMin = PS_MIN (xMin, x);
+            xMax = PS_MAX (xMax, x);
+            yMin = PS_MIN (yMin, y);
+            yMax = PS_MAX (yMax, y);
         }
     }
 
     if (stats) {
-        psMetadataAddS32(stats, PS_LIST_TAIL, "RANGE.XMIN", 0, "Minimum valid x value", xMin);
-        psMetadataAddS32(stats, PS_LIST_TAIL, "RANGE.XMAX", 0, "Maximum valid x value", xMax);
-        psMetadataAddS32(stats, PS_LIST_TAIL, "RANGE.YMIN", 0, "Minimum valid y value", yMin);
-        psMetadataAddS32(stats, PS_LIST_TAIL, "RANGE.YMAX", 0, "Maximum valid y value", yMax);
+        // XXX with multiple inputs (eg, output stacks -> exposure), these only represent the last input
+        psMetadataAddS32(stats, PS_LIST_TAIL, "RANGE.XMIN", PS_META_REPLACE, "Minimum valid x value", xMin);
+        psMetadataAddS32(stats, PS_LIST_TAIL, "RANGE.XMAX", PS_META_REPLACE, "Maximum valid x value", xMax);
+        psMetadataAddS32(stats, PS_LIST_TAIL, "RANGE.YMIN", PS_META_REPLACE, "Minimum valid y value", yMin);
+        psMetadataAddS32(stats, PS_LIST_TAIL, "RANGE.YMAX", PS_META_REPLACE, "Maximum valid y value", yMax);
     }
 

trunk/pswarp/src/pswarpTransformReadout.c

@@ -16 +16 @@
  * NOTE: in this function, the coordinates are transformed from the OUTPUT to the INPUT
  */
-bool pswarpTransformReadout(pmReadout *output, pmReadout *input, pmConfig *config)
+bool pswarpTransformReadout(pmReadout *output, pmReadout *input, pmConfig *config, bool backgroundWarp)
 {
     // XXX this implementation currently ignores the use of the region
@@ -52 +52 @@
     }
 
-    // XXX unused int nThreads = psMetadataLookupS32(&mdok, config->arguments, "NTHREADS"); // Number of threads
-    // XXX unused if (!mdok) {
-    // XXX unused     nThreads = 0;
-    // XXX unused }
     float poorFrac = psMetadataLookupF32(NULL, config->arguments, "POOR.FRAC"); ///< Flux fraction for "poor"
 
@@ -61 +57 @@
     // output coordinates to input coordinates
     pswarpMapGrid *grid = pswarpMapGridFromImage(input, output, nGridX, nGridY);
-    //    if (psMetadataLookupBool(NULL,config->arguments,"BACKGROUND_WARPING")) {
     
     // XXX optionally modify the grid based on this result and force the maxError < XXX
@@ -117 +112 @@
     psAssert (xGridMax < grid->nXpts, "xGridMax too big\n");
     psAssert (yGridMax < grid->nYpts, "yGridMax too big\n");
+
+    // fprintf (stderr, "warp %d,%d - %d,%d\n", xGridMin, yGridMin, xGridMax, yGridMax);
 
     // create jobs and supply them to the threads
@@ -132 +129 @@
             args->goodPixels = 0;
 
-            if (psMetadataLookupBool(NULL,config->arguments,"BACKGROUND_WARPING")) {
+            if (backgroundWarp) {
               args->background_warping = true;
               args->offset_x = psMetadataLookupF32(NULL,config->arguments,"BKG_WARP_XOFFSET");
               args->offset_y = psMetadataLookupF32(NULL,config->arguments,"BKG_WARP_YOFFSET");
             }
-
             
             // allocate a job
@@ -210 +206 @@
     psFree(interp);
 
-    if (goodPixels > 0 && psMetadataLookupBool(&mdok, recipe, "SOURCES")) {
-      if (!psMetadataLookupBool(NULL,config->arguments,"BACKGROUND_WARPING")) {
-        if (!pswarpTransformSources(output, input, config)) {
-          psError(psErrorCodeLast(), false, "Unable to interpolate image.");
-          return false;
+    if (goodPixels > 0 && !backgroundWarp && psMetadataLookupBool(&mdok, recipe, "SOURCES")) {
+        if (!pswarpTransformSources(output, input, config)) {
+            psError(psErrorCodeLast(), false, "Unable to transform sources.");
+            return false;
         }
-      }
     }
 

trunk/pswarp/src/pswarpTransformTile.c

@@ -89 +89 @@
 
     for (int y = yMin; y < yMax; y++) {
+
+      int yOut = y - outRow0; ///< Position on output image
+
         for (int x = xMin; x < xMax; x++) {
             // Only transform those pixels requested
@@ -94 +97 @@
                 continue;
             }
+
+            int xOut = x - outCol0;
+
+            // XXX the existing image may already have valid data -- probably should keep 
+            // the best, but for the moment, just keep the pixel which is not NAN
+            if (isfinite(outImageData[yOut][xOut])) continue;
+
             // pswarpMapApply converts the output coordinate (x,y) to the input coordinate.
             // both are in the parent frames of the input and output images.
@@ -127 +137 @@
             }
 
-            int xOut = x - outCol0, yOut = y - outRow0; ///< Position on output image
-
             if (outImageData) {
-                outImageData[yOut][xOut] = imageValue * jacobian;
+              // XXX TEST outImageData[yOut][xOut] = value;
+              outImageData[yOut][xOut] = imageValue * jacobian;
             }
             if (outVarData) {

trunk/pswarp/src/pswarpVersion.c

@@ -13 +13 @@
 #include <config.h>
 #endif
-
-#include <stdio.h>
-#include <pslib.h>
-#include <psmodules.h>
-#include <psphot.h>
-#include <ppStats.h>
 
 #include "pswarp.h"