Changeset 35434

Timestamp:

Apr 26, 2013, 9:07:23 AM (13 years ago)

Author:

eugene

Message:

major mods to allow multichip output

Location:

branches/eam_branches/ipp-20130419/pswarp

Files:

• doc/notes.20130406.txt (modified) (2 diffs)
• src/pswarp.c (modified) (4 diffs)
• src/pswarp.h (modified) (1 diff)
• src/pswarpDefineLayout.c (modified) (6 diffs)
• src/pswarpLoadAstrometry.c (modified) (7 diffs)
• src/pswarpLoop.c (modified) (13 diffs)
• src/pswarpLoopSkycell.c (added)
• src/pswarpMatchRange.c (modified) (2 diffs)
• src/pswarpOverlaps.c (modified) (5 diffs)
• src/pswarpParseCamera.c (modified) (4 diffs)
• src/pswarpTransformReadout.c (modified) (1 diff)

branches/eam_branches/ipp-20130419/pswarp/doc/notes.20130406.txt

@@ -10 +10 @@
                         uses pmFPAfileDefineFromArgs / pmFPAfileBindFromArgs for inputs
                         uses pswarpDefineSkycell for outputs
+
+                        pmFPAfileDefineSkycell is used to generate the output
+                        FPAs.  This forces the file to be FPA-level only, which is wrong. 
+                        Not certain of the impact of replacing with pmFPAfileDefineOutput
+
+                        pmFPAfileDefineSkycell:
+                                * appends -SKYCELL to the camera name
 
   * pswarpOptions     : check the recipe options, reconcile with cmd-line options
@@ -91 +98 @@
  * need to loop over output targets as well as input sources...
 
+*** looks like I can select the correct output chips and write them.
+    Currently, the output format is inherited from the output
+    astrometry description.  that means we get a MEF output file with
+    one chip per extension.  that seems poor, but I think we can fix
+    this with the filerules.
+

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

@@ -39 +39 @@
     }
 
-# if (0)
+# 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
@@ -46 +46 @@
     }
 
+# if (0)
     pswarpDumpOutput (config);
     exit (0);
 # endif
 
+    // load and warp
+    if (!pswarpLoop(config, statsFile->md)) {
+        pswarpCleanup(config, statsFile);
+    }
+    pswarpDumpOutput (config);
+    exit (0);
+
+# else
     // load the skycell layout information:
     //   read the output astrometry description
@@ -60 +69 @@
         pswarpCleanup(config, statsFile);
     }
-
     // load and warp
-    if (!pswarpLoop(config, statsFile->md)) {
+    if (!pswarpLoopSkycell(config, statsFile->md)) {
         pswarpCleanup(config, statsFile);
     }
@@ -69 +77 @@
         pswarpCleanup(config, statsFile);
     }
-
+# endif
     psLogMsg("pswarp", PS_LOG_INFO, "complete pswarp run: %f sec\n", psTimerMark("pswarp"));
     pswarpCleanup(config, statsFile);

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

@@ -194 +194 @@
 pswarpBounds *pswarpBoundsAlloc();
 
-bool pswarpLoadAstrometry (pmConfig *config, pmFPAfile *astrom, pmFPAfile *target);
+bool pswarpLoadAstrometry (pmFPAfile *target, pmFPAfile *astrom, pmConfig *config);
+
+bool pswarpTransformToTarget (pmFPA *output, pmReadout *input, pmConfig *config);
+bool pswarpMakePSF (pmFPA *output);
+bool pswarpUpdateStatistics (pmFPA *output, pmFPA *input, pmConfig *config);

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

@@ -25 +25 @@
         astrom = input;
     }
