Index: trunk/stac/Makefile.am
===================================================================
--- trunk/stac/Makefile.am	(revision 5745)
+++ trunk/stac/Makefile.am	(revision 5745)
@@ -0,0 +1,6 @@
+SUBDIRS = src/
+
+CLEANFILES = *~ core core.*
+
+EXTRA_DIST = \
+	autogen.sh
Index: trunk/stac/autogen.sh
===================================================================
--- trunk/stac/autogen.sh	(revision 5745)
+++ trunk/stac/autogen.sh	(revision 5745)
@@ -0,0 +1,103 @@
+#!/bin/sh
+# Run this to generate all the initial makefiles, etc.
+
+srcdir=`dirname $0`
+test -z "$srcdir" && srcdir=.
+
+ORIGDIR=`pwd`
+cd $srcdir
+
+PROJECT=stac
+TEST_TYPE=-f
+# change this to be a unique filename in the top level dir
+FILE=autogen.sh
+
+DIE=0
+
+LIBTOOLIZE=libtoolize
+ACLOCAL=aclocal
+AUTOHEADER=autoheader
+AUTOMAKE=automake
+AUTOCONF=autoconf
+
+#($LIBTOOLIZE --version) < /dev/null > /dev/null 2>&1 || {
+#        echo
+#        echo "You must have $LIBTOOlIZE installed to compile $PROJECT."
+#        echo "Download the appropriate package for your distribution,"
+#        echo "or get the source tarball at http://ftp.gnu.org/gnu/libtool/"
+#        DIE=1
+#}
+
+($ACLOCAL --version) < /dev/null > /dev/null 2>&1 || {
+        echo
+        echo "You must have $ACLOCAL installed to compile $PROJECT."
+        echo "Download the appropriate package for your distribution,"
+        echo "or get the source tarball at http://ftp.gnu.org/gnu/automake/"
+        DIE=1
+}
+
+#($AUTOHEADER --version) < /dev/null > /dev/null 2>&1 || {
+#        echo
+#        echo "You must have $AUTOHEADER installed to compile $PROJECT."
+#        echo "Download the appropriate package for your distribution,"
+#        echo "or get the source tarball at http://ftp.gnu.org/gnu/autoconf/"
+#        DIE=1
+#}
+
+($AUTOMAKE --version) < /dev/null > /dev/null 2>&1 || {
+        echo
+        echo "You must have $AUTOMAKE installed to compile $PROJECT."
+        echo "Download the appropriate package for your distribution,"
+        echo "or get the source tarball at http://ftp.gnu.org/gnu/automake/"
+        DIE=1
+}
+
+($AUTOCONF --version) < /dev/null > /dev/null 2>&1 || {
+        echo
+        echo "You must have $AUTOCONF installed to compile $PROJECT."
+        echo "Download the appropriate package for your distribution,"
+        echo "or get the source tarball at http://ftp.gnu.org/gnu/autoconf/"
+        DIE=1
+}
+
+if test "$DIE" -eq 1; then
+        exit 1
+fi
+
+test $TEST_TYPE $FILE || {
+        echo "You must run this script in the top-level $PROJECT directory"
+        exit 1
+}
+
+if test -z "$*"; then
+        echo "I am going to run ./configure with no arguments - if you wish "
+        echo "to pass any to it, please specify them on the $0 command line."
+fi
+
+#$LIBTOOLIZE --copy --force || echo "$LIBTOOlIZE failed"
+$ACLOCAL || echo "$ACLOCAL failed"
+#$AUTOHEADER || echo "$AUTOHEADER failed"
+$AUTOMAKE --add-missing --force-missing --copy || echo "$AUTOMAKE failed"
+$AUTOCONF || echo "$AUTOCONF failed"
+
+cd $ORIGDIR
+
+run_configure=true
+for arg in $*; do
+    case $arg in
+        --no-configure)
+            run_configure=false
+            ;;
+        *)
+            ;;
+    esac
+done
+
+if $run_configure; then
+    $srcdir/configure --enable-maintainer-mode "$@"
+    echo
+    echo "Now type 'make' to compile $PROJECT."
+else
+    echo
+    echo "Now run 'configure' and 'make' to compile $PROJECT."
+fi
Index: trunk/stac/configure.ac
===================================================================
--- trunk/stac/configure.ac	(revision 5745)
+++ trunk/stac/configure.ac	(revision 5745)
@@ -0,0 +1,26 @@
+AC_PREREQ(2.59)
+
+AC_INIT([stac], [0.0.1], [price@ifa.hawaii.edu])
+AC_CONFIG_SRCDIR([src])
+
+AM_INIT_AUTOMAKE([1.6 foreign dist-bzip2])
+dnl AM_CONFIG_HEADER([config.h])
+AM_MAINTAINER_MODE
+
+AC_LANG(C)
+AC_GNU_SOURCE
+AC_PROG_CC
+AC_PROG_INSTALL
+dnl AC_PROG_LIBTOOL
+
+PKG_CHECK_MODULES([PSLIB], [pslib >= 0.9.0]) 
+
+CFLAGS="${CFLAGS}"
+stac_CFLAGS="-Wall -Werror -std=c99 -O2 -DPS_NO_TRACE"
+AC_SUBST([stac_CFLAGS])
+
+AC_CONFIG_FILES([
+  Makefile
+  src/Makefile
+])
+AC_OUTPUT
Index: trunk/stac/src/Makefile
===================================================================
--- trunk/stac/src/Makefile	(revision 5743)
+++ 	(revision )
@@ -1,109 +1,0 @@
-SHELL = /bin/sh
-CC = gcc
-
-### 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 \
-	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 stacTime.o
-
-SHIFTSIZE = shiftSize.o stacWrite.o stacConfig.o stacRead.o stacCheckMemory.o stacSize.o
-
-FIXIMAGE = fixImage.o
-
-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)
-		$(CC) $(CFLAGS) -o $@ $(STAC) $(LDFLAGS) $(OPTFLAGS)
-
-shift:		$(SHIFT)
-		$(CC) $(CFLAGS) -o $@ $(SHIFT) $(LDFLAGS) $(OPTFLAGS)
-
-combine:	$(COMBINE)
-		$(CC) $(CFLAGS) -o $@ $(COMBINE) $(LDFLAGS) $(OPTFLAGS)
-
-shiftSize:	$(SHIFTSIZE)
-		$(CC) $(CFLAGS) -o $@ $(SHIFTSIZE) $(LDFLAGS) $(OPTFLAGS)
-
-calcGradient:	$(CALCGRADIENT)
-		$(CC) $(CFLAGS) -o $@ $(CALCGRADIENT) $(LDFLAGS) $(OPTFLAGS)
-
-sum:		$(SUM)
-		$(CC) $(CFLAGS) -o $@ $(SUM) $(LDFLAGS) $(OPTFLAGS)
-
-
-clean:
-		-$(RM) *.o gmon.* profile.txt
-
-empty:		clean
-		-$(RM) $(TARGETS) TAGS
-
-test:		stac test_0.fits test_1.fits test_2.fits test_3.fits
-		-$(RM) testout.fits
-		-$(RM) testout.fits.pre
-		-$(RM) chi2_*.fits
-		-$(RM) test_[0-3].fits.err
-		-$(RM) test_[0-3].fits.mask
-		-$(RM) test_[0-3].fits.shift.*
-		-$(RM) test_[0-3].fits.shiftrej
-		-$(RM) test_[0-3].fits.shifterr.*
-		-$(RM) test_[0-3].fits.grad
-		-$(RM) test_[0-3].fits.region
-		-$(RM) test_[0-3].fits.rejmap
-		-$(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-generate ; $(MAKE) test
-		-$(RM) $(TARGET)
-		$(MAKE) clean
-		export OPTFLAGS=-fbranch-use ; $(MAKE) $(TARGET)
-		-$(RM) *.da
-
-# Tags for emacs
-tags:
-		etags `find . -name \*.[ch] -print`
-
Index: trunk/stac/src/Makefile.am
===================================================================
--- trunk/stac/src/Makefile.am	(revision 5745)
+++ trunk/stac/src/Makefile.am	(revision 5745)
@@ -0,0 +1,107 @@
+# Compilation flags for Automake
+AM_CPPFLAGS = $(PSLIB_CFLAGS) $(stac_CFLAGS)
+AM_LDFLAGS = $(PSLIB_LIBS)
+
+# Programs to install in the "bin" directory
+bin_PROGRAMS = \
+	calcGradient \
+	combine \
+	shift \
+	shiftSize \
+	stac \
+	sum
+
+# Header files that don't need to be installed
+noinst_HEADERS = \
+	combineConfig.h \
+	stac.h \
+	stacConfig.h
+
+# List of sources for each of the programs
+stac_SOURCES = \
+	stac.c \
+	stacAreaOfInterest.c \
+	stacConfig.c \
+	stacCombine.c \
+	stacCheckMemory.c \
+	stacErrorImages.c \
+	stacHelp.c \
+	stacInvertMaps.c \
+	stacRead.c \
+	stacRejection.c \
+	stacScales.c \
+	stacSize.c \
+	stacTime.c \
+	stacTransform.c
+
+calcGradient_SOURCES = \
+	calcGradient.c \
+	stacRejection.c
+
+combine_SOURCES = \
+	combine.c \
+	combineConfig.c \
+	stacCheckMemory.c \
+	stacCombine.c \
+	stacErrorImages.c \
+	stacRead.c \
+	stacScales.c \
+	stacTime.c
+
+shift_SOURCES = \
+	shift.c \
+	stacCheckMemory.c \
+	stacInvertMaps.c \
+	stacRead.c \
+	stacSize.c \
+	stacTransform.c
+
+shiftSize_SOURCES = \
+	shiftSize.c \
+	stacCheckMemory.c \
+	stacConfig.c \
+	stacRead.c \
+	stacSize.c \
+	stacWrite.c
+
+sum_SOURCES = \
+	sum.c
+
+# Test programs
+TESTS = test
+
+# Clean target
+clean-local:
+	-$(RM) -f TAGS gmon.* profile.txt
+
+# Tags for emacs
+tags:
+	etags `find . -name \*.[ch] -print`
+
+test:	stac test_0.fits test_1.fits test_2.fits test_3.fits
+	-$(RM) testout.fits
+	-$(RM) testout.fits.pre
+	-$(RM) chi2_*.fits
+	-$(RM) test_[0-3].fits.err
+	-$(RM) test_[0-3].fits.mask
+	-$(RM) test_[0-3].fits.shift.*
+	-$(RM) test_[0-3].fits.shiftrej
+	-$(RM) test_[0-3].fits.shifterr.*
+	-$(RM) test_[0-3].fits.grad
+	-$(RM) test_[0-3].fits.region
+	-$(RM) test_[0-3].fits.rejmap
+	-$(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
+
+
Index: trunk/stac/src/stac.c
===================================================================
--- trunk/stac/src/stac.c	(revision 5743)
+++ trunk/stac/src/stac.c	(revision 5745)
@@ -7,5 +7,6 @@
 int main(int argc, char **argv)
 {
-    double startTime = getTime();
+    double startTime;
+    startTime = getTime();
 
 #if 0
@@ -155,5 +156,5 @@
         for (int i = 0; i < transformed->n; i++) {
             char shiftName[MAXCHAR];    // Filename of shifted image
-            sprintf(shiftName, "%s.shift", config->inputs->data[i]);
+            sprintf(shiftName, "%s.shift", (const char*)config->inputs->data[i]);
             psFits *shiftFile = psFitsOpen(shiftName, "w");
             image = psImageCopy(NULL, transformed->data[i], PS_TYPE_F32);
Index: trunk/stac/src/stacCheckMemory.c
===================================================================
--- trunk/stac/src/stacCheckMemory.c	(revision 5743)
+++ trunk/stac/src/stacCheckMemory.c	(revision 5745)
@@ -4,5 +4,5 @@
 
 
-#define LEAKS "leaks.dat"		// File to which to write leaks data
+#define LEAKS "leaks.dat"               // File to which to write leaks data
 
 
@@ -10,9 +10,9 @@
 {
     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));
+             "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;
 }
@@ -20,14 +20,14 @@
 
 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
+                       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));
+             "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));
 }
 
