Index: /trunk/stac/configure.ac
===================================================================
--- /trunk/stac/configure.ac	(revision 6886)
+++ /trunk/stac/configure.ac	(revision 6887)
@@ -26,4 +26,5 @@
 
 stac_CFLAGS="-Wall -Werror -std=c99 -DPS_NO_TRACE"
+#stac_CFLAGS="-Wall -Werror -std=c99"
 AC_SUBST([stac_CFLAGS])
 
Index: /trunk/stac/src/combine.c
===================================================================
--- /trunk/stac/src/combine.c	(revision 6886)
+++ /trunk/stac/src/combine.c	(revision 6887)
@@ -25,4 +25,5 @@
     // Make fake errors
     psArray *errors = psArrayAlloc(images->n);
+    errors->n = images->n;
     for (int i = 0; i < images->n; i++) {
         psImage *image = images->data[i];
@@ -44,4 +45,6 @@
     psVector *saturated = psVectorAlloc(images->n, PS_TYPE_F32); // Saturation limits
     psVector *bad = psVectorAlloc(images->n, PS_TYPE_F32); // Bad limits
+    saturated->n = images->n;
+    bad->n = images->n;
     for (int i = 0; i < images->n; i++) {
         saturated->data.F32[i] = (config->saturated - offsets->data.F32[i]) / scales->data.F32[i];
Index: /trunk/stac/src/combineConfig.c
===================================================================
--- /trunk/stac/src/combineConfig.c	(revision 6886)
+++ /trunk/stac/src/combineConfig.c	(revision 6887)
@@ -11,20 +11,20 @@
 {
     fprintf (stderr, "shift: shift an image, given the transformation\n"
-	     "Usage: %s [-h] [-v] [-g GAIN] [-r READNOISE] [-s SAT] [-b BAD] [-p FILE MAP] [-a APER] [-k SIGMAREJ] [-n NREJECT] OUT IN1 IN2...\n"
-	     "where\n"
-	     "\t-h           Help (this info)\n"
-	     "\t-v           Increase verbosity level\n"
-	     "\t-g GAIN      Gain in e/ADU (1.0)\n"
-	     "\t-r READNOISE Read noise in e (0.0)\n"
-	     "\t-s SAT       Saturation point (65535)\n"
-	     "\t-b BAD       Bad level (0)\n"
-	     "\t-p FILE MAP  Specify file containing star coordinates, with map\n"
-	     "\t-a APER      Aperture radius for photometry (3.0)\n"
-	     "\t-k SIGMAREJ  k-sigma rejection threshold (3.0)\n"
-	     "\t-n NREJECT   Number of rejection iterations (1)\n"
-	     "\tOUT          Output image\n"
-	     "\tIN1 IN2...   Input images (identical size)\n",
-	     programName
-	);
+             "Usage: %s [-h] [-v] [-g GAIN] [-r READNOISE] [-s SAT] [-b BAD] [-p FILE MAP] [-a APER] [-k SIGMAREJ] [-n NREJECT] OUT IN1 IN2...\n"
+             "where\n"
+             "\t-h           Help (this info)\n"
+             "\t-v           Increase verbosity level\n"
+             "\t-g GAIN      Gain in e/ADU (1.0)\n"
+             "\t-r READNOISE Read noise in e (0.0)\n"
+             "\t-s SAT       Saturation point (65535)\n"
+             "\t-b BAD       Bad level (0)\n"
+             "\t-p FILE MAP  Specify file containing star coordinates, with map\n"
+             "\t-a APER      Aperture radius for photometry (3.0)\n"
+             "\t-k SIGMAREJ  k-sigma rejection threshold (3.0)\n"
+             "\t-n NREJECT   Number of rejection iterations (1)\n"
+             "\tOUT          Output image\n"
+             "\tIN1 IN2...   Input images (identical size)\n",
+             programName
+        );
 }
 
@@ -35,16 +35,16 @@
 
     // Parameters with default values
-    config->verbose = 0;		// Verbosity level
-    config->gain = 1.0;			// Gain (e/ADU)
-    config->readnoise = 0.0;		// Read noise (e)
-    config->reject = 4.0;		// Rejection threshold (sigma)
-    config->nReject = 1;		// Number of rejection iterations
-    config->saturated = 65535.0;	// Saturation level
-    config->bad = 0.0;			// Bad level
-    config->outName = NULL;		// Output name
-    config->inNames = NULL;		// Input names;
-    config->starFile = NULL;		// Filename of file containing stars
-    config->starMap = NULL;		// Map for stars
-    config->aper = 3.0;			// Aperture for photometry
+    config->verbose = 0;                // Verbosity level
+    config->gain = 1.0;                 // Gain (e/ADU)
+    config->readnoise = 0.0;            // Read noise (e)
+    config->reject = 4.0;               // Rejection threshold (sigma)
+    config->nReject = 1;                // Number of rejection iterations
+    config->saturated = 65535.0;        // Saturation level
+    config->bad = 0.0;                  // Bad level
+    config->outName = NULL;             // Output name
+    config->inNames = NULL;             // Input names;
+    config->starFile = NULL;            // Filename of file containing stars
+    config->starMap = NULL;             // Map for stars
+    config->aper = 3.0;                 // Aperture for photometry
 
     return config;
@@ -62,5 +62,5 @@
     combineConfig *config = combineConfigAlloc(); // Configuration
 
-    const char *programName = argv[0];	// Program name
+    const char *programName = argv[0];  // Program name
 
     /* Variables for getopt */
@@ -72,74 +72,74 @@
     while ((opt = getopt(argc, argv, "hvg:r:s:b:p:a:k:n:")) != -1) {
         switch (opt) {
-	  case 'h':
-	    help(programName);
-	    exit(EXIT_SUCCESS);
-	  case 'v':
+          case 'h':
+            help(programName);
+            exit(EXIT_SUCCESS);
+          case 'v':
             config->verbose++;
             break;
-	  case 'r':
+          case 'r':
             if (sscanf(optarg, "%f", &config->readnoise) != 1) {
-		printf("Unable to read readnoise.\n");
-                help(programName);
-                exit(EXIT_FAILURE);
-	    }
-	    break;
-	  case 'g':
-	    if (sscanf(optarg, "%f", &config->gain) != 1) {
- 		printf("Unable to read gain.\n");
-		help(programName);
-		exit(EXIT_FAILURE);
-	    }
-	    break;
-	  case 's':
-            if (sscanf(optarg, "%f", &config->saturated) != 1) {
- 		printf("Unable to read saturation limit.\n");
-                help(programName);
-                exit(EXIT_FAILURE);
-	    }
-	    break;
-	  case 'b':
-            if (sscanf(optarg, "%f", &config->bad) != 1) {
- 		printf("Unable to read bad limit.\n");
-		help(programName);
-                exit(EXIT_FAILURE);
-	    }
-	    break;
-	  case 'p':
-	    if (argc < optind+1) {
- 		printf("Unable to read photometric files.\n");
-		help(programName);
-		exit(EXIT_FAILURE);
-	    }
-	    config->starFile = argv[optind-1];
-	    config->starMap = argv[optind++];
-	    // Note: incrementing optind, so I can read more than one parameter.
-	    break;
-	  case 'a':
-            if (sscanf(optarg, "%f", &config->aper) != 1) {
- 		printf("Unable to read aperture.\n");
+                printf("Unable to read readnoise.\n");
                 help(programName);
                 exit(EXIT_FAILURE);
             }
-	    break;
-	  case 'k':
-            if (sscanf(optarg, "%f", &config->reject) != 1) {
- 		printf("Unable to read rejection limit.\n");
+            break;
+          case 'g':
+            if (sscanf(optarg, "%f", &config->gain) != 1) {
+                printf("Unable to read gain.\n");
                 help(programName);
                 exit(EXIT_FAILURE);
             }
-	    break;
-	  case 'n':
-            if (sscanf(optarg, "%d", &config->nReject) != 1) {
-  		printf("Unable to read number of rejection iterations.\n");
-		help(programName);
+            break;
+          case 's':
+            if (sscanf(optarg, "%f", &config->saturated) != 1) {
+                printf("Unable to read saturation limit.\n");
+                help(programName);
                 exit(EXIT_FAILURE);
             }
-	    break;
-	  default:
-	    printf("Bad option: %c\n", opt);
-	    help(programName);
-	    exit(EXIT_FAILURE);
-	}
+            break;
+          case 'b':
+            if (sscanf(optarg, "%f", &config->bad) != 1) {
+                printf("Unable to read bad limit.\n");
+                help(programName);
+                exit(EXIT_FAILURE);
+            }
+            break;
+          case 'p':
+            if (argc < optind+1) {
+                printf("Unable to read photometric files.\n");
+                help(programName);
+                exit(EXIT_FAILURE);
+            }
+            config->starFile = argv[optind-1];
+            config->starMap = argv[optind++];
+            // Note: incrementing optind, so I can read more than one parameter.
+            break;
+          case 'a':
+            if (sscanf(optarg, "%f", &config->aper) != 1) {
+                printf("Unable to read aperture.\n");
+                help(programName);
+                exit(EXIT_FAILURE);
+            }
+            break;
+          case 'k':
+            if (sscanf(optarg, "%f", &config->reject) != 1) {
+                printf("Unable to read rejection limit.\n");
+                help(programName);
+                exit(EXIT_FAILURE);
+            }
+            break;
+          case 'n':
+            if (sscanf(optarg, "%d", &config->nReject) != 1) {
+                printf("Unable to read number of rejection iterations.\n");
+                help(programName);
+                exit(EXIT_FAILURE);
+            }
+            break;
+          default:
+            printf("Bad option: %c\n", opt);
+            help(programName);
+            exit(EXIT_FAILURE);
+        }
     }
 
@@ -151,9 +151,10 @@
         exit(EXIT_FAILURE);
     }
-    config->outName = argv[0];		// Output filename
+    config->outName = argv[0];          // Output filename
     config->inNames = psArrayAlloc(argc-1); // Input filenames
+    config->inNames->n = argc-1;
     for (int i = 1; i < argc; i++) {
-	config->inNames->data[i-1] = psAlloc(strlen(argv[i]));
-	strncpy(config->inNames->data[i-1], argv[i], strlen(argv[i]));
+        config->inNames->data[i-1] = psAlloc(strlen(argv[i]));
+        strncpy(config->inNames->data[i-1], argv[i], strlen(argv[i]));
     }
 
Index: /trunk/stac/src/shift.c
===================================================================
--- /trunk/stac/src/shift.c	(revision 6886)
+++ /trunk/stac/src/shift.c	(revision 6887)
@@ -126,4 +126,5 @@
     psArray *maps = psArrayAlloc(1);    // An array of one
     psArray *masks = psArrayAlloc(1);   // An array of one
+    images->n = maps->n = masks->n = 1;
     images->data[0] = image;
     maps->data[0] = map;
@@ -132,5 +133,13 @@
     // Get size
     if (outnx == 0 || outny == 0) {
-        stacSize(&outnx, &outny, NULL, NULL, images, maps);
+        psVector *xSize = psVectorAlloc(1, PS_TYPE_S32); // A vector of one
+        psVector *ySize = psVectorAlloc(1, PS_TYPE_S32); // A vector of one
+        xSize->n = 1;
+        ySize->n = 1;
+        xSize->data.S32[0] = image->numCols;
+        ySize->data.S32[0] = image->numRows;
+        stacSize(&outnx, &outny, NULL, NULL, xSize, ySize, maps);
+        psFree(xSize);
+        psFree(ySize);
     }
 
Index: /trunk/stac/src/shiftSize.c
===================================================================
--- /trunk/stac/src/shiftSize.c	(revision 6886)
+++ /trunk/stac/src/shiftSize.c	(revision 6887)
@@ -17,11 +17,11 @@
 {
     fprintf (stderr, "shiftSize: Calculate size for output combined image\n"
-	     "Usage: %s [-h] [-v] IN1 IN2...\n"
-	     "where\n"
-	     "\t-h           Help (this info)\n"
-	     "\t-v           Increase verbosity level\n"
-	     "\tIN1, IN2...  Input images, which have associated .map files.\n",
-	     name
-	);
+             "Usage: %s [-h] [-v] IN1 IN2...\n"
+             "where\n"
+             "\t-h           Help (this info)\n"
+             "\t-v           Increase verbosity level\n"
+             "\tIN1, IN2...  Input images, which have associated .map files.\n",
+             name
+        );
 }
 
