Index: trunk/stac/src/stacRejection.c
===================================================================
--- trunk/stac/src/stacRejection.c	(revision 3667)
+++ trunk/stac/src/stacRejection.c	(revision 5743)
@@ -74,4 +74,8 @@
     psTrace("stac.rejection", 1, "Mapping rejection masks back to source....\n");
 
+    // Vectors for calculating the mean gradient
+    psVector *grads = psVectorAlloc(nImages, PS_TYPE_F32); // Gradient for each image
+    psVector *gradsMask = psVectorAlloc(nImages, PS_TYPE_U8); // Mask for gradient vector
+
     // Transform rejection masks back to source
     psArray *inputRej = psArrayAlloc(nImages);
@@ -107,4 +111,5 @@
 	// calculate derivatives of the map, and use that as a buffer around the transformed position
 	// in the input image.
+	psStats *median = psStatsAlloc(PS_STAT_SAMPLE_MEDIAN);
 	for (int y = 0; y < nyInput; y++) {
 	    for (int x = 0; x < nxInput; x++) {
@@ -135,20 +140,28 @@
 				    (yPix >= 0) && (yPix <= ((psImage*)(inputs->data[j]))->numRows - 1)) {
 				    // Calculate the gradient
-				    meanGrads += stacGradient(inputs->data[j], xPix, yPix);
+				    grads->data.F32[j] = stacGradient(inputs->data[j], xPix, yPix);
+				    gradsMask->data.U8[j] = 0;
 				    numGrads++;
+				} else {
+				    gradsMask->data.U8[j] = 1; // Mask this one
 				}
+			    } else {
+				gradsMask->data.U8[j] = 1; // Mask this one
 			    }
 			}
 			if (numGrads > 0) {
-			    meanGrads /= (float)numGrads;
+			    (void)psVectorStats(median, grads, NULL, gradsMask, (psU8)1);
+			    meanGrads = median->sampleMedian;
 			} else {
-			    meanGrads = 0;
+			    meanGrads = 0.0;
 			}
 
 #ifdef TESTING
-			gradient->data.F32[y][x] = stacGradient(inputs->data[i], x, y) / meanGrads;
-#endif
-
-			if (stacGradient(inputs->data[i], x, y) < grad * meanGrads) {
+			//gradient->data.F32[y][x] = stacGradient(inputs->data[i], x, y) / meanGrads;
+			gradient->data.F32[y][x] = meanGrads;
+#endif
+
+			//if (stacGradient(inputs->data[i], x, y) < grad * meanGrads) {
+			if (meanGrads > grad) {
 			    mask->data.U8[y][x] = 1;
 			    nBad++;
@@ -173,4 +186,5 @@
 	    }
 	} // Iterating over pixels
+	psFree(median);
 
 #ifdef CRFLUX
@@ -194,13 +208,13 @@
 	psFits *rejmapFile = psFitsAlloc(rejmapName);
 	psFits *gradFile = psFitsAlloc(gradName);
-	if (!psFitsWriteImage(maskFile, NULL, mask, 0, NULL)) {
+	if (!psFitsWriteImage(maskFile, NULL, mask, 0)) {
 	    psErrorStackPrint(stderr, "Unable to write image: %s\n", maskName);
 	}
 	psTrace("stac", 1, "Mask image written to %s\n", maskName);
-	if (!psFitsWriteImage(rejmapFile, NULL, rejmap, 0, NULL)) {
+	if (!psFitsWriteImage(rejmapFile, NULL, rejmap, 0)) {
 	    psErrorStackPrint(stderr, "Unable to write image: %s\n", rejmapName);
 	}
 	psTrace("stac", 1, "Rejection map written to %s\n", rejmapName);
-	if (!psFitsWriteImage(gradFile, NULL, gradient, 0, NULL)) {
+	if (!psFitsWriteImage(gradFile, NULL, gradient, 0)) {
 	    psErrorStackPrint(stderr, "Unable to write image: %s\n", gradName);
 	}
@@ -218,4 +232,6 @@
     }
 
+    psFree(grads);
+    psFree(gradsMask);
 
     psFree(inCoords);