Index: /branches/czw_branch/20110406/psModules/src/imcombine/pmStack.c
===================================================================
--- /branches/czw_branch/20110406/psModules/src/imcombine/pmStack.c	(revision 31628)
+++ /branches/czw_branch/20110406/psModules/src/imcombine/pmStack.c	(revision 31629)
@@ -39,6 +39,8 @@
 # if (1)
 #define TESTING                         // Enable test output
-#define TEST_X 5745                       // x coordinate to examine
-#define TEST_Y 5331                       // y coordinate to examine
+/* #define TEST_X 5745                       // x coordinate to examine */
+/* #define TEST_Y 5331                       // y coordinate to examine */
+#define TEST_X 25
+#define TEST_Y 25
 #define TEST_RADIUS 0.5                 // Radius to examine
 # endif
@@ -102,4 +104,169 @@
     return;
 }
+
+// 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)));
+}
+
+static void KMMcalculate(const psVector *values,
+			 float *Punimodal,int *iter,
+			 float *mU, float *sU,
+			 float *pi1, float *m1, float *s1,
+			 float *pi2, float *m2, float *s2) {
+  double logL_bimodal = 0, logL_unimodal;
+  psVector *P1 = psVectorAlloc(values->n,PS_TYPE_F32);
+  psVector *P2 = psVectorAlloc(values->n,PS_TYPE_F32);
+  int i;
+/*   int debug = 0; */
+/*   if (fabs(values->data.F32[0] - -145.624954) < 1e-4) { */
+/*     debug = 1; */
+/*   } */
+  
+  // Calculate unimodal properties
+  *mU = 0;
+  *sU = 0;
+  logL_unimodal = 0;
+  for (i = 0; i < values->n; i++) { // Calculate mean
+    *mU += values->data.F32[i];
+  }
+  *mU /= values->n;
+  for (i = 0; i < values->n; i++) { // Calculate sigma
+    *sU += pow(values->data.F32[i] - *mU,2);
+  }
+  *sU = sqrt(*sU / values->n);
+  for (i = 0; i < values->n; i++) { // Calculate log likelihood
+    logL_unimodal += log(gaussian(values->data.F32[i],*mU,*sU));
+  }
+
+  // Do EM loop
+  float dL = 0;
+  float oldL = -999;
+  *iter = 0;
+  logL_bimodal = logL_unimodal;
+  *m1 = *mU - 3 * *sU;
+  *m2 = *mU + 3 * *sU;
+  *s1 = *sU / 2;
+  *s2 = *sU / 2;
+  *pi1 = 0.5;
+  *pi2 = 0.5;
+
+  float g1,g2,norm;
+  float w1,w2;
+
+  float KMM_TOLERANCE = 1e-6;
+  int KMM_MAX_ITERATIONS = 500;
+  float KMM_SMALL_NUMBER = 1e-5;
+  while (((dL > KMM_TOLERANCE)||(*iter < 3))&&(*iter < KMM_MAX_ITERATIONS)) {
+    *iter += 1;
+    dL = fabs(logL_bimodal - oldL);
+    oldL = logL_bimodal;
+/*     if (debug == 1) { */
+/*       fprintf(stderr,"KMM: %d %f %f %f %f (%f %f %f) (%f %f %f)\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
+      g1 = gaussian(values->data.F32[i],*m1,*s1);
+      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;
+    }
+    // Maximization/M-stage
+    logL_bimodal = 0;
+    w1 = 0;
+    w2 = 0;
+    for (i = 0; i < values->n; i++) { // Calculate log likelihood
+      if (!((*pi1 == 0)||(*pi2 == 0))) {
+	logL_bimodal += log(*pi1 * gaussian(values->data.F32[i],*m1,*s1) +
+			    *pi2 * gaussian(values->data.F32[i],*m2,*s2));
+      }
+    }
+    *m1 = 0;
+    *m2 = 0;
+    *s1 = 0;
+    *s2 = 0;
+    for (i = 0; i < values->n; i++) { // Calculate new means
+      *m1 += values->data.F32[i] * P1->data.F32[i];
+      *m2 += values->data.F32[i] * P2->data.F32[i];
+
+      w1 += P1->data.F32[i];
+      w2 += P2->data.F32[i];
+    }
+    *m1 /= w1;
+    *m2 /= w2;
+    for (i = 0; i < values->n; i++) { // Calculate new sigmas
+      *s1 += pow(values->data.F32[i] - *m1,2) * P1->data.F32[i];
+      *s2 += pow(values->data.F32[i] - *m2,2) * P2->data.F32[i];
+    }
+    *s1 = sqrt(*s1 / w1);
+    *s2 = sqrt(*s2 / w2);
+
+    *pi1 = w1 / values->n;
+    *pi2 = w2 / values->n;
+
+    if (!isfinite(*pi1)) { // finite checks
+      *pi1 = 0.0;
+    }
+    if (!isfinite(*pi2)) { // finite checks
+      *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;
+    }
+  } // End EM phase
+
+  // Calculate Punimodal
+  double lambda = -2.0 * (logL_unimodal - logL_bimodal);
+  int    df     = 2 + 2 * 1;
+  if (lambda > 0) {
+    *Punimodal = gsl_cdf_chisq_Q(lambda,df);
+  }
+  else {
+    *Punimodal = 1.0;
+  }
+  psFree(P1);
+  psFree(P2);
+}
+
+static void KMMFindPopular(const psVector *values, float *Punimodal, float *mean, float *sigma, float *pi) {
+  float KMM_MINIMUM_PVALUE = 0.05;
+  float mU,sU;
+  float pi1,m1,s1,pi2,m2,s2;
+  int iter;
+
+  KMMcalculate(values,Punimodal,&iter,
+	       &mU,&sU,
+	       &pi1,&m1,&s1,
+	       &pi2,&m2,&s2);
+
+  if (*Punimodal > KMM_MINIMUM_PVALUE) {
+    // Is unimodal
+    *mean = mU;
+    *sigma = sU;
+    *pi = 1.0;
+  }
+  else {
+    if (pi1 >= pi2) {
+      *mean = m1;
+      *sigma = s1;
+      *pi = pi1;
+    }
+    else {
+      *mean = m2;
+      *sigma = s2;
+      *pi = pi2;
+    }
+  }  
+}
+
+			       
 
 
@@ -287,5 +454,109 @@
     return;
 }
-
+#if (0)
+static void KMMRejectUnpopular(const psArray *inputs, int x, int y) {
+  float KMM_MINIMUM_PVALUE = 0.05;
+  float mU,sU;
+  float Punimodal,pi1,m1,s1,pi2,m2,s2;
+  int iter;
+  int j,k;
+
+  psVector *values = psVectorAlloc(inputs->n, PS_TYPE_F32);
+  k = 0;
+  for (j = 0; j < inputs->n; j++) {
+    pmStackData *data = inputs->data[j]; // Stack data of interest
+    if (!data) {
+      k++;
+      continue;
+    }
+    psImage *image = data->readout->image; // Image of interest
+    int xIn = x - data->readout->col0, yIn = y - data->readout->row0; // Coordinates on input readout
+    values->data.F32[j - k] = image->data.F32[yIn][xIn];
+  }
+  
+  KMMcalculate(values,&Punimodal,&iter,
+	       &mU,&sU,
+	       &pi1,&m1,&s1,
+	       &pi2,&m2,&s2);
+  //  fprintf(stderr,
+/*   psTrace("psModules.imcombine",3, */
+
+   CHECKPIX(x, y, 
+	  "KMM Unpopular Test: %d,%d: Puni: %g in %d",x,y,Punimodal,iter);  
+  if (Punimodal < KMM_MINIMUM_PVALUE) {
+    int i;
+    float g1,g2,norm;
+    float P1,P2;
+
+    for (i = 0; i < values->n; i++) { // Calculate probabilities for each mode
+      g1 = gaussian(values->data.F32[i],m1,s1);
+      g2 = gaussian(values->data.F32[i],m2,s2);
+      norm = (pi1 * g1 + pi2 * g2);
+      P1 = (pi1 * g1) / norm;
+      P2 = (pi2 * g2) / norm;
+
+      CHECKPIX(x, y, "KMM Unpopular Rejection: %d,%d: %f(%d): %d %f %f:(%f %f %f ) %f:(%f %f %f) rejection? %d %d\n",
+	       x, y,
+	       Punimodal,iter,
+	       i, values->data.F32[i],
+	       P1,m1,s1,pi1,
+	       P2,m2,s2,pi2,
+	       (pi1 > pi2)&&(P1 < P2),
+	       (pi1 < pi2)&&(P1 > P2));
+      if ((pi1 > pi2)&&(P1 < P2)) { // mode 1 is more popular, but this element belongs to mode 2
+	combineMarkReject(inputs,x,y,i);
+      }
+      if ((pi1 < pi2)&&(P1 > P2)) { // mode 2 is more popular, but this element belongs to mode 1
+	combineMarkReject(inputs,x,y,i);
+      }
+    }
+  }
+  psFree(values);
+  // else do nothing.
+}
+
+static void KMMRejectBright(const psArray *inputs, int x, int y) {
+  float KMM_MINIMUM_PVALUE = 0.05;
+  float mU,sU;
+  float Punimodal,pi1,m1,s1,pi2,m2,s2;
+  int iter;
+  int j;
+
+  psVector *values = psVectorAlloc(inputs->n, PS_TYPE_F32);
+  for (j = 0; j < inputs->n; j++) {
+    pmStackData *data = inputs->data[j]; // Stack data of interest
+    psImage *image = data->readout->image; // Image of interest
+    int xIn = x - data->readout->col0, yIn = y - data->readout->row0; // Coordinates on input readout
+    values->data.F32[j] = image->data.F32[yIn][xIn];
+  }
+  
+  KMMcalculate(values,&Punimodal,&iter,
+	       &mU,&sU,
+	       &pi1,&m1,&s1,
+	       &pi2,&m2,&s2);
+  if (Punimodal < KMM_MINIMUM_PVALUE) {
+    int i;
+    float g1,g2,norm;
+    float P1,P2;
+
+    for (i = 0; i < values->n; i++) { // Calculate probabilities for each mode
+      g1 = gaussian(values->data.F32[i],m1,s1);
+      g2 = gaussian(values->data.F32[i],m2,s2);
+      norm = (pi1 * g1 + pi2 * g2);
+      P1 = (pi1 * g1) / norm;
+      P2 = (pi2 * g2) / norm;
+
+      if ((m1 > m2)&&(P1 > P2)) { // m1 is larger, and this element belongs to mode 1
+	combineMarkReject(inputs,x,y,i);
+      }
+      if ((m1 < m2)&&(P1 < P2)) { // m2 is larger, and this element belongs to mode 2
+	combineMarkReject(inputs,x,y,i);
+      }
+    }
+  }
+  psFree(values);
+  // else do nothing.
+}
+#endif
 
 // Extract vectors for simple combination/rejection operations
@@ -430,5 +701,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],
@@ -571,4 +842,7 @@
     psVector *pixelLimits = buffer->limits; // Is the pixel suspect?
 
