Index: /trunk/pswarp/src/Makefile.am
===================================================================
--- /trunk/pswarp/src/Makefile.am	(revision 10956)
+++ /trunk/pswarp/src/Makefile.am	(revision 10957)
@@ -18,6 +18,9 @@
 pswarpDefine.c             \
 pswarpErrorCodes.c         \
+pswarpMapGrid.c            \
+pswarpMatchRange.c         \
 pswarpParseCamera.c   	   \
 pswarpTransformReadout.c   \
+pswarpTransformReadout_Opt.c \
 pswarpVersion.c            
 
Index: /trunk/pswarp/src/pswarp.h
===================================================================
--- /trunk/pswarp/src/pswarp.h	(revision 10956)
+++ /trunk/pswarp/src/pswarp.h	(revision 10957)
@@ -13,13 +13,20 @@
 # define PSWARP_RECIPE "PSWARP" // Name of the recipe to use
 
+// a single pswarpMap converts coordinates from one image to a second image 
+// the linear model is only valid over a limited range of pixels
 typedef struct {
-  double Xo, Xx, Xy;
-  double Yo, Yx, Yy;
+    double Xo, Xx, Xy;
+    double Yo, Yx, Yy;
+    int xo;
+    int yo;
 } pswarpMap;
 
+// the pswarpMapGrid carries a collection of pswarpMag structures representing the
+// local value of the pswarpMap at different locations in the image.
 typedef struct {
-  pswarpMap **grid;
-  int Nx, Ny;
-  int nXpix, nYpix;
+    pswarpMap ***maps;
+    int nXpts, nYpts;			// number of x,y samples in the grid
+    int nXpix, nYpix;			// x,y spacing in src image pixels of grid samples
+    int xMin,  yMin;			// coordinate of first grid sample
 } pswarpMapGrid;
 
@@ -31,3 +38,15 @@
 void pswarpCleanup (pmConfig *config);
 bool pswarpTransformReadout (pmReadout *output, pmReadout *input, pmConfig *config);
+bool pswarpTransformReadout_Opt (pmReadout *output, pmReadout *input, pmConfig *config);
 
+bool pswarpMatchRange (int *minX, int *minY, int *maxX, int *maxY, pmReadout *dest, pmReadout *src);
+
+pswarpMap *pswarpMapAlloc ();
+pswarpMapGrid *pswarpMapGridAlloc (int Nx, int Ny);
+
+pswarpMapGrid *pswarpMapGridFromImage (pmReadout *dest, pmReadout *src, int nXpix, int nYpix);
+bool pswarpMapGridSetGrid (pswarpMapGrid *grid, int ix, int iy, int *gridX, int *gridY);
+bool pswarpMapGridNextGrid (pswarpMapGrid *grid, int gridX, int gridY, int *nextX, int *nextY);
+double pswarpMapGridMaxError (pswarpMapGrid *grid);
+bool pswarpMapApply (double *outX, double *outY, pswarpMap *map, double inX, double inY);
+bool pswarpMapSetLocalModel (pswarpMap *map, pmReadout *dest, pmReadout *src, int ix, int iy);
Index: /trunk/pswarp/src/pswarpMapGrid.c
===================================================================
--- /trunk/pswarp/src/pswarpMapGrid.c	(revision 10956)
+++ /trunk/pswarp/src/pswarpMapGrid.c	(revision 10957)
@@ -1,5 +1,153 @@
 # include "pswarp.h"
 
