Changeset 35512

Timestamp:

May 5, 2013, 2:20:46 PM (13 years ago)

Author:

eugene

Message:

stats is basically working now

Location:

branches/eam_branches/ipp-20130419/pswarp

Files:

• doc/notes.20130406.txt (modified) (1 diff)
• src/Makefile.am (modified) (1 diff)
• src/pswarp.c (modified) (3 diffs)
• src/pswarp.h (modified) (2 diffs)
• src/pswarpCleanup.c (modified) (2 diffs)
• src/pswarpDefineLayout.c (modified) (2 diffs)
• src/pswarpLoop.c (modified) (9 diffs)
• src/pswarpMaskStats.c (modified) (2 diffs)
• src/pswarpOverlaps.c (modified) (4 diffs)
• src/pswarpPixelsLit.c (modified) (1 diff)
• src/pswarpStatsFile.c (modified) (1 diff)

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

@@ -104 +104 @@
     this with the filerules.
 
+** trouble running warp -> chip using astrometry from skycell (ok using astrom from warp)

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

@@ -50 +50 @@
         pswarpTransformSources.c        \
         pswarpTransformTile.c           \
+        pswarpUpdateStatistics.c \
         pswarpVersion.c                 \
         pswarpFiles.c

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

@@ -45 +45 @@
         pswarpCleanup(config, statsFile);
     }
-
-# if (0)
-    pswarpDumpOutput (config);
-    exit (0);
-# endif
-
     // load and warp
     if (!pswarpLoop(config, statsFile->md)) {
         pswarpCleanup(config, statsFile);
     }
-    pswarpDumpOutput (config);
-    exit (0);
 
 # else
@@ -78 +70 @@
     }
 # endif
+
     psLogMsg("pswarp", PS_LOG_INFO, "complete pswarp run: %f sec\n", psTimerMark("pswarp"));
     pswarpCleanup(config, statsFile);
@@ -83 +76 @@
 
 // NOTE: pswarpCleanup calls exit
+
+// pswarpDumpOutput (config);
+// exit (0);

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

@@ -189 +189 @@
 } pswarpBounds;
 
-bool pswarpFindOverlap (pmFPA *fpa, pswarpBounds *src, pswarpBounds *tgt);
+bool pswarpFindOverlap (pmFPA *input, pmFPA *output, pswarpBounds *src, pswarpBounds *tgt);
 psProjection *pswarpLocalFrame (pmFPA *fpa);
 pswarpBounds *pswarpMakeBounds (pmFPA *fpa, psProjection *frame);
@@ -199 +199 @@
 bool pswarpTransformToTarget (pmFPA *output, pmReadout *input, pmConfig *config);
 bool pswarpMakePSF (pmFPA *output);
-bool pswarpUpdateStatistics (pmFPA *output, pmFPA *input, pmConfig *config);
+bool pswarpUpdateStatistics (pmFPA *output, psMetadata *stats, pmFPA *input, pmFPA *astrom, pmConfig *config);

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

@@ -18 +18 @@
 
     // Ensure everything is written out, at every level
-    pswarpFileActivation(config, detectorFiles, true);
+    pswarpFileActivation(config, detectorFiles, false);
+    pswarpFileActivation(config, photFiles, false);
+    pswarpFileActivation(config, independentFiles, false);
     pswarpFileActivation(config, skycellFiles, true);
-    pswarpFileActivation(config, photFiles, true);
-    pswarpFileActivation(config, independentFiles, true);
 
