Index: trunk/stac/src/stacRejection.c
===================================================================
--- trunk/stac/src/stacRejection.c	(revision 2500)
+++ trunk/stac/src/stacRejection.c	(revision 2661)
@@ -2,4 +2,81 @@
 #include "pslib.h"
 #include "stac.h"
+
+
+#define MAX(x,y) ((x) > (y) ? (x) : (y))
+#define MIN(x,y) ((x) < (y) ? (x) : (y))
+
+float stacGradient(psImage *image,	// Input for which to measure the gradient
+		   int x, int y		// Coordinates at which to measure the gradient
+    )
+{
+    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
+
+    // Get limits
+    int xMin = MAX(x - 1, 0);
+    int xMax = MIN(x + 1, image->numCols - 1);
+    int yMin = MAX(y - 1, 0);
+    int yMax = MIN(y + 1, image->numRows - 1);
+    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;
+    }
+    psStats *stats = psStatsAlloc(PS_STAT_SAMPLE_MEDIAN);
+    (void)psVectorStats(stats, pixels, mask, 0);
+    float median = stats->sampleMedian;
+    psFree(stats);
+    psFree(pixels);
+    psFree(mask);
+    return median / image->data.F32[y][x];
+}
+
+#if 0
+float stacGradient(psImage *image,	// Input for which to measure the gradient
+		   int x, int y		// Coordinates at which to measure the gradient
+		   )
+{
+    float sum = 0.0; 			// The sum of surrounding pixels
+    float maxDiff = 0.0;
+    int num = 0;
+    // Get limits
+    int xMin = MAX(x - 1, 0);
+    int xMax = MIN(x + 1, image->numCols - 1);
+    int yMin = MAX(y - 1, 0);
+    int yMax = MIN(y + 1, image->numRows - 1);
+    for (int j = yMin; j <= yMax; j++) {
+	for (int i = xMin; i <= xMax; i++) {
+	    if (image->data.F32[j][i] - image->data.F32[y][x] < maxDiff) {
+		maxDiff = image->data.F32[j][i] - image->data.F32[y][x];
+	    }
+	    sum += image->data.F32[j][i];
+	    num++;
+	}
+    }
+    sum -= image->data.F32[y][x];
+    sum /= (float)(num-1);
+    sum -= image->data.F32[y][x];
+//    return sum / image->data.F32[y][x];
+    return maxDiff / image->data.F32[y][x];
+}   
+#endif
 
 
@@ -55,9 +132,25 @@
 	int nxInput = ((psImage*)(inputs->data[i]))->numCols;
 	int nyInput = ((psImage*)(inputs->data[i]))->numRows;
-	psImage *mask = psImageAlloc(nxInput, nyInput, PS_TYPE_F32); // The mask in the source frame
+	psImage *mask = psImageAlloc(nxInput, nyInput, PS_TYPE_U8); // The pixel mask for the input
+#ifdef TESTING
+	psImage *rejmap = psImageAlloc(nxInput, nyInput, PS_TYPE_F32); // The rejections in the source frame
+	psImage *grad = psImageAlloc(nxInput, nyInput, PS_TYPE_F32);	// The gradient image
+#endif
 	psImage *reject = rejected->data[i]; // Pull out the mask in the output frame
 	psPlaneTransform *map = maps->data[i]; // The map from input to output
 
 	psTrace("stac.rejection", 3, "Transforming rejection mask %d....\n", i);
+	int nBad = 0;			// Number of bad pixels
+
+#ifdef CRFLUX
+	// Set up CR output
+	FILE *crs = NULL;		// File for outputting details of rejected pixels
+	char crfile[MAXCHAR];	// Filename
+	sprintf(crfile,"%s.cr",config->inputs->data[i]);
+	if ((crs = fopen(crfile, "w")) == NULL) {
+	    psError("stac.rejection","Unable to open file for detailed output, %s\n");
+	    return NULL;
+	}
+#endif
 
 	// Transform the mask
@@ -67,29 +160,56 @@
 	for (int y = 0; y < nyInput; y++) {
 	    for (int x = 0; x < nxInput; x++) {
-    		inCoords->x = (double)x + 0.5;
-		inCoords->y = (double)y + 0.5;
+    		inCoords->x = (double)x;
+		inCoords->y = (double)y;
 		(void)psPlaneTransformApply(outCoords, map, inCoords);
-		mask->data.F32[y][x] = psImagePixelInterpolate(reject, outCoords->x, outCoords->y,
+		float maskVal = (float)psImagePixelInterpolate(reject, outCoords->x, outCoords->y,
 							       NULL, 0, 0.0, PS_INTERPOLATE_BILINEAR);
+#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++;
+#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));
+#endif
+		} else {
+		    mask->data.U8[y][x] = 0;
+		}
 
 	    }
-	}
+	} // Iterating over pixels
+
+#ifdef CRFLUX
+	// Close file used for detailed output
+	fclose(crs);
+#endif
+
+	psTrace("stac.rejection", 1, "%d pixels masked in image %d.\n", nBad, i);
+	// Clip the image, and convert to suitable mask format
 
 #ifdef TESTING
 	// Write error image out to check
-	char maskfile[MAXCHAR];	// Filename of mask image
+	char maskfile[MAXCHAR];		// Filename of mask image
+	char rejmapfile[MAXCHAR]; 	// Filename of rejection image
+	char gradfile[MAXCHAR];		// Filename of gradient image
 	sprintf(maskfile,"%s.mask",config->inputs->data[i]);
+	sprintf(rejmapfile,"%s.rejmap",config->inputs->data[i]);
+	sprintf(gradfile,"%s.grad",config->inputs->data[i]);
 	psImageWriteSection(mask,0,0,0,NULL,0,maskfile);
-#endif
-
-	// Clip the image, and convert to suitable mask format
-	(void)psImageClip(mask, config->frac, 0.0, config->frac, 1.0);
-	psImage *maskU8 = psImageCopy(NULL, mask, PS_TYPE_U8);
-	psFree(mask);
-    
+	psImageWriteSection(rejmap,0,0,0,NULL,0,rejmapfile);
+	psImageWriteSection(grad,0,0,0,NULL,0,gradfile);
+	psFree(rejmap);
+	psFree(grad);
+#endif
+
 	// Stuff into the array
-	inputRej->data[i] = maskU8;
-
-    }
+	inputRej->data[i] = mask;
+
+    }
+
 
     psFree(inCoords);