Index: /trunk/stac/src/stac.c
===================================================================
--- /trunk/stac/src/stac.c	(revision 2782)
+++ /trunk/stac/src/stac.c	(revision 2783)
@@ -36,5 +36,5 @@
     psTraceSetLevel("stac.rejection",10);
     psTraceSetLevel("stac.time",10);
-    
+    psTraceSetLevel("stac.area",10);
 
     // Set logging level
@@ -50,4 +50,6 @@
     psArray *maps = stacReadMaps(config);
 
+    psTrace("stac.time", 1, "I/O completed at %f seconds\n", getTime() - startTime);
+
     // Invert maps
     psArray *inverseMaps = stacInvertMaps(maps);
@@ -57,6 +59,7 @@
 
     // Transform inputs and errors
+    psArray *transformed = NULL;
     psArray *transformedErrors = NULL;
-    psArray *transformed = stacTransform(inputs, inverseMaps, errors, &transformedErrors, NULL, config);
+    (void)stacTransform(&transformed, &transformedErrors, inputs, inverseMaps, errors, NULL, NULL, config);
 
     psTrace("stac.time",1,"Transformation completed at %f seconds\n", getTime() - startTime);
@@ -64,5 +67,8 @@
     // Combine with rejection
     psArray *rejected = NULL;
-    psImage *combined = stacCombine(transformed, transformedErrors, 2,  &rejected, config);
+    psImage *combined = NULL;
+    (void)stacCombine(&combined, &rejected, transformed, transformedErrors, config->nReject, NULL, config);
+
+    psTrace("stac.time",1,"First combination completed at %f seconds\n", getTime() - startTime);
 
 #ifdef TESTING
@@ -75,5 +81,4 @@
     psArray *regions = psArrayAlloc(inputs->n); // Array of images denoting regions of interest
     for (int i = 0; i < inputs->n; i++) {
-	fprintf(stderr, "Finding area of interest for image %d\n", i);
 	regions->data[i] = stacAreaOfInterest(rejected->data[i], inverseMaps->data[i],
 					      ((psImage*)(inputs->data[i]))->numCols,
@@ -89,16 +94,26 @@
     psArray *rejectedSource = stacRejection(inputs, rejected, regions, maps, inverseMaps, config);
 
-    // Redo transformation with the mask
-    psFree(transformed);
-    psFree(transformedErrors);
-    psArray *transformedErrorsNew = NULL;
-    psArray *transformedNew = stacTransform(inputs, inverseMaps, errors, &transformedErrorsNew,
-					    rejectedSource, config);
+    // Get regions of interest in the output frame
+    psImage *combineRegion = psImageAlloc(config->outnx, config->outny, PS_TYPE_U8);
+    for (int y = 0; y < config->outny; y++) {
+	for (int x = 0; x < config->outnx; x++) {
+	    combineRegion->data.U8[y][x] = 0;
+	}
+    }
+    for (int i = 0; i < inputs->n; i++) {
+	psImage *region = stacAreaOfInterest(rejectedSource->data[i], maps->data[i], config->outnx,
+					     config->outny);
+	psBinaryOp(combineRegion, combineRegion, "+", region);
+	psFree(region);
+    }
+
+    // Redo transformation with the masks
+    (void)stacTransform(&transformed, &transformedErrors, inputs, inverseMaps, errors, rejectedSource,
+			combineRegion, config);
 
     // Combine the newly-transformed CR-free images, no rejection
-    psFree(combined);
     psFree(rejected);
     rejected = NULL;
-    combined = stacCombine(transformedNew, transformedErrorsNew, 0,  &rejected, config);
+    (void)stacCombine(&combined, &rejected, transformed, transformedErrors, 0, combineRegion, config);
 
     // Write output image
@@ -107,4 +122,5 @@
     // Free everything I've used
     stacConfigFree(config);
+    psFree(combineRegion);
     psFree(regions);
     psFree(inputs);
@@ -113,9 +129,9 @@
     psFree(errors);
     psFree(transformedErrors);
-    psFree(transformedNew);
+    psFree(transformed);
     psFree(rejectedSource);
     psFree(combined);
 
-    psTrace("stac.time",1,"Combination completed at %f seconds\n", getTime() - startTime);
+    psTrace("stac.time",1,"Final combination completed at %f seconds\n", getTime() - startTime);
 
    // Check memory for leaks, corruption
Index: /trunk/stac/src/stac.h
===================================================================
--- /trunk/stac/src/stac.h	(revision 2782)
+++ /trunk/stac/src/stac.h	(revision 2783)
@@ -28,4 +28,5 @@
 					// considered bad
     float grad;				// Multiplier of the gradient
+    int nReject;			// Number of rejection iterations
 } stacConfig;
 
@@ -67,11 +68,13 @@
 // stacTransform.c
 
-// Transform input images
-psArray *stacTransform(const psArray *images, // Array of images to be transformed
-		       const psArray *maps, // Array of polynomials that do the transformation
-		       const psArray *errors, // Array of error images
-		       psArray **outErrors, // Array of output errors
-		       const psArray *masks, // Masks of input images
-		       const stacConfig *config	// Configuration
+// 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 stacConfig *config // Configuration
     );
 
@@ -110,9 +113,11 @@
 
 // Combine the transformed images
-psImage *stacCombine(psArray *images,	// Array of transformed images
-		     psArray *errors,	// Array of transformed error images
-		     int nReject,	// Number of rejection iterations
-		     psArray **rejected, // Array of rejection masks
-		     stacConfig *config	// Configuration
+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
+		 stacConfig *config	// Configuration
     );
 
Index: /trunk/stac/src/stacAreaOfInterest.c
===================================================================
--- /trunk/stac/src/stacAreaOfInterest.c	(revision 2782)
+++ /trunk/stac/src/stacAreaOfInterest.c	(revision 2783)
@@ -89,4 +89,7 @@
     }
 
