Index: /trunk/psModules/src/imcombine/pmStack.c
===================================================================
--- /trunk/psModules/src/imcombine/pmStack.c	(revision 42090)
+++ /trunk/psModules/src/imcombine/pmStack.c	(revision 42091)
@@ -113,5 +113,5 @@
 // KMM functions to do bimodality rejection of pixels
 double gaussian(float x, float m, float s) {
-  return(pow(s * sqrt(2 * M_PI),-1) * exp(-0.5 * pow( (x - m) / s, 2)));
+    return(pow(s * sqrt(2 * M_PI),-1) * exp(-0.5 * pow( (x - m) / s, 2)));
 }
 
@@ -122,91 +122,91 @@
 			 float *pi2, float *m2, float *s2,
                          int xyrdebug) {
-  assert(values);
-  assert(values->type.type == PS_TYPE_F32);
+    assert(values);
+    assert(values->type.type == PS_TYPE_F32);
   
-  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 discrepant_index = -1;
-  double discrepant_value = 0;
+    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 discrepant_index = -1;
+    double discrepant_value = 0;
 /*   int debug = 0; */
   
-  // 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);
-
-    // 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;
-    }
+    // 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);
+
+	// 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);
-  for (i = 0; i < values->n; i++) { // Calculate log likelihood
-    logL_unimodal += log(gaussian(values->data.F32[i],*mU,*sU));
-  }
-
-  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;
-  float oldL = -999;
-  *iter = 0;
-  logL_bimodal = logL_unimodal;
-
-  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).
+    }
+    *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));
+    }
+
+    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;
+    float oldL = -999;
+    *iter = 0;
+    logL_bimodal = logL_unimodal;
+
+    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];
+	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;
-  }
+	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;
-
-  //MEH -- need to be double to help avoid 0 in norm
-  double g1,g2,norm;
-  float w1,w2;
-
-  // These should be options.
-  float KMM_TOLERANCE = 1e-6;
-  int KMM_MAX_ITERATIONS = 30;
-  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;
+    *pi1 = 0.5;
+    *pi2 = 0.5;
+
+    //MEH -- need to be double to help avoid 0 in norm
+    double g1,g2,norm;
+    float w1,w2;
+
+    // These should be options.
+    float KMM_TOLERANCE = 1e-6;
+    int KMM_MAX_ITERATIONS = 30;
+    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", */
@@ -216,124 +216,124 @@
 /*     } */
 
+	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
+	    g1 = gaussian(values->data.F32[i],*m1,*s1);
+	    g2 = gaussian(values->data.F32[i],*m2,*s2);
+	    norm = (*pi1 * g1 + *pi2 * g2);
+	    //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
+	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 and weights
+	    *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];
+
+	    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;
+	*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 -- 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
+
+    // Calculate Punimodal
+    double lambda = -2.0 * (logL_unimodal - logL_bimodal);
+    int    df     = 2 + 2 * 1; // I can't find my reference on this. 
+    if (lambda > 0) {
+	*Punimodal = gsl_cdf_chisq_Q(lambda,df);
+    }
+    else { // If lambda <= 0, then logL_unimodal > logL_bimodal, so Punimodal must be by definition 1.0
+	*Punimodal = 1.0;
+    }
+
     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); 
+	fprintf(stderr,"KMM calc Puni: %d %f %d %f\n", 
+		xyrdebug,lambda,df,*Punimodal); 
     } 
 
-    // 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);
-      //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
-    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 and weights
-      *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];
-
-      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;
-    *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 -- 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
-
-  // Calculate Punimodal
-  double lambda = -2.0 * (logL_unimodal - logL_bimodal);
-  int    df     = 2 + 2 * 1; // I can't find my reference on this. 
-  if (lambda > 0) {
-    *Punimodal = gsl_cdf_chisq_Q(lambda,df);
-  }
-  else { // If lambda <= 0, then logL_unimodal > logL_bimodal, so Punimodal must be by definition 1.0
-    *Punimodal = 1.0;
-  }
-
-  if (xyrdebug == 1) { 
-      fprintf(stderr,"KMM calc Puni: %d %f %d %f\n", 
-              xyrdebug,lambda,df,*Punimodal); 
-  } 
-
-  psFree(P1);
-  psFree(P2);
+    psFree(P1);
+    psFree(P2);
 }
 
 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;
