Index: trunk/stac/src/stacTransform.c
===================================================================
--- trunk/stac/src/stacTransform.c	(revision 2751)
+++ 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;
 }