-    if (!pswarpLoadAstrometry (config, input, astrom)) {
+    if (!pswarpLoadAstrometry (input, astrom, config)) {
         psError(PSWARP_ERR_CONFIG, false, "problem loading input astrometry\n");
         return false;
@@ -41 +41 @@
         return false;
     }
-    if (!pswarpLoadAstrometry (config, output, skycell)) {
+    if (!pswarpLoadAstrometry (output, skycell, config)) {
         psError(PSWARP_ERR_CONFIG, false, "problem loading output astrometry\n");
         return false;
@@ -50 +50 @@
 
     // find the R,D center for the skycell, use for common projection
-    psProjection *frame = pswarpLocalFrame (output->fpa);
+    psProjection *frame = pswarpLocalFrame (skycell->fpa);
 
     // generate Lmin,max, Mmin,max for both datasets
+    fprintf (stderr, "srcBounds:\n");
+    // pswarpBounds *srcBounds = pswarpMakeBounds (astrom->fpa, frame);
     pswarpBounds *srcBounds = pswarpMakeBounds (input->fpa, frame);
-    pswarpBounds *tgtBounds = pswarpMakeBounds (output->fpa, frame);
+
+    fprintf (stderr, "srcBounds:\n");
+    pswarpBounds *tgtBounds = pswarpMakeBounds (skycell->fpa, frame);
+
     pswarpFindOverlap (output->fpa, srcBounds, tgtBounds);
 
+    // XXX the loop below generates the output pixels. should this be deferred until 
+    // we actually need them?
     // Generate the output chips (pixels on output->fpa, concepts from skycell->fpa)
     pmFPAview *view = pmFPAviewAlloc(0);
@@ -63 +70 @@
     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; }
+        if (!chip->process) { 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; }
+            if (!cell->process) { continue; }
+
+            pmCell *refcell = pmFPAviewThisCell(view, skycell->fpa); ///< Target cell
 
             // we've got the output astrom header
-            pmHDU *hdu = pmHDUFromCell(cell); ///< HDU for source
+            pmHDU *hdu = pmHDUFromCell(refcell); ///< HDU for source
             if (!hdu || !hdu->header) {
                 psError(PM_ERR_PROG, false, "Unable to find header for sky cell.");
@@ -92 +101 @@
             // copy the image concepts from the skycell 
             bool status = false;
-            pmCell *refcell = pmFPAviewThisCell(view, skycell->fpa); ///< Target cell
             psMetadataItemSupplement(&status, cell->concepts, refcell->concepts, "CELL.XBIN");
             psMetadataItemSupplement(&status, cell->concepts, refcell->concepts, "CELL.YBIN");
@@ -101 +109 @@
             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;
+
+# if (0)
+            // ******* this chunk comes from pswarpLoop, but should be done when the output chip is generated
+            pmHDU *outHDU = pmHDUFromCell (cell);           ///< HDU for the output warped image
+            outHDU->header = psMetadataCopy(outHDU->header, hdu->header);
+            pswarpVersionHeader(outHDU->header);
+# endif
         }
     }

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

@@ -14 +14 @@
 // WCS and read the headers.  We place the resulting info on the target astrometry
 // containers
-bool pswarpLoadAstrometry (pmConfig *config, pmFPAfile *astrom, pmFPAfile *target) {
+bool pswarpLoadAstrometry (pmFPAfile *target, pmFPAfile *astrom, pmConfig *config) {
 
     pmChip *chip = NULL;
     pmCell *cell = NULL;
+
+    // XXX assert that astrom be WCS or CMF?
 
     // XXX set the type to be WCS
@@ -42 +44 @@
         while ((cell = pmFPAviewNextCell (view, astrom->fpa, 1)) != NULL) {
             psTrace ("pswarp", 4, "Cell %d: %x %x\n", view->cell, cell->file_exists, cell->process);
-            if (!cell->process || !cell->file_exists) { continue; }
+            if (!cell->process) { continue; }
             if (!pmFPAfileRead (astrom, view, config)) {
                 psError(psErrorCodeLast(), false, "failed READ at CELL %s", astrom->name);
@@ -89 +91 @@
     // apply the bilevel astrometry elements to the target
     if (bilevelAstrometry) {
-        if (!pmAstromReadBilevelMosaic(target->fpa, phu->header)) {
+        if (!pmAstromReadBilevelMosaic(astrom->fpa, phu->header)) {
             psError(psErrorCodeLast(), false, "Unable to read bilevel mosaic astrometry for skycell.");
             psFree(view);
@@ -95 +97 @@
         }
     }
+    // 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);
+    target->fpa->fromTPA = psMemIncrRefCounter (astrom->fpa->fromTPA);
+    target->fpa->toSky = psMemIncrRefCounter (astrom->fpa->toSky);
 
     while ((chip = pmFPAviewNextChip (view, astrom->fpa, 1)) != NULL) {
@@ -110 +116 @@
         // We read from the astrometry source into the target.
         pmChip *targetChip = pmFPAviewThisChip(view, target->fpa); ///< Chip in the output
+        if (!targetChip) continue; // only load astrometry into output chips which exist!
+
         if (bilevelAstrometry) {
-            if (!pmAstromReadBilevelChip(targetChip, hdu->header)) {
+            if (!pmAstromReadBilevelChip(chip, hdu->header)) {
                 psError(psErrorCodeLast(), false, "Unable to read bilevel chip astrometry for skycell.");
                 psFree(view);
@@ -118 +126 @@
         } else {
             // we use a default FPA pixel scale of 1.0
-            if (!pmAstromReadWCS(target->fpa, targetChip, hdu->header, 1.0)) {
+            if (!pmAstromReadWCS(target->fpa, chip, hdu->header, 1.0)) {
                 psError(psErrorCodeLast(), false, "Unable to read WCS astrometry for skycell.");
                 psFree(view);
@@ -124 +132 @@
             }
         }
+        targetChip->toFPA = psMemIncrRefCounter (chip->toFPA);
+        targetChip->fromFPA = psMemIncrRefCounter (chip->fromFPA);
     }
 

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

@@ -20 +20 @@
 bool pswarpLoop(pmConfig *config, psMetadata *stats)
 {
-    bool status;
-    bool mdok;                          // Status of MD lookup
-
     const char *skyCamera = psMetadataLookupStr(NULL, config->arguments,
                                                 "SKYCELL.CAMERA");  // Name of camera for skycell
@@ -29 +26 @@
 
     // load the recipe
+    bool status = false;
     psMetadata *recipe = psMetadataLookupPtr (&status, config->recipes, PSWARP_RECIPE);
     if (!recipe) {
@@ -35 +33 @@
     }
 
+    // XXX this should probably happen elsewhere (logically it belongs to file I/O, right?)
     if (!pswarpSetMaskBits(config)) {
         psError(psErrorCodeLast(), false, "failed to set mask bits");
         return NULL;
     }
-
-    // output mask bits
-    psImageMaskType maskValue = psMetadataLookupImageMask(&status, recipe, "MASK.OUTPUT");
-    psAssert (status, "MASK.OUTPUT was not defined");
 
     // select the input data sources
@@ -62 +57 @@
     }
 
-    // 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");
+    // select the input data sources
+    pmFPAfile *output = psMetadataLookupPtr(NULL, config->files, "PSWARP.OUTPUT");
     if (!output) {
         psError(PSWARP_ERR_CONFIG, true, "Can't find output 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.
-    // XXX why is this being done here, doesn't it duplicate the work in pswarpDefine.c??
+
+    // use the external astrometry source if supplied
+    pmFPAfile *skycell = psMetadataLookupPtr(NULL, config->files, "PSWARP.SKYCELL");
+    if (!skycell) {
+        psError(PSWARP_ERR_DATA, true, "Cannot find output astrometry.");
+        return false;
+    }
+
+    // XXX NOTE : 'output' was a readout, now we need to loop over all output targets
+
+    // ensure everyone is off except what we need below
     pswarpFileActivation(config, detectorFiles, false);
     pswarpFileActivation(config, photFiles, false);
@@ -83 +79 @@
     pswarpFileActivation(config, skycellFiles, false);
     
-    // pswarpFileActivation(config, skycellFiles, true);
-    // if (0) {
-    //   if (!pswarpIOChecksBefore(config)) {
-    //     psError(psErrorCodeLast(), false, "Unable to read files.");
-    //     goto DONE;
-    //   }
-    // }
-   
-    // Read the input astrometry
-    // XXX rather than use the activations here, this should just explicitly loop over the desired filerule
-    {
-        pmFPAfileActivate(config->files, true, "PSWARP.ASTROM");
-
-        pmChip *chip;
-        pmFPAview *view = pmFPAviewAlloc(0);
-        if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) {
-            psError(psErrorCodeLast(), false, "Unable to read files.");
-            goto DONE;
-        }
-        while ((chip = pmFPAviewNextChip (view, input->fpa, 1)) != NULL) {
-            psTrace ("pswarp", 4, "Chip %d: %x %x\n", view->chip, chip->file_exists, chip->process);
-            if (!chip->process || !chip->file_exists) { continue; }
-            if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) {
-                psError(psErrorCodeLast(), false, "Unable to read files.");
-                goto DONE;
-            }
-            pmCell *cell;
-            while ((cell = pmFPAviewNextCell (view, input->fpa, 1)) != NULL) {
-                psTrace ("pswarp", 4, "Cell %d: %x %x\n", view->cell, cell->file_exists, cell->process);
-                if (!cell->process || !cell->file_exists) { continue; }
-                if (!pmFPAfileIOChecks (config, view, PM_FPA_BEFORE) ||
-                    !pmFPAfileIOChecks (config, view, PM_FPA_AFTER)) {
-                    psError(psErrorCodeLast(), false, "Unable to read files.");
-                    goto DONE;
-                }
-            }
-            if (!pmFPAfileIOChecks (config, view, PM_FPA_AFTER)) {
-                psError(psErrorCodeLast(), false, "Unable to write files.");
-                goto DONE;
-            }
-        }
-        if (!pmFPAfileIOChecks (config, view, PM_FPA_AFTER)) {
-            psError(psErrorCodeLast(), false, "Unable to write files.");
-            goto DONE;
-        }
-        psFree(view);
-
-        pswarpFileActivation(config, detectorFiles, true);
-        pmFPAfileActivate(config->files, false, "PSWARP.ASTROM");
-    }
-
     // 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);
-            goto DONE;
-        }
-    }
-
-    psList *cells = psListAlloc(NULL);  // List of cells, for concepts averaging
+    // XXX NOTE : do not tackle source issues yet...
+    // XXX pmFPAfileActivate(config->files, true, "PSWARP.OUTPUT.SOURCES");
+
+    // Done with the detector side of things
+    pswarpFileActivation(config, detectorFiles, true);
+    pswarpFileActivation(config, independentFiles, true);
+
+    pmFPAview *view = pmFPAviewAlloc(0);
+
+    // XXX : need to generate a list of input cells for each output cell
+    // XXX : psList *cells = psListAlloc(NULL);  // List of cells, for concepts averaging
 
     // files associated with the science image
@@ -166 +98 @@
     }
 
-    // *** main transformation block
+    // *** main transformation block 
+    // *** this section loops over the input chips/cells and reads them one at a time
+    // *** the output chips/cells are filled where appropriate, but not yet written to disk 
     pmChip *chip;
     while ((chip = pmFPAviewNextChip (view, input->fpa, 1)) != NULL) {
@@ -176 +110 @@
         }
 
-        // read WCS data from the corresponding header
-        pmHDU *hdu = pmFPAviewThisHDU (view, astrom->fpa);
-
-        
-        if (bilevelAstrometry) {
-            if (!pmAstromReadBilevelChip (chip, hdu->header)) {
-                psError(psErrorCodeLast(), false, "Unable to read bilevel chip astrometry for input FPA.");
-                psFree(view);
-                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);
-                goto DONE;
-            }
-        }
-        
         pmCell *cell;
         while ((cell = pmFPAviewNextCell (view, input->fpa, 1)) != NULL) {
@@ -203 +118 @@
                 goto DONE;
             }
-
-            psListAdd(cells, PS_LIST_TAIL, cell);
 
             // process each of the readouts
@@ -217 +130 @@
                 }
 
+/// XXX ignore sources for now
+# if (0)
                 // Copy the detections from the astrometry carrier to the input, so they can be accessed by
                 // pswarpTransformReadout
@@ -224 +139 @@
                     psMetadataAddPtr(readout->analysis, PS_LIST_TAIL, "PSPHOT.DETECTIONS", PS_DATA_ARRAY, "Sources from input astrometry", detections);
                 }
-
-                pswarpTransformReadout(output, readout, config);
+# endif
+
+                pswarpTransformToTarget(output->fpa, readout, config);
                 
                 if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) {
@@ -242 +158 @@
         }
     }
-
-    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(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(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.");
-        goto DONE;
-    }
-
     if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) {
         psError(psErrorCodeLast(), false, "Unable to write files.");
@@ -352 +167 @@
     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.
+    // XXXX turn these on later .. pswarpUpdateStatistics (output);
+    // XXXX turn these on later .. pswarpMakePSF (output);
+
+    // XXX have not yet written the output to disk
+
+    psFree(view);
+
+ DONE:
+    return true;
+}
+
+# if (0)
+// once the output fpa elements have been built, loop over the fpa and generate stats
+// for each readout
+bool pswarpUpdateStatistics (pmFPA *output, pmFPA *input, pmConfig *config)  {
+
+    // output mask bits
+    psImageMaskType maskValue = psMetadataLookupImageMask(&status, recipe, "MASK.OUTPUT");
+    psAssert (status, "MASK.OUTPUT was not defined");
+
+    pmFPAview *view = pmFPAviewAlloc(0);
+    
+    pmChip *chip;
+    while ((chip = pmFPAviewNextChip (view, output, 1)) != NULL) {
+        psTrace ("pswarp", 4, "Chip %d: %x %x\n", view->chip, chip->file_exists, chip->process);
+        if (!chip->process || !chip->file_exists) { continue; }
+
+        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) {
+                    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;
+                }
+    
+                if (!pswarpPixelsLit(readout, stats, config)) {
+                    psError(psErrorCodeLast(), false, "Unable to calculate pixel regions.");
+                    psFree(view);
+                    goto DONE;
+                }
+                bool doStats = psMetadataLookupBool(&mdok,recipe,"MASK.STATS");
+                if (doStats) {
+                    if (!pswarpMaskStats(readout, 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);
+                    }
+                }
+            }
+
+            // XXX need to get the input Cells for this output cell (save on the analysis somewhere
+            // XXX skip this for now
+# if (0)
+            psArray *inputCells = psMetadataLookupPtr (NULL, cell->analysis, "something");
+            if (!pmConceptsAverageCells(cell, inputCells, NULL, NULL, false)) {
+                psError(psErrorCodeLast(), false, "Unable to average cell concepts.");
+                psFree(view);
+                goto DONE;
+            }
+# endif
+            // XXX Is this too ad-hoc?
+            psRegion *trimsec = psMetadataLookupPtr(NULL, cell->concepts, "CELL.TRIMSEC"); ///< Trim section
+            trimsec->x0 = trimsec->x1 = trimsec->y0 = trimsec->y1 = 0; ///< All pixels
+
+        }
+
+        // Perform statistics on the output image
+        if (stats) {
+            // XXX this should be done at each level...
+            if (!ppStatsFPA(stats, output->parent->parent->parent, view, maskValue, config)) {
+                psWarning("Unable to perform statistics on warped image.");
+            }
+        }
+    }
+    if (!psMetadataCopy(output->concepts, input->concepts)) {
+        psError(psErrorCodeLast(), false, "Unable to copy FPA concepts from input to output.");
+        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");
+    }
+
+    // need to update the chip and fpa level astrometry appropriate, see
+    // code in psastro
+    
+# if (0)
+    // 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);
+# endif
+
+# if (0)
+    if (!pmAstromWriteWCS(hdu->header, outFPA, outChip, WCS_NONLIN_TOL)) {
+        psError(psErrorCodeLast(), false, "Unable to generate WCS header.");
+        goto DONE;
+    }
+# endif
+    return true;
+}
+# endif
+
+# if (0)
+// 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.
+bool pswarpMakePSF (pmConfig *config) {
+
+    // 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 (psMetadataLookupBool(&mdok, recipe, "PSF")) {
@@ -419 +391 @@
         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;
 }
