Index: trunk/psModules/src/objects/pmDetections.c
===================================================================
--- trunk/psModules/src/objects/pmDetections.c	(revision 26450)
+++ trunk/psModules/src/objects/pmDetections.c	(revision 26893)
@@ -26,4 +26,8 @@
   psFree (detections->peaks);
   psFree (detections->oldPeaks);
+  psFree (detections->oldFootprints);
+
+  psFree (detections->newSources);
+  psFree (detections->allSources);
   return;
 }
@@ -35,8 +39,11 @@
     psMemSetDeallocator(detections, (psFreeFunc) pmDetectionsFree);
 
-    detections->footprints = NULL;
-    detections->peaks      = NULL;
-    detections->oldPeaks   = NULL;
-    detections->last       = 0;
+    detections->footprints    = NULL;
+    detections->peaks         = NULL;
+    detections->oldPeaks      = NULL;
+    detections->oldFootprints = NULL;
+    detections->newSources    = NULL;
+    detections->allSources    = NULL;
+    detections->last          = 0;
 
     return (detections);
Index: trunk/psModules/src/objects/pmDetections.h
===================================================================
--- trunk/psModules/src/objects/pmDetections.h	(revision 26450)
+++ trunk/psModules/src/objects/pmDetections.h	(revision 26893)
@@ -21,6 +21,9 @@
 typedef struct {
   psArray *footprints;        // collection of footprints in the image
+  psArray *oldFootprints;     // collection of footprints previously found
   psArray *peaks;             // collection of all peaks contained by the footprints
   psArray *oldPeaks;          // collection of all peaks previously found
+  psArray *newSources;        // collection of sources
+  psArray *allSources;        // collection of sources
   int last;
 } pmDetections;
Index: trunk/psModules/src/objects/pmFootprintCullPeaks.c
===================================================================
--- trunk/psModules/src/objects/pmFootprintCullPeaks.c	(revision 26450)
+++ trunk/psModules/src/objects/pmFootprintCullPeaks.c	(revision 26893)
@@ -29,11 +29,11 @@
   */
 
-# define IN_PEAK 1 
+# define IN_PEAK 1
 psErrorCode pmFootprintCullPeaks(const psImage *img, // the image wherein lives the footprint
-				 const psImage *weight,	// corresponding variance image
-				 pmFootprint *fp, // Footprint containing mortal peaks
-				 const float nsigma_delta, // how many sigma above local background a peak needs to be to survive
-				 const float fPadding, // fractional padding added to stdev since bright peaks have unreasonably high significance
-				 const float min_threshold) { // minimum permitted coll height
+                                 const psImage *weight, // corresponding variance image
+                                 pmFootprint *fp, // Footprint containing mortal peaks
+                                 const float nsigma_delta, // how many sigma above local background a peak needs to be to survive
+                                 const float fPadding, // fractional padding added to stdev since bright peaks have unreasonably high significance
+                                 const float min_threshold) { // minimum permitted coll height
     assert (img != NULL); assert (img->type.type == PS_TYPE_F32);
     assert (weight != NULL); assert (weight->type.type == PS_TYPE_F32);
@@ -42,8 +42,8 @@
 
     if (fp->peaks == NULL || fp->peaks->n == 0) { // nothing to do
-	return PS_ERR_NONE;
-    }
-
-    psRegion subRegion;			// desired subregion; 1 larger than bounding box (grr)
+        return PS_ERR_NONE;
+    }
+
+    psRegion subRegion;                 // desired subregion; 1 larger than bounding box (grr)
     subRegion.x0 = fp->bbox.x0; subRegion.x1 = fp->bbox.x1 + 1;
     subRegion.y0 = fp->bbox.y0; subRegion.y1 = fp->bbox.y1 + 1;
@@ -55,5 +55,5 @@
     psImage *idImg = psImageAlloc(subImg->numCols, subImg->numRows, PS_TYPE_S32);
 
-    // We need a psArray of peaks brighter than the current peak.  
+    // We need a psArray of peaks brighter than the current peak.
     // We reject peaks which either:
     // 1) are below the local threshold
