Changeset 35535

Timestamp:

May 8, 2013, 4:41:42 PM (13 years ago)

Author:

eugene

Message:

warping of background model appears to work

Location:

branches/eam_branches/ipp-20130419/pswarp/src

Files:

• Makefile.am (modified) (1 diff)
• pswarp.c (modified) (1 diff)
• pswarp.h (modified) (2 diffs)
• pswarpArguments.c (modified) (4 diffs)
• pswarpDefineLayout.c (modified) (6 diffs)
• pswarpLoadAstrometry.c (modified) (3 diffs)
• pswarpLoop.c (modified) (9 diffs)
• pswarpLoopBackground.c (modified) (3 diffs)
• pswarpMakePSF.c (modified) (4 diffs)
• pswarpMatchRange.c (modified) (1 diff)
• pswarpOptions.c (added)
• pswarpOverlaps.c (modified) (2 diffs)
• pswarpParseCamera.c (modified) (3 diffs)
• pswarpTransformReadout.c (modified) (4 diffs)
• pswarpUpdateMetadata.c (modified) (4 diffs)
• pswarpUpdateStatistics.c (modified) (1 diff)

branches/eam_branches/ipp-20130419/pswarp/src/Makefile.am

@@ -44 +44 @@
         pswarpMatchRange.c              \
         pswarpOverlaps.c                \
+        pswarpOptions.c                 \
         pswarpParseCamera.c             \
         pswarpPixelsLit.c               \

branches/eam_branches/ipp-20130419/pswarp/src/pswarp.c

@@ -39 +39 @@
     }
 
-# if (1)
-    // XXX I'm going to add the new code before removing the old code.
     // load the input & output astrometry, find the output overlaps, generate the output pixels
     if (!pswarpDefineLayout(config)) {
         pswarpCleanup(config, statsFile);
     }
-    // load and warp
+
+    // load input pixels & and warp
     if (!pswarpLoop(config, statsFile->md)) {
         pswarpCleanup(config, statsFile);
     }
 
-# else
-    // load the skycell layout information:
-    //   read the output astrometry description
-    //   generate the output pixels
-    //   convert header astrometry to pmFPA astrometry structures
-    if (!pswarpDefine(config)) {
-        pswarpCleanup(config, statsFile);
+    // load input pixels & and warp
+    if (!pswarpLoopBackground(config, statsFile->md)) {
+        pswarpCleanup(config, statsFile);
     }
-    if (!pswarpDefineBackground(config)) {
-        pswarpCleanup(config, statsFile);
-    }
-    // load and warp
-    if (!pswarpLoopSkycell(config, statsFile->md)) {
-        pswarpCleanup(config, statsFile);
-    }
-    if (!pswarpLoopBackground(config, statsFile->md)) {
-        fprintf(stderr,"Dying!\n");
-        pswarpCleanup(config, statsFile);
-    }
-# endif
 
+    // output and free 
+    // NOTE: pswarpCleanup calls exit
     psLogMsg("pswarp", PS_LOG_INFO, "complete pswarp run: %f sec\n", psTimerMark("pswarp"));
     pswarpCleanup(config, statsFile);
 }
-
-// NOTE: pswarpCleanup calls exit
-
-// pswarpDumpOutput (config);
-// exit (0);

branches/eam_branches/ipp-20130419/pswarp/src/pswarp.h

@@ -201 +201 @@
 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 *input, pmFPA *astrom, pmConfig *config);
+bool pswarpUpdateMetadata (pmFPA *output, pmFPA *skycell, pmFPA *input, pmFPA *astrom, pmConfig *config, bool fullImage);
 
 // XXX functions in pswarpParseCamera
@@ -216 +216 @@
 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);

branches/eam_branches/ipp-20130419/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);
     }
 
@@ -143 +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 pswarpOptions...\n");
-
-    return (config);
-}

branches/eam_branches/ipp-20130419/pswarp/src/pswarpDefineLayout.c

@@ -11 +11 @@
 # include "pswarp.h"
 
+bool pswarpGenerateOutputCell (pmChip *chip, pmCell *cell, pmCell *refcell, float xBin, float yBin);
+
 // load the astrometry header info and generate the tranformations 
 bool pswarpDefineLayout (pmConfig *config) {
@@ -31 +33 @@
         return false;
     }
+
+    // if a background model is supplied and the output is requested, then 
+    // this will exist.  if not, skip related features 
+    pmFPAfile *bkgModel = psMetadataLookupPtr(&status, config->files, "PSWARP.OUTPUT.BKGMODEL");
 
     // chips are not processed unless we have determined they overlap the inputs
@@ -54 +60 @@
             astrom = input;
         }
