Index: trunk/psLib/src/imageops/psImageGeomManip.c
===================================================================
--- trunk/psLib/src/imageops/psImageGeomManip.c	(revision 31435)
+++ trunk/psLib/src/imageops/psImageGeomManip.c	(revision 31446)
@@ -702,5 +702,5 @@
     } \
     /* note: output(i,j) = input(i+0.5-dx,j+0.5-dy) */ \
-    /* positive dx,dy moves pixel i-dx,j-dy to i,y */ \
+    /* positive dx,dy moves pixel i-dx,j-dy to i,j */ \
     /* also: pixel center is 0.5,0.5 */ \
     for (int row = 0; row < numRows; row++) { \
Index: trunk/psLib/src/imageops/psImageInterpolate.c
===================================================================
--- trunk/psLib/src/imageops/psImageInterpolate.c	(revision 31435)
+++ trunk/psLib/src/imageops/psImageInterpolate.c	(revision 31446)
@@ -31,4 +31,6 @@
 #include "psImageConvolve.h"
 
+# define IS_BILIN_SEPARABLE 0
+
 #include "psImageInterpolate.h"
 
@@ -48,4 +50,5 @@
 };
 
+# if (IS_BILIN_SEPARABLE)
 /// Generate a linear interpolation kernel
 static inline void interpolationKernelBilinear(psF32 *kernel, // Kernel vector to populate
@@ -56,4 +59,19 @@
     kernel[1] = frac;
 }
+# else
+/// Generate a linear interpolation kernel
+static inline void interpolationKernelBilinear(
+    psF32 kernel[kernelSizes[PS_INTERPOLATE_BILINEAR]][kernelSizes[PS_INTERPOLATE_BILINEAR]], // Kernel
+    float xFrac, float yFrac // Fraction of pixel
+    )
+{
+    // (xF*yF, (1-xF)(1-yF) = 1 - xF - yF + xFyF 
+
+    kernel[0][0] = 1.0 - xFrac - yFrac + xFrac*yFrac;
+    kernel[1][0] = yFrac - xFrac*yFrac;
+    kernel[0][1] = xFrac - xFrac*yFrac;
+    kernel[1][1] = xFrac*yFrac;
+}
+# endif
 
 #if 0
@@ -173,7 +191,9 @@
         // Nothing to pre-compute
         break;
+# if (IS_BILIN_SEPARABLE)
       case PS_INTERPOLATE_BILINEAR:
         interpolationKernelBilinear(kernel, frac);
         break;
+# endif
       case PS_INTERPOLATE_GAUSS:
         interpolationKernelGauss(kernel, kernelSizes[mode], 0.5, frac);
@@ -184,4 +204,7 @@
         interpolationKernelLanczos(kernel, kernelSizes[mode], frac);
         break;
+# if (!IS_BILIN_SEPARABLE)
+      case PS_INTERPOLATE_BILINEAR:
+# endif
       case PS_INTERPOLATE_BIQUADRATIC:  // 2D kernel
       default:
@@ -234,5 +257,7 @@
 
     switch (mode) {
+# if (IS_BILIN_SEPARABLE)
       case PS_INTERPOLATE_BILINEAR:
+# endif
       case PS_INTERPOLATE_GAUSS:
       case PS_INTERPOLATE_LANCZOS2:
@@ -258,4 +283,7 @@
         }
         break;
+# if (!IS_BILIN_SEPARABLE)
+      case PS_INTERPOLATE_BILINEAR:
+# endif
       case PS_INTERPOLATE_BIQUADRATIC:
         // 2D kernel, would cost too much memory to pre-calculate
@@ -360,4 +388,11 @@
     int xCentral = floor((X) - 0.5), yCentral = floor((Y) - 0.5); /* Central pixel */ \
     float xFrac = (X) - xCentral - 0.5, yFrac = (Y) - yCentral - 0.5; /* Fraction of pixel */ \