+# endif
+
+// once the output fpa elements have been built, loop over the fpa and generate stats
+// for each readout
+bool pswarpTransformToTarget (pmFPA *output, pmReadout *input, pmConfig *config)  {
+
+    pmFPAview *view = pmFPAviewAlloc(0);
+    
+    pmChip *chip;
+    while ((chip = pmFPAviewNextChip (view, output, 1)) != NULL) {
+        psTrace ("pswarp", 4, "Chip %d: %x %x\n", view->chip, chip->file_exists, chip->process);
+        if (!chip->process || !chip->file_exists) { continue; }
+
+        pmCell *cell;
+        while ((cell = pmFPAviewNextCell (view, output, 1)) != NULL) {
+            psTrace ("pswarp", 4, "Cell %d: %x %x\n", view->cell, cell->file_exists, cell->process);
+            if (!cell->process || !cell->file_exists) { continue; }
+
+            // process each of the readouts
+            pmReadout *readout;
+            while ((readout = pmFPAviewNextReadout(view, output, 1)) != NULL) {
+                pswarpTransformReadout (readout, input, config);
+            }
+        }
+    }
+    return true;
+}
+

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

@@ -98 +98 @@
 
     // demo forward and backward transformation
-# if (0)
+# if (1)
     srcPix->x = *minX;
     srcPix->y = *minY;
