Index: trunk/ppStack/src/ppStackPrepare.c
===================================================================
--- trunk/ppStack/src/ppStackPrepare.c	(revision 34800)
+++ trunk/ppStack/src/ppStackPrepare.c	(revision 35394)
@@ -253,5 +253,6 @@
     float maxFWHM = psMetadataLookupF32(&mdok, recipe, "PSF.INPUT.MAX"); // max allowed input fwhm
     float clipFWHMnSig = psMetadataLookupF32(&mdok, recipe, "PSF.INPUT.CLIP.NSIGMA"); // sigma clipping of inputs
-
+    float threshFWHM = psMetadataLookupF32(&mdok, recipe, "PSF.INPUT.THRESH"); // FWHM size that we refuse to clip below
+    float asymmetryFWHM = psMetadataLookupF32(&mdok, recipe, "PSF.INPUT.ASYMMETRY"); // max bimodal asymmetry
     psString log = psStringCopy("Input seeing FWHMs:\n"); // Log message
     bool havePSFs = false;                                // Do we have any PSFs?
@@ -286,5 +287,5 @@
             options->inputSeeing->data.F32[i] = sumFWHM / (float)numFWHM;
         }
-
+#ifdef OLD_REJECTION
 	// reject any input images which exceed the specified max FWHM
 	if (isfinite(maxFWHM) && (options->inputSeeing->data.F32[i] > maxFWHM)) {
@@ -292,5 +293,5 @@
             psLogMsg("ppStack", PS_LOG_INFO, "PSF FWHM for image %d is too large (%f vs %f maxFWHM) --- rejected.", i, options->inputSeeing->data.F32[i], maxFWHM);
 	}
-
+#endif
         psStringAppend(&log, "Input %d: %f\n", i, options->inputSeeing->data.F32[i]);
     }
@@ -300,4 +301,69 @@
     psFree(log);
 
+    // Do GMM test with two modes to decide where to put the break.
+    double Punimodal;
+    int m = 2;
+    psVector *modes = psVectorAlloc(num,PS_TYPE_F32);
+    psVector *means = psVectorAlloc(m,PS_TYPE_F32);
+    psVector *S     = psVectorAlloc(m,PS_TYPE_F32);
+    psVector *pi    = psVectorAlloc(m,PS_TYPE_F32);
+    psImage *P      = psImageAlloc(m,num,PS_TYPE_F32);
+    float limit;
+    if (!psMM1DClass(options->inputSeeing,
+		     options->inputSeeing->n,
+		     modes,means,
+		     S,pi,P,
+		     2,
+		     &Punimodal)) {
+      // Handle error here
+      psFree(modes);
+      psFree(means);
+      psFree(S);
+      psFree(pi);
+      psFree(P);
+      return(false);
+    }
+    fprintf(stderr,"means: %g %g\n",means->data.F32[0],means->data.F32[1]);
+    fprintf(stderr,"sigma: %g %g \n",S->data.F32[0],S->data.F32[1]);
+
+    fprintf(stderr,"pi:    %g %g\n",pi->data.F32[0],pi->data.F32[1]);
+    
+    //    unwrittenGMMfunction(options->inputSeeing, options->inputMask, 0xff,
+    //		 &Punimodal,
+    //		 &m1,&s1,&pi1,
+    //		 &m2,&s2,&pi2);
+
+    // Use Gaussian mixture model analysis of the FWHM distribution to decide an optional FWHM limit
+    if ((Punimodal > 0.5)|| // This distribution is best represented by a single mode
+	(num <= 4)) {       // Or we have a small number of inputs, making this statistic poor.
+      limit = maxFWHM;
+    }
+    else {                  // This is a bimodal distribution
+      if ((fabs(pi->data.F32[0] - pi->data.F32[1]) < asymmetryFWHM)||  // However, both modes are equally populated
+	  (pi->data.F32[1] > pi->data.F32[0])) {                       // Or the larger FWHM mode is more populated
+	limit = means->data.F32[1] + clipFWHMnSig * S->data.F32[1];
+      }
+      else {                                   // The smaller FWHM mode is more populated
+	limit = means->data.F32[0] + clipFWHMnSig * S->data.F32[0];
+      }
+    }
+    if (limit > maxFWHM)    { limit = maxFWHM; }    // We should not be larger than our max
+    if (limit < threshFWHM) { limit = threshFWHM; } // Nor smaller than our min
+
+    for (int i = 0; i < num; i++) {
+      if (options->inputSeeing->data.F32[i] > limit) {
+	options->inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PPSTACK_MASK_PSF;
+	psLogMsg("ppStack", PS_LOG_INFO,
+		 "PSF FWHM for image %d is too large (%f vs %f fwhmlimit) --- rejected",
+		 i, options->inputSeeing->data.F32[i],limit);
+      }
+    }
+    psFree(means);
+    psFree(modes);
+    psFree(S);
+    psFree(pi);
+    psFree(P);
+    
+#ifdef OLD_REJECTION
     // We should have the ability to filter the input list based on the seeing:
     // * reject above some max value and/or min value
@@ -323,5 +389,5 @@
         }
     }
-
+#endif
     // Generate target PSF
     if (options->convolve) {