@@ -67,68 +67,73 @@
     // The brightest peak is always safe; go through other peaks trying to cull them
     for (int i = 1; i < fp->peaks->n; i++) { // n.b. fp->peaks->n can change within the loop
-	const pmPeak *peak = fp->peaks->data[i];
-	int x = peak->x - subImg->col0;
-	int y = peak->y - subImg->row0;
-	//
-	// Find the level nsigma below the peak that must separate the peak
-	// from any of its friends
-	//
-	assert (x >= 0 && x < subImg->numCols && y >= 0 && y < subImg->numRows);
-
-	// const float stdev = sqrt(subWt->data.F32[y][x]);
-	// float threshold = subImg->data.F32[y][x] - nsigma_delta*stdev;
-
-	const float stdev = sqrt(subWt->data.F32[y][x]);
-	const float flux = subImg->data.F32[y][x];
-	const float fStdev = fabs(stdev/flux);
-	const float stdev_pad = fabs(flux * hypot(fStdev, fPadding));
-	// if flux is negative, careful with fStdev
-
-	float threshold = flux - nsigma_delta*stdev_pad;
-	
-	if (isnan(threshold) || threshold < min_threshold) {
-	    // min_threshold is assumed to be below the detection threshold,
-	    // so all the peaks are pmFootprint, and this isn't the brightest
+        const pmPeak *peak = fp->peaks->data[i];
+        int x = peak->x - subImg->col0;
+        int y = peak->y - subImg->row0;
+        //
+        // Find the level nsigma below the peak that must separate the peak
+        // from any of its friends
+        //
+        assert (x >= 0 && x < subImg->numCols && y >= 0 && y < subImg->numRows);
+
+        // const float stdev = sqrt(subWt->data.F32[y][x]);
+        // float threshold = subImg->data.F32[y][x] - nsigma_delta*stdev;
+
+        const float stdev = sqrt(subWt->data.F32[y][x]);
+        const float flux = subImg->data.F32[y][x];
+        const float fStdev = fabs(stdev/flux);
+        const float stdev_pad = fabs(flux * hypot(fStdev, fPadding));
+        // if flux is negative, careful with fStdev
+
+        float threshold = flux - nsigma_delta*stdev_pad;
+
+        if (isnan(threshold) || threshold < min_threshold) {
+            // min_threshold is assumed to be below the detection threshold,
+            // so all the peaks are pmFootprint, and this isn't the brightest
+            continue;
+        }
+
+        // XXX EAM : if stdev >= 0, i'm not sure how this can ever be true?
+        if (threshold > subImg->data.F32[y][x]) {
+            threshold = subImg->data.F32[y][x] - 10*FLT_EPSILON;
+        }
+
+        // XXX this is a bit expensive: psImageAlloc for every peak contained in this footprint
+        // perhaps this should alloc a single ID image above and pass it in to be set.
+
+        // XXX optionally use the faster pmFootprintsFind if the subimage size is large (eg, M31)
+
+        // XXX this is not quite there yet:
+        // pmFootprint *peakFootprint = NULL;
+        // int area = subImg->numCols * subImg->numRows;
+        // if (area > 30000) {
+
+        pmFootprint *peakFootprint = pmFootprintsFindAtPoint(subImg, threshold, brightPeaks, peak->y, peak->x);
+
+        // at this point brightPeaks only has the peaks brighter than the current
+
+        // XXX need to supply the image here
+        // we set the IDs to either 1 (in peak) or 0 (not in peak)
+        pmSetFootprintID (idImg, peakFootprint, IN_PEAK);
+        psFree(peakFootprint);
+
+        // If this peak has not already been assigned to a source, then we can look for any
+        // brighter peaks within its footprint. Check if any of the previous (brighter) peaks
+        // are within the footprint of this peak If so, the current peak is bogus; drop it.
+        bool keep = true;
+        for (int j = 0; keep && !peak->assigned && (j < brightPeaks->n); j++) {
+            const pmPeak *peak2 = fp->peaks->data[j];
+            int x2 = peak2->x - subImg->col0;
+            int y2 = peak2->y - subImg->row0;
+            if (idImg->data.S32[y2][x2] == IN_PEAK) {
+                // There's a brighter peak within the footprint above threshold; so cull our initial peak
+                keep = false;
+	    }
+        }
+        if (!keep) {
+	    psAssert (!peak->assigned, "logic error: trying to drop a previously-assigned peak");  // we should not drop any already assigned peaks.
 	    continue;
 	}
 
-	// XXX EAM : if stdev >= 0, i'm not sure how this can ever be true?
-	if (threshold > subImg->data.F32[y][x]) {
-	    threshold = subImg->data.F32[y][x] - 10*FLT_EPSILON;
-	}
-
-	// XXX this is a bit expensive: psImageAlloc for every peak contained in this footprint
-	// perhaps this should alloc a single ID image above and pass it in to be set.
-
-	// XXX optionally use the faster pmFootprintsFind if the subimage size is large (eg, M31)
-
-	// XXX this is not quite there yet:
-	// pmFootprint *peakFootprint = NULL;
-	// int area = subImg->numCols * subImg->numRows;
-	// if (area > 30000) {
-
-	pmFootprint *peakFootprint = pmFootprintsFindAtPoint(subImg, threshold, brightPeaks, peak->y, peak->x);
-
-	// at this point brightPeaks only has the peaks brighter than the current
-
-	// XXX need to supply the image here
-	// we set the IDs to either 1 (in peak) or 0 (not in peak)
-	pmSetFootprintID (idImg, peakFootprint, IN_PEAK);
-	psFree(peakFootprint);
-
-	// Check if any of the previous (brighter) peaks are within the footprint of this peak
-	// If so, the current peak is bogus; drop it.
-	bool keep = true;
-	for (int j = 0; keep && (j < brightPeaks->n); j++) {
-	    const pmPeak *peak2 = fp->peaks->data[j];
-	    int x2 = peak2->x - subImg->col0;
-	    int y2 = peak2->y - subImg->row0;
-	    if (idImg->data.S32[y2][x2] == IN_PEAK) 
-		// There's a brighter peak within the footprint above threshold; so cull our initial peak
-		keep = false;
-	}
-	if (!keep) continue;
-
-	psArrayAdd (brightPeaks, 128, fp->peaks->data[i]);
+        psArrayAdd (brightPeaks, 128, fp->peaks->data[i]);
     }
 
@@ -151,9 +156,9 @@
   */
 psErrorCode pmFootprintCullPeaks_OLD(const psImage *img, // the image wherein lives the footprint
-				 const psImage *weight,	// corresponding variance image
-				 pmFootprint *fp, // Footprint containing mortal peaks
-				 const float nsigma_delta, // how many sigma above local background a peak needs to be to survive
-				 const float fPadding, // fractional padding added to stdev since bright peaks have unreasonably high significance
-				 const float min_threshold) { // minimum permitted coll height
+                                 const psImage *weight, // corresponding variance image
+                                 pmFootprint *fp, // Footprint containing mortal peaks
+                                 const float nsigma_delta, // how many sigma above local background a peak needs to be to survive
+                                 const float fPadding, // fractional padding added to stdev since bright peaks have unreasonably high significance
+                                 const float min_threshold) { // minimum permitted coll height
     assert (img != NULL); assert (img->type.type == PS_TYPE_F32);
     assert (weight != NULL); assert (weight->type.type == PS_TYPE_F32);
@@ -162,8 +167,8 @@
 
     if (fp->peaks == NULL || fp->peaks->n == 0) { // nothing to do
-	return PS_ERR_NONE;
-    }
-
-    psRegion subRegion;			// desired subregion; 1 larger than bounding box (grr)
+        return PS_ERR_NONE;
+    }
+
+    psRegion subRegion;                 // desired subregion; 1 larger than bounding box (grr)
     subRegion.x0 = fp->bbox.x0; subRegion.x1 = fp->bbox.x1 + 1;
     subRegion.y0 = fp->bbox.y0; subRegion.y1 = fp->bbox.y1 + 1;
@@ -185,79 +190,79 @@
     //
     for (int i = 1; i < fp->peaks->n; i++) { // n.b. fp->peaks->n can change within the loop
-	const pmPeak *peak = fp->peaks->data[i];
-	int x = peak->x - subImg->col0;
-	int y = peak->y - subImg->row0;
-	//
-	// Find the level nsigma below the peak that must separate the peak
-	// from any of its friends
-	//
-	assert (x >= 0 && x < subImg->numCols && y >= 0 && y < subImg->numRows);
-
-	// stdev ~ sqrt(flux) : for bright regions (f ~ 1e6, sqrt(f) = 1e3 -- unrealistically small deviations)
-	// add additional padding in quadrature:
-
-	const float stdev = sqrt(subWt->data.F32[y][x]);
-	const float flux = subImg->data.F32[y][x];
-	const float fStdev = stdev/flux;
-	const float stdev_pad = flux * hypot(fStdev, fPadding);
-	float threshold = flux - nsigma_delta*stdev_pad;
-
-	if (isnan(threshold) || threshold < min_threshold) {
-#if 1	  // min_threshold is assumed to be below the detection threshold,
-	  // so all the peaks are pmFootprint, and this isn't the brightest
-	    // XXX mark peak to be dropped
-	    (void)psArrayRemoveIndex(fp->peaks, i);
-	    i--;			// we moved everything down one
-	    continue;
+        const pmPeak *peak = fp->peaks->data[i];
+        int x = peak->x - subImg->col0;
+        int y = peak->y - subImg->row0;
+        //
+        // Find the level nsigma below the peak that must separate the peak
+        // from any of its friends
+        //
+        assert (x >= 0 && x < subImg->numCols && y >= 0 && y < subImg->numRows);
+
+        // stdev ~ sqrt(flux) : for bright regions (f ~ 1e6, sqrt(f) = 1e3 -- unrealistically small deviations)
+        // add additional padding in quadrature:
+
+        const float stdev = sqrt(subWt->data.F32[y][x]);
+        const float flux = subImg->data.F32[y][x];
+        const float fStdev = stdev/flux;
+        const float stdev_pad = flux * hypot(fStdev, fPadding);
+        float threshold = flux - nsigma_delta*stdev_pad;
+
+        if (isnan(threshold) || threshold < min_threshold) {
+#if 1     // min_threshold is assumed to be below the detection threshold,
+          // so all the peaks are pmFootprint, and this isn't the brightest
+            // XXX mark peak to be dropped
+            (void)psArrayRemoveIndex(fp->peaks, i);
+            i--;                        // we moved everything down one
+            continue;
 #else
 #error n.b. We will be running LOTS of checks at this threshold, so only find the footprint once
-	    threshold = min_threshold;
+            threshold = min_threshold;
 #endif
-	}
-
-	// XXX EAM : if stdev >= 0, i'm not sure how this can ever be true?
-	if (threshold > subImg->data.F32[y][x]) {
-	    threshold = subImg->data.F32[y][x] - 10*FLT_EPSILON;
-	}
-
-	// XXX this is a bit expensive: psImageAlloc for every peak contained in this footprint
-	// perhaps this should alloc a single ID image above and pass it in to be set.
-
-	const int peak_id = 1;		// the ID for the peak of interest
-	brightPeaks->n = i;		// only stop at a peak brighter than we are
-
-	// XXX optionally use the faster pmFootprintsFind if the subimage size is large (eg, M31)
-
-	pmFootprint *peakFootprint = pmFootprintsFindAtPoint(subImg, threshold, brightPeaks, peak->y, peak->x);
-	brightPeaks->n = 0;		// don't double free
-	psImage *idImg = pmSetFootprintID(NULL, peakFootprint, peak_id);
-	psFree(peakFootprint);
-
-	// Check if any of the previous (brighter) peaks are within the footprint of this peak
-	// If so, the current peak is bogus; drop it.
-	int j;
-	for (j = 0; j < i; j++) {
-	    const pmPeak *peak2 = fp->peaks->data[j];
-	    int x2 = peak2->x - subImg->col0;
-	    int y2 = peak2->y - subImg->row0;
-	    const int peak2_id = idImg->data.S32[y2][x2]; // the ID for some other peak
-
-	    if (peak2_id == peak_id) {	// There's a brighter peak within the footprint above
-		;			// threshold; so cull our initial peak
-		(void)psArrayRemoveIndex(fp->peaks, i);
-		i--;			// we moved everything down one
-		break;
-	    }
-	}
-	if (j == i) {
-	    j++;
-	}
-
-	psFree(idImg);
+        }
+
+        // XXX EAM : if stdev >= 0, i'm not sure how this can ever be true?
+        if (threshold > subImg->data.F32[y][x]) {
+            threshold = subImg->data.F32[y][x] - 10*FLT_EPSILON;
+        }
+
+        // XXX this is a bit expensive: psImageAlloc for every peak contained in this footprint
+        // perhaps this should alloc a single ID image above and pass it in to be set.
+
+        const int peak_id = 1;          // the ID for the peak of interest
+        brightPeaks->n = i;             // only stop at a peak brighter than we are
+
+        // XXX optionally use the faster pmFootprintsFind if the subimage size is large (eg, M31)
+
+        pmFootprint *peakFootprint = pmFootprintsFindAtPoint(subImg, threshold, brightPeaks, peak->y, peak->x);
+        brightPeaks->n = 0;             // don't double free
+        psImage *idImg = pmSetFootprintID(NULL, peakFootprint, peak_id);
+        psFree(peakFootprint);
+
+        // Check if any of the previous (brighter) peaks are within the footprint of this peak
+        // If so, the current peak is bogus; drop it.
+        int j;
+        for (j = 0; j < i; j++) {
+            const pmPeak *peak2 = fp->peaks->data[j];
+            int x2 = peak2->x - subImg->col0;
+            int y2 = peak2->y - subImg->row0;
+            const int peak2_id = idImg->data.S32[y2][x2]; // the ID for some other peak
+
+            if (peak2_id == peak_id) {  // There's a brighter peak within the footprint above
+                ;                       // threshold; so cull our initial peak
+                (void)psArrayRemoveIndex(fp->peaks, i);
+                i--;                    // we moved everything down one
+                break;
+            }
+        }
+        if (j == i) {
+            j++;
+        }
+
+        psFree(idImg);
     }
 
     brightPeaks->n = 0; psFree(brightPeaks);
-    psFree((psImage *)subImg);
-    psFree((psImage *)subWt);
+    psFree(subImg);
+    psFree(subWt);
 
     return PS_ERR_NONE;
Index: trunk/psModules/src/objects/pmFootprintFindAtPoint.c
===================================================================
--- trunk/psModules/src/objects/pmFootprintFindAtPoint.c	(revision 26450)
+++ trunk/psModules/src/objects/pmFootprintFindAtPoint.c	(revision 26893)
@@ -29,5 +29,5 @@
  * so we set appropriate mask bits
  *
- * EAM : these function were confusingly using "startspan" and "spartspan"  
+ * EAM : these function were confusingly using "startspan" and "spartspan"
  * I've rationalized them all to 'startspan'
  */
@@ -36,18 +36,18 @@
 // An enum for what we should do with a Startspan
 //
-typedef enum {PM_SSPAN_DOWN = 0,	// scan down from this span
-	      PM_SSPAN_UP,		// scan up from this span
-	      PM_SSPAN_RESTART,		// restart scanning from this span
-	      PM_SSPAN_DONE		// this span is processed
-} PM_SSPAN_DIR;				// How to continue searching
+typedef enum {PM_SSPAN_DOWN = 0,        // scan down from this span
+              PM_SSPAN_UP,              // scan up from this span
+              PM_SSPAN_RESTART,         // restart scanning from this span
+              PM_SSPAN_DONE             // this span is processed
+} PM_SSPAN_DIR;                         // How to continue searching
 //
 // An enum for mask's pixel values.  We're looking for pixels that are above threshold, and
 // we keep extra book-keeping information in the PM_SSPAN_STOP plane.  It's simpler to be
-// able to check for 
+// able to check for
 //
 enum {
-    PM_SSPAN_INITIAL = 0x0,		// initial state of pixels.
-    PM_SSPAN_DETECTED = 0x1,		// we've seen this pixel
-    PM_SSPAN_STOP = 0x2			// you may stop searching when you see this pixel
+    PM_SSPAN_INITIAL = 0x0,             // initial state of pixels.
+    PM_SSPAN_DETECTED = 0x1,            // we've seen this pixel
+    PM_SSPAN_STOP = 0x2                 // you may stop searching when you see this pixel
 };
 //
@@ -55,17 +55,17 @@
 //
 typedef struct {
-    const pmSpan *span;			// save the pixel range
-    PM_SSPAN_DIR direction;		// How to continue searching
-    bool stop;				// should we stop searching?
+    const pmSpan *span;                 // save the pixel range
+    PM_SSPAN_DIR direction;             // How to continue searching
+    bool stop;                          // should we stop searching?
 } Startspan;
 
 static void startspanFree(Startspan *sspan) {
-    psFree((void *)sspan->span);
+    psFree(sspan->span);
 }
 
 static Startspan *
-StartspanAlloc(const pmSpan *span,	// The span in question
-	       psImage *mask,		// Pixels that we've already detected
-	       const PM_SSPAN_DIR dir	// Should we continue searching towards the top of the image?
+StartspanAlloc(const pmSpan *span,      // The span in question
+               psImage *mask,           // Pixels that we've already detected
+               const PM_SSPAN_DIR dir   // Should we continue searching towards the top of the image?
     ) {
     Startspan *sspan = psAlloc(sizeof(Startspan));
@@ -75,16 +75,16 @@
     sspan->direction = dir;
     sspan->stop = false;
-    
-    if (mask != NULL) {			// remember that we've detected these pixels
-	psImageMaskType *mpix = &mask->data.PS_TYPE_IMAGE_MASK_DATA[span->y - mask->row0][span->x0 - mask->col0];
-
-	for (int i = 0; i <= span->x1 - span->x0; i++) {
-	    mpix[i] |= PM_SSPAN_DETECTED;
-	    if (mpix[i] & PM_SSPAN_STOP) {
-		sspan->stop = true;
-	    }
-	}
-    }
-    
+
+    if (mask != NULL) {                 // remember that we've detected these pixels
+        psImageMaskType *mpix = &mask->data.PS_TYPE_IMAGE_MASK_DATA[span->y - mask->row0][span->x0 - mask->col0];
+
+        for (int i = 0; i <= span->x1 - span->x0; i++) {
+            mpix[i] |= PM_SSPAN_DETECTED;
+            if (mpix[i] & PM_SSPAN_STOP) {
+                sspan->stop = true;
+            }
+        }
+    }
+
     return sspan;
 }
@@ -94,22 +94,22 @@
 //
 static bool add_startspan(psArray *startspans, // the saved Startspans
-			  const pmSpan *sp, // the span in question
-			  psImage *mask, // mask of detected/stop pixels
-			  const PM_SSPAN_DIR dir) { // the desired direction to search
+                          const pmSpan *sp, // the span in question
+                          psImage *mask, // mask of detected/stop pixels
+                          const PM_SSPAN_DIR dir) { // the desired direction to search
     if (dir == PM_SSPAN_RESTART) {
-	if (add_startspan(startspans, sp, mask,  PM_SSPAN_UP) ||
-	    add_startspan(startspans, sp, NULL, PM_SSPAN_DOWN)) {
-	    return true;
-	}
+        if (add_startspan(startspans, sp, mask,  PM_SSPAN_UP) ||
+            add_startspan(startspans, sp, NULL, PM_SSPAN_DOWN)) {
+            return true;
+        }
     } else {
-	Startspan *sspan = StartspanAlloc(sp, mask, dir);
-	if (sspan->stop) {		// we detected a stop bit
-	    psFree(sspan);		// don't allocate new span
-
-	    return true;
-	} else {
-	    psArrayAdd(startspans, 1, sspan);
-	    psFree(sspan);		// as it's now owned by startspans
-	}
+        Startspan *sspan = StartspanAlloc(sp, mask, dir);
+        if (sspan->stop) {              // we detected a stop bit
+            psFree(sspan);              // don't allocate new span
+
+            return true;
+        } else {
+            psArrayAdd(startspans, 1, sspan);
+            psFree(sspan);              // as it's now owned by startspans
+        }
     }
 
@@ -127,44 +127,44 @@
  */
 static bool do_startspan(pmFootprint *fp, // the footprint that we're building
-			 const psImage *img, // the psImage we're working on
-			 psImage *mask, // the associated masks
-			 const float threshold,	// Threshold
-			 psArray *startspans) {	// specify which span to process next
-    bool F32 = false;			// is this an F32 image?
+                         const psImage *img, // the psImage we're working on
+                         psImage *mask, // the associated masks
+                         const float threshold, // Threshold
+                         psArray *startspans) { // specify which span to process next
+    bool F32 = false;                   // is this an F32 image?
     if (img->type.type == PS_TYPE_F32) {
-	F32 = true;
+        F32 = true;
     } else if (img->type.type == PS_TYPE_S32) {
-	F32 = false;
-    } else {				// N.b. You can't trivially add more cases here; F32 is just a bool
-	psError(PS_ERR_UNKNOWN, true, "Unsupported psImage type: %d", img->type.type);
-	return NULL;
-    }
-
-    psF32 *imgRowF32 = NULL;		// row pointer if F32
-    psS32 *imgRowS32 = NULL;		//  "   "   "  "  !F32
-    psImageMaskType *maskRow = NULL;		//  masks's row pointer
-    
+        F32 = false;
+    } else {                            // N.b. You can't trivially add more cases here; F32 is just a bool
+        psError(PS_ERR_UNKNOWN, true, "Unsupported psImage type: %d", img->type.type);
+        return NULL;
+    }
+
+    psF32 *imgRowF32 = NULL;            // row pointer if F32
+    psS32 *imgRowS32 = NULL;            //  "   "   "  "  !F32
+    psImageMaskType *maskRow = NULL;            //  masks's row pointer
+
     const int row0 = img->row0;
     const int col0 = img->col0;
     const int numRows = img->numRows;
     const int numCols = img->numCols;
-    
+
     /********************************************************************************************************/
-    
+
     Startspan *sspan = NULL;
     for (int i = 0; i < startspans->n; i++) {
-	sspan = startspans->data[i];
-	if (sspan->direction != PM_SSPAN_DONE) {
-	    break;
-	}
-	if (sspan->stop) {
-	    break;
-	}
+        sspan = startspans->data[i];
+        if (sspan->direction != PM_SSPAN_DONE) {
+            break;
+        }
+        if (sspan->stop) {
+            break;
+        }
     }
     if (sspan == NULL || sspan->direction == PM_SSPAN_DONE) { // no more Startspans to process
-	return false;
-    }
-    if (sspan->stop) {			// they don't want any more spans processed
-	return false;
+        return false;
+    }
+    if (sspan->stop) {                  // they don't want any more spans processed
+        return false;
     }
     /*
@@ -179,111 +179,111 @@
      * Go through image identifying objects
      */
-    int nx0, nx1 = -1;			// new values of x0, x1
+    int nx0, nx1 = -1;                  // new values of x0, x1
     const int di = (dir == PM_SSPAN_UP) ? 1 : -1; // how much i changes to get to the next row
-    bool stop = false;			// should I stop searching for spans?
+    bool stop = false;                  // should I stop searching for spans?
 
     for (int i = sspan->span->y -row0 + di; i < numRows && i >= 0; i += di) {
-	imgRowF32 = img->data.F32[i];	// only one of
-	imgRowS32 = img->data.S32[i];	//      these is valid!
-	maskRow = mask->data.PS_TYPE_IMAGE_MASK_DATA[i];
-	//
-	// Search left from the pixel diagonally to the left of (i - di, x0). If there's
-	// a connected span there it may need to grow up and/or down, so push it onto
-	// the stack for later consideration
-	//
-	nx0 = -1;
-	for (int j = x0 - 1; j >= -1; j--) {
-	    double pixVal = (j < 0) ? threshold - 100 : (F32 ? imgRowF32[j] : imgRowS32[j]);
-	    if ((maskRow[j] & PM_SSPAN_DETECTED) || pixVal < threshold) {
-		if (j < x0 - 1) {	// we found some pixels above threshold
-		    nx0 = j + 1;
-		}
-		break;
-	    }
-	}
-
-	if (nx0 < 0) {			// no span to the left
-	    nx1 = x0 - 1;		// we're going to resume searching at nx1 + 1
-	} else {
-	    //
-	    // Search right in leftmost span
-	    //
-	    //nx1 = 0;			// make gcc happy
-	    for (int j = nx0 + 1; j <= numCols; j++) {
-		double pixVal = (j >= numCols) ? threshold - 100 : (F32 ? imgRowF32[j] : imgRowS32[j]);
-		if ((maskRow[j] & PM_SSPAN_DETECTED) || pixVal < threshold) {
-		    nx1 = j - 1;
-		    break;
-		}
-	    }
-	    
-	    const pmSpan *sp = pmFootprintAddSpan(fp, i + row0, nx0 + col0, nx1 + col0);
-	    
-	    if (add_startspan(startspans, sp, mask, PM_SSPAN_RESTART)) {
-		stop = true;
-		break;
-	    }
-	}
-	//
-	// Now look for spans connected to the old span.  The first of these we'll
-	// simply process, but others will have to be deferred for later consideration.
-	//
-	// In fact, if the span overhangs to the right we'll have to defer the overhang
-	// until later too, as it too can grow in both directions
-	//
-	// Note that column numCols exists virtually, and always ends the last span; this
-	// is why we claim below that sx1 is always set
-	//
-	bool first = false;		// is this the first new span detected?
-	for (int j = nx1 + 1; j <= x1 + 1; j++) {
-	    double pixVal = (j >= numCols) ? threshold - 100 : (F32 ? imgRowF32[j] : imgRowS32[j]);
-	    if (!(maskRow[j] & PM_SSPAN_DETECTED) && pixVal >= threshold) {
-		int sx0 = j++;		// span that we're working on is sx0:sx1
-		int sx1 = -1;		// We know that if we got here, we'll also set sx1
-		for (; j <= numCols; j++) {
-		    double pixVal = (j >= numCols) ? threshold - 100 : (F32 ? imgRowF32[j] : imgRowS32[j]);
-		    if ((maskRow[j] & PM_SSPAN_DETECTED) || pixVal < threshold) { // end of span
-			sx1 = j;
-			break;
-		    }
-		}
-		assert (sx1 >= 0);
-
-		const pmSpan *sp;
-		if (first) {
-		    if (sx1 <= x1) {
-			sp = pmFootprintAddSpan(fp, i + row0, sx0 + col0, sx1 + col0 - 1);
-			if (add_startspan(startspans, sp, mask, PM_SSPAN_DONE)) {
-			    stop = true;
-			    break;
-			}
-		    } else {		// overhangs to right
-			sp = pmFootprintAddSpan(fp, i + row0, sx0 + col0, x1 + col0);
-			if (add_startspan(startspans, sp, mask, PM_SSPAN_DONE)) {
-			    stop = true;
-			    break;
-			}
-			sp = pmFootprintAddSpan(fp, i + row0, x1 + 1 + col0, sx1 + col0 - 1);
-			if (add_startspan(startspans, sp, mask, PM_SSPAN_RESTART)) {
-			    stop = true;
-			    break;
-			}
-		    }
-		    first = false;
-		} else {
-		    sp = pmFootprintAddSpan(fp, i + row0, sx0 + col0, sx1 + col0 - 1);
-		    if (add_startspan(startspans, sp, mask, PM_SSPAN_RESTART)) {
-			stop = true;
-			break;
-		    }
-		}
-	    }
-	}
-
-	if (stop || first == false) {	// we're done
-	    break;
-	}
-
-	x0 = nx0; x1 = nx1;
+        imgRowF32 = img->data.F32[i];   // only one of
+        imgRowS32 = img->data.S32[i];   //      these is valid!
+        maskRow = mask->data.PS_TYPE_IMAGE_MASK_DATA[i];
+        //
+        // Search left from the pixel diagonally to the left of (i - di, x0). If there's
+        // a connected span there it may need to grow up and/or down, so push it onto
+        // the stack for later consideration
+        //
+        nx0 = -1;
+        for (int j = x0 - 1; j >= -1; j--) {
+            double pixVal = (j < 0) ? threshold - 100 : (F32 ? imgRowF32[j] : imgRowS32[j]);
+            if ((maskRow[j] & PM_SSPAN_DETECTED) || pixVal < threshold) {
+                if (j < x0 - 1) {       // we found some pixels above threshold
+                    nx0 = j + 1;
+                }
+                break;
+            }
+        }
+
+        if (nx0 < 0) {                  // no span to the left
+            nx1 = x0 - 1;               // we're going to resume searching at nx1 + 1
+        } else {
+            //
+            // Search right in leftmost span
+            //
+            //nx1 = 0;                  // make gcc happy
+            for (int j = nx0 + 1; j <= numCols; j++) {
+                double pixVal = (j >= numCols) ? threshold - 100 : (F32 ? imgRowF32[j] : imgRowS32[j]);
+                if ((maskRow[j] & PM_SSPAN_DETECTED) || pixVal < threshold) {
+                    nx1 = j - 1;
+                    break;
+                }
+            }
+
+            const pmSpan *sp = pmFootprintAddSpan(fp, i + row0, nx0 + col0, nx1 + col0);
+
+            if (add_startspan(startspans, sp, mask, PM_SSPAN_RESTART)) {
+                stop = true;
+                break;
+            }
+        }
+        //
+        // Now look for spans connected to the old span.  The first of these we'll
+        // simply process, but others will have to be deferred for later consideration.
+        //
+        // In fact, if the span overhangs to the right we'll have to defer the overhang
+        // until later too, as it too can grow in both directions
+        //
+        // Note that column numCols exists virtually, and always ends the last span; this
+        // is why we claim below that sx1 is always set
+        //
+        bool first = false;             // is this the first new span detected?
+        for (int j = nx1 + 1; j <= x1 + 1; j++) {
+            double pixVal = (j >= numCols) ? threshold - 100 : (F32 ? imgRowF32[j] : imgRowS32[j]);
+            if (!(maskRow[j] & PM_SSPAN_DETECTED) && pixVal >= threshold) {
+                int sx0 = j++;          // span that we're working on is sx0:sx1
+                int sx1 = -1;           // We know that if we got here, we'll also set sx1
+                for (; j <= numCols; j++) {
+                    double pixVal = (j >= numCols) ? threshold - 100 : (F32 ? imgRowF32[j] : imgRowS32[j]);
+                    if ((maskRow[j] & PM_SSPAN_DETECTED) || pixVal < threshold) { // end of span
+                        sx1 = j;
+                        break;
+                    }
+                }
+                assert (sx1 >= 0);
+
+                const pmSpan *sp;
+                if (first) {
+                    if (sx1 <= x1) {
+                        sp = pmFootprintAddSpan(fp, i + row0, sx0 + col0, sx1 + col0 - 1);
+                        if (add_startspan(startspans, sp, mask, PM_SSPAN_DONE)) {
+                            stop = true;
+                            break;
+                        }
+                    } else {            // overhangs to right
+                        sp = pmFootprintAddSpan(fp, i + row0, sx0 + col0, x1 + col0);
+                        if (add_startspan(startspans, sp, mask, PM_SSPAN_DONE)) {
+                            stop = true;
+                            break;
+                        }
+                        sp = pmFootprintAddSpan(fp, i + row0, x1 + 1 + col0, sx1 + col0 - 1);
+                        if (add_startspan(startspans, sp, mask, PM_SSPAN_RESTART)) {
+                            stop = true;
+                            break;
+                        }
+                    }
+                    first = false;
+                } else {
+                    sp = pmFootprintAddSpan(fp, i + row0, sx0 + col0, sx1 + col0 - 1);
+                    if (add_startspan(startspans, sp, mask, PM_SSPAN_RESTART)) {
+                        stop = true;
+                        break;
+                    }
+                }
+            }
+        }
+
+        if (stop || first == false) {   // we're done
+            break;
+        }
+
+        x0 = nx0; x1 = nx1;
     }
     /*
@@ -309,22 +309,22 @@
  */
 pmFootprint *
-pmFootprintsFindAtPoint(const psImage *img,	// image to search
-		       const float threshold,	// Threshold
-		       const psArray *peaks, // array of peaks; finding one terminates search for footprint
-		       int row, int col) { // starting position (in img's parent's coordinate system)
+pmFootprintsFindAtPoint(const psImage *img,     // image to search
+                       const float threshold,   // Threshold
+                       const psArray *peaks, // array of peaks; finding one terminates search for footprint
+                       int row, int col) { // starting position (in img's parent's coordinate system)
    assert(img != NULL);
 
-   bool F32 = false;			// is this an F32 image?
+   bool F32 = false;                    // is this an F32 image?
    if (img->type.type == PS_TYPE_F32) {
        F32 = true;
    } else if (img->type.type == PS_TYPE_S32) {
        F32 = false;
-   } else {				// N.b. You can't trivially add more cases here; F32 is just a bool
+   } else {                             // N.b. You can't trivially add more cases here; F32 is just a bool
        psError(PS_ERR_UNKNOWN, true, "Unsupported psImage type: %d", img->type.type);
        return NULL;
    }
-   psF32 *imgRowF32 = NULL;		// row pointer if F32
-   psS32 *imgRowS32 = NULL;		//  "   "   "  "  !F32
-   
+   psF32 *imgRowF32 = NULL;             // row pointer if F32
+   psS32 *imgRowS32 = NULL;             //  "   "   "  "  !F32
+
    const int row0 = img->row0;
    const int col0 = img->col0;
@@ -339,5 +339,5 @@
         psError(PS_ERR_BAD_PARAMETER_VALUE, true,
                 "row/col == (%d, %d) are out of bounds [%d--%d, %d--%d]",
-		row + row0, col + col0, row0, row0 + numRows - 1, col0, col0 + numCols - 1);
+                row + row0, col + col0, row0, row0 + numRows - 1, col0, col0 + numCols - 1);
        return NULL;
    }
@@ -347,5 +347,5 @@
        return pmFootprintAlloc(0, img);
    }
-   
+
    pmFootprint *fp = pmFootprintAlloc(1 + img->numRows/10, img);
 /*
@@ -364,6 +364,6 @@
    if (peaks != NULL) {
        for (int i = 0; i < peaks->n; i++) {
-	   pmPeak *peak = peaks->data[i];
-	   mask->data.PS_TYPE_IMAGE_MASK_DATA[peak->y - mask->row0][peak->x - mask->col0] |= PM_SSPAN_STOP;
+           pmPeak *peak = peaks->data[i];
+           mask->data.PS_TYPE_IMAGE_MASK_DATA[peak->y - mask->row0][peak->x - mask->col0] |= PM_SSPAN_STOP;
        }
    }
@@ -372,22 +372,22 @@
  */
    psArray *startspans = psArrayAllocEmpty(1); // spans where we have to restart the search
-   
-   imgRowF32 = img->data.F32[row];	// only one of
-   imgRowS32 = img->data.S32[row];	//      these is valid!
+
+   imgRowF32 = img->data.F32[row];      // only one of
+   imgRowS32 = img->data.S32[row];      //      these is valid!
    psImageMaskType *maskRow = mask->data.PS_TYPE_IMAGE_MASK_DATA[row];
    {
        int i;
        for (i = col; i >= 0; i--) {
-	   pixVal = F32 ? imgRowF32[i] : imgRowS32[i];
-	   if ((maskRow[i] & PM_SSPAN_DETECTED) || pixVal < threshold) {
-	       break;
-	   }
+           pixVal = F32 ? imgRowF32[i] : imgRowS32[i];
+           if ((maskRow[i] & PM_SSPAN_DETECTED) || pixVal < threshold) {
+               break;
+           }
        }
        int i0 = i;
        for (i = col; i < numCols; i++) {
-	   pixVal = F32 ? imgRowF32[i] : imgRowS32[i];
-	   if ((maskRow[i] & PM_SSPAN_DETECTED) || pixVal < threshold) {
-	       break;
-	   }
+           pixVal = F32 ? imgRowF32[i] : imgRowS32[i];
+           if ((maskRow[i] & PM_SSPAN_DETECTED) || pixVal < threshold) {
+               break;
+           }
        }
        int i1 = i;
@@ -404,6 +404,6 @@
     */
    psFree(mask);
-   psFree(startspans);			// restores the image pixel
-
-   return fp;				// pmFootprint really
+   psFree(startspans);                  // restores the image pixel
+
+   return fp;                           // pmFootprint really
 }
Index: trunk/psModules/src/objects/pmObjects.c
===================================================================
--- trunk/psModules/src/objects/pmObjects.c	(revision 26450)
+++ 	(revision )
@@ -1,25 +1,0 @@
-/** @file  pmObjects.c
- *
- *  This file will ...
- *
- *  @author GLG, MHPCC
- *  @author EAM, IfA: significant modifications.
- *
- *  @version $Revision: 1.13 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2006-09-15 09:49:01 $
- *
- *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
- *
- */
-
-#ifdef HAVE_CONFIG_H
-#include <config.h>
-#endif
-
-#include <stdio.h>
-#include <math.h>
-#include <string.h>
-#include <pslib.h>
-#include "pmObjects.h"
-#include "pmModelGroup.h"
-
Index: trunk/psModules/src/objects/pmObjects.h
===================================================================
--- trunk/psModules/src/objects/pmObjects.h	(revision 26450)
+++ 	(revision )
@@ -1,68 +1,0 @@
-/* @file  pmObjects.h
- *
- * The process of finding, measuring, and classifying astronomical sources on
- * images is one of the critical tasks of the IPP or any astronomical software
- * system. This file will define structures and functions related to the task
- * of source detection and measurement. The elements defined in this section
- * are generally low-level components which can be connected together to
- * construct a complete object measurement suite.
- *
- * @author GLG, MHPCC
- *
- * @version $Revision: 1.8 $ $Name: not supported by cvs2svn $
- * @date $Date: 2007-01-24 02:54:15 $
- * Copyright 2004 Maui High Performance Computing Center, University of Hawaii
- */
-
-#ifndef PM_OBJECTS_H
-#define PM_OBJECTS_H
-
-/// @addtogroup Objects Object Detection / Analysis Functions
-/// @{
-
-#include <stdio.h>
-#include <math.h>
-#include <pslib.h>
-
-/**
- *
- * This function converts the source classification into the closest available
- * approximation to the Dophot classification scheme:
- * XXX EAM : fix this to use current source classification scheme
- *
- * PM_SOURCE_DEFECT: 8
- * PM_SOURCE_SATURATED: 8
- * PM_SOURCE_SATSTAR: 10
- * PM_SOURCE_PSFSTAR: 1
- * PM_SOURCE_GOODSTAR: 1
- * PM_SOURCE_POOR_FIT_PSF: 7
- * PM_SOURCE_FAIL_FIT_PSF: 4
- * PM_SOURCE_FAINTSTAR: 4
- * PM_SOURCE_GALAXY: 2
- * PM_SOURCE_FAINT_GALAXY: 2
- * PM_SOURCE_DROP_GALAXY: 2
- * PM_SOURCE_FAIL_FIT_GAL: 2
- * PM_SOURCE_POOR_FIT_GAL: 2
- * PM_SOURCE_OTHER: ?
- *
- */
-int pmSourceDophotType(
-    pmSource *source                    ///< Add comment.
-);
-
-
-/** pmSourceSextractType()
- *
- * This function converts the source classification into the closest available
- * approximation to the Sextractor classification scheme. the correspondence is
- * not yet defined (TBD) .
- *
- * XXX: Must code this.
- *
- */
-int pmSourceSextractType(
-    pmSource *source                    ///< Add comment.
-);
-
-/// @}
-#endif
Index: trunk/psModules/src/objects/pmPSF.c
===================================================================
--- trunk/psModules/src/objects/pmPSF.c	(revision 26450)
+++ trunk/psModules/src/objects/pmPSF.c	(revision 26893)
@@ -42,4 +42,7 @@
 #include "pmErrorCodes.h"
 
+
+#define MAX_AXIS_RATIO 20.0             // Maximum axis ratio for PSF model
+
 /*****************************************************************************/
 /* FUNCTION IMPLEMENTATION - PUBLIC                                          */
@@ -405,2 +408,23 @@
     return psf;
 }
+
+
+float pmPSFtoFWHM(const pmPSF *psf, float x, float y)
+{
+    PS_ASSERT_PTR_NON_NULL(psf, NAN);
+
+    pmModel *model = pmModelFromPSFforXY(psf, x, y, 1.0); // Model of source
+    if (!model) {
+        psError(PS_ERR_UNKNOWN, false, "Unable to determine PSF model at %f,%f\n", x, y);
+        return NAN;
+    }
+    psF32 *params = model->params->data.F32; // Model parameters
+    psEllipseAxes axes = pmPSF_ModelToAxes(params, MAX_AXIS_RATIO); // Ellipse axes
+
+    // Curiously, the minor axis can be larger than the major axis, so need to check.
+    float fwhm = 2.355 * PS_MAX(axes.minor, axes.major); // FWHM, converted from sigma
+
+    psFree(model);
+
+    return fwhm;
+}
Index: trunk/psModules/src/objects/pmPSF.h
===================================================================
--- trunk/psModules/src/objects/pmPSF.h	(revision 26450)
+++ trunk/psModules/src/objects/pmPSF.h	(revision 26893)
@@ -111,4 +111,11 @@
 psEllipseAxes pmPSF_ModelToAxes (psF32 *modelPar, double maxAR);
 
+/// Calculate FWHM value from a PSF
+float pmPSFtoFWHM(
+    const pmPSF *psf,                   // PSF of interest
+    float x, float y                    // Position of interest
+    );
+
+
 /// @}
 # endif
Index: trunk/psModules/src/objects/pmPeaks.c
===================================================================
--- trunk/psModules/src/objects/pmPeaks.c	(revision 26450)
+++ trunk/psModules/src/objects/pmPeaks.c	(revision 26893)
@@ -416,5 +416,5 @@
 
         } else {
-            psError(PS_ERR_UNKNOWN, true, "peak specified valid column range.");
+            psLogMsg ("psModules.objects", 5, "peak specified outside valid column range.");
         }
     }
