Index: trunk/stac/src/stacCombine.c
===================================================================
--- trunk/stac/src/stacCombine.c	(revision 2500)
+++ trunk/stac/src/stacCombine.c	(revision 2661)
@@ -16,6 +16,7 @@
     psStats *stats = psStatsAlloc(PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV);
     (void)psVectorStats(stats, values, masks, 1);
+    float mean = stats->sampleMean;
     psFree(stats);
-    return stats->sampleMean;
+    return mean;
 #else
     // Instead, do it ourselves
@@ -45,7 +46,8 @@
 {
     psStats *stats = psStatsAlloc(PS_STAT_SAMPLE_MEDIAN);
-    (void)psVectorStats(stats, values, masks, 1);
+    (void)psVectorStats(stats, values, masks, 0);
+    float median = stats->sampleMedian;
     psFree(stats);
-    return stats->sampleMedian;
+    return median;
 }
 
@@ -113,5 +115,10 @@
 	}
     }
-    
+
+#ifdef TESTING
+    // chi^2 image
+    psImage *chi2 = psImageAlloc(numCols, numRows, PS_TYPE_F32);
+    static int iteration = 0;		// Number of times function has been called
+#endif
     
     for (int y = 0; y < numRows; y++) {
@@ -131,13 +138,26 @@
 	    }
 	    
-	    float average = stacCombineMean(pixels, deltas, mask); // Combined value
+	    float average = stacCombineMedian(pixels, deltas, mask); // Combined value
+
+#ifdef TESTING
+	    // Calculate chi^2
+	    chi2->data.F32[y][x] = 0.0;
+	    int numGoodPix = 0;
+	    for (int i = 0; i < nImages; i++) {
+		if (mask->data.U8[i]) {
+		    chi2->data.F32[y][x] += SQUARE((pixels->data.F32[i] - average) / deltas->data.F32[i]);
+		    numGoodPix++;
+		}
+	    }
+	    chi2->data.F32[y][x] /= (float)numGoodPix;
+#endif
 	    
 	    // Rejection iterations
 	    for (int rejNum = 0; rejNum < nReject; rejNum++) {
 		float max = 0.0;
-		int maxIndex = 0;
+		int maxIndex = -1;
 		for (int i = 0; i < nImages; i++) {
-		    if (mask->data.U8[i] && (ABS(pixels->data.F32[i] - average) / deltas->data.F32[i] > max)) {
-			max = ABS(pixels->data.F32[i] - average) / deltas->data.F32[i];
+		    if (mask->data.U8[i] && ((pixels->data.F32[i] - average) / deltas->data.F32[i] > max)) {
+			max = (pixels->data.F32[i] - average) / deltas->data.F32[i];
 			maxIndex = i;
 		    }
@@ -152,5 +172,5 @@
 	    } // Rejection iterations
 	    
-	    combined->data.F32[y][x] = average;
+	    combined->data.F32[y][x] = stacCombineMean(pixels, deltas, mask);
 	}
     } // Iterating over output pixels
@@ -160,9 +180,15 @@
     if (nReject > 0) {
 	for (int i = 0; i < nImages; i++) {
-	    char rejfile[MAXCHAR];		// Filename of rejection image
+	    char rejfile[MAXCHAR];	// Filename of rejection image
 	    sprintf(rejfile,"%s.shiftrej",config->inputs->data[i]);
 	    psImageWriteSection((psImage*)((*rejected)->data[i]),0,0,0,NULL,0,rejfile);
 	}
     }
+    // Write chi^2 image
+    iteration++;
+    char chifile[MAXCHAR];		// Filename of chi^2 image
+    sprintf(chifile,"chi2_%d.fits",iteration);
+    psImageWriteSection(chi2,0,0,0,NULL,0,chifile);
+    psFree(chi2);
 #endif