Index: trunk/psLib/src/image/psImage.c
===================================================================
--- trunk/psLib/src/image/psImage.c	(revision 1205)
+++ trunk/psLib/src/image/psImage.c	(revision 1261)
@@ -9,6 +9,6 @@
  *  @author Ross Harman, MHPCC
  *
- *  @version $Revision: 1.33 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2004-07-09 21:48:07 $
+ *  @version $Revision: 1.34 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2004-07-22 20:09:04 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -400,2 +400,236 @@
 }
 
+/*****************************************************************************
+p_psImagePixelInterpolation(image, x, y): this routine takes as input an
+image and coordinates (x, y) and produces as output the corresponding pixel
+value at the those coordinates.  For fractional corrdinates (x, y), 2-D
+linear interpolation is performed on the image.
+ *****************************************************************************/
+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); \
+    } \
+    break
+
+    switch (input->type.type) {
+        PSIMAGE_PIXEL_INTERPOLATE_CASE(U8);
+        PSIMAGE_PIXEL_INTERPOLATE_CASE(U16);
+        PSIMAGE_PIXEL_INTERPOLATE_CASE(U32);
+        PSIMAGE_PIXEL_INTERPOLATE_CASE(U64);
+        PSIMAGE_PIXEL_INTERPOLATE_CASE(S8);
+        PSIMAGE_PIXEL_INTERPOLATE_CASE(S16);
+        PSIMAGE_PIXEL_INTERPOLATE_CASE(S32);
+        PSIMAGE_PIXEL_INTERPOLATE_CASE(S64);
+        PSIMAGE_PIXEL_INTERPOLATE_CASE(F32);
+        PSIMAGE_PIXEL_INTERPOLATE_CASE(F64);
+        PSIMAGE_PIXEL_INTERPOLATE_CASE(C32);
+        PSIMAGE_PIXEL_INTERPOLATE_CASE(C64);
+    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]; \
+}
+
+#define PSIMAGE_PIXEL_INTERPOLATE_FLAT_COMPLEX(TYPE) \
+inline psC64 p_psImagePixelInterpolateFLAT_##TYPE(const psImage *input, \
+        float x, \
+        float y, \
+        psC64 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)
+PSIMAGE_PIXEL_INTERPOLATE_FLAT_COMPLEX(C32)
+PSIMAGE_PIXEL_INTERPOLATE_FLAT_COMPLEX(C64)
+
+#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); \
+}
+
+#define PSIMAGE_PIXEL_INTERPOLATE_BILINEAR_COMPLEX(TYPE) \
+inline psC64 p_psImagePixelInterpolateBILINEAR_##TYPE(const psImage *input, \
+        float x, \
+        float y, \
+        psC64 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; \
+    psC64 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 * (currentRow[intFloorX+1] - \
+                             currentRow[intFloorX]); \
+        } \
+        return(pixel); \
+    } \
+    nextRow = input->data.TYPE[intFloorY+1]; \
+    if (intFloorX == lastX) { \
+        pixel = currentRow[intFloorX]; \
+        if (intFloorY < lastY) { \
+            pixel+= fracX * (nextRow[intFloorX] - \
+                             currentRow[intFloorX]); \
+        } \
+        return(pixel); \
+    } \
+    \
+    rx = currentRow[intFloorX] + \
+         fracX * (currentRow[intFloorX+1] - \
+                  currentRow[intFloorX]); \
+    \
+    pixel = rx + fracY * ((psF64)nextRow[intFloorX] + \
+                          fracX * (nextRow[intFloorX+1] - \
+                                   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)
+PSIMAGE_PIXEL_INTERPOLATE_BILINEAR_COMPLEX(C32)
+PSIMAGE_PIXEL_INTERPOLATE_BILINEAR_COMPLEX(C64)
+
+