Index: trunk/stac/src/stacRejection.c
===================================================================
--- trunk/stac/src/stacRejection.c	(revision 2661)
+++ 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