Index: /trunk/stac/src/Makefile
===================================================================
--- /trunk/stac/src/Makefile	(revision 2669)
+++ /trunk/stac/src/Makefile	(revision 2670)
@@ -1,5 +1,5 @@
 SHELL = /bin/sh
 CC = gcc
-CFLAGS += -O2 -g -std=c99 -I/home/mithrandir/price/psLib3/psLib/include -DTESTING -DCRFLUX
+CFLAGS += -O2 -g -std=c99 -I/home/mithrandir/price/psLib3/psLib/include -DCRFLUX
 PSLIB += -L/home/mithrandir/price/psLib3/psLib/lib/ -lpslib -lgsl -lgslcblas -lfftw3f -lsla -lcfitsio -lm
 LDFLAGS += $(PSLIB)
@@ -21,4 +21,8 @@
 		$(CC) $(CFLAGS) -o $@ calcGrad.o $(LDFLAGS) $(OPTFLAGS)
 
+median:		median.o
+		$(CC) $(CFLAGS) -o $@ median.o $(LDFLAGS) $(OPTFLAGS)
+
+
 clean:
 		-$(RM) *.o gmon.* profile.txt
@@ -37,5 +41,5 @@
 		-$(RM) test_[0-3].fits.shifterr.*
 		-$(RM) leaks.dat
-		./stac -v -k 2.5 -f 0.3 -G 0.6 testout.fits test_0.fits test_1.fits test_2.fits test_3.fits
+		./stac -v -k 3 -f 0.7 -G 0.77 testout.fits test_0.fits test_1.fits test_2.fits test_3.fits
 
 # Run profiling.
Index: /trunk/stac/src/stac.c
===================================================================
--- /trunk/stac/src/stac.c	(revision 2669)
+++ /trunk/stac/src/stac.c	(revision 2670)
@@ -58,5 +58,5 @@
 
     // Transform rejected pixels to source frame
-    psArray *rejectedSource = stacRejection(inputs, transformed, combined, maps, rejected, config);
+    psArray *rejectedSource = stacRejection(inputs, transformed, combined, maps, inverseMaps, rejected, config);
 
     // Redo transformation with the mask
Index: /trunk/stac/src/stac.h
===================================================================
--- /trunk/stac/src/stac.h	(revision 2669)
+++ /trunk/stac/src/stac.h	(revision 2670)
@@ -140,4 +140,5 @@
 		       psImage *combined, // Combined image
 		       psArray *maps,	// Maps from input to transformed image
+		       psArray *inverseMaps, // Maps from transformed to input image
 		       psArray *rejected, // Rejected images
 		       stacConfig *config // Configuration
Index: /trunk/stac/src/stacCombine.c
===================================================================
--- /trunk/stac/src/stacCombine.c	(revision 2669)
+++ /trunk/stac/src/stacCombine.c	(revision 2670)
@@ -14,5 +14,5 @@
 #if 0
     // Would like to use psVectorStats, but it doesn't have errors built in yet
-    psStats *stats = psStatsAlloc(PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV);
+    psStats *stats = psStatsAlloc(PS_STAT_SAMPLE_MEAN);
     (void)psVectorStats(stats, values, masks, 1);
     float mean = stats->sampleMean;
