Index: trunk/stac/src/Makefile
===================================================================
--- trunk/stac/src/Makefile	(revision 3828)
+++ trunk/stac/src/Makefile	(revision 5743)
@@ -1,15 +1,26 @@
 SHELL = /bin/sh
 CC = gcc
-CFLAGS += -g -std=c99 -I/home/mithrandir/price/psLib/include -D_GNU_SOURCE # -DTESTING -DCRFLUX
-PSLIB += -L/home/mithrandir/price/psLib/lib/ -lpslib -lgsl -lgslcblas -lfftw3f -lsla -lcfitsio -lm -lxml2 -lmysqlclient
-LDFLAGS += $(PSLIB)
 
+### psLib
+PSLIB_INCLUDE := $(shell pslib-config --cflags)
+PSLIB_LINK := $(shell pslib-config --libs)
+
+### Additional flags for diagnostics, testing, etc
+#PAP_CFLAGS = -DTESTING
+#PAP_CFLAGS = -DCRFLUX
+
+### Build flags
+CFLAGS += -g -O2 -std=c99 -Werror -D_GNU_SOURCE -DPS_NO_TRACE $(PSLIB_INCLUDE) $(PAP_CFLAGS)
+LDFLAGS += $(PSLIB_LINK)
+
+### Ingredients for each program
 STAC = stac.o stacConfig.o stacRead.o stacErrorImages.o stacTransform.o stacCheckMemory.o stacHelp.o \
-	stacCombine.o stacInvertMaps.o stacRejection.o stacAreaOfInterest.o stacSize.o stacScales.o time.o
+	stacCombine.o stacInvertMaps.o stacRejection.o stacAreaOfInterest.o stacSize.o stacScales.o \
+	stacTime.o
 
 SHIFT = shift.o stacRead.o stacTransform.o stacCheckMemory.o stacInvertMaps.o stacSize.o
 
 COMBINE = combine.o combineConfig.o stacRead.o stacErrorImages.o stacCombine.o stacScales.o \
-	stacCheckMemory.o time.o
+	stacCheckMemory.o stacTime.o
 
 SHIFTSIZE = shiftSize.o stacWrite.o stacConfig.o stacRead.o stacCheckMemory.o stacSize.o
@@ -17,10 +28,19 @@
 FIXIMAGE = fixImage.o
 
-TARGETS = stac shift combine shiftSize fixImage
+CALCGRADIENT = calcGradient.o stacRejection.o
 
+SUM = sum.o
+
+### List of targets
+TARGETS = stac shift combine shiftSize calcGradient sum
+
+
+### Build recipes
 .PHONY: tags clean empty test profile optimise all
 
 .c.o:
 		$(CC) -c $(CFLAGS) $(OPTFLAGS) $<
+
+all:		$(TARGETS)
 
 stac:		$(STAC)
@@ -36,8 +56,10 @@
 		$(CC) $(CFLAGS) -o $@ $(SHIFTSIZE) $(LDFLAGS) $(OPTFLAGS)
 
-fixImage:	$(FIXIMAGE)
-		$(CC) $(CFLAGS) -o $@ $(FIXIMAGE) $(LDFLAGS) $(OPTFLAGS)
+calcGradient:	$(CALCGRADIENT)
+		$(CC) $(CFLAGS) -o $@ $(CALCGRADIENT) $(LDFLAGS) $(OPTFLAGS)
 
-all:		$(TARGETS)
+sum:		$(SUM)
+		$(CC) $(CFLAGS) -o $@ $(SUM) $(LDFLAGS) $(OPTFLAGS)
+
 
 clean:
Index: trunk/stac/src/calcGradient.c
===================================================================
--- trunk/stac/src/calcGradient.c	(revision 5743)
+++ trunk/stac/src/calcGradient.c	(revision 5743)
@@ -0,0 +1,55 @@
+#include <stdio.h>
+#include "pslib.h"
+#include "stac.h"
+
+int main(int argc, char *argv[])
+{
+    if (argc != 3) {
+        printf("Usage: %s FITS_IN FITS_OUT\n", argv[0]);
+        exit(EXIT_FAILURE);
+    }
+    const char *inName = argv[1];       // Input filename
+    const char *outName = argv[2];      // Output filename
+
+    psFits *inFile = psFitsOpen(inName, "r"); // File to read
+    psRegion imageRegion = {0, 0, 0, 0}; // Region of image to read
+
+    // We only read PHUs --- not mucking around with extensions for now
+    psImage *image = psFitsReadImage(NULL, inFile, imageRegion, 0);
+    if (image == NULL) {
+        psErrorStackPrint(stderr,"Fatal error: Unable to read %s\n", inName);
+        exit(EXIT_FAILURE);
+    }
+    psTrace("calcGradients.read", 4, "Image %s is %dx%d\n", inName, image->numCols, image->numRows);
+    // Convert to 32-bit floating point, in necessary
+    if (image->type.type != PS_TYPE_F32) {
+        psTrace("calcGradients.read", 3, "Converting %s to floating point in memory....", inName);
+        psImage *temp = psImageCopy(NULL, image, PS_TYPE_F32);
+        psFree(image);
+        image = temp;
+    } else {
+        int numNaN = psImageClipNaN(image, FLT_MAX);
+        if (numNaN) {
+            psTrace("calcGradients.read", 5, "Clipped %d NaN pixels.\n", numNaN);
+        }
+    }
+    psFitsClose(inFile);
+
+    psImage *grad = psImageAlloc(image->numCols, image->numRows, PS_TYPE_F32);
+    for (int y = 0; y < image->numRows; y++) {
+        for (int x = 0; x < image->numCols; x++) {
+            grad->data.F32[y][x] = stacGradient(image, x, y);
+        }
+    }
+
+    psFree(image);
+
+    psFits *outFile = psFitsOpen(outName, "w");
+    if (!psFitsWriteImage(outFile, NULL, grad, 0)) {
+        psErrorStackPrint(stderr, "Unable to write image: %s\n", outName);
+    }
+    psFitsClose(outFile);
+
+    return EXIT_SUCCESS;
+}
+
Index: trunk/stac/src/combine.c
===================================================================
--- trunk/stac/src/combine.c	(revision 3828)
+++ trunk/stac/src/combine.c	(revision 5743)
@@ -21,21 +21,22 @@
     combineConfig *config = combineParseConfig(argc, argv); // Configuration
 
-    psArray *images = stacReadImages(config->inNames); // The images
+    psArray *headers = NULL;            // Array of headers
+    psArray *images = stacReadImages(&headers, config->inNames); // The images
     // Make fake errors
     psArray *errors = psArrayAlloc(images->n);
     for (int i = 0; i < images->n; i++) {
-	psImage *image = images->data[i];
-	psImage *err = psImageAlloc(image->numCols, image->numRows, PS_TYPE_F32);
-	for (int y = 0; y < image->numRows; y++) {
-	    for (int x = 0; x < image->numCols; x++) {
-		err->data.F32[y][x] = 1.0;
-	    }
-	}
-	errors->data[i] = err;
+        psImage *image = images->data[i];
+        psImage *err = psImageAlloc(image->numCols, image->numRows, PS_TYPE_F32);
+        for (int y = 0; y < image->numRows; y++) {
+            for (int x = 0; x < image->numCols; x++) {
+                err->data.F32[y][x] = 1.0;
+            }
+        }
+        errors->data[i] = err;
     }
 
     // Calculate scales between images
-    psVector *scales = NULL;		// Relative scales between images
-    psVector *offsets = NULL;		// Offsets between images
+    psVector *scales = NULL;            // Relative scales between images
+    psVector *offsets = NULL;           // Offsets between images
     (void)stacScales(&scales, &offsets, images, config->starFile, config->starMap, 0.0, 0.0, config->aper);
 
@@ -44,6 +45,6 @@
     psVector *bad = psVectorAlloc(images->n, PS_TYPE_F32); // Bad limits
     for (int i = 0; i < images->n; i++) {
-	saturated->data.F32[i] = (config->saturated - offsets->data.F32[i]) / scales->data.F32[i];
-	bad->data.F32[i] = (config->bad - offsets->data.F32[i]) / scales->data.F32[i];
+        saturated->data.F32[i] = (config->saturated - offsets->data.F32[i]) / scales->data.F32[i];
+        bad->data.F32[i] = (config->bad - offsets->data.F32[i]) / scales->data.F32[i];
     }
 
@@ -52,16 +53,17 @@
 
     // Combine the images
