Index: trunk/psphot/src/psphotSourceStats.c
===================================================================
--- trunk/psphot/src/psphotSourceStats.c	(revision 41359)
+++ trunk/psphot/src/psphotSourceStats.c	(revision 42842)
@@ -120,4 +120,5 @@
         source->moments = pmMomentsAlloc();
 
+	// XXX can this ever be set at this point? (we just allocated the source and moments)
         if (source->mode & PM_SOURCE_MODE_MOMENTS_FAILURE) {
 	    fprintf (stderr, "moment failure\n");
@@ -144,4 +145,7 @@
 
     if (setWindow) {
+	// identify sources on lines, drop from MomentsWindow analysis
+	psphotFindFeatures (readout, sources);
+
         if (!psphotSetMomentsWindow(recipe, readout->analysis, sources, maskVal)) {
             psError(PS_ERR_UNEXPECTED_NULL, false, "Failed to determine Moments Window!");
@@ -503,4 +507,8 @@
     }
 
+    // use clipped stats or robust stats to define the clump
+    bool PSF_CLUMP_USE_CLIPPED_STATS = psMetadataLookupF32(&status, recipe, "PSF_CLUMP_USE_CLIPPED_STATS");
+    if (!status) { PSF_CLUMP_USE_CLIPPED_STATS = true; }
+
     // when we set the window, we are not attempting to measure spatial variations; we can use a somewhat higher S/N limit
     // since we are using all sources (true?)
@@ -539,11 +547,23 @@
     sources = psArraySort (sources, pmSourceSortByFlux);
 
+    // generate an array of the sources which we want to use for this analysis
+    // XXX move max source number to config
+    psArray *sourceSubset = psArrayAllocEmpty (100);
+    for (int i = 0; (i < sources->n) && (i < 400); i++) {
+
+	pmSource *source = sources->data[i];
+
+	// skip faint sources for moments measurement
+	if (sqrt(source->peak->detValue) < MIN_SN) continue;
+
+	// skip sources on lines
+	if (source->mode2 & PM_SOURCE_MODE2_ON_LINE) continue;
+	psArrayAdd (sourceSubset, 100, source);
+    }
+
     // loop over radii:
     for (int i = 0; i < nsigma; i++) {
-
-        // XXX move max source number to config
-        for (int j = 0; (j < sources->n) && (j < 400); j++) {
-
-            pmSource *source = sources->data[j];
+        for (int j = 0; j < sourceSubset->n; j++) {
+            pmSource *source = sourceSubset->data[j];
             psAssert (source->moments, "force moments to exist");
             source->moments->nPixels = 0;
@@ -565,6 +585,6 @@
 	  sprintf (name, "moments.v%d.dat", i);
 	  FILE *fout = fopen (name, "w");
-	  for (int j = 0; j < sources->n; j++) {
-            pmSource *source = sources->data[j];
+	  for (int j = 0; j < sourceSubset->n; j++) {
+            pmSource *source = sourceSubset->data[j];
             psAssert (source->moments, "force moments to exist");
             source->moments->nPixels = 0;
@@ -581,8 +601,8 @@
 
         // determine the PSF parameters from the source moment values
-        pmPSFClump psfClump = pmSourcePSFClump (NULL, NULL, sources, PSF_SN_LIM, PSF_CLUMP_GRID_SCALE, MOMENTS_SX_MAX, MOMENTS_SY_MAX, MOMENTS_SX_MIN, MOMENTS_SY_MIN, MOMENTS_AR_MAX);
+        pmPSFClump psfClump = pmSourcePSFClump (NULL, NULL, sourceSubset, PSF_SN_LIM, PSF_CLUMP_GRID_SCALE, MOMENTS_SX_MAX, MOMENTS_SY_MAX, MOMENTS_SX_MIN, MOMENTS_SY_MIN, MOMENTS_AR_MAX, PSF_CLUMP_USE_CLIPPED_STATS);
         psLogMsg ("psphot", 3, "sigma guess (pix) %.1f, nStars: %d of %d in clump, nSigma: %5.2f, X,  Y: %f, %f (%f, %f)\n", sigma[i], psfClump.nStars, psfClump.nTotal, psfClump.nSigma, psfClump.X, psfClump.Y, sqrt(psfClump.X) / sigma[i], sqrt(psfClump.Y) / sigma[i]);
 
-	Rmin[i] = pmSourceMinKronRadius(sources, PSF_SN_LIM);
+	Rmin[i] = pmSourceMinKronRadius(sourceSubset, PSF_SN_LIM);
 
 #if 0
@@ -600,5 +620,5 @@
         psMetadataAddF32 (regionMD, PS_LIST_TAIL, "PSF.CLUMP.DY", PS_META_REPLACE, "psf clump center", psfClump.dY);
 	if (pmVisualTestLevel("psphot.moments.full", 2)) {
-	    psphotVisualPlotMoments (recipe, analysis, sources);
+	    psphotVisualPlotMoments (recipe, analysis, sourceSubset);
 	}
 #endif
@@ -608,4 +628,6 @@
         Sout[i] = (Nout[i] == 0) ? NAN : sqrt(0.5*(psfClump.X + psfClump.Y)) / sigma[i];
     }
+
+    psFree (sourceSubset);
 
     // we are looking for sigma for which Sout = 0.65 (or some other value)