Index: /trunk/stac/src/Makefile
===================================================================
--- /trunk/stac/src/Makefile	(revision 3374)
+++ /trunk/stac/src/Makefile	(revision 3375)
@@ -1,12 +1,15 @@
 SHELL = /bin/sh
 CC = gcc
-CFLAGS += -g -std=c99 -I/home/mithrandir/price/psLib3/psLib/include -DCRFLUX -DTESTING
-PSLIB += -L/home/mithrandir/price/psLib3/psLib/lib/ -lpslib -lgsl -lgslcblas -lfftw3f -lsla -lcfitsio -lm
+CFLAGS += -g -std=c99 -I/home/mithrandir/price/psLib/include -DCRFLUX -DTESTING
+PSLIB += -L/home/mithrandir/price/psLib/lib/ -lpslib -lgsl -lgslcblas -lfftw3f -lsla -lcfitsio -lm -lxml2
 LDFLAGS += $(PSLIB)
 
 STAC = stac.o stacConfig.o stacRead.o stacErrorImages.o stacTransform.o stacCheckMemory.o \
-	stacCombine.o stacInvertMaps.o stacRejection.o psPlaneTransform.o stacAreaOfInterest.o stacSize.o
+	stacCombine.o stacInvertMaps.o stacRejection.o stacAreaOfInterest.o stacSize.o
 
-FITTRANS = fitTrans.o fitTransConfig.o fitTransRead.o stacCheckMemory.o fitTransSolve.o psPlaneTransform.o
+SHIFT = shift.o stacConfig.o stacRead.o stacErrorImages.o stacTransform.o stacCheckMemory.o \
+	stacCombine.o stacInvertMaps.o stacSize.o
+
+FITTRANS = fitTrans.o fitTransConfig.o fitTransRead.o stacCheckMemory.o fitTransSolve.o
 
 .PHONY: tags clean empty test profile optimise
@@ -17,4 +20,7 @@
 stac:		$(STAC)
 		$(CC) $(CFLAGS) -o $@ $(STAC) $(LDFLAGS) $(OPTFLAGS)
+
+shift:		$(SHIFT)
+		$(CC) $(CFLAGS) -o $@ $(SHIFT) $(LDFLAGS) $(OPTFLAGS)
 
 calcGrad:	calcGrad.o
Index: /trunk/stac/src/psPlaneTransform.c
===================================================================
--- /trunk/stac/src/psPlaneTransform.c	(revision 3374)
+++ /trunk/stac/src/psPlaneTransform.c	(revision 3375)
@@ -1,30 +1,0 @@
-#include <stdio.h>
-#include "pslib.h"
-#include "stac.h"
-
-#define PS_INT_CHECK_NON_NEGATIVE(NAME, RVAL) \
-if (NAME < 0) { \
-    psError("pslib.badValue", \
-            "Error: %s is less than 0.", #NAME); \
-    return(RVAL); \
-}
-
-static void planeTransformFree(psPlaneTransform *pt)
-{
-    psFree(pt->x);
-    psFree(pt->y);
-}
-
-psPlaneTransform* psPlaneTransformAlloc(psS32 n1, psS32 n2)
-{
-    PS_INT_CHECK_NON_NEGATIVE(n1, NULL);
-    PS_INT_CHECK_NON_NEGATIVE(n2, NULL);
-
-    psPlaneTransform *pt = psAlloc(sizeof(psPlaneTransform));
-    pt->x = psDPolynomial2DAlloc(n1, n2, PS_POLYNOMIAL_ORD);
-    pt->y = psDPolynomial2DAlloc(n1, n2, PS_POLYNOMIAL_ORD);
-
-    p_psMemSetDeallocator(pt, (psFreeFcn) planeTransformFree);
-    return(pt);
-}
-
Index: /trunk/stac/src/stac.c
===================================================================
--- /trunk/stac/src/stac.c	(revision 3374)
+++ /trunk/stac/src/stac.c	(revision 3375)
@@ -77,7 +77,12 @@
 
 #ifdef TESTING
-    char prefile[MAXCHAR];		// Filename of precombined image
-    sprintf(prefile,"%s.pre",config->output);
-    psImageWriteSection(combined,0,0,0,NULL,0,prefile);
+    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);
+    }
+    psTrace("stac", 1, "Pre-combined image written to %s\n", preName);
+    psFree(preFile);
 #endif
 