-  float pi1,m1,s1,pi2,m2,s2;
-  int iter;
-
-  assert(values);
-  assert(values->type.type == PS_TYPE_F32);
+    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,
-	       &pi1,&m1,&s1,
-	       &pi2,&m2,&s2,xyrdebug);
+    KMMcalculate(values,Punimodal,&iter,
+		 &mU,&sU,
+		 &pi1,&m1,&s1,
+		 &pi2,&m2,&s2,xyrdebug);
 /*   fprintf(stdout,"%g %g : %g %g %g : %g %g %g : %g %d\t", */
 /* 	  mU,sU, */
@@ -347,23 +347,23 @@
 /*   } */
 /*   fprintf(stdout,"\n"); */
-  if (*Punimodal > KMM_MINIMUM_PVALUE) {
-    // Is unimodal
-    *mean = mU;
-    *sigma = sU;
-    *pi = 1.0;
-  }
-  else {
-    // Is bimodal. Select most popular mode.
-    if (pi1 >= pi2) {
-      *mean = m1;
-      *sigma = s1;
-      *pi = pi1;
+    if (*Punimodal > KMM_MINIMUM_PVALUE) {
+	// Is unimodal
+	*mean = mU;
+	*sigma = sU;
+	*pi = 1.0;
     }
     else {
-      *mean = m2;
-      *sigma = s2;
-      *pi = pi2;
-    }
-  }  
+	// Is bimodal. Select most popular mode.
+	if (pi1 >= pi2) {
+	    *mean = m1;
+	    *sigma = s1;
+	    *pi = pi1;
+	}
+	else {
+	    *mean = m2;
+	    *sigma = s2;
+	    *pi = pi2;
+	}
+    }  
 }
 
@@ -381,5 +381,5 @@
                                     const psVector *exps,      // Exposure times to combine
                                     const psVector *weights // Weights to apply
-                                    )
+    )
 {
     assert(mean);
@@ -400,4 +400,6 @@
     //     variance_combination = variance_individual / N
     // which makes sense --- the standard deviation of the combination is reduced by a factor of sqrt(N).
+    // NOTE: in 2012.07.13, the variance calculation was changed without justification to the variance
+    // appropriate to a weighted mean, while the pixel mean was kept as the average unweighted by pixel variance
 
     float sumValueWeight = 0.0;         // Sum of the value multiplied by the weight
@@ -411,6 +413,6 @@
         sumWeight += weights->data.F32[i];
         if (variances) {
-	  //            sumVarianceWeight += variances->data.F32[i] * PS_SQR(weights->data.F32[i]);
-	  sumVarianceWeight += 1 / variances->data.F32[i];
+	    //            sumVarianceWeight += variances->data.F32[i] * PS_SQR(weights->data.F32[i]);
+	    sumVarianceWeight += 1 / variances->data.F32[i];
         }
         if (exps) {
@@ -427,6 +429,6 @@
     *mean = sumValueWeight / sumWeight;
     if (var) {
-      //*var = sumVarianceWeight / PS_SQR(sumWeight);
-      *var = 1 / sumVarianceWeight;
+	//*var = sumVarianceWeight / PS_SQR(sumWeight);
+	*var = 1 / sumVarianceWeight;
     }
     if (exp) {
@@ -445,5 +447,5 @@
                                    const psVector *values, // Values to combine
                                    psVector *sortBuffer // Buffer for sorting
-                                   )
+    )
 {
     assert(values);
@@ -496,5 +498,5 @@
                                         float frac, // Fraction to discard
                                         psVector *sortBuffer // Buffer for sorting
-                                        )
+    )
 {
     int numGood = values->n;            // Number of good values
@@ -527,5 +529,5 @@
                                       int x, int y, // Pixel
                                       int source // Source image index
-                                      )
+    )
 {
     CHECKPIX(x, y, "Marking image %d, pixel %d,%d for inspection\n", source, x, y);
@@ -545,5 +547,5 @@
                                      int x, int y, // Pixel
                                      int source // Source image index
-                                     )
+    )
 {
     CHECKPIX(x, y, "Marking pixel image %d, pixel %d,%d for rejection\n", source, x, y);
@@ -561,60 +563,60 @@
 // least popular mode.
 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];
-  }
+    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);
-
-   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.
+    KMMcalculate(values,&Punimodal,&iter,
+		 &mU,&sU,
+		 &pi1,&m1,&s1,
+		 &pi2,&m2,&s2);
+
+    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.
 }
 
@@ -622,44 +624,44 @@
 // 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;