-    psImage *combined = NULL;		// Combined image
-    psArray *rejected = NULL;		// Array of rejection masks
+    psImage *combined = NULL;           // Combined image
+    psArray *rejected = NULL;           // Array of rejection masks
     stacCombine(&combined, &rejected, images, errors, config->nReject, NULL, saturated, bad, config->reject);
 
-    psFits *outFile = psFitsAlloc(config->outName);
-    if (!psFitsWriteImage(outFile, NULL, combined, 0, NULL)) {
-	psErrorStackPrint(stderr, "Unable to write image: %s\n", config->outName);
+    psFits *outFile = psFitsOpen(config->outName, "r");
+    if (!psFitsWriteImage(outFile, headers->data[0], combined, 0)) {
+        psErrorStackPrint(stderr, "Unable to write image: %s\n", config->outName);
     }
     psTrace("combine", 1, "Combined image written to %s\n", config->outName);
-    psFree(outFile);
+    psFitsClose(outFile);
 
     // Clean up
+    psFree(headers);
     psFree(combined);
     psFree(rejected);
Index: trunk/stac/src/fixImage.c
===================================================================
--- trunk/stac/src/fixImage.c	(revision 3828)
+++ 	(revision )
@@ -1,56 +1,0 @@
-#include <stdio.h>
-#include "pslib.h"
-
-void help(const char *programName
-    )
-{
-    fprintf (stderr, "fixImage: Fix an image (floating point with no NaNs)\n"
-	     "Usage: %s IN OUT\n"
-	     "where\n"
-	     "\tIN           Input image\n"
-	     "\tOUT          Output image\n",
-	     programName
-	);
-}
-
-
-int main(int argc, char **argv)
-{
-    if (argc != 3) {
-	help(argv[0]);
-	exit(EXIT_FAILURE);
-    }
-    const char *inName = argv[1];	// Input filename
-    const char *outName = argv[2];	// Output filename
-
-    psFits *imageFile = psFitsAlloc(inName);
-    psRegion *imageRegion = psRegionAlloc(0, 0, 0, 0); // Region of image to read
-    psImage *image = psFitsReadImage(NULL, imageFile, *imageRegion, 0);
-    if (image == NULL) {
-	psErrorStackPrint(stderr, "Fatal error: Unable to read %s\n", inName);
-	exit(EXIT_FAILURE);
-    }
-    psTrace("shift", 4, "Image %s is %dx%d\n", inName, image->numCols, image->numRows);
-    // Convert to 32-bit floating point, in necessary
-    if (image->type.type != PS_TYPE_F32) {
-	psTrace("shift", 5, "Converting %s to floating point in memory....", inName);
-	psImage *temp = psImageCopy(NULL, image, PS_TYPE_F32);
-	psFree(image);
-	image = temp;
-    }
-    psFree(imageFile);
-    psFree(imageRegion);
-
-    int numPix = psImageClipNaN(image, 0.0);
-    printf("%d NaN pixels clipped.\n", numPix);
-
-    // Write out transformed image
-    psFits *outFile = psFitsAlloc(outName);
-    if (!psFitsWriteImage(outFile, NULL, image, 0, NULL)) {
-	psErrorStackPrint(stderr, "Unable to write image: %s\n", outName);
-    }
-    psTrace("shift", 1, "Fixed image written to %s\n", outName);
-    psFree(outFile);
-    psFree(image);
-
-}
Index: trunk/stac/src/shift.c
===================================================================
--- trunk/stac/src/shift.c	(revision 3828)
+++ trunk/stac/src/shift.c	(revision 5743)
@@ -11,13 +11,13 @@
 {
     fprintf (stderr, "shift: shift an image, given the transformation\n"
-	     "Usage: %s [-h] [-v] [-s NX NY] IN OUT\n"
-	     "where\n"
-	     "\t-h           Help (this info)\n"
-	     "\t-v           Increase verbosity level\n"
-	     "\t-s NX NY     Size of output image\n"
-	     "\tIN           Input image, which has associated .map file\n"
-	     "\tOUT          Output image\n",
-	     programName
-	);
+             "Usage: %s [-h] [-v] [-s NX NY] IN OUT\n"
+             "where\n"
+             "\t-h           Help (this info)\n"
+             "\t-v           Increase verbosity level\n"
+             "\t-s NX NY     Size of output image\n"
+             "\tIN           Input image, which has associated .map file\n"
+             "\tOUT          Output image\n",
+             programName
+        );
 }
 
@@ -38,9 +38,10 @@
 
     // Parameters with default values
-    int outnx = 0, outny = 0;		// Size of output image
-    int verbose = 0;			// Verbosity level
+    int outnx = 0, outny = 0;           // Size of output image
+    int verbose = 0;                    // Verbosity level
+    float bad = 0.0;                    // Bad level
 
     // Parse command line
-    const char *programName = argv[0];	// Program name
+    const char *programName = argv[0];  // Program name
     /* Variables for getopt */
     int opt;   /* Option, from getopt */
@@ -49,13 +50,13 @@
 
     /* Parse command-line arguments using getopt */
-    while ((opt = getopt(argc, argv, "hvs:")) != -1) {
+    while ((opt = getopt(argc, argv, "hvb:s:")) != -1) {
         switch (opt) {
-	  case 'h':
-	    help(programName);
-	    exit(EXIT_SUCCESS);
-	  case 'v':
+          case 'h':
+            help(programName);
+            exit(EXIT_SUCCESS);
+          case 'v':
             verbose++;
             break;
-	  case 's':
+          case 's':
             if ((sscanf(argv[optind-1], "%d", &outnx) != 1) || (sscanf(argv[optind++], "%d", &outny) != 1)) {
                 // Note: incrementing optind, so I can read more than one parameter.
@@ -63,9 +64,15 @@
                 exit(EXIT_FAILURE);
             }
-	    break;
-	  default:
-	    help(programName);
-	    exit(EXIT_FAILURE);
-	}
+            break;
+          case 'b':
+            if (sscanf(optarg, "%f", &bad) != 1) {
+                help(programName);
+                exit(EXIT_FAILURE);
+            }
+            break;
+          default:
+            help(programName);
+            exit(EXIT_FAILURE);
+        }
     }
 
@@ -77,39 +84,53 @@
         exit(EXIT_FAILURE);
     }
-    const char *inName = argv[0];	// Input filename
-    const char *outName = argv[1];	// Output filename
+    const char *inName = argv[0];       // Input filename
+    const char *outName = argv[1];      // Output filename
 
-    psFits *imageFile = psFitsAlloc(inName);
-    psRegion *imageRegion = psRegionAlloc(0, 0, 0, 0); // Region of image to read
-    psImage *image = psFitsReadImage(NULL, imageFile, *imageRegion, 0);
+    psFits *imageFile = psFitsOpen(inName, "r");
+    psRegion imageRegion = {0, 0, 0, 0}; // Region of image to read
+    psMetadata *header = psFitsReadHeader(NULL, imageFile); // FITS header
+    psImage *image = psFitsReadImage(NULL, imageFile, imageRegion, 0);
     if (image == NULL) {
-	psErrorStackPrint(stderr, "Fatal error: Unable to read %s\n", inName);
-	exit(EXIT_FAILURE);
+        psErrorStackPrint(stderr, "Fatal error: Unable to read %s\n", inName);
+        exit(EXIT_FAILURE);
     }
     psTrace("shift", 4, "Image %s is %dx%d\n", inName, image->numCols, image->numRows);
     // Convert to 32-bit floating point, in necessary
     if (image->type.type != PS_TYPE_F32) {
-	psTrace("shift", 5, "Converting %s to floating point in memory....", inName);
-	psImage *temp = psImageCopy(NULL, image, PS_TYPE_F32);
-	psFree(image);
-	image = temp;
+        psTrace("shift", 5, "Converting %s to floating point in memory....", inName);
+        psImage *temp = psImageCopy(NULL, image, PS_TYPE_F32);
+        psFree(image);
+        image = temp;
     }
-    psFree(imageFile);
-    psFree(imageRegion);
+    psFitsClose(imageFile);
+
+    // Generate masks
+    psImage *mask = psImageAlloc(image->numCols, image->numRows, PS_TYPE_U8);
+    for (int y = 0; y < image->numRows; y++) {
+        for (int x = 0; x < image->numCols; x++) {
+            if (image->data.F32[y][x] <= bad) {
+                mask->data.U8[y][x] = 1;
+            } else {
+                mask->data.U8[y][x] = 0;
+            }
+        }
+    }
 
     // Read map
