Index: /trunk/pswarp/src/pswarp.h
===================================================================
--- /trunk/pswarp/src/pswarp.h	(revision 12522)
+++ /trunk/pswarp/src/pswarp.h	(revision 12523)
@@ -48,6 +48,6 @@
 pswarpMapGrid *pswarpMapGridFromImage (pmReadout *dest, pmReadout *src, int nXpix, int nYpix);
 bool pswarpMapGridSetGrid (pswarpMapGrid *grid, int ix, int iy, int *gridX, int *gridY);
-bool pswarpMapGridNextGrid_X (pswarpMapGrid *grid, int gridX, int *nextX);
-bool pswarpMapGridNextGrid_Y (pswarpMapGrid *grid, int gridY, int *nextY);
+int pswarpMapGridNextGrid_X (pswarpMapGrid *grid, int gridX);
+int pswarpMapGridNextGrid_Y (pswarpMapGrid *grid, int gridY);
 double pswarpMapGridMaxError (pswarpMapGrid *grid);
 bool pswarpMapApply (double *outX, double *outY, pswarpMap *map, double inX, double inY);
Index: /trunk/pswarp/src/pswarpDataLoad.c
===================================================================
--- /trunk/pswarp/src/pswarpDataLoad.c	(revision 12522)
+++ /trunk/pswarp/src/pswarpDataLoad.c	(revision 12523)
@@ -9,23 +9,8 @@
 bool pswarpDataLoad (pmConfig *config) {
 
-    bool status;
     pmChip *chip;
     pmCell *cell;
     pmReadout *readout;
     pmFPAview *view;
-
-    // select the current recipe
-    psMetadata *recipe  = psMetadataLookupPtr (NULL, config->recipes, PSASTRO_RECIPE);
-    if (!recipe) {
-	psError(PSWARP_ERR_CONFIG, false, "Can't find PSASTRO recipe needed for pixel scale!\n");
-	return false;
-    }
-
-    // physical pixel scale in microns per pixel
-    double pixelScale = psMetadataLookupF32 (&status, recipe, "PSASTRO.PIXEL.SCALE");
-    if (!status) {
-	psError(PS_ERR_IO, false, "Failed to lookup pixel scale"); 
-	return false; 
-    } 
 
     // select the input data sources
@@ -61,18 +46,18 @@
     view = pmFPAviewAlloc (0);
 
+    // XXX need to read only the headers for the skycell
+    // XXX these pmAstromReadBilevel functions seem to be broken
+
     // find the FPA phu
-    // XXX wrap the test below into a function
-    // XXX need to read only the headers for the skycell
-    // XXX need to optionally load the astrometry datafile
     bool bilevelAstrometry = false;
     pmHDU *phu = pmFPAviewThisPHU (view, astrom->fpa);
     if (phu) {
-      char *ctype = psMetadataLookupStr (NULL, phu->header, "CTYPE1");
-      if (ctype) {
-	bilevelAstrometry = !strcmp (&ctype[4], "-DIS");
-      }
+	char *ctype = psMetadataLookupStr (NULL, phu->header, "CTYPE1");
+	if (ctype) {
+	    bilevelAstrometry = !strcmp (&ctype[4], "-DIS");
+	}
     }
     if (bilevelAstrometry) {
-      pmAstromReadBilevelMosaic (input->fpa, phu->header, pixelScale);
+	pmAstromReadBilevelMosaic (input->fpa, phu->header);
     } 
 
@@ -88,7 +73,8 @@
 	pmHDU *hdu = pmFPAviewThisHDU (view, astrom->fpa);
 	if (bilevelAstrometry) {
-	  pmAstromReadBilevelChip (chip, hdu->header); 
+	    pmAstromReadBilevelChip (chip, hdu->header); 
 	} else {
-	  pmAstromReadWCS (input->fpa, chip, hdu->header, pixelScale);
+	    // we use a default FPA pixel scale of 1.0
+	    pmAstromReadWCS (input->fpa, chip, hdu->header, 1.0);
 	}
 
Index: /trunk/pswarp/src/pswarpDefine.c
===================================================================
--- /trunk/pswarp/src/pswarpDefine.c	(revision 12522)
+++ /trunk/pswarp/src/pswarpDefine.c	(revision 12523)
@@ -4,22 +4,6 @@
 bool pswarpDefine (pmConfig *config) {
 
-    bool status;
-    psMetadata *recipe = NULL;
-
-    // load the pixel scale from the PSASTRO recipe
-    recipe  = psMetadataLookupPtr (NULL, config->recipes, PSASTRO_RECIPE);
-    if (!recipe) {
-	psError(PSWARP_ERR_CONFIG, false, "Can't find PSASTRO recipe needed for pixel scale!\n");
-	return false;
-    }
-    // physical pixel scale in microns per pixel
-    double pixelScale = psMetadataLookupF32 (&status, recipe, "PSASTRO.PIXEL.SCALE");
-    if (!status) {
-	psError(PS_ERR_IO, false, "Failed to lookup pixel scale"); 
-	return false; 
-    } 
-
     // load the PSWARP recipe
-    recipe  = psMetadataLookupPtr (NULL, config->recipes, PSWARP_RECIPE);
+    psMetadata *recipe = psMetadataLookupPtr (NULL, config->recipes, PSWARP_RECIPE);
     if (!recipe) {
 	psError(PSWARP_ERR_CONFIG, true, "Can't find PSWARP recipe!\n");
@@ -46,14 +30,15 @@
     bool bilevelAstrometry = false;
     if (phu) {
-      char *ctype = psMetadataLookupStr (NULL, phu->header, "CTYPE1");
-      if (ctype) {
-	bilevelAstrometry = !strcmp (&ctype[4], "-DIS");
-      }
+	char *ctype = psMetadataLookupStr (NULL, phu->header, "CTYPE1");
+	if (ctype) {
+	    bilevelAstrometry = !strcmp (&ctype[4], "-DIS");
+	}
     }
     if (bilevelAstrometry) {
-      pmAstromReadBilevelMosaic (skycell->fpa, phu->header, pixelScale);
-      pmAstromReadBilevelChip (chip, hdu->header); 
+	pmAstromReadBilevelMosaic (skycell->fpa, phu->header);
+	pmAstromReadBilevelChip (chip, hdu->header); 
     } else {
-      pmAstromReadWCS (skycell->fpa, chip, hdu->header, pixelScale);
+	// we use a default FPA pixel scale of 1.0
+	pmAstromReadWCS (skycell->fpa, chip, hdu->header, 1.0);
     }
 
Index: /trunk/pswarp/src/pswarpMapGrid.c
===================================================================
--- /trunk/pswarp/src/pswarpMapGrid.c	(revision 12522)
+++ /trunk/pswarp/src/pswarpMapGrid.c	(revision 12523)
@@ -1,6 +1,8 @@
 # include "pswarp.h"
 
-// construct a grid with superpixel spacing of nXpix, nYpix
-// XXX for the moment, ignore readout->cell->chip offsets
+// pswarpMapGridFromImage builds a set (a grid) of locally-linear maps which convert the source
+// coordinates (src) to destination coordinates (dest).  we construct a grid with superpixel
+// spacing of nXpix, nYpix.  The transformation for each grid cell is valid for the superpixel.
+// The grid over-fills the source image so ever source image pixel is guaranteed to have a map.
 pswarpMapGrid *pswarpMapGridFromImage (pmReadout *dest, pmReadout *src, int nXpix, int nYpix) {
 
@@ -8,18 +10,22 @@
     int j, nj;
 
-    // XXX I was trying to match the grids too closely
-    // split the difference of the remainder
-    // int xMin = 0.5*(src->image->numCols % nXpix);
-    // int yMin = 0.5*(src->image->numRows % nYpix);
-    int xMin = 0;
-    int yMin = 0;
-
-    int nXpts = src->image->numCols / nXpix + 1;
-    if (src->image->numCols % nXpix) nXpts ++;
-    int nYpts = src->image->numRows / nYpix + 1;
-    if (src->image->numRows % nYpix) nYpts ++;
-
+    // start counting from the center of the superpixels
+    int xMin = 0.5*nXpix;
+    int yMin = 0.5*nYpix;
+
+    // the map is defined for coordinates in the image parent frame.
+    int Nx = src->image->numCols + src->image->col0;
+    int Ny = src->image->numRows + src->image->row0;
+
+    // allocate an extra superpixel to carry the remainder
+    int nXpts = Nx / nXpix;
+    int nYpts = Ny / nYpix;
+    if (Nx % nXpix) nXpts ++;
+    if (Ny % nYpix) nYpts ++;
+
+    // create the grid of maps
     pswarpMapGrid *grid = pswarpMapGridAlloc (nXpts, nYpts);
 
+    // measure the map for the center of each superpixel
     for (ni = 0, i = xMin; ni < nXpts; i += nXpix, ni++) {
 	for (nj = 0, j = yMin; nj < nYpts; j += nYpix, nj++) {
@@ -35,24 +41,30 @@
 }
 
+// set the grid coordinate (gridX,gridY) for the given source image coordinate (ix,iy)
 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;
+    *gridX = 0.5 + (ix - grid->xMin) / (double) grid->nXpix;
+    *gridY = 0.5 + (iy - grid->yMin) / (double) grid->nYpix;
     return true;
 }
 
-bool pswarpMapGridNextGrid_X (pswarpMapGrid *grid, int gridX, int *nextX) {
-
-    *nextX = gridX*grid->nXpix + grid->xMin + 0.5*grid->nXpix;
-    return true;
-}
-
-bool pswarpMapGridNextGrid_Y (pswarpMapGrid *grid, int gridY, int *nextY) {
-
-    *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
+// given the specified grid coordinate (gridX), return the x-coordinate for the source image
+// corresponding to the next grid cell
+int pswarpMapGridNextGrid_X (pswarpMapGrid *grid, int gridX) {
+
+    int nextX = (gridX + 0.5)*grid->nXpix + grid->xMin;
+    return nextX;
+}
+
+// given the specified grid coordinate (gridY), return the y-coordinate for the source image
+// corresponding to the next grid cell
+int pswarpMapGridNextGrid_Y (pswarpMapGrid *grid, int gridY) {
+
+    int nextY = (gridY + 0.5)*grid->nYpix + grid->yMin;
+    return nextY;
+}
+
+// measure the max error accumulated in applying one grid point to its neighbors
+// XXX double-check this
 double pswarpMapGridMaxError (pswarpMapGrid *grid) {
 
@@ -76,4 +88,5 @@
 }
 
+// given the source coordinate (inX,inY), return the destination coordinate (outX,outY)
 bool pswarpMapApply (double *outX, double *outY, pswarpMap *map, double inX, double inY) {
 
@@ -84,5 +97,7 @@
 }
 
-// determine the map for the given pixel from src to dest. pixel is in src coords
+// determine the (linear) map for the given pixel (ix,iy) from source image (src) to the destination image (dest)
+// pixel is in src coords. input and output pixel coordinates are in the parent frame of the image (Note that the
+// astrometric transformations are supplied for the parent image coordinate frame.
 bool pswarpMapSetLocalModel (pswarpMap *map, pmReadout *dest, pmReadout *src, int ix, int iy) {
 
Index: /trunk/pswarp/src/pswarpMatchRange.c
===================================================================
--- /trunk/pswarp/src/pswarpMatchRange.c	(revision 12522)
+++ /trunk/pswarp/src/pswarpMatchRange.c	(revision 12523)
@@ -2,17 +2,20 @@
 
 # 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); \
+    psPlaneTransformApply(destFP, chipDest->toFPA, destPix); \
+    psPlaneTransformApply (destTP, fpaDest->toTPA, destFP); \
+    psDeproject (sky, destTP, fpaDest->toSky); \
+    psProject (srcTP, sky, fpaSrc->toSky); \
+    psPlaneTransformApply (srcFP, fpaSrc->fromTPA, srcTP); \
+    psPlaneTransformApply (srcPix, chipSrc->fromFPA, srcFP); \
+    fprintf (stderr, "%6.1f,%6.1f -> ", srcPix->x, srcPix->y); \
     *minX = PS_MIN (*minX, srcPix->x); \
     *minY = PS_MIN (*minY, srcPix->y); \
     *maxX = PS_MAX (*maxX, srcPix->x); \
-    *maxY = PS_MAX (*maxY, srcPix->y);
+    *maxY = PS_MAX (*maxY, srcPix->y); \
+    fprintf (stderr, "%4d,%4d - %4d,%4d\n", *minX, *minY, *maxX, *maxY);
 
 // we are warping from src to dest.  find the max overlapping pixels in the INPUT (src)
-// coordinate frame
+// coordinate frame.  NOTE: these are in the parent pixel frame since the astrometric
+// transformation refers to the parent frame
 bool pswarpMatchRange (int *minX, int *minY, int *maxX, int *maxY, pmReadout *dest, pmReadout *src) {
 
@@ -30,58 +33,87 @@
     pmFPA *fpaDest = chipDest->parent;
 
-    *minX = src->image->numCols;
-    *minY = src->image->numRows;
-    *maxX = 0;
-    *maxY = 0;
+    *minX = src->image->numCols + src->image->col0 - 1;
+    *minY = src->image->numRows + src->image->row0 - 1;
+    *maxX = src->image->col0;
+    *maxY = src->image->row0;
+    fprintf (stderr, "%4d,%4d - %4d,%4d\n", *minX, *minY, *maxX, *maxY);
 
     // XXX save these as static for speed?
     psPlane *srcPix = psPlaneAlloc();
+    psPlane *srcFP  = psPlaneAlloc();
+    psPlane *srcTP  = psPlaneAlloc();
+
     psPlane *destPix = psPlaneAlloc();
+    psPlane *destFP  = psPlaneAlloc();
+    psPlane *destTP  = psPlaneAlloc();
 
-    psPlane *FP = psPlaneAlloc();
-    psPlane *TP = psPlaneAlloc();
     psSphere *sky = psSphereAlloc();
 
-    destPix->x = 0;
-    destPix->y = 0;
+    destPix->x = dest->image->col0;
+    destPix->y = dest->image->row0;
     TEST_MINMAX;
     
-    destPix->x = dest->image->numCols;
-    destPix->y = 0;
+    destPix->x = dest->image->col0 + dest->image->numCols;
+    destPix->y = dest->image->row0;
     TEST_MINMAX;
     
-    destPix->x = dest->image->numCols;
-    destPix->y = dest->image->numRows;
+    destPix->x = dest->image->col0 + dest->image->numCols;
+    destPix->y = dest->image->row0 + dest->image->numRows;
     TEST_MINMAX;
     
-    destPix->x = 0;
-    destPix->y = dest->image->numRows;
+    destPix->x = dest->image->col0;
+    destPix->y = dest->image->row0 + dest->image->numRows;
     TEST_MINMAX;
     
-    destPix->x = 0.5*dest->image->numCols;
-    destPix->y = 0.0;
+    destPix->x = dest->image->col0 + 0.5*dest->image->numCols;
+    destPix->y = dest->image->row0;
     TEST_MINMAX;
     
-    destPix->x = 0.0;
-    destPix->y = 0.5*dest->image->numRows;
+    destPix->x = dest->image->col0;
+    destPix->y = dest->image->row0 + 0.5*dest->image->numRows;
     TEST_MINMAX;
     
-    destPix->x = 0.5*dest->image->numCols;
-    destPix->y = dest->image->numRows;
+    destPix->x = dest->image->col0 + 0.5*dest->image->numCols;
+    destPix->y = dest->image->row0;
     TEST_MINMAX;
 
-    destPix->x = dest->image->numCols;
-    destPix->y = 0.5*dest->image->numRows;
+    destPix->x = dest->image->col0;
+    destPix->y = dest->image->row0 + 0.5*dest->image->numRows;
     TEST_MINMAX;
 
-    *minX = PS_MAX (*minX, 0);
-    *minY = PS_MAX (*minY, 0);
-    *maxX = PS_MIN (*maxX, src->image->numCols);
-    *maxY = PS_MIN (*maxY, src->image->numRows);
+    *minX = PS_MAX (*minX, src->image->col0);
+    *minY = PS_MAX (*minY, src->image->row0);
+    *maxX = PS_MIN (*maxX, src->image->numCols + src->image->col0);
+    *maxY = PS_MIN (*maxY, src->image->numRows + src->image->row0);
+
+    // demo forward and backward transformation
+    srcPix->x = *minX;
+    srcPix->y = *minY;
+    psPlaneTransformApply(srcFP, chipSrc->toFPA, srcPix); 
+    psPlaneTransformApply (srcTP, fpaSrc->toTPA, srcFP); 
+    psDeproject (sky, srcTP, fpaSrc->toSky); 
+    psProject (destTP, sky, fpaDest->toSky); 
+    psPlaneTransformApply (destFP, fpaDest->fromTPA, destTP); 
+    psPlaneTransformApply (destPix, chipDest->fromFPA, destFP); 
+
+    fprintf (stderr, "%f,%f -> %f,%f ", srcPix->x, srcPix->y, destPix->x, destPix->y);
+
+    psPlaneTransformApply(destFP, chipDest->toFPA, destPix); 
+    psPlaneTransformApply (destTP, fpaDest->toTPA, destFP); 
+    psDeproject (sky, destTP, fpaDest->toSky); 
+    psProject (srcTP, sky, fpaSrc->toSky); 
+    psPlaneTransformApply (srcFP, fpaSrc->fromTPA, srcTP); 
+    psPlaneTransformApply (srcPix, chipSrc->fromFPA, srcFP); 
+
+    fprintf (stderr, "-> %f,%f ", srcPix->x, srcPix->y);
 
     psFree (srcPix);
+    psFree (srcFP);
+    psFree (srcTP);
+
     psFree (destPix);
-    psFree (FP);
-    psFree (TP);
+    psFree (destFP);
+    psFree (destTP);
+
     psFree (sky);
 
Index: /trunk/pswarp/src/pswarpTransformReadout_Opt.c
===================================================================
--- /trunk/pswarp/src/pswarpTransformReadout_Opt.c	(revision 12522)
+++ /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
-
-