-  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];
-  }
+    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.
+    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 // End if(0) to prevent KMMReject{Unpopular|Bright} from being defined.
@@ -682,5 +684,5 @@
                            psImageMaskType badMaskBits, // Value to mask as 'bad'
                            psImageMaskType suspectMaskBits // Value to mask as 'suspect'
-                           )
+    )
 {
     // Rudimentary error checking
@@ -795,5 +797,5 @@
     // set the mask bits if nGoodBits[i] > f*numGood
     {
-	# define SUSPECT_FRACTION 0.65
+# define SUSPECT_FRACTION 0.65
 	*goodMask = 0x0000;
 	psImageMaskType value = 0x0001;
@@ -840,5 +842,5 @@
 			  int nminpix,         // Minimum number of input per pixel
                           float invTotalWeight    // Inverse of total weight for all inputs
-                          )
+    )
 {
     psVector *pixelData = buffer->pixels; // Values for the pixel of interest
@@ -944,5 +946,5 @@
                         bool useVariance, // Use variance for rejection when combining?
                         bool safe    // Combine safely?
-                        )
+    )
 {
     if (iter <= 0) {
@@ -968,5 +970,5 @@
     bool useKMM = false;
     if (num >= KMM_MINIMUM_INPUTS) {
-      useKMM = true;
+	useKMM = true;
     }
     
@@ -986,7 +988,7 @@
             }
 # 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);
+	    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);
 	}
         for (int i = 0; i < num; i++) {
@@ -994,12 +996,12 @@
         }
         for (int i = 0; i < num; i++) {
-	  // 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]);
-	  }
+	    // 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", 
@@ -1151,12 +1153,12 @@
           default: {
               if (useVariance) {
-		float median;
-		if ((useKMM)&&(Punimodal < 0.05)) {
-		  median = KMMmean;
-		}
-		else {
-                  median = combinationWeightedOlympic(pixelData, pixelWeights,
-						      olympic, buffer->sort); // Median for stack
-		}
+		  float median;
+		  if ((useKMM)&&(Punimodal < 0.05)) {
+		      median = KMMmean;
+		  }
+		  else {
+		      median = combinationWeightedOlympic(pixelData, pixelWeights,
+							  olympic, buffer->sort); // Median for stack
+		  }
 		
 		  CHECKPIX(x, y, "Flag with variance pixel %d,%d: median = %f\n", x, y, median);
@@ -1170,7 +1172,7 @@
                       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);
+			       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]) {
@@ -1435,92 +1437,124 @@
 
 bool pmStackSimpleMedianCombine(
-				pmReadout *combined,
-				psArray *input) {
-  int num = input->n;
-  //  int numCols, numRows;
-  int minInputCols, maxInputCols, minInputRows, maxInputRows; // Smallest and largest values to combine
-  int xSize, ySize;                   // Size of the output image
-
-  psArray *stack = psArrayAlloc(num); // Stack of readouts  
-  for (int i = 0; i < num; i++) {
-    //    pmStackData *data = input->data[i]; // Stack data for this input
-    pmReadout *ro = input->data[i]; // data->readout;  // Readout of interest
-    if (!ro) {
-      continue;
-    }
-    stack->data[i] = psMemIncrRefCounter(ro);
-  }    
-
-  if (!pmReadoutStackValidate(&minInputCols, &maxInputCols, &minInputRows, &maxInputRows, &xSize, &ySize,
-			      stack)) {
-    psError(psErrorCodeLast(), false, "Input stack is not valid.");
+    pmReadout *combined,
+    psArray *input) {
+    int num = input->n;
+    //  int numCols, numRows;
+    int minInputCols, maxInputCols, minInputRows, maxInputRows; // Smallest and largest values to combine
+    int xSize, ySize;                   // Size of the output image
+
+    psArray *stack = psArrayAlloc(num); // Stack of readouts  
+    for (int i = 0; i < num; i++) {
+	//    pmStackData *data = input->data[i]; // Stack data for this input
+	pmReadout *ro = input->data[i]; // data->readout;  // Readout of interest
+	if (!ro) {
+	    continue;
+	}
+	stack->data[i] = psMemIncrRefCounter(ro);
+    }    
+
+    if (!pmReadoutStackValidate(&minInputCols, &maxInputCols, &minInputRows, &maxInputRows, &xSize, &ySize,
+				stack)) {
+	psError(psErrorCodeLast(), false, "Input stack is not valid.");
+	psFree(stack);
+	return false;
+    }
+
+    psVector *pixelData = psVectorAlloc(input->n,PS_TYPE_F32);
+    psVector *pixelMask = psVectorAlloc(input->n,PS_TYPE_VECTOR_MASK);
+    psStats  *stats     = psStatsAlloc(PS_STAT_ROBUST_MEDIAN);
+
+    for (int y = minInputRows; y < maxInputRows; y++) {
+	for (int x = minInputCols; x < maxInputCols; x++) {
+	    for (int i = 0; i < input->n; i++) {
+		pmReadout *ro  = stack->data[i];
+		psImage *image = ro->image;
+		pixelMask->data.PS_TYPE_VECTOR_MASK_DATA[i] = 0;
+		pixelData->data.F32[i] = image->data.F32[y][x];
+		if (isfinite(image->data.F32[y][x])&&
+		    (fabs(image->data.F32[y][x]) < 1e5)) {
+		    pixelMask->data.PS_TYPE_VECTOR_MASK_DATA[i] = 0;
+		}
+		else {
+		    pixelMask->data.PS_TYPE_VECTOR_MASK_DATA[i] = 1;
+		}
+#if (0)
+		if ((x == 59)&&(y > 40)&&(y < 50)) {
+		    fprintf(stderr,"%d %d %d %d %g\n",
+			    x,y,i,pixelMask->data.PS_TYPE_VECTOR_MASK_DATA[i],pixelData->data.F32[i]);
+		}
+#endif
+	    }
+	    if (!psVectorStats(stats,pixelData,NULL,pixelMask,1)) {
+		psError(PS_ERR_UNKNOWN, false, "Unable to calculate median");
+		psFree(stats);
+		psFree(pixelData);
+		psFree(pixelMask);
+		psFree(stack);
+		return(false);
+	    }
+	    combined->image->data.F32[y][x] = stats->robustMedian;
+	    if (combined->mask) {
+		combined->mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] = 0;
+	    }
+#if (0)
+	    if ((x == 59)&&(y > 40)&&(y < 50)) {
+		fprintf(stderr,"%d %d %d %d %g\n",
+			x,y,-1,combined->mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x],
+			combined->image->data.F32[y][x]);
+	    }
+#endif
+	}
+    }
+  
+    psFree(stats);
+    psFree(pixelData);
+    psFree(pixelMask);
     psFree(stack);
-    return false;
-  }
-
-  psVector *pixelData = psVectorAlloc(input->n,PS_TYPE_F32);
-  psVector *pixelMask = psVectorAlloc(input->n,PS_TYPE_VECTOR_MASK);
-  psStats  *stats     = psStatsAlloc(PS_STAT_ROBUST_MEDIAN);
-
-  for (int y = minInputRows; y < maxInputRows; y++) {
-    for (int x = minInputCols; x < maxInputCols; x++) {
-      for (int i = 0; i < input->n; i++) {
-	pmReadout *ro  = stack->data[i];
-	psImage *image = ro->image;
-	pixelMask->data.PS_TYPE_VECTOR_MASK_DATA[i] = 0;
-	pixelData->data.F32[i] = image->data.F32[y][x];
-	if (isfinite(image->data.F32[y][x])&&
-	    (fabs(image->data.F32[y][x]) < 1e5)) {
-	  pixelMask->data.PS_TYPE_VECTOR_MASK_DATA[i] = 0;
-	}
-	else {
-	  pixelMask->data.PS_TYPE_VECTOR_MASK_DATA[i] = 1;
-	}
-#if (0)
-      if ((x == 59)&&(y > 40)&&(y < 50)) {
-	  fprintf(stderr,"%d %d %d %d %g\n",
-		  x,y,i,pixelMask->data.PS_TYPE_VECTOR_MASK_DATA[i],pixelData->data.F32[i]);
-	}
-#endif
-      }
-      if (!psVectorStats(stats,pixelData,NULL,pixelMask,1)) {
-	psError(PS_ERR_UNKNOWN, false, "Unable to calculate median");
-	psFree(stats);
-	psFree(pixelData);
-	psFree(pixelMask);
-	psFree(stack);
-	return(false);
-      }
-      combined->image->data.F32[y][x] = stats->robustMedian;
-      if (combined->mask) {
-	combined->mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] = 0;
-      }
-#if (0)
-      if ((x == 59)&&(y > 40)&&(y < 50)) {
-	fprintf(stderr,"%d %d %d %d %g\n",
-		x,y,-1,combined->mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x],
-		combined->image->data.F32[y][x]);
-      }
-#endif
-    }
-  }
-  
-  psFree(stats);
-  psFree(pixelData);
-  psFree(pixelMask);
-  psFree(stack);
-  return (true);
-}
-
-# define MIN_GOOD_PERCENTILE 2
+    return (true);
+}
+
+# define SUSPECT_FRACTION 0.65
+
+// Comparison and swap functions for sorting values directly
+#define SORT_VV_COMPARE(A,B) (pixelData->data.F32[A] < pixelData->data.F32[B])
+#define SORT_VV_SWAP(TYPE,A,B) {					\
+	if (A != B) {							\
+	    psF32 tempVal = pixelData->data.F32[A];			\
+	    pixelData->data.F32[A] = pixelData->data.F32[B];		\
+	    pixelData->data.F32[B] = tempVal;				\
+	    psF32 tempVar = pixelVariances->data.F32[A];		\
+	    pixelVariances->data.F32[A] = pixelVariances->data.F32[B];	\
+	    pixelVariances->data.F32[B] = tempVar;			\
+	    if (expTime) {						\
+		psF32 tempExp = expTime->data.F32[A];			\
+		expTime->data.F32[A] = expTime->data.F32[B];		\
+		expTime->data.F32[B] = tempExp;				\
+	    }								\
+	}								\
+    }
+
+// this macro uses the macros above which assume pixelData, pixelVariances, expTime
+#define SORT_VALUES(NVALUES) { PSSORT(NVALUES, SORT_VV_COMPARE, SORT_VV_SWAP, F32); }
+
+#define ESCAPE	{							\
+  combined->image->data.F32[y][x] = NAN;				\
+  combined->variance->data.F32[y][x] = NAN;				\
+  combined->mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] = blankMaskBits;	\
+  if (expmaps) {							\
+    expmaps->image->data.F32[y][x] = 0.0;				\
+    expmaps->variance->data.F32[y][x] = 0.0;				\
+    expmaps->mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] = 0;		\
+  } continue; }
+
 bool pmStackCombineByPercentile(
-    pmReadout *combined,
-    psArray *stackData,
-    psF64 minRange,
-    psF64 maxRange,
+    pmReadout *combined,		// output stacked readout
+    pmReadout *expmaps,			// output exposure map information
+    psArray *stackData,			// input exposures
+    psF64 rejectFraction,               // outlier fraction of pixels to reject
     psImageMaskType badMaskBits, 	// treat these bits as 'bad'
     psImageMaskType suspectMaskBits,	// treat these bits as 'suspect'
     psImageMaskType blankMaskBits       // use this mask value for pixels missing input data (distinguish between Ninput = 0 and Ngood = 0?)
-) {
+    ) {
     int minInputCols, maxInputCols, minInputRows, maxInputRows; // Smallest and largest values to combine
     int xSize, ySize;                   // Size of the output image
@@ -1530,8 +1564,10 @@
     for (int i = 0; i < stackData->n; i++) {
         pmStackData *data = stackData->data[i]; // Stack data for this input
-        pmReadout *ro = data->readout;  // Readout of interest
-        stackReadouts->data[i] = psMemIncrRefCounter(ro);
-    }
-
+	stackReadouts->data[i] = NULL;
+	if (!data) continue;
+	pmReadout *ro = data->readout;  // Readout of interest
+	if (!ro) continue;
+	stackReadouts->data[i] = psMemIncrRefCounter(ro); // need to bump the counter since the free below will decrement
+    }
     if (!pmReadoutStackValidate(&minInputCols, &maxInputCols, &minInputRows, &maxInputRows, &xSize, &ySize, stackReadouts)) {
 	psError(psErrorCodeLast(), false, "Input stack is not valid.");
@@ -1541,6 +1577,19 @@
     psFree(stackReadouts);
 
-    psVector *pixelData = psVectorAlloc(stackData->n, PS_TYPE_F32);
-    psVector *pixelMask = psVectorAlloc(stackData->n, PS_TYPE_VECTOR_MASK);
+    // make sure the output readout matches the inputs, and set to blank by default
+    pmReadoutUpdateSize(combined, minInputCols, minInputRows, xSize, ySize, true, true, blankMaskBits);
+    if (expmaps) {
+	// if we are generating the expmaps, update to match the images, set blank mask areas to 0
+	pmReadoutUpdateSize(expmaps, minInputCols, minInputRows, xSize, ySize, expmaps->mask != NULL, expmaps->variance != NULL, 0);
+    }
+   
+    psVector *pixelData      = psVectorAlloc(stackData->n, PS_TYPE_F32);
+    psVector *pixelVariances = psVectorAlloc(stackData->n, PS_TYPE_F32);
+
+    // if we are asking for the exptime maps, generate a storage vector expTime
+    psVector *expTime        = expmaps && expmaps->image ? psVectorAlloc(stackData->n, PS_TYPE_F32) : NULL;
+
+    int nGoodBits[16]; // accumulate the good pixel bits here for fuzzy logic
+    psAssert (sizeof(psImageMaskType) == 2, "psImageMaskType is not the expected size");
 
     for (int y = minInputRows; y < maxInputRows; y++) {
@@ -1548,91 +1597,203 @@
 
 	    int nGood = 0; 
+	    memset (nGoodBits, 0, 16*sizeof(int));
 	    for (int i = 0; i < stackData->n; i++) {
 
 		pmStackData *data = stackData->data[i]; // Stack data for this input
+		if (!data) continue;
 		pmReadout *ro = data->readout;  // Readout of interest
-		psImage *image = ro->image;
-		psImage *mask = ro->mask;
+		if (!ro) continue;
+
+		psAssert (ro->mask,     "must must exist, but does not");
+		psAssert (ro->variance, "variance must exist, but does not");
+
+		psImage *image    = ro->image;
+		psImage *variance = ro->variance;
+		psImage *mask     = ro->mask;
 
 		int xIn = x - data->readout->col0;
 		int yIn = y - data->readout->row0; // Coordinates on input readout
 
+		// skip obviously bad input data
 		if (!isfinite(image->data.F32[yIn][xIn])) continue;
-		if (fabs(image->data.F32[yIn][xIn]) > 1e5) continue;
-		// XXX need to test the input mask as well here
-		if (mask->data.PS_TYPE_IMAGE_MASK_DATA[yIn][xIn] & badMaskBits);
-
-		// XXX save the count of suspect bits
-
-# if (0)
+		if (mask->data.PS_TYPE_IMAGE_MASK_DATA[yIn][xIn] & badMaskBits) continue;
+		if (mask->data.PS_TYPE_IMAGE_MASK_DATA[yIn][xIn] & suspectMaskBits) continue;
+
 		// count the number of times a given mask bit is set in the input pixels.
 		// NOTE: since we have explicitly skipped the pixels with any bad bits, these are only
 		// the suspect bits (nGoodBits is a bit of a misnomer: it is more like 'nSuspectBitsForGoodInputs'
-		if (1) { 
-		  psImageMaskType value = 0x0001;
-		  for (int nbit = 0; nbit < 16; nbit ++) {
+		// NOTE: skip the full bit-by-bit check if we know the mask byte is empty
+		psImageMaskType value = 0x0001;
+		for (int nbit = 0; mask->data.PS_TYPE_IMAGE_MASK_DATA[yIn][xIn] && (nbit < 16); nbit ++) {
 		    if (mask->data.PS_TYPE_IMAGE_MASK_DATA[yIn][xIn] & value) {
-		      // XXX nGoodBits[nbit] ++;
+			nGoodBits[nbit] ++;
 		    }
 		    value <<= 1;
-		  }
 		}
+
+		// accumulate pixel data and variance values:
+		pixelData->data.F32[nGood] = image->data.F32[yIn][xIn];
+		pixelVariances->data.F32[nGood] = variance->data.F32[yIn][xIn];
+
+		// accumulate exposure times if required
+		if (expTime)  { expTime->data.F32[nGood] = data->exp; }
+		nGood ++;
+	    }
+# define MIN_GOOD_PERCENTILE 3
+
+	    if (nGood < MIN_GOOD_PERCENTILE) ESCAPE;
+	    
+	    pixelData->n = nGood;
+
+	    // sort pixelData, pixelVariance, expTime (if it exists)
+	    SORT_VALUES (pixelData->n);
+
+	    // we now have a sorted vector of values.  We can make a very coarse outlier
+	    // cut based on the median and interquartile range.  This will let us reject
+	    // some extreme outliers that bias the signal otherwise.
+
+	    int Nlo = 0;     
+	    int Nhi = nGood; 
+
+	    if (nGood >= 5) {
+	      int midPoint = nGood / 2;
+	      float rawMedian = (nGood % 2) ? pixelData->data.F32[midPoint] : 0.5*(pixelData->data.F32[midPoint] + pixelData->data.F32[midPoint-1]);
+	      
+	      // XXX measure interquartile range
+	      int P25 = 0.25*nGood;
+	      int P75 = 0.75*nGood;
+	      float rawSigma = 0.74*(pixelData->data.F32[P75] - pixelData->data.F32[P25]);
+	      float minThresh = rawMedian - 5.0*rawSigma;
+	      float maxThresh = rawMedian + 5.0*rawSigma;
+
+	      // find the entries which are in the range
+	      // these should be safe since minThresh and maxThresh are guaranteed to contain the median
+	      while (pixelData->data.F32[Nlo    ] < minThresh) { Nlo ++; }
+	      while (pixelData->data.F32[Nhi - 1] > maxThresh) { Nhi --; }
+	    }
+
+	    // if we did not clip above (either nGood < 5 or no rejection), Nlo will be 0, Nhi will be nGood
+	    // In either case, Nlo is the offset of the first unclipped point.
+	    int nGoodClip = Nhi - Nlo;
+
+	    int isTest = false;
+	    isTest = isTest || ((x == 4743) && (y == 2903));
+	    isTest = isTest || ((x == 4953) && (y == 2919));
+	    isTest = isTest || ((x == 4751) && (y == 2725));
+	    isTest = false;
+
+	    if (isTest) {
+	      char testname[256];
+	      snprintf (testname, 256, "testpix.%04d.%04d.txt", (int) x, (int) y);
+	      FILE *f = fopen (testname, "w");
+	      int fd = fileno (f);
+
+	      fprintf (f, "# nGood: %d, Nlo: %d, Nhi: %d, nGoodClip: %d\n", nGood, Nlo, Nhi, nGoodClip);
+	      p_psVectorPrint (fd, pixelData, "pixelData");
+	      fclose (f);
+	    }
+	    // rather than define a min and max value,
+	    // what we really want is a symmetric selection about the middle,
+	    // or the output is biased.  If N % 2 = 1, then 
+	    
+	    // we are going to exclude rejectFraction*nGood measurements.  But the
+	    // rejection needs to be symmetric
+
+	    // 'AT_LEAST' means we reject at least 'rejectFraction' of values, but it could
+	    // be a higher percentage for smaller numbers of inputs.
+# define AT_LEAST 1
+# define MIN_GOOD_PERCENTILE 3
+# if (AT_LEAST)
+	    int Ns = MIN(MAX(1, 0.5*rejectFraction * nGoodClip), nGoodClip) + Nlo;
+	    int Ne = nGoodClip - Ns + Nlo;
+	    int Npt = Ne - Ns;
+# else
+	    int Ns = MIN(MAX(0, 0.5*rejectFraction * nGoodClip), nGoodClip) + Nlo;
+	    int Ne = nGoodClip - Ns + Nlo;
+	    int Npt = Ne - Ns;
 # endif
-
-		pixelData->data.F32[nGood] = image->data.F32[yIn][xIn];
-		nGood ++;
-	    }
-	    pixelData->n = nGood;
-	    psVectorSortInPlace (pixelData);
-
-	    if (nGood < MIN_GOOD_PERCENTILE) {
-		combined->image->data.F32[y][x] = NAN;
-		combined->variance->data.F32[y][x] = NAN;
-		combined->mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] = 1;
-		continue;
-	    }
-
-# if (0)
-# define SUSPECT_FRACTION 0.65
-	    // set the mask bits if nGoodBits[i] > f*numGood
-	    if (0) {
-	      psImageMaskType value = 0x0001;
-	      for (int nbit = 0; nbit < 16; nbit ++) {
-		if (nGoodBits[nbit] > SUSPECT_FRACTION*numGood) {
-		  *goodMask |= value;
+	    if ((Npt < 1) || (nGood < MIN_GOOD_PERCENTILE)) ESCAPE;
+
+	    // Set the given (suspect) mask bit if nGoodBits[i] > f*nGood in other words
+	    // if more than 65% of the good inputs had one of these bits set, then we
+	    // should set that bit in the output mask.  Note that this analysis counts the
+	    // mask bits of pixels rejected by the clipping above.
+	    psImageMaskType value = 0x0001;
+	    psImageMaskType outputMask = 0x0000;
+	    for (int nbit = 0; nbit < 16; nbit ++) {
+		if (nGoodBits[nbit] > SUSPECT_FRACTION*nGood) {
+		    outputMask |= value;
 		}
 		value <<= 1;
-	      }
-	    }
-# endif
-
-	    int Ns = MIN(MAX(0, minRange * nGood),nGood); // e.g., 0.1 * 50 = 5
-	    int Ne = MIN(MAX(0, maxRange * nGood),nGood); // e.g., 0.9 * 50 = 45 
-	    int Npt = Ne - Ns;
+	    }
 
 	    float sum = 0.0;
+	    float varSum = 0.0;
 	    for (int n = Ns; n < Ne; n++) {
 		sum += pixelData->data.F32[n];
+		varSum += pixelVariances->data.F32[n];
 	    }
 	    float mean = sum / (float) Npt;
-
-	    float varSum = 0.0;
-	    for (int n = Ns; n < Ne; n++) {
-		varSum += SQ(pixelData->data.F32[n] - mean);
-	    }
+	    float varValue = varSum / (float) (nGoodClip*nGoodClip);
+
+	    // alternative: calculate the stdev of the pixel values
+	    // float varSum = 0.0;
+	    // for (int n = Ns; n < Ne; n++) {
+	    // 	varSum += SQ(pixelData->data.F32[n] - mean);
+	    // }
 	    // variance on the mean (stdev / sqrt(N))^2
-	    float varValue = varSum / (Npt - 1) / Npt;
-
-	    // XXX this is probably not the correct variance to save
-	    // the predicted variance should be the 1 / sum (1/var)
+
+	    // the reported variance values can be extremely high / wrong.
+	    // if we have enough measurements, let's just use the interquartile range
+	    // of the data to estimate the per-pixel variance.  NOTE: this is not valid
+	    // if the inputs have been significantly smoothed.  In that case we need
+	    // to include the covariance explicitly.  But this algorithm should be used
+	    // without convolution.
+	    // XXX How do we choose the cutoff here?
+	    if (nGoodClip >= 9) {
+	      // Measure interquartile range
+	      int P25 = 0.25*nGoodClip + Nlo;
+	      int P75 = 0.75*nGoodClip + Nlo;
+	      float rawSigma = 0.74*(pixelData->data.F32[P75] - pixelData->data.F32[P25]);
+	      varValue = PS_MIN(9e7, PS_SQR(rawSigma) / (float) nGoodClip); // sigma_mean = sigma_meas / sqrt(Nmeas) -> var_mean = var_meas / Nmeas
+	      // XXX the upper limit of 9e7 is set to match the output
+	      // format (STK_UNIONS) which can only represent values up to that limit.
+	      // perhaps it would be better to saturate the output image in psFits
+	      // rather than here.  
+	    }
+
+	    // Note: since we are calculating the average of a subset of a sorted
+	    // list of values, the denominator should not be the number of measurements
+	    // in the calculation above (Npt): in the extreme case of a median, we would
+	    // have a single value (Npt = 1), but the variance of a median is only ~1.4 x
+	    // the variance of the average / sqrt(N).  We should use nGood (the total number
+	    // of values in the sorted list), but the variance should be scaled by a factor
+	    // which depends on the fraction of values included.  
+
+	    // this coefficient varies between 1.4 (for pure median) and 1.05 for 68% range.
 
 	    combined->image->data.F32[y][x] = mean;
 	    combined->variance->data.F32[y][x] = varValue;
-	    combined->mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] = 0;
-	}
-    }
-  
+	    combined->mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] = outputMask;
+
+	    // The exposure time of interest should be the total number of values, after
+	    // rejection of known bad measurements, not the sorted and clipped number.
+	    // Note that if we were to take the median, the relevant exposure time would
+	    // still be the total of all inputs, not the single exposure for which the
+	    // median was generated.
+
+	    if (expTime) { 
+		float sum = 0.0;
+		for (int n = Nlo; n < Nhi; n++) {
+		    sum += expTime->data.F32[n];
+		}
+		expmaps->image->data.F32[y][x] = sum;
+	    }
+	    if (expmaps) { expmaps->mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] = nGoodClip; }
+	}
+    }
     psFree(pixelData);
-    psFree(pixelMask);
+    psFree(pixelVariances);
+    psFree(expTime);
 
     return (true);
Index: /trunk/psModules/src/imcombine/pmStack.h
===================================================================
--- /trunk/psModules/src/imcombine/pmStack.h	(revision 42090)
+++ /trunk/psModules/src/imcombine/pmStack.h	(revision 42091)
@@ -49,7 +49,7 @@
 bool pmStackCombineByPercentile(
     pmReadout *combined,
+    pmReadout *expmaps,
     psArray *input,
-    psF64 minRange,
-    psF64 maxRange,
+    psF64 rejectFraction,
     psImageMaskType badMaskBits, 	// treat these bits as 'bad'
     psImageMaskType suspectMaskBits,	// treat these bits as 'suspect'