@@ -41,8 +41,8 @@
 
 
-    int verbose = 0;			// Verbosity level
+    int verbose = 0;                    // Verbosity level
 
     // Parse command line
-    const char *programName = argv[0];	// Program name
+    const char *programName = argv[0];  // Program name
     /* Variables for getopt */
     int opt;   /* Option, from getopt */
@@ -53,13 +53,13 @@
     while ((opt = getopt(argc, argv, "hv")) != -1) {
         switch (opt) {
-	  case 'h':
-	    help(programName);
-	    exit(EXIT_SUCCESS);
-	  case 'v':
+          case 'h':
+            help(programName);
+            exit(EXIT_SUCCESS);
+          case 'v':
             verbose++;
             break;
-	  default:
-	    help(programName);
-	}
+          default:
+            help(programName);
+        }
     }
 
@@ -73,11 +73,34 @@
 
     psArray *inputs = psArrayAlloc(argc); // Input filenames
+    inputs->n = argc;
     for (int i = 0; i < argc; i++) {
-	inputs->data[i] = psAlloc(strlen(argv[i]));
-	strncpy(inputs->data[i], argv[i], strlen(argv[i]));
+        inputs->data[i] = psAlloc(strlen(argv[i]));
+        strncpy(inputs->data[i], argv[i], strlen(argv[i]));
+        psTrace("stac.size", 8, "Input file: %s\n", inputs->data[i]);
     }
 
     // Read input files
-    psArray *images = stacReadImages(NULL, inputs);
+    psVector *xSize = psVectorAlloc(inputs->n, PS_TYPE_S32);
+    psVector *ySize = psVectorAlloc(inputs->n, PS_TYPE_S32);
+    xSize->n = inputs->n;
+    ySize->n = inputs->n;
+    for (int i = 0; i < inputs->n; i++) {
+        const char *filename = inputs->data[i]; // Name of FITS file
+        psFits *file = psFitsOpen(filename, "r"); // FITS file handle
+        psMetadata *header = psFitsReadHeader(NULL, file); // FITS header
+        psFitsClose(file);
+        bool mdok = true;               // Status of MD lookup
+        xSize->data.S32[i] = psMetadataLookupS32(&mdok, header, "NAXIS1");
+        if (!mdok) {
+            psError(PS_ERR_IO, true, "Unable to find NAXIS1 in %s\n", filename);
+            exit(EXIT_FAILURE);
+        }
+        ySize->data.S32[i] = psMetadataLookupS32(&mdok, header, "NAXIS2");
+        if (!mdok) {
+            psError(PS_ERR_IO, true, "Unable to find NAXIS2 in %s\n", filename);
+            exit(EXIT_FAILURE);
+        }
+        psFree(header);
+    }
 
     // Read maps
@@ -85,7 +108,9 @@
 
     // Get size
-    int outnx, outny;			// Size of combined image
-    float xMapDiff, yMapDiff;		// Difference to apply to maps
-    stacSize(&outnx, &outny, &xMapDiff, &yMapDiff, images, maps);
+    int outnx, outny;                   // Size of combined image
+    float xMapDiff, yMapDiff;           // Difference to apply to maps
+    stacSize(&outnx, &outny, &xMapDiff, &yMapDiff, xSize, ySize, maps);
+    psFree(xSize);
+    psFree(ySize);
 
     printf("%d %d\n", outnx, outny);
@@ -95,5 +120,4 @@
 
     psFree(inputs);
-    psFree(images);
     psFree(maps);
 
Index: /trunk/stac/src/stac.c
===================================================================
--- /trunk/stac/src/stac.c	(revision 6886)
+++ /trunk/stac/src/stac.c	(revision 6887)
@@ -45,4 +45,5 @@
     // Generate masks
     psArray *masks = psArrayAlloc(inputs->n);