-    if (!pswarpIOChecksAfter(config)) {
+    pmFPAfile *output = psMetadataLookupPtr(NULL, config->files, "PSWARP.OUTPUT");
+    if (!output) {
+        psError(PSWARP_ERR_CONFIG, false, "Can't find output data!\n");
+        pmFPAfileFreeSetStrict(false);
+        goto DONE;
+    }
+
+    pmFPAview *view = pmFPAviewAlloc(0);
+    if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) {
+        psError(psErrorCodeLast(), false, "Unable to read files.");
+        pmFPAfileFreeSetStrict(false);
+        goto DONE;
+    }
+    pmChip *chip;
+    while ((chip = pmFPAviewNextChip (view, output->fpa, 1)) != NULL) {
+        psTrace ("pswarp", 4, "Chip %d: %x %x\n", view->chip, chip->file_exists, chip->process);
+        if (!chip->process || !chip->file_exists) { continue; }
+        if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) {
+            psError(psErrorCodeLast(), false, "Unable to read files.");
+            pmFPAfileFreeSetStrict(false);
+            goto DONE;
+        }
+        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 (!pmFPAfileIOChecks (config, view, PM_FPA_BEFORE) ||
+                !pmFPAfileIOChecks (config, view, PM_FPA_AFTER)) {
+                psError(psErrorCodeLast(), false, "Unable to read files.");
+                pmFPAfileFreeSetStrict(false);
+                goto DONE;
+            }
+        }
+        if (!pmFPAfileIOChecks (config, view, PM_FPA_AFTER)) {
+            psError(psErrorCodeLast(), false, "Unable to write files.");
+            pmFPAfileFreeSetStrict(false);
+            goto DONE;
+        }
+    }
+    if (!pmFPAfileIOChecks (config, view, PM_FPA_AFTER)) {
         psError(psErrorCodeLast(), false, "Unable to write files.");
-        exitValue = pswarpExitCode(exitValue);
         pmFPAfileFreeSetStrict(false);
+        goto DONE;
     }
+    psFree(view);
 
     if (!pswarpStatsFileSave (config, statsFile)) {
@@ -56 +96 @@
     psLibFinalize();
 
+DONE:
     exitValue = pswarpExitCode(exitValue);
     exit (exitValue);

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

@@ -53 +53 @@
 
     // generate Lmin,max, Mmin,max for both datasets
-    fprintf (stderr, "srcBounds:\n");
-    // pswarpBounds *srcBounds = pswarpMakeBounds (astrom->fpa, frame);
     pswarpBounds *srcBounds = pswarpMakeBounds (input->fpa, frame);
-
-    fprintf (stderr, "srcBounds:\n");
     pswarpBounds *tgtBounds = pswarpMakeBounds (skycell->fpa, frame);
 
-    pswarpFindOverlap (output->fpa, srcBounds, tgtBounds);
+    // find the output (tgt) chips which overlap the input (src) chips
+    pswarpFindOverlap (input->fpa, output->fpa, srcBounds, tgtBounds);
 
-    // XXX the loop below generates the output pixels. should this be deferred until 
-    // we actually need them?
+    // The loop below generates the output pixels. XXX 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);
+
+    // generate the output hdu and contents bits
+    pmFPAAddSourceFromFormat(output->fpa, output->format);
 
     pmChip *chip;
@@ -115 +116 @@
             cell->file_exists = true;
 
-# if (0)
-            // ******* this chunk comes from pswarpLoop, but should be done when the output chip is generated
+            // copy the basic headers across from astrom ref to output
             pmHDU *outHDU = pmHDUFromCell (cell);           ///< HDU for the output warped image
             outHDU->header = psMetadataCopy(outHDU->header, hdu->header);
             pswarpVersionHeader(outHDU->header);
-# endif
         }
     }
-
-    // XXX not sure what this does...
-    pmFPAviewReset (view);
-    pmFPAAddSourceFromView(output->fpa, view, output->format);
-
     psFree (view);
     return true;

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

@@ -36 +36 @@
     if (!pswarpSetMaskBits(config)) {
         psError(psErrorCodeLast(), false, "failed to set mask bits");
-        return NULL;
+        return false;
     }
 
@@ -83 +83 @@
     // XXX pmFPAfileActivate(config->files, true, "PSWARP.OUTPUT.SOURCES");
 
-    // Done with the detector side of things
+    // the loops below load the input data
     pswarpFileActivation(config, detectorFiles, true);
-    pswarpFileActivation(config, independentFiles, true);
+
+    // pswarpFileActivation(config, independentFiles, true);
 
     pmFPAview *view = pmFPAviewAlloc(0);
@@ -95 +96 @@
     if (!pmFPAfileIOChecks (config, view, PM_FPA_BEFORE)) {
         psError(psErrorCodeLast(), false, "Unable to read files.");
-        goto DONE;
+        goto FAIL;
     }
 