@@ -36,6 +36,6 @@
 void stacCheckMemory(void)
 {
-    psMemBlock **leaks = NULL;		// List of leaks
-    FILE *leakFile;			// File to write leaks to
+    psMemBlock **leaks = NULL;          // List of leaks
+    FILE *leakFile;                     // File to write leaks to
 
     psTrace("stac.checkMemory", 1, "Checking for memory problems....\n");
@@ -44,6 +44,6 @@
 
     if ((leakFile = fopen(LEAKS, "w")) == NULL) {
-	fprintf(stderr, "Unable to open leaks file, %s\n", LEAKS);
-	return;
+        fprintf(stderr, "Unable to open leaks file, %s\n", LEAKS);
+        return;
     }
 
@@ -51,11 +51,12 @@
     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);
+        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
+    int nCorrupted;                     // Number of corrupted memory blocks
+    nCorrupted = psMemCheckCorruption(false);
     psTrace("stac.checkMemory", 1, "%d memory blocks corrupted.\n", nCorrupted);
 
Index: trunk/stac/src/stacCombine.c
===================================================================
--- trunk/stac/src/stacCombine.c	(revision 5743)
+++ trunk/stac/src/stacCombine.c	(revision 5745)
@@ -8,7 +8,7 @@
 #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, 1 = don't use, 0 = use