-    char mapName[MAXCHAR];		// Name of map file
+    char mapName[MAXCHAR];              // Name of map file
     sprintf(mapName, "%s.map", inName);
     psPlaneTransform *map = stacReadMap(mapName);
 
     // Functions work on array of images, so:
-    psArray *images = psArrayAlloc(1); // An array of one
-    psArray *maps = psArrayAlloc(1);// An array of one
+    psArray *images = psArrayAlloc(1);  // An array of one
+    psArray *maps = psArrayAlloc(1);    // An array of one
+    psArray *masks = psArrayAlloc(1);   // An array of one
     images->data[0] = image;
     maps->data[0] = map;
+    masks->data[0] = mask;
 
     // Get size
     if (outnx == 0 || outny == 0) {
-	stacSize(&outnx, &outny, NULL, NULL, images, maps);
+        stacSize(&outnx, &outny, NULL, NULL, images, maps);
     }
 
@@ -119,17 +140,29 @@
     // Transform inputs and errors
     psArray *transformed = NULL;
-    (void)stacTransform(&transformed, NULL, images, inverseMaps, NULL, NULL, NULL, NULL, NULL, outnx, outny);
+    (void)stacTransform(&transformed, NULL, images, inverseMaps, NULL, masks, NULL, NULL, NULL, outnx, outny);
+
+    // Update FITS header appropriately
+    psMetadataAddS32(header, PS_LIST_TAIL, "NAXIS1", PS_META_REPLACE, "Number of pixels in x",
+                     (int)((psImage*)transformed->data[0])->numCols);
+    psMetadataAddS32(header, PS_LIST_TAIL, "NAXIS2", PS_META_REPLACE, "Number of pixels in y",
+                     (int)((psImage*)transformed->data[0])->numRows);
+    psMetadataAddS32(header, PS_LIST_TAIL, "BITPIX", PS_META_REPLACE, "Bits per pixel", -32);
 
     // Write out transformed image
-    psFits *outFile = psFitsAlloc(outName);
-    if (!psFitsWriteImage(outFile, NULL, transformed->data[0], 0, NULL)) {
-	psErrorStackPrint(stderr, "Unable to write image: %s\n", outName);
+    psFits *outFile = psFitsOpen(outName, "w");
+    int numPix = psImageClipNaN(transformed->data[0], 0.0);
+    if (numPix > 0) {
+        psTrace("stac", 3, "Clipping %d NaN pixels to zero.\n", numPix);
+    }
+    if (!psFitsWriteImage(outFile, header, transformed->data[0], 0)) {
+        psErrorStackPrint(stderr, "Unable to write image: %s\n", outName);
     }
     psTrace("shift", 1, "Transformed image written to %s\n", outName);
-    psFree(outFile);
+    psFitsClose(outFile);
 
     // Free everything I've used
     psFree(images);
     psFree(maps);
+    psFree(masks);
     psFree(inverseMaps);
     psFree(transformed);
Index: trunk/stac/src/shiftSize.c
===================================================================
--- trunk/stac/src/shiftSize.c	(revision 3828)
+++ trunk/stac/src/shiftSize.c	(revision 5743)
@@ -79,5 +79,5 @@
 
     // Read input files
-    psArray *images = stacReadImages(inputs);
+    psArray *images = stacReadImages(NULL, inputs);
 
     // Read maps
Index: trunk/stac/src/stac.c
===================================================================
--- trunk/stac/src/stac.c	(revision 3828)
+++ trunk/stac/src/stac.c	(revision 5743)
@@ -18,4 +18,5 @@
     // Set trace levels
     psTraceSetLevel(".",0);
+    psTraceSetLevel("stac",10);
     psTraceSetLevel("stac.checkMemory",10);
     psTraceSetLevel("stac.config",10);
@@ -38,5 +39,23 @@
 
     // Read input files
-    psArray *inputs = stacReadImages(config->inputs);
+    psArray *headers = NULL;            // Array of headers
+    psArray *inputs = stacReadImages(&headers, config->inputs);
+
+    // Generate masks
+    psArray *masks = psArrayAlloc(inputs->n);
+    for (int i = 0; i < inputs->n; i++) {
+        psImage *image = inputs->data[i]; // Image for which to get mask
+        psImage *mask = psImageAlloc(image->numCols, image->numRows, PS_TYPE_U8);
+        for (int y = 0; y < image->numRows; y++) {
+            for (int x = 0; x < image->numCols; x++) {
+                if (image->data.F32[y][x] <= config->bad) {
+                    mask->data.U8[y][x] = 1;
+                } else {
+                    mask->data.U8[y][x] = 0;
+                }
+            }
+        }
+        masks->data[i] = mask;
+    }
 
     // Read maps
@@ -47,5 +66,22 @@
     // Get size, if not input
     if (config->outnx == 0 || config->outny == 0) {
-	stacSize(&config->outnx, &config->outny, &config->xMapDiff, &config->yMapDiff, inputs, maps);
+        stacSize(&config->outnx, &config->outny, &config->xMapDiff, &config->yMapDiff, inputs, maps);
+    }
+
+    // Fix up the headers
+    for (int i = 0; i < headers->n; i++) {
+        psMetadata *header = headers->data[i];
+        psMetadataAddS32(header, PS_LIST_TAIL, "NAXIS1", PS_META_REPLACE, "Number of pixels in x",
+                         config->outnx);
+        psMetadataAddS32(header, PS_LIST_TAIL, "NAXIS2", PS_META_REPLACE, "Number of pixels in y",
+                         (int)config->outny);
+        psMetadataAddS32(header, PS_LIST_TAIL, "BITPIX", PS_META_REPLACE, "Bits per pixel", -32);
+#if 0
+        bool mdok = true;               // Result of MD lookup
+        float crpix1 = psMetadataLookupF32(&mdok, header, "CRPIX1");
+        if (mdok && !isnan(crpix1)) {
+            psMetadataAddF32(header, PS_LIST_TAIL, "CRPIX1", PS_META_REPLACE,
+                             "WCS Coordinate reference pixel", crpix1 -
+#endif
     }
 
@@ -58,12 +94,12 @@
     // Write error images out to check
     for (int i = 0; i < inputs->n; i++) {
-	char errName[MAXCHAR];		// Filename of error image
-	sprintf(errName,"%s.err",config->inputs->data[i]);
-	psFits *errorFile = psFitsAlloc(errName);
-	if (!psFitsWriteImage(errorFile, NULL, errors->data[i], 0, NULL)) {
-	    psErrorStackPrint(stderr, "Unable to write image: %s\n", errName);
-	}
-	psTrace("stac", 1, "Error image written to %s\n", errName);
-	psFree(errorFile);
+        char errName[MAXCHAR];          // Filename of error image
+        sprintf(errName,"%s.err",config->inputs->data[i]);
+        psFits *errorFile = psFitsOpen(errName, "w");
+        if (!psFitsWriteImage(errorFile, NULL, errors->data[i], 0)) {
+            psErrorStackPrint(stderr, "Unable to write image: %s\n", errName);
+        }
+        psTrace("stac", 1, "Error image written to %s\n", errName);
+        psFitsClose(errorFile);
     }
 #endif
@@ -72,6 +108,6 @@
     psArray *transformed = NULL;
     psArray *transformedErrors = NULL;
-    (void)stacTransform(&transformed, &transformedErrors, inputs, inverseMaps, errors, NULL, NULL, NULL, NULL,
-			config->outnx, config->outny);
+    (void)stacTransform(&transformed, &transformedErrors, inputs, inverseMaps, errors, masks, NULL, NULL, NULL,
+                        config->outnx, config->outny);
     psTrace("stac.time",1,"Transformation completed at %f seconds\n", getTime() - startTime);
 
@@ -79,39 +115,55 @@
     // Write transformed images out to check
     for (int i = 0; i < inputs->n; i++) {
-	char shiftName[MAXCHAR];	// Filename of shift image
-	char errName[MAXCHAR];		// Filename of error image
-	sprintf(shiftName,"%s.shift1",config->inputs->data[i]);
-	sprintf(errName,"%s.shifterr1",config->inputs->data[i]);
-	psFits *shiftFile = psFitsAlloc(shiftName);
-	psFits *errFile = psFitsAlloc(errName);
-	psImage *trans = transformed->data[i]; // Transformed image
-	psImage *transErr = transformedErrors->data[i];	// Transformed error image
-	if (!psFitsWriteImage(shiftFile, NULL, trans, 0, NULL)) {
-	    psErrorStackPrint(stderr, "Unable to write image: %s\n", shiftName);
-	}
-	psTrace("stac", 1, "Shifted image written to %s\n", shiftName);
-	if (!psFitsWriteImage(errFile, NULL, transErr, 0, NULL)) {
-	    psErrorStackPrint(stderr, "Unable to write image: %s\n", errName);
-	}
-	psTrace("stac", 1, "Shifted error image written to %s\n", errName);
-	psFree(shiftFile);
-	psFree(errFile);
-    }
-#endif
-
-
-
-    psVector *scales = NULL;		// Relative scales between images
-    psVector *offsets = NULL;		// Offsets between images
+        char shiftName[MAXCHAR];        // Filename of shift image
+        char errName[MAXCHAR];          // Filename of error image
+        sprintf(shiftName,"%s.shift1",config->inputs->data[i]);
+        sprintf(errName,"%s.shifterr1",config->inputs->data[i]);
+        psFits *shiftFile = psFitsOpen(shiftName, "w");
+        psFits *errFile = psFitsOpen(errName, "w");
+        psImage *trans = transformed->data[i]; // Transformed image
+        psImage *transErr = transformedErrors->data[i]; // Transformed error image
+        if (!psFitsWriteImage(shiftFile, NULL, trans, 0)) {
+            psErrorStackPrint(stderr, "Unable to write image: %s\n", shiftName);
+        }
+        psTrace("stac", 1, "Shifted image written to %s\n", shiftName);
+        if (!psFitsWriteImage(errFile, NULL, transErr, 0)) {
+            psErrorStackPrint(stderr, "Unable to write image: %s\n", errName);
+        }
+        psTrace("stac", 1, "Shifted error image written to %s\n", errName);
+        psFitsClose(shiftFile);
+        psFitsClose(errFile);
+    }
+#endif
+
+    psVector *scales = NULL;            // Relative scales between images
+    psVector *offsets = NULL;           // Offsets between images
     (void)stacScales(&scales, &offsets, transformed, config->starFile, config->starMapFile, config->xMapDiff,
-		     config->yMapDiff, config->aper);
+                     config->yMapDiff, config->aper);
     // Rescale the images
     (void)stacRescale(transformed, transformedErrors, NULL, scales, offsets);
     // Set the saturation and bad values
     psVector *saturated = psVectorAlloc(transformed->n, PS_TYPE_F32); // Saturation limits
-    psVector *bad = psVectorAlloc(transformed->n, PS_TYPE_F32);	// Bad limits
+    psVector *bad = psVectorAlloc(transformed->n, PS_TYPE_F32); // Bad limits
     for (int i = 0; i < transformed->n; i++) {
-	saturated->data.F32[i] = (config->saturated - offsets->data.F32[i]) / scales->data.F32[i];
-	bad->data.F32[i] = (config->bad - offsets->data.F32[i]) / scales->data.F32[i];
+        saturated->data.F32[i] = (config->saturated - offsets->data.F32[i]) / scales->data.F32[i];
+        bad->data.F32[i] = (config->bad - offsets->data.F32[i]) / scales->data.F32[i];
+    }
+
+    // Save shifted images
+    if (config->saveShifts) {
+        psImage *image = NULL;          // Copy of image, to remove NAN for output
+        for (int i = 0; i < transformed->n; i++) {
+            char shiftName[MAXCHAR];    // Filename of shifted image
+            sprintf(shiftName, "%s.shift", config->inputs->data[i]);
+            psFits *shiftFile = psFitsOpen(shiftName, "w");
+            image = psImageCopy(NULL, transformed->data[i], PS_TYPE_F32);
+            (void)psImageClipNaN(image, 0.0);
+            if (!psFitsWriteImage(shiftFile, headers->data[i], image, 0)) {
+                psErrorStackPrint(stderr, "Unable to write image: %s\n", shiftName);
+            }
+            psTrace("stac", 1, "Shifted image %d written to %s\n", i, shiftName);
+            psFitsClose(shiftFile);
+            psFree(image);
+        }
     }
 
@@ -120,5 +172,5 @@
     psImage *combined = NULL;
     (void)stacCombine(&combined, &rejected, transformed, transformedErrors, config->nReject, NULL, saturated,
-		      bad, config->reject);
+                      bad, config->reject);
 
     psTrace("stac.time",1,"First combination completed at %f seconds\n", getTime() - startTime);