@@ -107 +108 @@
         if (!pmFPAfileIOChecks (config, view, PM_FPA_BEFORE)) {
             psError(psErrorCodeLast(), false, "Unable to read files.");
-            goto DONE;
+            goto FAIL;
         }
 
@@ -116 +117 @@
             if (!pmFPAfileIOChecks (config, view, PM_FPA_BEFORE)) {
                 psError(psErrorCodeLast(), false, "Unable to read files.");
-                goto DONE;
+                goto FAIL;
             }
 
@@ -124 +125 @@
                 if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) {
                     psError(psErrorCodeLast(), false, "Unable to read files.");
-                    goto DONE;
+                    goto FAIL;
                 }
                 if (!readout->data_exists) {
@@ -145 +146 @@
                 if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) {
                     psError(psErrorCodeLast(), false, "Unable to write files.");
-                    goto DONE;
+                    goto FAIL;
                 }
             }
             if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) {
                 psError(psErrorCodeLast(), false, "Unable to write files.");
-                goto DONE;
+                goto FAIL;
             }
         }
         if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) {
             psError(psErrorCodeLast(), false, "Unable to write files.");
-            goto DONE;
+            goto FAIL;
         }
     }
     if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) {
         psError(psErrorCodeLast(), false, "Unable to write files.");
-        goto DONE;
+        goto FAIL;
     }
 
@@ -167 +168 @@
     pswarpFileActivation(config, independentFiles, false);
 
-    // XXXX turn these on later .. pswarpUpdateStatistics (output);
+    if (!pswarpUpdateStatistics (output->fpa, stats, input->fpa, astrom->fpa, config)) {
+      psError(psErrorCodeLast(), false, "problem generating statistics.");
+      goto FAIL;
+    }
+
     // XXXX turn these on later .. pswarpMakePSF (output);
-
     // XXX have not yet written the output to disk
-
     psFree(view);
-
- DONE:
     return true;
+
+ FAIL:
+    psFree (view);
+    return false;
 }
 
-# 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");
+bool pswarpTransformToTarget (pmFPA *output, pmReadout *input, pmConfig *config)  {
 
     pmFPAview *view = pmFPAviewAlloc(0);
@@ -195 +195 @@
 
         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")) {
-        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);
-    }
-    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);

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

@@ -31 +31 @@
                               dynamicMaskVal,advisoryMaskVal)) {
     psError(PS_ERR_UNKNOWN, false, "Unable to calculate masks for readout.");
-    return(false);
+    return false;
   }
   psMetadataAddS32(stats, PS_LIST_TAIL,"MASKFRAC_NPIX", 0,
@@ -43 +43 @@
   psMetadataAddF32(stats,PS_LIST_TAIL, "MASKFRAC_ADVISORY", 0,
                    "Fraction of pixels masked as an advisory", (float) Npix_advisory / Npix_valid);
-  return(true);
+  return true;
 }

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

@@ -45 +45 @@
         if (!chip->process || !chip->file_exists) { 
             // this ensures one entry per chip (regardless of existence)
-            // XXX maybe not? pswarpBoundsAppend (bounds, NAN, NAN, NAN, NAN);
+            pswarpBoundsAppend (bounds, NAN, NAN, NAN, NAN);
             continue; 
         }
@@ -83 +83 @@
             }
         }
-        fprintf (stderr, "bounds %d : %f - %f, %f - %f\n", view->chip, Pmin, Pmax, Qmin, Qmax);
+        // fprintf (stderr, "bounds %d : %f - %f, %f - %f\n", view->chip, Pmin, Pmax, Qmin, Qmax);
         pswarpBoundsAppend (bounds, Pmin, Pmax, Qmin, Qmax);
     }
@@ -164 +164 @@
 }
 