@@ -501,5 +501,5 @@
                 }
             } else {
-                psError(PS_ERR_UNKNOWN, true, "peak specified outside valid column range.");
+		psLogMsg ("psModules.objects", 5, "peak specified outside valid column range.");
             }
 
@@ -545,5 +545,5 @@
             }
         } else {
-            psError(PS_ERR_UNKNOWN, true, "peak specified outside valid column range.");
+            psLogMsg ("psModules.objects", 5, "peak specified outside valid column range.");
         }
     }
Index: trunk/psModules/src/objects/pmPhotObj.c
===================================================================
--- trunk/psModules/src/objects/pmPhotObj.c	(revision 26893)
+++ trunk/psModules/src/objects/pmPhotObj.c	(revision 26893)
@@ -0,0 +1,39 @@
+/** @file  pmPhotObj.c
+ *
+ *  @author EAM, IfA
+ *
+ *  @version $Revision: 1.13 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2006-09-15 09:49:01 $
+ * Copyright 2009 Institute for Astronomy, University of Hawaii
+ */
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include <stdio.h>
+#include <math.h>
+#include <string.h>
+#include <pslib.h>
+#include "pmPhotObj.h"
+
+static void pmPhotObjFree (pmPhotObj *tmp)
+{
+    if (!tmp) return;
+    psFree(tmp->sources);
+}
+
+
+pmPhotObj *pmPhotObjAlloc() 
+{
+    static int id = 1;
+    pmPhotObj *obj = (pmPhotObj *) psAlloc(sizeof(pmPhotObj));
+    psMemSetDeallocator(obj, (psFreeFunc) pmPhotObjFree);
+
+    *(int *)&obj->id = id; // cast away the const to set the id.
+    id++;
+
+    obj->sources = NULL;
+    return obj;
+}
+
Index: trunk/psModules/src/objects/pmPhotObj.h
===================================================================
--- trunk/psModules/src/objects/pmPhotObj.h	(revision 26893)
+++ trunk/psModules/src/objects/pmPhotObj.h	(revision 26893)
@@ -0,0 +1,41 @@
+/* @file  pmPhotObj.h
+ *
+ * @author EAM, IfA
+ *
+ * @version $Revision: $
+ * @date $Date: 2009-02-16 22:30:50 $
+ * Copyright 2009 Institute for Astronomy, University of Hawaii
+ */
+
+# ifndef PM_PHOT_OBJ_H
+# define PM_PHOT_OBJ_H
+
+#include <pslib.h>
+#include "pmPeaks.h"
+#include "pmModel.h"
+#include "pmMoments.h"
+
+/// @addtogroup Objects Object Detection / Analysis Functions
+/// @{
+
+
+/** pmPhotObj data structure
+ *
+ *  A Photometry Object is a connected set of source measurements with a common
+ *  connection.  Each source that comprises the object represents the detection of some
+ *  astronomical object on a single image.  The object represents the related collection
+ *  of measurements.  The fits are coupled in some way.  For example, they may all have
+ *  the same position, but independent fluxes.  Or, they may have a common set of
+ *  positions and shape parameters.  Or the position in each image may be related by a
+ *  function.
+ *
+ *  XXX is thre any info that is neaded for each object beyond that carried by the sources
+ *  (besides ID)?
+ */
+typedef struct {
+  int seq;                            ///< ID for output (generated on write OR set on read)
+  psArray *sources;
+} pmPhotObj;
+
+/// @}
+# endif /* PM_PHOT_OBJ_H */
Index: trunk/psModules/src/objects/pmSource.c
===================================================================
--- trunk/psModules/src/objects/pmSource.c	(revision 26450)
+++ trunk/psModules/src/objects/pmSource.c	(revision 26893)
@@ -3,12 +3,10 @@
  *  Functions to define and manipulate sources on images
  *
