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