Index: trunk/psLib/src/imageops/psImageUnbin.c
===================================================================
--- trunk/psLib/src/imageops/psImageUnbin.c	(revision 7766)
+++ trunk/psLib/src/imageops/psImageUnbin.c	(revision 10136)
@@ -50,4 +50,5 @@
                 for (int ix = Xs; (ix < Xe) && (ix < Nx); ix++) {
                     vOut[iy][ix] = V;
+                    //assert (fabs(V - psImageUnbinPixel(ix, iy, in, DX, DY, dx, dy)) < 1e-5*fabs(V));
                     V += dV;
                 }
@@ -161,2 +162,47 @@
     return out;
 }
+
+/************************************************************************************************************/
+/*
+ * Get the value of a single unbinned pixel from the binned representation
+ *
+ * N.b. This code only works for the central part of the image; the edge
+ * cases should be added
+ */
+double psImageUnbinPixel(const int ix, const int iy, // desired Unbinned point
+                         const psImage *in, // binned image
+                         const int DX, const int DY,  //!< Scaling factors in x and y
+                         const int dx, const int dy)   //!< Overhang
+{
+    PS_ASSERT_IMAGE_NON_NULL(in, NAN);
+    assert (in->type.type == PS_TYPE_F32);
+    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);
+    /*
+     * 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) {
+        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]
+}