Index: /trunk/stac/src/Makefile
===================================================================
--- /trunk/stac/src/Makefile	(revision 2660)
+++ /trunk/stac/src/Makefile	(revision 2661)
@@ -1,7 +1,7 @@
 SHELL = /bin/sh
 CC = gcc
-CFLAGS = -O2 -g -std=c99 -I/home/mithrandir/price/psLib3/psLib/include # -DTESTING
-PSLIB = -L/home/mithrandir/price/psLib3/psLib/lib/ -lpslib -lgsl -lgslcblas -lfftw3f -lsla -lcfitsio -lm
-LDFLAGS = $(PSLIB)
+CFLAGS += -O2 -g -std=c99 -I/home/mithrandir/price/psLib3/psLib/include -DTESTING -DCRFLUX
+PSLIB += -L/home/mithrandir/price/psLib3/psLib/lib/ -lpslib -lgsl -lgslcblas -lfftw3f -lsla -lcfitsio -lm
+LDFLAGS += $(PSLIB)
 
 OBJECTS = stac.o stacConfig.o stacRead.o stacErrorImages.o stacTransform.o stacCheckMemory.o \
@@ -18,4 +18,7 @@
 		$(CC) $(CFLAGS) -o $@ $(OBJECTS) $(LDFLAGS) $(OPTFLAGS)
 
+calcGrad:	calcGrad.o
+		$(CC) $(CFLAGS) -o $@ calcGrad.o $(LDFLAGS) $(OPTFLAGS)
+
 clean:
 		-$(RM) *.o gmon.* profile.txt
@@ -27,11 +30,12 @@
 		-$(RM) testout.fits
 		-$(RM) testout.fits.pre
+		-$(RM) chi2_*.fits
 		-$(RM) test_[0-3].fits.err
 		-$(RM) test_[0-3].fits.mask
-		-$(RM) test_[0-3].fits.shift.?
+		-$(RM) test_[0-3].fits.shift.*
 		-$(RM) test_[0-3].fits.shiftrej
 		-$(RM) test_[0-3].fits.shifterr.*
 		-$(RM) leaks.dat
-		./stac -v testout.fits test_0.fits test_1.fits test_2.fits test_3.fits
+		./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
 
 # Run profiling.
@@ -49,8 +53,8 @@
 optimise:	clean
 		-$(RM) $(TARGET)
-		export OPTFLAGS=-fprofile-arcs ; $(MAKE) test
+		export OPTFLAGS=-fprofile-generate ; $(MAKE) test
 		-$(RM) $(TARGET)
 		$(MAKE) clean
-		export OPTFLAGS=-fbranch-probabilities ; $(MAKE) $(TARGET)
+		export OPTFLAGS=-fbranch-use ; $(MAKE) $(TARGET)
 		-$(RM) *.da
 
Index: /trunk/stac/src/stac.h
===================================================================
--- /trunk/stac/src/stac.h	(revision 2660)
+++ /trunk/stac/src/stac.h	(revision 2661)
@@ -27,4 +27,5 @@
     float frac;				// Fraction of input pixel that must be masked before the pixel is
 					// considered bad
+    float grad;				// Multiplier of the gradient
 } stacConfig;
 
@@ -129,4 +130,9 @@
 // stacRejection.c
 