@@ -128,24 +180,24 @@
     // Write rejection images out to check
     for (int i = 0; i < rejected->n; i++) {
-	char rejName[MAXCHAR];	// Filename of rejection image
-	sprintf(rejName, "%s.shiftrej", config->inputs->data[i]);
-	
-	psFits *rejFile = psFitsAlloc(rejName);
-	if (!psFitsWriteImage(rejFile, NULL, rejected->data[i], 0, NULL)) {
-	    psErrorStackPrint(stderr, "Unable to write image: %s\n", rejName);
-	}
-	psTrace("stac", 1, "Rejection image written to %s\n", rejName);
-	psFree(rejFile);
+        char rejName[MAXCHAR];  // Filename of rejection image
+        sprintf(rejName, "%s.shiftrej", config->inputs->data[i]);
+
+        psFits *rejFile = psFitsOpen(rejName, "w");
+        if (!psFitsWriteImage(rejFile, NULL, rejected->data[i], 0)) {
+            psErrorStackPrint(stderr, "Unable to write image: %s\n", rejName);
+        }
+        psTrace("stac", 1, "Rejection image written to %s\n", rejName);
+        psFitsClose(rejFile);
     }
 
     // Write out pre-combined image
-    char preName[MAXCHAR];		// Filename of precombined image
+    char preName[MAXCHAR];              // Filename of precombined image
     sprintf(preName, "%s.pre", config->output);
-    psFits *preFile = psFitsAlloc(preName);
-    if (!psFitsWriteImage(preFile, NULL, combined, 0, NULL)) {
-	psErrorStackPrint(stderr, "Unable to write image: %s\n", preName);
+    psFits *preFile = psFitsOpen(preName, "w");
+    if (!psFitsWriteImage(preFile, NULL, combined, 0)) {
+        psErrorStackPrint(stderr, "Unable to write image: %s\n", preName);
     }
     psTrace("stac", 1, "Pre-combined image written to %s\n", preName);
-    psFree(preFile);
+    psFitsClose(preFile);
 #endif
 
@@ -153,16 +205,16 @@
     psArray *regions = psArrayAlloc(inputs->n); // Array of images denoting regions of interest
     for (int i = 0; i < inputs->n; i++) {
-	regions->data[i] = stacAreaOfInterest(rejected->data[i], inverseMaps->data[i],
-					      ((psImage*)(inputs->data[i]))->numCols,
-					      ((psImage*)(inputs->data[i]))->numRows);
-#ifdef TESTING
-	char regionName[MAXCHAR];	// Filename of region image
-	sprintf(regionName,"%s.region",config->inputs->data[i]);
-	psFits *regionFile = psFitsAlloc(regionName);
-	if (!psFitsWriteImage(regionFile, NULL, regions->data[i], 0, NULL)) {
-	    psErrorStackPrint(stderr, "Unable to write image: %s\n", regionName);
-	}
-	psTrace("stac", 1, "Region image written to %s\n", regionName);
-	psFree(regionFile);
+        regions->data[i] = stacAreaOfInterest(rejected->data[i], inverseMaps->data[i],
+                                              ((psImage*)(inputs->data[i]))->numCols,
+                                              ((psImage*)(inputs->data[i]))->numRows);
+#ifdef TESTING
+        char regionName[MAXCHAR];       // Filename of region image
+        sprintf(regionName,"%s.region",config->inputs->data[i]);
+        psFits *regionFile = psFitsOpen(regionName, "w");
+        if (!psFitsWriteImage(regionFile, NULL, regions->data[i], 0)) {
+            psErrorStackPrint(stderr, "Unable to write image: %s\n", regionName);
+        }
+        psTrace("stac", 1, "Region image written to %s\n", regionName);
+        psFitsClose(regionFile);
 #endif
     }
@@ -170,23 +222,34 @@
     // Transform rejected pixels to source frame
     psArray *rejectedSource = stacRejection(inputs, rejected, regions, maps, inverseMaps, config->frac,
-					    config->grad, config->inputs);
+                                            config->grad, config->inputs);
 
     // Get regions of interest in the output frame
     psImage *combineRegion = psImageAlloc(config->outnx, config->outny, PS_TYPE_U8);
     for (int y = 0; y < config->outny; y++) {
-	for (int x = 0; x < config->outnx; x++) {
-	    combineRegion->data.U8[y][x] = 0;
-	}
-    }
-    for (int i = 0; i < inputs->n; i++) {
-	psImage *region = stacAreaOfInterest(rejectedSource->data[i], maps->data[i], config->outnx,
-					     config->outny);
-	psBinaryOp(combineRegion, combineRegion, "+", region);
-	psFree(region);
+        for (int x = 0; x < config->outnx; x++) {
+            combineRegion->data.U8[y][x] = 0;
+        }
+    }
+    for (int i = 0; i < inputs->n; i++) {
+        psImage *region = stacAreaOfInterest(rejectedSource->data[i], maps->data[i], config->outnx,
+                                             config->outny);
+        psBinaryOp(combineRegion, combineRegion, "+", region);
+        psFree(region);
+
+        // Do OR of masks
+        psImage *mask = masks->data[i]; // Mask of input image (bad columns etc)
+        psImage *cr = rejectedSource->data[i]; // Mask of CRs
+        for (int y = 0; y < mask->numRows; y++) {
+            for (int x = 0; x < mask->numCols; x++) {
+                if (cr->data.U8[y][x] > 0) {
+                    mask->data.U8[y][x] = 1;
+                }
+            }
+        }
     }
 
     // Redo transformation with the masks and scales/offsets
