Index: trunk/pswarp/src/pswarpTransformReadout_Opt.c
===================================================================
--- trunk/pswarp/src/pswarpTransformReadout_Opt.c	(revision 12505)
+++ trunk/pswarp/src/pswarpTransformReadout_Opt.c	(revision 12523)
@@ -5,5 +5,5 @@
 
     int minX, minY, maxX, maxY;
-    int gridXo, gridX, gridY, nextGridX, nextGridY;
+    int gridXo, gridX, gridY, nextGridXo, nextGridX, nextGridY;
     pswarpMap *map = NULL;
 
@@ -18,4 +18,11 @@
     psF32 **outData = output->image->data.F32;
 
+    // we need to apply the offset to convert parent coordinates to child coordinates for
+    // psImagePixelInterpolate below
+    int inCol0 = input->image->col0;
+    int inRow0 = input->image->row0;
+    int outCol0 = output->image->col0;
+    int outRow0 = output->image->row0;
+
     // we might want to do the rectangular regions outside of the selection independently
     // psImageInit (output->image, NAN);
@@ -24,4 +31,6 @@
     pswarpMatchRange (&minX, &minY, &maxX, &maxY, input, output);
 
+    // pswarpMapGridFromImage builds a set of locally-linear maps which convert the
+    // output coordinates to input coordinates
     pswarpMapGrid *grid = pswarpMapGridFromImage (input, output, 128, 128);
 
@@ -31,20 +40,30 @@
 
     pswarpMapGridSetGrid (grid, minX, minY, &gridX, &gridY); 
-    pswarpMapGridNextGrid_Y (grid, gridY, &nextGridY);
+    nextGridY = pswarpMapGridNextGrid_Y (grid, gridY);
     map = grid->maps[gridX][gridY];
 
-    // Iterate over the output image pixels
+    assert ((int)(minX - outCol0) >= 0);
+    assert ((int)(maxX - outCol0) <= output->image->numCols);
+    assert ((int)(minY - outRow0) >= 0);
+    assert ((int)(maxY - outRow0) <= output->image->numRows);
+
+    gridXo = gridX;
+    nextGridXo = pswarpMapGridNextGrid_X (grid, gridX);
+
+    // Iterate over the output image pixels (parent frame)
     for (int y = minY; y < maxY; y++) {
 	if (y >= nextGridY) {
 	    gridY ++;
-	    pswarpMapGridNextGrid_Y (grid, gridY, &nextGridY);
+	    nextGridY = pswarpMapGridNextGrid_Y (grid, gridY);
 	    map = grid->maps[gridX][gridY];
 	}
 
-	gridXo = gridX;
+	gridX = gridXo;
+	nextGridX = nextGridXo;
+	map = grid->maps[gridX][gridY];
 	for (int x = minX; x < maxX; x++) {
 	    if (x >= nextGridX) {
 		gridX ++;
-		pswarpMapGridNextGrid_X (grid, gridX, &nextGridX);
+		nextGridX = pswarpMapGridNextGrid_X (grid, gridX);
 		map = grid->maps[gridX][gridY];
 	    }
@@ -53,16 +72,19 @@
 	    if (region && region->data.U8[y][x]) continue;
 
-	    // XXX double check this 1/2 pixel offset
-	    // XXX subtract 0.5,0.5 from result?
+	    // pswarpMapApply converts the output coordinate (x,y) to the input coordinate.
+	    // both are in the parent frames of the input and output images.
 	    pswarpMapApply (&inPix->x, &inPix->y, map, x + 0.5, y + 0.5);
+
+	    if (inPix->x - inCol0 < 0) continue;
+	    if (inPix->x - inCol0 >= inImage->numCols) continue;
+	    if (inPix->y - inRow0 < 0) continue;
+	    if (inPix->y - inRow0 >= inImage->numRows) continue;
 
 	    // XXX get interpolation method from the recipe
 	    // XXX include mask
 	    // XXX apply scale and offset?
-	    outData[y][x] = (psF32)psImagePixelInterpolate(inImage, inPix->x, inPix->y, NULL, 1, NAN, PS_INTERPOLATE_BILINEAR);
+	    // psImagePixelInterpolate determines the value at pixel coordinate (x,y) in child coordinates 
+	    outData[y-outRow0][x-outCol0] = (psF32)psImagePixelInterpolate(inImage, inPix->x - inCol0, inPix->y - inRow0, NULL, 1, NAN, PS_INTERPOLATE_BILINEAR);
 	}
-	gridX = gridXo;
-	pswarpMapGridNextGrid_X (grid, gridX, &nextGridX);
-	map = grid->maps[gridX][gridY];
     }
 
@@ -71,17 +93,2 @@
     return true;
 }
-
-# if (0)
-    if (error) {
-	// Error is actually the variance
-	outError->data.F32[y][x] = (psF32)p_psImageErrorInterpolateBILINEAR_F32(error,
-										detector->x,
-										detector->y,
-										mask, 1, NAN);
-    }
-    if (error) {
-	outError->data.F32[y][x] = outError->data.F32[y][x] / SQUARE(scale);
-    }
-# endif
-
-