Index: trunk/pswarp/src/pswarpTransformReadout_Opt.c
===================================================================
--- trunk/pswarp/src/pswarpTransformReadout_Opt.c	(revision 12528)
+++ trunk/pswarp/src/pswarpTransformReadout_Opt.c	(revision 12744)
@@ -1,4 +1,3 @@
 # include "pswarp.h"
-psImageInterpolateMode psImageInterpolateModeFromString (char *name);
 
 // NOTE: in this function, the coordinates are transformed from the OUTPUT to the INPUT
@@ -11,5 +10,5 @@
     psImageInterpolateMode interpolationMode = PS_INTERPOLATE_NONE;
 
-    // XXX this implementation currently ignores the use of the region 
+    // XXX this implementation currently ignores the use of the region
     psImage *region = NULL;
 
@@ -19,8 +18,8 @@
     psMetadata *recipe  = psMetadataLookupPtr (NULL, config->recipes, PSWARP_RECIPE);
     if (!recipe) {
-	psError(PSWARP_ERR_CONFIG, true, "Can't find PSWARP recipe!\n");
-	return false;
+        psError(PSWARP_ERR_CONFIG, true, "Can't find PSWARP recipe!\n");
+        return false;
     }
-    
+
     int nGridX = psMetadataLookupS32 (&status, recipe, "GRID.NX");
     if (!status) nGridX = 128;
@@ -30,15 +29,15 @@
     char *name = psMetadataLookupStr (&status, recipe, "INTERPOLATION.MODE");
     if (!name) {
-	interpolationMode = PS_INTERPOLATE_BILINEAR;
-	psLogMsg ("pswarp", 3, "defaulting to bilinear interpolation\n");
+        interpolationMode = PS_INTERPOLATE_BILINEAR;
+        psLogMsg ("pswarp", 3, "defaulting to bilinear interpolation\n");
     } else {
-	interpolationMode = psImageInterpolateModeFromString (name);
-	if (interpolationMode == PS_INTERPOLATE_NONE) {
-	    interpolationMode = PS_INTERPOLATE_BILINEAR;
-	    psLogMsg ("pswarp", 3, "unknown interpolation mode %s, defaulting to bilinear interpolation\n", name);
-	}
+        interpolationMode = psImageInterpolateModeFromString (name);
+        if (interpolationMode == PS_INTERPOLATE_NONE) {
+            interpolationMode = PS_INTERPOLATE_BILINEAR;
+            psLogMsg ("pswarp", 3, "unknown interpolation mode %s, defaulting to bilinear interpolation\n", name);
+        }
     }
 
-    psPlane *inPix = psPlaneAlloc();	// Coordinates on the input detector
+    psPlane *inPix = psPlaneAlloc();    // Coordinates on the input detector
     psImage *inImage = input->image;
     psF32 **outData = output->image->data.F32;
@@ -65,5 +64,5 @@
     psLogMsg ("pswarp", 3, "maximum error using this grid sampling: %f\n", maxError);
 
-    pswarpMapGridSetGrid (grid, minX, minY, &gridX, &gridY); 
+    pswarpMapGridSetGrid (grid, minX, minY, &gridX, &gridY);
     nextGridY = pswarpMapGridNextGrid_Y (grid, gridY);
     map = grid->maps[gridX][gridY];
@@ -77,41 +76,53 @@
     nextGridXo = pswarpMapGridNextGrid_X (grid, gridX);
 
+    psImageInterpolateOptions *interp = psImageInterpolateOptionsAlloc(PS_INTERPOLATE_BILINEAR, inImage,
+                                                                       NULL, NULL, 0, NAN, NAN, 0, 0, 0.0);
+
     // Iterate over the output image pixels (parent frame)
     for (int y = minY; y < maxY; y++) {
-	if (y >= nextGridY) {
-	    gridY ++;
-	    nextGridY = pswarpMapGridNextGrid_Y (grid, gridY);
-	    map = grid->maps[gridX][gridY];
-	}
+        if (y >= nextGridY) {
+            gridY ++;
+            nextGridY = pswarpMapGridNextGrid_Y (grid, gridY);
+            map = grid->maps[gridX][gridY];
+        }
 
-	gridX = gridXo;
-	nextGridX = nextGridXo;
-	map = grid->maps[gridX][gridY];
-	for (int x = minX; x < maxX; x++) {
-	    if (x >= nextGridX) {
-		gridX ++;
-		nextGridX = pswarpMapGridNextGrid_X (grid, gridX);
-		map = grid->maps[gridX][gridY];
-	    }
+        gridX = gridXo;
+        nextGridX = nextGridXo;
+        map = grid->maps[gridX][gridY];
+        for (int x = minX; x < maxX; x++) {
+            if (x >= nextGridX) {
+                gridX ++;
+                nextGridX = pswarpMapGridNextGrid_X (grid, gridX);
+                map = grid->maps[gridX][gridY];
+            }
 
-	    // Only transform those pixels requested
-	    if (region && region->data.U8[y][x]) continue;
+            // Only transform those pixels requested
+            if (region && region->data.U8[y][x]) continue;
 
-	    // 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);
+            // 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;
+            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 include mask
-	    // XXX apply scale and offset?
-	    // 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, interpolationMode);
-	}
+            // XXX include mask
+            // XXX apply scale and offset?
+            // psImagePixelInterpolate determines the value at pixel coordinate (x,y) in child coordinates
+            double value;
+            if (!psImageInterpolate(&value, NULL, NULL, inPix->x - inCol0, inPix->y - inRow0, interp)) {
+                psError(PS_ERR_UNKNOWN, false, "Unable to interpolate image.");
+                psFree(interp);
+                psFree (inPix);
+                psFree (grid);
+                return false;
+            }
+            outData[y-outRow0][x-outCol0] = value;
+        }
     }
 
+    psFree(interp);
     psFree (inPix);
     psFree (grid);