-    (void)stacTransform(&transformed, &transformedErrors, inputs, inverseMaps, errors, rejectedSource,
-			combineRegion, scales, offsets, config->outnx, config->outny);
+    (void)stacTransform(&transformed, &transformedErrors, inputs, inverseMaps, errors, masks,
+                        combineRegion, scales, offsets, config->outnx, config->outny);
 
     // Combine the newly-transformed CR-free images, no rejection
@@ -194,27 +257,17 @@
     rejected = NULL;
     (void)stacCombine(&combined, &rejected, transformed, transformedErrors, 0, combineRegion, saturated,
-		      bad, config->reject);
+                      bad, config->reject);
 
     // Write output image
-    psFits *outFile = psFitsAlloc(config->output);
-    if (!psFitsWriteImage(outFile, NULL, combined, 0, NULL)) {
-	psErrorStackPrint(stderr, "Unable to write image: %s\n", config->output);
+    psFits *outFile = psFitsOpen(config->output, "w");
+    int numPix = psImageClipNaN(combined, 0.0);
+    if (numPix > 0) {
+        psTrace("stac", 3, "Clipping %d NaN pixels to zero.\n", numPix);
+    }
+    if (!psFitsWriteImage(outFile, headers->data[0], combined, 0)) {
+        psErrorStackPrint(stderr, "Unable to write image: %s\n", config->output);
     }
     psTrace("stac", 1, "Combined image written to %s\n", config->output);
-    psFree(outFile);
-
-    // Save shifted images
-    if (config->saveShifts) {
-	for (int i = 0; i < transformed->n; i++) {
-	    char shiftName[MAXCHAR];	// Filename of shifted image
-	    sprintf(shiftName, "%s.shift", config->inputs->data[i]);
-	    psFits *shiftFile = psFitsAlloc(shiftName);
-	    if (!psFitsWriteImage(shiftFile, NULL, transformed->data[i], 0, NULL)) {
-		psErrorStackPrint(stderr, "Unable to write image: %s\n", shiftName);
-	    }
-	    psTrace("stac", 1, "Shifted image %d written to %s\n", i, shiftName);
-	    psFree(shiftFile);
-	}
-    }
+    psFitsClose(outFile);
 
     // Free everything I've used
@@ -228,8 +281,11 @@
     psFree(inverseMaps);
     psFree(errors);
+    psFree(masks);
     psFree(transformedErrors);
     psFree(transformed);
     psFree(rejectedSource);
     psFree(combined);
+    psFree(saturated);
+    psFree(bad);
 
     psTrace("stac.time",1,"Final combination completed at %f seconds\n", getTime() - startTime);
Index: trunk/stac/src/stac.h
===================================================================
--- trunk/stac/src/stac.h	(revision 3828)
+++ trunk/stac/src/stac.h	(revision 5743)
@@ -18,5 +18,6 @@
 
 // Read the input files and return an array of images
-psArray *stacReadImages(psArray *filenames // The file names of the images
+psArray *stacReadImages(psArray **headers, // The image headers, to be returned
+			psArray *filenames // The file names of the images
     );
 
@@ -72,5 +73,5 @@
 
 // Print out a memblock when it's allocated --- this function used as a callback
-psMemoryId stacMemPrint(const psMemBlock *ptr);
+psMemId stacMemPrint(const psMemBlock *ptr);
 
 
Index: trunk/stac/src/stacAreaOfInterest.c
===================================================================
--- trunk/stac/src/stacAreaOfInterest.c	(revision 3828)
+++ trunk/stac/src/stacAreaOfInterest.c	(revision 5743)
@@ -13,6 +13,6 @@
     )
 {
-    psDPolynomial2D *xPoly = transform->x;
-    psDPolynomial2D *yPoly = transform->y;
+    psPolynomial2D *xPoly = transform->x;
+    psPolynomial2D *yPoly = transform->y;
 
     assert(xPoly->type == PS_POLYNOMIAL_ORD);
Index: trunk/stac/src/stacCheckMemory.c
===================================================================
--- trunk/stac/src/stacCheckMemory.c	(revision 3828)
+++ trunk/stac/src/stacCheckMemory.c	(revision 5743)
@@ -7,5 +7,5 @@
 
 
-psMemoryId stacMemPrint(const psMemBlock *ptr)
+psMemId stacMemPrint(const psMemBlock *ptr)
 {
     psLogMsg("stac.memoryPrint", PS_LOG_INFO,
@@ -41,5 +41,5 @@
     psTrace("stac.checkMemory", 1, "Checking for memory problems....\n");
 
-    (void)psMemProblemCallbackSet(stacMemoryProblem); // Set callback for corruption
+    (void)psMemProblemCallbackSet((psMemProblemCallback)stacMemoryProblem); // Set callback for corruption
 
     if ((leakFile = fopen(LEAKS, "w")) == NULL) {
Index: trunk/stac/src/stacCombine.c
===================================================================
--- trunk/stac/src/stacCombine.c	(revision 3828)
+++ trunk/stac/src/stacCombine.c	(revision 5743)
@@ -10,5 +10,5 @@
 float stacCombineMean(psVector *values,	// Values for which to take the mean
 		      psVector *errors,	// Errors in the values
-		      psVector *masks	// Masks for the values, 0 = don't use, 1 = use
+		      psVector *masks	// Masks for the values, 1 = don't use, 0 = use
     )
 {
@@ -149,6 +149,10 @@
 		    }
 		}
-	    
+		
 		float average = stacCombineMean(pixels, deltas, mask); // Combined value
+
+		// We set the value BEFORE the rejection iteration because not all pixels that we reject
+		// here will be rejected in the final cut.
+		(*combined)->data.F32[y][x] = average;
 
 #ifdef TESTING
@@ -188,6 +192,4 @@
 		}
 	    
-		(*combined)->data.F32[y][x] = average;
-
 	    } // Pixels of interest
 
@@ -201,5 +203,5 @@
     sprintf(chiName,"chi2_%d.fits",iteration);
     psFits *chiFile = psFitsAlloc(chiName);
-    if (!psFitsWriteImage(chiFile, NULL, chi2 , 0, NULL)) {
+    if (!psFitsWriteImage(chiFile, NULL, chi2 , 0)) {
 	psErrorStackPrint(stderr, "Unable to write image: %s\n", chiName);
     }
Index: trunk/stac/src/stacConfig.c
===================================================================
--- trunk/stac/src/stacConfig.c	(revision 3828)
+++ trunk/stac/src/stacConfig.c	(revision 5743)
@@ -26,7 +26,7 @@
     config->saturated = 65535.0;	// Saturation level
     config->bad = 0.0;			// Bad level
-    config->reject = 3.0;
-    config->frac = 0.5;
-    config->grad = 0.6;
+    config->reject = 2.5;
+    config->frac = 0.45;
+    config->grad = 0.7;
     config->nReject = 2;
 
Index: trunk/stac/src/stacInvertMaps.c
===================================================================
--- trunk/stac/src/stacInvertMaps.c	(revision 3828)
+++ trunk/stac/src/stacInvertMaps.c	(revision 5743)
@@ -27,10 +27,10 @@
 	assert(oldMap->x->nX == oldMap->x->nY && oldMap->y->nX == oldMap->y->nY &&
 	       oldMap->x->nX == oldMap->y->nX);
-	int order = oldMap->x->nX;	// Polynomial order
+	int order = oldMap->x->nX + 1;	// Polynomial order
 	psTrace("stac.invertMaps", 4, "Generating order %d polynomial inverse transformation.\n", order);
-	psPlaneTransform *newMap = psPlaneTransformAlloc(order, order);	// Inverted map
+	psPlaneTransform *newMap = psPlaneTransformAlloc(order + 1, order + 1);	// Inverted map
 
 	// Create fake polynomial to use in evaluation
-	psDPolynomial2D *fakePoly = psDPolynomial2DAlloc(order, order, PS_POLYNOMIAL_ORD);
+	psPolynomial2D *fakePoly = psPolynomial2DAlloc(order, order, PS_POLYNOMIAL_ORD);
 	for (int i = 0; i < order; i++) {
 	    for (int j = 0; j < order; j++) {
@@ -82,5 +82,5 @@
 
 		    fakePoly->mask[i][j] = 0;
-		    double ijPoly = psDPolynomial2DEval(fakePoly, xIn->data.F32[g], yIn->data.F32[g]);
+		    double ijPoly = psPolynomial2DEval(fakePoly, xIn->data.F32[g], yIn->data.F32[g]);
 		    fakePoly->mask[i][j] = 1;
 
@@ -89,5 +89,5 @@
 			    
 			    fakePoly->mask[m][n] = 0;
-			    double mnPoly = psDPolynomial2DEval(fakePoly, xIn->data.F32[g], yIn->data.F32[g]);
+			    double mnPoly = psPolynomial2DEval(fakePoly, xIn->data.F32[g], yIn->data.F32[g]);
 			    fakePoly->mask[m][n] = 1;
 
Index: trunk/stac/src/stacRead.c
===================================================================
--- trunk/stac/src/stacRead.c	(revision 3828)
+++ trunk/stac/src/stacRead.c	(revision 5743)
@@ -1,45 +1,53 @@
 #include <stdio.h>
 #include <string.h>
+#include <assert.h>
 #include "pslib.h"
 #include "stac.h"
 
-#define BUFFER 100			// Size of buffer for incrementally reading coordinates
-
-psArray *stacReadImages(psArray *filenames // The file names of the images
+#define BUFFER 100                      // Size of buffer for incrementally reading coordinates
+
+psArray *stacReadImages(psArray **headers, // The image headers, to be returned
+    psArray *filenames // The file names of the images
     )
 {
-    int nFiles = filenames->n;		// The number of input files
+    int nFiles = filenames->n;          // The number of input files
     psArray *images = psArrayAlloc(nFiles); // The input files, to be returned
-    psRegion *imageRegion = psRegionAlloc(0, 0, 0, 0); // Region of image to read
+    assert(!headers || ! *headers || (*headers)->n == nFiles);
+    if (headers && ! *headers) {
+        *headers = psArrayAlloc(nFiles);
+    }
+
+    psRegion imageRegion = {0, 0, 0, 0}; // Region of image to read
 
     psTrace("stac.read.images", 1, "Reading input images....\n");
     for (int i = 0; i < nFiles; i++) {
-	psTrace("stac.read.images", 2, "Reading input image %s....\n",filenames->data[i]);
-	psFits *imageFile = psFitsAlloc(filenames->data[i]);
-	// We only read PHUs --- not mucking around with extensions for now
-	psImage *image = psFitsReadImage(NULL, imageFile, *imageRegion, 0);
-	if (image == NULL) {
-	    psErrorStackPrint(stderr,"Fatal error: Unable to read %s\n",filenames->data[i]);
-	    exit(EXIT_FAILURE);
-	}
-	psTrace("stac.read.images", 4, "Image %s is %dx%d\n", filenames->data[i], image->numCols,
-		image->numRows);
-	// Convert to 32-bit floating point, in necessary
-	if (image->type.type != PS_TYPE_F32) {
-	    psTrace("stac.read.images", 3, "Converting %s to floating point in memory....",
-		    filenames->data[i]);
-	    images->data[i] = psImageCopy(NULL, image, PS_TYPE_F32);
-	    psFree(image);
-	} else {
-	    int numNaN = psImageClipNaN(image, FLT_MAX);
-	    if (numNaN) {
-		psTrace("stac.read.images", 5, "Clipped %d NaN pixels.\n", numNaN);
-	    }
-	    images->data[i] = image;
-	}
-	psFree(imageFile);
+        psTrace("stac.read.images", 2, "Reading input image %s....\n",filenames->data[i]);
+        psFits *imageFile = psFitsOpen(filenames->data[i], "r");
+        // We only read PHUs --- not mucking around with extensions for now
+        if (headers) {
+            (*headers)->data[i] = psFitsReadHeader(NULL, imageFile);
+        }
+        psImage *image = psFitsReadImage(NULL, imageFile, imageRegion, 0);
+        if (image == NULL) {
+            psErrorStackPrint(stderr,"Fatal error: Unable to read %s\n",filenames->data[i]);
+            exit(EXIT_FAILURE);
+        }
+        psFitsClose(imageFile);
+        psTrace("stac.read.images", 4, "Image %s is %dx%d\n", filenames->data[i], image->numCols,
+                image->numRows);
+        // Convert to 32-bit floating point, in necessary
+        if (image->type.type != PS_TYPE_F32) {
+            psTrace("stac.read.images", 3, "Converting %s to floating point in memory....\n",
+                    filenames->data[i]);
+            images->data[i] = psImageCopy(NULL, image, PS_TYPE_F32);
+            psFree(image);
+        }
+        int numNaN = psImageClipNaN(image, 0.0);
+        if (numNaN) {
+            psTrace("stac.read.images", 5, "Clipped %d NaN pixels to zero.\n", numNaN);
+        }
+        images->data[i] = image;
     }
     psTrace("stac.read.images",1,"%d input images read.\n",nFiles);
-    psFree(imageRegion);
 
     return images;
@@ -50,6 +58,6 @@
     FILE *file = fopen(filename, "r");
     if (file == NULL) {
-	psLogMsg("stac.read.coords", PS_LOG_ERROR, "Cannot open coordinate file, %s\n", filename);
-	return NULL;
+        psLogMsg("stac.read.coords", PS_LOG_ERROR, "Cannot open coordinate file, %s\n", filename);
+        return NULL;
     }
 
@@ -57,15 +65,15 @@
 
     psArray *coords = psArrayAlloc(BUFFER); // The array of coordinates to be returned
-    float x, y;				// Coordinates to read
-    int num = 0;			// Number of coordinates read
+    float x, y;                         // Coordinates to read
+    int num = 0;                        // Number of coordinates read
     while (fscanf(file, "%f %f\n", &x, &y) == 2) {
-	psPlane *coord = psPlaneAlloc();// A coordinate
-	coord->x = x;
-	coord->y = y;
-	coords->data[num] = coord;
-	num++;
-	if (num % BUFFER) {
-	    coords = psArrayRealloc(coords, num + BUFFER);
-	}
+        psPlane *coord = psPlaneAlloc();// A coordinate
+        coord->x = x;
+        coord->y = y;
+        coords->data[num] = coord;
+        num++;
+        if (num % BUFFER) {
+            coords = psArrayRealloc(coords, num + BUFFER);
+        }
     }
     coords->n = num;
@@ -82,86 +90,86 @@
     FILE *mapfp = fopen(filename, "r");
     if (mapfp == NULL) {
-	psLogMsg("stac.read.map", PS_LOG_ERROR, "Cannot open map file, %s\n", filename);
-	return NULL;
+        psLogMsg("stac.read.map", PS_LOG_ERROR, "Cannot open map file, %s\n", filename);
+        return NULL;
     }
     // Read the file
     psTrace("stac.read.map", 5, "Reading map file %s....\n", filename);
-    
+
     // Format is now:
     // order
     // x coefficients
     // y coefficients
-    // 
+    //
     // where the order is 1 for linear, 2 for quadratic, 3 for cubic.
     // and the coefficients are read by the following pseudo-code:
-    // 
-    // 	for (int k = 0; k < order + 1; k++)
-    // 	    for (int j = 0; j < k; j++)
-    // 		int i = k - j;
+    //
+    //  for (int k = 0; k < order + 1; k++)
+    //      for (int j = 0; j < k; j++)
+    //          int i = k - j;
     //              read coefficient of x^i y^j
-    // 
+    //
     // This produces the ordering:
     // x^0 y^0, x^1 y^0, x^0 y^1, x^2 y^0, x^1 y^1, x^0 y^2, x^3 y^0, x^2 y^1, x^1 y^2, x^0 y^3
-    // 
+    //
     // This is, of course, for third order polynomials.
     // For lower orders, the list is truncated at the appropriate level.
 
-    int order = 0;			// Polynomial order
+    int order = 0;                      // Polynomial order
     if (fscanf(mapfp, "%d", &order) != 1) {
-	psLogMsg("stac.read.map", PS_LOG_ERROR, "Unable to read map file %s\n", filename);
-	fclose(mapfp);
-	return NULL;
-    }
-    
+        psLogMsg("stac.read.map", PS_LOG_ERROR, "Unable to read map file %s\n", filename);
+        fclose(mapfp);
+        return NULL;
+    }
+
     psTrace("stac.read.map", 5, "Polynomial order: %d\n", order);
-    
+
     psPlaneTransform *map = psPlaneTransformAlloc(order + 1, order + 1); // The transformation
     // Set coefficient masks
     for (int i = 0; i < order + 1; i++) {
-	for (int j = 0; j < order + 1; j++) {
-	    if (i + j > order + 1) {
-		map->x->mask[i][j] = 1;
-		map->y->mask[i][j] = 1;
-	    } else {
-		map->x->mask[i][j] = 0;
-		map->y->mask[i][j] = 0;
-	    }
-	}
-    }
-    
+        for (int j = 0; j < order + 1; j++) {
+            if (i + j > order + 1) {
+                map->x->mask[i][j] = 1;
+                map->y->mask[i][j] = 1;
+            } else {
+                map->x->mask[i][j] = 0;
+                map->y->mask[i][j] = 0;
+            }
+        }
+    }
+
     // Read x coefficients
     psTrace("stac.read.map", 7, "x' = \n");
     for (int k = 0; k < order + 1; k++) {
-	for (int j = 0; j < k + 1; j++) {
-	    int i = k - j;
-	    if (fscanf(mapfp, "%lf", &map->x->coeff[i][j]) != 1) {
-		psLogMsg("stac.read.map", PS_LOG_ERROR, "Unable to read map file %s\n", filename);
-		fclose(mapfp);
-		psFree(map);
-		return NULL;
-	    }
-	    psTrace("stac.read.map", 7, "      %f x^%d y^%d\n", map->x->coeff[i][j], i, j);
-	}
+        for (int j = 0; j < k + 1; j++) {
+            int i = k - j;
+            if (fscanf(mapfp, "%lf", &map->x->coeff[i][j]) != 1) {
+                psLogMsg("stac.read.map", PS_LOG_ERROR, "Unable to read map file %s\n", filename);
+                fclose(mapfp);
+                psFree(map);
+                return NULL;
+            }
+            psTrace("stac.read.map", 7, "      %f x^%d y^%d\n", map->x->coeff[i][j], i, j);
+        }
     }
     // Read y coefficients
     psTrace("stac.read.maps", 7, "y' = \n");
     for (int k = 0; k < order + 1; k++) {
-	for (int j = 0; j < k + 1; j++) {
-	    int i = k - j;
-	    if (fscanf(mapfp, "%lf", &map->y->coeff[i][j]) != 1) {
-		psLogMsg("stac.read.map", PS_LOG_ERROR, "Unable to read map file %s\n", filename);
-		fclose(mapfp);
-		psFree(map);
-		return NULL;
-	    }
-	    psTrace("stac.read.map", 7, "      %f x^%d y^%d\n", map->y->coeff[i][j], i, j);
-	}
-    }
-    
+        for (int j = 0; j < k + 1; j++) {
+            int i = k - j;
+            if (fscanf(mapfp, "%lf", &map->y->coeff[i][j]) != 1) {
+                psLogMsg("stac.read.map", PS_LOG_ERROR, "Unable to read map file %s\n", filename);
+                fclose(mapfp);
+                psFree(map);
+                return NULL;
+            }
+            psTrace("stac.read.map", 7, "      %f x^%d y^%d\n", map->y->coeff[i][j], i, j);
+        }
+    }
+
     fclose(mapfp);
 
     return map;
 }
-    
+
 
 
@@ -169,20 +177,20 @@
     )
 {
-    int nFiles = filenames->n;		// The number of input files
+    int nFiles = filenames->n;          // The number of input files
     psArray *maps = psArrayAlloc(nFiles); // The maps, to be returned
-    char mapfile[MAXCHAR];		// Filename of map
-    
+    char mapfile[MAXCHAR];              // Filename of map
+
     psTrace("stac.read.maps",1,"Reading maps....\n");
     for (int i = 0; i < nFiles; i++) {
-	if (strlen(filenames->data[i]) > MAXCHAR - 4) {
-	    psLogMsg("stac.read.maps",PS_LOG_ERROR,"Filename %s is too long.\n",filenames->data[i]);
-	    exit(EXIT_FAILURE);
-	}
-	// Read the file
-	sprintf(mapfile,"%s.map",filenames->data[i]);
-	maps->data[i] = stacReadMap(mapfile);
-	if (maps->data[i] == NULL) {
-	    psLogMsg("stac.read.maps", PS_LOG_ERROR, "Unable to read map: %s\n", mapfile);
-	}
+        if (strlen(filenames->data[i]) > MAXCHAR - 4) {
+            psLogMsg("stac.read.maps",PS_LOG_ERROR,"Filename %s is too long.\n",filenames->data[i]);
+            exit(EXIT_FAILURE);
+        }
+        // Read the file
+        sprintf(mapfile,"%s.map",filenames->data[i]);
+        maps->data[i] = stacReadMap(mapfile);
+        if (maps->data[i] == NULL) {
+            psLogMsg("stac.read.maps", PS_LOG_ERROR, "Unable to read map: %s\n", mapfile);
+        }
 
     }
Index: trunk/stac/src/stacRejection.c
===================================================================
--- trunk/stac/src/stacRejection.c	(revision 3828)
+++ 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);
Index: trunk/stac/src/stacScales.c
===================================================================
--- trunk/stac/src/stacScales.c	(revision 3828)
+++ trunk/stac/src/stacScales.c	(revision 5743)
@@ -14,9 +14,9 @@
     assert(image->type.type == PS_TYPE_F32);
 
-// Will use robust median instead of sampling --- it's supposed to be fast.
 #if 1
     int size = image->numCols * image->numRows;	// Number of pixels in image
     int numSamples = size / sample;	// Number of samples in image
     psVector *values = psVectorAlloc(numSamples + 1, PS_TYPE_F32); // Vector containing sub-sample
+    psVector *mask = psVectorAlloc(numSamples + 1, PS_TYPE_U8);	// Mask for sample
 
     int offset = 0;			// Offset from start of the row
@@ -27,4 +27,9 @@
 	while (col < image->numCols) {
 	    values->data.F32[index] = image->data.F32[row][col];
+	    if (isnan(values->data.F32[index])) {
+		mask->data.U8[index] = 1;
+	    } else {
+		mask->data.U8[index] = 0;
+	    }
 	    col += sample;
 	    index++;
@@ -34,9 +39,11 @@
 
     psStats *stats = psStatsAlloc(PS_STAT_SAMPLE_MEDIAN);
-    stats = psVectorStats(stats, values, NULL, NULL, 0);
+    stats = psVectorStats(stats, values, NULL, mask, 1);
     float median = stats->sampleMedian;
     psFree(stats);
     psFree(values);
+    psFree(mask);
 #else
+    // Will use robust median instead of sampling --- it's supposed to be fast.
     psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN); // Using a clipped mean because median is SLOW
     stats->clipSigma = 3.0;
@@ -162,12 +169,16 @@
 			}
 		    }