+float stacCombineMean(psVector *values, // Values for which to take the mean
+                      psVector *errors, // Errors in the values
+                      psVector *masks   // Masks for the values, 1 = don't use, 0 = use
     )
 {
@@ -26,14 +26,14 @@
     int num = values->n;
     for (int i = 0; i < num; i++) {
-	if (! masks->data.U8[i]) {
-	    // "error" here is the variance
-	    sum += values->data.F32[i] / errors->data.F32[i];
-	    weights += 1.0 / errors->data.F32[i];
-	}
+        if (! masks->data.U8[i]) {
+            // "error" here is the variance
+            sum += values->data.F32[i] / errors->data.F32[i];
+            weights += 1.0 / errors->data.F32[i];
+        }
     }
     if (weights > 0.0) {
-	return (float)(sum/weights);
+        return (float)(sum/weights);
     } else {
-	return NAN;
+        return NAN;
     }
 #endif
@@ -43,6 +43,6 @@
 
 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
+                        psVector *errors, // Errors in the values
+                        psVector *masks // Masks for the values, 0 = don't use, 1 = use
     )
 {
@@ -56,13 +56,13 @@
 
 
-bool stacCombine(psImage **combined,	// The combined image for output
-		 psArray **rejected,	// Array of rejection masks
-		 psArray *images,	// Array of transformed images
-		 psArray *errors,	// Array of transformed error images
-		 int nReject,		// Number of rejection iterations
-		 psImage *region,	// Region to combine
-		 psVector *saturated,	// Saturation limits for each image
-		 psVector *bad,		// Bad pixel limits for each image
-		 float reject		// Rejection (k-sigma)
+bool stacCombine(psImage **combined,    // The combined image for output
+                 psArray **rejected,    // Array of rejection masks
+                 psArray *images,       // Array of transformed images
+                 psArray *errors,       // Array of transformed error images
+                 int nReject,           // Number of rejection iterations
+                 psImage *region,       // Region to combine
+                 psVector *saturated,   // Saturation limits for each image
+                 psVector *bad,         // Bad pixel limits for each image
+                 float reject           // Rejection (k-sigma)
     )
 {
@@ -76,25 +76,25 @@
     assert(bad->n == images->n);
 
-    int nImages = images->n;		// Number of images
-    int numRows = ((psImage*)images->data[0])->numRows;	// Image size
+    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
-    
+
     // 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
-
-	assert(image->numCols == numCols && image->numRows == numRows);
-	assert(error->numCols == numCols && error->numRows == numRows);
+        psImage *image = (psImage *)images->data[i]; // The image
+        psImage *error = (psImage *)errors->data[i]; // The error image
+
+        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));
+    assert(!*combined || ((*combined)->numRows == numRows && (*combined)->numCols == numCols));
     if (*combined == NULL) {
-	*combined = psImageAlloc(numCols, numRows, PS_TYPE_F32); // Combined image
+        *combined = psImageAlloc(numCols, numRows, PS_TYPE_F32); // Combined image
     }
 
     // Check area of interest
-    assert(!region || (region->numRows == numRows) && (region->numCols == numCols));
+    assert(!region || (region->numRows == numRows && region->numCols == numCols));
 
     psTrace("stac.combine", 1, "Combining images....\n");
@@ -106,20 +106,20 @@
     // Set up rejection masks
     if (nReject > 0) {
-	if (*rejected == NULL) {
-	    // Allocate the rejection masks, if required
-	    *rejected = psArrayAlloc(nImages);
-	} else {
-	    assert((*rejected)->n != nImages);
-	}
-
-	// Create and initialise rejection masks
-	for (int i = 0; i < nImages; i++) {
-	    (*rejected)->data[i] = psImageAlloc(numCols, numRows, PS_TYPE_U8);
-	    for (int r = 0; r < numRows; r++) {
-		for (int c = 0; c < numCols; c++) {
-		    ((psImage*)((*rejected)->data[i]))->data.U8[r][c] = 0;
-		}
-	    }
-	}
+        if (*rejected == NULL) {
+            // Allocate the rejection masks, if required
+            *rejected = psArrayAlloc(nImages);
+        } else {
+            assert((*rejected)->n != nImages);
+        }
+
+        // Create and initialise rejection masks
+        for (int i = 0; i < nImages; i++) {
+            (*rejected)->data[i] = psImageAlloc(numCols, numRows, PS_TYPE_U8);
+            for (int r = 0; r < numRows; r++) {
+                for (int c = 0; c < numCols; c++) {
+                    ((psImage*)((*rejected)->data[i]))->data.U8[r][c] = 0;
+                }
+            }
+        }
     }
 
@@ -127,72 +127,72 @@
     // chi^2 image
     psImage *chi2 = psImageAlloc(numCols, numRows, PS_TYPE_F32);
-    static int iteration = 0;		// Number of times function has been called
-#endif
-    
+    static int iteration = 0;           // Number of times function has been called
+#endif
+
     for (int y = 0; y < numRows; y++) {
-	for (int x = 0; x < numCols; x++) {
-
-	    // Only combine those pixels requested
-	    if (!region || (region && region->data.U8[y][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]; // This is the variance!
-		    pixels->data.F32[i] = pixel;
-		    deltas->data.F32[i] = delta;
-		    if ((pixel >= saturated->data.F32[i]) || (pixel <= bad->data.F32[i]) ||
-			(! isfinite(pixel)) || (! isfinite(delta))) {
-			mask->data.U8[i] = (psU8)1; // Don't use!
-		    } else {
-			mask->data.U8[i] = (psU8)0; // Use.
-		    }
-		}
-		
-		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;
+        for (int x = 0; x < numCols; x++) {
+
+            // Only combine those pixels requested
+            if (!region || (region && region->data.U8[y][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]; // This is the variance!
+                    pixels->data.F32[i] = pixel;
+                    deltas->data.F32[i] = delta;
+                    if ((pixel >= saturated->data.F32[i]) || (pixel <= bad->data.F32[i]) ||
+                        (! isfinite(pixel)) || (! isfinite(delta))) {
+                        mask->data.U8[i] = (psU8)1; // Don't use!
+                    } else {
+                        mask->data.U8[i] = (psU8)0; // Use.
+                    }
+                }
+
+                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
-		// Calculate chi^2
-		chi2->data.F32[y][x] = 0.0;
-		int numGoodPix = 0;
-		for (int i = 0; i < nImages; i++) {
-		    if (! mask->data.U8[i]) {
-			chi2->data.F32[y][x] += SQUARE(pixels->data.F32[i] - average) / deltas->data.F32[i];
-			numGoodPix++;
-		    }
-		}
-		chi2->data.F32[y][x] /= (float)numGoodPix;
-#endif
-		
-		// Rejection iterations
-		bool keepGoing = true;	// Keep going with rejection?
-		for (int rejNum = 0; (rejNum < nReject) && keepGoing; rejNum++) {
-		    float max = 0.0;	// Maximum deviation
-		    int maxIndex = -1;	// Index of the maximum deviation
-		    for (int i = 0; i < nImages; i++) {
-			if (!mask->data.U8[i] &&
-			    ((pixels->data.F32[i] - average) / sqrtf(deltas->data.F32[i]) > max)) {
-			    max = (pixels->data.F32[i] - average) / sqrtf(deltas->data.F32[i]);
-			    maxIndex = i;
-			}
-		    }
-		    // Reject the pixel with the maximum deviation
-		    if (max > reject) {
-			mask->data.U8[maxIndex] = 1;
-			((psImage*)((*rejected)->data[maxIndex]))->data.U8[y][x] += 1;
-			// Re-do combination following rejection
-			average = stacCombineMean(pixels, deltas, mask);
-		    } else {
-			keepGoing = false;
-		    }
-		}
-	    
-	    } // Pixels of interest
-
-	}
+                // Calculate chi^2
+                chi2->data.F32[y][x] = 0.0;
+                int numGoodPix = 0;
+                for (int i = 0; i < nImages; i++) {
+                    if (! mask->data.U8[i]) {
+                        chi2->data.F32[y][x] += SQUARE(pixels->data.F32[i] - average) / deltas->data.F32[i];
+                        numGoodPix++;
+                    }
+                }
+                chi2->data.F32[y][x] /= (float)numGoodPix;
+#endif
+
+                // Rejection iterations
+                bool keepGoing = true;  // Keep going with rejection?
+                for (int rejNum = 0; (rejNum < nReject) && keepGoing; rejNum++) {
+                    float max = 0.0;    // Maximum deviation
+                    int maxIndex = -1;  // Index of the maximum deviation
+                    for (int i = 0; i < nImages; i++) {
+                        if (!mask->data.U8[i] &&
+                            ((pixels->data.F32[i] - average) / sqrtf(deltas->data.F32[i]) > max)) {
+                            max = (pixels->data.F32[i] - average) / sqrtf(deltas->data.F32[i]);
+                            maxIndex = i;
+                        }
+                    }
+                    // Reject the pixel with the maximum deviation
+                    if (max > reject) {
+                        mask->data.U8[maxIndex] = 1;
+                        ((psImage*)((*rejected)->data[maxIndex]))->data.U8[y][x] += 1;
+                        // Re-do combination following rejection
+                        average = stacCombineMean(pixels, deltas, mask);
+                    } else {
+                        keepGoing = false;
+                    }
+                }
+
+            } // Pixels of interest
+
+        }
     } // Iterating over output pixels
 
@@ -200,9 +200,9 @@
     // Write chi^2 image
     iteration++;
-    char chiName[MAXCHAR];		// Filename of chi^2 image
+    char chiName[MAXCHAR];              // Filename of chi^2 image
     sprintf(chiName,"chi2_%d.fits",iteration);
     psFits *chiFile = psFitsAlloc(chiName);
     if (!psFitsWriteImage(chiFile, NULL, chi2 , 0)) {
-	psErrorStackPrint(stderr, "Unable to write image: %s\n", chiName);
+        psErrorStackPrint(stderr, "Unable to write image: %s\n", chiName);
     }
     psTrace("stac.combine", 1, "Chi^2 image written to %s\n", chiName);
Index: trunk/stac/src/stacInvertMaps.c
===================================================================
--- trunk/stac/src/stacInvertMaps.c	(revision 5743)
+++ trunk/stac/src/stacInvertMaps.c	(revision 5745)
@@ -9,8 +9,8 @@
 
 psArray *stacInvertMaps(const psArray *maps, // Array of maps to invert
-			int outnx, int outny // Size of output image
+                        int outnx, int outny // Size of output image
     )
 {
-    int nMaps = maps->n;		// Number of maps
+    int nMaps = maps->n;                // Number of maps
     psArray *inverted = psArrayAlloc(nMaps); // Array of inverted maps for output
 
@@ -23,176 +23,176 @@
     for (int mapNum = 0; mapNum < nMaps; mapNum++) {
 
-	psPlaneTransform *oldMap = (psPlaneTransform*)maps->data[mapNum]; // Uninverted map
-	// Check input
-	assert(oldMap->x->nX == oldMap->x->nY && oldMap->y->nX == oldMap->y->nY &&
-	       oldMap->x->nX == oldMap->y->nX);
-	int order = oldMap->x->nX + 1;	// Polynomial order
-	psTrace("stac.invertMaps", 4, "Generating order %d polynomial inverse transformation.\n", order);
-	psPlaneTransform *newMap = psPlaneTransformAlloc(order + 1, order + 1);	// Inverted map
-
-	// Create fake polynomial to use in evaluation
-	psPolynomial2D *fakePoly = psPolynomial2DAlloc(order, order, PS_POLYNOMIAL_ORD);
-	for (int i = 0; i < order; i++) {
-	    for (int j = 0; j < order; j++) {
-		fakePoly->coeff[i][j] = 1.0; // Set all coeffecients to 1
-		fakePoly->mask[i][j] = 1; // Mask all coefficients; unmask to evaluate
-	    }
-	}
-
-	// A grid of xin,yin --> xout,yout
-	psVector *xIn = psVectorAlloc(NUM_GRID * NUM_GRID, PS_TYPE_F32);
-	psVector *yIn = psVectorAlloc(NUM_GRID * NUM_GRID, PS_TYPE_F32);
-	psVector *xOut = psVectorAlloc(NUM_GRID * NUM_GRID, PS_TYPE_F32);
-	psVector *yOut = psVectorAlloc(NUM_GRID * NUM_GRID, PS_TYPE_F32);
-
-	// Create grid of points
-	for (int yint = 0; yint < NUM_GRID; yint++) {
-	    inCoord->y = (float)(yint * outny) / (float)(NUM_GRID - 1);
-	    for (int xint = 0; xint < NUM_GRID; xint++) {
-		inCoord->x = (float)(xint * outnx) / (float)(NUM_GRID - 1);
-
-		(void)psPlaneTransformApply(outCoord, oldMap, inCoord);
-
-		xOut->data.F32[yint*NUM_GRID + xint] = inCoord->x;
-		yOut->data.F32[yint*NUM_GRID + xint] = inCoord->y;
-		xIn->data.F32[yint*NUM_GRID + xint] = outCoord->x;
-		yIn->data.F32[yint*NUM_GRID + xint] = outCoord->y;
-	    }
-	}
-
-	// Initialise the matrix and vectors
-	int nCoeff = order * (order + 1) / 2; // Number of polynomial coefficients
-	psImage *matrix = psImageAlloc(nCoeff, nCoeff, PS_TYPE_F64); // Matrix for solution
-	psVector *xVector = psVectorAlloc(nCoeff, PS_TYPE_F64); // Vector for solution in x
-	psVector *yVector = psVectorAlloc(nCoeff, PS_TYPE_F64); // Vector for solution in y
-	for (int i = 0; i < nCoeff; i++) {
-	    for (int j = 0; j < nCoeff; j++) {
-		matrix->data.F64[i][j] = 0.0;
-	    }
-	    xVector->data.F64[i] = 0.0;
-	    yVector->data.F64[i] = 0.0;
-	}
-
-	// Iterate over the grid points
-	for (int g = 0; g < NUM_GRID*NUM_GRID; g++) {
-
-	    // Iterate over the polynomial coefficients, accumulating the matrix and vectors
-	    for (int i = 0, ijIndex = 0; i < order; i++) {
-		for (int j = 0; j < order - i; j++, ijIndex++) {
-
-		    fakePoly->mask[i][j] = 0;
-		    double ijPoly = psPolynomial2DEval(fakePoly, xIn->data.F32[g], yIn->data.F32[g]);
-		    fakePoly->mask[i][j] = 1;
-
-		    for (int m = 0, mnIndex = 0; m < order; m++) {
-			for (int n = 0; n < order - m; n++, mnIndex++) {
-			    
-			    fakePoly->mask[m][n] = 0;
-			    double mnPoly = psPolynomial2DEval(fakePoly, xIn->data.F32[g], yIn->data.F32[g]);
-			    fakePoly->mask[m][n] = 1;
-
-			    matrix->data.F64[ijIndex][mnIndex] += ijPoly * mnPoly;
-			}
-		    }
-		    
-		    xVector->data.F64[ijIndex] += ijPoly * (double)xOut->data.F32[g];
-		    yVector->data.F64[ijIndex] += ijPoly * (double)yOut->data.F32[g];
-		}
-	    } // Iterating over coefficients
-	} // Iterating over grid points
-
-	// 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
-	psVector *xSolution = psMatrixLUSolve(NULL, luMatrix, xVector, permutation); // Solution in x
-	psVector *ySolution = psMatrixLUSolve(NULL, luMatrix, yVector, permutation); // Solution in y
-
-	// Stuff coefficients into transformation
-	for (int i = 0, ijIndex = 0; i < order; i++) {
-	    for (int j = 0; j < order - i; j++, ijIndex++) {
-		newMap->x->coeff[i][j] = xSolution->data.F64[ijIndex];
-		newMap->y->coeff[i][j] = ySolution->data.F64[ijIndex];
-	    }
-	}
-	inverted->data[mapNum] = newMap;
+        psPlaneTransform *oldMap = (psPlaneTransform*)maps->data[mapNum]; // Uninverted map
+        // Check input
+        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);
+        psPlaneTransform *newMap = psPlaneTransformAlloc(order, order); // Inverted map
+
+        // Create fake polynomial to use in evaluation
+        psPolynomial2D *fakePoly = psPolynomial2DAlloc(order, order, PS_POLYNOMIAL_ORD);
+        for (int i = 0; i <= order; i++) {
+            for (int j = 0; j <= order; j++) {
+                fakePoly->coeff[i][j] = 1.0; // Set all coeffecients to 1
+                fakePoly->mask[i][j] = 1; // Mask all coefficients; unmask to evaluate
+            }
+        }
+
+        // A grid of xin,yin --> xout,yout
+        psVector *xIn = psVectorAlloc(NUM_GRID * NUM_GRID, PS_TYPE_F32);
+        psVector *yIn = psVectorAlloc(NUM_GRID * NUM_GRID, PS_TYPE_F32);
+        psVector *xOut = psVectorAlloc(NUM_GRID * NUM_GRID, PS_TYPE_F32);
+        psVector *yOut = psVectorAlloc(NUM_GRID * NUM_GRID, PS_TYPE_F32);
+
+        // Create grid of points
+        for (int yint = 0; yint < NUM_GRID; yint++) {
+            inCoord->y = (float)(yint * outny) / (float)(NUM_GRID - 1);
+            for (int xint = 0; xint < NUM_GRID; xint++) {
+                inCoord->x = (float)(xint * outnx) / (float)(NUM_GRID - 1);
+
+                (void)psPlaneTransformApply(outCoord, oldMap, inCoord);
+
+                xOut->data.F32[yint*NUM_GRID + xint] = inCoord->x;
+                yOut->data.F32[yint*NUM_GRID + xint] = inCoord->y;
+                xIn->data.F32[yint*NUM_GRID + xint] = outCoord->x;
+                yIn->data.F32[yint*NUM_GRID + xint] = outCoord->y;
+            }
+        }
+
+        // Initialise the matrix and vectors
+        int nCoeff = (order + 1) * (order + 2) / 2; // Number of polynomial coefficients
+        psImage *matrix = psImageAlloc(nCoeff, nCoeff, PS_TYPE_F64); // Matrix for solution
+        psVector *xVector = psVectorAlloc(nCoeff, PS_TYPE_F64); // Vector for solution in x
+        psVector *yVector = psVectorAlloc(nCoeff, PS_TYPE_F64); // Vector for solution in y
+        for (int i = 0; i < nCoeff; i++) {
+            for (int j = 0; j < nCoeff; j++) {
+                matrix->data.F64[i][j] = 0.0;
+            }
+            xVector->data.F64[i] = 0.0;
+            yVector->data.F64[i] = 0.0;
+        }
+
+        // Iterate over the grid points
+        for (int g = 0; g < NUM_GRID*NUM_GRID; g++) {
+
+            // Iterate over the polynomial coefficients, accumulating the matrix and vectors
+            for (int i = 0, ijIndex = 0; i <= order; i++) {
+                for (int j = 0; j <= order - i; j++, ijIndex++) {
+
+                    fakePoly->mask[i][j] = 0;
+                    double ijPoly = psPolynomial2DEval(fakePoly, xIn->data.F32[g], yIn->data.F32[g]);
+                    fakePoly->mask[i][j] = 1;
+
+                    for (int m = 0, mnIndex = 0; m <= order; m++) {
+                        for (int n = 0; n <= order - m; n++, mnIndex++) {
+
+                            fakePoly->mask[m][n] = 0;
+                            double mnPoly = psPolynomial2DEval(fakePoly, xIn->data.F32[g], yIn->data.F32[g]);
+                            fakePoly->mask[m][n] = 1;
+
+                            matrix->data.F64[ijIndex][mnIndex] += ijPoly * mnPoly;
+                        }
+                    }
+
+                    xVector->data.F64[ijIndex] += ijPoly * (double)xOut->data.F32[g];
+                    yVector->data.F64[ijIndex] += ijPoly * (double)yOut->data.F32[g];
+                }
+            } // Iterating over coefficients
+        } // Iterating over grid points
+
+        // 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
+        psVector *xSolution = psMatrixLUSolve(NULL, luMatrix, xVector, permutation); // Solution in x
+        psVector *ySolution = psMatrixLUSolve(NULL, luMatrix, yVector, permutation); // Solution in y
+
+        // Stuff coefficients into transformation
+        for (int i = 0, ijIndex = 0; i <= order; i++) {
+            for (int j = 0; j <= order - i; j++, ijIndex++) {
+                newMap->x->coeff[i][j] = xSolution->data.F64[ijIndex];
+                newMap->y->coeff[i][j] = ySolution->data.F64[ijIndex];
+            }
+        }
+        inverted->data[mapNum] = newMap;
 
 #ifdef TESTING
-	// Print x coefficients
-	psTrace("stac.invertMaps", 7, "x' = \n");
-	for (int i = 0; i < order; i++) {
-	    for (int j = 0; j < order - i; j++) {
-		psTrace("stac.invertMaps", 7, "      %f x^%d y^%d\n", newMap->x->coeff[i][j], i, j);
-	    }
-	}
-	// Print y coefficients
-	psTrace("stac.invertMaps", 7, "y' = \n");
-	for (int i = 0; i < order; i++) {
-	    for (int j = 0; j < order - i; j++) {
-		psTrace("stac.invertMaps", 7, "      %f x^%d y^%d\n", newMap->y->coeff[i][j], i, j);
-	    }
-	}
+        // Print x coefficients
+        psTrace("stac.invertMaps", 7, "x' = \n");
+        for (int i = 0; i <= order; i++) {
+            for (int j = 0; j <= order - i; j++) {
+                psTrace("stac.invertMaps", 7, "      %f x^%d y^%d\n", newMap->x->coeff[i][j], i, j);
+            }
+        }
+        // Print y coefficients
+        psTrace("stac.invertMaps", 7, "y' = \n");
+        for (int i = 0; i <= order; i++) {
+            for (int j = 0; j <= order - i; j++) {
+                psTrace("stac.invertMaps", 7, "      %f x^%d y^%d\n", newMap->y->coeff[i][j], i, j);
+            }
+        }
 #endif
 
 #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);
+        // 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
 
-	psFree(permutation);
-	psFree(luMatrix);
-	psFree(matrix);
-	psFree(xVector);
-	psFree(yVector);
-	psFree(xSolution);
-	psFree(ySolution);
-	psFree(fakePoly);
-	psFree(xIn);
-	psFree(yIn);
-	psFree(xOut);
-	psFree(yOut);
+        psFree(permutation);
+        psFree(luMatrix);
+        psFree(matrix);
+        psFree(xVector);
+        psFree(yVector);
+        psFree(xSolution);
+        psFree(ySolution);
+        psFree(fakePoly);
+        psFree(xIn);
+        psFree(yIn);
+        psFree(xOut);
+        psFree(yOut);
     }
 
     psFree(inCoord);
     psFree(outCoord);
-		
+
 
 #if 0
-	// Can't handle higher order than linear yet
-	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
-
-	// 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);
+        // Can't handle higher order than linear yet
+        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
+
+        // 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);
 #endif
 
Index: trunk/stac/src/stacRead.c
===================================================================
--- trunk/stac/src/stacRead.c	(revision 5743)
+++ trunk/stac/src/stacRead.c	(revision 5745)
@@ -124,8 +124,8 @@
     psTrace("stac.read.map", 5, "Polynomial order: %d\n", order);
 
-    psPlaneTransform *map = psPlaneTransformAlloc(order + 1, order + 1); // The transformation
+    psPlaneTransform *map = psPlaneTransformAlloc(order, order); // The transformation
     // Set coefficient masks
-    for (int i = 0; i < order + 1; i++) {
-        for (int j = 0; j < order + 1; j++) {
+    for (int i = 0; i <= order; i++) {
+        for (int j = 0; j <= order; j++) {
             if (i + j > order + 1) {
                 map->x->mask[i][j] = 1;
@@ -140,6 +140,6 @@
     // 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++) {
+    for (int k = 0; k <= order; k++) {
+        for (int j = 0; j <= k; j++) {
             int i = k - j;
             if (fscanf(mapfp, "%lf", &map->x->coeff[i][j]) != 1) {
@@ -154,6 +154,6 @@
     // 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++) {
+    for (int k = 0; k <= order; k++) {
+        for (int j = 0; j <= k; j++) {
             int i = k - j;
             if (fscanf(mapfp, "%lf", &map->y->coeff[i][j]) != 1) {
@@ -188,5 +188,5 @@
         }
         // Read the file
-        sprintf(mapfile,"%s.map",filenames->data[i]);
+        sprintf(mapfile, "%s.map", (const char*)filenames->data[i]);
         maps->data[i] = stacReadMap(mapfile);
         if (maps->data[i] == NULL) {
Index: trunk/stac/src/stacScales.c
===================================================================
--- trunk/stac/src/stacScales.c	(revision 5743)
+++ trunk/stac/src/stacScales.c	(revision 5745)
@@ -6,8 +6,8 @@
 #define SQUARE(x) ((x)*(x))
 
-#define SAMPLE 10			// Subsample rate for images
+#define SAMPLE 10                       // Subsample rate for images
 
 float stacBackground(const psImage *image, // Image for which to get the background
-		     int sample		// Sample in increments of this value
+                     int sample         // Sample in increments of this value
     )
 {
@@ -15,25 +15,25 @@
 
 #if 1
-    int size = image->numCols * image->numRows;	// Number of pixels in image
-    int numSamples = size / sample;	// Number of samples in image
+    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
-    int index = 0;			// Sample number
+    psVector *mask = psVectorAlloc(numSamples + 1, PS_TYPE_U8); // Mask for sample
+
+    int offset = 0;                     // Offset from start of the row
+    int index = 0;                      // Sample number
     for (int row = 0; row < image->numRows; row++) {
-	// I'd cast this as a "for", but this makes it a bit easier to understand.
-	int col = offset;
-	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++;
-	}
-	offset = col - image->numCols;
+        // I'd cast this as a "for", but this makes it a bit easier to understand.
+        int col = offset;
+        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++;
+        }
+        offset = col - image->numCols;
     }
 
@@ -58,11 +58,11 @@
 
 
-bool stacScales(psVector **scalesPtr,	// Scales to return
-		psVector **offsetsPtr,	// Offsets to return
-		const psArray *images,	// Images on which to measure the scales and offsets
-		const char *starFile,	// File containing coordinates to photometer
-		const char *starMapFile, // Map for coodinates to the common output frame
-		float xMapDiff, float yMapDiff, // Difference from the map to apply (due to size difference)
-		float aper		// Aperture to use for photometry (radius)
+bool stacScales(psVector **scalesPtr,   // Scales to return
+                psVector **offsetsPtr,  // Offsets to return
+                const psArray *images,  // Images on which to measure the scales and offsets
+                const char *starFile,   // File containing coordinates to photometer
+                const char *starMapFile, // Map for coodinates to the common output frame
+                float xMapDiff, float yMapDiff, // Difference from the map to apply (due to size difference)
+                float aper              // Aperture to use for photometry (radius)
     )
 {
@@ -71,6 +71,6 @@
     assert(images);
     for (int i = 0; i < images->n; i++) {
-	psImage *image = images->data[i];
-	assert(image->type.type == PS_TYPE_F32);
+        psImage *image = images->data[i];
+        assert(image->type.type == PS_TYPE_F32);
     }
 
@@ -78,20 +78,20 @@
     psVector *offsets = NULL; // Offsets between images (ADU)
     if (*scalesPtr) {
-	scales = *scalesPtr;
-	assert(scales);
-	assert(scales->n == images->n);
-	assert(scales->type.type == PS_TYPE_F32);
+        scales = *scalesPtr;
+        assert(scales);
+        assert(scales->n == images->n);
+        assert(scales->type.type == PS_TYPE_F32);
     } else {
-	scales = psVectorAlloc(images->n, PS_TYPE_F32);
-	*scalesPtr = scales;
+        scales = psVectorAlloc(images->n, PS_TYPE_F32);
+        *scalesPtr = scales;
     }
     if (*offsetsPtr) {
-	offsets = *offsetsPtr;
-	assert(offsets);
-	assert(offsets->n == images->n);
-	assert(offsets->type.type == PS_TYPE_F32);
+        offsets = *offsetsPtr;
+        assert(offsets);
+        assert(offsets->n == images->n);
+        assert(offsets->type.type == PS_TYPE_F32);
     } else {
-	offsets = psVectorAlloc(images->n, PS_TYPE_F32);
-	*offsetsPtr = offsets;
+        offsets = psVectorAlloc(images->n, PS_TYPE_F32);
+        *offsetsPtr = offsets;
     }
 
@@ -100,140 +100,138 @@
     double startTime = getTime();
     for (int i = 0; i < images->n; i++) {
-	offsets->data.F32[i] = stacBackground(images->data[i], SAMPLE);
-	psTrace("stac.scales", 5, "Background in image %d is %f\n", i, offsets->data.F32[i]);
-	double time = getTime();
-	psTrace("stac.scales", 10, "Took %f sec\n", time - startTime);
-	startTime = time;
+        offsets->data.F32[i] = stacBackground(images->data[i], SAMPLE);
+        psTrace("stac.scales", 5, "Background in image %d is %f\n", i, offsets->data.F32[i]);
+        double time = getTime();
+        psTrace("stac.scales", 10, "Took %f sec\n", time - startTime);
+        startTime = time;
     }
 
     // Now the scales
     if (starFile == NULL || starMapFile == NULL) {
-	psLogMsg("stac.scales", PS_LOG_INFO,
-		 "No coordinates available to set scales --- assuming all are identical.\n");
-	for (int i = 0; i < images->n; i++) {
-	    scales->data.F32[i] = 1.0;
-	    psTrace("stac.scales", 5, "Scale for image %d is %f\n", i, scales->data.F32[i]);
-	}
+        psLogMsg("stac.scales", PS_LOG_INFO,
+                 "No coordinates available to set scales --- assuming all are identical.\n");
+        for (int i = 0; i < images->n; i++) {
+            scales->data.F32[i] = 1.0;
+            psTrace("stac.scales", 5, "Scale for image %d is %f\n", i, scales->data.F32[i]);
+        }
     } else {
-	// Read star coordinates and map
-	psArray *starCoords = stacReadCoords(starFile); // Array of star coordinates
-	psPlaneTransform *starMap = stacReadMap(starMapFile); // Transformation for star coordinates
-
-	// Transform the stellar positions to match the transformed reference frame
-	psArray *starCoordsTransformed = psArrayAlloc(starCoords->n); // Transformed positions
-	// Fix up difference between map and output frame
-	starMap->x->coeff[0][0] -= xMapDiff;
-	starMap->y->coeff[0][0] -= yMapDiff;
-	for (int i = 0; i < starCoords->n; i++) {
-	    starCoordsTransformed->data[i] = psPlaneTransformApply(NULL, starMap, starCoords->data[i]);
-	}
-	psFree(starMap);
-
-	psArray *stars = psArrayAlloc(images->n); // Array of stellar photometry vectors
-	psArray *masks = psArrayAlloc(images->n); // Array of masks for stars
-
-	// Set scales relative to the first image
-	scales->data.F32[0] = 1.0;
-
-	// Iterate over images
-	for (int i = 0; i < scales->n; i++) {
-	    psImage *image = images->data[i]; // The image we're working with
-
-	    float background = offsets->data.F32[i]; // Background in image
-
-	    // Do photometry on transformed image i
-	    psTrace("stac.scales", 3, "Doing photometry on image %d....\n", i);
-	    psVector *photometry = psVectorAlloc(starCoords->n, PS_TYPE_F32); // Photometry of the stars
-	    psVector *mask = psVectorAlloc(starCoords->n, PS_TYPE_U8); // Mask for the photometry
-	    for (int j = 0; j < starCoordsTransformed->n; j++) {
-		psPlane *coords = starCoordsTransformed->data[j]; // The coordinates of the star
-		
-		if (coords->x < aper || coords->y < aper ||
-		    coords->x + aper > image->numCols - 1 ||
-		    coords->y + aper > image->numRows) {
-		    mask->data.U8[j] = 1;
-		} else {
-		    // Sum flux within the aperture
-		    float sum = 0.0;
-		    int numPix = 0;
-		    float aper2 = SQUARE(aper);
-		    for (int y = (int)floorf(coords->y - aper);
-			 y <= (int)ceilf(coords->y + aper); y++) {
-			for (int x = (int)floorf(coords->x - aper);
-			     x <= (int)ceilf(coords->x + aper); x++) {
-			    if (SQUARE((float)x + 0.5 - coords->x) + SQUARE((float)y + 0.5 - coords->y) <=
-				aper2) {
-				sum += image->data.F32[y][x];
-				numPix++;
-			    }
-			}
-		    }
-		    // 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;
-		    }
-		}
-	    }
-	    stars->data[i] = photometry;
-	    masks->data[i] = mask;
-	}
-	psFree(starCoords);
-	psFree(starCoordsTransformed);
-
-	// Get the scales
-	psVector *ref = stars->data[0];	// The reference photometry
-	psVector *refMask = masks->data[0]; // The reference mask
+        // Read star coordinates and map
+        psArray *starCoords = stacReadCoords(starFile); // Array of star coordinates
+        psPlaneTransform *starMap = stacReadMap(starMapFile); // Transformation for star coordinates
+
+        // Transform the stellar positions to match the transformed reference frame
+        psArray *starCoordsTransformed = psArrayAlloc(starCoords->n); // Transformed positions
+        // Fix up difference between map and output frame
+        starMap->x->coeff[0][0] -= xMapDiff;
+        starMap->y->coeff[0][0] -= yMapDiff;
+        for (int i = 0; i < starCoords->n; i++) {
+            starCoordsTransformed->data[i] = psPlaneTransformApply(NULL, starMap, starCoords->data[i]);
+        }
+        psFree(starMap);
+
+        psArray *stars = psArrayAlloc(images->n); // Array of stellar photometry vectors
+        psArray *masks = psArrayAlloc(images->n); // Array of masks for stars
+
+        // Set scales relative to the first image
+        scales->data.F32[0] = 1.0;
+
+        // Iterate over images
+        for (int i = 0; i < scales->n; i++) {
+            psImage *image = images->data[i]; // The image we're working with
+
+            // Do photometry on transformed image i
+            psTrace("stac.scales", 3, "Doing photometry on image %d....\n", i);
+            psVector *photometry = psVectorAlloc(starCoords->n, PS_TYPE_F32); // Photometry of the stars
+            psVector *mask = psVectorAlloc(starCoords->n, PS_TYPE_U8); // Mask for the photometry
+            for (int j = 0; j < starCoordsTransformed->n; j++) {
+                psPlane *coords = starCoordsTransformed->data[j]; // The coordinates of the star
+
+                if (coords->x < aper || coords->y < aper ||
+                    coords->x + aper > image->numCols - 1 ||
+                    coords->y + aper > image->numRows) {
+                    mask->data.U8[j] = 1;
+                } else {
+                    // Sum flux within the aperture
+                    float sum = 0.0;
+                    int numPix = 0;
+                    float aper2 = SQUARE(aper);
+                    for (int y = (int)floorf(coords->y - aper);
+                         y <= (int)ceilf(coords->y + aper); y++) {
+                        for (int x = (int)floorf(coords->x - aper);
+                             x <= (int)ceilf(coords->x + aper); x++) {
+                            if (SQUARE((float)x + 0.5 - coords->x) + SQUARE((float)y + 0.5 - coords->y) <=
+                                aper2) {
+                                sum += image->data.F32[y][x];
+                                numPix++;
+                            }
+                        }
+                    }
+                    // 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;
+                    }
+                }
+            }
+            stars->data[i] = photometry;
+            masks->data[i] = mask;
+        }
+        psFree(starCoords);
+        psFree(starCoordsTransformed);
+
+        // Get the scales
+        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;
+        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 *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);
+        psStats *stats = psStatsAlloc(PS_STAT_SAMPLE_MEDIAN); // Statistics
+#endif
+        for (int i = 1; i < scales->n; i++) {
+            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);
+            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;
+            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);
-
-	psFree(stars);
-	psFree(masks);
+            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);
+
+        psFree(stars);
+        psFree(masks);
     }
 
@@ -242,9 +240,9 @@
 
 
-bool stacRescale(psArray *images,	// Images to rescale
-		 psArray *errImages,	// Variance images to rescale
-		 const psImage *mask,	// Mask for pixels to scale
-		 const psVector *scales,// Scales for images
-		 const psVector *offsets // Offsets for images
+bool stacRescale(psArray *images,       // Images to rescale
+                 psArray *errImages,    // Variance images to rescale
+                 const psImage *mask,   // Mask for pixels to scale
+                 const psVector *scales,// Scales for images
+                 const psVector *offsets // Offsets for images
     )
 {
@@ -258,36 +256,36 @@
     assert(offsets->type.type == PS_TYPE_F32);
     for (int i = 0; i < images->n; i++) {
-	psImage *image = images->data[i]; // Image of interest
-	assert(image->type.type == PS_TYPE_F32);
-	assert(scales->data.F32[i] != 0);
-	if (mask) {
-	    assert(mask->type.type == PS_TYPE_U8);
-	    assert(image->numCols == mask->numCols && image->numRows == mask->numRows);
-	}
-	if (errImages) {
-	    psImage *errImage = errImages->data[i];
-	    assert(errImage->type.type == PS_TYPE_F32);
-	    assert(errImage->numCols == image->numCols && errImage->numRows == image->numRows);
-	}
-    }
-
-    for (int i = 0; i < images->n; i++) {
-	psImage *image = images->data[i]; // Image to rescale
-	psImage *errImage = NULL;	// Variance image to rescale
-	if (errImages) {
-	    errImage = errImages->data[i];
-	}
-	float scale = scales->data.F32[i]; // Scale to use
-	float offset = offsets->data.F32[i]; // Offset to use
-	for (int y = 0; y < image->numRows; y++) {
-	    for (int x = 0; x < image->numCols; x++) {
-		if (!mask || mask->data.F32[y][x]) {
-		    image->data.F32[y][x] = (image->data.F32[y][x] - offset) / scale;
-		    if (errImage) {
-			errImage->data.F32[y][x] /= SQUARE(scale);
-		    }
-		}
-	    }
-	}
+        psImage *image = images->data[i]; // Image of interest
+        assert(image->type.type == PS_TYPE_F32);
+        assert(scales->data.F32[i] != 0);
+        if (mask) {
+            assert(mask->type.type == PS_TYPE_U8);
+            assert(image->numCols == mask->numCols && image->numRows == mask->numRows);
+        }
+        if (errImages) {
+            psImage *errImage = errImages->data[i];
+            assert(errImage->type.type == PS_TYPE_F32);
+            assert(errImage->numCols == image->numCols && errImage->numRows == image->numRows);
+        }
+    }
+
+    for (int i = 0; i < images->n; i++) {
+        psImage *image = images->data[i]; // Image to rescale
+        psImage *errImage = NULL;       // Variance image to rescale
+        if (errImages) {
+            errImage = errImages->data[i];
+        }
+        float scale = scales->data.F32[i]; // Scale to use
+        float offset = offsets->data.F32[i]; // Offset to use
+        for (int y = 0; y < image->numRows; y++) {
+            for (int x = 0; x < image->numCols; x++) {
+                if (!mask || mask->data.F32[y][x]) {
+                    image->data.F32[y][x] = (image->data.F32[y][x] - offset) / scale;
+                    if (errImage) {
+                        errImage->data.F32[y][x] /= SQUARE(scale);
+                    }
+                }
+            }
+        }
     }
 
Index: trunk/stac/src/stacTransform.c
===================================================================
--- trunk/stac/src/stacTransform.c	(revision 5743)
+++ trunk/stac/src/stacTransform.c	(revision 5745)
@@ -10,108 +10,108 @@
 // 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)
+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 */ 
+    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 = 0.0;
+    psF32 V01 = 0.0;
+    psF32 V10 = 0.0;
+    psF32 V11 = 0.0;
+    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; 
+            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 = 0.0;
+    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 = 0.0;
+    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;
 }
 
 
 
-bool stacTransform(psArray **outputs,	// Transformed images for output
-		   psArray **outErrors, // Transformed error images for output
-		   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
-		   const psArray *masks, // Masks of input images
-		   const psImage *region, // Region of interest for transformation
-		   const psVector *scales, // Relative scales
-		   const psVector *offsets, // Relative offsets
-		   int outnx, int outny	// Size of output images
+bool stacTransform(psArray **outputs,   // Transformed images for output
+                   psArray **outErrors, // Transformed error images for output
+                   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
+                   const psArray *masks, // Masks of input images
+                   const psImage *region, // Region of interest for transformation
+                   const psVector *scales, // Relative scales
+                   const psVector *offsets, // Relative offsets
+                   int outnx, int outny // Size of output images
     )
 {
-    int nImages = images->n;		// Number of images
+    int nImages = images->n;            // Number of images
 
     // Check input sizes
@@ -126,10 +126,10 @@
     assert(!*outputs || (*outputs)->n == nImages);
     if (*outputs == NULL) {
-	*outputs = psArrayAlloc(nImages);
-	psTrace("stac.transform", 5, "Allocating space for transformed images, %dx%d\n", outnx, outny);
-	for (int i = 0; i < nImages; i++) {
-	    (*outputs)->data[i] = psImageAlloc(outnx, outny, PS_TYPE_F32);
-	    psImageInit((*outputs)->data[i], 0.0);
-	}
+        *outputs = psArrayAlloc(nImages);
+        psTrace("stac.transform", 5, "Allocating space for transformed images, %dx%d\n", outnx, outny);
+        for (int i = 0; i < nImages; i++) {
+            (*outputs)->data[i] = psImageAlloc(outnx, outny, PS_TYPE_F32);
+            psImageInit((*outputs)->data[i], 0.0);
+        }
     }
 
@@ -137,9 +137,9 @@
     assert(!errors || ! *outErrors || errors->n == (*outErrors)->n);
     if (errors && (*outErrors == NULL)) {
-	*outErrors = psArrayAlloc(errors->n);
-	psTrace("stac.transform", 5, "Allocating space for transformed error images, %dx%d\n", outnx, outny);
-	for (int i = 0; i < nImages; i++) {
-	    (*outErrors)->data[i] = psImageAlloc(outnx, outny, PS_TYPE_F32);
-	}
+        *outErrors = psArrayAlloc(errors->n);
+        psTrace("stac.transform", 5, "Allocating space for transformed error images, %dx%d\n", outnx, outny);
+        for (int i = 0; i < nImages; i++) {
+            (*outErrors)->data[i] = psImageAlloc(outnx, outny, PS_TYPE_F32);
+        }
     }
 
@@ -147,9 +147,9 @@
     assert(!masks || masks->n == nImages);
     if (masks != NULL) {
-	for (int i = 0; i < nImages; i++) {
-	    psImage *image = images->data[i];
-	    psImage *mask = masks->data[i];
-	    assert(mask->numRows == image->numRows && mask->numCols == image->numCols);
-	}
+        for (int i = 0; i < nImages; i++) {
+            psImage *image = images->data[i];
+            psImage *mask = masks->data[i];
+            assert(mask->numRows == image->numRows && mask->numCols == image->numCols);
+        }
     }
 
@@ -162,60 +162,60 @@
     // 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 *outImage = (*outputs)->data[n]; // The output image
-	psImage *error = NULL; // The error image
-	psImage *outError = NULL; // The output error image
-	if (errors) {
-	    error = errors->data[n];
-	    outError = (*outErrors)->data[n];
-	}
-	float offset = 0.0;		// Relative offset
-	float scale = 1.0;		// Relative scale
-	if (offsets) {
-	    offset = offsets->data.F32[n];
-	}
-	if (scales) {
-	    scale = scales->data.F32[n];
-	}
-
-	// Mask
-	psImage *mask = NULL;
-	if (masks != NULL) {
-	    mask = masks->data[n];
-	}
-
-	// Iterate over the output image pixels
-	for (int y = 0; y < outny; y++) {
-	    for (int x = 0; x < outnx; x++) {
-		// Only transform those pixels requested
-		if (!region || (region && region->data.U8[y][x])) {
-		    // Transform!
-		    sky->x = (double)x + 0.5;
-		    sky->y = (double)y + 0.5;
-		    (void)psPlaneTransformApply(detector, map, sky);
-		    
-		    // Change PS_INTERPOLATE_BILINEAR to best available technique.
-		    outImage->data.F32[y][x] = (psF32)psImagePixelInterpolate(image, detector->x,
-									      detector->y, mask, 1, NAN,
-									      PS_INTERPOLATE_BILINEAR);
-		    if (error) {
-			// Error is actually the variance
-			outError->data.F32[y][x] = (psF32)p_psImageErrorInterpolateBILINEAR_F32(error,
-												detector->x,
-												detector->y,
-												mask, 1, NAN);
-		    }
-
-		    outImage->data.F32[y][x] = (outImage->data.F32[y][x] - offset) / scale;
-		    outImage->data.F32[y][x] = outImage->data.F32[y][x] / SQUARE(scale);
-
-		} // Pixels of interest
-
-	    }
-	} // Iterating over output pixels
+        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 *outImage = (*outputs)->data[n]; // The output image
+        psImage *error = NULL; // The error image
+        psImage *outError = NULL; // The output error image
+        if (errors) {
+            error = errors->data[n];
+            outError = (*outErrors)->data[n];
+        }
+        float offset = 0.0;             // Relative offset
+        float scale = 1.0;              // Relative scale
+        if (offsets) {
+            offset = offsets->data.F32[n];
+        }
+        if (scales) {
+            scale = scales->data.F32[n];
+        }
+
+        // Mask
+        psImage *mask = NULL;
+        if (masks != NULL) {
+            mask = masks->data[n];
+        }
+
+        // Iterate over the output image pixels
+        for (int y = 0; y < outny; y++) {
+            for (int x = 0; x < outnx; x++) {
+                // Only transform those pixels requested
+                if (!region || (region && region->data.U8[y][x])) {
+                    // Transform!
+                    sky->x = (double)x + 0.5;
+                    sky->y = (double)y + 0.5;
+                    (void)psPlaneTransformApply(detector, map, sky);
+
+                    // Change PS_INTERPOLATE_BILINEAR to best available technique.
+                    outImage->data.F32[y][x] = (psF32)psImagePixelInterpolate(image, detector->x,
+                                                                              detector->y, mask, 1, NAN,
+                                                                              PS_INTERPOLATE_BILINEAR);
+                    if (error) {
+                        // Error is actually the variance
+                        outError->data.F32[y][x] = (psF32)p_psImageErrorInterpolateBILINEAR_F32(error,
+                                                                                                detector->x,
+                                                                                                detector->y,
+                                                                                                mask, 1, NAN);
+                    }
+
+                    outImage->data.F32[y][x] = (outImage->data.F32[y][x] - offset) / scale;
+                    outImage->data.F32[y][x] = outImage->data.F32[y][x] / SQUARE(scale);
+
+                } // Pixels of interest
+
+            }
+        } // Iterating over output pixels
 
     } // Iterating over images
Index: trunk/stac/src/stacWrite.c
===================================================================
--- trunk/stac/src/stacWrite.c	(revision 5743)
+++ trunk/stac/src/stacWrite.c	(revision 5745)
@@ -4,6 +4,6 @@
 #include "stac.h"
 
-bool stacWriteMap(const char *mapName,	// Filename to write to
-		  psPlaneTransform *map	// Map to write
+bool stacWriteMap(const char *mapName,  // Filename to write to
+                  psPlaneTransform *map // Map to write
     )
 {
@@ -12,42 +12,42 @@
     FILE *mapFile = fopen(mapName, "w");
     if (!mapFile) {
-	fprintf(stderr, "Unable to open map file: %s\n", mapName);
-	return false;
+        fprintf(stderr, "Unable to open map file: %s\n", mapName);
+        return false;
     }
 
-    psPolynomial2D *xMap = map->x;	// x transform
-    psPolynomial2D *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;	// The polynomial order
+    int order = xMap->nX;       // The polynomial order
     fprintf(mapFile, "%d\n", order);
-    
+
     // x coefficients
-    for (int k = 0; k < order + 1; k++) {
-	for (int j = 0; j < k + 1; j++) {
-	    int i = k - j;
-	    if (xMap->mask[i][j]) {
-		fprintf(mapFile, "0.0 ");
-	    } else {
-		fprintf(mapFile, "%g ", xMap->coeff[i][j]);
-	    }
-	}
+    for (int k = 0; k <= order; k++) {
+        for (int j = 0; j <= k; j++) {
+            int i = k - j;
+            if (xMap->mask[i][j]) {
+                fprintf(mapFile, "0.0 ");
+            } else {
+                fprintf(mapFile, "%g ", xMap->coeff[i][j]);
+            }
+        }
     }
     fprintf(mapFile, "\n");
 
     // y coefficients
-    for (int k = 0; k < order + 1; k++) {
-	for (int j = 0; j < k + 1; j++) {
-	    int i = k - j;
-	    if (yMap->mask[i][j]) {
-		fprintf(mapFile, "0.0 ");
-	    } else {
-		fprintf(mapFile, "%g ", yMap->coeff[i][j]);
-	    }
-	}
+    for (int k = 0; k <= order; k++) {
+        for (int j = 0; j <= k; j++) {
+            int i = k - j;
+            if (yMap->mask[i][j]) {
+                fprintf(mapFile, "0.0 ");
+            } else {
+                fprintf(mapFile, "%g ", yMap->coeff[i][j]);
+            }
+        }
     }
     fprintf(mapFile, "\n");
-    
+
     fclose(mapFile);
 
@@ -57,5 +57,5 @@
 
 bool stacWriteMaps(const psArray *names, // Filenames of the input images (will add ".map")
-		   const psArray *maps	// Maps to write
+                   const psArray *maps  // Maps to write
     )
 {
@@ -65,9 +65,9 @@
 
     for (int i = 0; i < names->n; i++) {
- 	char mapName[MAXCHAR];		// Filename of error image
-	sprintf(mapName, "%s.map", names->data[i]);
-	if (!stacWriteMap(mapName, maps->data[i])) {
-	    return false;
-	}
+        char mapName[MAXCHAR];          // Filename of error image
+        sprintf(mapName, "%s.map", (const char*)names->data[i]);
+        if (!stacWriteMap(mapName, maps->data[i])) {
+            return false;
+        }
     }
 
Index: trunk/stac/src/sum.c
===================================================================
--- trunk/stac/src/sum.c	(revision 5743)
+++ trunk/stac/src/sum.c	(revision 5745)
@@ -8,5 +8,4 @@
 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);