Index: trunk/psphot/src/psphotBasicDeblend.c
===================================================================
--- trunk/psphot/src/psphotBasicDeblend.c	(revision 5704)
+++ trunk/psphot/src/psphotBasicDeblend.c	(revision 5718)
@@ -1,74 +1,134 @@
 # include "psphot.h"
 
-// try PSF models and select best option
+bool psphotBasicDeblend (psArray *sources, psMetadata *config, psStats *sky) { 
 
-pmPSF *psphotBasicDeblend (psMetadata *config, psArray *sources, psStats *skystats) 
-{ 
+    int N;
+    bool status;
+    float threshold;
+    pmSource *source, *testSource;
 
-    // source analysis is done in S/N order (brightest first)
+    FILE *f = fopen ("deblend.dat", "w");
+    if (f == NULL) psAbort ("psphot", "can't open deblend.dat output file");
+    int Nblend = 0;
+
+    psTimerStart ("psphot");
+
+    // this should be added to the PM_SOURCE flags:
+    int PM_SOURCE_BLEND = PM_SOURCE_OTHER + 1;
+
+    float FRACTION = psMetadataLookupF32 (&status, config, "DEBLEND_PEAK_FRACTION");
+    float NSIGMA   = psMetadataLookupF32 (&status, config, "DEBLEND_SKY_NSIGMA");
+    float minThreshold = NSIGMA*sky->sampleStdev;
+    fprintf (stderr, "min threshold: %f\n", minThreshold);
+
+    // we need sources spatially-sorted to find overlaps
     sources = psArraySort (sources, psphotSortByY);
 
+    // source analysis is done in peak order (brightest first)
+    // we use an index for this so the spatial sorting is kept
     psVector *SN = psVectorAlloc (sources->n, PS_DATA_F32);
     for (int i = 0; i < SN->n; i++) {
-      source = sources->data[i];
-      SN->data.F32[i] = source->moments->SN;
+	source = sources->data[i];
+	SN->data.F32[i] = source->peak->counts;
     }
     psVector *index = psVectorSortIndex (NULL, SN);
-    
+    // this results in an index of increasing SN
+
+    // temporary array for overlapping objects we find
     psArray *overlap = psArrayAlloc (100);
 
-    for (int i = 0; i < sources->n; i++) {
-      N = index->data.U32[i];
-      source = sources->data[N];
-      overlap->n = 0;
+    // XXX make sure this results in decreasing, not increasing, SN
+    for (int i = sources->n - 1; i >= 0; i--) {
+	N = index->data.U32[i];
+	source = sources->data[N];
 
-      // search backwards for overlapping sources
-      for (int j = N - 1; j >= 0; j--) {
-	testSource = sources->data[j];
-	if (testSource->moments->x < source->pixels->col0) continue;
-	if (testSource->moments->x >= source->pixels->col0 + source->pixels->numCols) continue;
-	if (testSource->moments->y < source->pixels->row0) break;
-	if (testSource->moments->y >= source->pixels->row0 + source->pixels->numRows) {
-	  fprintf (stderr, "warning: invalid condition\n");
-	  continue;
+	if (source->type == PM_SOURCE_BLEND) continue;
+
+	overlap->n = 0;
+
+	// search backwards for overlapping sources
+	for (int j = N - 1; j >= 0; j--) {
+	    testSource = sources->data[j];
+	    if (testSource->peak->x <  source->pixels->col0) continue;
+	    if (testSource->peak->x >= source->pixels->col0 + source->pixels->numCols) continue;
+	    if (testSource->peak->y <  source->pixels->row0) break;
+	    if (testSource->peak->y >= source->pixels->row0 + source->pixels->numRows) {
+		fprintf (stderr, "warning: invalid condition\n");
+		continue;
+	    }
+	    psArrayAdd (overlap, 100, testSource);
 	}
-	psArrayAdd (overlap, 100, testSource);
-      }
 
-      // search forwards for overlapping sources
-      for (int j = N + 1; j < sources->n; j++) {
-	testSource = sources->data[j];
-	if (testSource->moments->x < source->pixels->col0) continue;
-	if (testSource->moments->x >= source->pixels->col0 + source->pixels->numCols) continue;
-	if (testSource->moments->y < source->pixels->row0) {
-	  fprintf (stderr, "warning: invalid condition\n");
-	  continue;
+	// search forwards for overlapping sources
+	for (int j = N + 1; j < sources->n; j++) {
+	    testSource = sources->data[j];
+	    if (testSource->peak->x <  source->pixels->col0) continue;
+	    if (testSource->peak->x >= source->pixels->col0 + source->pixels->numCols) continue;
+	    if (testSource->peak->y <  source->pixels->row0) {
+		fprintf (stderr, "warning: invalid condition\n");
+		continue;
+	    }
+	    if (testSource->peak->y >= source->pixels->row0 + source->pixels->numRows) break;
+	    psArrayAdd (overlap, 100, testSource);
 	}
-	if (testSource->moments->y >= source->pixels->row0 + source->pixels->numRows) break;
-	psArrayAdd (overlap, 100, testSource);
-      }
 
-      // generate source contour (1/4 peak counts)
-      if (overlap->n > 0) {
-	// XXX EAM : make the 1/4 user-defined.
-	threshold = 0.25 * (source->moments->Peak - source->moments->Sky) + source->moments->Sky); 
-	psArray *contour = pmSourceContour_EAM (source->pixels, source->moments->x, source->moments->x, threshold);
+	// generate source contour (1/4 peak counts)
+	if (overlap->n > 0) {
+	    // XXX EAM : make the 1/4 user-defined.
+	    // XXX keep threshold from dropping too low (N*sky.sigma)
+	    threshold = FRACTION * (source->peak->counts - source->moments->Sky) + source->moments->Sky; 
+	    threshold = PS_MAX (threshold, minThreshold);
 
-	// XXX EAM : the order is wrong: need to check j, j+1 contour points
-	psVector *xv = contour->data[0];
-	psVector *yv = contour->data[1];
-	for (int j = 0; j < xv->n; j++) {
-	  for (int k = 0; k < overlap->n; k++) {
-	    testSource = overlap->data[k];
-	    if (fabs(yv->data.F32[j] - testSource->moments->y) > 0.5) continue;
-	    if (xv->data.F32[j] > testSource->moments->x) continue;
-	    if (xv->data.F32[j] < testSource->moments->x) continue;
-	    // XXX EAM : mark source with flag (blended)
-	  }
-	}  
+	    // XXX EAM : should the contour input coordinate be in parent or subimage coords? parent, for now
+	    psArray *contour = pmSourceContour_EAM (source->pixels, source->peak->x, source->peak->y, threshold);
+	    if (contour == NULL) {
+		fprintf (stderr, "below threshold? invalid peak?\n");
+		continue;
+	    }
+
+	    // the source contour consists of two vectors, xv and yv.  the contour is 
+	    // a series of coordinate pairs, (xv[i],yv[i]) & (xv[i+1],yv[i+1]).  both
+	    // coordinate pairs have the same yv[] value, with xv[i] corresponding to
+	    // the left boundary and xv[i+1] corresponding to the right boundary
+
+	    psVector *xv = contour->data[0];
+	    psVector *yv = contour->data[1];
+	    for (int k = 0; k < overlap->n; k++) {
+		testSource = overlap->data[k];
+		if (testSource->peak->counts > source->peak->counts) continue;
+		for (int j = 0; j < xv->n; j+=2) {
+		    if (fabs(yv->data.F32[j] - testSource->peak->y) > 0.5) continue;
+		    if (xv->data.F32[j+0] > testSource->peak->x) break;
+		    if (xv->data.F32[j+1] < testSource->peak->x) break;
+
+		    int xp0 = source->moments->x - source->pixels->col0;
+		    int xp1 = source->peak->x - source->pixels->col0;
+		    int xp2 = testSource->moments->x - source->pixels->col0;
+		    int xp3 = testSource->peak->x - source->pixels->col0;
+
+		    int yp0 = source->moments->y - source->pixels->row0;
+		    int yp1 = source->peak->y - source->pixels->row0;
+		    int yp2 = testSource->moments->y - testSource->pixels->row0;
+		    int yp3 = testSource->peak->y - testSource->pixels->row0;
+
+		    fprintf (f, "%d %d (%f, %f) :  %d %d (%f, %f)  vs %f\n", 
+			     source->peak->x, source->peak->y, 
+			     source->pixels->data.F32[yp0][xp0], source->pixels->data.F32[yp1][xp1], 
+			     testSource->peak->x, testSource->peak->y, 
+			     testSource->pixels->data.F32[yp2][xp2], testSource->pixels->data.F32[yp3][xp3], threshold
+			     );
+
+		    testSource->type = PM_SOURCE_BLEND;
+		    Nblend ++;
+		    j = xv->n;
+		}
+	    }  
+	}
     }
+    fclose (f);
+    psLogMsg ("psphot.deblend", 3, "identified %d blended objects (%f)\n", Nblend, psTimerMark ("psphot"));
+    return true;
+}
 
-    psLogMsg ("psphot.pspsf", 3, "selected psf model %s, ApResid: %f +/- %f\n", modelName, psf->ApResid, psf->dApResid);
 
-    return (psf);
-}
+