@@ -25,5 +25,5 @@
     int num = values->n;
     for (int i = 0; i < num; i++) {
-	if (masks->data.U8[i]) {
+	if (! masks->data.U8[i]) {
 	    sum += values->data.F32[i] / SQUARE(errors->data.F32[i]);
 	    weights += 1.0 / SQUARE(errors->data.F32[i]);
@@ -46,5 +46,5 @@
 {
     psStats *stats = psStatsAlloc(PS_STAT_SAMPLE_MEDIAN);
-    (void)psVectorStats(stats, values, masks, 0);
+    (void)psVectorStats(stats, values, masks, 1);
     float median = stats->sampleMedian;
     psFree(stats);
@@ -132,11 +132,11 @@
 		deltas->data.F32[i] = delta;
 		if ((pixel >= saturated) || (pixel <= bad) || (! isfinite(pixel)) || (! isfinite(delta))) {
-		    mask->data.U8[i] = (psU8)0; // Don't use!
+		    mask->data.U8[i] = (psU8)1; // Don't use!
 		} else {
-		    mask->data.U8[i] = (psU8)1; // Use.
+		    mask->data.U8[i] = (psU8)0; // Use.
 		}
 	    }
 	    
-	    float average = stacCombineMedian(pixels, deltas, mask); // Combined value
+	    float average = stacCombineMean(pixels, deltas, mask); // Combined value
 
 #ifdef TESTING
@@ -145,5 +145,5 @@
 	    int numGoodPix = 0;
 	    for (int i = 0; i < nImages; i++) {
-		if (mask->data.U8[i]) {
+		if (! mask->data.U8[i]) {
 		    chi2->data.F32[y][x] += SQUARE((pixels->data.F32[i] - average) / deltas->data.F32[i]);
 		    numGoodPix++;
@@ -158,5 +158,5 @@
 		int maxIndex = -1;
 		for (int i = 0; i < nImages; i++) {
-		    if (mask->data.U8[i] && ((pixels->data.F32[i] - average) / deltas->data.F32[i] > max)) {
+		    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;
@@ -165,5 +165,5 @@
 		// Reject the pixel with the maximum deviation
 		if (max > reject) {    
-		    mask->data.U8[maxIndex] = 0;
+		    mask->data.U8[maxIndex] = 1;
 		    ((psImage*)((*rejected)->data[maxIndex]))->data.F32[y][x] += 1.0;
 		    // Re-do combination following rejection
Index: /trunk/stac/src/stacRejection.c
===================================================================
--- /trunk/stac/src/stacRejection.c	(revision 2669)
+++ /trunk/stac/src/stacRejection.c	(revision 2670)
@@ -11,8 +11,7 @@
     )
 {
-    float sum = 0.0; 			// The sum of surrounding pixels
     int num = 0;
-    psVector *pixels = psVectorAlloc(9, PS_TYPE_F32); // Array of pixels
-    psVector *mask = psVectorAlloc(9, PS_TYPE_U8); // Corresponding mask
+    psVector *pixels = psVectorAlloc(8, PS_TYPE_F32); // Array of pixels
+    psVector *mask = psVectorAlloc(8, PS_TYPE_U8); // Corresponding mask
 
     // Get limits
@@ -23,24 +22,17 @@
     for (int j = yMin; j <= yMax; j++) {
 	for (int i = xMin; i <= xMax; i++) {
-	    sum += image->data.F32[j][i];
-	    pixels->data.F32[num] = image->data.F32[j][i];
-	    mask->data.U8[num] = 1;
-	    num++;
-	}
-    }
-#if 0
-    sum -= image->data.F32[y][x];
-    sum /= (float)(num-1);
-//    sum -= image->data.F32[y][x];
-    return sum / image->data.F32[y][x];
-#endif 
-
-    // Fill out the array
-    for (int i = num; i < 9; i++) {
-	pixels->data.F32[i] = 0;
-	mask->data.U8[i] = 0;
-    }
+	    if ((i != x) && (j != y)) {
+		pixels->data.F32[num] = image->data.F32[j][i];
+		mask->data.U8[num] = 0;
+		num++;
+	    }
+	}
+    }
+    pixels->n = num;
+    mask->n = num;
+
+    // Get the median
     psStats *stats = psStatsAlloc(PS_STAT_SAMPLE_MEDIAN);
-    (void)psVectorStats(stats, pixels, mask, 0);
+    (void)psVectorStats(stats, pixels, mask, 1);
     float median = stats->sampleMedian;
     psFree(stats);
@@ -85,4 +77,5 @@
 		       psImage *combined, // Combined image
 		       psArray *maps,	// Maps from input to transformed image
+		       psArray *inverseMaps, // Maps from transformed to input image
 		       psArray *rejected, // Rejected images
 		       stacConfig *config // Configuration
@@ -167,13 +160,42 @@
 #ifdef TESTING
 		rejmap->data.F32[y][x] = maskVal;
-		grad->data.F32[y][x] = stacGradient(inputs->data[i], x, y);
-#endif
-		// Check threshold, then check gradient
-		if ((maskVal > config->frac) && (stacGradient(inputs->data[i], x, y) < config->grad)) {
-		    mask->data.U8[y][x] = 1;
-		    nBad++;
+#endif
+
+		// Check gradient
+		float meanGrads = 0.0;
+		int numGrads = 0;
+		for (int j = 0; j < nImages; j++) {
+		    if (i != j) {
+			// Get coordinates for that image
+			(void)psPlaneTransformApply(inCoords, inverseMaps->data[j], outCoords);
+			int xPix = (int)(inCoords->x + 0.5);
+			int yPix = (int)(inCoords->y + 0.5);
+			if ((xPix >= 0) && (xPix <= ((psImage*)(inputs->data[j]))->numCols - 1) &&
+			    (yPix >= 0) && (yPix <= ((psImage*)(inputs->data[j]))->numRows - 1)) {
+			    // Calculate the gradient
+			    meanGrads += stacGradient(inputs->data[j], xPix, yPix);
+			    numGrads++;
+			}
+		    }
+		}
+		if (numGrads > 0) {
+		    meanGrads /= (float)numGrads;
+		} else {
+		    meanGrads = 0;
+		}
+		
+#ifdef TESTING
+		grad->data.F32[y][x] = stacGradient(inputs->data[i], x, y) / meanGrads;
+#endif
+		if (maskVal > config->frac) {
+
+		    if ((maskVal > config->frac) &&
+			(stacGradient(inputs->data[i], x, y)) < config->grad * meanGrads) {
+			mask->data.U8[y][x] = 1;
+			nBad++;
 #ifdef CRFLUX
-		    fprintf(crs, "%d %d --> %f %f %f\n", x, y, ((psImage*)(inputs->data[i]))->data.F32[y][x],
-			    maskVal, stacGradient(inputs->data[i], x, y));
+			fprintf(crs, "%d %d --> %f %f %f\n", x, y, ((psImage*)(inputs->data[i]))->data.F32[y][x],
+				maskVal, stacGradient(inputs->data[i], x, y));
+		    }
 #endif
 		} else {
@@ -189,5 +211,5 @@
 #endif
 
-	psTrace("stac.rejection", 1, "%d pixels masked in image %d.\n", nBad, i);
+	psTrace("stac.rejection", 1, "%d pixels masked in image %d\n", nBad, i);
 	// Clip the image, and convert to suitable mask format