+    masks->n = inputs->n;
     for (int i = 0; i < inputs->n; i++) {
         psImage *image = inputs->data[i]; // Image for which to get mask
@@ -67,5 +68,16 @@
     // Get size, if not input
     if (config->outnx == 0 || config->outny == 0) {
-        stacSize(&config->outnx, &config->outny, &config->xMapDiff, &config->yMapDiff, inputs, maps);
+        psVector *xSize = psVectorAlloc(inputs->n, PS_TYPE_S32); // Sizes of images in x
+        psVector *ySize = psVectorAlloc(inputs->n, PS_TYPE_S32); // Sizes of images in y
+        xSize->n = inputs->n;
+        ySize->n = inputs->n;
+        for (int i = 0; i < inputs->n; i++) {
+            psImage *image = inputs->data[i]; // The i-th image
+            xSize->data.S32[i] = image->numCols;
+            ySize->data.S32[i] = image->numRows;
+        }
+        stacSize(&config->outnx, &config->outny, &config->xMapDiff, &config->yMapDiff, xSize, ySize, maps);
+        psFree(xSize);
+        psFree(ySize);
     }
 
@@ -96,7 +108,7 @@
     for (int i = 0; i < inputs->n; i++) {
         char errName[MAXCHAR];          // Filename of error image
-        sprintf(errName,"%s.err",config->inputs->data[i]);
+        sprintf(errName,"%s.err", (char*)config->inputs->data[i]);
         psFits *errorFile = psFitsOpen(errName, "w");
-        if (!psFitsWriteImage(errorFile, NULL, errors->data[i], 0)) {
+        if (!psFitsWriteImage(errorFile, NULL, errors->data[i], 0, NULL)) {
             psErrorStackPrint(stderr, "Unable to write image: %s\n", errName);
         }
@@ -118,15 +130,15 @@
         char shiftName[MAXCHAR];        // Filename of shift image
         char errName[MAXCHAR];          // Filename of error image
-        sprintf(shiftName,"%s.shift1",config->inputs->data[i]);
-        sprintf(errName,"%s.shifterr1",config->inputs->data[i]);
+        sprintf(shiftName,"%s.shift1", (char*)config->inputs->data[i]);
+        sprintf(errName,"%s.shifterr1",(char*)config->inputs->data[i]);
         psFits *shiftFile = psFitsOpen(shiftName, "w");
         psFits *errFile = psFitsOpen(errName, "w");
         psImage *trans = transformed->data[i]; // Transformed image
         psImage *transErr = transformedErrors->data[i]; // Transformed error image
-        if (!psFitsWriteImage(shiftFile, NULL, trans, 0)) {
+        if (!psFitsWriteImage(shiftFile, NULL, trans, 0, NULL)) {
             psErrorStackPrint(stderr, "Unable to write image: %s\n", shiftName);
         }
         psTrace("stac", 1, "Shifted image written to %s\n", shiftName);
-        if (!psFitsWriteImage(errFile, NULL, transErr, 0)) {
+        if (!psFitsWriteImage(errFile, NULL, transErr, 0, NULL)) {
             psErrorStackPrint(stderr, "Unable to write image: %s\n", errName);
         }
@@ -146,4 +158,6 @@
     psVector *saturated = psVectorAlloc(transformed->n, PS_TYPE_F32); // Saturation limits
     psVector *bad = psVectorAlloc(transformed->n, PS_TYPE_F32); // Bad limits
+    saturated->n = transformed->n;
+    bad->n = transformed->n;
     for (int i = 0; i < transformed->n; i++) {
         saturated->data.F32[i] = (config->saturated - offsets->data.F32[i]) / scales->data.F32[i];
@@ -182,8 +196,8 @@
     for (int i = 0; i < rejected->n; i++) {
         char rejName[MAXCHAR];  // Filename of rejection image
-        sprintf(rejName, "%s.shiftrej", config->inputs->data[i]);
+        sprintf(rejName, "%s.shiftrej", (char*)config->inputs->data[i]);
 
         psFits *rejFile = psFitsOpen(rejName, "w");
-        if (!psFitsWriteImage(rejFile, NULL, rejected->data[i], 0)) {
+        if (!psFitsWriteImage(rejFile, NULL, rejected->data[i], 0, NULL)) {
             psErrorStackPrint(stderr, "Unable to write image: %s\n", rejName);
         }
@@ -196,5 +210,5 @@
     sprintf(preName, "%s.pre", config->output);
     psFits *preFile = psFitsOpen(preName, "w");
-    if (!psFitsWriteImage(preFile, NULL, combined, 0)) {
+    if (!psFitsWriteImage(preFile, NULL, combined, 0, NULL)) {
         psErrorStackPrint(stderr, "Unable to write image: %s\n", preName);
     }
@@ -205,4 +219,5 @@
     // Get regions of interest in the source frame
     psArray *regions = psArrayAlloc(inputs->n); // Array of images denoting regions of interest
+    regions->n = inputs->n;
     for (int i = 0; i < inputs->n; i++) {
         regions->data[i] = stacAreaOfInterest(rejected->data[i], inverseMaps->data[i],
@@ -211,7 +226,7 @@
 #ifdef TESTING
         char regionName[MAXCHAR];       // Filename of region image
-        sprintf(regionName,"%s.region",config->inputs->data[i]);
+        sprintf(regionName,"%s.region",(char*)config->inputs->data[i]);
         psFits *regionFile = psFitsOpen(regionName, "w");
-        if (!psFitsWriteImage(regionFile, NULL, regions->data[i], 0)) {
+        if (!psFitsWriteImage(regionFile, NULL, regions->data[i], 0, NULL)) {
             psErrorStackPrint(stderr, "Unable to write image: %s\n", regionName);
         }
Index: /trunk/stac/src/stac.h
===================================================================
--- /trunk/stac/src/stac.h	(revision 6886)
+++ /trunk/stac/src/stac.h	(revision 6887)
@@ -19,5 +19,5 @@
 // Read the input files and return an array of images
 psArray *stacReadImages(psArray **headers, // The image headers, to be returned
-			psArray *filenames // The file names of the images
+                        psArray *filenames // The file names of the images
     );
 
@@ -38,6 +38,6 @@
 // Calculate the error images
 psArray *stacErrorImages(psArray *inputs, // Array of input images
-			 float gain,	// Gain, in e/ADU
-			 float rn	// Read noise, in e
+                         float gain,    // Gain, in e/ADU
+                         float rn       // Read noise, in e
     );
 
@@ -47,14 +47,14 @@
 
 // Transform input images, return success
-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
     );
 
@@ -68,6 +68,6 @@
 // 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
+                       const char *file, ///< the file in which the problem originated
+                       psS32 lineno     ///< the line number in which the problem originated
     );
 
@@ -81,25 +81,25 @@
 
 // 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
+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
+                        psVector *errors, // Errors in the values
+                        psVector *masks // Masks for the values, 0 = don't use, 1 = use
     );
 
 // Combine the transformed images
-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)
     );
 
@@ -110,5 +110,5 @@
 // Invert an array of maps
 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
     );
 
@@ -118,17 +118,17 @@
 
 // Return the relative gradient for a given pixel
-float stacGradient(psImage *image,	// Input for which to measure the gradient
-		   int x, int y		// Coordinates at which to measure the gradient
+float stacGradient(psImage *image,      // Input for which to measure the gradient
+                   int x, int y         // Coordinates at which to measure the gradient
     );
 
 // Transform the rejection masks back to the source frame
-psArray *stacRejection(psArray *inputs,	// Input images
-		       psArray *rejected, // Rejected images
-		       psArray *regions, // Regions of interest
-		       psArray *maps,	// Maps from input to transformed image
-		       psArray *inverseMaps, // Maps from transformed to input image
-		       float frac,	// Fraction of pixel rejected before looking more carefully
-		       float grad,	// Gradient limit for rejection
-		       psArray *names	// Names of original images (only for writing out when TESTING)
+psArray *stacRejection(psArray *inputs, // Input images
+                       psArray *rejected, // Rejected images
+                       psArray *regions, // Regions of interest
+                       psArray *maps,   // Maps from input to transformed image
+                       psArray *inverseMaps, // Maps from transformed to input image
+                       float frac,      // Fraction of pixel rejected before looking more carefully
+                       float grad,      // Gradient limit for rejection
+                       psArray *names   // Names of original images (only for writing out when TESTING)
     );
 
@@ -139,6 +139,6 @@
 // Returns the change in x' and y' for a change in x and y of 1
 bool stacPlaneTransformDeriv(psPlane *out, // Output derivative, assumed already allocated
-			     psPlaneTransform *transform, // The transform for which to obtain the derivative
-			     psPlane *in // The position at which to obtain the derivative
+                             psPlaneTransform *transform, // The transform for which to obtain the derivative
+                             psPlane *in // The position at which to obtain the derivative
     );
 
@@ -146,6 +146,6 @@
 // Basically, uses derivatives of the transformation to work out what pixels in the source might be of interest
 psImage *stacAreaOfInterest(psImage *mask, // Mask that is to be used to determine the area of interest
-			    psPlaneTransform *map, // Map from the mask to the area of interest
-			    int nxOut, int nyOut // Size of the area of interest
+                            psPlaneTransform *map, // Map from the mask to the area of interest
+                            int nxOut, int nyOut // Size of the area of interest
     );
 
@@ -155,8 +155,9 @@
 
 // Calculate the size of the output image
-bool stacSize(int *outnx, int *outny,	// Size of output image (returned)
-	      float *xMapDiff, float *yMapDiff, // Difference applied to maps
-	      const psArray *images,	// Input images
-	      psArray *maps		// Transformation maps
+bool stacSize(int *outnx, int *outny,   // Size of output image (returned)
+              float *xMapDiff, float *yMapDiff, // Difference applied to maps
+              const psVector *xSizes,   // Sizes of images in x
+              const psVector *ySizes,   // Sizes of images in y
+              psArray *maps             // Transformation maps
     );
 
@@ -167,24 +168,24 @@
 // Calculate the background in an image
 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
     );
 
 
 // Calculate the relative scales and offsets between the images
-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)
     );
 
 // Rescale images
-bool stacRescale(psArray *images,	// Images to rescale
-		 psArray *errImages,	// Variance images to rescale
-		 const psImage *mask,	// Mask indicating which 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 indicating which pixels to scale
+                 const psVector *scales,// Scales for images
+                 const psVector *offsets // Offsets for images
     );
 
@@ -196,11 +197,11 @@
 
 // Write map out
-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
     );
 
 // Write multiple maps
 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
     );
 
Index: /trunk/stac/src/stacCombine.c
===================================================================
--- /trunk/stac/src/stacCombine.c	(revision 6886)
+++ /trunk/stac/src/stacCombine.c	(revision 6887)
@@ -103,4 +103,7 @@
     psVector *deltas = psVectorAlloc(nImages, PS_TYPE_F32); // Will hold the errors in the statistics step
     psVector *mask = psVectorAlloc(nImages, PS_TYPE_U8); // Will hold the mask in the statistics step
+    pixels->n = nImages;
+    deltas->n = nImages;
+    mask->n = nImages;
 
     // Set up rejection masks
@@ -109,4 +112,5 @@
             // Allocate the rejection masks, if required
             *rejected = psArrayAlloc(nImages);
+            (*rejected)->n = nImages;
         } else {
             assert((*rejected)->n != nImages);
@@ -202,10 +206,10 @@
     char chiName[MAXCHAR];              // Filename of chi^2 image
     sprintf(chiName,"chi2_%d.fits",iteration);
-    psFits *chiFile = psFitsAlloc(chiName);
-    if (!psFitsWriteImage(chiFile, NULL, chi2 , 0)) {
+    psFits *chiFile = psFitsOpen(chiName, "w");
+    if (!psFitsWriteImage(chiFile, NULL, chi2 , 0, NULL)) {
         psErrorStackPrint(stderr, "Unable to write image: %s\n", chiName);
     }
     psTrace("stac.combine", 1, "Chi^2 image written to %s\n", chiName);
-    psFree(chiFile);
+    psFitsClose(chiFile);
     psFree(chi2);
 #endif
Index: /trunk/stac/src/stacConfig.c
===================================================================
--- /trunk/stac/src/stacConfig.c	(revision 6886)
+++ /trunk/stac/src/stacConfig.c	(revision 6887)
@@ -24,6 +24,6 @@
     config->outnx = 0;
     config->outny = 0;
-    config->saturated = 65535.0;	// Saturation level
-    config->bad = 0.0;			// Bad level
+    config->saturated = 65535.0;        // Saturation level
+    config->bad = 0.0;                  // Bad level
     config->reject = 2.5;
     config->frac = 0.45;
@@ -39,5 +39,5 @@
     // Free the vectors, if necessary
     if (config->inputs) {
-	psFree(config->inputs);
+        psFree(config->inputs);
     }
     // Free everything
@@ -46,10 +46,10 @@
 
 
-stacConfig *stacParseConfig(int argc,	// Number of command-line arguments
-			    char **argv // Command-line arguments
+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
+    const char *programName = argv[0];  // Program name
 
     /* Variables for getopt */
@@ -61,23 +61,23 @@
     while ((opt = getopt(argc, argv, "hvSg:r:o:s:b:k:n:f:G:p:a:")) != -1) {
         switch (opt) {
-	  case 'h':
-	    help(programName);
-	    exit(EXIT_SUCCESS);
-	  case 'v':
+          case 'h':
+            help(programName);
+            exit(EXIT_SUCCESS);
+          case 'v':
             config->verbose++;
             break;
-	  case 'g':
-	    if (sscanf(optarg, "%f", &config->gain) != 1) {
-		help(programName);
-		exit(EXIT_FAILURE);
-	    }
-	    break;
-	  case 'r':
-	    if (sscanf(optarg, "%f", &config->readnoise) != 1) {
-		help(programName);
-		exit(EXIT_FAILURE);
-	    }
-	    break;
-	  case 'o':
+          case 'g':
+            if (sscanf(optarg, "%f", &config->gain) != 1) {
+                help(programName);
+                exit(EXIT_FAILURE);
+            }
+            break;
+          case 'r':
+            if (sscanf(optarg, "%f", &config->readnoise) != 1) {
+                help(programName);
+                exit(EXIT_FAILURE);
+            }
+            break;
+          case 'o':
             if ((sscanf(argv[optind-1], "%d", &config->outnx) != 1) ||
                 (sscanf(argv[optind++], "%d", &config->outny) != 1)) {
@@ -86,62 +86,62 @@
                 exit(EXIT_FAILURE);
             }
-	    break;
-	  case 's':
-	    if (sscanf(optarg, "%f", &config->saturated) != 1) {
-		help(programName);
-		exit(EXIT_FAILURE);
-	    }
-	    break;
-	  case 'b':
-	    if (sscanf(optarg, "%f", &config->bad) != 1) {
-		help(programName);
-		exit(EXIT_FAILURE);
-	    }
-	    break;
-	  case 'p':
-	    if (argc < optind+1) {
-		help(programName);
-		exit(EXIT_FAILURE);
-	    }
-	    config->starFile = argv[optind-1];
-	    config->starMapFile = argv[optind++];
-	    // Note: incrementing optind, so I can read more than one parameter.
-	    break;
-	  case 'a':
-	    if (sscanf(optarg, "%f", &config->aper) != 1) {
-		help(programName);
-		exit(EXIT_FAILURE);
-	    }
-	    break;
-	  case 'k':
-	    if (sscanf(optarg, "%f", &config->reject) != 1) {
-		help(programName);
-		exit(EXIT_FAILURE);
-	    }
-	    break;
-	  case 'n':
-	    if (sscanf(optarg, "%d", &config->nReject) != 1) {
-		help(programName);
-		exit(EXIT_FAILURE);
-	    }
-	    break;
-	  case 'f':
-	    if (sscanf(optarg, "%f", &config->frac) != 1) {
-		help(programName);
-		exit(EXIT_FAILURE);
-	    }
-	    break;
-	  case 'G':
-	    if (sscanf(optarg, "%f", &config->grad) != 1) {
-		help(programName);
-		exit(EXIT_FAILURE);
-	    }
-	    break;
-	  case 'S':
-	    config->saveShifts = true;
-	    break;
-	  default:
-	    help(programName);
-	}
+            break;
+          case 's':
+            if (sscanf(optarg, "%f", &config->saturated) != 1) {
+                help(programName);
+                exit(EXIT_FAILURE);
+            }
+            break;
+          case 'b':
+            if (sscanf(optarg, "%f", &config->bad) != 1) {
+                help(programName);
+                exit(EXIT_FAILURE);
+            }
+            break;
+          case 'p':
+            if (argc < optind+1) {
+                help(programName);
+                exit(EXIT_FAILURE);
+            }
+            config->starFile = argv[optind-1];
+            config->starMapFile = argv[optind++];
+            // Note: incrementing optind, so I can read more than one parameter.
+            break;
+          case 'a':
+            if (sscanf(optarg, "%f", &config->aper) != 1) {
+                help(programName);
+                exit(EXIT_FAILURE);
+            }
+            break;
+          case 'k':
+            if (sscanf(optarg, "%f", &config->reject) != 1) {
+                help(programName);
+                exit(EXIT_FAILURE);
+            }
+            break;
+          case 'n':
+            if (sscanf(optarg, "%d", &config->nReject) != 1) {
+                help(programName);
+                exit(EXIT_FAILURE);
+            }
+            break;
+          case 'f':
+            if (sscanf(optarg, "%f", &config->frac) != 1) {
+                help(programName);
+                exit(EXIT_FAILURE);
+            }
+            break;
+          case 'G':
+            if (sscanf(optarg, "%f", &config->grad) != 1) {
+                help(programName);
+                exit(EXIT_FAILURE);
+            }
+            break;
+          case 'S':
+            config->saveShifts = true;
+            break;
+          default:
+            help(programName);
+        }
     }
 
@@ -154,10 +154,11 @@
     }
 
-    config->output = argv[0];		// Output file
+    config->output = argv[0];           // Output file
     // Get the input files
-    config->inputs = psArrayAlloc(argc-1);
+    config->inputs = psArrayAlloc(argc - 1);
+    config->inputs->n = 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]));
+        config->inputs->data[i-1] = psAlloc(strlen(argv[i]));
+        strncpy(config->inputs->data[i-1], argv[i], strlen(argv[i]));
     }
 