+    // KMM values;
+    float Punimodal,KMMmean,KMMsigma,KMMpi;
+    
     // Set up rejection limits
     float rej2 = PS_SQR(rej); // Rejection level squared
@@ -577,4 +851,10 @@
         // Using squared rejection limit because it's cheaper than sqrts
         double sumWeights = 0.0;
+
+	// 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);
+	
         for (int i = 0; i < num; i++) {
             sumWeights += pixelWeights->data.F32[i];
@@ -582,5 +862,6 @@
         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(sys * pixelData->data.F32[i]); */
+	  float sysVar = PS_SQR(KMMsigma);
 	    CHECKPIX(x, y, "Variance %d (%d), pixel %d,%d: %f %f %f\n", 
 		     i, pixelSources->data.U16[i], x, y, 
@@ -731,9 +1012,16 @@
           default: {
               if (useVariance) {
-                  float median = combinationWeightedOlympic(pixelData, pixelWeights,
-                                                            olympic, buffer->sort); // Median for stack
+		float median;
+		if (Punimodal > 0.05) {
+                  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
 		      CHECKPIX(x, y, "Testing input %d for pixel %d,%d: %f\n", j, x, y, diff);
@@ -742,4 +1030,9 @@
                       // pixelVariances includes the rejection limit, from above
                       float diff2 = PS_SQR(diff); // Square difference
+		      CHECKPIX(x,y, "Input %d, pixel %d,%d (%" PRIu16 "): %f %f (%f) %f %f :: %f %f %f %f\n",
+			      i, x, y, pixelSources->data.U16[j], pixelData->data.F32[j], pixelVariances->data.F32[j],
+			      1.0, pixelWeights->data.F32[j], 1.0,
+			      pixelLimits->data.F32[j], diff2, diff2 / pixelLimits->data.F32[j],worst);
+
                       if (diff2 > pixelLimits->data.F32[j]) {
                           float dev = diff2 / pixelLimits->data.F32[j]; // Deviation
@@ -980,220 +1273,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)));
-}
-
-static void KMMcalculate(const psVector *values,
-			 float *Punimodal,int *iter,
-			 float *pi1, float *m1, float *s1,
-			 float *pi2, float *m2, float *s2) {
-  double logL_bimodal = 0, logL_unimodal;
-  float mU,sU;
-  psVector *P1 = psVectorAlloc(values->n,PS_TYPE_F32);
-  psVector *P2 = psVectorAlloc(values->n,PS_TYPE_F32);
-  int i;
-  
-  // Calculate unimodal properties
-  mU = 0;
-  sU = 0;
-  logL_unimodal = 0;
-  for (i = 0; i < values->n; i++) { // Calculate mean
-    mU += values->data.F32[i];
-  }
-  mU /= values->n;
-  for (i = 0; i < values->n; i++) { // Calculate sigma
-    sU += pow(values->data.F32[i] - mU,2);
-  }
-  sU = sqrt(sU / values->n);
-  for (i = 0; i < values->n; i++) { // Calculate log likelihood
-    logL_unimodal += log(gaussian(values->data.F32[i],mU,sU));
-  }
-
-  // Do EM loop
-  float dL = 0;
-  float oldL = -999;
-  *iter = 0;
-  logL_bimodal = logL_unimodal;
-  *m1 = mU - 3 * sU;
-  *m2 = mU + 3 * sU;
-  *s1 = sU / 2;
-  *s2 = sU / 2;
-  *pi1 = 0.5;
-  *pi2 = 0.5;
-
-  float g1,g2,norm;
-  float w1,w2;
-
-  float KMM_TOLERANCE = 1e-6;
-  int KMM_MAX_ITERATIONS = 500;
-  float KMM_SMALL_NUMBER = 1e-5;
-  while (((dL > KMM_TOLERANCE)||(*iter < 3))&&(*iter < KMM_MAX_ITERATIONS)) {
-    *iter += 1;
-    dL = fabs(logL_bimodal - oldL);
-    oldL = logL_bimodal;
-
-    // Expectation/P-stage
-    for (i = 0; i < values->n; i++) { // Calculate probabilities for each mode
-      g1 = gaussian(values->data.F32[i],*m1,*s1);
-      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;
-    }
-    // Maximization/M-stage
-    logL_bimodal = 0;
-    w1 = 0;
-    w2 = 0;
-    for (i = 0; i < values->n; i++) { // Calculate log likelihood
-      if (!((*pi1 == 0)||(*pi2 == 0))) {
-	logL_bimodal += log(*pi1 * gaussian(values->data.F32[i],*m1,*s1) +
-			    *pi2 * gaussian(values->data.F32[i],*m2,*s2));
-      }
-    }
-    *m1 = 0;
-    *m2 = 0;
-    *s1 = 0;
-    *s2 = 0;
-    for (i = 0; i < values->n; i++) { // Calculate new means
-      *m1 += values->data.F32[i] * P1->data.F32[i];
-      *m2 += values->data.F32[i] * P2->data.F32[i];
-
-      w1 += P1->data.F32[i];
-      w2 += P2->data.F32[i];
-    }
-    *m1 /= w1;
-    *m2 /= w2;
-    for (i = 0; i < values->n; i++) { // Calculate new sigmas
-      *s1 += pow(values->data.F32[i] - *m1,2) * P1->data.F32[i];
-      *s2 += pow(values->data.F32[i] - *m2,2) * P2->data.F32[i];
-    }
-    *s1 = sqrt(*s1 / w1);
-    *s2 = sqrt(*s2 / w2);
-
-    *pi1 = w1 / values->n;
-    *pi2 = w2 / values->n;
-
-    if (!isfinite(*pi1)) { // finite checks
-      *pi1 = 0.0;
-    }
-    if (!isfinite(*pi2)) { // finite checks
-      *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;
-    }
-  } // End EM phase
-
-  // Calculate Punimodal
-  double lambda = -2.0 * (logL_unimodal - logL_bimodal);
-  int    df     = 2 + 2 * 1;
-  if (lambda > 0) {
-    *Punimodal = gsl_cdf_chisq_Q(lambda,df);
-  }
-  else {
-    *Punimodal = 1.0;
-  }  
-}
-
-static void KMMRejectUnpopular(const psArray *inputs, int x, int y) {
-  float KMM_MINIMUM_PVALUE = 0.05;
-  float Punimodal,pi1,m1,s1,pi2,m2,s2;
-  int iter;
-  int j;
-
-  psVector *values = psVectorAlloc(inputs->n, PS_TYPE_F32);
-  for (j = 0; j < inputs->n; j++) {
-    pmStackData *data = inputs->data[j]; // Stack data of interest
-    psImage *image = data->readout->image; // Image of interest
-    int xIn = x - data->readout->col0, yIn = y - data->readout->row0; // Coordinates on input readout
-    values->data.F32[j] = image->data.F32[yIn][xIn];
-  }
-  
-  KMMcalculate(values,&Punimodal,&iter,
-	       &pi1,&m1,&s1,
-	       &pi2,&m2,&s2);
-  //  fprintf(stderr,
-  psTrace("psModules.imcombine",3,
-	  "KMM Unpopular Test: %d,%d: Puni: %f in %d",x,y,Punimodal,iter);
-  //  CHECKPIX(x, y, "
-  
-  if (Punimodal < KMM_MINIMUM_PVALUE) {
-    int i;
-    float g1,g2,norm;
-    float P1,P2;
-
-    for (i = 0; i < values->n; i++) { // Calculate probabilities for each mode
-      g1 = gaussian(values->data.F32[i],m1,s1);
-      g2 = gaussian(values->data.F32[i],m2,s2);
-      norm = (pi1 * g1 + pi2 * g2);
-      P1 = (pi1 * g1) / norm;
-      P2 = (pi2 * g2) / norm;
-
-      CHECKPIX(x, y, "KMM Unpopular Rejection: %d,%d: %d %f %f:(%f %f %f ) %f:(%f %f %f) rejection? %d %d\n",
-	       x, y, i, values->data.F32[i],
-	       P1,m1,s1,pi1,
-	       P2,m2,s2,pi2,
-	       (pi1 > pi2)&&(P1 < P2),
-	       (pi1 < pi2)&&(P1 > P2));
-      if ((pi1 > pi2)&&(P1 < P2)) { // mode 1 is more popular, but this element belongs to mode 2
-	combineMarkReject(inputs,x,y,i);
-      }
-      if ((pi1 < pi2)&&(P1 > P2)) { // mode 2 is more popular, but this element belongs to mode 1
-	combineMarkReject(inputs,x,y,i);
-      }
-    }
-  }
-  psFree(values);
-  // else do nothing.
-}
-
-static void KMMRejectBright(const psArray *inputs, int x, int y) {
-  float KMM_MINIMUM_PVALUE = 0.05;
-  float Punimodal,pi1,m1,s1,pi2,m2,s2;
-  int iter;
-  int j;
-
-  psVector *values = psVectorAlloc(inputs->n, PS_TYPE_F32);
-  for (j = 0; j < inputs->n; j++) {
-    pmStackData *data = inputs->data[j]; // Stack data of interest
-    psImage *image = data->readout->image; // Image of interest
-    int xIn = x - data->readout->col0, yIn = y - data->readout->row0; // Coordinates on input readout
-    values->data.F32[j] = image->data.F32[yIn][xIn];
-  }
-  
-  KMMcalculate(values,&Punimodal,&iter,
-	       &pi1,&m1,&s1,
-	       &pi2,&m2,&s2);
-  if (Punimodal < KMM_MINIMUM_PVALUE) {
-    int i;
-    float g1,g2,norm;
-    float P1,P2;
-
-    for (i = 0; i < values->n; i++) { // Calculate probabilities for each mode
-      g1 = gaussian(values->data.F32[i],m1,s1);
-      g2 = gaussian(values->data.F32[i],m2,s2);
-      norm = (pi1 * g1 + pi2 * g2);
-      P1 = (pi1 * g1) / norm;
-      P2 = (pi2 * g2) / norm;
-
-      if ((m1 > m2)&&(P1 > P2)) { // m1 is larger, and this element belongs to mode 1
-	combineMarkReject(inputs,x,y,i);
-      }
-      if ((m1 < m2)&&(P1 < P2)) { // m2 is larger, and this element belongs to mode 2
-	combineMarkReject(inputs,x,y,i);
-      }
-    }
-  }
-  // else do nothing.
-}
-			       
-			       
-    
-  
   
   
@@ -1353,9 +1430,9 @@
     }
 
-    // Pre-reject inputs using KMM bimodality test.
-    if (1)  {
-/*       KMMRejectUnpopular(input,x,y); */
-      rejection = true;
-    }
+/*     // Pre-reject inputs using KMM bimodality test. */
+/*     if (1)  { */
+/* /\*       KMMRejectUnpopular(input,x,y); *\/ */
+/*       rejection = true; */
+/*     } */
    
     // Set up rejection list
@@ -1369,14 +1446,14 @@
         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);
-
-	    // Pre-reject inputs using KMM bimodality test.
-	  if (1)  {
-	    KMMRejectUnpopular(input,x,y);
-/* 	    rejection = true; */
-	  }
-	  else {
-	    KMMRejectBright(input,x,y);
-	  }
+/* 	    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. */
+/* 	  if (1)  { */
+/* 	    KMMRejectUnpopular(input,x,y); */
+/* /\* 	    rejection = true; *\/ */
+/* 	  } */
+/* 	  else { */
+/* 	    KMMRejectBright(input,x,y); */
+/* 	  } */
 
 #ifdef TESTING