Changeset 34411

Timestamp:

Sep 6, 2012, 1:24:39 PM (14 years ago)

Author:

watersc1

Message:

Broken background warping. Most of the code is in place, still trying to figure out why the results don't look good.

Location:

branches/czw_branch/20120906/pswarp/src

Files:

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

branches/czw_branch/20120906/pswarp/src/pswarp.c

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

branches/czw_branch/20120906/pswarp/src/pswarp.h

@@ -81 +81 @@
 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);

branches/czw_branch/20120906/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 = psMetadataLookupBool(&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 = psMetadataLookupBool(&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");

branches/czw_branch/20120906/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;
+}

branches/czw_branch/20120906/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
 };

branches/czw_branch/20120906/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)) {
@@ -196 +198 @@
         }
 
+        
         pmCell *cell;
         while ((cell = pmFPAviewNextCell (view, input->fpa, 1)) != NULL) {
@@ -227 +230 @@
 
                 pswarpTransformReadout(output, readout, config);
-
+                
                 if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) {
                     psError(psErrorCodeLast(), false, "Unable to write files.");
@@ -255 +258 @@
         goto DONE;
     }
-
+    
     pmCell *outCell = output->parent;   ///< Output cell
     pmChip *outChip = outCell->parent;  ///< Output chip
@@ -450 +453 @@
     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) {
+          psTrace ("pswarp", 4, "AChip %d: %x %x\n", view->chip, chip->file_exists, chip->process);
+            if (!chip->process || !chip->file_exists) { continue; }
+            if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) {
+                psError(psErrorCodeLast(), false, "Unable to read files.");
+                goto DONE;
+            }
+            pmCell *cell;
+            while ((cell = pmFPAviewNextCell (view, input->fpa, 1)) != NULL) {
+                psTrace ("pswarp", 4, "ACell %d: %x %x\n", view->cell, cell->file_exists, cell->process);
+                if (!cell->process || !cell->file_exists) { continue; }
+                if (!pmFPAfileIOChecks (config, view, PM_FPA_BEFORE) ||
+                    !pmFPAfileIOChecks (config, view, PM_FPA_AFTER)) {
+                    psError(psErrorCodeLast(), false, "Unable to read files.");
+                    goto DONE;
+                }
+            }
+            if (!pmFPAfileIOChecks (config, view, PM_FPA_AFTER)) {
+                psError(psErrorCodeLast(), false, "Unable to write files.");
+                goto DONE;
+            }
+        }
+        if (!pmFPAfileIOChecks (config, view, PM_FPA_AFTER)) {
+            psError(psErrorCodeLast(), false, "Unable to write files.");
+            goto DONE;
+        }
+        psFree(view);
+
+        pswarpFileActivation(config, detectorFiles, true);
+        pmFPAfileActivate(config->files, false, "PSWARP.ASTROM");
+    }
+
+/*     // Turn on the source output --- we need to get rid of these so that we can measure the PSF */
+/*     pmFPAfileActivate(config->files, true, "PSWARP.OUTPUT.SOURCES"); */
+
+    // 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);
+    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;
+        }
+
+        // read WCS data from the corresponding header
+        pmHDU *hdu = pmFPAviewThisHDU (view, astrom->fpa);
+        psMetadataAddS32(config->arguments,PS_LIST_TAIL,
+                         "INTERPOLATION.MODE",PS_META_REPLACE,"", PS_INTERPOLATE_BILINEAR); ///< Mode
+        psMetadataAddS32(config->arguments,PS_LIST_TAIL,
+                         "INTERPOLATION.NUM",PS_META_REPLACE,"", 64); ///< Mode
+        
+        if (bilevelAstrometry) {
+            if (!pmAstromReadBilevelChip (chip, hdu->header)) {
+                psError(psErrorCodeLast(), false, "Unable to read bilevel chip astrometry for input FPA.");
+                psFree(view);
+                psFree(stats);
+                goto DONE;
+            }
+        } else {
+            // we use a default FPA pixel scale of 1.0
+            if (!pmAstromReadWCS (input->fpa, chip, hdu->header, 1.0)) {
+                psError(psErrorCodeLast(), false, "Unable to read WCS astrometry for input FPA.");
+                psFree(view);
+                psFree(stats);
+                goto DONE;
+            }
+        }
+#define CORRECT_INPUT_WCS 1
+#if CORRECT_INPUT_WCS
+        // Correct the input WCS
+        pmAstromWCS *WCS = pmAstromWCSfromHeader(hdu->header);
+
+        double cd1f = 1.0 * 400;
+        double cd2f = 1.0 * 400;
+
+        WCS->crpix1 = WCS->crpix1 / cd1f + 0.5;
+        WCS->crpix2 = WCS->crpix2 / cd2f - 0.5;
+
+        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);
+        // End WCS work.
+#endif
+        
+        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;
+                }
+
+                // Copy the detections from the astrometry carrier to the input, so they can be accessed by
+                // pswarpTransformReadout
+                pmReadout *astromRO = pmFPAviewThisReadout(view, astrom->fpa); // Readout for astrometry
+                pmDetections *detections = psMetadataLookupPtr(&mdok, astromRO->analysis, "PSPHOT.DETECTIONS"); // Sources from astrometry
+                if (detections) {
+                    psMetadataAddPtr(readout->analysis, PS_LIST_TAIL, "PSPHOT.DETECTIONS", PS_DATA_ARRAY, "Sources from input astrometry", detections);
+                }
+
+                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, "", 2);
+/*              psMetadataAddS32(config->arguments,PS_LIST_TAIL, "GRID.NX",PS_META_REPLACE,"",100); ///< Number of grid points in x */
+/*              psMetadataAddS32(config->arguments,PS_LIST_TAIL, "GRID.NY",PS_META_REPLACE,"",100); ///< Number of grid points in y */
+
+                pswarpTransformReadout(output, readout, config);
+                
+                if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) {
+                    psError(psErrorCodeLast(), false, "Unable to write files.");
+                    goto DONE;
+                }
+            }
+            if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) {
+                psError(psErrorCodeLast(), false, "Unable to write files.");
+                goto DONE;
+            }
+        }
+        if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) {
+            psError(psErrorCodeLast(), false, "Unable to write files.");
+            goto DONE;
+        }
+    }
+
+    if (!output->data_exists) {
+        psWarning("No overlap between input and skycell.");
+/*         if (stats) { */
+/*             psMetadataAddS32(stats, PS_LIST_TAIL, "QUALITY", PS_META_REPLACE, */
+/*                              "No overlap between input and skycell", PSWARP_ERR_NO_OVERLAP); */
+/*         } */
+        psphotFilesActivate(config, false);
+        psFree(cells);
+        psFree(view);
+        goto DONE;
+    }
+    
+    pmCell *outCell = output->parent;   ///< Output cell
+    pmChip *outChip = outCell->parent;  ///< Output chip
+    pmFPA *outFPA = outChip->parent;    ///< Output FP
+
+    if (!pswarpPixelsLit(output, stats, config)) {
+        psError(psErrorCodeLast(), false, "Unable to calculate pixel regions.");
+        psFree(cells);
+        psFree(view);
+        goto DONE;
+    }
+/*     bool doStats = psMetadataLookupBool(&mdok,recipe,"MASK.STATS"); */
+/*     if (doStats) { */
+/*       if (!pswarpMaskStats(output, stats, config)) { */
+/*      psError(psErrorCodeLast(), false, "Unable to calculate mask stats."); */
+/*      psFree(cells); */
+/*      psFree(view); */
+/*      goto DONE; */
+/*       } */
+/*     } */
+/*     // Set covariance matrix for output */
+/*     { */
+/*         psList *covariances = psMetadataLookupPtr(&mdok, output->analysis, */
+/*                                                   PSWARP_ANALYSIS_COVARIANCES); // Covariance matrices */
+/*         psAssert(covariances, "Should be there"); */
+/*         psArray *covars = psListToArray(covariances); // Array of covariance matrices */
+/*         psKernel *covar = psImageCovarianceAverage(covars); */
+/*         psFree(covars); */
+/*         psMetadataRemoveKey(output->analysis, PSWARP_ANALYSIS_COVARIANCES); */
+
+/*         // Correct covariance matrix scale for the mean (square root of the) Jacobian */
+/*         double jacobian = psMetadataLookupF64(NULL, output->analysis, PSWARP_ANALYSIS_JACOBIAN); // Jacobian */
+/*         int goodPixels = psMetadataLookupS32(NULL, output->analysis, PSWARP_ANALYSIS_GOODPIX);   // Good pixels */
+/*         jacobian /= goodPixels; */
+/*         output->covariance = psImageCovarianceScale(covar, jacobian); */
+/*         psFree(covar); */
+
+/*         if (output->variance) { */
+/*             psImageCovarianceTransfer(output->variance, output->covariance); */
+/*         } */
+/*     } */
+
+/*     if (!pmConceptsAverageCells(outCell, cells, NULL, NULL, false)) { */
+/*         psError(psErrorCodeLast(), false, "Unable to average cell concepts."); */
+/*         psFree(stats); */
+/*         psFree(cells); */
+/*         psFree(view); */
+/*         goto DONE; */
+/*     } */
+/*     psFree(cells); */
+
+    psRegion *trimsec = psMetadataLookupPtr(NULL, outCell->concepts, "CELL.TRIMSEC"); ///< Trim section
+    trimsec->x0 = trimsec->x1 = trimsec->y0 = trimsec->y1 = 0; ///< All pixels
+
+    if (!psMetadataCopy(outFPA->concepts, input->fpa->concepts)) {
+        psError(psErrorCodeLast(), false, "Unable to copy FPA concepts from input to output.");
+        psFree(stats);
+        psFree(view);
+        goto DONE;
+    }
+
+    // Update ZP from the astrometry
+/*     { */
+/*         psMetadataItem *item = psMetadataLookup(outFPA->concepts, "FPA.ZP"); */
+/*         item->data.F32 = psMetadataLookupF32(NULL, astrom->fpa->concepts, "FPA.ZP"); */
+/*     } */
+
+    pmHDU *hdu = outFPA->hdu;           ///< HDU for the output warped image
+
+    // Copy header from target
+    {
+        pmFPAview *skyView = pmFPAviewAlloc(0); ///< View into skycell
+        skyView->chip = skyView->cell = 0;
+        pmCell *cell = pmFPAfileThisCell(config->files, skyView, "PSWARP.SKYCELL"); // Skycell cell
+        psFree(skyView);
+        pmHDU *skyHDU = pmHDUFromCell(cell); ///< HDU
+        if (!skyHDU) {
+            psError(PSWARP_ERR_DATA, false, "Unable to find skycell HDU.");
+            psFree(view);
+            goto DONE;
+        }
+        hdu->header = psMetadataCopy(hdu->header, skyHDU->header);
+    }
+
+    pswarpVersionHeader(hdu->header);
+    
+    if (!pmAstromWriteWCS(hdu->header, outFPA, outChip, WCS_NONLIN_TOL)) {
+        psError(psErrorCodeLast(), false, "Unable to generate WCS header.");
+        psFree(stats);
+        goto DONE;
+    }
+
+    if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) {
+        psError(psErrorCodeLast(), false, "Unable to write files.");
+        goto DONE;
+    }
+
+    // Done with the detector side of things
+    pswarpFileActivation(config, detectorFiles, false);
+    pswarpFileActivation(config, independentFiles, false);
+
+
+    // We need a new PSF model for the warped frame.  It would be good to generate this analytically, but
+    // that's going to be tricky.  We have a list of sources, so we use those to redetermine the PSF model.
+
+/*     if (psMetadataLookupBool(&mdok, recipe, "PSF")) { */
+/*         pswarpFileActivation(config, photFiles, true); */
+/*         if (!pswarpIOChecksBefore(config)) { */
+/*             psError(psErrorCodeLast(), false, "Unable to read files."); */
+/*             goto DONE; */
+/*         } */
+
+/*         // supply the readout and fpa of interest to psphot */
+/*         pmFPAfile *photFile = psMetadataLookupPtr(NULL, config->files, "PSPHOT.INPUT"); */
+/*         pmFPACopy(photFile->fpa, outFPA); */
+
+/*         pmFPAview *view = pmFPAviewAlloc(0); ///< View into skycell */
+/*         view->chip = view->cell = view->readout = 0; */
+
+/*         // grab the sources of interest from the storage location (pmFPAfile PSPHOT.INPUT.CMF) */
+/*         psArray *sources = psphotLoadPSFSources (config, view); */
+/*         if (!sources) { */
+/*             psError(psErrorCodeLast(), false, "No sources supplied to measure PSF"); */
+/*             goto DONE; */
+/*         } */
+
+/*         pmModelClassSetLimits(PM_MODEL_LIMITS_STRICT); */
+
+/*         // measure the PSF using these sources */
+/*         if (!psphotReadoutFindPSF(config, view, "PSPHOT.INPUT", sources)) { */
+/*             // This is likely a data quality issue */
+/*             // XXX Split into multiple cases using error codes? */
+/*             psErrorStackPrint(stderr, "Unable to determine PSF"); */
+/*             psWarning("Unable to determine PSF --- suspect bad data quality."); */
+/*             if (stats && psMetadataLookupS32(NULL, stats, "QUALITY") == 0) { */
+/*                 psMetadataAddS32(stats, PS_LIST_TAIL, "QUALITY", PS_META_REPLACE, */
+/*                                  "Unable to determine PSF", psErrorCodeLast()); */
+/*             } */
+/*             psErrorClear(); */
+/*             psphotFilesActivate(config, false); */
+/*         } */
+
+/*         // Ensure seeing is carried over */
+/*         pmChip *photChip = pmFPAviewThisChip(view, photFile->fpa);                 // Chip with seeing */
+/*         psMetadataItem *item = psMetadataLookup(outChip->concepts, "CHIP.SEEING"); // Concept with seeing */
+/*         item->data.F32 = psMetadataLookupF32(NULL, photChip->concepts, "CHIP.SEEING"); */
+
+/* // XXX EAM : put this in a visualization function */
+/* #if (TESTING) */
+/*         { */
+/*             #define PSF_SIZE 20         ///< Half-size of PSF */
+/*             #define PSF_FLUX 10000      ///< Central flux for PSF */
+/*             pmChip *photChip = pmFPAviewThisChip(view, photFile->fpa); */
+/*             pmPSF *psf = psMetadataLookupPtr(NULL, photChip->analysis, "PSPHOT.PSF"); */
+/*             psImage *image = psImageAlloc(2 * PSF_SIZE + 1, 2 * PSF_SIZE + 1, PS_TYPE_F32); */
+/*             psImageInit(image, 0); */
+/*             pmModel *model = pmModelFromPSFforXY(psf, PSF_SIZE, PSF_SIZE, PSF_FLUX); */
+/*             pmModelAdd(image, NULL, model, PM_MODEL_OP_FULL, 0); */
+/*             psFree(model); */
+/*             psFits *fits = psFitsOpen("psf.fits", "w"); */
+/*             psFitsWriteImage(fits, NULL, image, 0, NULL); */
+/*             psFitsClose(fits); */
+/*             psFree(image); */
+/*         } */
+/* #endif */
+
+/*         psFree(view); */
+/*     } */
+
+/*     // Perform statistics on the output image */
+/*     if (stats) { */
+/*         if (!ppStatsFPA(stats, output->parent->parent->parent, view, maskValue, config)) { */
+/*             psWarning("Unable to perform statistics on warped image."); */
+/*         } */
+/*     } */
+    
+
+    // Add MD5 information for readout
+    const char *chipName = psMetadataLookupStr(NULL, output->parent->parent->concepts, "CHIP.NAME");
+    const char *cellName = psMetadataLookupStr(NULL, output->parent->concepts, "CELL.NAME");
+    psString headerName = NULL; ///< Header name for MD5
+    psStringAppend(&headerName, "MD5_%s_%s_%d", chipName, cellName, view->readout);
+    psVector *md5 = psImageMD5(output->image); ///< md5 hash
+    psString md5string = psMD5toString(md5); ///< String
+    psFree(md5);
+    psMetadataAddStr(hdu->header, PS_LIST_TAIL, headerName, PS_META_REPLACE,
+                     "Image MD5", md5string);
+    psFree(md5string);
+    psFree(headerName);
+    psFree(view);
+
+ DONE:
+
+    return true;
+}

branches/czw_branch/20120906/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

branches/czw_branch/20120906/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
 

branches/czw_branch/20120906/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;