Index: trunk/psLib/src/imageops/psImageStats.c
===================================================================
--- trunk/psLib/src/imageops/psImageStats.c	(revision 1234)
+++ trunk/psLib/src/imageops/psImageStats.c	(revision 1251)
@@ -13,6 +13,5 @@
 #include "psImage.h"
 #include "psFunctions.h"
-#include "float.h"
-#include <math.h>
+#include "psImageStats.h"
 
 /// This routine must determine the various statistics for the image.
@@ -377,64 +376,138 @@
 linear interpolation is performed on the image.
  *****************************************************************************/
-float psImagePixelInterpolation(psImage *input,
-                                float x,
-                                float y)
-{
-    float floorX = 0.0;
-    float floorY = 0.0;
-    float fracX = 0.0;
-    float fracY = 0.0;
-    int intFloorX = 0;
-    int intFloorY = 0;
-    float x1 = 0.0;
-    float x2 = 0.0;
-    float pixel = 0.0;
-
-    if ((x < 0.0) ||
-            (x > ((float) input->numRows-1)) ||
-            (y < 0.0) ||
-            (y > ((float) input->numCols-1))) {
-        psAbort(__func__,
-                "Fractional coordinates (%f %f) outside image range.", x, y);
-    }
-
-    floorX = floorf(x);
-    intFloorX = (int) floorX;
-    fracX = x - floorX;
-
-    floorY = floorf(y);
-    intFloorY = (int) floorY;
-    fracY = y - floorY;
-
-    if (intFloorX == (input->numRows-1)) {
-        pixel = input->data.F32[intFloorX][intFloorY];
-        if (intFloorY < (input->numCols-1)) {
-            pixel+= fracY * (input->data.F32[intFloorX][intFloorY+1] -
-                             input->data.F32[intFloorX][intFloorY]);
-        }
-        return(pixel);
-    } else if (intFloorX == (input->numCols-1)) {
-        if (intFloorX < (input->numRows-1)) {
-            pixel+= fracX * (input->data.F32[intFloorX+1][intFloorY] -
-                             input->data.F32[intFloorX][intFloorY]);
-        }
-        return(pixel);
-    }
-
-
-    x1 = input->data.F32[intFloorX][intFloorY];
-    x1+= fracY * (input->data.F32[intFloorX][intFloorY+1] -
-                  input->data.F32[intFloorX][intFloorY]);
-
-    x2 = input->data.F32[intFloorX+1][intFloorY];
-    x2+= fracY * (input->data.F32[intFloorX+1][intFloorY+1] -
-                  input->data.F32[intFloorX+1][intFloorY]);
-    pixel = x1;
-    pixel+= fracX * (x2 - x1);
-
-    return(pixel);
-}
-
-
-
-
+psF32 psImagePixelInterpolate(
+    const psImage *input,
+    float x,
+    float y,
+    psF32 unexposedValue,
+    psImageInterpolateMode mode)
+{
+
+    if (input == NULL) {
+        psError(__func__,"Image can not be NULL.");
+        return unexposedValue;
+    }
+
+    #define PSIMAGE_PIXEL_INTERPOLATE_CASE(TYPE) \
+case PS_TYPE_##TYPE: \
+    switch (mode) { \
+    case PS_INTERPOLATE_FLAT: \
+        return p_psImagePixelInterpolateFLAT_##TYPE(input,x,y,unexposedValue); \
+        break; \
+    case PS_INTERPOLATE_BILINEAR: \
+        return p_psImagePixelInterpolateBILINEAR_##TYPE(input,x,y,unexposedValue); \
+        break; \
+    default: \
+        psError(__func__,"Unsupported interpolation mode (#%d)",mode); \
+        return unexposedValue; \
+    } \
+    break
+
+    switch (input->type.type) {
+        PSIMAGE_PIXEL_INTERPOLATE_CASE(F32);
+    default:
+        psError(__func__,"Unsupported image datatype (%d)",input->type.type);
+    }
+
+    return unexposedValue;
+}
+
+#define PSIMAGE_PIXEL_INTERPOLATE_FLAT(TYPE) \
+inline psF64 p_psImagePixelInterpolateFLAT_##TYPE(const psImage *input, \
+        float x, \
+        float y, \
+        psF64 unexposedValue) \
+{ \
+    int intX = (int) round((psF64)(x) - 0.5); \
+    int intY = (int) round((psF64)(y) - 0.5); \
+    int lastX = input->numCols - 1; \
+    int lastY = input->numRows - 1; \
+    \
+    if ((intX < 0) || \
+            (intX > lastX) || \
+            (intY < 0) || \
+            (intY > lastY)) { \
+        return unexposedValue; \
+    } \
+    \
+    return input->data.TYPE[intY][intX]; \
+}
+
+PSIMAGE_PIXEL_INTERPOLATE_FLAT(U8)
+PSIMAGE_PIXEL_INTERPOLATE_FLAT(U16)
+PSIMAGE_PIXEL_INTERPOLATE_FLAT(U32)
+PSIMAGE_PIXEL_INTERPOLATE_FLAT(U64)
+PSIMAGE_PIXEL_INTERPOLATE_FLAT(S8)
+PSIMAGE_PIXEL_INTERPOLATE_FLAT(S16)
+PSIMAGE_PIXEL_INTERPOLATE_FLAT(S32)
+PSIMAGE_PIXEL_INTERPOLATE_FLAT(S64)
+PSIMAGE_PIXEL_INTERPOLATE_FLAT(F32)
+PSIMAGE_PIXEL_INTERPOLATE_FLAT(F64)
+
+#define PSIMAGE_PIXEL_INTERPOLATE_BILINEAR(TYPE) \
+inline psF64 p_psImagePixelInterpolateBILINEAR_##TYPE(const psImage *input, \
+        float x, \
+        float y, \
+        psF64 unexposedValue) \
+{ \
+    double floorX = floor((psF64)(x) - 0.5); \
+    double floorY = floor((psF64)(y) - 0.5); \
+    double fracX = x - 0.5 - floorX; \
+    double fracY = y - 0.5 - floorY; \
+    int intFloorX = (int) floorX; \
+    int intFloorY = (int) floorY; \
+    int lastX = input->numCols - 1; \
+    int lastY = input->numRows - 1; \
+    double rx = 0.0; \
+    psF64 pixel = 0.0; \
+    ps##TYPE* currentRow; \
+    ps##TYPE* nextRow; \
+    \
+    if ((intFloorX < 0) || \
+            (intFloorX > lastX) || \
+            (intFloorY < 0) || \
+            (intFloorY > lastY)) { \
+        return unexposedValue; \
+    } \
+    \
+    currentRow = input->data.TYPE[intFloorY]; \
+    if (intFloorY == lastY) { \
+        pixel = currentRow[intFloorX]; \
+        if (intFloorX < lastX) { \
+            pixel+= fracY * ((psF64)currentRow[intFloorX+1] - \
+                             (psF64)currentRow[intFloorX]); \
+        } \
+        return(pixel); \
+    } \
+    nextRow = input->data.TYPE[intFloorY+1]; \
+    if (intFloorX == lastX) { \
+        pixel = currentRow[intFloorX]; \
+        if (intFloorY < lastY) { \
+            pixel+= fracX * ((psF64)nextRow[intFloorX] - \
+                             (psF64)currentRow[intFloorX]); \
+        } \
+        return(pixel); \
+    } \
+    \
+    rx = currentRow[intFloorX] + \
+         fracX * ((psF64)currentRow[intFloorX+1] - \
+                  (psF64)currentRow[intFloorX]); \
+    \
+    pixel = rx + fracY * ((psF64)nextRow[intFloorX] + \
+                          fracX * ((psF64)nextRow[intFloorX+1] - \
+                                   (psF64)nextRow[intFloorX]) - rx); \
+    \
+    return(pixel); \
+}
+
+PSIMAGE_PIXEL_INTERPOLATE_BILINEAR(U8)
+PSIMAGE_PIXEL_INTERPOLATE_BILINEAR(U16)
+PSIMAGE_PIXEL_INTERPOLATE_BILINEAR(U32)
+PSIMAGE_PIXEL_INTERPOLATE_BILINEAR(U64)
+PSIMAGE_PIXEL_INTERPOLATE_BILINEAR(S8)
+PSIMAGE_PIXEL_INTERPOLATE_BILINEAR(S16)
+PSIMAGE_PIXEL_INTERPOLATE_BILINEAR(S32)
+PSIMAGE_PIXEL_INTERPOLATE_BILINEAR(S64)
+PSIMAGE_PIXEL_INTERPOLATE_BILINEAR(F32)
+PSIMAGE_PIXEL_INTERPOLATE_BILINEAR(F64)
+