Changeset 34800 for trunk/pswarp/src

Timestamp:

Dec 11, 2012, 2:04:31 PM (14 years ago)

Author:

watersc1

Message:

Merge from my branch including background restoration and stack stage projection cell binned images.

Location:

trunk

Files:

• . (modified) (1 prop)
• pswarp/src/pswarp.c (modified) (2 diffs)
• pswarp/src/pswarp.h (modified) (2 diffs)
• pswarp/src/pswarpArguments.c (modified) (4 diffs)
• pswarp/src/pswarpDefine.c (modified) (3 diffs)
• pswarp/src/pswarpFileNames.h (modified) (2 diffs)
• pswarp/src/pswarpLoop.c (modified) (5 diffs)
• pswarp/src/pswarpParseCamera.c (modified) (2 diffs)
• pswarp/src/pswarpTransformReadout.c (modified) (6 diffs)
• pswarp/src/pswarpTransformSources.c (modified) (1 diff)
• pswarp/src/pswarpTransformTile.c (modified) (2 diffs)

trunk/pswarp/src/pswarp.c

@@ -46 +46 @@
         goto DIE;
     }
+    if (psMetadataLookupBool(NULL,config->arguments,"BACKGROUND.MODEL")) {
+      if (!pswarpDefineBackground(config)) {
+        goto DIE;
+      }
+    }
 
     // Open the statistics file
@@ -66 +71 @@
     if (!pswarpLoop(config, stats)) {
         goto DIE;
+    }
+    if (psMetadataLookupBool(NULL,config->arguments,"BACKGROUND.MODEL")) {
+      if (!pswarpLoopBackground(config, stats)) {
+        fprintf(stderr,"Dying!\n");
+        goto DIE;
+      }
     }
 

trunk/pswarp/src/pswarp.h

@@ -62 +62 @@
     psImage *region;
 
+    /** values which are needed to control the background model warping. */
+    bool background_warping;
+    double offset_x;
+    double offset_y;
+  
     /** input values for this tile */
     int gridX;
@@ -81 +86 @@
 bool pswarpParseCamera (pmConfig *config);
 bool pswarpDefine (pmConfig *config);
+bool pswarpDefineBackground (pmConfig *config);
 bool pswarpLoop (pmConfig *config, psMetadata *stats);
+bool pswarpLoopBackground (pmConfig *config, psMetadata *stats);
 psExit pswarpExitCode(psExit exitValue);
 bool pswarpTransformReadout (pmReadout *output, pmReadout *input, pmConfig *config);

trunk/pswarp/src/pswarpArguments.c

@@ -86 +86 @@
     pmConfigFileSetsMD (config->arguments, &argc, argv, "MASK", "-mask", "-masklist");
     pmConfigFileSetsMD (config->arguments, &argc, argv, "VARIANCE", "-variance", "-variancelist");
-
+    pmConfigFileSetsMD (config->arguments, &argc, argv, "BACKGROUND", "-background", "-bkglist");
+    
     if (argc != 3) {
         usage();
@@ -164 +165 @@
     }
 
+    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,
@@ -176 +194 @@
                      "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,
@@ -189 +211 @@
                      "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");

trunk/pswarp/src/pswarpDefine.c

@@ -41 +41 @@
         return false;
     }
-
+    
     // open the full skycell file; no need to defer different depths. only load the header data
     pmFPAview *view = pmFPAviewAlloc (0);
@@ -86 +86 @@
         psMetadataItemSupplement(&status, target->concepts, source->concepts, "CELL.X0");
         psMetadataItemSupplement(&status, target->concepts, source->concepts, "CELL.Y0");
+
+        
     }
 
@@ -123 +125 @@
         }
     }
-
+    
     view->chip = view->cell = view->readout = -1;
     pmFPAAddSourceFromView(output->fpa, view, output->format);
+
 
     psFree (view);
     return true;
 }
+
+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/pswarpFileNames.h

@@ -12 +12 @@
   "PSWARP.MASK",
   "PSWARP.VARIANCE",
+  "PSWARP.BKGMODEL",
   NULL
 };
@@ -20 +21 @@
   "PSWARP.OUTPUT.MASK",
   "PSWARP.OUTPUT.VARIANCE",
+  "PSWARP.OUTPUT.BKGMODEL",
   NULL
 };

trunk/pswarp/src/pswarpLoop.c

@@ -179 +179 @@
         // read WCS data from the corresponding header
         pmHDU *hdu = pmFPAviewThisHDU (view, astrom->fpa);
+
+        
         if (bilevelAstrometry) {
             if (!pmAstromReadBilevelChip (chip, hdu->header)) {
@@ -195 +197 @@
             }
         }
