Index: /trunk/psLib/src/imageops/psImageUnbin.c
===================================================================
--- /trunk/psLib/src/imageops/psImageUnbin.c	(revision 11298)
+++ /trunk/psLib/src/imageops/psImageUnbin.c	(revision 11299)
@@ -12,4 +12,5 @@
 // interpolate from model to background (~bresenham linear interpolation)
 // XXX this code skips the initial pixels?
+// XXX this code may be off by 0.5,0.5 pixels
 psImage *psImageUnbin(psImage *out, const psImage *in, int DX, int DY, int dx, int dy)
 {
@@ -18,6 +19,4 @@
     PS_ASSERT_INT_POSITIVE(DX, NULL);
     PS_ASSERT_INT_POSITIVE(DY, NULL);
-    PS_ASSERT_INT_NONNEGATIVE(dx, NULL);
-    PS_ASSERT_INT_NONNEGATIVE(dy, NULL);
     PS_ASSERT_INT_LESS_THAN_OR_EQUAL(dx, DX, NULL);
     PS_ASSERT_INT_LESS_THAN_OR_EQUAL(dy, DY, NULL);
@@ -31,4 +30,5 @@
     psF32 **vOut = out->data.F32;
 
+    // loop over all pixels with
     for (int Iy = 0; Iy < ny-1; Iy ++) {
         for (int Ix = 0; Ix < nx-1; Ix ++) {
@@ -42,8 +42,8 @@
             // (Xs,Ys) : (Xe,Ye) : binned pixel centers in unbinned coords
             // corresponding to (Ix,Iy), (Ix+1,Iy+1)
-            int Xs = (Ix + 1)*DX - dx;
-            int Ys = (Iy + 1)*DY - dy;
-            int Xe = Xs + DX;
-            int Ye = Ys + DY;
+            int Xs = PS_MAX (0, PS_MIN (Nx, (Ix + 0.5)*DX - dx));
+            int Ys = PS_MAX (0, PS_MIN (Ny, (Iy + 0.5)*DY - dy));
+            int Xe = PS_MAX (0, PS_MIN (Nx, Xs + DX));
+            int Ye = PS_MAX (0, PS_MIN (Ny, Ys + DY));
 
             for (int iy = Ys; (iy < Ye) && (iy < Ny); iy++) {
@@ -62,10 +62,10 @@
 
     // side pixels
-    int Xs = DX - dx;
-    int Xe = nx*DX - dx;
+    int Xs = PS_MAX (0, PS_MIN (Nx, ( 0 + 0.5)*DX - dx));
+    int Xe = PS_MAX (0, PS_MIN (Nx, (nx - 0.5)*DX - dx));
     for (int Iy = 0; Iy < ny - 1; Iy++) {
 
-        int Ys = (Iy + 1)*DY - dy;
-        int Ye = Ys + DY;
+        int Ys = PS_MAX (0, PS_MIN (Ny, (Iy + 0.5)*DY - dy));
+        int Ye = PS_MAX (0, PS_MIN (Ny, Ys + DY));
 
         // leading edge
@@ -95,10 +95,10 @@
 
     // top and bottom pixels
-    int Ys = DY - dy;
-    int Ye = ny*DY - dy;
+    int Ys = PS_MAX (0, PS_MIN (Ny, ( 0 + 0.5)*DY - dy));
+    int Ye = PS_MAX (0, PS_MIN (Ny, (ny - 0.5)*DY - dy));
     for (int Ix = 0; Ix < nx - 1; Ix++) {
 
-        int Xs = (Ix + 1)*DX - dx;
-        int Xe = Xs + DX;
+        int Xs = PS_MAX (0, PS_MIN (Nx, (Ix + 0.5)*DX - dx));
+        int Xe = PS_MAX (0, PS_MIN (Nx, Xs + DX));
 
         // top edge
@@ -183,6 +183,4 @@
     PS_ASSERT_INT_POSITIVE(DX, NAN);
     PS_ASSERT_INT_POSITIVE(DY, NAN);
-    PS_ASSERT_INT_NONNEGATIVE(dx, NAN);
-    PS_ASSERT_INT_NONNEGATIVE(dy, NAN);
     PS_ASSERT_INT_LESS_THAN_OR_EQUAL(dx, DX, NAN);
     PS_ASSERT_INT_LESS_THAN_OR_EQUAL(dy, DY, NAN);
@@ -190,23 +188,24 @@
      * Find which binned pixel we're in
      */
-    const int Ix = (ix + dx)/DX - 1; // index of binned pixel
-    const int Iy = (iy + dy)/DY - 1;
-
-    if (Ix < 0 || Ix >= in->numCols - 1 || Iy < 0 || Iy >= in->numRows - 1) {
+    const int Xs = (ix + dx)/DX; // index of binned pixel
+    const int Ys = (iy + dy)/DY;
+    if (Xs < 0 || Xs >= in->numCols || Ys < 0 || Ys >= in->numRows) {
         psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Point (%d,%d) lies outside binned image", ix, iy);
         return NAN;
     }
-
-    double V00 = in->data.F32[Iy+0][Ix+0];
-    double V01 = in->data.F32[Iy+0][Ix+1];
-    double V10 = in->data.F32[Iy+1][Ix+0];
-    double V11 = in->data.F32[Iy+1][Ix+1];
-
-    const int Xs = (Ix + 1)*DX - dx; // centre of bottom left corner of binned pixel
-    const int Ys = (Iy + 1)*DY - dy; // (i.e. [Iy][Ix]) in unbinned coordinates
-
-    const double Vxs = (V10 - V00)*(iy - Ys)/DY + V00;
-    const double Vxe = (V11 - V01)*(iy - Ys)/DY + V01;
-
-    return Vxs + (Vxe - Vxs)*(ix - Xs)/DX; // value at [iy][ix]
+    const int Xe = PS_MIN (Xs + 1, in->numCols - 1);
+    const int Ye = PS_MIN (Ys + 1, in->numRows - 1);
+
+    double V00 = in->data.F32[Ys][Xs];
+    double V01 = in->data.F32[Ys][Xe];
+    double V10 = in->data.F32[Ye][Xs];
+    double V11 = in->data.F32[Ye][Xe];
+
+    const int xs = Xs*DX - dx; // centre of bottom left corner of binned pixel
+    const int ys = Ys*DY - dy; // (i.e. [Xs][Ys]) in unbinned coordinates
+
+    const double Vxs = (V10 - V00)*(iy - ys)/DY + V00;
+    const double Vxe = (V11 - V01)*(iy - ys)/DY + V01;
+
+    return Vxs + (Vxe - Vxs)*(ix - xs)/DX; // value at [iy][ix]
 }