Index: trunk/psphot/src/psphotEnsemblePSF.c
===================================================================
--- trunk/psphot/src/psphotEnsemblePSF.c	(revision 5704)
+++ trunk/psphot/src/psphotEnsemblePSF.c	(revision 5718)
@@ -13,4 +13,7 @@
     sources = psArraySort (sources, psphotSortByY);
 
+    // this should be added to the PM_SOURCE flags:
+    int PM_SOURCE_BLEND = PM_SOURCE_OTHER + 1;
+
     float OUTER_RADIUS     = psMetadataLookupF32 (&status, config, "SKY_OUTER_RADIUS");
     float PSF_FIT_NSIGMA   = psMetadataLookupF32 (&status, config, "PSF_FIT_NSIGMA");
@@ -23,6 +26,17 @@
     // pre-calculate all model pixels
     psArray  *models = psArrayAlloc (sources->n);
+    psVector *index  = psVectorAlloc (sources->n, PS_TYPE_U32);
+    models->n = 0;
+    index->n = 0;
+
     for (int i = 0; i < sources->n; i++) {
 	pmSource *inSource = sources->data[i];
+
+	// skip non-astronomical objects (very likely defects)
+	// XXX EAM : should we try these anyway?
+	if (inSource->type == PM_SOURCE_BLEND) continue;
+	if (inSource->type == PM_SOURCE_DEFECT) continue; 
+	if (inSource->type == PM_SOURCE_SATURATED) continue;
+
 	pmSource *otSource = pmSourceAlloc ();
 
@@ -72,5 +86,6 @@
 	// save source in list
 	otSource->modelPSF = model;
-	models->data[i] = otSource;
+	index->data.U32[models->n] = i;
+	psArrayAdd (models, 100, otSource);
     }
     psLogMsg ("psphot.emsemble", 4, "built models: %f (%d objects)\n", psTimerMark ("psphot"), sources->n);
@@ -79,7 +94,8 @@
 
     // fill out the sparse matrix
-    psSparse *sparse = psSparseAlloc (sources->n, 100);
-    for (int i = 0; i < sources->n; i++) {
-	pmSource *Fi = sources->data[i];
+    psSparse *sparse = psSparseAlloc (models->n, 100);
+    for (int i = 0; i < models->n; i++) {
+	int N = index->data.U32[i];
+	pmSource *Fi = sources->data[N];
 	pmSource *Mi = models->data[i];
 
@@ -93,13 +109,12 @@
 
 	// loop over all other stars following this one
-	for (int j = i + 1; j < sources->n; j++) {
+	for (int j = i + 1; j < models->n; j++) {
 	    pmSource *Mj = models->data[j];
 
-	    // XXX double check that this is working!  should not break here
-	    check me here;
-	    if (Mi->pixels->col0 + Mi->pixels->numCols < Mj->pixels->col0) break;
+	    // skip over disjoint source images, break after last possible overlap
+	    if (Mi->pixels->row0 + Mi->pixels->numRows < Mj->pixels->row0) break;
+	    if (Mj->pixels->row0 + Mj->pixels->numRows < Mi->pixels->row0) continue;
+	    if (Mi->pixels->col0 + Mi->pixels->numCols < Mj->pixels->col0) continue;
 	    if (Mj->pixels->col0 + Mj->pixels->numCols < Mi->pixels->col0) continue;
-	    if (Mi->pixels->row0 + Mi->pixels->numRows < Mj->pixels->row0) continue;
-	    if (Mj->pixels->row0 + Mj->pixels->numRows < Mi->pixels->row0) continue;
 	    
 	    // got an overlap; calculate cross-product and add to output array
@@ -115,6 +130,7 @@
 
     // adjust models, set sources and subtract
-    for (int i = 0; i < sources->n; i++) {
-	pmSource *Fi = sources->data[i];
+    for (int i = 0; i < models->n; i++) {
+	int N = index->data.U32[i];
+	pmSource *Fi = sources->data[N];
 	pmSource *Mi = models->data[i];