- *  @author GLG, MHPCC
- *  @author EAM, IfA: significant modifications.
+ *  @author EAM, IfA
+ *  @author GLG, MHPCC (initial code base)
  *
  *  @version $Revision: 1.70 $ $Name: not supported by cvs2svn $
  *  @date $Date: 2009-02-16 22:29:59 $
- *
- *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
- *
+ *  Copyright 2009 Institute for Astronomy, University of Hawaii
  */
 
@@ -277,5 +275,5 @@
 // psphot-specific function which applies the recipe values
 // only apply selection to sources within specified region
-pmPSFClump pmSourcePSFClump(psRegion *region, psArray *sources, psMetadata *recipe)
+pmPSFClump pmSourcePSFClump(psImage **savedImage, psRegion *region, psArray *sources, float PSF_SN_LIM, float PSF_CLUMP_GRID_SCALE, psF32 SX_MAX, psF32 SY_MAX, psF32 AR_MAX)
 {
     psTrace("psModules.objects", 10, "---- begin ----\n");
@@ -287,33 +285,7 @@
 
     PS_ASSERT_PTR_NON_NULL(sources, errorClump);
-    PS_ASSERT_PTR_NON_NULL(recipe, errorClump);
-
-    bool status = true;                 // Status of MD lookup
-    float PSF_SN_LIM = psMetadataLookupF32(&status, recipe, "PSF_SN_LIM");
-    if (!status) {
-        PSF_SN_LIM = 0;
-    }
-    float PSF_CLUMP_GRID_SCALE = psMetadataLookupF32(&status, recipe, "PSF_CLUMP_GRID_SCALE");
-    if (!status) {
-        PSF_CLUMP_GRID_SCALE = 0.1;
-    }
 
     // find the sigmaX, sigmaY clump
     {
-        psF32 SX_MAX = psMetadataLookupF32(&status, recipe, "MOMENTS_SX_MAX");
-        if (!status) {
-            psWarning("MOMENTS_SX_MAX not set in recipe");
-            SX_MAX = 10.0;
-        }
-        psF32 SY_MAX = psMetadataLookupF32(&status, recipe, "MOMENTS_SY_MAX");
-        if (!status) {
-            psWarning("MOMENTS_SY_MAX not set in recipe");
-            SY_MAX = 10.0;
-        }
-        psF32 AR_MAX = psMetadataLookupF32(&status, recipe, "MOMENTS_AR_MAX");
-        if (!status) {
-            psWarning("MOMENTS_AR_MAX not set in recipe");
-            AR_MAX =  3.0;
-        }
         psF32 AR_MIN = 1.0 / AR_MAX;
 
@@ -401,9 +373,6 @@
 	psfClump.nSigma = stats->sampleStdev;
 
-        const bool keep_psf_clump = psMetadataLookupBool(NULL, recipe, "KEEP_PSF_CLUMP");
-        if (keep_psf_clump)
-        {
-            psMetadataAdd(recipe, PS_LIST_TAIL,
-                          "PSF_CLUMP", PS_DATA_IMAGE, "Image of PSF coefficients", splane);
+	if (savedImage) {
+	    *savedImage = psMemIncrRefCounter(splane);
         }
         psFree (splane);
@@ -411,5 +380,12 @@
 
         // if we failed to find a valid peak, return the empty clump (failure signal)
-        if (!peaks || !peaks->n) {
+	if (peaks == NULL) {
+            psError(PS_ERR_UNKNOWN, false, "failure in peak analysis for PSF clump.\n");
+	    psFree (peaks);
+            return emptyClump;
+	}
+
+        if (peaks->n == 0)
+        {
             psLogMsg ("psphot", 3, "failed to find a peak in the PSF clump image\n");
             if (nValid == 0) {
@@ -418,4 +394,5 @@
                 psLogMsg ("psphot", 3, "no significant peak\n");
             }
+	    psFree (peaks);
             return (emptyClump);
         }
@@ -524,10 +501,9 @@
 *****************************************************************************/
 
-bool pmSourceRoughClass(psRegion *region, psArray *sources, psMetadata *recipe, pmPSFClump clump, psImageMaskType maskSat)
+bool pmSourceRoughClass(psRegion *region, psArray *sources, float PSF_SN_LIM, float PSF_CLUMP_NSIGMA, pmPSFClump clump, psImageMaskType maskSat)
 {
     psTrace("psModules.objects", 10, "---- begin ----");
 
     PS_ASSERT_PTR_NON_NULL(sources, false);
-    PS_ASSERT_PTR_NON_NULL(recipe, false);
 
     int Nsat     = 0;
@@ -542,13 +518,4 @@
     psVector *starsn_peaks = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
     psVector *starsn_moments = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
-
-    // get basic parameters, or set defaults
-    bool status;
-    float PSF_SN_LIM = psMetadataLookupF32 (&status, recipe, "PSF_SN_LIM");
-    if (!status) PSF_SN_LIM = 20.0;
-    float PSF_CLUMP_NSIGMA = psMetadataLookupF32 (&status, recipe, "PSF_CLUMP_NSIGMA");
-    if (!status) PSF_CLUMP_NSIGMA = 1.5;
-
-    // float INNER_RADIUS = psMetadataLookupF32 (&status, recipe, "SKY_INNER_RADIUS");
 
     pmSourceMode noMoments = PM_SOURCE_MODE_MOMENTS_FAILURE | PM_SOURCE_MODE_SKYVAR_FAILURE | PM_SOURCE_MODE_SKY_FAILURE | PM_SOURCE_MODE_BELOW_MOMENTS_SN;
Index: trunk/psModules/src/objects/pmSource.h
===================================================================
--- trunk/psModules/src/objects/pmSource.h	(revision 26450)
+++ trunk/psModules/src/objects/pmSource.h	(revision 26893)
@@ -176,13 +176,15 @@
  *
  * The return value indicates the success (TRUE) of the operation.
- *
- * XXX: Limit the S/N of the candidate sources (part of Metadata)? (TBD).
- * XXX: Save the clump parameters on the Metadata (TBD)
- *
- */
+ */
+
 pmPSFClump pmSourcePSFClump(
+    psImage **savedImage, 
     psRegion *region,                   ///< restrict measurement to specified region
     psArray *source,                    ///< The input pmSource
-    psMetadata *metadata                ///< Contains classification parameters
+    float PSF_SN_LIM, 
+    float PSF_CLUMP_GRID_SCALE, 
+    psF32 SX_MAX, 
+    psF32 SY_MAX, 
+    psF32 AR_MAX
 );
 
@@ -200,5 +202,6 @@
     psRegion *region,                   ///< restrict measurement to specified region
     psArray *sources,                    ///< The input pmSources
-    psMetadata *metadata,               ///< Contains classification parameters
+    float PSF_SN_LIM,			 ///< min S/N for source to be used for PSF model
+    float PSF_CLUMP_NSIGMA,		 ///< size of region around peak of clump for PSF stars
     pmPSFClump clump,                   ///< Statistics about the PSF clump
     psImageMaskType maskSat             ///< Mask value for saturated pixels
@@ -220,5 +223,6 @@
     float radius,     ///< Use a circle of pixels around the peak
     float sigma,      ///< size of Gaussian window function (<= 0.0 -> skip window)
-    float minSN	      ///< minimum pixel significance
+    float minSN,	      ///< minimum pixel significance
+    psImageMaskType maskVal
 );
 
Index: trunk/psModules/src/objects/pmSourceIO.c
===================================================================
--- trunk/psModules/src/objects/pmSourceIO.c	(revision 26450)
+++ trunk/psModules/src/objects/pmSourceIO.c	(revision 26893)
@@ -40,4 +40,5 @@
 #include "pmPSF.h"
 #include "pmModel.h"
+#include "pmDetections.h"
 #include "pmSource.h"
 #include "pmModelClass.h"
@@ -344,10 +345,17 @@
 
     // if sources is NULL, write out an empty table
-    // input / output sources are stored on the readout->analysis as "PSPHOT.SOURCES" -- a better name might be something like PM_SOURCE_DATA
-    psArray *sources = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.SOURCES");
+    // input / output sources are stored on the readout->analysis as "PSPHOT.DETECTIONS"
+
+    psArray *sources = NULL;
+    pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS");
+    if (detections) {
+	sources = detections->allSources;
+    }
     if (!sources) {
+	detections = pmDetectionsAlloc();
         sources = psArrayAlloc(0);
-        psMetadataAddArray(readout->analysis, PS_LIST_TAIL, "PSPHOT.SOURCES", PS_META_REPLACE, "Blank array of sources", sources);
-        psFree(sources); // Held onto by the metadata, so we can continue to use
+	detections->allSources = sources;
+        psMetadataAddPtr(readout->analysis, PS_LIST_TAIL, "PSPHOT.DETECTIONS", PS_DATA_UNKNOWN | PS_META_REPLACE, "Blank array of sources", detections);
+        psFree(detections); // Held onto by the metadata, so we can continue to use
     }
 
@@ -1031,6 +1039,9 @@
     }
     readout->data_exists = true;
-    status = psMetadataAdd (readout->analysis, PS_LIST_TAIL, "PSPHOT.SOURCES", PS_DATA_ARRAY, "input sources", sources);
-    psFree (sources);
+
+    pmDetections *detections = pmDetectionsAlloc();
+    detections->allSources = sources;
+    status = psMetadataAdd (readout->analysis, PS_LIST_TAIL, "PSPHOT.DETECTIONS", PS_DATA_ARRAY, "input sources", detections);
+    psFree (detections);
     return true;
 }
@@ -1124,9 +1135,10 @@
     bool status;
 
-    // select the psf of interest
-    pmPSF *psf = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.SOURCES");
-    if (!psf) return false;
-    return true;
-}
-
-
+    // select the detections of interest
+    pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS");
+    if (!detections) return false;
+    if (!detections->allSources) return false;
+    return true;
+}
+
+
Index: trunk/psModules/src/objects/pmSourceIO_CMF_PS1_V2.c
===================================================================
--- trunk/psModules/src/objects/pmSourceIO_CMF_PS1_V2.c	(revision 26450)
+++ trunk/psModules/src/objects/pmSourceIO_CMF_PS1_V2.c	(revision 26893)
@@ -304,4 +304,5 @@
         source->peak->dx   = dPAR[PM_PAR_XPOS];
         source->peak->dy   = dPAR[PM_PAR_YPOS];
