Index: /branches/czw_branch/20110406/psModules/src/imcombine/pmStack.c
===================================================================
--- /branches/czw_branch/20110406/psModules/src/imcombine/pmStack.c	(revision 31630)
+++ /branches/czw_branch/20110406/psModules/src/imcombine/pmStack.c	(revision 31631)
@@ -37,5 +37,5 @@
 
 
-# if (1)
+# if (0)
 #define TESTING                         // Enable test output
 /* #define TEST_X 5745                       // x coordinate to examine */
@@ -106,5 +106,4 @@
 
 // KMM functions to do bimodality rejection of pixels
-
 float gaussian(float x, float m, float s) {
   return(pow(s * sqrt(2 * M_PI),-1) * exp(-0.5 * pow( (x - m) / s, 2)));
@@ -116,4 +115,7 @@
 			 float *pi1, float *m1, float *s1,
 			 float *pi2, float *m2, float *s2) {
+  assert(values);
+  assert(values->type.type == PS_TYPE_F32);
+  
   double logL_bimodal = 0, logL_unimodal;
   psVector *P1 = psVectorAlloc(values->n,PS_TYPE_F32);
@@ -121,7 +123,4 @@
   int i;
 /*   int debug = 0; */
-/*   if (fabs(values->data.F32[0] - -145.624954) < 1e-4) { */
-/*     debug = 1; */
-/*   } */
   
   // Calculate unimodal properties
@@ -156,4 +155,5 @@
   float w1,w2;
 
+  // These should be options.
   float KMM_TOLERANCE = 1e-6;
   int KMM_MAX_ITERATIONS = 500;
@@ -191,5 +191,5 @@
     *s1 = 0;
     *s2 = 0;
-    for (i = 0; i < values->n; i++) { // Calculate new means
+    for (i = 0; i < values->n; i++) { // Calculate new means and weights
       *m1 += values->data.F32[i] * P1->data.F32[i];
       *m2 += values->data.F32[i] * P2->data.F32[i];
@@ -226,9 +226,9 @@
   // Calculate Punimodal
   double lambda = -2.0 * (logL_unimodal - logL_bimodal);
-  int    df     = 2 + 2 * 1;
+  int    df     = 2 + 2 * 1; // I can't find my reference on this. 
   if (lambda > 0) {
     *Punimodal = gsl_cdf_chisq_Q(lambda,df);
   }
-  else {
+  else { // If lambda <= 0, then logL_unimodal > logL_bimodal, so Punimodal must be by definition 1.0
     *Punimodal = 1.0;
   }
@@ -238,9 +238,12 @@
 
 static void KMMFindPopular(const psVector *values, float *Punimodal, float *mean, float *sigma, float *pi) {
-  float KMM_MINIMUM_PVALUE = 0.05;
+  float KMM_MINIMUM_PVALUE = 0.05; // Should be an option.
   float mU,sU;
   float pi1,m1,s1,pi2,m2,s2;
   int iter;
 
+  assert(values);
+  assert(values->type.type == PS_TYPE_F32);
+  
   KMMcalculate(values,Punimodal,&iter,
 	       &mU,&sU,
@@ -255,4 +258,5 @@
   }
   else {
+    // Is bimodal. Select most popular mode.
     if (pi1 >= pi2) {
       *mean = m1;
@@ -455,4 +459,7 @@
 }
 #if (0)
+// Currently unused function to reject inputs for a given pixel(x,y) based on the selecting the
+// most popular mode after running the KMM test, and rejecting all inputs that belong to the
+// least popular mode.
 static void KMMRejectUnpopular(const psArray *inputs, int x, int y) {
   float KMM_MINIMUM_PVALUE = 0.05;
@@ -479,6 +486,4 @@
 	       &pi1,&m1,&s1,
 	       &pi2,&m2,&s2);
-  //  fprintf(stderr,
-/*   psTrace("psModules.imcombine",3, */
 
    CHECKPIX(x, y, 
@@ -516,4 +521,6 @@
 }
 
+// Currently unused function to reject inputs for a given pixel(x,y) based on the selecting the
+// faintest mode as determined by the KMM test, and rejecting all inputs that belong to the brighest.
 static void KMMRejectBright(const psArray *inputs, int x, int y) {
   float KMM_MINIMUM_PVALUE = 0.05;
@@ -558,5 +565,5 @@
   // else do nothing.
 }
-#endif
+#endif // End if(0) to prevent KMMReject{Unpopular|Bright} from being defined.
 
 // Extract vectors for simple combination/rejection operations
@@ -701,5 +708,5 @@
     if (PS_SQR(x - TEST_X) + PS_SQR(y - TEST_Y) <= PS_SQR(TEST_RADIUS)) {
         for (int i = 0; i < numGood; i++) {
-	  fprintf(stderr,"Input %d, pixel %d,%d (%" PRIu16 "): %f %f (%f) %f %f %d %x %x -> %x %x\n",
+	    fprintf(stderr,"Input %d, pixel %d,%d (%" PRIu16 "): %f %f (%f) %f %f %d %x %x -> %x %x\n",
                     i, x, y, pixelSources->data.U16[i], pixelData->data.F32[i], pixelVariances->data.F32[i],
                     addVariance->data.F32[i], pixelWeights->data.F32[i], pixelExps->data.F32[i],
@@ -843,5 +850,10 @@
 
     // KMM values;
-    float Punimodal,KMMmean,KMMsigma,KMMpi;
+    float Punimodal = 1.0,KMMmean = NAN,KMMsigma = NAN,KMMpi = NAN;
+    int KMM_MINIMUM_INPUTS = 6;
+    bool useKMM = false;
+    if (num >= KMM_MINIMUM_INPUTS) {
+      useKMM = true;
+    }
     
     // Set up rejection limits
@@ -853,15 +865,23 @@
 
 	// Determine the systematic error from the most popular population in the sample
-	KMMFindPopular(pixelData,&Punimodal,&KMMmean,&KMMsigma,&KMMpi);
-	CHECKPIX(x,y,"KMM Popularity Contest: (%d,%d) Puni: %g Mean: %f Sigma %f Pi: %f\n",
-		 x,y,Punimodal,KMMmean,KMMsigma,KMMpi);
-	
+	// This should probably be an option
+	if (useKMM) {
+	  KMMFindPopular(pixelData,&Punimodal,&KMMmean,&KMMsigma,&KMMpi);
+	  CHECKPIX(x,y,"KMM Popularity Contest: (%d,%d) Puni: %g Mean: %f Sigma %f Pi: %f\n",
+		   x,y,Punimodal,KMMmean,KMMsigma,KMMpi);
+	}
         for (int i = 0; i < num; i++) {
             sumWeights += pixelWeights->data.F32[i];
         }
         for (int i = 0; i < num; i++) {
-            // Systematic error contributes to the rejection level
-/*             float sysVar = PS_SQR(sys * pixelData->data.F32[i]); */
-	  float sysVar = PS_SQR(KMMsigma);
+	  // Systematic error contributes to the rejection level
+	  float sysVar;
+	  if (useKMM) { // If we can trust KMM results, set the systematic variance
+	    sysVar = PS_SQR(KMMsigma);
+	  }
+	  else { // Otherwise, use the 10% systematic variance we've done in the past.
+	    sysVar = PS_SQR(sys * pixelData->data.F32[i]);
+	  }
+
 	    CHECKPIX(x, y, "Variance %d (%d), pixel %d,%d: %f %f %f\n", 
 		     i, pixelSources->data.U16[i], x, y, 
@@ -1013,16 +1033,16 @@
               if (useVariance) {
 		float median;
-		if (Punimodal > 0.05) {
+		if ((useKMM)&&(Punimodal < 0.05)) {
+		  median = KMMmean;
+		}
+		else {
                   median = combinationWeightedOlympic(pixelData, pixelWeights,
 						      olympic, buffer->sort); // Median for stack
 		}
-		else {
-		  median = KMMmean;
-		}
+		
 		  CHECKPIX(x, y, "Flag with variance pixel %d,%d: median = %f\n", x, y, median);
                   float worst = -INFINITY; // Largest deviation
                   for (int j = 0; j < num; j++) {
-		    
-                      float diff = pixelData->data.F32[j] - median; // Difference from expected
+		      float diff = pixelData->data.F32[j] - median; // Difference from expected
 		      CHECKPIX(x, y, "Testing input %d for pixel %d,%d: %f\n", j, x, y, diff);
 
@@ -1273,6 +1293,4 @@
 }
 
-  
-  
 
 //////////////////////////////////////////////////////////////////////////////////////////////////////////////
@@ -1445,6 +1463,5 @@
     for (int y = minInputRows; y < maxInputRows; y++) {
         for (int x = minInputCols; x < maxInputCols; x++) {
-
-/* 	    CHECKPIX(x, y, "Combining pixel %d,%d: %x %x %f %f %f %f %d %d %d\n", x, y, badMaskBits, blankMaskBits, iter, rej, sys, olympic, useVariance, safe, rejection); */
+	    CHECKPIX(x, y, "Combining pixel %d,%d: %x %x %f %f %f %f %d %d %d\n", x, y, badMaskBits, blankMaskBits, iter, rej, sys, olympic, useVariance, safe, rejection);
 
 /* 	    // Pre-reject inputs using KMM bimodality test. */