+    psTrace("stac.area", 3, "Finding area of interest....\n");
+
+
     // Coordinates for the transformations
     psPlane *inFrame = psAlloc(sizeof(psPlane));
Index: /trunk/stac/src/stacCombine.c
===================================================================
--- /trunk/stac/src/stacCombine.c	(revision 2782)
+++ /trunk/stac/src/stacCombine.c	(revision 2783)
@@ -26,6 +26,7 @@
     for (int i = 0; i < num; i++) {
 	if (! masks->data.U8[i]) {
-	    sum += values->data.F32[i] / SQUARE(errors->data.F32[i]);
-	    weights += 1.0 / SQUARE(errors->data.F32[i]);
+	    // "error" here is the variance
+	    sum += values->data.F32[i] / errors->data.F32[i];
+	    weights += 1.0 / errors->data.F32[i];
 	}
     }
@@ -53,9 +54,12 @@
 
 
-psImage *stacCombine(psArray *images,	// Array of transformed images
-		     psArray *errors,	// Array of transformed error images
-		     int nReject,	// Number of rejection iterations
-		     psArray **rejected, // Array of rejection masks
-		     stacConfig *config	// Configuration
+
+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
+		 stacConfig *config	// Configuration
     )
 {
@@ -76,5 +80,5 @@
 		    "Image size mismatch on error image %d\nExpected %dx%d, got %dx%d\n",
 		    i, numCols, numRows, image->numCols, image->numRows);
-	    return NULL;
+	    return false;
 	}
 	if ((error->numCols != numCols) || (error->numRows != numRows)) {
@@ -82,16 +86,29 @@
 		    "Image size mismatch on error image %d\nExpected %dx%d, got %dx%d\n",
 		    i, numCols, numRows, error->numCols, error->numRows);