-		    psTrace("stac.scales", 9, "Star at %f,%f --> %f\n", coords->x, coords->y, sum);
 		    // Subtract background, renormalise to account for circular aperture
 		    if (numPix > 0 && sum > 0) {
 			sum -= offsets->data.F32[i] * (float)numPix;
 			photometry->data.F32[j] = sum * M_PI * aper2 / (float)numPix;
+			if (photometry->data.F32[j] > 0 && finite(photometry->data.F32[j])) {
+			    mask->data.U8[j] = 1;
+			    psTrace("stac.scales", 8, "Star at %f,%f --> %f\n", coords->x, coords->y, sum);
+			} else {
+			    mask->data.U8[j] = 0;
+			}
+		    } else {
 			mask->data.U8[j] = 0;
-		    } else {
-			mask->data.U8[j] = 1;
 		    }
 		}
@@ -182,18 +193,42 @@
 	psVector *ref = stars->data[0];	// The reference photometry
 	psVector *refMask = masks->data[0]; // The reference mask
+#if 0
 	psStats *stats = psStatsAlloc(PS_STAT_CLIPPED_MEAN); // Statistics
 	stats->clipSigma = 2.5;
 	stats->clipIter = 3;
+#else
+	psStats *stats = psStatsAlloc(PS_STAT_SAMPLE_MEDIAN); // Statistics
+#endif
 	for (int i = 1; i < scales->n; i++) {
-	    psVector *compare = stars->data[i];	// The comparison photometry
-	    psVector *compareMask = masks->data[i]; // The comparison mask
-
-	    compare = (psVector*)psBinaryOp(compare, compare, "/", ref);
-	    compareMask = (psVector*)psBinaryOp(compareMask, compareMask, "+", refMask);
-
-	    stats = psVectorStats(stats, compare, NULL, compareMask, 3); // Use maskVal of 3 to catch 1 and 2
-
+	    psVector *input = stars->data[i];	// The comparison photometry
+	    psVector *inputMask = masks->data[i]; // The comparison mask
+
+	    psVector *compare = (psVector*)psBinaryOp(NULL, input, "/", ref);
+	    psVector *compareMask = (psVector*)psBinaryOp(NULL, inputMask, "*", refMask);
+	    (void)psBinaryOp(compareMask, psScalarAlloc(1, PS_TYPE_U8), "-", compareMask);
+
+#if 0
+	    psTrace("stac.scales", 9, "Getting scale for image %d...\n", i);
+	    for (int j = 0; j < compare->n; j++) {
+		if (compareMask->data.U8[j] & 1) {
+		    psTrace("stac.scales", 9, "Bad star %d: %f %f --> %f %d\n", j, input->data.F32[j],
+			    ref->data.F32[j], compare->data.F32[j], compareMask->data.U8[j]);
+		} else {
+		    psTrace("stac.scales", 9, "Good star %d: %f %f --> %f %d\n", j, input->data.F32[j],
+			    ref->data.F32[j], compare->data.F32[j], compareMask->data.U8[j]);
+		}
+	    }
+#endif
+
+	    psVectorStats(stats, compare, NULL, compareMask, 1);
+
+#if 0
 	    scales->data.F32[i] = stats->clippedMean;
+#else
+	    scales->data.F32[i] = stats->sampleMedian;
+#endif
 	    psTrace("stac.scales", 5, "Scale for image %d is %f\n", i, scales->data.F32[i]);
+	    psFree(compare);
+	    psFree(compareMask);
 	}
 	psFree(stats);
