Index: trunk/stac/src/stacCombine.c
===================================================================
--- trunk/stac/src/stacCombine.c	(revision 2764)
+++ 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