@@ -89,7 +94,12 @@
 					      ((psImage*)(inputs->data[i]))->numRows);
 #ifdef TESTING
-	char regionFile[MAXCHAR];	// Filename of region image
-	sprintf(regionFile,"%s.region",config->inputs->data[i]);
-	psImageWriteSection(regions->data[i], 0, 0, 0, NULL, 0, regionFile);
+	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);
 #endif
     }
@@ -122,5 +132,10 @@
 
     // Write output image
-    psImageWriteSection(combined,0,0,0,NULL,0,config->output);
+    psFits *outFile = psFitsAlloc(config->output);
+    if (!psFitsWriteImage(outFile, NULL, combined, 0, NULL)) {
+	psErrorStackPrint(stderr, "Unable to write image: %s\n", config->output);
+    }
+    psTrace("stac", 1, "Combined image written to %s\n", config->output);
+    psFree(outFile);
 
     // Free everything I've used
Index: /trunk/stac/src/stacAreaOfInterest.c
===================================================================
--- /trunk/stac/src/stacAreaOfInterest.c	(revision 3374)
+++ /trunk/stac/src/stacAreaOfInterest.c	(revision 3375)
@@ -1,3 +1,4 @@
 #include <stdio.h>
+#include <assert.h>
 #include "pslib.h"
 #include "stac.h"
@@ -15,8 +16,6 @@
     psDPolynomial2D *yPoly = transform->y;
 
-    if (xPoly->type != PS_POLYNOMIAL_ORD || yPoly->type != PS_POLYNOMIAL_ORD) {
-	psError("stac.area.deriv", "Transformation polynomials aren't ordinary polynomials.\n");
-	return false;
-    }
+    assert(xPoly->type == PS_POLYNOMIAL_ORD);
+    assert(yPoly->type == PS_POLYNOMIAL_ORD);
 
     // Initialise derivatives