Index: trunk/stac/src/stacTime.c
===================================================================
--- trunk/stac/src/stacTime.c	(revision 5743)
+++ trunk/stac/src/stacTime.c	(revision 5743)
@@ -0,0 +1,12 @@
+#include <stdio.h>
+#include <sys/time.h>
+
+double getTime(void)
+/* Gets the current time.  Got this from Nick Kaiser's fetchpix.c */
+{
+    struct timeval tv;
+    gettimeofday(&tv, NULL);
+    return(tv.tv_sec + 1.e-6 * tv.tv_usec);
+}
+
+
Index: trunk/stac/src/stacTransform.c
===================================================================
--- trunk/stac/src/stacTransform.c	(revision 3828)
+++ trunk/stac/src/stacTransform.c	(revision 5743)
@@ -130,4 +130,5 @@
 	for (int i = 0; i < nImages; i++) {
 	    (*outputs)->data[i] = psImageAlloc(outnx, outny, PS_TYPE_F32);
+	    psImageInit((*outputs)->data[i], 0.0);
 	}
     }
@@ -200,5 +201,5 @@
 		    // Change PS_INTERPOLATE_BILINEAR to best available technique.
 		    outImage->data.F32[y][x] = (psF32)psImagePixelInterpolate(image, detector->x,
-									      detector->y, mask, 1, 0.0,
+									      detector->y, mask, 1, NAN,
 									      PS_INTERPOLATE_BILINEAR);
 		    if (error) {
@@ -207,5 +208,5 @@
 												detector->x,
 												detector->y,
-												mask, 1, 0.0);
+												mask, 1, NAN);
 		    }
 