-	    return NULL;
-	}
+	    return false;
+	}
+    }
+
+    // Check combined image
+    if (*combined == NULL) {
+	*combined = psImageAlloc(numCols, numRows, PS_TYPE_F32); // Combined image
+    } else if (((*combined)->numRows != numRows) || ((*combined)->numCols != numCols)) {
+	psError("stac.combine", "Size of combined image (%dx%d) does not match inputs (%dx%d)\n",
+		(*combined)->numCols, (*combined)->numRows, numCols, numRows);
+	return false;
+    }
+
+    // Check area of interest
+    if (region && ((region->numRows != numRows) || (region->numCols != numCols))) {
+	psError("stac.combine", "Size of area of interest (%dx%d) does not match inputs (%dx%d)\n",
+		region->numCols, region->numRows, numCols, numRows);
+	return false;
     }
 
     psTrace("stac.combine", 1, "Combining images....\n");
-    psTrace("stac.combine", 3, "Saturation: %f Bad: %f\n", saturated, bad);
 
     psVector *pixels = psVectorAlloc(nImages, PS_TYPE_F32); // Will hold the pixels in the statistics step
     psVector *deltas = psVectorAlloc(nImages, PS_TYPE_F32); // Will hold the errors in the statistics step
-    psVector *mask = psVectorAlloc(nImages, PS_TYPE_U8); // Mask bad pixels
-    psImage *combined = psImageAlloc(numCols, numRows, PS_TYPE_F32); // Combined image
-
+    psVector *mask = psVectorAlloc(nImages, PS_TYPE_U8); // Will hold the mask in the statistics step
 
     // Set up rejection masks
@@ -105,7 +122,8 @@
 	    return NULL;
 	}
+
+	// Create and initialise rejection masks
 	for (int i = 0; i < nImages; i++) {
-	    (*rejected)->data[i] = psImageAlloc(numCols, numRows, PS_TYPE_U8); // Create rejection mask
-	    // Zero out the mask
+	    (*rejected)->data[i] = psImageAlloc(numCols, numRows, PS_TYPE_U8);
 	    for (int r = 0; r < numRows; r++) {
 		for (int c = 0; c < numCols; c++) {
@@ -124,56 +142,63 @@
     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];
-		pixels->data.F32[i] = pixel;
-		deltas->data.F32[i] = delta;
-		if ((pixel >= saturated) || (pixel <= bad) || (! isfinite(pixel)) || (! isfinite(delta))) {
-		    mask->data.U8[i] = (psU8)1; // Don't use!
-		} else {
-		    mask->data.U8[i] = (psU8)0; // Use.
-		}
-	    }
+		// 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) || (pixel <= bad) || (! 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
+		float average = stacCombineMean(pixels, deltas, mask); // Combined value
 
 #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
+		// 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;
+		    }
+		}
 	    
-	    // 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) / deltas->data.F32[i] > max)) {
-			max = (pixels->data.F32[i] - average) / deltas->data.F32[i];
-			maxIndex = i;
-		    }
-		}
-		// Reject the pixel with the maximum deviation
-		if (max > reject) {
-		    mask->data.U8[maxIndex] = 1;
-		    ((psImage*)((*rejected)->data[maxIndex]))->data.U8[y][x] += 1.0;
-		    // Re-do combination following rejection
-		    average = stacCombineMean(pixels, deltas, mask);
-		} else {
-		    keepGoing = false;
-		}
-	    }
-	    
-	    combined->data.F32[y][x] = average;
+		(*combined)->data.F32[y][x] = average;
+
+	    } // Pixels of interest
+
 	}
     } // Iterating over output pixels
