Index: trunk/pswarp/src/pswarpTransformReadout_Opt.c
===================================================================
--- trunk/pswarp/src/pswarpTransformReadout_Opt.c	(revision 10954)
+++ trunk/pswarp/src/pswarpTransformReadout_Opt.c	(revision 10957)
@@ -1,33 +1,20 @@
 # include "pswarp.h"
 
-bool pswarpTransformReadout (pmReadout *output, pmReadout *input, pmConfig *config) {
+// NOTE: in this function, the coordinates are transformed from the OUTPUT to the INPUT
+bool pswarpTransformReadout_Opt (pmReadout *output, pmReadout *input, pmConfig *config) {
+
+    int minX, minY, maxX, maxY;
+    int gridX, gridY, nextGridX, nextGridY;
+    pswarpMap *map = NULL;
 
     // XXX this implementation currently ignores the use of the region 
     psImage *region = NULL;
-    pmCell *cell = NULL;
 
     // select the current recipe
-    psMetadata *recipe  = psMetadataLookupPtr (NULL, config->recipes, PSPHOT_RECIPE);
-
-    int outNx = output->image->numCols;
-    int outNy = output->image->numRows;
+    // psMetadata *recipe  = psMetadataLookupPtr (NULL, config->recipes, PSPHOT_RECIPE);
 
     psPlane *inPix = psPlaneAlloc();	// Coordinates on the input detector
-    psPlane *outPix = psPlaneAlloc();	// Coordinates on the output detector
-
-    psPlane *FP = psPlaneAlloc();	// Coordinates on the focal plane
-    psPlane *TP = psPlaneAlloc();	// Coordinates on the tangent plane
-    psSphere *sky = psSphereAlloc();	// Coordinates on the sky
-
-    cell = input->parent;
-    pmChip *chipInput = cell->parent;
-    pmFPA *fpaInput = chipInput->parent;
-
-    cell = output->parent;
-    pmChip *chipOutput = cell->parent;
-    pmFPA *fpaOutput = chipOutput->parent;
-
+    psImage *inImage = input->image;
     psF32 **outData = output->image->data.F32;
-    psImage *inImage = input->image;
 
     // we might want to do the rectangular regions outside of the selection independently
@@ -37,50 +24,56 @@
     pswarpMatchRange (&minX, &minY, &maxX, &maxY, input, output);
 
+    pswarpMapGrid *grid = pswarpMapGridFromImage (input, output, 128, 128);
+
+    // XXX need to modify the grid based on this result and force the maxError < XXX
+    double maxError = pswarpMapGridMaxError (grid);
+    fprintf (stderr, "maximum error using this grid sampling: %f\n", maxError);
+
+    pswarpMapGridSetGrid (grid, minX, minY, &gridX, &gridY); 
+    pswarpMapGridNextGrid (grid, gridX, gridY, &nextGridX, &nextGridY);
+    map = grid->maps[gridX][gridY];
+
     // Iterate over the output image pixels
-    for (int y = 0; y < outNy; y++) {
-	for (int x = 0; x < outNx; x++) {
+    for (int y = minY; y < maxY; y++) {
+	if (y >= nextGridY) {
+	    gridY ++;
+	    pswarpMapGridNextGrid (grid, gridX, gridY, &nextGridX, &nextGridY);
+	    map = grid->maps[gridX][gridY];
+	}
+	for (int x = minX; x < maxX; x++) {
+	    if (x >= nextGridX) {
+		gridX ++;
+		pswarpMapGridNextGrid (grid, gridX, gridY, &nextGridX, &nextGridY);
+		map = grid->maps[gridX][gridY];
+	    }
+
 	    // Only transform those pixels requested
-	    if (!region || (region && region->data.U8[y][x])) {
-		// Transform!
+	    if (region && region->data.U8[y][x]) continue;
 
-		// XXX double check this 1/2 pixel offset
-		outPix->x = (double)x + 0.5;
-		outPix->y = (double)y + 0.5;
+	    // XXX double check this 1/2 pixel offset
+	    // XXX subtract 0.5,0.5 from result?
+	    pswarpMapApply (&inPix->x, &inPix->y, map, x + 0.5, y + 0.5);
 
-		psPlaneTransformApply(FP, chipOutput->toFPA, outPix);
-		psPlaneTransformApply (TP, fpaOutput->toTPA, FP);
-		psDeproject (sky, TP, fpaOutput->toSky);
-		
-		psProject (TP, sky, fpaInput->toSky);
-		psPlaneTransformApply (FP, fpaInput->fromTPA, TP);
-		psPlaneTransformApply (inPix, chipInput->fromFPA, FP);
-
-		// XXX get interpolation method from the recipe
-		outImage->data.F32[y][x] = (psF32)psImagePixelInterpolate(image, inPix->x, inPix->y, mask, 1, NAN, PS_INTERPOLATE_BILINEAR);
-
-# 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);
-		}
-# endif
-
-		outImage->data.F32[y][x] = (outImage->data.F32[y][x] - offset) / scale;
-# if (0)
-		if (error) {
-		    outError->data.F32[y][x] = outError->data.F32[y][x] / SQUARE(scale);
-		}
-# endif
-
-	    } // Pixels of interest
-
+	    // 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);
 	}
-    } // Iterating over output pixels
-
+    }
     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
+