Index: trunk/pswarp/src/pswarpTransformReadout_Opt.c
===================================================================
--- trunk/pswarp/src/pswarpTransformReadout_Opt.c	(revision 10958)
+++ trunk/pswarp/src/pswarpTransformReadout_Opt.c	(revision 12505)
@@ -5,5 +5,5 @@
 
     int minX, minY, maxX, maxY;
-    int gridX, gridY, nextGridX, nextGridY;
+    int gridXo, gridX, gridY, nextGridX, nextGridY;
     pswarpMap *map = NULL;
 
@@ -19,5 +19,5 @@
 
     // we might want to do the rectangular regions outside of the selection independently
-    psImageInit (output->image, NAN);
+    // psImageInit (output->image, NAN);
 
     // find the output pixel range
@@ -31,5 +31,5 @@
 
     pswarpMapGridSetGrid (grid, minX, minY, &gridX, &gridY); 
-    pswarpMapGridNextGrid (grid, gridX, gridY, &nextGridX, &nextGridY);
+    pswarpMapGridNextGrid_Y (grid, gridY, &nextGridY);
     map = grid->maps[gridX][gridY];
 
@@ -38,11 +38,13 @@
 	if (y >= nextGridY) {
 	    gridY ++;
-	    pswarpMapGridNextGrid (grid, gridX, gridY, &nextGridX, &nextGridY);
+	    pswarpMapGridNextGrid_Y (grid, gridY, &nextGridY);
 	    map = grid->maps[gridX][gridY];
 	}
+
+	gridXo = gridX;
 	for (int x = minX; x < maxX; x++) {
 	    if (x >= nextGridX) {
 		gridX ++;
-		pswarpMapGridNextGrid (grid, gridX, gridY, &nextGridX, &nextGridY);
+		pswarpMapGridNextGrid_X (grid, gridX, &nextGridX);
 		map = grid->maps[gridX][gridY];
 	    }
@@ -60,4 +62,7 @@
 	    outData[y][x] = (psF32)psImagePixelInterpolate(inImage, inPix->x, inPix->y, NULL, 1, NAN, PS_INTERPOLATE_BILINEAR);
 	}
+	gridX = gridXo;
+	pswarpMapGridNextGrid_X (grid, gridX, &nextGridX);
+	map = grid->maps[gridX][gridY];
     }