-static pswarpMapFree (pswarpMap *map) {
+// construct a grid with superpixel spacing of nXpix, nYpix
+// XXX for the moment, ignore readout->cell->chip offsets
+pswarpMapGrid *pswarpMapGridFromImage (pmReadout *dest, pmReadout *src, int nXpix, int nYpix) {
+
+    // split the difference of the remainder
+    int xMin = 0.5*(src->image->numCols % nXpix);
+    int yMin = 0.5*(src->image->numRows % nYpix);
+
+    pswarpMapGrid *grid = pswarpMapGridAlloc (src->image->numCols / nXpix, src->image->numRows / nYpix);
+
+    int ni = 0;
+    int nj = 0;
+    for (int i = xMin; i < src->image->numCols; i += nXpix, ni++) {
+	for (int j = yMin; j < src->image->numRows; j += nYpix, nj++) {
+	    pswarpMapSetLocalModel (grid->maps[ni][nj], dest, src, i, j);
+	}
+    }
+	    
+    grid->nXpix = nXpix;
+    grid->nYpix = nYpix;
+    grid->xMin = xMin;
+    grid->yMin = yMin;
+    return grid;
+}
+
+bool pswarpMapGridSetGrid (pswarpMapGrid *grid, int ix, int iy, int *gridX, int *gridY) {
+
+    *gridX = (ix - grid->xMin + 0.5*grid->nXpix) / grid->nXpix;
+    *gridY = (iy - grid->yMin + 0.5*grid->nYpix) / grid->nYpix;
+    return true;
+}
+
+bool pswarpMapGridNextGrid (pswarpMapGrid *grid, int gridX, int gridY, int *nextX, int *nextY) {
+
+    *nextX = gridX*grid->nXpix + grid->xMin + 0.5*grid->nXpix;
+    *nextY = gridY*grid->nYpix + grid->yMin + 0.5*grid->nYpix;
+    return true;
+}
+
+// measure the max error accumulated in appling one grid point to its neighbors
+double pswarpMapGridMaxError (pswarpMapGrid *grid) {
+
+    double xRaw, yRaw;
+    double xRef, yRef;
+    double maxError = 0;
+
+    for (int i = 0; i < grid->nXpts - 1; i++) {
+	for (int j = 0; j < grid->nYpts - 1; j++) {
+
+	    // measure the output coordinates for the next grid position using the current grid map
+	    // compare with the coordinates measured using the next grid map
+	    pswarpMapApply (&xRaw, &yRaw, grid->maps[i][j], grid->maps[i][j]->xo + grid->nXpix, grid->maps[i][j]->yo);
+	    pswarpMapApply (&xRef, &yRef, grid->maps[i+1][j], grid->maps[i][j]->xo + grid->nXpix, grid->maps[i][j]->yo);
+
+	    double posError = hypot (xRaw-xRef, yRaw-yRef);
+	    maxError = PS_MAX (maxError, posError);
+	}
+    }	
+    return maxError;
+}
+
+bool pswarpMapApply (double *outX, double *outY, pswarpMap *map, double inX, double inY) {
+
+    *outX = map->Xo + map->Xx*inX + map->Xy*inY;
+    *outY = map->Yo + map->Yx*inX + map->Yy*inY;
+
+    return true;
+}
+
+// determine the map for the given pixel from src to dest. pixel is in src coords
+bool pswarpMapSetLocalModel (pswarpMap *map, pmReadout *dest, pmReadout *src, int ix, int iy) {
+
+    pmCell *cell = NULL;
+
+    cell = src->parent;
+    pmChip *chipSrc = cell->parent;
+    pmFPA *fpaSrc = chipSrc->parent;
+
+    cell = dest->parent;
+    pmChip *chipDest = cell->parent;
+    pmFPA *fpaDest = chipDest->parent;
+
+    // XXX save these as static for speed?
+    psPlane *offset = psPlaneAlloc();
+
+    psPlane *FP = psPlaneAlloc();
+    psPlane *TP = psPlaneAlloc();
+    psSphere *sky = psSphereAlloc();
+
+    psPlane *V00 = psPlaneAlloc();
+    psPlane *V10 = psPlaneAlloc();
+    psPlane *V01 = psPlaneAlloc();
+
+    // XXX need to include readout->cell->chip offsets
+
+    // V(0,0) position
+    offset->x = ix;
+    offset->y = iy;
+    psPlaneTransformApply(FP, chipSrc->toFPA, offset);
+    psPlaneTransformApply (TP, fpaSrc->toTPA, FP);
+    psDeproject (sky, TP, fpaSrc->toSky);
+    psProject (TP, sky, fpaDest->toSky);
+    psPlaneTransformApply (FP, fpaDest->fromTPA, TP);
+    psPlaneTransformApply (V00, chipDest->fromFPA, FP);
+    
+    // V(1,0) position
+    offset->x = ix + 1;
+    offset->y = iy;
+    psPlaneTransformApply(FP, chipSrc->toFPA, offset);
+    psPlaneTransformApply (TP, fpaSrc->toTPA, FP);
+    psDeproject (sky, TP, fpaSrc->toSky);
+    psProject (TP, sky, fpaDest->toSky);
+    psPlaneTransformApply (FP, fpaDest->fromTPA, TP);
+    psPlaneTransformApply (V10, chipDest->fromFPA, FP);
+    
+    // V(0,1) position
+    offset->x = ix;
+    offset->y = iy + 1;
+    psPlaneTransformApply(FP, chipSrc->toFPA, offset);
+    psPlaneTransformApply (TP, fpaSrc->toTPA, FP);
+    psDeproject (sky, TP, fpaSrc->toSky);
+    psProject (TP, sky, fpaDest->toSky);
+    psPlaneTransformApply (FP, fpaDest->fromTPA, TP);
+    psPlaneTransformApply (V01, chipDest->fromFPA, FP);
+
+    map->Xx = V10->x - V00->x;
+    map->Xy = V01->x - V00->x;
+    map->Xo = V00->x + map->Xx*ix + map->Xy*iy;
+
+    map->Yx = V10->y - V00->y;
+    map->Yy = V01->y - V00->y;
+    map->Yo = V00->y + map->Yx*ix + map->Yy*iy;
+ 
+    map->xo = ix;
+    map->yo = iy;
+
+    psFree (offset);
+    psFree (FP);
+    psFree (TP);
+    psFree (sky);
+
+    psFree (V00);
+    psFree (V10);
+    psFree (V01);
+
+    return true;
+}
+
+static void pswarpMapFree (pswarpMap *map) {
   return;
 }
@@ -13,3 +161,38 @@
 }
 