@@ -64,16 +63,7 @@
 {
     // Input checks
-    if (mask == NULL) {
-	psError("stac.area", "Input mask is NULL for some weird reason.\n");
-	return NULL;
-    }
-    if (map == NULL) {
-	psError("stac.area", "Input map is NULL for some weird reason.\n");
-	return NULL;
-    }
-    if (mask->type.type != PS_TYPE_U8) {
-	psError("stac.area", "Input mask must be of type U8.\n");
-	return NULL;
-    }
+    assert(mask);
+    assert(map);
+    assert(mask->type.type == PS_TYPE_U8);
     
     // Number of pixels in x and y
Index: /trunk/stac/src/stacCheckMemory.c
===================================================================
--- /trunk/stac/src/stacCheckMemory.c	(revision 3374)
+++ /trunk/stac/src/stacCheckMemory.c	(revision 3375)
@@ -44,9 +44,9 @@
 
     if ((leakFile = fopen(LEAKS, "w")) == NULL) {
-	psError("stac.checkMemory", "Unable to open leaks file, %s\n",LEAKS);
+	fprintf(stderr, "Unable to open leaks file, %s\n", LEAKS);
 	return;
     }
 
-    int nLeaks = psMemCheckLeaks(0, &leaks, leakFile); // Number of leaks
+    int nLeaks = psMemCheckLeaks(0, &leaks, leakFile, false); // Number of leaks
     psTrace("stac.checkMemory", 1, "%d leaks found.\n", nLeaks);
     for (int i = 0; i < nLeaks; i++) {
Index: /trunk/stac/src/stacCombine.c
===================================================================
--- /trunk/stac/src/stacCombine.c	(revision 3374)
+++ /trunk/stac/src/stacCombine.c	(revision 3375)
@@ -1,3 +1,4 @@
 #include <stdio.h>
+#include <assert.h>
 #include <math.h>
 #include "pslib.h"
@@ -15,5 +16,5 @@
     // Would like to use psVectorStats, but it doesn't have errors built in yet
     psStats *stats = psStatsAlloc(PS_STAT_SAMPLE_MEAN);
-    (void)psVectorStats(stats, values, masks, 1);
+    (void)psVectorStats(stats, values, NULL, masks, 1);
     float mean = stats->sampleMean;
     psFree(stats);
@@ -47,5 +48,5 @@
 {
     psStats *stats = psStatsAlloc(PS_STAT_SAMPLE_MEDIAN);
-    (void)psVectorStats(stats, values, masks, 1);
+    (void)psVectorStats(stats, values, NULL, masks, 1);
     float median = stats->sampleMedian;
     psFree(stats);
@@ -76,33 +77,16 @@
 	psImage *error = (psImage *)errors->data[i]; // The error image
 
-	if ((image->numCols != numCols) || (image->numRows != numRows)) {
-	    psError("stac.combine",
-		    "Image size mismatch on error image %d\nExpected %dx%d, got %dx%d\n",
-		    i, numCols, numRows, image->numCols, image->numRows);
-	    return false;
-	}
-	if ((error->numCols != numCols) || (error->numRows != numRows)) {
-	    psError("stac.combine",
-		    "Image size mismatch on error image %d\nExpected %dx%d, got %dx%d\n",
-		    i, numCols, numRows, error->numCols, error->numRows);
-	    return false;
-	}
+	assert(image->numCols == numCols && image->numRows == numRows);
+	assert(error->numCols == numCols && error->numRows == numRows);
     }
 
     // Check combined image
+    assert(!*combined || ((*combined)->numRows == numRows) && ((*combined)->numCols == numCols));
     if (*combined == NULL) {
 	*combined = psImageAlloc(numCols, numRows, PS_TYPE_F32); // Combined image
-    } else if (((*combined)->numRows != numRows) || ((*combined)->numCols != numCols)) {
-	psError("stac.combine", "Size of combined image (%dx%d) does not match inputs (%dx%d)\n",
-		(*combined)->numCols, (*combined)->numRows, numCols, numRows);
-	return false;
     }
 
     // Check area of interest
-    if (region && ((region->numRows != numRows) || (region->numCols != numCols))) {
-	psError("stac.combine", "Size of area of interest (%dx%d) does not match inputs (%dx%d)\n",
-		region->numCols, region->numRows, numCols, numRows);
-	return false;
-    }
+    assert(!region || (region->numRows == numRows) && (region->numCols == numCols));
 
     psTrace("stac.combine", 1, "Combining images....\n");
@@ -117,8 +101,6 @@
 	    // Allocate the rejection masks, if required
 	    *rejected = psArrayAlloc(nImages);
-	} else if ((*rejected)->n != nImages) {
-	    psError("stac.combine", "Number of rejection masks (%d) does not match number of input"
-		    "images (%d).\n", (*rejected)->n, nImages);
-	    return NULL;
+	} else {
+	    assert((*rejected)->n != nImages);
 	}
 
@@ -208,14 +190,25 @@
     if (nReject > 0) {
 	for (int i = 0; i < nImages; i++) {
-	    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);
+	    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);
 	}
     }
     // 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);
+    char chiName[MAXCHAR];		// Filename of chi^2 image
+    sprintf(chiName,"chi2_%d.fits",iteration);
+    psFits *chiFile = psFitsAlloc(chiName);
+    if (!psFitsWriteImage(chiFile, NULL, chi2 , 0, NULL)) {
+	psErrorStackPrint(stderr, "Unable to write image: %s\n", chiName);
+    }
+    psTrace("stac", 1, "Chi^2 image written to %s\n", chiName);
+    psFree(chiFile);
     psFree(chi2);
 #endif