+        if (astrom->camera != input->camera) {
+            // XXX this may cause an error for output / skycell
+            psError(PSWARP_ERR_DATA, true, "Input camera and astrometry camera do not match.");
+            return false;
+        }
         if (!pswarpLoadAstrometry (input, astrom, config)) {
             psError(PSWARP_ERR_CONFIG, false, "problem loading input astrometry\n");
@@ -81 +92 @@
     pmFPAview *view = pmFPAviewAlloc(0);
 
+    // check input astrometry mode
+    bool bilevelAstrometry = psMetadataLookupBool (NULL, skycell->fpa->analysis, "ASTROMETRY.BILEVEL");
+    
     // generate the output hdu and contents bits
     pmFPAAddSourceFromFormat(output->fpa, output->format);
+
+    // XXX this is a bit muddled : we are setting up the astrometry elements for the output
+    // background model here instead of in pswarpLoadAstrometry, where it would be more natural
+    // the only explanation I have is that 1) I need to have the binning factors for the input image, but 
+    // 2) at this point I don't yet have the input background models loaded and 3) the source
+    // for the binning factor of the input and output background models is different
+    // (header->CCDSUMi for input, recipe / pmFPAfile.xbin.ybin for the other)
+
+    // if we have an output background model, generate hdu, etc and supply astrometry
+    if (bkgModel) {
+        pmFPAAddSourceFromFormat(bkgModel->fpa, bkgModel->format);
+
+        if (bilevelAstrometry) {
+            // top-level elements of the bkgModel astrometry
+            bkgModel->fpa->toTPA   = psMemIncrRefCounter (skycell->fpa->toTPA);
+            bkgModel->fpa->fromTPA = psMemIncrRefCounter (skycell->fpa->fromTPA);
+            bkgModel->fpa->toSky   = psMemIncrRefCounter (skycell->fpa->toSky);
+        }
+    }
 
     pmChip *chip;
@@ -95 +128 @@
             pmCell *refcell = pmFPAviewThisCell(view, skycell->fpa); ///< Target cell
 
-            // we've got the output astrom header
-            pmHDU *hdu = pmHDUFromCell(refcell); ///< HDU for source
-            if (!hdu || !hdu->header) {
-                psError(PM_ERR_PROG, false, "Unable to find header for sky cell.");
+            if (!pswarpGenerateOutputCell (chip, cell, refcell, output->xBin, output->yBin)){
+                psError(PSWARP_ERR_DATA, false, "failed to generate output cell");
                 psFree(view);
                 return false;
             }
-            int numCols = psMetadataLookupS32(&status, hdu->header, "NAXIS1"); ///< Number of columns
-            int numRows = psMetadataLookupS32(&status, 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;
+
+            if (bkgModel) {
+                pmChip *bkgChip = pmFPAviewThisChip(view, bkgModel->fpa); ///< Target cell
+                pmCell *bkgCell = pmFPAviewThisCell(view, bkgModel->fpa); ///< Target cell
+                if (!pswarpGenerateOutputCell (bkgChip, bkgCell, refcell, bkgModel->xBin, bkgModel->yBin)) {
+                    psError(PSWARP_ERR_DATA, false, "failed to generate output cell");
+                    psFree(view);
+                    return false;
+                }
+                if (!pswarpModifyChipAstrom (config, view, bkgChip, skycell, bilevelAstrometry, bkgModel->xBin, bkgModel->yBin)) {
+                    psError(PSWARP_ERR_DATA, false, "problem with output astrometry");
+                    psFree(view);
+                    return false;
+                }
             }
-                
-            // generate the pixels for output->fpa
-            pmReadout *readout = pmReadoutAlloc(cell); ///< output readout
-            readout->image = psImageAlloc(numCols, numRows, PS_TYPE_F32);
-
-            psImageInit(readout->image, NAN);
-            psFree(readout);                // Drop reference
-                
-            // copy the image concepts from the skycell 
-            bool status = false;
-            psMetadataItemSupplement(&status, cell->concepts, refcell->concepts, "CELL.XBIN");
-            psMetadataItemSupplement(&status, cell->concepts, refcell->concepts, "CELL.YBIN");
-            psMetadataItemSupplement(&status, cell->concepts, refcell->concepts, "CELL.XSIZE");
-            psMetadataItemSupplement(&status, cell->concepts, refcell->concepts, "CELL.YSIZE");
-            psMetadataItemSupplement(&status, cell->concepts, refcell->concepts, "CELL.XPARITY");
-            psMetadataItemSupplement(&status, cell->concepts, refcell->concepts, "CELL.YPARITY");
-            psMetadataItemSupplement(&status, cell->concepts, refcell->concepts, "CELL.X0");
-            psMetadataItemSupplement(&status, cell->concepts, refcell->concepts, "CELL.Y0");
-
-            chip->data_exists = true;
-            chip->file_exists = true;
-            cell->data_exists = true;
-            cell->file_exists = true;
-
-            // copy the basic headers across from astrom ref to output
-            pmHDU *outHDU = pmHDUFromCell (cell);           ///< HDU for the output warped image
-            outHDU->header = psMetadataCopy(outHDU->header, hdu->header);
-            pswarpVersionHeader(outHDU->header);
         }
     }
@@ -143 +154 @@
 }
 