@@ -167,5 +168,5 @@
     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, "\t%s\n", config->inputs->data[i]);
     }
     psTrace("stac.config", 9, "Output file is %s\n",config->output);
Index: /trunk/stac/src/stacErrorImages.c
===================================================================
--- /trunk/stac/src/stacErrorImages.c	(revision 6886)
+++ /trunk/stac/src/stacErrorImages.c	(revision 6887)
@@ -4,11 +4,12 @@
 
 psArray *stacErrorImages(psArray *inputs, // Array of input images
-			 float gain,	// Gain, in e/ADU
-			 float rn	// Read noise, in e
+                         float gain,    // Gain, in e/ADU
+                         float rn       // Read noise, in e
     )
 {
-    float invGain = 1.0/gain;		// Inverse square root of gain
-    rn /= gain;				// Read noise in ADU
+    float invGain = 1.0/gain;           // Inverse square root of gain
+    rn /= gain;                         // Read noise in ADU
     psArray *errors = psArrayAlloc(inputs->n);
+    errors->n = inputs->n;
 
     psTrace("stac.errors", 1, "Calculating error images....\n");
@@ -16,21 +17,21 @@
     // Iterate over the input images
     for (int i = 0; i < inputs->n; i++) {
-	psTrace("stac.errors",5,"Working on image #%d\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
+        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;
-	    }
-	}
+        // 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;
+        // Put image onto the array
+        errors->data[i] = error;
     }
 