+    bool xExact = fabsf(xFrac) < FLT_EPSILON, yExact = fabsf(yFrac) < FLT_EPSILON; /* Are shifts exact? */ \
+
+// Set up the kernel parameters; defines some useful values
+// EAM : I'm unsure of the right answer here
+#define INTERPOLATE_SETUP_EAM(X, Y)				 \
+    int xCentral = floor((X)), yCentral = floor((Y)); /* Central pixel */ \
+    float xFrac = (X) - xCentral, yFrac = (Y) - yCentral; /* Fraction of pixel */ \
     bool xExact = fabsf(xFrac) < FLT_EPSILON, yExact = fabsf(yFrac) < FLT_EPSILON; /* Are shifts exact? */ \
 
@@ -722,5 +757,5 @@
 }
 
-// Interpolation engine for (separable) interpolation kernels
+// Interpolation engine for (non-separable) interpolation kernels
 static psImageInterpolateStatus interpolateKernel(double *imageValue, double *varianceValue,
                                                   psImageMaskType *maskValue, float x, float y,
@@ -747,6 +782,25 @@
     // Get the appropriate kernels
     psF32 kernel[size][size];
+
+# if (IS_BILIN_SEPARABLE)
     psAssert(mode == PS_INTERPOLATE_BIQUADRATIC, "Mode is %x", mode);
     interpolationKernelBiquadratic(kernel, xFrac, yFrac);
+# else
+    psAssert((mode == PS_INTERPOLATE_BIQUADRATIC) || (mode == PS_INTERPOLATE_BILINEAR), "Mode is %x", mode);
+    if (mode == PS_INTERPOLATE_BIQUADRATIC) {
+	interpolationKernelBiquadratic(kernel, xFrac, yFrac);
+    } else {
+	interpolationKernelBilinear(kernel, xFrac, yFrac);
+    }
+# endif
+
+# if (0)
+    for (int i = 0; i < size; i++) {
+	for (int j = 0; j < size; j++) {
+	    fprintf (stderr, "%f ", kernel[i][j]);
+	}
+	fprintf (stderr, "\n");
+    }
+# endif
 
     float sumKernel = 0.0;              // Sum of the kernel
@@ -908,4 +962,7 @@
         return interpolateKernel(imageValue, varianceValue, maskValue, x, y, interp);
       case PS_INTERPOLATE_BILINEAR:
+# if (!IS_BILIN_SEPARABLE)
+        return interpolateKernel(imageValue, varianceValue, maskValue, x, y, interp);
+# endif
       case PS_INTERPOLATE_GAUSS:
       case PS_INTERPOLATE_LANCZOS2:
Index: trunk/psLib/src/imageops/psImageMaskOps.c
===================================================================
--- trunk/psLib/src/imageops/psImageMaskOps.c	(revision 31435)
+++ trunk/psLib/src/imageops/psImageMaskOps.c	(revision 31446)
@@ -146,6 +146,6 @@
     for (int iy = 0; iy < image->numRows; iy++) { \
         for (int ix = 0; ix < image->numCols; ix++) { \
-            dx = ix + image->col0 - x; \
-            dy = iy + image->row0 - y; \
+            dx = ix + 0.5 + image->col0 - x; \
+            dy = iy + 0.5 + image->row0 - y; \
             r2 = PS_SQR(dx) + PS_SQR(dy); \
             if (r2 <= R2) { \
@@ -197,6 +197,6 @@
     for (int iy = 0; iy < image->numRows; iy++) { \
         for (int ix = 0; ix < image->numCols; ix++) { \
-            dx = ix + image->col0 - x; \
-            dy = iy + image->row0 - y; \
+            dx = ix + 0.5 + image->col0 - x; \
+            dy = iy + 0.5 + image->row0 - y; \
             r2 = PS_SQR(dx) + PS_SQR(dy); \
             if (r2 > R2) { \