Index: /trunk/stac/src/stacErrorImages.c
===================================================================
--- /trunk/stac/src/stacErrorImages.c	(revision 3374)
+++ /trunk/stac/src/stacErrorImages.c	(revision 3375)
@@ -35,7 +35,12 @@
 #ifdef TESTING
 	// Write error image out to check
-	char errfile[MAXCHAR];		// Filename of error image
-	sprintf(errfile,"%s.err",config->inputs->data[i]);
-	psImageWriteSection(error,0,0,0,NULL,0,errfile);
+	char errName[MAXCHAR];		// Filename of error image
+	sprintf(errName,"%s.err",config->inputs->data[i]);
+	psFits *errorFile = psFitsAlloc(errName);
+	if (!psFitsWriteImage(errorFile, NULL, error, 0, NULL)) {
+	    psErrorStackPrint(stderr, "Unable to write image: %s\n", errName);
+	}
+	psTrace("stac", 1, "Error image written to %s\n", errName);
+	psFree(errorFile);
 #endif
 
Index: /trunk/stac/src/stacInvertMaps.c
===================================================================
--- /trunk/stac/src/stacInvertMaps.c	(revision 3374)
+++ /trunk/stac/src/stacInvertMaps.c	(revision 3375)
@@ -1,3 +1,4 @@
 #include <stdio.h>
+#include <assert.h>
 #include "pslib.h"
 #include "stac.h"
@@ -24,16 +25,6 @@
 	psPlaneTransform *oldMap = (psPlaneTransform*)maps->data[mapNum]; // Uninverted map
 	// Check input
-	if (oldMap->x->nX != oldMap->x->nY) {
-	    psError("stac.invertMaps", "Polynomial order in x and y don't match for map %d\n", mapNum);
-	    return NULL;
-	}
-	if (oldMap->y->nX != oldMap->y->nY) {
-	    psError("stac.invertMaps", "Polynomial order in x and y don't match for map %d\n", mapNum);
-	    return NULL;
-	}
-	if (oldMap->x->nX != oldMap->y->nX) {
-	    psError("stac.invertMaps", "Polynomial order in x and y don't match for map %d\n", mapNum);
-	    return NULL;
-	}
+	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
 	psTrace("stac.invertMaps", 4, "Generating order %d polynomial inverse transformation.\n", order);
@@ -91,5 +82,5 @@
 
 		    fakePoly->mask[i][j] = 0;
-		    double ijPoly = psDPolynomial2DEval(xIn->data.F32[g], yIn->data.F32[g], fakePoly);
+		    double ijPoly = psDPolynomial2DEval(fakePoly, xIn->data.F32[g], yIn->data.F32[g]);
 		    fakePoly->mask[i][j] = 1;
 
@@ -98,5 +89,5 @@
 			    
 			    fakePoly->mask[m][n] = 0;
-			    double mnPoly = psDPolynomial2DEval(xIn->data.F32[g], yIn->data.F32[g], fakePoly);
+			    double mnPoly = psDPolynomial2DEval(fakePoly, xIn->data.F32[g], yIn->data.F32[g]);
 			    fakePoly->mask[m][n] = 1;
 
@@ -113,5 +104,5 @@
 	// Solution via LU Decomposition
 	psVector *permutation = psVectorAlloc(nCoeff, PS_TYPE_F64); // Permutation vector for LU Decomposition
-	psImage *luMatrix = psMatrixLUD(NULL, permutation, matrix); // LU decomposed matrix
+	psImage *luMatrix = psMatrixLUD(NULL, &permutation, matrix); // LU decomposed matrix
 	psVector *xSolution = psMatrixLUSolve(NULL, luMatrix, xVector, permutation); // Solution in x
 	psVector *ySolution = psMatrixLUSolve(NULL, luMatrix, yVector, permutation); // Solution in y
@@ -178,14 +169,8 @@
 #if 0
 	// Can't handle higher order than linear yet