Index: /trunk/stac/src/stacInvertMaps.c
===================================================================
--- /trunk/stac/src/stacInvertMaps.c	(revision 6886)
+++ /trunk/stac/src/stacInvertMaps.c	(revision 6887)
@@ -14,4 +14,5 @@
     int nMaps = maps->n;                // Number of maps
     psArray *inverted = psArrayAlloc(nMaps); // Array of inverted maps for output
+    inverted->n = nMaps;
 
     psTrace("stac.invertMaps", 1, "Inverting maps....\n");
@@ -45,4 +46,8 @@
         psVector *xOut = psVectorAlloc(NUM_GRID * NUM_GRID, PS_TYPE_F32);
         psVector *yOut = psVectorAlloc(NUM_GRID * NUM_GRID, PS_TYPE_F32);
+        xIn->n = NUM_GRID * NUM_GRID;
+        yIn->n = NUM_GRID * NUM_GRID;
+        xOut->n = NUM_GRID * NUM_GRID;
+        yOut->n = NUM_GRID * NUM_GRID;
 
         // Create grid of points
@@ -66,4 +71,6 @@
         psVector *xVector = psVectorAlloc(nCoeff, PS_TYPE_F64); // Vector for solution in x
         psVector *yVector = psVectorAlloc(nCoeff, PS_TYPE_F64); // Vector for solution in y