Index: /trunk/stac/src/stacConfig.c
===================================================================
--- /trunk/stac/src/stacConfig.c	(revision 2782)
+++ /trunk/stac/src/stacConfig.c	(revision 2783)
@@ -9,5 +9,5 @@
 {
     fprintf (stderr, "STAC: Simultaneous Telescope Array Combination\n"
-	     "Usage: %s [-h] [-v] [-g GAIN] [-r RN] [-o NX NY] [-s SAT] [-b BAD] [-k REJ] [-k FRAC] [-G GRAD] OUT IN1 IN2...\n"
+	     "Usage: %s [-h] [-v] [-g GAIN] [-r RN] [-o NX NY] [-s SAT] [-b BAD] [-k REJ] [-n NREJECT] [-k FRAC] [-G GRAD] OUT IN1 IN2...\n"
 	     "where\n"
 	     "\t-h           Help (this info)\n"
@@ -19,4 +19,5 @@
 	     "\t-b BAD       Bad level (0)\n"
 	     "\t-k REJ       Rejection level (k-sigma; 3.5)\n"
+	     "\t-n NREJECT   Number of rejection iterations (2)\n"
 	     "\t-f FRAC      Fraction of pixel to be marked before considered bad (0.5)\n"
 	     "\t-d GRAD      Gradient threshold for pixel to be marked (0.0)\n"
@@ -45,4 +46,5 @@
     config->frac = 0.5;
     config->grad = 0.4;
+    config->nReject = 2;
 
     return config;
@@ -74,5 +76,5 @@
 
     /* Parse command-line arguments using getopt */
-    while ((opt = getopt(argc, argv, "hvg:r:o:s:b:k:f:G:")) != -1) {
+    while ((opt = getopt(argc, argv, "hvg:r:o:s:b:k:n:f:G:")) != -1) {
         switch (opt) {
 	  case 'h':
@@ -115,4 +117,9 @@
 	  case 'k':
 	    if (sscanf(optarg, "%f", &config->reject) != 1) {
+		help(programName);
+	    }
+	    break;
+	  case 'n':
+	    if (sscanf(optarg, "%d", &config->nReject) != 1) {
 		help(programName);
 	    }
Index: /trunk/stac/src/stacTransform.c
===================================================================
--- /trunk/stac/src/stacTransform.c	(revision 2782)
+++ /trunk/stac/src/stacTransform.c	(revision 2783)
@@ -105,10 +105,12 @@
 
 
-psArray *stacTransform(const psArray *images, // Array of images to be transformed
-		       const psArray *maps, // Array of polynomials that do the transformation
-		       const psArray *errors, // Array of error images to be transformed
-		       psArray **outErrors, // Transformed error images for output
-		       const psArray *masks, // Masks of input images
-		       const stacConfig *config	// Configuration
+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 stacConfig *config // Configuration
     )
 {
@@ -121,26 +123,44 @@
 	psError("stac.transform", "Number of maps (%d) does not match number of images (%d).\n",
 		maps->n, images->n);
-	return NULL;
+	return false;
     }
     if (errors && (images->n != errors->n)) {
 	psError("stac.transform", "Number of error images (%d) does not match number of images (%d).\n",
 		errors->n, images->n);
-	return NULL;
-    }
-
+	return false;
+    }
+
+    // Allocate the output images if required, otherwise check the number
+    if (*outputs == NULL) {
+	*outputs = psArrayAlloc(nImages);
+	psTrace("stac.transform", 5, "Allocating space for transformed images, %dx%d\n", nx, ny);
+	for (int i = 0; i < nImages; i++) {
+	    (*outputs)->data[i] = psImageAlloc(nx, ny, PS_TYPE_F32);
+	}
+    } else if ((*outputs)->n != nImages) {
+	psError("stac.transform", "Number of output images (%d) does not match number of input images "
+		"(%d).\n", (*outputs)->n, nImages);
+	return false;
+    }
+
+    // Allocate the output error images, if required, otherwise check the number
     if (errors && (*outErrors == NULL)) {
-	// Allocate the output error images, if required
 	*outErrors = psArrayAlloc(errors->n);
+	psTrace("stac.transform", 5, "Allocating space for transformed error images, %dx%d\n", nx, ny);
+	for (int i = 0; i < nImages; i++) {
+	    (*outErrors)->data[i] = psImageAlloc(nx, ny, PS_TYPE_F32);
+	}
     } else if (errors->n != (*outErrors)->n) {
 	psError("stac.transform", "Number of error output images (%d) does not match number of input"
 		"images (%d).\n", (*outErrors)->n, errors->n);
-	return NULL;
-    }
-
+	return false;
+    }
+
+    // Check the masks, if specified
     if (masks != NULL) {
 	if (masks->n != nImages) {
 	    psError("stac.transform", "Number of masks (%d) does not match number of input images (%d).\n",
 		    masks->n, nImages);
-	    return NULL;
+	    return false;
 	}
 	for (int i = 0; i < nImages; i++) {
@@ -151,15 +171,14 @@
 			 "Size of input mask (%dx%d) does not match that of input image (%dx%d).\n",
 			 mask->numCols, mask->numRows, image->numCols, image->numRows);
-		return NULL;
+		return false;
 	    }
 	}
     }
 
-    // Arrays of images for return
-    psArray *outImages = psArrayAlloc(nImages);	// Output images
 
     // Stuff for the transformations
     psPlane *detector = psAlloc(sizeof(psPlane)); // Coordinates on the detector
     psPlane *sky = psAlloc(sizeof(psPlane)); // Coordinates on the sky
+
 
     // Iterate over the images
@@ -171,9 +190,6 @@
 	psPlaneTransform *map = maps->data[n]; // The map
 	psImage *error = errors->data[n]; // The error image
-
-	// Initialise the output images
-	psImage *outImage = psImageAlloc(nx, ny, PS_TYPE_F32);
-	psImage *outError = psImageAlloc(nx, ny, PS_TYPE_F32);
-	psTrace("stac.transform", 5, "Allocating space for transformed image, %dx%d\n", nx, ny);
+	psImage *outImage = (*outputs)->data[n]; // The output image
+	psImage *outError = (*outErrors)->data[n]; // The output error image
 
 #if 0
@@ -193,29 +209,26 @@
 	}
 