-	if (((psPlaneTransform*)maps->data[i])->x->nX != 2 ||
-	    ((psPlaneTransform*)maps->data[i])->x->nY != 2 ||
-	    ((psPlaneTransform*)maps->data[i])->y->nX != 2 ||
-	    ((psPlaneTransform*)maps->data[i])->y->nY != 2) {
-	    psError("stac.invertMaps",
-		    "STAC cannot currently support orders other than linear.\n");
-	    psFree(inverted);
-	    return NULL;
-	}
-
+	assert(((psPlaneTransform*)maps->data[i])->x->nX == 2 &&
+	       ((psPlaneTransform*)maps->data[i])->x->nY == 2 &&
+	       ((psPlaneTransform*)maps->data[i])->y->nX == 2 &&
+	       ((psPlaneTransform*)maps->data[i])->y->nY == 2);
 	psPlaneTransform *newMap = psPlaneTransformAlloc(2, 2); // Inverted map
 	psPlaneTransform *oldMap = (psPlaneTransform*)maps->data[i]; // Uninverted map
Index: /trunk/stac/src/stacRead.c
===================================================================
--- /trunk/stac/src/stacRead.c	(revision 3374)
+++ /trunk/stac/src/stacRead.c	(revision 3375)
@@ -8,23 +8,34 @@
     psArray *filenames = config->inputs;// The file names
     int nFiles = filenames->n;		// The number of input files
-    psArray *files = psArrayAlloc(nFiles); // The input files, to be returned
+    psArray *images = psArrayAlloc(nFiles); // The input files, to be returned
+    psRegion *imageRegion = psRegionAlloc(0, 0, 0, 0); // Region of image to read
 
-    psTrace("stac.read.images",1,"Reading input images....\n");
+    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
-	psTrace("stac.read.images",2,"Reading input image %s....\n",filenames->data[i]);
-	files->data[i] = psImageReadSection(NULL, 0, 0, 0, 0, 0, NULL, 0, filenames->data[i]);
-	if (files->data[i] == NULL) {
+	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);
 	}
-	// Big assumption: input images are 32-bit.
-	if (((psImage*)(files->data[i]))->type.type != PS_TYPE_F32) {
-	    psError("stac.read.images", "Image %d is not 32-bit\n",i);
+	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 {
+	    images->data[i] = image;
 	}
+	psFree(imageFile);
     }
     psTrace("stac.read.images",1,"%d input images read.\n",nFiles);
+    psFree(imageRegion);
 
-    return files;
+    return images;
 }
 
Index: /trunk/stac/src/stacRejection.c
===================================================================
--- /trunk/stac/src/stacRejection.c	(revision 3374)
+++ /trunk/stac/src/stacRejection.c	(revision 3375)
@@ -1,3 +1,4 @@
 #include <stdio.h>
+#include <assert.h>
 #include "pslib.h"
 #include "stac.h"
@@ -34,5 +35,5 @@
     // Get the median
     psStats *stats = psStatsAlloc(PS_STAT_SAMPLE_MEDIAN);
-    (void)psVectorStats(stats, pixels, mask, 1);
+    (void)psVectorStats(stats, pixels, NULL, mask, 1);
     float median = stats->sampleMedian;
     psFree(stats);
@@ -53,33 +54,13 @@
 
     // Check inputs
