Index: trunk/psModules/src/imcombine/pmStack.c
===================================================================
--- trunk/psModules/src/imcombine/pmStack.c	(revision 35455)
+++ trunk/psModules/src/imcombine/pmStack.c	(revision 36651)
@@ -43,7 +43,11 @@
 // #define TEST_X 972
 // #define TEST_Y 3213
-#define TEST_X 3289
-#define TEST_Y 4810
-#define TEST_RADIUS 0.5                 // Radius to examine
+//#define TEST_X 3289
+//#define TEST_Y 4810
+//#define TEST_RADIUS 0.5                 // Radius to examine
+//MEH -- streak-like junk md04s065i
+#define TEST_X 1129 
+#define TEST_Y 4256
+#define TEST_RADIUS 2.0                 // Radius to examine
 # endif
 
@@ -108,5 +112,5 @@
 
 // KMM functions to do bimodality rejection of pixels
-float gaussian(float x, float m, float s) {
+double gaussian(float x, float m, float s) {
   return(pow(s * sqrt(2 * M_PI),-1) * exp(-0.5 * pow( (x - m) / s, 2)));
 }
@@ -116,5 +120,6 @@
 			 float *mU, float *sU,
 			 float *pi1, float *m1, float *s1,
-			 float *pi2, float *m2, float *s2) {
+			 float *pi2, float *m2, float *s2,
+                         int xyrdebug) {
   assert(values);
   assert(values->type.type == PS_TYPE_F32);
@@ -151,4 +156,10 @@
   }
 
+  if (xyrdebug == 1) { 
+      fprintf(stderr,"KMM uni: %d %f %d (%f %f)\n", 
+              xyrdebug,logL_unimodal,discrepant_index, 
+              *mU,*sU); 
+  } 
+
   // Do EM loop
   float dL = 0;
@@ -186,5 +197,6 @@
   *pi2 = 0.5;
 
-  float g1,g2,norm;
+  //MEH -- need to be double to help avoid 0 in norm
+  double g1,g2,norm;
   float w1,w2;
 
@@ -203,4 +215,12 @@
 /* 	      *m2,*s2,*pi2); */
 /*     } */
+
+    if (xyrdebug == 1) { 
+        fprintf(stderr,"KMM EM iter: %d %f %f %f %f (%f %f %e) (%f %f %e)\n", 
+                *iter,logL_unimodal,logL_bimodal,oldL,dL, 
+                *m1,*s1,*pi1, 
+                *m2,*s2,*pi2); 
+    } 
+
     // Expectation/P-stage
     for (i = 0; i < values->n; i++) { // Calculate probabilities for each mode
@@ -208,6 +228,18 @@
       g2 = gaussian(values->data.F32[i],*m2,*s2);
       norm = (*pi1 * g1 + *pi2 * g2);
-      P1->data.F32[i] = (*pi1 * g1) / norm;
-      P2->data.F32[i] = (*pi2 * g2) / norm;
+      //MEH -- must protect denom from norm=0
+      if (norm > 0) {
+	  P1->data.F32[i] = (*pi1 * g1) / norm;
+	  P2->data.F32[i] = (*pi2 * g2) / norm;
+      } else {
+	  P1->data.F32[i] = 0.0;
+	  P2->data.F32[i] = 0.0;
+      }	 
+ 
+      if (xyrdebug == 1) {
+          fprintf(stderr,"KMM EM-P loop: %d %d %le %le %le\n",
+                         *iter,i,norm,g1,g2);
+      }
+
     }
     // Maximization/M-stage
@@ -231,4 +263,10 @@
       w1 += P1->data.F32[i];
       w2 += P2->data.F32[i];
+
+      if (xyrdebug == 1) {
+          fprintf(stderr,"KMM EM-M loop: %d %d (%f %f %f %e) (%f %f %f %e)\n",
+                         *iter,i,*m1,values->data.F32[i],w1,P1->data.F32[i],*m2,values->data.F32[i],w2,P2->data.F32[i]);
+      }
+
     }
     *m1 /= w1;
@@ -250,10 +288,18 @@
       *pi2 = 0.0;
     }
-    if (*s1 == 0) { // sigma may not be zero
-      *s1 = KMM_SMALL_NUMBER * *m1;
-    }
-    if (*s2 == 0) { // sigma may not be zero
-      *s2 = KMM_SMALL_NUMBER * *m2;
-    }
+    if (*s1 == 0) { // sigma may not be zero -- MEH -- nor <0 and need additive offset if m~0
+      *s1 = fabsf(KMM_SMALL_NUMBER * *m1) + KMM_SMALL_NUMBER;
+    }
+    if (*s2 == 0) { // sigma may not be zero 
+      *s2 = fabsf(KMM_SMALL_NUMBER * *m2) + KMM_SMALL_NUMBER;
+    }
+
+    if (xyrdebug == 1) { 
+        fprintf(stderr,"KMM EM end: %d %f %f %f %f (%f %e %e %f) (%f %e %e %f)\n", 
+                *iter,logL_unimodal,logL_bimodal,oldL,dL, 
+                *m1,*s1,*pi1,w1, 
+                *m2,*s2,*pi2,w2); 
+    } 
+
   } // End EM phase
 
@@ -267,9 +313,15 @@
     *Punimodal = 1.0;
   }
+
+  if (xyrdebug == 1) { 
+      fprintf(stderr,"KMM calc Puni: %d %f %d %f\n", 
+              xyrdebug,lambda,df,*Punimodal); 
+  } 
+
   psFree(P1);
   psFree(P2);
 }
 
-static void KMMFindPopular(const psVector *values, float *Punimodal, float *mean, float *sigma, float *pi) {
+static void KMMFindPopular(const psVector *values, float *Punimodal, float *mean, float *sigma, float *pi, int xyrdebug) {
   float KMM_MINIMUM_PVALUE = 0.05; // Should be an option.
   float mU,sU;
@@ -283,5 +335,5 @@
 	       &mU,&sU,
 	       &pi1,&m1,&s1,
-	       &pi2,&m2,&s2);
+	       &pi2,&m2,&s2,xyrdebug);
 /*   fprintf(stdout,"%g %g : %g %g %g : %g %g %g : %g %d\t", */
 /* 	  mU,sU, */
@@ -898,4 +950,6 @@
     psVector *pixelSuspects = buffer->suspects; // Is the pixel suspect?
     psVector *pixelLimits = buffer->limits; // Is the pixel suspect?
+    //MEH -- adding a debug option for TESTING xyr position but could be better..
+    int xyrdebug = 0;
 
     // KMM values;
@@ -917,5 +971,10 @@
 	// This should probably be an option
 	if (useKMM) {
-	  KMMFindPopular(pixelData,&Punimodal,&KMMmean,&KMMsigma,&KMMpi);
+#ifdef TESTING
+            if (PS_SQR(x - TEST_X) + PS_SQR(y - TEST_Y) <= PS_SQR(TEST_RADIUS)) {
+                xyrdebug = 1;
+            }
+# endif
+	  KMMFindPopular(pixelData,&Punimodal,&KMMmean,&KMMsigma,&KMMpi,xyrdebug);
 	  CHECKPIX(x,y,"KMM Popularity Contest: (%d,%d) Puni: %g Mean: %f Sigma %f Pi: %f\n",
 		   x,y,Punimodal,KMMmean,KMMsigma,KMMpi);