+// Return the relative gradient for a given pixel
+float stacGradient(psImage *image,	// Input for which to measure the gradient
+		   int x, int y		// Coordinates at which to measure the gradient
+    );
+
 // Transform the rejection masks back to the source frame
 psArray *stacRejection(psArray *inputs,	// Input images
Index: /trunk/stac/src/stacCombine.c
===================================================================
--- /trunk/stac/src/stacCombine.c	(revision 2660)
+++ /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
 
Index: /trunk/stac/src/stacConfig.c
===================================================================
--- /trunk/stac/src/stacConfig.c	(revision 2660)
+++ /trunk/stac/src/stacConfig.c	(revision 2661)
@@ -9,5 +9,5 @@
 {
     fprintf (stderr, "STAC: Simultaneous Telescope Array Combination\n"
-	     "Usage: %s [-h] [-v] [-g GAIN] [-r RN] [-o NX NY] [-s SAT] [-b BAD] [-k REJ] [-k FRAC] OUT IN1 IN2...\n"
+	     "Usage: %s [-h] [-v] [-g GAIN] [-r RN] [-o NX NY] [-s SAT] [-b BAD] [-k REJ] [-k FRAC] [-G GRAD] OUT IN1 IN2...\n"
 	     "where\n"
 	     "\t-h           Help (this info)\n"
@@ -19,5 +19,6 @@
 	     "\t-b BAD       Bad level (0)\n"
 	     "\t-k REJ       Rejection level (k-sigma; 3.5)\n"
-	     "\t-f FRAC      Fraction of pixel to be marked before considered bad (0.8)\n"
+	     "\t-f FRAC      Fraction of pixel to be marked before considered bad (0.5)\n"
+	     "\t-d GRAD      Gradient threshold for pixel to be marked (0.0)\n"
 	     "\tOUT          Output image\n"
 	     "\tIN1, IN2...  Input images, which have associated .map files.\n",
@@ -37,10 +38,11 @@
     config->inputs = NULL;
     config->output = NULL;
-    config->outnx = 512;
-    config->outny = 512;
+    config->outnx = 1024;
+    config->outny = 1024;
     config->saturated = 65536.0;
     config->bad = 0.0;
-    config->reject = 2.75;
+    config->reject = 3.5;
     config->frac = 0.5;
+    config->grad = 0.4;
 
     return config;
@@ -72,5 +74,5 @@
 
     /* Parse command-line arguments using getopt */
-    while ((opt = getopt(argc, argv, "hvg:r:o:s:b:k:f:")) != -1) {
+    while ((opt = getopt(argc, argv, "hvg:r:o:s:b:k:f:G:")) != -1) {
         switch (opt) {
 	  case 'h':
@@ -121,4 +123,9 @@
 	    }
 	    break;
+	  case 'G':
+	    if (sscanf(optarg, "%f", &config->grad) != 1) {
+		help(programName);
+	    }
+	    break;
 	  default:
 	    help(programName);
Index: /trunk/stac/src/stacRejection.c
===================================================================
--- /trunk/stac/src/stacRejection.c	(revision 2660)
+++ /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);
Index: /trunk/stac/src/stacTransform.c
===================================================================
--- /trunk/stac/src/stacTransform.c	(revision 2660)
+++ /trunk/stac/src/stacTransform.c	(revision 2661)
@@ -20,5 +20,5 @@
 						   unsigned int maskVal,
 						   psF64 unexposedValue)
-{ 
+{
     double floorX = floor((psF64)(x) - 0.5); 
     double floorY = floor((psF64)(y) - 0.5); 
@@ -206,10 +206,11 @@
 		(void)psPlaneTransformApply(detector, map, sky);
 
-		outImage->data.F32[y][x] = (psF32)psImagePixelInterpolate(image, detector->x, detector->y,
-									  mask, 1, 0.0,
+		outImage->data.F32[y][x] = (psF32)psImagePixelInterpolate(image, detector->x + 0.5,
+									  detector->y + 0.5, mask, 1, 0.0,
 									  PS_INTERPOLATE_BILINEAR);
-		outError->data.F32[y][x] = (psF32)p_psImageErrorInterpolateBILINEAR_F32(error, detector->x,
-											detector->y, mask, 1,
-											0.0);
+		outError->data.F32[y][x] = (psF32)p_psImageErrorInterpolateBILINEAR_F32(error,
+											detector->x + 0.5,
+											detector->y + 0.5,
+											mask, 1, 0.0);
 		outError->data.F32[y][x] = sqrtf(outError->data.F32[y][x]);
 	    }