-    if (inputs->n != rejected->n) {
-	psError("stac.rejection", "Number of input images (%d) and rejection masks (%d) does not match.\n",
-		inputs->n, rejected->n);
-	return NULL;
-    }
-    if (inputs->n != regions->n) {
-	psError("stac.rejection", "Number of input images (%d) and region masks (%d) does not match.\n",
-		inputs->n, regions->n);
-	return NULL;
-    }
-    if (inputs->n != maps->n) {
-	psError("stac.rejection", "Number of input images (%d) and transformation maps (%d) does "
-		"not match.\n", inputs->n, maps->n);
-	return NULL;
-    }
-    if (inputs->n != inverseMaps->n) {
-	psError("stac.rejection", "Number of input images (%d) and inverse transformation maps (%d) does "
-		"not match.\n", inputs->n, inverseMaps->n);
-	return NULL;
-    }
+    assert(inputs->n == rejected->n);
+    assert(inputs->n == regions->n);
+    assert(inputs->n == maps->n);
+    assert(inputs->n == inverseMaps->n);
+
     for (int i = 0; i < nImages; i++) {
-	if ((((psImage*)(inputs->data[i]))->numRows != ((psImage*)(regions->data[i]))->numRows) ||
-	    (((psImage*)(inputs->data[i]))->numCols != ((psImage*)(regions->data[i]))->numCols)) {
-	    psError("stac.rejection",
-		    "Sizes of input image (%dx%d) and region mask (%dx%d) for input %d do not match.\n",
-		    ((psImage*)(inputs->data[i]))->numCols, ((psImage*)(inputs->data[i]))->numRows,
-		    ((psImage*)(regions->data[i]))->numCols, ((psImage*)(regions->data[i]))->numRows, i);
-	    return NULL;
-	}
+	psImage *input = inputs->data[i];
+	psImage *region = regions->data[i];
+	assert(input->numRows == region->numRows && input->numCols == region->numCols);
     }
 
@@ -114,5 +95,5 @@
 	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");
+	    fprintf(stderr, "Unable to open file for detailed output, %s\n");
 	    return NULL;
 	}
@@ -198,13 +179,29 @@
 #ifdef TESTING
 	// Write error image out to check
-	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);
-	psImageWriteSection(rejmap,0,0,0,NULL,0,rejmapfile);
-	psImageWriteSection(grad,0,0,0,NULL,0,gradfile);
+	char maskName[MAXCHAR];		// Filename of mask image
+	char rejmapName[MAXCHAR]; 	// Filename of rejection image
+	char gradName[MAXCHAR];		// Filename of gradient image
+	sprintf(maskName,"%s.mask",config->inputs->data[i]);
+	sprintf(rejmapName,"%s.rejmap",config->inputs->data[i]);
+	sprintf(gradName,"%s.grad",config->inputs->data[i]);
+
+	psFits *maskFile = psFitsAlloc(maskName);
+	psFits *rejmapFile = psFitsAlloc(rejmapName);
+	psFits *gradFile = psFitsAlloc(gradName);
+	if (!psFitsWriteImage(maskFile, NULL, mask, 0, NULL)) {
+	    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)) {
+	    psErrorStackPrint(stderr, "Unable to write image: %s\n", rejmapName);
+	}
+	psTrace("stac", 1, "Rejection map written to %s\n", rejmapName);
+	if (!psFitsWriteImage(gradFile, NULL, grad, 0, NULL)) {
+	    psErrorStackPrint(stderr, "Unable to write image: %s\n", gradName);
+	}
+	psTrace("stac", 1, "Gradient image written to %s\n", gradName);
+	psFree(maskFile);
+	psFree(rejmapFile);
+	psFree(gradFile);
 	psFree(rejmap);
 	psFree(grad);
Index: /trunk/stac/src/stacSize.c
===================================================================
--- /trunk/stac/src/stacSize.c	(revision 3374)
+++ /trunk/stac/src/stacSize.c	(revision 3375)
@@ -1,3 +1,4 @@
 #include <stdio.h>
+#include <assert.h>
 #include "pslib.h"
 #include "stac.h"