-// Lists of file rules used at different stages
-# include "pswarpFileNames.h"
-
-static int Nout = 0;
-
-// XX function for tests
-bool pswarpDumpOutput (pmConfig *config) {
-
-    pmFPAfile *output = psMetadataLookupPtr(NULL, config->files, "PSWARP.OUTPUT");
-    if (!output) {
-        psError(PSWARP_ERR_CONFIG, false, "Can't find output data!\n");
-        return false;
-    }
-
-    pswarpFileActivation(config, detectorFiles, false);
-    pswarpFileActivation(config, photFiles, false);
-    pswarpFileActivation(config, independentFiles, false);
-    pswarpFileActivation(config, skycellFiles, false);
-    pmFPAfileActivate(config->files, true, "PSWARP.OUTPUT");
-
-    pmFPAview *view = pmFPAviewAlloc(0);
-
-    pmChip *chip;
-    while ((chip = pmFPAviewNextChip (view, output->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; }
-
-        pmCell *cell;
-        while ((cell = pmFPAviewNextCell (view, output->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; }
-
-            // process each of the readouts
-            pmReadout *readout;
-            while ((readout = pmFPAviewNextReadout(view, output->fpa, 1)) != NULL) {
-                if (!readout->data_exists) {
-                    continue;
-                }
-
-                char name[64];
-                snprintf (name, 64, "dumpwarp.%02d.fits", Nout);
-                Nout ++;
+bool pswarpGenerateOutputCell (pmChip *chip, pmCell *cell, pmCell *refcell, float xBin, float yBin) {
+
+    bool status = false;
+
+    // we've got the output astrom header
+    pmHDU *hdu = pmHDUFromCell(refcell); ///< HDU for source
+    if (!hdu || !hdu->header) {
+        psError(PM_ERR_PROG, false, "Unable to find header for sky cell.");
+        return false;
+    }
+    int numCols = psMetadataLookupS32(&status, hdu->header, "NAXIS1") / xBin; ///< Number of output columns
+    int numRows = psMetadataLookupS32(&status, hdu->header, "NAXIS2") / yBin; ///< Number of output rows
+    if ((numCols == 0) || (numRows == 0)) {
+        psError(PSWARP_ERR_DATA, false, "output cell has invalid dimensions %d x %d", numCols, numRows);
+        return false;
+    }
                 
-                psphotSaveImage (NULL, readout->image, name);
-            }
-        }
-        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);
-
-DONE:
-
-    pswarpFileActivation(config, detectorFiles, true);
-    pmFPAfileActivate(config->files, false, "PSWARP.ASTROM");
+    // generate the pixels for output->fpa
+    pmReadout *readout = pmReadoutAlloc(cell); ///< output readout
+    readout->image = psImageAlloc(numCols, numRows, PS_TYPE_F32);
+
+    psImageInit(readout->image, NAN);
+    psFree(readout);                // Drop reference
+                
+    // copy the image concepts from the skycell 
+    psMetadataItemSupplement(&status, cell->concepts, refcell->concepts, "CELL.XBIN");
+    psMetadataItemSupplement(&status, cell->concepts, refcell->concepts, "CELL.YBIN");
+    psMetadataItemSupplement(&status, cell->concepts, refcell->concepts, "CELL.XSIZE");
+    psMetadataItemSupplement(&status, cell->concepts, refcell->concepts, "CELL.YSIZE");
+    psMetadataItemSupplement(&status, cell->concepts, refcell->concepts, "CELL.XPARITY");
+    psMetadataItemSupplement(&status, cell->concepts, refcell->concepts, "CELL.YPARITY");
+    psMetadataItemSupplement(&status, cell->concepts, refcell->concepts, "CELL.X0");
+    psMetadataItemSupplement(&status, cell->concepts, refcell->concepts, "CELL.Y0");
+
+    chip->data_exists = true;
+    chip->file_exists = true;
+    chip->process     = true;
+
+    cell->data_exists = true;
+    cell->file_exists = true;
+    cell->process     = true;
+
+    // copy the basic headers across from astrom ref to output
+    pmHDU *outHDU = pmHDUFromCell (cell);           ///< HDU for the output warped image
+    outHDU->header = psMetadataCopy(outHDU->header, hdu->header);
+    pswarpVersionHeader(outHDU->header);
     return true;
-}    
-
+}

branches/eam_branches/ipp-20130419/pswarp/src/pswarpLoadAstrometry.c

@@ -25 +25 @@
     astrom->type = PM_FPA_FILE_WCS;
 
-    // Read the pmFPA headers
+    // **** read in all of the headers from the astrometry source
+
     pmFPAview *view = pmFPAviewAlloc(0);
 
@@ -76 +77 @@
     }
 
-    // the section below converts the header astrometry info (in astrom) to the pmFPA
-    // astrometry structures, saving them on the target fpa structure.
+    // **** convert the header astrometry info (in astrom) to the pmFPA astrometry
+    // structures, saving them on the target fpa structure.
 
     pmFPAviewReset (view);
@@ -97 +98 @@
         }
     }