+        source->peak->SN   = sqrt(source->peak->flux); // XXX a proxy: various functions sort by peak S/N
 
         source->pixWeight = psMetadataLookupF32 (&status, row, "PSF_QF");
Index: trunk/psModules/src/objects/pmSourceIO_MatchedRefs.c
===================================================================
--- trunk/psModules/src/objects/pmSourceIO_MatchedRefs.c	(revision 26450)
+++ trunk/psModules/src/objects/pmSourceIO_MatchedRefs.c	(revision 26893)
@@ -163,4 +163,5 @@
 
     // try find the MATCHED_REFS extension.  if non-existent, note that we tried, and move on.
+    // It is not an error to lack this entry -- psFitsMoveExtNameClean does not raise an error
     if (!psFitsMoveExtNameClean (fits, "MATCHED_REFS")) {
 	psMetadataAddBool (fpa->analysis, PS_LIST_TAIL, "READ.REFMATCH", PS_META_REPLACE, "attempted to read MATCHED_REFS", true);
Index: trunk/psModules/src/objects/pmSourceMoments.c
===================================================================
--- trunk/psModules/src/objects/pmSourceMoments.c	(revision 26450)
+++ trunk/psModules/src/objects/pmSourceMoments.c	(revision 26893)
@@ -54,5 +54,5 @@
 # define VALID_RADIUS(X,Y,RAD2) (((RAD2) >= (PS_SQR(X) + PS_SQR(Y))) ? 1 : 0)
 
-bool pmSourceMoments(pmSource *source, psF32 radius, psF32 sigma, psF32 minSN)
+bool pmSourceMoments(pmSource *source, psF32 radius, psF32 sigma, psF32 minSN, psImageMaskType maskVal)
 {
     PS_ASSERT_PTR_NON_NULL(source, false);
@@ -114,5 +114,5 @@
         for (psS32 col = 0; col < source->pixels->numCols ; col++, vPix++, vWgt++) {
             if (vMsk) {
-                if (*vMsk) {
+                if (*vMsk & maskVal) {
                     vMsk++;
                     continue;
@@ -135,5 +135,5 @@
 	    // stars.
             if (PS_SQR(pDiff) < minSN2*wDiff) continue;
-            if (pDiff < 0) continue; // XXX : MWV says I should include < 0.0 valued points...
+            // if (pDiff < 0) continue; // XXX : MWV says I should include < 0.0 valued points...
 
 	    // Apply a Gaussian window function.  Be careful with the window function.  S/N
@@ -226,5 +226,5 @@
 	for (psS32 col = 0; col < source->pixels->numCols ; col++, vPix++, vWgt++) {
 	    if (vMsk) {
-		if (*vMsk) {
+		if (*vMsk & maskVal) {
 		    vMsk++;
 		    continue;
@@ -249,5 +249,5 @@
 	    // stars.
 	    if (PS_SQR(pDiff) < minSN2*wDiff) continue;
-	    if (pDiff < 0) continue;
+	    // if (pDiff < 0) continue;
 
 	    // Apply a Gaussian window function.  Be careful with the window function.  S/N
@@ -315,10 +315,10 @@
     source->moments->Myyyy = YYYY/Sum;
 
-    if (source->moments->Mxx < 0) {
-	fprintf (stderr, "error: neg second moment??\n");
-    }
-    if (source->moments->Myy < 0) {
-	fprintf (stderr, "error: neg second moment??\n");
-    }
+    // if (source->moments->Mxx < 0) {
+    // 	fprintf (stderr, "error: neg second moment??\n");
+    // }
+    // if (source->moments->Myy < 0) {
+    // 	fprintf (stderr, "error: neg second moment??\n");
+    // }
 
     psTrace ("psModules.objects", 4, "Mxx: %f  Mxy: %f  Myy: %f  Mxxx: %f  Mxxy: %f  Mxyy: %f  Myyy: %f  Mxxxx: %f  Mxxxy: %f  Mxxyy: %f  Mxyyy: %f  Mxyyy: %f\n",
Index: trunk/psModules/src/objects/pmSourcePlotApResid.c
===================================================================
--- trunk/psModules/src/objects/pmSourcePlotApResid.c	(revision 26450)
+++ trunk/psModules/src/objects/pmSourcePlotApResid.c	(revision 26893)
@@ -34,4 +34,5 @@
 #include "pmPSF.h"
 #include "pmModel.h"
+#include "pmDetections.h"
 #include "pmSource.h"
 #include "pmSourcePlots.h"
@@ -53,4 +54,5 @@
     PS_ASSERT_PTR_NON_NULL(layout, false);
 
+    bool status;
     Graphdata graphdata;
     KapaSection section;
@@ -61,7 +63,9 @@
     pmReadout  *readout = pmFPAfileThisReadout (config->files, view, "PSPHOT.INPUT");
 
-    psArray *sources = psMetadataLookupPtr (NULL, readout->analysis, "PSPHOT.SOURCES");
-    if (sources == NULL)
-        return false;
+    pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS");
+    if (detections == NULL) return false;
+
+    psArray *sources = detections->allSources;
+    if (sources == NULL) return false;
 
     int kapa = pmKapaOpen (false);
Index: trunk/psModules/src/objects/pmSourcePlotMoments.c
===================================================================
--- trunk/psModules/src/objects/pmSourcePlotMoments.c	(revision 26450)
+++ trunk/psModules/src/objects/pmSourcePlotMoments.c	(revision 26893)
@@ -37,4 +37,5 @@
 #include "pmPSF.h"
 #include "pmModel.h"
+#include "pmDetections.h"
 #include "pmSource.h"
 #include "pmSourcePlots.h"
@@ -54,4 +55,5 @@
     PS_ASSERT_PTR_NON_NULL(layout, false);
 
+    bool status;
     Graphdata graphdata;
     KapaSection section;
@@ -62,7 +64,9 @@
     pmReadout  *readout = pmFPAfileThisReadout (config->files, view, "PSPHOT.INPUT");
 
-    psArray *sources = psMetadataLookupPtr (NULL, readout->analysis, "PSPHOT.SOURCES");
-    if (sources == NULL)
-        return false;
+    pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS");
+    if (detections == NULL) return false;
+
+    psArray *sources = detections->allSources;
+    if (sources == NULL) return false;
 
     int kapa = pmKapaOpen (false);
Index: trunk/psModules/src/objects/pmSourcePlotPSFModel.c
===================================================================
--- trunk/psModules/src/objects/pmSourcePlotPSFModel.c	(revision 26450)
+++ trunk/psModules/src/objects/pmSourcePlotPSFModel.c	(revision 26893)
@@ -37,4 +37,5 @@
 #include "pmPSF.h"
 #include "pmModel.h"
+#include "pmDetections.h"
 #include "pmSource.h"
 #include "pmSourcePlots.h"
@@ -56,4 +57,5 @@
     PS_ASSERT_PTR_NON_NULL(layout, false);
 
+    bool status;
     Graphdata graphdata;
     KapaSection section;
@@ -64,7 +66,9 @@
     pmReadout  *readout = pmFPAfileThisReadout (config->files, view, "PSPHOT.INPUT");
 
-    psArray *sources = psMetadataLookupPtr (NULL, readout->analysis, "PSPHOT.SOURCES");
-    if (sources == NULL)
-        return false;
+    pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS");
+    if (detections == NULL) return false;
+
+    psArray *sources = detections->allSources;
+    if (sources == NULL) return false;
 
     int kapa = pmKapaOpen (false);