@@ -9,8 +10,5 @@
 {
     int nImages = images->n;		// Number of input images
-    if (nImages != maps->n) {
-	psError("stac.size", "Number of images (%d) and maps (%d) do not match.\n", nImages, maps->n);
-	return false;
-    }
+    assert(maps->n == nImages);
 
     psPlane *inCoord = psAlloc(sizeof(psPlane)); // Input coordinates
Index: /trunk/stac/src/stacTransform.c
===================================================================
--- /trunk/stac/src/stacTransform.c	(revision 3374)
+++ /trunk/stac/src/stacTransform.c	(revision 3375)
@@ -1,3 +1,4 @@
 #include <stdio.h>
+#include <assert.h>
 #include "pslib.h"
 #include "stac.h"
@@ -120,16 +121,9 @@
 
     // Check input sizes
-    if (images->n != maps->n) {
-	psError("stac.transform", "Number of maps (%d) does not match number of images (%d).\n",
-		maps->n, images->n);
-	return false;
-    }
-    if (errors && (images->n != errors->n)) {
-	psError("stac.transform", "Number of error images (%d) does not match number of images (%d).\n",
-		errors->n, images->n);
-	return false;
-    }
+    assert(images->n == maps->n);
+    assert(!errors || (images->n == errors->n));
 
     // Allocate the output images if required, otherwise check the number
+    assert(!*outputs || (*outputs)->n == nImages);
     if (*outputs == NULL) {
 	*outputs = psArrayAlloc(nImages);
@@ -138,11 +132,8 @@
 	    (*outputs)->data[i] = psImageAlloc(nx, ny, PS_TYPE_F32);
 	}
-    } else if ((*outputs)->n != nImages) {
-	psError("stac.transform", "Number of output images (%d) does not match number of input images "
-		"(%d).\n", (*outputs)->n, nImages);
-	return false;
     }
 
     // Allocate the output error images, if required, otherwise check the number
+    assert(!errors || ! *outErrors || errors->n == (*outErrors)->n);
     if (errors && (*outErrors == NULL)) {
 	*outErrors = psArrayAlloc(errors->n);
@@ -151,26 +142,13 @@
 	    (*outErrors)->data[i] = psImageAlloc(nx, ny, PS_TYPE_F32);
 	}
-    } else if (errors->n != (*outErrors)->n) {
-	psError("stac.transform", "Number of error output images (%d) does not match number of input"
-		"images (%d).\n", (*outErrors)->n, errors->n);
-	return false;
     }
 
     // Check the masks, if specified
+    assert(!masks || masks->n == nImages);
     if (masks != NULL) {
-	if (masks->n != nImages) {
-	    psError("stac.transform", "Number of masks (%d) does not match number of input images (%d).\n",
-		    masks->n, nImages);
-	    return false;
-	}
 	for (int i = 0; i < nImages; i++) {
 	    psImage *image = images->data[i];
 	    psImage *mask = masks->data[i];
-	    if ((mask->numRows != image->numRows) || (mask->numCols != image->numCols)) {
-		psError ("stac.transform",
-			 "Size of input mask (%dx%d) does not match that of input image (%dx%d).\n",
-			 mask->numCols, mask->numRows, image->numCols, image->numRows);
-		return false;
-	    }
+	    assert(mask->numRows == image->numRows && mask->numCols == image->numCols);
 	}
     }
@@ -236,12 +214,23 @@
 #ifdef TESTING
 	// Write error image out to check
-	char shiftfile[MAXCHAR];	// Filename of shift image
-	char errfile[MAXCHAR];		// Filename of error image
-	sprintf(shiftfile,"%s.shift.%d",config->inputs->data[n],numTransforms);
-	sprintf(errfile,"%s.shifterr.%d",config->inputs->data[n],numTransforms);
+	char shiftName[MAXCHAR];	// Filename of shift image
+	char errName[MAXCHAR];		// Filename of error image
+	sprintf(shiftName,"%s.shift.%d",config->inputs->data[n],numTransforms);
+	sprintf(errName,"%s.shifterr.%d",config->inputs->data[n],numTransforms);
 	psTrace("stac.transform.test", 6,
 		"Output files have size: %dx%d\n",outImage->numCols,outImage->numRows);
-	psImageWriteSection(outImage,0,0,0,NULL,0,shiftfile);
-	psImageWriteSection(outError,0,0,0,NULL,0,errfile);
+
+	psFits *shiftFile = psFitsAlloc(shiftName);
+	psFits *errFile = psFitsAlloc(errName);
+	if (!psFitsWriteImage(shiftFile, NULL, outImage, 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, outError, 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