Index: /trunk/stac/src/Makefile
===================================================================
--- /trunk/stac/src/Makefile	(revision 2500)
+++ /trunk/stac/src/Makefile	(revision 2500)
@@ -0,0 +1,60 @@
+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)
+
+OBJECTS = stac.o stacConfig.o stacRead.o stacErrorImages.o stacTransform.o stacCheckMemory.o \
+	stacCombine.o stacInvertMaps.o stacRejection.o psPlaneTransform.o
+
+TARGET = stac
+
+.PHONY: tags clean empty test profile optimise
+
+.c.o:
+		$(CC) -c $(CFLAGS) $(OPTFLAGS) $<
+
+$(TARGET):	$(OBJECTS)
+		$(CC) $(CFLAGS) -o $@ $(OBJECTS) $(LDFLAGS) $(OPTFLAGS)
+
+clean:
+		-$(RM) *.o gmon.* profile.txt
+
+empty:		clean
+		-$(RM) $(TARGET) TAGS
+
+test:		stac test_0.fits test_1.fits test_2.fits test_3.fits
+		-$(RM) testout.fits
+		-$(RM) testout.fits.pre
+		-$(RM) test_[0-3].fits.err
+		-$(RM) test_[0-3].fits.mask
+		-$(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
+
+# Run profiling.
+profile:	CFLAGS += -pg
+profile:	empty $(TARGET)
+		$(RM) gmon.*
+		for ((i = 0; i < 10; i++)) \
+		do \
+			$(MAKE) test; \
+			mv -f gmon.out gmon.$$i; \
+		done
+		gprof -p -q -l $(TARGET) gmon.* > profile.txt
+
+# Do gcc optimisation
+optimise:	clean
+		-$(RM) $(TARGET)
+		export OPTFLAGS=-fprofile-arcs ; $(MAKE) test
+		-$(RM) $(TARGET)
+		$(MAKE) clean
+		export OPTFLAGS=-fbranch-probabilities ; $(MAKE) $(TARGET)
+		-$(RM) *.da
+
+# Tags for emacs
+tags:
+		etags `find . -name \*.[ch] -print`
+
Index: /trunk/stac/src/psPlaneTransform.c
===================================================================
--- /trunk/stac/src/psPlaneTransform.c	(revision 2500)
+++ /trunk/stac/src/psPlaneTransform.c	(revision 2500)
@@ -0,0 +1,30 @@
+#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 2500)
+++ /trunk/stac/src/stac.c	(revision 2500)
@@ -0,0 +1,91 @@
+#include <stdio.h>
+#include "pslib.h"
+#include "stac.h"
+
+
+int main(int argc, char **argv)
+{
+
+#if 0
+    psMemAllocateCallbackSet(stacMemPrint);
+    psMemAllocateCallbackSetID(185);
+    psMemFreeCallbackSet(stacMemPrint);
+    psMemFreeCallbackSetID(185);
+#endif
+
+    // Set trace levels
+    psTraceSetLevel(".",0);
+    psTraceSetLevel("stac.checkMemory",10);
+    psTraceSetLevel("stac.config",10);
+    psTraceSetLevel("stac.read",10);
+    psTraceSetLevel("stac.invertMaps",10);
+    psTraceSetLevel("stac.errors",10);
+    psTraceSetLevel("stac.transform",10);
+    psTraceSetLevel("stac.combine",10);
+    psTraceSetLevel("stac.rejection",10);
+
+    // Set logging level
+    psLogSetLevel(9);
+
+    // Parse command line
+    stacConfig *config = stacParseConfig(argc, argv);
+
+    // Read input files
+    psArray *inputs = stacReadImages(config);
+
+    // Read maps
+    psArray *maps = stacReadMaps(config);
+
+    // Invert maps
+    psArray *inverseMaps = stacInvertMaps(maps);
+
+    // Generate errors
+    psArray *errors = stacErrorImages(inputs, config);
+
+    // Transform inputs and errors
+    psArray *transformedErrors = NULL;
+    psArray *transformed = stacTransform(inputs, inverseMaps, errors, &transformedErrors, NULL, config);
+
+    // Combine with rejection
+    psArray *rejected = NULL;
+    psImage *combined = stacCombine(transformed, transformedErrors, 1,  &rejected, config);
+
+#ifdef TESTING
+    char prefile[MAXCHAR];		// Filename of precombined image
+    sprintf(prefile,"%s.pre",config->output);
+    psImageWriteSection(combined,0,0,0,NULL,0,prefile);
+#endif
+
+    // Transform rejected pixels to source frame
+    psArray *rejectedSource = stacRejection(inputs, transformed, combined, maps, rejected, config);
+
+    // Redo transformation with the mask
+    psFree(transformed);
+    psFree(transformedErrors);
+    psArray *transformedErrorsNew = NULL;
+    psArray *transformedNew = stacTransform(inputs, inverseMaps, errors, &transformedErrorsNew,
+					    rejectedSource, config);
+
+    // Combine the newly-transformed CR-free images, no rejection
+    psFree(combined);
+    psFree(rejected);
+    rejected = NULL;
+    combined = stacCombine(transformedNew, transformedErrorsNew, 0,  &rejected, config);
+
+    // Write output image
+    psImageWriteSection(combined,0,0,0,NULL,0,config->output);
+
+    // Free everything I've used
+    stacConfigFree(config);
+    psFree(inputs);
+    psFree(maps);
+    psFree(inverseMaps);
+    psFree(errors);
+    psFree(transformedErrors);
+    psFree(transformedNew);
+    psFree(rejectedSource);
+    psFree(combined);
+
+    // Check memory for leaks, corruption
+    stacCheckMemory();
+}
Index: /trunk/stac/src/stac.h
===================================================================
--- /trunk/stac/src/stac.h	(revision 2500)
+++ /trunk/stac/src/stac.h	(revision 2500)
@@ -0,0 +1,141 @@
+#ifndef STAC_H
+#define STAC_H
+
+#include "pslib.h"
+
+#define abs(x) ((x) >= 0 ? (x) : (-(x)))
+
+#define MAXCHAR 80
+
+// Temporary implementation
+psPlaneTransform *psPlaneTransformAlloc(psS32 n1, psS32 n2);
+
+/************************************************************************************************************/
+
+// stacConfig.c
+
+// Configuration options
+typedef struct {
+    int verbose;			// Verbosity level
+    float gain, readnoise;		// Gain and readnoise for detectors
+    psArray *inputs;			// Input file names
+    const char *output;			// Output file name
+    int outnx, outny;			// Size of output image
+    float saturated;			// Saturation level
+    float bad;				// Bad level
+    float reject;			// Rejection level
+    float frac;				// Fraction of input pixel that must be masked before the pixel is
+					// considered bad
+} stacConfig;
+
+// Allocator
+stacConfig *stacConfigAlloc(void);
+// Deallocator
+void stacConfigFree(stacConfig *config);
+
+// Help message
+void help(const char *name);
+
+// Parse the command line and return config
+stacConfig *stacParseConfig(int argc,	// Number of command-line arguments
+			    char **argv // Command-line arguments
+    );
+
+
+/************************************************************************************************************/
+
+// stacRead.c
+
+// Read the input files and return an array of images
+psArray *stacReadImages(stacConfig *config);
+
+// Read the map files and return an array of transformations
+psArray *stacReadMaps(stacConfig *config);
+
+/************************************************************************************************************/
+
+// stacErrorImages.c
+
+// Calculate the error images
+psArray *stacErrorImages(psArray *inputs, // Array of input images
+			 stacConfig *config // Configuration details
+    );
+
+/************************************************************************************************************/
+
+// stacTransform.c
+
+// Transform input images
+psArray *stacTransform(const psArray *images, // Array of images to be transformed
+		       const psArray *maps, // Array of polynomials that do the transformation
+		       const psArray *errors, // Array of error images
+		       psArray **outErrors, // Array of output errors
+		       const psArray *masks, // Masks of input images
+		       const stacConfig *config	// Configuration
+    );
+
+/************************************************************************************************************/
+
+// stacCheckMemory.c
+
+// Check memory
+void stacCheckMemory(void);
+
+// Print out the problem
+void stacMemoryProblem(const psMemBlock* ptr, ///< the pointer to the problematic memory block.
+		       const char *file, ///< the file in which the problem originated
+		       psS32 lineno	///< the line number in which the problem originated
+    );
+
+// Print out a memblock when it's allocated --- this function used as a callback
+psMemoryId stacMemPrint(const psMemBlock *ptr);
+
+
+/************************************************************************************************************/
+
+// stacCombine.c
+
+// Get the mean for a bunch of values
+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
+    );
+
+// Get the median for a bunch of values
+float stacCombineMedian(psVector *values, // Values for which to take the median
+			psVector *errors, // Errors in the values
+			psVector *masks	// Masks for the values, 0 = don't use, 1 = use
+    );
+
+// Combine the transformed images
+psImage *stacCombine(psArray *images,	// Array of transformed images
+		     psArray *errors,	// Array of transformed error images
+		     int nReject,	// Number of rejection iterations
+		     psArray **rejected, // Array of rejection masks
+		     stacConfig *config	// Configuration
+    );
+
+/************************************************************************************************************/
+
+// stacInvertMaps.c
+
+// Invert an array of maps
+psArray *stacInvertMaps(const psArray *maps // Array of maps to invert
+    );
+
+
+/************************************************************************************************************/
+
+// stacRejection.c
+
+// Transform the rejection masks back to the source frame
+psArray *stacRejection(psArray *inputs,	// Input images
+		       psArray *transformed, // Transformed images
+		       psImage *combined, // Combined image
+		       psArray *maps,	// Maps from input to transformed image
+		       psArray *rejected, // Rejected images
+		       stacConfig *config // Configuration
+    );
+
+
+#endif
Index: /trunk/stac/src/stacCheckMemory.c
===================================================================
--- /trunk/stac/src/stacCheckMemory.c	(revision 2500)
+++ /trunk/stac/src/stacCheckMemory.c	(revision 2500)
@@ -0,0 +1,62 @@
+#include <stdio.h>
+#include "pslib.h"
+#include "stac.h"
+
+
+#define LEAKS "leaks.dat"		// File to which to write leaks data
+
+
+psMemoryId stacMemPrint(const psMemBlock *ptr)
+{
+    psLogMsg("stac.memoryPrint", PS_LOG_INFO,
+	     "Memory block %d:\n"
+	     "\tFile %s, line %d, size %d\n"
+	     "\tPosts: %x %x %x\n",
+	     ptr->id, ptr->file, ptr->lineno, ptr->userMemorySize, ptr->startblock, ptr->endblock,
+	     *(void**)((int8_t *)(ptr + 1) + ptr->userMemorySize));
+    return 0;
+}
+
+
+void stacMemoryProblem(const psMemBlock* ptr, ///< the pointer to the problematic memory block.
+		       const char *file, ///< the file in which the problem originated
+		       psS32 lineno	///< the line number in which the problem originated
+    )
+{
+    psLogMsg("stac.checkMemory.corruption", PS_LOG_WARN,
+	     "Memory corruption detected in memBlock %d\n"
+	     "\tFile %s, line %d, size %d\n"
+	     "\tPosts: %x %x %x\n",
+	     ptr->id, file, lineno, ptr->userMemorySize, ptr->startblock, ptr->endblock,
+	     (ptr + 1 + ptr->userMemorySize));
+}
+
+
+
+void stacCheckMemory(void)
+{
+    psMemBlock **leaks = NULL;		// List of leaks
+    FILE *leakFile;			// File to write leaks to
+
+    psTrace("stac.checkMemory", 1, "Checking for memory problems....\n");
+
+    (void)psMemProblemCallbackSet(stacMemoryProblem); // Set callback for corruption
+
+    if ((leakFile = fopen(LEAKS, "w")) == NULL) {
+	psError("stac.checkMemory", "Unable to open leaks file, %s\n",LEAKS);
+	return;
+    }
+
+    int nLeaks = psMemCheckLeaks(0, &leaks, leakFile); // Number of leaks
+    psTrace("stac.checkMemory", 1, "%d leaks found.\n", nLeaks);
+    for (int i = 0; i < nLeaks; i++) {
+	psLogMsg("stac.checkMemory.leaks", PS_LOG_WARN,
+		 "Memory leak detection: memBlock %d\n"
+		 "\tFile %s, line %d, size %d\n",
+		 leaks[i]->id, leaks[i]->file, leaks[i]->lineno, leaks[i]->userMemorySize);
+    }
+
+    int nCorrupted = psMemCheckCorruption(false); // Number of corrupted
+    psTrace("stac.checkMemory", 1, "%d memory blocks corrupted.\n", nCorrupted);
+
+}
Index: /trunk/stac/src/stacCombine.c
===================================================================
--- /trunk/stac/src/stacCombine.c	(revision 2500)
+++ /trunk/stac/src/stacCombine.c	(revision 2500)
@@ -0,0 +1,176 @@
+#include <stdio.h>
+#include <math.h>
+#include "pslib.h"
+#include "stac.h"
+
+#define SQUARE(x) ((x)*(x))
+#define ABS(x) ((x) >= 0 ? (x) : -(x))
+
+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
+    )
+{
+#if 0
+    // Would like to use psVectorStats, but it doesn't have errors built in yet
+    psStats *stats = psStatsAlloc(PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV);
+    (void)psVectorStats(stats, values, masks, 1);
+    psFree(stats);
+    return stats->sampleMean;
+#else
+    // Instead, do it ourselves
+    double sum = 0.0;
+    double weights = 0.0;
+    int num = values->n;
+    for (int i = 0; i < num; i++) {
+	if (masks->data.U8[i]) {
+	    sum += values->data.F32[i] / SQUARE(errors->data.F32[i]);
+	    weights += 1.0 / SQUARE(errors->data.F32[i]);
+	}
+    }
+    if (weights > 0.0) {
+	return (float)(sum/weights);
+    } else {
+	return NAN;
+    }
+#endif
+
+}
+
+
+float stacCombineMedian(psVector *values, // Values for which to take the median
+			psVector *errors, // Errors in the values
+			psVector *masks	// Masks for the values, 0 = don't use, 1 = use
+    )
+{
+    psStats *stats = psStatsAlloc(PS_STAT_SAMPLE_MEDIAN);
+    (void)psVectorStats(stats, values, masks, 1);
+    psFree(stats);
+    return stats->sampleMedian;
+}
+
+
+psImage *stacCombine(psArray *images,	// Array of transformed images
+		     psArray *errors,	// Array of transformed error images
+		     int nReject,	// Number of rejection iterations
+		     psArray **rejected, // Array of rejection masks
+		     stacConfig *config	// Configuration
+    )
+{
+    int nImages = images->n;		// Number of images
+    int numRows = ((psImage*)images->data[0])->numRows;	// Image size
+    int numCols = ((psImage*)images->data[0])->numCols; // Image size
+    float saturated = config->saturated;// Saturation limit
+    float bad = config->bad;		// Bad pixel limit
+    float reject = config->reject;	// Rejection (k-sigma)
+    
+    // Check dimensions for consistency
+    for (int i = 0; i < nImages; i++) {
+	psImage *image = (psImage *)images->data[i]; // The image
+	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 NULL;
+	}
+	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 NULL;
+	}
+    }
+
+    psTrace("stac.combine", 1, "Combining images....\n");
+    psTrace("stac.combine", 3, "Saturation: %f Bad: %f\n", saturated, bad);
+
+    psVector *pixels = psVectorAlloc(nImages, PS_TYPE_F32); // Will hold the pixels in the statistics step
+    psVector *deltas = psVectorAlloc(nImages, PS_TYPE_F32); // Will hold the errors in the statistics step
+    psVector *mask = psVectorAlloc(nImages, PS_TYPE_U8); // Mask bad pixels
+    psImage *combined = psImageAlloc(numCols, numRows, PS_TYPE_F32); // Combined image
+
+
+    // Set up rejection masks
+    if (nReject > 0) {
+	if (*rejected == NULL) {
+	    // 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;
+	}
+	for (int i = 0; i < nImages; i++) {
+	    (*rejected)->data[i] = psImageAlloc(numCols, numRows, PS_TYPE_F32); // Create rejection mask
+	    // Zero out the mask
+	    for (int r = 0; r < numRows; r++) {
+		for (int c = 0; c < numCols; c++) {
+		    ((psImage*)((*rejected)->data[i]))->data.F32[r][c] = 0.0;
+		}
+	    }
+	}
+    }
+    
+    
+    for (int y = 0; y < numRows; y++) {
+	for (int x = 0; x < numCols; x++) {
+	    
+	    // Export pixels into the vector and get stats
+	    for (int i = 0; i < nImages; i++) {
+		float pixel = ((psImage*)images->data[i])->data.F32[y][x];
+		float delta = ((psImage*)errors->data[i])->data.F32[y][x];
+		pixels->data.F32[i] = pixel;
+		deltas->data.F32[i] = delta;
+		if ((pixel >= saturated) || (pixel <= bad) || (! isfinite(pixel)) || (! isfinite(delta))) {
+		    mask->data.U8[i] = (psU8)0; // Don't use!
+		} else {
+		    mask->data.U8[i] = (psU8)1; // Use.
+		}
+	    }
+	    
+	    float average = stacCombineMean(pixels, deltas, mask); // Combined value
+	    
+	    // Rejection iterations
+	    for (int rejNum = 0; rejNum < nReject; rejNum++) {
+		float max = 0.0;
+		int maxIndex = 0;
+		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];
+			maxIndex = i;
+		    }
+		}
+		// Reject the pixel with the maximum deviation
+		if (max > reject) {    
+		    mask->data.U8[maxIndex] = 0;
+		    ((psImage*)((*rejected)->data[maxIndex]))->data.F32[y][x] += 1.0;
+		    // Re-do combination following rejection
+		    average = stacCombineMean(pixels, deltas, mask);
+		}
+	    } // Rejection iterations
+	    
+	    combined->data.F32[y][x] = average;
+	}
+    } // Iterating over output pixels
+
+#ifdef TESTING
+    // Write rejection images out to check
+    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);
+	}
+    }
+#endif
+
+
+    // Clean up
+    psFree(pixels);
+    psFree(deltas);
+    psFree(mask);
+
+    return combined;
+}
Index: /trunk/stac/src/stacConfig.c
===================================================================
--- /trunk/stac/src/stacConfig.c	(revision 2500)
+++ /trunk/stac/src/stacConfig.c	(revision 2500)
@@ -0,0 +1,156 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <unistd.h>
+#include <getopt.h>
+#include <string.h>
+#include "stac.h"
+
+void help(const char *name)
+{
+    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"
+	     "where\n"
+	     "\t-h           Help (this info)\n"
+	     "\t-v           Increase verbosity level\n"
+	     "\t-g           Gain (e/ADU) for detectors\n"
+	     "\t-r           Read noise (e) for detectors\n"
+	     "\t-o NX NY     Size of output image (512, 512)\n"
+	     "\t-s SAT       Saturation level (65536)\n"
+	     "\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"
+	     "\tOUT          Output image\n"
+	     "\tIN1, IN2...  Input images, which have associated .map files.\n",
+	     name
+	);
+}
+
+
+stacConfig *stacConfigAlloc(void)
+{
+    // Allocate memory
+    stacConfig *config = (stacConfig*) psAlloc(sizeof(stacConfig));
+    // Set defaults
+    config->verbose = 0;
+    config->gain = 1.0;
+    config->readnoise = 0.0;
+    config->inputs = NULL;
+    config->output = NULL;
+    config->outnx = 512;
+    config->outny = 512;
+    config->saturated = 65536.0;
+    config->bad = 0.0;
+    config->reject = 2.75;
+    config->frac = 0.5;
+
+    return config;
+}
+
+
+void stacConfigFree(stacConfig *config)
+{
+    // Free the vector, if necessary
+    if (config->inputs) {
+	psFree(config->inputs);
+    }
+    // Free everything
+    psFree(config);
+}
+
+
+stacConfig *stacParseConfig(int argc,	// Number of command-line arguments
+			    char **argv // Command-line arguments
+    )
+{
+    stacConfig *config = stacConfigAlloc(); // Configuration values
+    const char *programName = argv[0];	// Program name
+
+    /* Variables for getopt */
+    int opt;   /* Option, from getopt */
+    extern char *optarg;   /* Argument accompanying switch */
+    extern int optind;   /* getopt variables */
+
+    /* Parse command-line arguments using getopt */
+    while ((opt = getopt(argc, argv, "hvg:r:o:s:b:k:f:")) != -1) {
+        switch (opt) {
+	  case 'h':
+	    help(programName);
+	    exit(EXIT_SUCCESS);
+	  case 'v':
+            config->verbose++;
+            break;
+	  case 'g':
+	    if (sscanf(optarg, "%f", &config->gain) != 1) {
+		help(programName);
+	    }
+	    break;
+	  case 'r':
+	    if (sscanf(optarg, "%f", &config->readnoise) != 1) {
+		help(programName);
+	    }
+	    break;
+	  case 'o':
+            if ((sscanf(argv[optind-1], "%d", &config->outnx) != 1) ||
+                (sscanf(argv[optind++], "%d", &config->outny) != 1)) {
+                /*
+                  Note: incrementing optind, so I can read more than
+                  one parameter.
+                */
+                help(programName);
+                exit(EXIT_FAILURE);
+            }
+	    break;
+	  case 's':
+	    if (sscanf(optarg, "%f", &config->saturated) != 1) {
+		help(programName);
+	    }
+	    break;
+	  case 'b':
+	    if (sscanf(optarg, "%f", &config->bad) != 1) {
+		help(programName);
+	    }
+	    break;
+	  case 'k':
+	    if (sscanf(optarg, "%f", &config->reject) != 1) {
+		help(programName);
+	    }
+	    break;
+	  case 'f':
+	    if (sscanf(optarg, "%f", &config->frac) != 1) {
+		help(programName);
+	    }
+	    break;
+	  default:
+	    help(programName);
+	}
+    }
+
+    /* Get the remaining, mandatory values */
+    argc -= optind;
+    argv += optind;
+    if (argc < 2) {
+        help(programName);
+        exit(EXIT_FAILURE);
+    }
+
+    config->output = argv[0];		// Output file
+    // Get the input files
+    config->inputs = psArrayAlloc(argc-1);
+    for (int i = 1; i < argc; i++) {
+	config->inputs->data[i-1] = psAlloc(strlen(argv[i]));
+	strncpy(config->inputs->data[i-1], argv[i], strlen(argv[i]));
+    }
+
+    // Debugging output
+    psTrace("stac.config", 8, "Parsed command line for configuration\n");
+    psTrace("stac.config", 9, "Verbosity level: %d\n",config->verbose);
+    psTrace("stac.config", 9, "%d inputs:\n",config->inputs->n);
+    for (int i = 0; i < config->inputs->n; i++) {
+	psTrace("stac.config", 9, "\t%s\n", config->inputs->data[i]);
+    }
+    psTrace("stac.config", 9, "Output file is %s\n",config->output);
+    psTrace("stac.config", 9, "Output file size is %dx%d\n", config->outnx, config->outny);
+
+    return config;
+}
+
Index: /trunk/stac/src/stacErrorImages.c
===================================================================
--- /trunk/stac/src/stacErrorImages.c	(revision 2500)
+++ /trunk/stac/src/stacErrorImages.c	(revision 2500)
@@ -0,0 +1,45 @@
+#include <stdio.h>
+#include "pslib.h"
+#include "stac.h"
+
+psArray *stacErrorImages(psArray *inputs, // Array of input images
+			 stacConfig *config // Configuration details
+    )
+{
+    float invGain = 1.0/config->gain;	// Inverse square root of gain
+    float rn = config->readnoise/config->gain; // Read noise in ADU
+    psArray *errors = psArrayAlloc(inputs->n);
+
+    psTrace("stac.errors", 1, "Calculating error images....\n");
+
+    // Iterate over the input images
+    for (int i = 0; i < inputs->n; i++) {
+	psTrace("stac.errors",5,"Working on image #%d\n",i);
+
+	psImage *image = inputs->data[i]; // Pull out the image of interest
+	int numRows = image->numRows;	// Number of rows
+	int numCols = image->numCols;	// Number of columns	
+	psImage *error = psImageAlloc(numCols, numRows, PS_TYPE_F32); // The error image
+
+	// Iterate over the pixels
+	for (int r = 0; r < numRows; r++) {
+	    for (int c = 0; c < numCols; c++) {
+		// We actually calculate the variance
+		error->data.F32[r][c] = image->data.F32[r][c]*invGain + rn*rn;
+	    }
+	}
+
+	// Put image onto the array
+	errors->data[i] = error;
+ 
+#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);
+#endif
+
+    }
+
+    return errors;
+}
Index: /trunk/stac/src/stacInvertMaps.c
===================================================================
--- /trunk/stac/src/stacInvertMaps.c	(revision 2500)
+++ /trunk/stac/src/stacInvertMaps.c	(revision 2500)
@@ -0,0 +1,72 @@
+#include <stdio.h>
+#include "pslib.h"
+#include "stac.h"
+
+
+#define MAX(x,y) ((x) > (y) ? (x) : (y))
+
+psArray *stacInvertMaps(const psArray *maps // Array of maps to invert
+    )
+{
+    int nMaps = maps->n;		// Number of maps
+    psArray *inverted = psArrayAlloc(nMaps); // Array of inverted maps for output
+
+    psTrace("stac.invertMaps", 1, "Inverting maps....\n");
+
+    for (int i = 0; i < nMaps; i++) {
+	// 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;
+	}
+
+	psPlaneTransform *newMap = psPlaneTransformAlloc(2, 2); // Inverted map
+	psPlaneTransform *oldMap = (psPlaneTransform*)maps->data[i]; // Uninverted map
+
+	// Now, simply do a 2x2 matrix inversion
+
+	double a = oldMap->x->coeff[1][0];
+	double b = oldMap->x->coeff[0][1];
+	double c = oldMap->y->coeff[1][0];
+	double d = oldMap->y->coeff[0][1];
+	double e = oldMap->x->coeff[0][0];
+	double f = oldMap->y->coeff[0][0];
+
+	double invDet = 1.0 / (a * d - b * c); // Inverse of the determinant
+
+	// Not entirely sure why this works, but it appears to do so.......................................
+	newMap->x->coeff[1][0] = invDet * a;
+	newMap->x->coeff[0][1] = - invDet * b;
+	newMap->y->coeff[1][0] = - invDet * c;
+	newMap->y->coeff[0][1] = invDet * d;
+
+	newMap->x->coeff[0][0] = - invDet * (d * e + c * f);
+	newMap->y->coeff[0][0] = - invDet * (b * e + a * f);
+
+#ifdef TESTING
+	// Go forward then backward
+	double x = 123.4;
+	double y = 432.1;
+	psPlane *oldCoords = psAlloc(sizeof(psPlane));
+	oldCoords->x = x;
+	oldCoords->y = y;
+	psPlane *newCoords = psPlaneTransformApply(NULL, oldMap, oldCoords);
+	psTrace("stac.invertMaps.test", 5, "%f,%f --> %f,%f\n", x, y, newCoords->x, newCoords->y);
+	(void)psPlaneTransformApply(oldCoords, newMap, newCoords);
+	psTrace("stac.invertMaps.test", 5, "--------> %f,%f\n", oldCoords->x, oldCoords->y);
+	psFree(newCoords);
+	psFree(oldCoords);
+#endif
+
+	inverted->data[i] = newMap;	// Stuff into the array
+    }
+
+    return inverted;
+}
+
+
Index: /trunk/stac/src/stacRead.c
===================================================================
--- /trunk/stac/src/stacRead.c	(revision 2500)
+++ /trunk/stac/src/stacRead.c	(revision 2500)
@@ -0,0 +1,114 @@
+#include <stdio.h>
+#include <string.h>
+#include "pslib.h"
+#include "stac.h"
+
+#define ORDER 2
+
+psArray *stacReadImages(stacConfig *config)
+{
+    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
+
+    psTrace("stac.read.images",1,"Reading input images....\n");
+    for (int i = 0; i < nFiles; 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) {
+	    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",1,"%d input images read.\n",nFiles);
+
+    return files;
+}
+
+
+psArray *stacReadMaps(stacConfig *config)
+{
+    psArray *filenames = config->inputs;// The file names
+    int nFiles = filenames->n;		// The number of input files
+    psArray *maps = psArrayAlloc(nFiles); // The maps, to be returned
+    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);
+	}
+	// Open the file
+	sprintf(mapfile,"%s.map",filenames->data[i]);
+	FILE *mapfp = fopen(mapfile,"r");
+	if (mapfp == NULL) {
+	    psLogMsg("stac.read.maps", PS_LOG_ERROR, "Cannot open map file, %s\n",mapfile);
+	    exit(EXIT_FAILURE);
+	}
+	// Read the file
+	psTrace("stac.read.maps",5,"Reading map file %s....\n",mapfile);
+	psPlaneTransform *map = psPlaneTransformAlloc(ORDER,ORDER);
+
+	// Format is as following:
+	// A0  A1
+	// B00 B01
+	// B10 B11
+	//
+	// where:
+	// x_sky = A0 + B00*x + B10*y
+	// y_sky = A1 + B01*x + B11*y
+
+	double a0, a1, b00, b10, b01, b11;
+
+	if (fscanf(mapfp,"%lf %lf %lf %lf %lf %lf", &a0, &a1, &b00, &b01, &b10, &b11) != 6) {
+	    psLogMsg("stac.read.maps", PS_LOG_ERROR, "Couldn't read all 6 transformation values from %s\n",
+		     mapfile);
+	    exit(EXIT_FAILURE);
+	}
+	fclose(mapfp);
+	map->x->coeff[0][0] = a0;
+	map->x->coeff[1][0] = b00;
+	map->x->coeff[0][1] = b10;
+	map->y->coeff[0][0] = a1;
+	map->y->coeff[1][0] = b01;
+	map->y->coeff[0][1] = b11;
+
+	// Mask out the cross-terms
+	map->x->mask[1][1] = 0;
+	map->y->coeff[1][1] = 0;
+	// Set the cross-terms to zero, just in case...
+	map->x->coeff[1][1] = 0.0;
+	map->y->coeff[1][1] = 0.0;
+
+	// Plug it into the array
+	maps->data[i] = map;
+
+#ifdef TESTING
+	// Check psLib's answer is the same as I would get with old version
+	double detx = 256.0, dety = 256.0;
+	// My code
+	double myx = a0 + b00*detx + b10*dety;
+	double myy = a1 + b01*detx + b11*dety;
+	// psLib
+	psPlane *det = psAlloc(sizeof(psPlane));
+	det->x = 256.0;
+	det->y = 256.0;
+	psPlane *pslib = psPlaneTransformApply(NULL, map, det);
+
+	psTrace("stac.read.maps.test",0,"me: %lf,%lf\n",myx,myy);
+	psTrace("stac.read.maps.test",0,"ps: %lf,%lf\n",pslib->x,pslib->y);
+	psFree(det);
+	psFree(pslib);
+#endif
+
+    }
+    psTrace("stac.read.maps",1,"%d maps read.\n",nFiles);
+
+    return maps;
+}
Index: /trunk/stac/src/stacRejection.c
===================================================================
--- /trunk/stac/src/stacRejection.c	(revision 2500)
+++ /trunk/stac/src/stacRejection.c	(revision 2500)
@@ -0,0 +1,99 @@
+#include <stdio.h>
+#include "pslib.h"
+#include "stac.h"
+
+
+psArray *stacRejection(psArray *inputs,	// Input images
+		       psArray *transformed, // Transformed images
+		       psImage *combined, // Combined image
+		       psArray *maps,	// Maps from input to transformed image
+		       psArray *rejected, // Rejected images
+		       stacConfig *config // Configuration
+    )
+{
+    int nImages = inputs->n;		// Number of input images
+
+    // Check inputs
+    if (inputs->n != transformed->n) {
+	psError("stac.rejection", "Number of input (%d) and transformed (%d) images does not match.\n",
+		inputs->n, transformed->n);
+	return NULL;
+    }
+    for (int i = 0; i < nImages; i++) {
+	if ((((psImage*)(transformed->data[i]))->numRows != combined->numRows) ||
+	    (((psImage*)(transformed->data[i]))->numCols != combined->numCols)) {
+	    psError("stac.rejection",
+		    "Sizes of transformed (%dx%d) and combined (%dx%d) images do not match.\n",
+		    ((psImage*)(transformed->data[i]))->numCols, ((psImage*)(transformed->data[i]))->numRows,
+		    combined->numCols, combined->numRows);
+	    return NULL;
+	}
+	if ((((psImage*)(transformed->data[i]))->numRows != ((psImage*)(rejected->data[i]))->numRows) ||
+	    (((psImage*)(transformed->data[i]))->numCols != ((psImage*)(rejected->data[i]))->numCols)) {
+	    psError("stac.rejection",
+		    "Sizes of transformed (%dx%d) and rejected (%dx%d) images do not match.\n",
+		    ((psImage*)(transformed->data[i]))->numCols, ((psImage*)(transformed->data[i]))->numRows,
+		    ((psImage*)(rejected->data[i]))->numCols, ((psImage*)(rejected->data[i]))->numRows);
+	    return NULL;
+	}
+    }
+
+    // Size of output images
+    int nxOutput = combined->numCols;
+    int nyOutput = combined->numRows;
+
+    // Stuff for the transformations
+    psPlane *inCoords = psAlloc(sizeof(psPlane)); // Coordinates on the input
+    psPlane *outCoords = psAlloc(sizeof(psPlane)); // Coordinates on the output
+
+    psTrace("stac.rejection", 1, "Mapping rejection masks back to source....\n");
+
+    // Transform rejection masks back to source
+    psArray *inputRej = psArrayAlloc(nImages);
+    for (int i = 0; i < nImages; i++) {
+	// Size of input image
+	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 *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);
+
+	// Transform the mask
+	// Optimisation option is to only transform the pixels that have been rejected in the output,
+	// calculate derivatives of the map, and use that as a buffer around the transformed position
+	// in the input image.
+	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;
+		(void)psPlaneTransformApply(outCoords, map, inCoords);
+		mask->data.F32[y][x] = psImagePixelInterpolate(reject, outCoords->x, outCoords->y,
+							       NULL, 0, 0.0, PS_INTERPOLATE_BILINEAR);
+
+	    }
+	}
+
+#ifdef TESTING
+	// Write error image out to check
+	char maskfile[MAXCHAR];	// Filename of mask image
+	sprintf(maskfile,"%s.mask",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);
+    
+	// Stuff into the array
+	inputRej->data[i] = maskU8;
+
+    }
+
+    psFree(inCoords);
+    psFree(outCoords);
+
+    return inputRej;
+}
Index: /trunk/stac/src/stacTransform.c
===================================================================
--- /trunk/stac/src/stacTransform.c	(revision 2500)
+++ /trunk/stac/src/stacTransform.c	(revision 2500)
@@ -0,0 +1,242 @@
+#include <stdio.h>
+#include "pslib.h"
+#include "stac.h"
+
+#define SQUARE(x) ((x)*(x))
+#define MIN(x,y) (((x) > (y)) ? (y) : (x))
+#define MAX(x,y) (((x) > (y)) ? (x) : (y))
+
+#define MAXCHAR 80
+
+static int numTransforms = 0;		// Number of transformations performed
+
+
+// Hacked the original ps_ImagePixelInterpolateBILINEAR_F32 to add variances
+// i.e., to square the fractions when combining.
+inline psF64 p_psImageErrorInterpolateBILINEAR_F32(const psImage* input, 
+						   float x, 
+						   float y, 
+						   const psImage* mask, 
+						   unsigned int maskVal,
+						   psF64 unexposedValue)
+{ 
+    double floorX = floor((psF64)(x) - 0.5); 
+    double floorY = floor((psF64)(y) - 0.5); 
+    psF64 fracX = x - 0.5 - floorX; 
+    psF64 fracY = y - 0.5 - floorY; 
+    int intFloorX = (int) floorX; 
+    int intFloorY = (int) floorY; 
+    int lastX = input->numCols - 1; 
+    int lastY = input->numRows - 1; 
+    psF32 V00; 
+    psF32 V01; 
+    psF32 V10; 
+    psF32 V11; 
+    bool valid00 = false; 
+    bool valid01 = false; 
+    bool valid10 = false; 
+    bool valid11 = false; 
+    
+    if (intFloorY >= 0 && intFloorY <= lastY) { 
+        if (intFloorX >= 0 && intFloorX <= lastX) { 
+            V00 = input->data.F32[intFloorY][intFloorX]; 
+	    valid00 = true;
+        } 
+        if (intFloorX >= -1 && intFloorX < lastX) { 
+            V10 = input->data.F32[intFloorY][intFloorX+1]; 
+	    valid10 = true;
+        } 
+    } 
+    if (intFloorY >= -1 && intFloorY < lastY) { 
+        if (intFloorX >= 0 && intFloorX <= lastX) { 
+            V01 = input->data.F32[intFloorY+1][intFloorX]; 
+	    valid01 = true;
+        } 
+        if (intFloorX >= -1 && intFloorX < lastX) { 
+            V11 = input->data.F32[intFloorY+1][intFloorX+1]; 
+	    valid11 = true;
+        } 
+    } 
+    
+    /* cover likely case of all pixels being valid more efficiently */  
+    if (valid00 && valid10 && valid01 && valid11) { 
+        /* formula from the ADD */ 
+        return V00*SQUARE((1.0-fracX)*(1.0-fracY)) + V10*SQUARE(fracX*(1.0-fracY)) +
+	    V01*SQUARE(fracY*(1.0-fracX)) + V11*SQUARE(fracX*fracY);
+    } 
+    
+    /* OK, at least one pixel is not valid - need to do it piecemeal */ 
+    
+    psF64 V0; 
+    bool valid0 = true; 
+    if (valid00 && valid10) { 
+        V0 = V00*SQUARE(1-fracX)+V10*SQUARE(fracX); 
+    } else if (valid00) { 
+        V0 = V00; 
+    } else if (valid10) { 
+        V0 = V10; 
+    } else { 
+        valid0 = false; 
+    } 
+    
+    psF64 V1; 
+    bool valid1 = true; 
+    if (valid01 && valid11) { 
+        V1 = V01*SQUARE(1-fracX)+V11*SQUARE(fracX); 
+    } else if (valid01) { 
+        V1 = V01; 
+    } else if (valid11) { 
+        V1 = V11; 
+    } else { 
+        valid1 = false; 
+    } 
+    
+    if (valid0 && valid1) { 
+        return ( V0*SQUARE(1-fracY) + V1*SQUARE(fracY) ); 
+    } else if (valid0) { 
+        return V0; 
+    } else if (valid1) { 
+        return V1; 
+    } 
+    
+    return unexposedValue; 
+}
+
+
+
+psArray *stacTransform(const psArray *images, // Array of images to be transformed
+		       const psArray *maps, // Array of polynomials that do the transformation
+		       const psArray *errors, // Array of error images to be transformed
+		       psArray **outErrors, // Transformed error images for output
+		       const psArray *masks, // Masks of input images
+		       const stacConfig *config	// Configuration
+    )
+{
+    int nImages = images->n;		// Number of images
+    int nx = config->outnx, ny = config->outny; // Size of output images
+    numTransforms++;
+
+    // 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 NULL;
+    }
+    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 NULL;
+    }
+
+    if (errors && (*outErrors == NULL)) {
+	// Allocate the output error images, if required
+	*outErrors = psArrayAlloc(errors->n);
+    } 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 NULL;
+    }
+
+    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 NULL;
+	}
+	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 NULL;
+	    }
+	}
+    }
+
+    // Arrays of images for return
+    psArray *outImages = psArrayAlloc(nImages);	// Output images
+
+    // Stuff for the transformations
+    psPlane *detector = psAlloc(sizeof(psPlane)); // Coordinates on the detector
+    psPlane *sky = psAlloc(sizeof(psPlane)); // Coordinates on the sky
+
+    // Iterate over the images
+    for (int n = 0; n < nImages; n++) {
+	psTrace("stac.transform", 1, "Transforming image %d....\n",n);
+
+	// Pull out the various stuff we're working on
+	psImage *image = images->data[n]; // The input image
+	psPlaneTransform *map = maps->data[n]; // The map
+	psImage *error = errors->data[n]; // The error image
+
+	// Initialise the output images
+	psImage *outImage = psImageAlloc(nx, ny, PS_TYPE_F32);
+	psImage *outError = psImageAlloc(nx, ny, PS_TYPE_F32);
+	psTrace("stac.transform", 5, "Allocating space for transformed image, %dx%d\n", nx, ny);
+	for (int y = 0; y < ny; y++) {
+	    for (int x = 0; x < nx; x++) {
+		outImage->data.F32[y][x] = 0.0;
+		outError->data.F32[y][x] = 0.0;
+	    }
+	}
+
+	// Mask
+	psImage *mask = NULL;
+	if (masks != NULL) {
+	    mask = masks->data[n];
+	}
+
+#if 0
+	// Calculate scales; could be adapted for higher order than linear by moving this into the
+	// pixel-dependent section and calculating derivatives
+	double xscale = 0.5*sqrt(SQUARE(map->x->coeff[0][1]) + SQUARE(map->x->coeff[1][0]));
+	double yscale = 0.5*sqrt(SQUARE(map->y->coeff[0][1]) + SQUARE(map->y->coeff[1][0]));
+	double areascale = 0.25 / xscale / yscale;
+#endif
+
+	// Iterate over the output image pixels
+	for (int y = 0; y < ny; y++) {
+	    for (int x = 0; x < nx; x++) {
+		// Transform!
+		// Could no doubt make this quicker by calculating derivatives and incrementing by those.
+		sky->x = (double)x + 0.5;
+		sky->y = (double)y + 0.5;
+		(void)psPlaneTransformApply(detector, map, sky);
+
+		outImage->data.F32[y][x] = (psF32)psImagePixelInterpolate(image, detector->x, detector->y,
+									  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] = sqrtf(outError->data.F32[y][x]);
+	    }
+	} // Iterating over output pixels
+
+#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);
+	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);
+#endif
+
+	// Plug the new images into arrays
+	outImages->data[n] = outImage;
+	(*outErrors)->data[n] = outError;
+
+    } // Iterating over images
+
+    // Done with transformations
+    psFree(detector);
+    psFree(sky);
+
+    return outImages;
+}
+