-
+        
         pmCell *cell;
         while ((cell = pmFPAviewNextCell (view, input->fpa, 1)) != NULL) {
@@ -227 +229 @@
 
                 pswarpTransformReadout(output, readout, config);
-
+                
                 if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) {
                     psError(psErrorCodeLast(), false, "Unable to write files.");
@@ -255 +257 @@
         goto DONE;
     }
-
+    
     pmCell *outCell = output->parent;   ///< Output cell
     pmChip *outChip = outCell->parent;  ///< Output chip
@@ -450 +452 @@
     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/pswarpParseCamera.c

@@ -84 +84 @@
     }
 
+    pmFPAfile *inBackground = defineInputFile(config, NULL, "PSWARP.BKGMODEL", "BACKGROUND",
+                                              PM_FPA_FILE_IMAGE);
+    if (!inBackground) {
+      psLogMsg("pswarp", PS_LOG_INFO, "No background models supplied");
+    }
+    
     // The input skycell is a required argument: it defines the output image
     // XXX we may need a different skycell structure here
@@ -134 +140 @@
     }
 
+    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

trunk/pswarp/src/pswarpTransformReadout.c

@@ -29 +29 @@
     psImageInterpolateMode interpolationMode = psMetadataLookupS32(NULL, config->arguments,
                                                                    "INTERPOLATION.MODE"); ///< Mode for interp
+
     int numKernels = psMetadataLookupS32(NULL, config->arguments, "INTERPOLATION.NUM"); ///< Number of kernels
 
@@ -60 +61 @@
     // 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
     double maxError = pswarpMapGridMaxError(grid); // Maximum (positional) error from using grid
@@ -130 +132 @@
             args->goodPixels = 0;
 
+            if (psMetadataLookupBool(NULL,config->arguments,"BACKGROUND_WARPING")) {
+              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
             psThreadJob *job = psThreadJobAlloc ("PSWARP_TRANSFORM_TILE");
@@ -171 +180 @@
         } else {
             pswarpTransformTileArgs *args = job->args->data[0];
-            // fprintf (stderr, "finished job %d,%d, Nargs: %ld\n", args->gridX, args->gridY, job->args->n);
             goodPixels += args->goodPixels;
             xMin = PS_MIN(args->xMin, xMin);
@@ -183 +191 @@
                 jacobian += args->jacobian * args->goodPixels;
             }
+
+            
         }
         psFree(job);
@@ -201 +211 @@
 
     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;
+          psError(psErrorCodeLast(), false, "Unable to interpolate image.");
+          return false;
         }
+      }
     }
 

trunk/pswarp/src/pswarpTransformSources.c

@@ -41 +41 @@
     if (!outDetections) {
         outDetections = pmDetectionsAlloc();
-        psMetadataAddPtr(output->analysis, PS_LIST_TAIL, "PSPHOT.DETECTIONS", PS_DATA_ARRAY, "Warped sources", outDetections);
+        psMetadataAddPtr(output->analysis, PS_LIST_TAIL, "PSPHOT.DETECTIONS", PS_DATA_ARRAY , "Warped sources", outDetections);
     }
     psArray *outSources = outDetections->allSources;

trunk/pswarp/src/pswarpTransformTile.c

@@ -46 +46 @@
     args->jacobian = NAN;
 
+    args->background_warping = false;
+    args->offset_x = 0.0;
+    args->offset_y = 0.0;
+    
     return args;
 }
@@ -83 +87 @@
     // Iterate over the output image pixels (parent frame)
     long goodPixels = 0;                ///< Number of input pixels landing on the output image
+
     for (int y = yMin; y < yMax; y++) {
         for (int x = xMin; x < xMax; x++) {
-
             // Only transform those pixels requested
             if (region && region->data.PS_TYPE_IMAGE_MASK_DATA[y][x]) {
                 continue;
             }
-
             // pswarpMapApply converts the output coordinate (x,y) to the input coordinate.
             // both are in the parent frames of the input and output images.
             double xIn, yIn;            // Input pixel coordinates
             pswarpMapApply(&xIn, &yIn, map, x + 0.5, y + 0.5);
+
+            // This needs to use a more reliable method to do this offset and limiting
+            //      if (args->interp->mode == 8) {
+            if (args->background_warping) {
+              //              double xOffset = 177.0 / 400.0; // (modelsize * modelbinning - xsize) / 2.0
+              //              double yOffset = 166.0 / 400.0; // (modelsize * modelbinning - ysize) / 2.0
+              xIn += args->offset_x;
+              yIn += args->offset_y;
+
+              if ((xIn > inNumCols - args->offset_x)||
+                  (yIn > inNumRows - args->offset_y)||
+                  (xIn < args->offset_x)||
+                  (yIn < args->offset_y)) {
+                continue;
+              }
+            }
+            
             if (xIn < 0 || xIn >= inNumCols || yIn < 0 || yIn >= inNumRows) {
                 continue;