-bool pswarpFindOverlap (pmFPA *fpa, pswarpBounds *src, pswarpBounds *tgt) {
+bool pswarpFindOverlap (pmFPA *input, pmFPA *output, pswarpBounds *src, pswarpBounds *tgt) {
 
     // we have the source and target bounds.  loop over all output elements and check if any of the source
     // elements overlap it.
 
-    assert (fpa->chips->n == tgt->Pmin->n);
-
-    for (int j = 0; j < fpa->chips->n; j++) {
-        pmChip *chip = fpa->chips->data[j];
+    // also add the list if input (src) cells which land on each output (tgt) chip we do this
+    // by only identifying the input and output chips and adding all input chip cells to each
+    // output cell
+
+    assert (output->chips->n == tgt->Pmin->n);
+    assert (input->chips->n == src->Pmin->n);
+
+    for (int j = 0; j < output->chips->n; j++) {
+        pmChip *chip = output->chips->data[j];
 
         // we have src bounds
@@ -180 +185 @@
         float Qmax = tgt->Qmax->data.F32[j];
 
-        bool valid = false;
-        for (int i = 0; !valid && (i < src->Pmin->n); i++) {
+        psArray *inputChips = psArrayAllocEmpty(8);
+
+        bool hasOverlap = false;
+        for (int i = 0; i < src->Pmin->n; i++) {
             // overlaps in P?
-            bool Pvalid = (Pmin < src->Pmax->data.F32[i]) && (Pmax > src->Pmin->data.F32[i]);
-            bool Qvalid = (Qmin < src->Qmax->data.F32[i]) && (Qmax > src->Qmin->data.F32[i]);
-            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;
+            if (!isfinite(src->Pmin->data.F32[i])) continue;
+            if (!isfinite(src->Pmax->data.F32[i])) continue;
+            if (!isfinite(src->Qmin->data.F32[i])) continue;
+            if (!isfinite(src->Qmax->data.F32[i])) continue;
+
+            if (Pmin > src->Pmax->data.F32[i]) continue;
+            if (Pmax < src->Pmin->data.F32[i]) continue;
+
+            if (Qmin > src->Qmax->data.F32[i]) continue;
+            if (Qmax < src->Qmin->data.F32[i]) continue;
+
+            hasOverlap = true;
+            psArrayAdd (inputChips, 1, input->chips->data[i]);
+        }
+        chip->process = hasOverlap;
+        if (hasOverlap) {
+            // get the complete list of input cells for this set of input chips
+            // pmConceptsAverageCells needs a psList (not a psArray)
+            psList *inputCells = psListAlloc(NULL);
+            for (int nChip = 0; nChip < inputChips->n; nChip++) {
+                pmChip *inputChip = inputChips->data[nChip];
+                for (int nCell = 0; nCell < inputChip->cells->n; nCell++) {
+                    fprintf (stderr, "input %d, %d : output %d\n", nChip, nCell, j);
+                    psListAdd (inputCells, PS_LIST_TAIL, inputChip->cells->data[nCell]);
+                }
+            }
+            
+            // place the inputCells on each of the output cells 
+            for (int nCell = 0; nCell < chip->cells->n; nCell++) {
+                pmCell *outputCell = chip->cells->data[nCell];
+                psMetadataAddPtr(outputCell->analysis, PS_LIST_TAIL, "INPUT.CELLS", PS_DATA_LIST , "input cells touching this output cell", inputCells);
+            }
+            psFree (inputCells);
+        }
+        psFree (inputChips);
     }
     return true;

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

@@ -64 +64 @@
     for (int y = 0; y < numRows; y++) {
         for (int x = 0; x < numCols; x++) {
-            if (!(mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & maskValue)) {
-                if (y > yMax) {
-                    yMax = y;
-                }
-                if (y < yMin) {
-                    yMin = y;
-                }
-                if (x > xMax) {
-                    xMax = x;
-                }
-                if (x < xMin) {
-                    xMin = x;
-                }
-            }
+            if (mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & maskValue) { continue; }
+            xMin = PS_MIN (xMin, x);
+            xMax = PS_MAX (xMax, x);
+            yMin = PS_MIN (yMin, y);
+            yMax = PS_MAX (yMax, y);
         }
     }
 
+    // XXX need to add these in a container appropriate to the element
     if (stats) {
         psMetadataAddS32(stats, PS_LIST_TAIL, "RANGE.XMIN", 0, "Minimum valid x value", xMin);

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

@@ -37 +37 @@
 
     pswarpStatsFile *statsFile = pswarpStatsFileAlloc ();
-    statsFile->name = name;
+    statsFile->name = psMemIncrRefCounter(name);
 
     if (!mdok) return statsFile; // XXX this is probably a config error, but treat as ok for now