+        xVector->n = nCoeff;
+        yVector->n = nCoeff;
         for (int i = 0; i < nCoeff; i++) {
             for (int j = 0; j < nCoeff; j++) {
@@ -104,4 +111,5 @@
         // Solution via LU Decomposition
         psVector *permutation = psVectorAlloc(nCoeff, PS_TYPE_F64); // Permutation vector for LU Decomposition
+        permutation->n = nCoeff;
         psImage *luMatrix = psMatrixLUD(NULL, &permutation, matrix); // LU decomposed matrix
         psVector *xSolution = psMatrixLUSolve(NULL, luMatrix, xVector, permutation); // Solution in x
Index: /trunk/stac/src/stacRead.c
===================================================================
--- /trunk/stac/src/stacRead.c	(revision 6886)
+++ /trunk/stac/src/stacRead.c	(revision 6887)
@@ -13,7 +13,9 @@
     int nFiles = filenames->n;          // The number of input files
     psArray *images = psArrayAlloc(nFiles); // The input files, to be returned
+    images->n = nFiles;
     assert(!headers || ! *headers || (*headers)->n == nFiles);
     if (headers && ! *headers) {
         *headers = psArrayAlloc(nFiles);
+        (*headers)->n = nFiles;
     }
 
@@ -71,16 +73,14 @@
     psArray *coords = psArrayAlloc(BUFFER); // The array of coordinates to be returned
     float x, y;                         // Coordinates to read
-    int num = 0;                        // Number of coordinates read
     while (fscanf(file, "%f %f\n", &x, &y) == 2) {
         psPlane *coord = psPlaneAlloc();// A coordinate
         coord->x = x;
         coord->y = y;
-        coords->data[num] = coord;
-        num++;
-        if (num % BUFFER) {
-            coords = psArrayRealloc(coords, num + BUFFER);
-        }
-    }
-    coords->n = num;
+        coords->data[coords->n] = coord;
+        coords->n++;
+        if (coords->n % BUFFER == 0) {
+            coords = psArrayRealloc(coords, coords->n + BUFFER);
+        }
+    }
 
     psTrace("stac.read.coords", 5, "%d coordinates read.\n", num);
@@ -184,4 +184,5 @@
     int nFiles = filenames->n;          // The number of input files
     psArray *maps = psArrayAlloc(nFiles); // The maps, to be returned
+    maps->n = nFiles;
     char mapfile[MAXCHAR];              // Filename of map
 
Index: /trunk/stac/src/stacRejection.c
===================================================================
--- /trunk/stac/src/stacRejection.c	(revision 6886)
+++ /trunk/stac/src/stacRejection.c	(revision 6887)
@@ -8,6 +8,6 @@
 #define MIN(x,y) ((x) < (y) ? (x) : (y))
 
-float stacGradient(psImage *image,	// Input for which to measure the gradient
-		   int x, int y		// Coordinates at which to measure the gradient
+float stacGradient(psImage *image,      // Input for which to measure the gradient
+                   int x, int y         // Coordinates at which to measure the gradient
     )
 {
@@ -22,11 +22,11 @@
     int yMax = MIN(y + 1, image->numRows - 1);
     for (int j = yMin; j <= yMax; j++) {
-	for (int i = xMin; i <= xMax; i++) {
-	    if ((i != x) && (j != y)) {
-		pixels->data.F32[num] = image->data.F32[j][i];
-		mask->data.U8[num] = 0;
-		num++;
-	    }
-	}
+        for (int i = xMin; i <= xMax; i++) {
+            if ((i != x) && (j != y)) {
+                pixels->data.F32[num] = image->data.F32[j][i];
+                mask->data.U8[num] = 0;
+                num++;
+            }
+        }
     }
     pixels->n = num;
@@ -43,15 +43,15 @@
 }
 
-psArray *stacRejection(psArray *inputs,	// Input images
-		       psArray *rejected, // Rejected images
-		       psArray *regions, // Regions of interest
-		       psArray *maps,	// Maps from input to transformed image
-		       psArray *inverseMaps, // Maps from transformed to input image
-		       float frac,	// Fraction of pixel rejected before looking more carefully
-		       float grad,	// Gradient limit for rejection
-		       psArray *names	// Names of original images (only for writing out when TESTING)
+psArray *stacRejection(psArray *inputs, // Input images
+                       psArray *rejected, // Rejected images
+                       psArray *regions, // Regions of interest
+                       psArray *maps,   // Maps from input to transformed image
+                       psArray *inverseMaps, // Maps from transformed to input image
+                       float frac,      // Fraction of pixel rejected before looking more carefully
+                       float grad,      // Gradient limit for rejection
+                       psArray *names   // Names of original images (only for writing out when TESTING)
     )
 {
-    int nImages = inputs->n;		// Number of input images
+    int nImages = inputs->n;            // Number of input images
 
     // Check inputs
@@ -63,7 +63,7 @@
 
     for (int i = 0; i < nImages; i++) {
-	psImage *input = inputs->data[i];
-	psImage *region = regions->data[i];
-	assert(input->numRows == region->numRows && input->numCols == region->numCols);
+        psImage *input = inputs->data[i];
+        psImage *region = regions->data[i];
+        assert(input->numRows == region->numRows && input->numCols == region->numCols);
     }
 
@@ -77,156 +77,159 @@
     psVector *grads = psVectorAlloc(nImages, PS_TYPE_F32); // Gradient for each image
     psVector *gradsMask = psVectorAlloc(nImages, PS_TYPE_U8); // Mask for gradient vector
+    grads->n = nImages;
+    gradsMask->n = nImages;
 
     // Transform rejection masks back to source
     psArray *inputRej = psArrayAlloc(nImages);
+    inputRej->n = 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_U8); // The pixel mask for the input
-#ifdef TESTING
-	psImage *rejmap = psImageAlloc(nxInput, nyInput, PS_TYPE_F32); // The rejections in the source frame
-	psImage *gradient = psImageAlloc(nxInput, nyInput, PS_TYPE_F32);	// The gradient image
-#endif
-	psImage *reject = rejected->data[i]; // Pull out the mask in the output frame
-	psPlaneTransform *map = maps->data[i]; // The map from input to output
-	psImage *region = regions->data[i]; // The region of interest for this image
-
-	psTrace("stac.rejection", 3, "Transforming rejection mask %d....\n", i);
-	int nBad = 0;			// Number of bad pixels
+        // 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_U8); // The pixel mask for the input
+#ifdef TESTING
+        psImage *rejmap = psImageAlloc(nxInput, nyInput, PS_TYPE_F32); // The rejections in the source frame
+        psImage *gradient = psImageAlloc(nxInput, nyInput, PS_TYPE_F32);        // The gradient image
+#endif
+        psImage *reject = rejected->data[i]; // Pull out the mask in the output frame
+        psPlaneTransform *map = maps->data[i]; // The map from input to output
+        psImage *region = regions->data[i]; // The region of interest for this image
+
+        psTrace("stac.rejection", 3, "Transforming rejection mask %d....\n", i);
+        int nBad = 0;                   // Number of bad pixels
 
 #ifdef CRFLUX
-	// Set up CR output
-	FILE *crs = NULL;		// File for outputting details of rejected pixels
-	char crfile[MAXCHAR];	// Filename
-	sprintf(crfile,"%s.cr",names->data[i]);
-	if ((crs = fopen(crfile, "w")) == NULL) {
-	    fprintf(stderr, "Unable to open file for detailed output, %s\n");
-	    return NULL;
-	}
-#endif
-
-	// 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.
-	psStats *median = psStatsAlloc(PS_STAT_SAMPLE_MEDIAN);
-	for (int y = 0; y < nyInput; y++) {
-	    for (int x = 0; x < nxInput; x++) {
-
-		// Only transform pixels of interest
-		if (region->data.U8[y][x]) {
-
-		    inCoords->x = (double)x + 0.5;
-		    inCoords->y = (double)y + 0.5;
-		    (void)psPlaneTransformApply(outCoords, map, inCoords);
-		    float maskVal = (float)psImagePixelInterpolate(reject, outCoords->x, outCoords->y,
-								   NULL, 0, 0.0, PS_INTERPOLATE_BILINEAR);
-#ifdef TESTING
-		    rejmap->data.F32[y][x] = maskVal;
-#endif
-
-		    if (maskVal > frac) {
-			// Calculate mean gradient on other images
-			float meanGrads = 0.0;
-			int numGrads = 0;
-			for (int j = 0; j < nImages; j++) {
-			    if (i != j) {
-				// Get coordinates for that image
-				(void)psPlaneTransformApply(inCoords, inverseMaps->data[j], outCoords);
-				int xPix = (int)(inCoords->x + 0.5);
-				int yPix = (int)(inCoords->y + 0.5);
-				if ((xPix >= 0) && (xPix <= ((psImage*)(inputs->data[j]))->numCols - 1) &&
-				    (yPix >= 0) && (yPix <= ((psImage*)(inputs->data[j]))->numRows - 1)) {
-				    // Calculate the gradient
-				    grads->data.F32[j] = stacGradient(inputs->data[j], xPix, yPix);
-				    gradsMask->data.U8[j] = 0;
-				    numGrads++;
-				} else {
-				    gradsMask->data.U8[j] = 1; // Mask this one
-				}
-			    } else {
-				gradsMask->data.U8[j] = 1; // Mask this one
-			    }
-			}
-			if (numGrads > 0) {
-			    (void)psVectorStats(median, grads, NULL, gradsMask, (psU8)1);
-			    meanGrads = median->sampleMedian;
-			} else {
-			    meanGrads = 0.0;
-			}
-
-#ifdef TESTING
-			//gradient->data.F32[y][x] = stacGradient(inputs->data[i], x, y) / meanGrads;
-			gradient->data.F32[y][x] = meanGrads;
-#endif
-
-			//if (stacGradient(inputs->data[i], x, y) < grad * meanGrads) {
-			if (meanGrads > grad) {
-			    mask->data.U8[y][x] = 1;
-			    nBad++;
+        // Set up CR output
+        FILE *crs = NULL;               // File for outputting details of rejected pixels
+        char crfile[MAXCHAR];   // Filename
+        sprintf(crfile,"%s.cr",names->data[i]);
+        if ((crs = fopen(crfile, "w")) == NULL) {
+            fprintf(stderr, "Unable to open file for detailed output, %s\n");
+            return NULL;
+        }
+#endif
+
+        // 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.
+        psStats *median = psStatsAlloc(PS_STAT_SAMPLE_MEDIAN);
+        for (int y = 0; y < nyInput; y++) {
+            for (int x = 0; x < nxInput; x++) {
+
+                // Only transform pixels of interest
+                if (region->data.U8[y][x]) {
+
+                    inCoords->x = (double)x + 0.5;
+                    inCoords->y = (double)y + 0.5;
+                    (void)psPlaneTransformApply(outCoords, map, inCoords);
+                    float maskVal = (float)psImagePixelInterpolate(reject, outCoords->x, outCoords->y,
+                                                                   NULL, 0, 0.0, PS_INTERPOLATE_BILINEAR);
+#ifdef TESTING
+                    rejmap->data.F32[y][x] = maskVal;
+#endif
+
+                    if (maskVal > frac) {
+                        // Calculate mean gradient on other images
+                        float meanGrads = 0.0;
+                        int numGrads = 0;
+                        for (int j = 0; j < nImages; j++) {
+                            if (i != j) {
+                                // Get coordinates for that image
+                                (void)psPlaneTransformApply(inCoords, inverseMaps->data[j], outCoords);
+                                int xPix = (int)(inCoords->x + 0.5);
+                                int yPix = (int)(inCoords->y + 0.5);
+                                if ((xPix >= 0) && (xPix <= ((psImage*)(inputs->data[j]))->numCols - 1) &&
+                                    (yPix >= 0) && (yPix <= ((psImage*)(inputs->data[j]))->numRows - 1)) {
+                                    // Calculate the gradient
+                                    grads->data.F32[j] = stacGradient(inputs->data[j], xPix, yPix);
+                                    gradsMask->data.U8[j] = 0;
+                                    numGrads++;
+                                } else {
+                                    gradsMask->data.U8[j] = 1; // Mask this one
+                                }
+                            } else {
+                                gradsMask->data.U8[j] = 1; // Mask this one
+                            }
+                        }
+                        if (numGrads > 0) {
+                            (void)psVectorStats(median, grads, NULL, gradsMask, (psU8)1);
+                            meanGrads = median->sampleMedian;
+                        } else {
+                            meanGrads = 0.0;
+                        }
+
+#ifdef TESTING
+                        //gradient->data.F32[y][x] = stacGradient(inputs->data[i], x, y) / meanGrads;
+                        gradient->data.F32[y][x] = meanGrads;
+#endif
+
+                        //if (stacGradient(inputs->data[i], x, y) < grad * meanGrads) {
+                        if (meanGrads > grad) {
+                            mask->data.U8[y][x] = 1;
+                            nBad++;
 
 #ifdef CRFLUX
-			    fprintf(crs, "%d %d --> %f %f %f\n", x, y,
-				    ((psImage*)(inputs->data[i]))->data.F32[y][x], maskVal,
-				    stacGradient(inputs->data[i], x, y));
-#endif
-
-			} else {
-			    mask->data.U8[y][x] = 0;
-			} // Gradient threshold
-		    } else {
-			mask->data.U8[y][x] = 0;
-		    } // Rejection threshold
-
-		} else {
-		    mask->data.U8[y][x] = 0;
-		} // Only touching pixels of interest
-
-	    }
-	} // Iterating over pixels
-	psFree(median);
+                            fprintf(crs, "%d %d --> %f %f %f\n", x, y,
+                                    ((psImage*)(inputs->data[i]))->data.F32[y][x], maskVal,
+                                    stacGradient(inputs->data[i], x, y));
+#endif
+
+                        } else {
+                            mask->data.U8[y][x] = 0;
+                        } // Gradient threshold
+                    } else {
+                        mask->data.U8[y][x] = 0;
+                    } // Rejection threshold
+
+                } else {
+                    mask->data.U8[y][x] = 0;
+                } // Only touching pixels of interest
+
+            }
+        } // Iterating over pixels
+        psFree(median);
 
 #ifdef CRFLUX
-	// Close file used for detailed output
-	fclose(crs);
-#endif
-
-	psTrace("stac.rejection", 1, "%d pixels masked in image %d\n", nBad, i);
-	// Clip the image, and convert to suitable mask format
-
-#ifdef TESTING
-	// Write error image out to check
-	char maskName[MAXCHAR];		// Filename of mask image
-	char rejmapName[MAXCHAR]; 	// Filename of rejection image
-	char gradName[MAXCHAR];		// Filename of gradient image
-	sprintf(maskName, "%s.mask", names->data[i]);
-	sprintf(rejmapName, "%s.rejmap", names->data[i]);
-	sprintf(gradName, "%s.grad", names->data[i]);
-
-	psFits *maskFile = psFitsAlloc(maskName);
-	psFits *rejmapFile = psFitsAlloc(rejmapName);
-	psFits *gradFile = psFitsAlloc(gradName);
-	if (!psFitsWriteImage(maskFile, NULL, mask, 0)) {
-	    psErrorStackPrint(stderr, "Unable to write image: %s\n", maskName);
-	}
-	psTrace("stac", 1, "Mask image written to %s\n", maskName);
-	if (!psFitsWriteImage(rejmapFile, NULL, rejmap, 0)) {
-	    psErrorStackPrint(stderr, "Unable to write image: %s\n", rejmapName);
-	}
-	psTrace("stac", 1, "Rejection map written to %s\n", rejmapName);
-	if (!psFitsWriteImage(gradFile, NULL, gradient, 0)) {
-	    psErrorStackPrint(stderr, "Unable to write image: %s\n", gradName);
-	}
-	psTrace("stac", 1, "Gradient image written to %s\n", gradName);
-	psFree(maskFile);
-	psFree(rejmapFile);
-	psFree(gradFile);
-	psFree(rejmap);
-	psFree(gradient);
-#endif
-
-	// Stuff into the array
-	inputRej->data[i] = mask;
+        // Close file used for detailed output
+        fclose(crs);
+#endif
+
+        psTrace("stac.rejection", 1, "%d pixels masked in image %d\n", nBad, i);
+        // Clip the image, and convert to suitable mask format
+
+#ifdef TESTING
+        // Write error image out to check
+        char maskName[MAXCHAR];         // Filename of mask image
+        char rejmapName[MAXCHAR];       // Filename of rejection image
+        char gradName[MAXCHAR];         // Filename of gradient image
+        sprintf(maskName, "%s.mask", (char*)names->data[i]);
+        sprintf(rejmapName, "%s.rejmap", (char*)names->data[i]);
+        sprintf(gradName, "%s.grad", (char*)names->data[i]);
+
+        psFits *maskFile = psFitsOpen(maskName, "w");
+        psFits *rejmapFile = psFitsOpen(rejmapName, "w");
+        psFits *gradFile = psFitsOpen(gradName, "w");
+        if (!psFitsWriteImage(maskFile, NULL, mask, 0, NULL)) {
+            psErrorStackPrint(stderr, "Unable to write image: %s\n", maskName);
+        }
+        psTrace("stac", 1, "Mask image written to %s\n", maskName);
+        if (!psFitsWriteImage(rejmapFile, NULL, rejmap, 0, NULL)) {
+            psErrorStackPrint(stderr, "Unable to write image: %s\n", rejmapName);
+        }
+        psTrace("stac", 1, "Rejection map written to %s\n", rejmapName);
+        if (!psFitsWriteImage(gradFile, NULL, gradient, 0, NULL)) {
+            psErrorStackPrint(stderr, "Unable to write image: %s\n", gradName);
+        }
+        psTrace("stac", 1, "Gradient image written to %s\n", gradName);
+        psFitsClose(maskFile);
+        psFitsClose(rejmapFile);
+        psFitsClose(gradFile);
+        psFree(rejmap);
+        psFree(gradient);
+#endif
+
+        // Stuff into the array
+        inputRej->data[i] = mask;
 
     }
Index: /trunk/stac/src/stacScales.c
===================================================================
--- /trunk/stac/src/stacScales.c	(revision 6886)
+++ /trunk/stac/src/stacScales.c	(revision 6887)
@@ -19,4 +19,6 @@
     psVector *values = psVectorAlloc(numSamples + 1, PS_TYPE_F32); // Vector containing sub-sample
     psVector *mask = psVectorAlloc(numSamples + 1, PS_TYPE_U8); // Mask for sample
+    values->n = numSamples + 1;
+    mask->n = numSamples + 1;
 
     int offset = 0;                     // Offset from start of the row
@@ -84,4 +86,5 @@
     } else {
         scales = psVectorAlloc(images->n, PS_TYPE_F32);
+        scales->n = images->n;
         *scalesPtr = scales;
     }
@@ -93,4 +96,5 @@
     } else {
         offsets = psVectorAlloc(images->n, PS_TYPE_F32);
+        offsets->n = images->n;
         *offsetsPtr = offsets;
     }