@@ -115 +115 @@
     psPlaneTransformApply (srcFP, fpaSrc->fromTPA, srcTP); 
     psPlaneTransformApply (srcPix, chipSrc->fromFPA, srcFP); 
-    fprintf (stderr, "-> %f,%f ", srcPix->x, srcPix->y);
+    fprintf (stderr, "-> %f,%f\n", srcPix->x, srcPix->y);
 # endif
 

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);
     }
@@ -109 +110 @@
     pmFPAview *view = pmFPAviewAlloc (0);
 
+    // we are looping over the chips of the OUTPUT fpa.  chip->file_exists is false for all chips. 
+    // at this point, I cannot exclude any of these chips (I do not know which output chips will not be created).
+    // NOTE : I either have to use SKYCELL fpa and ensure SKYCELL->chip->toFPA exists or grab hdu from SKYCELL
+
     pmChip *chip = NULL;
     while ((chip = pmFPAviewNextChip (view, fpa, 1)) != NULL) {
         psTrace ("pswarp", 4, "Chip %d: %x %x\n", view->chip, chip->file_exists, chip->process);
-        if (!chip->process || !chip->file_exists) { continue; }
+        // if (!chip->file_exists) { continue; } // this fpa refers to the output file, so no chip exist
+        assert (chip->toFPA);
+        chip->process = true; // start with process true, we will de-activate later (pswarpFindOverlaps)
 
         // we've got the output astrom header
@@ -153 +160 @@
     double Dmid = stats->sampleMean;
 
-    psProjection *frame = psProjectionAlloc (Rmid, Dmid, 1.0, 1.0, PS_PROJ_TAN);
+    psProjection *frame = psProjectionAlloc (Rmid, Dmid, PS_RAD_DEG/3600.0, PS_RAD_DEG/3600.0, PS_PROJ_TAN);
     return frame;
 }
