Index: trunk/psLib/src/imageops/psImageInterpolate.c
===================================================================
--- trunk/psLib/src/imageops/psImageInterpolate.c	(revision 34724)
+++ trunk/psLib/src/imageops/psImageInterpolate.c	(revision 34800)
@@ -31,5 +31,5 @@
 #include "psImageConvolve.h"
 
-# define IS_BILIN_SEPARABLE 0
+# define IS_BILIN_SEPARABLE 1
 
 #include "psImageInterpolate.h"
@@ -189,4 +189,5 @@
     switch (mode) {
       case PS_INTERPOLATE_FLAT:
+      case PS_INTERPOLATE_BILINEAR_SIMPLE:
         // Nothing to pre-compute
         break;
@@ -287,6 +288,8 @@
 # endif
       case PS_INTERPOLATE_BIQUADRATIC:
+      case PS_INTERPOLATE_BILINEAR_SIMPLE:
         // 2D kernel, would cost too much memory to pre-calculate
         break;
+
       default:
         psAbort("Unsupported interpolation mode: %x", mode);
@@ -937,4 +940,59 @@
 
 
+psImageInterpolateStatus interpolateJustWork(double *imageValue, double *varianceValue,
+					     psImageMaskType *maskValue, float x, float y,
+					     const psImageInterpolation *interp) {
+  float u1,u2;
+  int xl,xh;
+  int yl,yh;
+  const psImage *image = interp->image; // Image to interpolate
+
+  xl = floor(x);
+  if (xl < 0) { xl = 0; }
+  if (xl >= image->numCols) {xl = image->numCols - 1; }
+  xh = xl + 1;
+  if (xh >= image->numCols) {xh = image->numCols - 1; }
+
+  yl = floor(y);
+  if (yl < 0) { yl = 0; }
+  if (yl >= image->numRows) {yl = image->numRows - 1; }
+  yh = yl + 1;
+  if (yh >= image->numRows) {yh = image->numRows - 1; }
+
+  if (imageValue && image) {
+    if (image->data.F32[yl][xl] == image->data.F32[yl][xh]) {
+      u1 = image->data.F32[yl][xh];
+    }
+    else {
+      u1 = (xh - x) * image->data.F32[yl][xl] + (x - xl) * image->data.F32[yl][xh];
+    }
+    if (image->data.F32[yh][xl] == image->data.F32[yh][xh]) {
+      u2 = image->data.F32[yh][xh];
+    }
+    else {
+      u2 = (xh - x) * image->data.F32[yh][xl] + (x - xl) * image->data.F32[yh][xh];
+    }
+    if (u1 == u2) {
+      *imageValue = u1;
+    }
+    else {
+      *imageValue = (yh - y) * u1 + (y - yl) * u2;
+    }
+    if ((floor(x) >= image->numCols)||
+	(floor(y) >= image->numRows)) {
+      *imageValue = NAN;
+    }
+  }
+#if (0)
+  if (*imageValue == 0.0) {
+    fprintf(stderr,"IJK: Zero!: %g %g [%d %d %d %d] %g %g (%g %g %g %g)\n",
+	    x,y,xl,xh,yl,yh,u1,u2,
+	    image->data.F32[yl][xl],image->data.F32[yl][xh],image->data.F32[yh][xl],image->data.F32[yh][xh]
+	    );
+  }
+#endif
+  return PS_INTERPOLATE_STATUS_GOOD;
+}
+  
 
 psImageInterpolateStatus psImageInterpolate(double *imageValue, double *varianceValue,
@@ -969,7 +1027,9 @@
       case PS_INTERPOLATE_BIQUADRATIC:
         return interpolateKernel(imageValue, varianceValue, maskValue, x, y, interp);
+      case PS_INTERPOLATE_BILINEAR_SIMPLE:
+	return interpolateJustWork(imageValue, varianceValue, maskValue, x, y, interp);
       case PS_INTERPOLATE_BILINEAR:
 # if (!IS_BILIN_SEPARABLE)
-        return interpolateKernel(imageValue, varianceValue, maskValue, x, y, interp);
+	return interpolateKernel(imageValue, varianceValue, maskValue, x, y, interp);
 # endif
       case PS_INTERPOLATE_GAUSS:
@@ -1019,4 +1079,5 @@
     switch (mode) {
       case PS_INTERPOLATE_FLAT:
+      case PS_INTERPOLATE_BILINEAR_SIMPLE:
         // No smearing by design
         return 1.0;
@@ -1083,4 +1144,5 @@
     if (!strcasecmp(name, "FLAT"))     return PS_INTERPOLATE_FLAT;
     if (!strcasecmp(name, "BILINEAR")) return PS_INTERPOLATE_BILINEAR;
+    if (!strcasecmp(name, "SIMPLEBILINEAR")) return PS_INTERPOLATE_BILINEAR_SIMPLE;
     if (!strcasecmp(name, "BIQUADRATIC")) return PS_INTERPOLATE_BIQUADRATIC;
     if (!strcasecmp(name, "GAUSS"))    return PS_INTERPOLATE_GAUSS;
Index: trunk/psLib/src/imageops/psImageInterpolate.h
===================================================================
--- trunk/psLib/src/imageops/psImageInterpolate.h	(revision 34724)
+++ trunk/psLib/src/imageops/psImageInterpolate.h	(revision 34800)
@@ -33,4 +33,5 @@
     PS_INTERPOLATE_LANCZOS3,            ///< Sinc interpolation with 6x6 pixel kernel
     PS_INTERPOLATE_LANCZOS4,            ///< Sinc interpolation with 8x8 pixel kernel
+    PS_INTERPOLATE_BILINEAR_SIMPLE,     ///< Simple manual bilinear interpolation
 } psImageInterpolateMode;
 
Index: trunk/psLib/src/imageops/psImagePixelManip.c
===================================================================
--- trunk/psLib/src/imageops/psImagePixelManip.c	(revision 34724)
+++ trunk/psLib/src/imageops/psImagePixelManip.c	(revision 34800)
@@ -216,15 +216,54 @@
     }
 
-
+#define psImageOverlayLoopClean(DATATYPE,OP) {	 \
+      for (int row=y0;row<imageRowLimit;row++) {	    \
+	ps##DATATYPE* imageRow = image->data.DATATYPE[row];	   \
+	ps##DATATYPE* overlayRow = overlay->data.DATATYPE[row-y0]; \
+	for (int col=x0;col<imageColLimit;col++) {		   \
+	  if (!isfinite(imageRow[col])) {			   \
+	    imageRow[col] OP overlayRow[col-x0];		   \
+	  }							   \
+	  else if (!isfinite(overlayRow[col-x0])) {		   \
+	    imageRow[col] = imageRow[col];			   \
+	  }							   \
+	  else {						   \
+	    imageRow[col] = overlayRow[col-x0];			   \
+	  }								\
+	}								\
+      }									\
+      pixelsOverlaid += (imageRowLimit - y0) * (imageColLimit - x0);	\
+    }
+
+#define psImageOverlayLoopMask(DATATYPE) {		    \
+      for (int row=y0;row<imageRowLimit;row++) {	    \
+	ps##DATATYPE* imageRow = image->data.DATATYPE[row];	   \
+	ps##DATATYPE* overlayRow = overlay->data.DATATYPE[row-y0]; \
+	for (int col=x0;col<imageColLimit;col++) {		   \
+	  if (!isfinite(imageRow[col])) {			   \
+	    imageRow[col] = overlayRow[col-x0];		   \
+	  }							   \
+	  else if (!isfinite(overlayRow[col-x0])) {		   \
+	    imageRow[col] = imageRow[col];			   \
+	  }							   \
+	  else {						   \
+	    imageRow[col] &= overlayRow[col-x0];		   \
+	  }								\
+	}								\
+      }									\
+      pixelsOverlaid += (imageRowLimit - y0) * (imageColLimit - x0);	\
+    }
+	
+    
+    
     #define psImageOverlayLoop(DATATYPE,OP) { \
         for (int row=y0;row<imageRowLimit;row++) { \
             ps##DATATYPE* imageRow = image->data.DATATYPE[row]; \
             ps##DATATYPE* overlayRow = overlay->data.DATATYPE[row-y0]; \
-            for (int col=x0;col<imageColLimit;col++) { \
-                imageRow[col] OP overlayRow[col-x0]; \
-            } \
-        } \
+            for (int col=x0;col<imageColLimit;col++) {		       \
+                imageRow[col] OP overlayRow[col-x0];		       \
+	      }							       \
+        }							       \
         pixelsOverlaid += (imageRowLimit - y0) * (imageColLimit - x0); \