-#if 0
-	// Calculate scales; could be adapted for higher order than linear by moving this into the
-	// pixel-dependent section and calculating derivatives
-	double xscale = 0.5*sqrt(SQUARE(map->x->coeff[0][1]) + SQUARE(map->x->coeff[1][0]));
-	double yscale = 0.5*sqrt(SQUARE(map->y->coeff[0][1]) + SQUARE(map->y->coeff[1][0]));
-	double areascale = 0.25 / xscale / yscale;
-#endif
-
 	// Iterate over the output image pixels
 	for (int y = 0; y < ny; y++) {
 	    for (int x = 0; x < nx; x++) {
-		// Transform!
-		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 + 0.5,
-									  detector->y + 0.5, mask, 1, 0.0,
-									  PS_INTERPOLATE_BILINEAR);
-		outError->data.F32[y][x] = (psF32)p_psImageErrorInterpolateBILINEAR_F32(error,
-											detector->x + 0.5,
-											detector->y + 0.5,
-											mask, 1, 0.0);
-		outError->data.F32[y][x] = sqrtf(outError->data.F32[y][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 + 0.5,
+									      detector->y + 0.5, mask, 1, 0.0,
+									      PS_INTERPOLATE_BILINEAR);
+		    // Error is actually the variance
+		    outError->data.F32[y][x] = (psF32)p_psImageErrorInterpolateBILINEAR_F32(error,
+											    detector->x + 0.5,
+											    detector->y + 0.5,
+											    mask, 1, 0.0);
+		} // Pixels of interest
+
 	    }
 	} // Iterating over output pixels
@@ -233,8 +246,4 @@
 #endif
 
-	// Plug the new images into arrays
-	outImages->data[n] = outImage;
-	(*outErrors)->data[n] = outError;
-
     } // Iterating over images
 
@@ -243,5 +252,5 @@
     psFree(sky);
 
-    return outImages;
+    return true;
 }