@@ -162 +169 @@
     // elements overlap it.
 
-    pmFPAview *view = pmFPAviewAlloc (0);
+    assert (fpa->chips->n == tgt->Pmin->n);
 
-    pmChip *chip = NULL;
-    int Nout = 0;
-    while ((chip = pmFPAviewNextChip (view, 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; }
+    for (int j = 0; j < fpa->chips->n; j++) {
+        pmChip *chip = fpa->chips->data[j];
 
         // we have src bounds
-        float Pmin = tgt->Pmin->data.F32[Nout];
-        float Pmax = tgt->Pmax->data.F32[Nout];
-        float Qmin = tgt->Qmin->data.F32[Nout];
-        float Qmax = tgt->Qmax->data.F32[Nout];
+        float Pmin = tgt->Pmin->data.F32[j];
+        float Pmax = tgt->Pmax->data.F32[j];
+        float Qmin = tgt->Qmin->data.F32[j];
+        float Qmax = tgt->Qmax->data.F32[j];
 
         bool valid = false;
@@ -183 +187 @@
             valid = Pvalid && Qvalid;
         }
+        // if (valid) {
+        //     fprintf (stderr, "chip %d is valid\n", j);
+        // } else {
+        //     fprintf (stderr, "chip %d is NOT valid\n", j);
+        // }
         chip->process = valid;
     }

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

@@ -128 +128 @@
 
     // The output skycell
-    pmFPAfile *output = pmFPAfileDefineSkycell(config, NULL, "PSWARP.OUTPUT");
+    // XXX The output needs to be in a more generic format, not just skycell.
+    // I need to understand how badly things depend on SKYCELL being a special format
+    // (is it just pswarp, or is it everything downstream as well??)
+    // XXX for a test, generate a generic output (based on output target astrometry)
+    // XXX pmFPAfile *output = pmFPAfileDefineSkycell(config, NULL, "PSWARP.OUTPUT");
+    pmFPAfile *output = pmFPAfileDefineOutputForFormat(config, NULL, "PSWARP.OUTPUT", skycell->cameraName, skycell->formatName);
     if (!output) {
         psError(psErrorCodeLast(), false, "Failed to build FPA from PSWARP.OUTPUT");
@@ -135 +140 @@
     output->save = true;
 
+    // XXX for a test, skip these
+# if (0)
     pmFPAfile *outMask = pmFPAfileDefineSkycell(config, output->fpa, "PSWARP.OUTPUT.MASK");
     if (!outMask) {
@@ -159 +166 @@
         outSources->save = true;
     }
+# endif
 
     if (inBackground && psMetadataLookupBool(&mdok, recipe, "BACKGROUND.MODEL")) {
@@ -172 +180 @@
     }
     
-    if (psMetadataLookupBool(&mdok, recipe, "PSF")) {
+    /// XXX TEST 
+    if (false && psMetadataLookupBool(&mdok, recipe, "PSF")) {
         // This file, PSPHOT.INPUT, is just used as a carrier; output files (eg, PSPHOT.RESID) are defined by
         // psphotDefineFiles

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

@@ -117 +117 @@
     psAssert (xGridMax < grid->nXpts, "xGridMax too big\n");
     psAssert (yGridMax < grid->nYpts, "yGridMax too big\n");
+
+    fprintf (stderr, "warp %d,%d - %d,%d\n", xGridMin, yGridMin, xGridMax, yGridMax);
 
     // create jobs and supply them to the threads