Index: trunk/stac/src/stacRejection.c
===================================================================
--- trunk/stac/src/stacRejection.c	(revision 2751)
+++ trunk/stac/src/stacRejection.c	(revision 2764)
@@ -42,41 +42,9 @@
 }
 
-#if 0
-float stacGradient(psImage *image,	// Input for which to measure the gradient
-		   int x, int y		// Coordinates at which to measure the gradient
-		   )
-{
-    float sum = 0.0; 			// The sum of surrounding pixels
-    float maxDiff = 0.0;
-    int num = 0;
-    // Get limits
-    int xMin = MAX(x - 1, 0);
-    int xMax = MIN(x + 1, image->numCols - 1);
-    int yMin = MAX(y - 1, 0);
-    int yMax = MIN(y + 1, image->numRows - 1);
-    for (int j = yMin; j <= yMax; j++) {
-	for (int i = xMin; i <= xMax; i++) {
-	    if (image->data.F32[j][i] - image->data.F32[y][x] < maxDiff) {
-		maxDiff = image->data.F32[j][i] - image->data.F32[y][x];
-	    }
-	    sum += image->data.F32[j][i];
-	    num++;
-	}
-    }
-    sum -= image->data.F32[y][x];
-    sum /= (float)(num-1);
-    sum -= image->data.F32[y][x];
-//    return sum / image->data.F32[y][x];
-    return maxDiff / image->data.F32[y][x];
-}   
-#endif
-
-
 psArray *stacRejection(psArray *inputs,	// Input images
-		       psArray *transformed, // Transformed images
-		       psImage *combined, // Combined image
+		       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
-		       psArray *rejected, // Rejected images
 		       stacConfig *config // Configuration
     )
@@ -85,31 +53,34 @@
 
     // Check inputs
-    if (inputs->n != transformed->n) {
-	psError("stac.rejection", "Number of input (%d) and transformed (%d) images does not match.\n",
-		inputs->n, transformed->n);
+    if (inputs->n != rejected->n) {
+	psError("stac.rejection", "Number of input images (%d) and rejection masks (%d) does not match.\n",
+		inputs->n, rejected->n);
+	return NULL;
+    }
+    if (inputs->n != regions->n) {
+	psError("stac.rejection", "Number of input images (%d) and region masks (%d) does not match.\n",
+		inputs->n, regions->n);
+	return NULL;
+    }
+    if (inputs->n != maps->n) {
+	psError("stac.rejection", "Number of input images (%d) and transformation maps (%d) does "
+		"not match.\n", inputs->n, maps->n);
+	return NULL;
+    }
+    if (inputs->n != inverseMaps->n) {
+	psError("stac.rejection", "Number of input images (%d) and inverse transformation maps (%d) does "
+		"not match.\n", inputs->n, inverseMaps->n);
 	return NULL;
     }
     for (int i = 0; i < nImages; i++) {
-	if ((((psImage*)(transformed->data[i]))->numRows != combined->numRows) ||
-	    (((psImage*)(transformed->data[i]))->numCols != combined->numCols)) {
+	if ((((psImage*)(inputs->data[i]))->numRows != ((psImage*)(regions->data[i]))->numRows) ||
+	    (((psImage*)(inputs->data[i]))->numCols != ((psImage*)(regions->data[i]))->numCols)) {
 	    psError("stac.rejection",
-		    "Sizes of transformed (%dx%d) and combined (%dx%d) images do not match.\n",
-		    ((psImage*)(transformed->data[i]))->numCols, ((psImage*)(transformed->data[i]))->numRows,
-		    combined->numCols, combined->numRows);
+		    "Sizes of input image (%dx%d) and region mask (%dx%d) for input %d do not match.\n",
+		    ((psImage*)(inputs->data[i]))->numCols, ((psImage*)(inputs->data[i]))->numRows,
+		    ((psImage*)(regions->data[i]))->numCols, ((psImage*)(regions->data[i]))->numRows, i);
 	    return NULL;
 	}
-	if ((((psImage*)(transformed->data[i]))->numRows != ((psImage*)(rejected->data[i]))->numRows) ||
-	    (((psImage*)(transformed->data[i]))->numCols != ((psImage*)(rejected->data[i]))->numCols)) {
-	    psError("stac.rejection",
-		    "Sizes of transformed (%dx%d) and rejected (%dx%d) images do not match.\n",
-		    ((psImage*)(transformed->data[i]))->numCols, ((psImage*)(transformed->data[i]))->numRows,
-		    ((psImage*)(rejected->data[i]))->numCols, ((psImage*)(rejected->data[i]))->numRows);
-	    return NULL;
-	}
-    }
-
-    // Size of output images
-    int nxOutput = combined->numCols;
-    int nyOutput = combined->numRows;
+    }
 
     // Stuff for the transformations
@@ -132,4 +103,5 @@
 	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);
@@ -153,45 +125,40 @@
 	for (int y = 0; y < nyInput; y++) {
 	    for (int x = 0; x < nxInput; x++) {
-    		inCoords->x = (double)x;
-		inCoords->y = (double)y;
-		(void)psPlaneTransformApply(outCoords, map, inCoords);
-		float maskVal = (float)psImagePixelInterpolate(reject, outCoords->x, outCoords->y,
-							       NULL, 0, 0.0, PS_INTERPOLATE_BILINEAR);
+
+		// Only transform pixels of interest
+		if (region->data.U8[y][x]) {
+
+		    inCoords->x = (double)x;
+		    inCoords->y = (double)y;
+		    (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
-
-#ifdef TESTING
-		// Check gradient
-		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
-			    meanGrads += stacGradient(inputs->data[j], xPix, yPix);
-			    numGrads++;
+		    rejmap->data.F32[y][x] = maskVal;
+
+		    // Calculate gradient
+		    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
+				meanGrads += stacGradient(inputs->data[j], xPix, yPix);
+				numGrads++;
+			    }
 			}
 		    }
-		}
-		if (numGrads > 0) {
-		    meanGrads /= (float)numGrads;
-		} else {
-		    meanGrads = 0;
-		}
-#else
-		float meanGrads = 10000.0;
-#endif
-
-		
-#ifdef TESTING
-		grad->data.F32[y][x] = stacGradient(inputs->data[i], x, y) / meanGrads;
-#endif
-		if (maskVal > config->frac) {
+		    if (numGrads > 0) {
+			meanGrads /= (float)numGrads;
+		    } else {
+			meanGrads = 0;
+		    }
+		    grad->data.F32[y][x] = stacGradient(inputs->data[i], x, y) / meanGrads;
+#endif
 
 		    if ((maskVal > config->frac) &&
@@ -200,11 +167,13 @@
 			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
+			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;
 		    }
-		} else {
-		    mask->data.U8[y][x] = 0;
-		}
+
+		} // Only touching pixels of interest
 
 	    }