+    psMetadataAddBool (astrom->fpa->analysis, PS_LIST_TAIL, "ASTROMETRY.BILEVEL", PS_META_REPLACE, "bilevel astrometry?", bilevelAstrometry);
+
     // for pswarpLocalFrame, I need transformations and HDUs on a single fpa (so set on astrom as well as target)
     // target->fpa->toTPA = psMemIncrRefCounter (astrom->fpa->toTPA);

branches/eam_branches/ipp-20130419/pswarp/src/pswarpLoop.c

@@ -1 +1 @@
 /** @file pswarpLoop.c
  *
- *  @brief mail processing loop for pswarp
+ *  @brief main processing loop for pswarp
  *  @ingroup pswarp
  *
@@ -13 +13 @@
 #include "pswarpFileNames.h"            // Lists of file rules used at different stages
 
-// 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)
 {
-    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
@@ -26 +22 @@
     if (!recipe) {
         psError(PSWARP_ERR_CONFIG, false, "missing recipe %s", PSWARP_RECIPE);
-        return false;
-    }
-
-    // XXX move this to pswarpParseCamera?  has to come after input mask header is read
-    if (!pswarpSetMaskBits(config)) {
-        psError(psErrorCodeLast(), false, "failed to set mask bits");
         return false;
     }
@@ -66 +56 @@
         }
 
-        fprintf (stderr, "loading file %s\n", input->filename);
-
-        // use the external astrometry source if supplied
+        // select the input data sources
         pmFPAfile *astrom = pmFPAfileSelectSingle(config->files, "PSWARP.ASTROM", i);
         if (!astrom) {
@@ -74 +62 @@
         }
 
-        if (astrom->camera != input->camera) {
-            psError(PSWARP_ERR_DATA, true, "Input camera and astrometry camera do not match.");
-            return false;
-        }
-
-        // ensure everyone is off except what we need below
-        pswarpFileActivation(config, detectorFiles, false);
+        // ensure everyone is off except what we need below 
+        // XXX double-check activation issues in general
+        pswarpFileActivation(config, detectorFiles, true);
         pswarpFileActivation(config, photFiles, false);
         pswarpFileActivation(config, independentFiles, false);
         pswarpFileActivation(config, skycellFiles, false);
     
-        // the loops below load the input data
-        pswarpFileActivation(config, detectorFiles, true);
-
-        // need to read in these (but only if CMF, right?)
+        // 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");
-
-        // pswarpFileActivation(config, independentFiles, true);
 
         pmFPAviewReset (view);
@@ -174 +154 @@
             goto FAIL;
         }
-        if (!pswarpUpdateMetadata (output->fpa, input->fpa, astrom->fpa, config)) {
+        if (!pswarpUpdateMetadata (output->fpa, skycell->fpa, input->fpa, astrom->fpa, config, true)) {
             psError(psErrorCodeLast(), false, "problem generating statistics.");
             goto FAIL;
@@ -184 +164 @@
       goto FAIL;
     }
-
-    // XXX
-    // pswarpDumpOutput (config);
-    // exit (0);
 
     psFree(view);
@@ -196 +172 @@
     return false;
 }
-
-// static int Nout = 0;
 
 // once the output fpa elements have been built, loop over the fpa and generate stats
@@ -219 +193 @@
             while ((readout = pmFPAviewNextReadout(view, output, 1)) != NULL) {
                 pswarpTransformReadout (readout, input, config);
-                // char name[64];
-                // snprintf (name, 64, "dumpwarp.%02d.fits", Nout);
-                // Nout ++;
-                // 
-                // psphotSaveImage (NULL, readout->image, name);
             }
         }

branches/eam_branches/ipp-20130419/pswarp/src/pswarpLoopBackground.c

@@ -11 +11 @@
 
 #include "pswarp.h"
-#include <ppStats.h>
 #include "pswarpFileNames.h"            // Lists of file rules used at different stages
 
@@ -22 +21 @@
     bool status;
     bool mdok;                          // Status of MD lookup
-    const char *skyCamera = psMetadataLookupStr(NULL, config->arguments,
-                                                "SKYCELL.CAMERA");  // Name of camera for skycell
+
     if (!psMetadataLookupBool(NULL,config->arguments,"BACKGROUND.MODEL")) return true;
 
-    pmConfigCamerasCull(config, skyCamera);
-    pmConfigRecipesCull(config, "PSWARP,PPSTATS,PSPHOT,PSASTRO,MASKS,JPEG");
+    int xGrid = psMetadataLookupS32(&mdok,config->arguments,"BKG.XGRID");
+    int yGrid = psMetadataLookupS32(&mdok,config->arguments,"BKG.YGRID");
+
+    // WHAT IS THIS?? should it be disabled at the end?
+    psMetadataAddS32(config->arguments, PS_LIST_TAIL, "INTERPOLATION.MODE", PS_META_REPLACE, "", 8);
 
     // load the recipe
@@ -36 +37 @@
     }
 
-    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");
+    pmFPAfile *output = psMetadataLookupPtr(&status, config->files, "PSWARP.OUTPUT.BKGMODEL");
+    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.");
-        return false;
-    }
-
-    // select the output readout
+    pmFPAfile *skycell = psMetadataLookupPtr(&status, config->files, "PSWARP.SKYCELL");
+    if (!skycell) {
+        psError(PSWARP_ERR_DATA, true, "Cannot find output astrometry.");
+        return false;
+    }
+
     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;
-    }
+
+    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;
+    }
+
+    // loop over this section once per input group
+    for (int i = 0; i < nInputs; i++) {
+
+        // select the input data source : we are reading & transforming pixels for the background model
+        pmFPAfile *input = pmFPAfileSelectSingle(config->files, "PSWARP.BKGMODEL", 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;
+        }
+
+        // check input astrometry mode
+        bool bilevelAstrometry = psMetadataLookupBool (NULL, astrom->fpa->analysis, "ASTROMETRY.BILEVEL");
+        if (bilevelAstrometry) {
+            // top-level elements of the input model astrometry
+            input->fpa->toTPA = psMemIncrRefCounter (astrom->fpa->toTPA);
+            input->fpa->fromTPA = psMemIncrRefCounter (astrom->fpa->fromTPA);
+            input->fpa->toSky = psMemIncrRefCounter (astrom->fpa->toSky);
+        }
+
+        // we only need to load the background models here 
+        pswarpFileActivation(config, detectorFiles, false);
+        pswarpFileActivation(config, photFiles, false);
+        pswarpFileActivation(config, independentFiles, false);
+        pswarpFileActivation(config, skycellFiles, false);
+
+        // we are only reading the background models, not writing anything
+        pmFPAfileActivate(config->files, true, "PSWARP.BKGMODEL");
+
+        pmFPAviewReset (view);
+
+        pmFPAview *view = pmFPAviewAlloc(0);
+        if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) {
+            psError(psErrorCodeLast(), false, "Unable to read files.");
+            goto FAIL;
+        }
+
+        pmChip *chip;
+        while ((chip = pmFPAviewNextChip (view, input->fpa, 1)) != NULL) {
+            if (!chip->process || !chip->file_exists) { continue; }
+            if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) {
+                psError(psErrorCodeLast(), false, "Unable to read files.");
+                goto FAIL;
+            }
+
+            // adjust output astrometry
+            double xBin = NAN;
+            double yBin = NAN;
+            pswarpGetInputScales (&xBin, &yBin, config, view, chip);
+            double fluxRatio = xBin * yBin / (float) (xGrid * yGrid);
+            
+            if (!pswarpModifyChipAstrom (config, view, chip, astrom, bilevelAstrometry, xBin, yBin)) {
+                psError(psErrorCodeLast(), false, "failed to set BKGMODEL astrometry.");
+                goto FAIL;
+            }
+
+            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)) {
+                    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;
+                    }
+
+                    // re-normalize the BKGMODEL pixels by modified astrometry
+                    for (int x = 0; x < readout->image->numCols; x++) {
+                        for (int y = 0; y < readout->image->numRows; y++) {
+                            readout->image->data.F32[y][x] *= fluxRatio;
+                        }
+                    }
+                
+                    // transform the actual BKGMODEL pixels from readout to output
+                    // XXX replace the use of this temporary metadata item with an argument?
+                    // XXX move the source transform out of pswarpTransformReadout?
+                    psMetadataAddBool(config->arguments,PS_LIST_TAIL, "BACKGROUND_WARPING", PS_META_REPLACE, "", true);
+                    pswarpTransformToTarget (output->fpa, readout, config);
+                    psMetadataAddBool(config->arguments,PS_LIST_TAIL, "BACKGROUND_WARPING", PS_META_REPLACE, "", false);
+
+                    if (!pmFPAfileIOChecks (config, view, PM_FPA_AFTER)) { // Readout
+                        psError(psErrorCodeLast(), false, "Unable to close files.");
+                        goto FAIL;
+                    }
+                }
+
+                if (!pmFPAfileIOChecks (config, view, PM_FPA_AFTER)) { // Cell
+                    psError(psErrorCodeLast(), false, "Unable to close files.");
+                    goto FAIL;
+                }
+            }
+
+            if (!pmFPAfileIOChecks (config, view, PM_FPA_AFTER)) { // Chip
+                psError(psErrorCodeLast(), false, "Unable to close files.");
+                goto FAIL;
+            }
+        }
+
+        if (!pmFPAfileIOChecks (config, view, PM_FPA_AFTER)) { // FPA
+            psError(psErrorCodeLast(), false, "Unable to close files.");
+            goto FAIL;
+        }
+
+        if (false && !pswarpUpdateMetadata (output->fpa, skycell->fpa, input->fpa, astrom->fpa, config, false)) {
+            psError(psErrorCodeLast(), false, "problem generating statistics.");
+            goto FAIL;
+        }
+    }
+
     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;
+    return true;
+
+FAIL:
+    psFree (view);
+    return false;
+}
+
+// **** input == (small) BKGMODEL, astrom == (fullsize) ASTROM
+// use astrometry from input ref astrometry and BINNING info from input model to modify the astrometry 
+// elements of input chip.  scale is the area of an input pixel is raw image pixels
+bool pswarpGetInputScales (double *xBin, double *yBin, pmConfig *config, pmFPAview *view, pmChip *chip) {
+
+    // Pull information from the header of the background files so we can use it to set values.
+    pmHDU *hdu = pmFPAviewThisHDU(view, chip->parent);
+    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," "));
+    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:
+    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"));
+
+    *xBin = CCDSUM1;
+    *yBin = CCDSUM2;
 
     return true;
 }
+
+// **** input == (small) BKGMODEL, astrom == (fullsize) ASTROM
+// use astrometry from input ref astrometry and BINNING info from input model to modify the astrometry 
+// elements of input chip.  scale is the area of an input pixel is raw image pixels
+bool pswarpModifyChipAstrom (pmConfig *config, pmFPAview *view, pmChip *chip, pmFPAfile *astrom, bool bilevelAstrometry, double xBin, double yBin) {
+
+    // read WCS data from the corresponding header 
+    pmHDU *hdu = pmFPAviewThisHDU (view, astrom->fpa);
+
+    // generate a WCS structure from the ASTROM header keywords
+    pmAstromWCS *WCS = pmAstromWCSfromHeader(hdu->header);
+
+    // re-scale the terms of the WCS to match the BKGMODEL scale
+    double cd1f = xBin;
+    double cd2f = yBin;
+
+    WCS->cdelt1 *= cd1f;
+    WCS->cdelt2 *= cd2f;
+    WCS->crpix1 = WCS->crpix1 / cd1f;
+    WCS->crpix2 = WCS->crpix2 / cd2f;
+
+    // WCS->trans->x->nX/nY
+    // adjust the scale of the WCS coeffs
+    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);
+        }
+    }
+
+    // this is a bit crude : use a locally temporary header to store the WCS information
+    // which is then applied to the input chip (why not modify the input chip header?)
+    psMetadata *tmpHeader = psMetadataAlloc();
+
+    // write the modified WCS to the header
+    // XXX this should probably not modify the astrom header
+    pmAstromWCStoHeader (tmpHeader, WCS);
+        
+    // use the modified header to update the WCS elements of the BKGMODEL fpa/chip structures
+    if (bilevelAstrometry) {
+        if (!pmAstromReadBilevelChip (chip, tmpHeader)) {
+            psError(psErrorCodeLast(), false, "Unable to read bilevel chip astrometry for input FPA.");
+            return false;
+        }
+    } else {
+        // we use a default FPA pixel scale of 1.0
+        if (!pmAstromReadWCS (chip->parent, chip, tmpHeader, 1.0)) {
+            psError(psErrorCodeLast(), false, "Unable to read WCS astrometry for input FPA.");
+            return false;
+        }
+    }
+
+    return true;
+}

branches/eam_branches/ipp-20130419/pswarp/src/pswarpMakePSF.c

@@ -11 +11 @@
 
 #include "pswarp.h"
-
-#define WCS_NONLIN_TOL 0.001            // Non-linear tolerance for header WCS
-#define TESTING 0                       // Testing output?
 
 // We need a new PSF model for the warped frame.  It would be good to generate this analytically, but
@@ -34 +31 @@
     }
 
-    // XXX What was being read here??  skip this for now
-    if (false && !pswarpIOChecksBefore(config)) {
-        psError(psErrorCodeLast(), false, "Unable to read files.");
-        return false;
-    }
-
+    // XXX move above the loop
+    pmModelClassSetLimits(PM_MODEL_LIMITS_STRICT);
+                
     // supply the readout and fpa of interest to psphot
     pmFPAfile *photFile = psMetadataLookupPtr(&status, config->files, "PSPHOT.INPUT");
@@ -71 +65 @@
                 }
 
-                // XXX move above the loop
-                pmModelClassSetLimits(PM_MODEL_LIMITS_STRICT);
-                
                 // measure the PSF using these sources
                 if (!psphotReadoutFindPSF(config, view, "PSPHOT.INPUT", sources)) {
@@ -96 +87 @@
         item->data.F32 = psMetadataLookupF32(&status, chip->concepts, "CHIP.SEEING");
     }
-
-    // XXX test function found below
-    // pswarpMakePSF_test ();
 
     psFree(view);

branches/eam_branches/ipp-20130419/pswarp/src/pswarpMatchRange.c

@@ -98 +98 @@
 
     // demo forward and backward transformation
-# if (1)
+# if (0)
     srcPix->x = *minX;
     srcPix->y = *minY;

branches/eam_branches/ipp-20130419/pswarp/src/pswarpOverlaps.c

@@ -83 +83 @@
             }
         }
-        // fprintf (stderr, "bounds %d : %f - %f, %f - %f\n", view->chip, Pmin, Pmax, Qmin, Qmax);
         pswarpBoundsAppend (bounds, Pmin, Pmax, Qmin, Qmax);
     }
@@ -215 +214 @@
                 pmChip *inputChip = inputChips->data[nChip];
                 for (int nCell = 0; nCell < inputChip->cells->n; nCell++) {
-                    fprintf (stderr, "input %d, %d : output %d\n", nChip, nCell, j);
                     psArrayAdd (inputCells, 8, inputChip->cells->data[nCell]);
                 }

branches/eam_branches/ipp-20130419/pswarp/src/pswarpParseCamera.c

@@ -38 +38 @@
     }
 
+    // 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])
 
@@ -95 +101 @@
     // 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");
-      outBackground->xBin = psMetadataLookupS32(&status, recipe, "BKG.XGRID");
-      outBackground->yBin = psMetadataLookupS32(&status, recipe, "BKG.YGRID");
-      
-      if (!outBackground) {
-        psError(psErrorCodeLast(), false, "Failed to build FPA from PSWARP.OUTPUT.BKGMODEL");
-        return false;
-      }
-      outBackground->save = true;
+        // pmFPAfile *outBackground = pmFPAfileDefineSkycell(config, NULL, "PSWARP.OUTPUT.BKGMODEL");
+        pmFPAfile *outBackground = pmFPAfileDefineOutputForFormat(config, NULL, "PSWARP.OUTPUT.BKGMODEL", skycell->cameraName, skycell->formatName);
+        outBackground->xBin = psMetadataLookupS32(&status, recipe, "BKG.XGRID");
+        outBackground->yBin = psMetadataLookupS32(&status, recipe, "BKG.YGRID");
+        
+        if (!outBackground) {
+            psError(psErrorCodeLast(), false, "Failed to build FPA from PSWARP.OUTPUT.BKGMODEL");
+            return false;
+        }
+        outBackground->save = true;
     }
 
@@ -152 +159 @@
     }
  
+    // 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;

branches/eam_branches/ipp-20130419/pswarp/src/pswarpTransformReadout.c

@@ -48 +48 @@
     psAssert(mdok, "MASK.INPUT was not defined");
 
+    bool backgroundWarp = psMetadataLookupBool(NULL,config->arguments,"BACKGROUND_WARPING");
+
     if (!input->covariance) {
         input->covariance = psImageCovarianceNone();
     }
 
-    // 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 +59 @@
     // 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
@@ -118 +115 @@
     psAssert (yGridMax < grid->nYpts, "yGridMax too big\n");
 
-    fprintf (stderr, "warp %d,%d - %d,%d\n", xGridMin, yGridMin, xGridMax, yGridMax);
+    // fprintf (stderr, "warp %d,%d - %d,%d\n", xGridMin, yGridMin, xGridMax, yGridMax);
 
     // create jobs and supply them to the threads
@@ -134 +131 @@
             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

branches/eam_branches/ipp-20130419/pswarp/src/pswarpUpdateMetadata.c

@@ -10 +10 @@
 
 #include "pswarp.h"
+#define WCS_NONLIN_TOL 0.001            // Non-linear tolerance for header WCS
 
-// once the output fpa elements have been built, loop over the fpa and generate stats
-// for each readout
-bool pswarpUpdateMetadata (pmFPA *output, pmFPA *input, pmFPA *astrom, pmConfig *config)  {
+// once the output fpa elements have been built, loop over the fpa and generate metadata
+bool pswarpUpdateMetadata (pmFPA *output, pmFPA *skycell, pmFPA *input, pmFPA *astrom, pmConfig *config, bool fullImage)  {
 
     pmFPAview *view = pmFPAviewAlloc(0);
     
+    // check output astrometry mode
+    bool bilevelAstrometry = psMetadataLookupBool (NULL, skycell->analysis, "ASTROMETRY.BILEVEL");
+
     pmChip *chip;
     while ((chip = pmFPAviewNextChip (view, output, 1)) != NULL) {
@@ -39 +42 @@
                 if (!readout->data_exists) continue;
     
-                // Set covariance matrix for output
-# if (0)
-                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);
+                if (fullImage) {
+                    // Set covariance matrix for output
+                    bool status = false;
+                    psList *covariances = psMetadataLookupPtr(&status, readout->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(readout->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);
+                    // Correct covariance matrix scale for the mean (square root of the) Jacobian
+                    double jacobian = psMetadataLookupF64(NULL, readout->analysis, PSWARP_ANALYSIS_JACOBIAN); // Jacobian
+                    int goodPixels = psMetadataLookupS32(NULL, readout->analysis, PSWARP_ANALYSIS_GOODPIX);   // Good pixels
+                    jacobian /= goodPixels;
+                    readout->covariance = psImageCovarianceScale(covar, jacobian);
+                    psFree(covar);
 
-                if (output->variance) {
-                  psImageCovarianceTransfer(output->variance, output->covariance);
+                    if (readout->variance) {
+                        psImageCovarianceTransfer(readout->variance, readout->covariance);
+                    }
                 }
-# endif
 
-                {
-                  // Add MD5 information for readout
-                  psString headerName = NULL; ///< Header name for MD5
-                  psVector *md5 = psImageMD5(readout->image); ///< md5 hash
-                  psString md5string = psMD5toString(md5); ///< String
-
-                  pmHDU *hdu = pmHDUFromReadout(readout);
-                  psStringAppend(&headerName, "MD5_%s_%s_%d", chipName, cellName, view->readout);
-                  psMetadataAddStr(hdu->header, PS_LIST_TAIL, headerName, PS_META_REPLACE, "Image MD5", md5string);
-                  psFree(md5);
-                  psFree(md5string);
-                  psFree(headerName);
-                }
+                // Add MD5 information for readout
+                psString headerName = NULL; ///< Header name for MD5
+                psVector *md5 = psImageMD5(readout->image); ///< md5 hash
+                psString md5string = psMD5toString(md5); ///< String
+                
+                pmHDU *hdu = pmHDUFromReadout(readout);
+                psStringAppend(&headerName, "MD5_%s_%s_%d", chipName, cellName, view->readout);
+                psMetadataAddStr(hdu->header, PS_LIST_TAIL, headerName, PS_META_REPLACE, "Image MD5", md5string);
+                psFree(md5);
+                psFree(md5string);
+                psFree(headerName);
             }
 
-            psList *inputCells = psMetadataLookupPtr (NULL, cell->analysis, "INPUT.CELLS");
-            if (!pmConceptsAverageCells(cell, inputCells, NULL, NULL, false)) {
-                psError(psErrorCodeLast(), false, "Unable to average cell concepts.");
-                psFree(view);
-                return false;
+            if (fullImage) {
+                psList *inputCells = psMetadataLookupPtr (NULL, cell->analysis, "INPUT.CELLS");
+                if (!pmConceptsAverageCells(cell, inputCells, NULL, NULL, false)) {
+                    psError(psErrorCodeLast(), false, "Unable to average cell concepts.");
+                    psFree(view);
+                    return false;
+                }
             }
 
@@ -85 +89 @@
             psRegion *trimsec = psMetadataLookupPtr(NULL, cell->concepts, "CELL.TRIMSEC"); ///< Trim section
             trimsec->x0 = trimsec->x1 = trimsec->y0 = trimsec->y1 = 0; ///< All pixels
+        }
 
+        // update astrometry headers
+        pmHDU *hdu = pmHDUFromChip(chip); ///< HDU for source
+        if (!hdu || !hdu->header) {
+            psError(PM_ERR_PROG, false, "Unable to find header for output.");
+            psFree(view);
+            return false;
+        }
+
+        if (bilevelAstrometry) {
+            if (!pmAstromWriteBilevelChip(hdu->header, chip, WCS_NONLIN_TOL)) {
+                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 (!pmAstromWriteWCS(hdu->header, output, chip, WCS_NONLIN_TOL)) {
+                psError(psErrorCodeLast(), false, "Unable to read WCS astrometry for skycell.");
+                psFree(view);
+                return false;
+            }
         }
     }
@@ -93 +119 @@
         return false;
     }
-
+    
     // Update ZP from the astrometry
-    {
+    if (fullImage) {
         psMetadataItem *item = psMetadataLookup(output->concepts, "FPA.ZP");
         item->data.F32 = psMetadataLookupF32(NULL, astrom->concepts, "FPA.ZP");
     }
 
-    // need to update the chip and fpa level astrometry appropriately, see
-    // code in psastro
-    
-# if (0)
-    if (!pmAstromWriteWCS(hdu->header, outFPA, outChip, WCS_NONLIN_TOL)) {
-        psError(psErrorCodeLast(), false, "Unable to generate WCS header.");
-        return false;
+    // apply the bilevel astrometry elements to the target
+    if (bilevelAstrometry) {
+        pmHDU *phu = pmFPAviewThisPHU(view, output); ///< Astrometry PHU
+        if (!phu->header) {
+            phu->header = psMetadataAlloc ();
+        }
+        if (!pmAstromWriteBilevelMosaic(phu->header, output, WCS_NONLIN_TOL)) {
+            psError(psErrorCodeLast(), false, "Unable to write bilevel mosaic astrometry for skycell.");
+            psFree(view);
+            return false;
+        }
     }
-# endif
 
     return true;

branches/eam_branches/ipp-20130419/pswarp/src/pswarpUpdateStatistics.c

@@ -1 @@
-/** @file pswarpLoop.c
+/** @file pswarpUpdateStatistics.c
  *
- *  @brief mail processing loop for pswarp
+ *  @brief generate output statistics
  *  @ingroup pswarp
  *
  *  @author IfA
- *  @version $Revision: 1.38 $ $Name: not supported by cvs2svn $
  *  @date $Date: 2009-02-13 21:54:32 $
  *  Copyright 2009 Institute for Astronomy, University of Hawaii