-    }
+      }
 
     #define psImageOverlayLoopDivide(DATATYPE,BADVALUE) { \
@@ -273,9 +312,12 @@
         psImageOverlayLoop(DATATYPE,*=); \
         break; \
+    case 'E': \
+      psImageOverlayLoopClean(DATATYPE,=); \
+      break;				   \
     case '/': \
         psImageOverlayLoopDivide(DATATYPE,BADVALUE); \
         break; \
     case '=': \
-        psImageOverlaySetLoop(DATATYPE); \
+        psImageOverlaySetLoop(DATATYPE);		\
         break; \
     default: \
@@ -286,14 +328,45 @@
     } \
     break;
+    #define psImageOverlayCaseInt(DATATYPE,BADVALUE) \
+case PS_TYPE_##DATATYPE: \
+    switch (*op) { \
+    case '+': \
+        psImageOverlayLoop(DATATYPE,+=); \
+        break; \
+    case '-': \
+        psImageOverlayLoop(DATATYPE,-=); \
+        break; \
+    case '*': \
+        psImageOverlayLoop(DATATYPE,*=); \
+        break; \
+    case 'E': \
+      psImageOverlayLoopClean(DATATYPE,=); \
+      break;				   \
+    case 'M':				   \
+      psImageOverlayLoopMask(DATATYPE);	   \
+      break;				   \
+    case '/': \
+        psImageOverlayLoopDivide(DATATYPE,BADVALUE); \
+        break; \
+    case '=': \
+        psImageOverlaySetLoop(DATATYPE);		\
+        break; \
+    default: \
+        psError(PS_ERR_BAD_PARAMETER_VALUE, true, \
+                _("Specified operation, '%s', is not supported."), \
+                op); \
+        return pixelsOverlaid; \
+    } \
+    break;
 
     switch (type) {
-        psImageOverlayCase(U8, 0);
-        psImageOverlayCase(U16,0);
-        psImageOverlayCase(U32,0);       // Not a requirement
-        psImageOverlayCase(U64,0);       // Not a requirement
-        psImageOverlayCase(S8, 0);
-        psImageOverlayCase(S16,0);
-        psImageOverlayCase(S32,0);       // Not a requirement
-        psImageOverlayCase(S64,0);       // Not a requirement
+        psImageOverlayCaseInt(U8, 0);
+        psImageOverlayCaseInt(U16,0);
+        psImageOverlayCaseInt(U32,0);       // Not a requirement
+        psImageOverlayCaseInt(U64,0);       // Not a requirement
+        psImageOverlayCaseInt(S8, 0);
+        psImageOverlayCaseInt(S16,0);
+        psImageOverlayCaseInt(S32,0);       // Not a requirement
+        psImageOverlayCaseInt(S64,0);       // Not a requirement
         psImageOverlayCase(F32,NAN);
         psImageOverlayCase(F64,NAN);