+static void pswarpMapGridFree (pswarpMapGrid *grid) {
+
+    if (grid == NULL) return;
+    if (grid->maps == NULL) return;
+
+    for (int i = 0; i < grid->nXpts; i++) {
+	for (int j = 0; j < grid->nYpts; j++) {
+	    psFree (grid->maps[i][j]);
+	}
+	psFree (grid->maps[i]);
+    }
+    psFree (grid->maps);
+    return;
+}
+
+pswarpMapGrid *pswarpMapGridAlloc (int nXpts, int nYpts) {
+
+  pswarpMapGrid *grid = (pswarpMapGrid *) psAlloc (sizeof(pswarpMapGrid));
+  psMemSetDeallocator(grid, (psFreeFunc) pswarpMapGridFree);
+
+  grid->maps = psAlloc (nXpts*sizeof(void **));
+  for (int i = 0; i < nXpts; i++) {
+      grid->maps[i] = psAlloc (nYpts*sizeof(void *));
+      for (int j = 0; j < nYpts; j++) {
+	  grid->maps[i][j] = pswarpMapAlloc();
+      }
+  }
+  grid->nXpts = nXpts;
+  grid->nYpts = nYpts;
+
+  grid->nXpix = 0;
+  grid->nYpix = 0;
   
+  return grid;
+}
+
Index: /trunk/pswarp/src/pswarpMatchRange.c
===================================================================
--- /trunk/pswarp/src/pswarpMatchRange.c	(revision 10957)
+++ /trunk/pswarp/src/pswarpMatchRange.c	(revision 10957)
@@ -0,0 +1,79 @@
+# include "pswarp.h"
+
+# define TEST_MINMAX \
+    psPlaneTransformApply(FP, chipDest->toFPA, destPix); \
+    psPlaneTransformApply (TP, fpaDest->toTPA, FP); \
+    psDeproject (sky, TP, fpaDest->toSky); \
+    psProject (TP, sky, fpaSrc->toSky); \
+    psPlaneTransformApply (FP, fpaSrc->fromTPA, TP); \
+    psPlaneTransformApply (srcPix, chipSrc->fromFPA, FP); \
+    *minX = PS_MAX (*minX, srcPix->x); \
+    *minY = PS_MAX (*minY, srcPix->y); \
+    *maxX = PS_MIN (*maxX, srcPix->x); \
+    *maxY = PS_MIN (*maxY, srcPix->y);
+
+// we are warping from src to dest.  find the max overlapping pixels in the INPUT (src)
+// coordinate frame
+bool pswarpMatchRange (int *minX, int *minY, int *maxX, int *maxY, pmReadout *dest, pmReadout *src) {
+
+    // transform input corners and edge centers to output coords
+    // find max overlapping region of output image
+    
+    pmCell *cell = NULL;
+
+    cell = src->parent;
+    pmChip *chipSrc = cell->parent;
+    pmFPA *fpaSrc = chipSrc->parent;
+
+    cell = dest->parent;
+    pmChip *chipDest = cell->parent;
+    pmFPA *fpaDest = chipDest->parent;
+
+    *minX = 0;
+    *minY = 0;
+    *maxX = src->image->numCols;
+    *maxY = src->image->numRows;
+
+    // XXX save these as static for speed?
+    psPlane *srcPix = psPlaneAlloc();
+    psPlane *destPix = psPlaneAlloc();
+
+    psPlane *FP = psPlaneAlloc();
+    psPlane *TP = psPlaneAlloc();
+    psSphere *sky = psSphereAlloc();
+
+    destPix->x = 0;
+    destPix->y = 0;
+    TEST_MINMAX;
+    
+    destPix->x = dest->image->numCols;
+    destPix->y = 0;
+    TEST_MINMAX;
+    
+    destPix->x = dest->image->numCols;
+    destPix->y = dest->image->numRows;
+    TEST_MINMAX;
+    
+    destPix->x = 0;
+    destPix->y = dest->image->numRows;
+    TEST_MINMAX;
+    
+    destPix->x = 0.5*dest->image->numCols;
+    destPix->y = 0.0;
+    TEST_MINMAX;
+    
+    destPix->x = 0.0;
+    destPix->y = 0.5*dest->image->numRows;
+    TEST_MINMAX;
+    
+    destPix->x = 0.5*dest->image->numCols;
+    destPix->y = dest->image->numRows;
+    TEST_MINMAX;
+
+    destPix->x = dest->image->numCols;
+    destPix->y = 0.5*dest->image->numRows;
+    TEST_MINMAX;
+
+    return true;
+}
+
Index: /trunk/pswarp/src/pswarpTransformReadout.c
===================================================================
--- /trunk/pswarp/src/pswarpTransformReadout.c	(revision 10956)
+++ /trunk/pswarp/src/pswarpTransformReadout.c	(revision 10957)
@@ -35,45 +35,37 @@
 	for (int x = 0; x < outNx; x++) {
 	    // 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
+	    outPix->x = (double)x + 0.5;
+	    outPix->y = (double)y + 0.5;
 
-		psPlaneTransformApply(FP, chipOutput->toFPA, outPix);
-		psPlaneTransformApply (TP, fpaOutput->toTPA, FP);
-		psDeproject (sky, TP, fpaOutput->toSky);
+	    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);
+	    psProject (TP, sky, fpaInput->toSky);
+	    psPlaneTransformApply (FP, fpaInput->fromTPA, TP);
+	    psPlaneTransformApply (inPix, chipInput->fromFPA, FP);
 
-		// XXX get interpolation method from the recipe
-		outData[y][x] = (psF32)psImagePixelInterpolate(inImage, inPix->x, inPix->y, NULL, 1, NAN, PS_INTERPOLATE_BILINEAR);
-		// modify zero and scale?
-
-		// outData[y][x] = (psF32)psImagePixelInterpolate(inImage, 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);
-		}
-		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
+	    outData[y][x] = (psF32)psImagePixelInterpolate(inImage, inPix->x, inPix->y, NULL, 1, NAN, PS_INTERPOLATE_BILINEAR);
+	    // outData[y][x] = (psF32)psImagePixelInterpolate(inImage, inPix->x, inPix->y, mask, 1, NAN, PS_INTERPOLATE_BILINEAR);
+	    // modify zero and scale?
 	}
-    } // 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
Index: /trunk/pswarp/src/pswarpTransformReadout_Opt.c
===================================================================
--- /trunk/pswarp/src/pswarpTransformReadout_Opt.c	(revision 10956)
+++ /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
+