Index: trunk/stac/src/stacWrite.c
===================================================================
--- trunk/stac/src/stacWrite.c	(revision 3828)
+++ trunk/stac/src/stacWrite.c	(revision 5743)
@@ -16,10 +16,10 @@
     }
 
-    psDPolynomial2D *xMap = map->x;	// x transform
-    psDPolynomial2D *yMap = map->y;	// y transform
+    psPolynomial2D *xMap = map->x;	// x transform
+    psPolynomial2D *yMap = map->y;	// y transform
 
     // A crucial limitation of the current system --- the order of each polynomial must be the same
     assert(xMap->nX == xMap->nY && yMap->nX == yMap->nY && xMap->nX == yMap->nX);
-    int order = xMap->nX - 1;	// The polynomial order
+    int order = xMap->nX;	// The polynomial order
     fprintf(mapFile, "%d\n", order);
     
Index: trunk/stac/src/sum.c
===================================================================
--- trunk/stac/src/sum.c	(revision 5743)
+++ trunk/stac/src/sum.c	(revision 5743)
@@ -0,0 +1,83 @@
+// Brain-dead simple stack of multiple frames
+// Makes no effort to zap bad pixels, or rescale images.
+// Assumes that images are already registered
+
+#include <stdio.h>
+#include "pslib.h"
+
+int main(int argc, char *argv[])
+{
+    const char *programName = argv[0];  // Program name
+    const char *outputName = argv[1];   // Output file name
+    psArray *inputNames = psArrayAlloc(argc-2);
+    for (int i = 2; i < argc; i++) {
+        inputNames->data[i-2] = psStringCopy(argv[i]);
+    }
+
+    int numCols = 0;                    // Number of columns
+    int numRows = 0;                    // Number of rows
+    psImage *output = NULL;             // Output image
+    psImage *scale = NULL;              // Scale image
+    psMetadata *header = NULL;          // FITS header for output image
+    for (int i = 0; i < inputNames->n; i++) {
+        psFits *fits = psFitsOpen(inputNames->data[i], "r"); // FITS file
+        psRegion readRegion = {0, 0, 0, 0}; // Region to read
+        psImage *image = psFitsReadImage(NULL, fits, readRegion, 0);
+        if (numCols == 0 && numRows == 0) {
+            // First run through --- set the size, initialise the output, read the header
+            numCols = image->numCols;
+            numRows = image->numRows;
+            output = psImageAlloc(numCols, numRows, PS_TYPE_F32);
+            scale = psImageAlloc(numCols, numRows, PS_TYPE_U8);
+            for (int y = 0; y < numCols; y++) {
+                for (int x = 0; x < numRows; x++) {
+                    output->data.F32[y][x] = 0.0;
+                    scale->data.U8[y][x] = 0;
+                }
+            }
+            header = psFitsReadHeader(NULL, fits);
+        } else if (image->numCols != numCols || image->numRows != numRows) {
+            psError(PS_ERR_IO, true, "Input images have different dimensions.  %s is %dx%d, but should be "
+                    "%dx%d\n", inputNames->data[i], image->numCols, image->numRows, numCols, numRows);
+            return EXIT_FAILURE;
+        }
+        psFitsClose(fits);
+
+        // Convert the type, if necessary
+        if (image->type.type != PS_TYPE_F32) {
+            psImage *temp = psImageCopy(NULL, image, PS_TYPE_F32);
+            psFree(image);
+            image = temp;
+        }
+
+        // Add in
+        for (int y = 0; y < numCols; y++) {
+            for (int x = 0; x < numRows; x++) {
+                if (image->data.F32[y][x] > 0) {
+                    output->data.F32[y][x] += image->data.F32[y][x];
+                    scale->data.U8[y][x]++;
+                }
+            }
+        }
+
+        psFree(image);
+    }
+    psFree(inputNames);
+
+    for (int y = 0; y < numCols; y++) {
+        for (int x = 0; x < numRows; x++) {
+            if (scale->data.U8[y][x] > 0) {
+                output->data.F32[y][x] /= (float)scale->data.U8[y][x];
+            } else {
+                output->data.F32[y][x] = 0.0;
+            }
+        }
+    }
+
+    psFits *fits = psFitsOpen(outputName, "w");
+    (void)psFitsWriteImage(fits, header, output, 0);
+    psFitsClose(fits);
+    psFree(output);
+    // Pau.
+    return EXIT_SUCCESS;
+}