Index: /tags/ipp-20110622/psModules/src/imcombine/pmStack.c
===================================================================
--- /tags/ipp-20110622/psModules/src/imcombine/pmStack.c	(revision 31721)
+++ /tags/ipp-20110622/psModules/src/imcombine/pmStack.c	(revision 31722)
@@ -37,6 +37,6 @@
 
 
-//# if (0)
-#define TESTING                         // Enable test output
+# if (0)
+#define TESTING                         Enable test output
 /* #define TEST_X 5745                       // x coordinate to examine */
 /* #define TEST_Y 5331                       // y coordinate to examine */
@@ -44,5 +44,5 @@
 #define TEST_Y 25
 #define TEST_RADIUS 0.5                 // Radius to examine
-//# endif
+# endif
 
 # ifdef TESTING
@@ -122,4 +122,6 @@
   psVector *P2 = psVectorAlloc(values->n,PS_TYPE_F32);
   int i;
+  int discrepant_index = -1;
+  double discrepant_value = 0;
 /*   int debug = 0; */
   
@@ -134,4 +136,11 @@
   for (i = 0; i < values->n; i++) { // Calculate sigma
     *sU += pow(values->data.F32[i] - *mU,2);
+
+    // Attempt to guess better starting values
+    if (pow(values->data.F32[i] - *mU,2) > discrepant_value) {
+      discrepant_value = pow(values->data.F32[i] - *mU,2);
+      discrepant_index = i;
+    }
+    
   }
   *sU = sqrt(*sU / values->n);
@@ -145,8 +154,31 @@
   *iter = 0;
   logL_bimodal = logL_unimodal;
-  *m1 = *mU - 3 * *sU;
-  *m2 = *mU + 3 * *sU;
-  *s1 = *sU / 2;
-  *s2 = *sU / 2;
+
+  if (discrepant_index == -1) {
+    *m1 = *mU - 3 * *sU;
+    *m2 = *mU + 3 * *sU;
+    *s1 = *sU / 2;
+    *s2 = *sU / 2;
+  }
+  else {
+    // This is an attempt to speed up convergence. Find the largest contributor to sigma, and set one mean
+    // to that value.  Set the other mean to the mean of all other points with this one removed.  Next,
+    // set the sigmas to be equal to each other.  Take the value of sigma to be such that a point equidistant
+    // to the initial values of the two modes is equally not believed by either mode (2.5 sigma away).
+    
+    discrepant_value = values->data.F32[discrepant_index];
+    
+    if (discrepant_value >  *mU) {
+      *m1 = ((*mU * values->n) - discrepant_value) / (values->n - 1);
+      *m2 = discrepant_value;
+    }
+    else {
+      *m1 = discrepant_value;
+      *m2 = ((*mU * values->n) - discrepant_value) / (values->n - 1);
+    }
+    *s1 = fabs((*m1 - *m2) / 5);
+    *s2 = *s1;
+  }
+    
   *pi1 = 0.5;
   *pi2 = 0.5;
@@ -157,5 +189,5 @@
   // These should be options.
   float KMM_TOLERANCE = 1e-6;
-  int KMM_MAX_ITERATIONS = 500;
+  int KMM_MAX_ITERATIONS = 30;
   float KMM_SMALL_NUMBER = 1e-5;
   while (((dL > KMM_TOLERANCE)||(*iter < 3))&&(*iter < KMM_MAX_ITERATIONS)) {
@@ -250,5 +282,15 @@
 	       &pi1,&m1,&s1,
 	       &pi2,&m2,&s2);
-
+/*   fprintf(stdout,"%g %g : %g %g %g : %g %g %g : %g %d\t", */
+/* 	  mU,sU, */
+/* 	  m1,s1,pi1, */
+/* 	  m2,s2,pi2, */
+/* 	  *Punimodal,iter); */
+/*   if (iter > 3) { */
+/*     for (int i = 0; i < values->n; i++) { */
+/*       fprintf(stdout," %f ",values->data.F32[i]); */
+/*     } */
+/*   } */
+/*   fprintf(stdout,"\n"); */
   if (*Punimodal > KMM_MINIMUM_PVALUE) {
     // Is unimodal