@@ -122,4 +126,5 @@
         // Transform the stellar positions to match the transformed reference frame
         psArray *starCoordsTransformed = psArrayAlloc(starCoords->n); // Transformed positions
+        starCoordsTransformed->n = starCoords->n;
         // Fix up difference between map and output frame
         starMap->x->coeff[0][0] -= xMapDiff;
@@ -132,4 +137,6 @@
         psArray *stars = psArrayAlloc(images->n); // Array of stellar photometry vectors
         psArray *masks = psArrayAlloc(images->n); // Array of masks for stars
+        stars->n = images->n;
+        masks->n = images->n;
 
         // Set scales relative to the first image
@@ -144,4 +151,6 @@
             psVector *photometry = psVectorAlloc(starCoords->n, PS_TYPE_F32); // Photometry of the stars
             psVector *mask = psVectorAlloc(starCoords->n, PS_TYPE_U8); // Mask for the photometry
+            photometry->n = starCoords->n;
+            mask->n = starCoords->n;
             for (int j = 0; j < starCoordsTransformed->n; j++) {
                 psPlane *coords = starCoordsTransformed->data[j]; // The coordinates of the star
Index: /trunk/stac/src/stacSize.c
===================================================================
--- /trunk/stac/src/stacSize.c	(revision 6886)
+++ /trunk/stac/src/stacSize.c	(revision 6887)
@@ -4,13 +4,17 @@
 #include "stac.h"
 
-bool stacSize(int *outnx, int *outny,	// Size of output image (returned)
-	      float *xMapDiff, float *yMapDiff, // Difference applied to maps
-	      const psArray *images,	// Input images
-	      psArray *maps		// Transformation maps
+bool stacSize(int *outnx, int *outny,   // Size of output image (returned)
+              float *xMapDiff, float *yMapDiff, // Difference applied to maps
+              const psVector *xSizes,   // Sizes of images in x
+              const psVector *ySizes,   // Sizes of images in y
+              psArray *maps             // Transformation maps
     )
 {
-    int nImages = images->n;		// Number of input images
+    int nImages = xSizes->n;            // Number of input images
+    assert(ySizes->n == nImages);
+    assert(xSizes->type.type == PS_TYPE_S32);
+    assert(ySizes->type.type == PS_TYPE_S32);
     assert(maps->n == nImages);
-    assert(outnx && outny);		// Must be able to write to these
+    assert(outnx && outny);             // Must be able to write to these
 
     psPlane *inCoord = psAlloc(sizeof(psPlane)); // Input coordinates
@@ -24,83 +28,84 @@
 
     for (int i = 0; i < nImages; i++) {
-	psImage *image = images->data[i]; // The image
-	psPlaneTransform *map = maps->data[i]; // The map
+        int xSize = xSizes->data.S32[i];// Size of image in x
+        int ySize = ySizes->data.S32[i];// Size of image in y
+        psPlaneTransform *map = maps->data[i]; // The map
 
-	psTrace("stac.size", 4, "Image %d:\n", i);
+        psTrace("stac.size", 4, "Image %d:\n", i);
 
-	// Lower left corner
-	inCoord->x = 0.0;
-	inCoord->y = 0.0;
-	(void)psPlaneTransformApply(outCoord, map, inCoord);
-	if (outCoord->x < xMin) {
-	    xMin = outCoord->x;
-	} else if (outCoord->x > xMax) {
-	    xMax = outCoord->x;
-	}
-	if (outCoord->y < yMin) {
-	    yMin = outCoord->y;
-	} else if (outCoord->y > yMax) {
-	    yMax = outCoord->y;
-	}
-	psTrace("stac.size", 4, "Lower left: %f %f\n", outCoord->x, outCoord->y);
+        // Lower left corner
+        inCoord->x = 0.0;
+        inCoord->y = 0.0;
+        (void)psPlaneTransformApply(outCoord, map, inCoord);
+        if (outCoord->x < xMin) {
+            xMin = outCoord->x;
+        } else if (outCoord->x > xMax) {
+            xMax = outCoord->x;
+        }
+        if (outCoord->y < yMin) {
+            yMin = outCoord->y;
+        } else if (outCoord->y > yMax) {
+            yMax = outCoord->y;
+        }
+        psTrace("stac.size", 4, "Lower left: %f %f\n", outCoord->x, outCoord->y);
 
-	// Lower right corner
-	inCoord->x = image->numCols;
-	(void)psPlaneTransformApply(outCoord, map, inCoord);
-	if (outCoord->x < xMin) {
-	    xMin = outCoord->x;
-	} else if (outCoord->x > xMax) {
-	    xMax = outCoord->x;
-	}
-	if (outCoord->y < yMin) {
-	    yMin = outCoord->y;
-	} else if (outCoord->y > yMax) {
-	    yMax = outCoord->y;
-	}
-	psTrace("stac.size", 4, "Lower right: %f %f\n", outCoord->x, outCoord->y);
+        // Lower right corner
+        inCoord->x = xSize;
+        (void)psPlaneTransformApply(outCoord, map, inCoord);
+        if (outCoord->x < xMin) {
+            xMin = outCoord->x;
+        } else if (outCoord->x > xMax) {
+            xMax = outCoord->x;
+        }
+        if (outCoord->y < yMin) {
+            yMin = outCoord->y;
+        } else if (outCoord->y > yMax) {
+            yMax = outCoord->y;
+        }
+        psTrace("stac.size", 4, "Lower right: %f %f\n", outCoord->x, outCoord->y);
 
-	// Upper right corner
-	inCoord->x = 0.0;
-	inCoord->y = image->numRows;
-	(void)psPlaneTransformApply(outCoord, map, inCoord);
-	if (outCoord->x < xMin) {
-	    xMin = outCoord->x;
-	} else if (outCoord->x > xMax) {
-	    xMax = outCoord->x;
-	}
-	if (outCoord->y < yMin) {
-	    yMin = outCoord->y;
-	} else if (outCoord->y > yMax) {
-	    yMax = outCoord->y;
-	}
-	psTrace("stac.size", 4, "Upper right: %f %f\n", outCoord->x, outCoord->y);
+        // Upper right corner
+        inCoord->x = 0.0;
+        inCoord->y = ySize;
+        (void)psPlaneTransformApply(outCoord, map, inCoord);
+        if (outCoord->x < xMin) {
+            xMin = outCoord->x;
+        } else if (outCoord->x > xMax) {
+            xMax = outCoord->x;
+        }
+        if (outCoord->y < yMin) {
+            yMin = outCoord->y;
+        } else if (outCoord->y > yMax) {
+            yMax = outCoord->y;
+        }
+        psTrace("stac.size", 4, "Upper right: %f %f\n", outCoord->x, outCoord->y);
 
-	// Upper left corner
-	inCoord->x = image->numCols;
-	(void)psPlaneTransformApply(outCoord, map, inCoord);
-	if (outCoord->x < xMin) {
-	    xMin = outCoord->x;
-	} else if (outCoord->x > xMax) {
-	    xMax = outCoord->x;
-	}
-	if (outCoord->y < yMin) {
-	    yMin = outCoord->y;
-	} else if (outCoord->y > yMax) {
-	    yMax = outCoord->y;
-	}
-	psTrace("stac.size", 4, "Upper left: %f %f\n", outCoord->x, outCoord->y);
+        // Upper left corner
+        inCoord->x = xSize;
+        (void)psPlaneTransformApply(outCoord, map, inCoord);
+        if (outCoord->x < xMin) {
+            xMin = outCoord->x;
+        } else if (outCoord->x > xMax) {
+            xMax = outCoord->x;
+        }
+        if (outCoord->y < yMin) {
+            yMin = outCoord->y;
+        } else if (outCoord->y > yMax) {
+            yMax = outCoord->y;
+        }
+        psTrace("stac.size", 4, "Upper left: %f %f\n", outCoord->x, outCoord->y);
     }
 
     // Tweak the maps to account for the offset
     if (xMapDiff != NULL) {
-	*xMapDiff = floorf(xMin);
+        *xMapDiff = floorf(xMin);
     }
     if (yMapDiff != NULL) {
-	*yMapDiff = floorf(yMin);
+        *yMapDiff = floorf(yMin);
     }
     for (int i = 0; i < nImages; i++) {
-	psPlaneTransform *map = maps->data[i]; // The map of interest
-	map->x->coeff[0][0] -= floorf(xMin);
-	map->y->coeff[0][0] -= floorf(yMin);
+        psPlaneTransform *map = maps->data[i]; // The map of interest
+        map->x->coeff[0][0] -= floorf(xMin);
+        map->y->coeff[0][0] -= floorf(yMin);
     }
 
Index: /trunk/stac/src/stacTransform.c
===================================================================
--- /trunk/stac/src/stacTransform.c	(revision 6886)
+++ /trunk/stac/src/stacTransform.c	(revision 6887)
@@ -127,4 +127,5 @@
     if (*outputs == NULL) {
         *outputs = psArrayAlloc(nImages);
+        (*outputs)->n = nImages;
         psTrace("stac.transform", 5, "Allocating space for transformed images, %dx%d\n", outnx, outny);
         for (int i = 0; i < nImages; i++) {
@@ -138,4 +139,5 @@
     if (errors && (*outErrors == NULL)) {
         *outErrors = psArrayAlloc(errors->n);
+        (*outErrors)->n = errors->n;
         psTrace("stac.transform", 5, "Allocating space for transformed error images, %dx%d\n", outnx, outny);
         for (int i = 0; i < nImages; i++) {
Index: /trunk/stac/src/sum.c
===================================================================
--- /trunk/stac/src/sum.c	(revision 6886)
+++ /trunk/stac/src/sum.c	(revision 6887)
@@ -10,4 +10,5 @@
     const char *outputName = argv[1];   // Output file name
     psArray *inputNames = psArrayAlloc(argc-2);
+    inputNames->n = argc - 2;
     for (int i = 2; i < argc; i++) {
         inputNames->data[i-2] = psStringCopy(argv[i]);