Index: /branches/rel10_ifa/psModules/src/pslib/Makefile.am
===================================================================
--- /branches/rel10_ifa/psModules/src/pslib/Makefile.am	(revision 6825)
+++ /branches/rel10_ifa/psModules/src/pslib/Makefile.am	(revision 6826)
@@ -8,4 +8,5 @@
     psImageJpeg.c \
     psImageFlip.c \
+    psImageUnbin.c \
     psLine.c \
     psMetadataItemParse.c \
@@ -20,4 +21,5 @@
     psImageJpeg.h \
     psImageFlip.h \
+    psImageUnbin.h \
     psLine.h \
     psMetadataItemParse.h \
Index: /branches/rel10_ifa/psModules/src/pslib/psImageUnbin.c
===================================================================
--- /branches/rel10_ifa/psModules/src/pslib/psImageUnbin.c	(revision 6826)
+++ /branches/rel10_ifa/psModules/src/pslib/psImageUnbin.c	(revision 6826)
@@ -0,0 +1,269 @@
+#include <stdio.h>
+#include "pslib.h"
+#include "psImageUnbin.h"
+
+// interpolate from model to background (~bresenham linear interpolation)
+// XXX this code skips the initial pixels?
+psImage *psImageUnbin (psImage *out, psImage *in, int DX, int DY, int dx, int dy)
+{
+
+    int nx = in->numCols;
+    int ny = in->numRows;
+
+    if (out == NULL)
+        psAbort ("psImageUnbin", "please supply output image");
+
+    int Nx = out->numCols;
+    int Ny = out->numRows;
+
+    psF32 **vIn = in->data.F32;
+    psF32 **vOut = out->data.F32;
+
+    for (int Iy = 0; Iy < ny-1; Iy ++) {
+        for (int Ix = 0; Ix < nx-1; Ix ++) {
+
+            float V00 = vIn[Iy+0][Ix+0];
+            float V01 = vIn[Iy+0][Ix+1];
+            float V10 = vIn[Iy+1][Ix+0];
+            float V11 = vIn[Iy+1][Ix+1];
+
+            // a single binned pixel quad
+            // (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;
+
+            for (int iy = Ys; (iy < Ye) && (iy < Ny); iy++) {
+                float Vxs = (V10 - V00)*(iy - Ys) / DY + V00;
+                float Vxe = (V11 - V01)*(iy - Ys) / DY + V01;
+                float dV = (Vxe - Vxs) / DX;
+                float V  = Vxs;
+                for (int ix = Xs; (ix < Xe) && (ix < Nx); ix++) {
+                    vOut[iy][ix] = V;
+                    V += dV;
+                }
+            }
+        }
+    }
+
+    // side pixels
+    int Xs = DX - dx;
+    int Xe = nx*DX - dx;
+    for (int Iy = 0; Iy < ny - 1; Iy++) {
+
+        int Ys = (Iy + 1)*DY - dy;
+        int Ye = Ys + DY;
+
+        // leading edge
+        float V0 = vIn[Iy+0][0];
+        float V1 = vIn[Iy+1][0];
+        float dV = (V1 - V0) / DY;
+        float V = V0;
+        for (int iy = Ys; (iy < Ye) && (iy < Ny); iy++) {
+            for (int ix = 0; ix < Xs; ix++) {
+                vOut[iy][ix] = V;
+            }
+            V += dV;
+        }
+
+        // trailing edge
+        V0 = vIn[Iy+0][nx-1];
+        V1 = vIn[Iy+1][nx-1];
+        dV = (V1 - V0) / DY;
+        V = V0;
+        for (int iy = Ys; (iy < Ye) && (iy < Ny); iy++) {
+            for (int ix = Xe; ix < Nx; ix++) {
+                vOut[iy][ix] = V;
+            }
+            V += dV;
+        }
+    }
+
+    // top and bottom pixels
+    int Ys = DY - dy;
+    int Ye = ny*DY - dy;
+    for (int Ix = 0; Ix < nx - 1; Ix++) {
+
+        int Xs = (Ix + 1)*DX - dx;
+        int Xe = Xs + DX;
+
+        // top edge
+        float V0 = vIn[0][Ix+0];
+        float V1 = vIn[0][Ix+1];
+        float dV = (V1 - V0) / DX;
+        for (int iy = 0; iy < Ys; iy++) {
+            float V = V0;
+            for (int ix = Xs; (ix < Xe) && (ix < Nx); ix++) {
+                vOut[iy][ix] = V;
+                V += dV;
+            }
+        }
+
+        // bottom edge
+        V0 = vIn[ny-1][Ix+0];
+        V1 = vIn[ny-1][Ix+1];
+        dV = (V1 - V0) / DX;
+        for (int iy = Ye; iy < Ny; iy++) {
+            float V = V0;
+            for (int ix = Xs; (ix < Xe) && (ix < Nx); ix++) {
+                vOut[iy][ix] = V;
+                V += dV;
+            }
+        }
+    }
+
+    // the four corners
+    {
+        float V;
+        // 0,0
+        V = vIn[0][0];
+        for (int iy = 0; iy < DY - dy; iy++)
+        {
+            for (int ix = 0; ix < DX - dx; ix++) {
+                vOut[iy][ix] = V;
+            }
+        }
+        // Nx,0
+        V = vIn[0][nx-1];
+        for (int iy = 0; iy < DY - dy; iy++)
+        {
+            for (int ix = nx*DX - dx; ix < Nx; ix++) {
+                vOut[iy][ix] = V;
+            }
+        }
+        // 0,Ny
+        V = vIn[ny-1][0];
+        for (int iy = ny*DY - dy; iy < Ny; iy++)
+        {
+            for (int ix = 0; ix < DX - dx; ix++) {
+                vOut[iy][ix] = V;
+            }
+        }
+        // Nx,Ny
+        V = vIn[nx-1][ny-1];
+        for (int iy = ny*DY - dy; iy < Ny; iy++)
+        {
+            for (int ix = nx*DX - dx; ix < Nx; ix++) {
+                vOut[iy][ix] = V;
+            }
+        }
+    }
+
+    return out;
+}
+
+// random number seed to select a fraction of the image pixels
+static psRandom *rnd = NULL;
+
+// use no more than MAX_SAMPLE_PIXELS pixels for each median box
+static int MAX_SAMPLE_PIXELS;
+
+void psImageClippedStatsInit (int nMax)
+{
+
+    MAX_SAMPLE_PIXELS = nMax;
+
+    if (!rnd) {
+        rnd = psRandomAlloc (PS_RANDOM_TAUS, 0);
+    }
+    return;
+}
+
+void psImageClippedStatsCleanup ()
+{
+    psFree (rnd);
+    rnd = NULL;
+    return;
+}
+
+double psImageClippedStats (psImage *image, psImage *mask, psU8 maskValue, double fmin, double fmax)
+{
+    double value;
+    int nx = image->numCols;
+    int ny = image->numRows;
+
+    if (nx*ny <= 0)
+        return 0.0;
+
+    int Nsubset = PS_MIN (MAX_SAMPLE_PIXELS, nx*ny);
+    int Npixels = nx*ny;
+
+    psVector *values = psVectorAlloc (Nsubset, PS_TYPE_F32);
+
+    float min = values->data.F32[0];
+    float max = values->data.F32[0];
+
+    int n = 0;
+    for (int i = 0; i < Nsubset; i++) {
+        double frnd = psRandomUniform (rnd);
+        int pixel = Npixels * frnd;
+        int ix = pixel % nx;
+        int iy = pixel / nx;
+
+        if (mask->data.U8[iy][ix] & maskValue)
+            continue;
+
+        value = image->data.F32[iy][ix];
+        min = PS_MIN (value, min);
+        max = PS_MIN (value, max);
+        values->data.F32[n] = value;
+        n++;
+    }
+    values->n = n;
+
+    int imin = fmin * n;
+    int imax = fmax * n;
+    int npts = imax - imin + 1;
+
+    fsort (values->data.F32, n);
+
+    value = 0;
+    for (int i = imin; (i <= imax) && (i < n); i++) {
+        value += values->data.F32[i];
+    }
+    value = value / npts;
+
+    // XXX correct for selection bias??
+
+    psFree (values);
+    return value;
+}
+
+void fsort (float *value, int N)
+{
+
+    int l,j,ir,i;
+    float temp;
+
+    if (N < 2)
+        return;
+    l = N >> 1;
+    ir = N - 1;
+    for (;;) {
+        if (l > 0) {
+            temp = value[--l];
+        } else {
+            temp = value[ir];
+            value[ir] = value[0];
+            if (--ir == 0) {
+                value[0] = temp;
+                return;
+            }
+        }
+        i = l;
+        j = (l << 1) + 1;
+        while (j <= ir) {
+            if (j < ir && value[j] < value[j+1])
+                ++j;
+            if (temp < value[j]) {
+                value[i]=value[j];
+                j += (i=j) + 1;
+            } else
+                j = ir + 1;
+        }
+        value[i] = temp;
+    }
+}
+
Index: /branches/rel10_ifa/psModules/src/pslib/psImageUnbin.h
===================================================================
--- /branches/rel10_ifa/psModules/src/pslib/psImageUnbin.h	(revision 6826)
+++ /branches/rel10_ifa/psModules/src/pslib/psImageUnbin.h	(revision 6826)
@@ -0,0 +1,18 @@
+#ifndef PS_IMAGE_UNBIN_H
+#define PS_IMAGE_UNBIN_H
+
+#include "pslib.h"
+
+// This needs to be considered more carefully
+psImage *psImageUnbin (psImage *out, psImage *in, int DX, int DY, int dx, int dy);
+
+void psImageClippedStatsInit (int nMax);
+void psImageClippedStatsCleanup ();
+
+// my temporary image stats function; seems to be much faster than psLib???
+double psImageClippedStats (psImage *image, psImage *mask, psU8 maskValue, double fmin, double fmax);
+
+// my temporary sort which is based on N.R.; seems to be faster than psLib sort???
+void fsort (float *value, int N);
+
+#endif
Index: /branches/rel10_ifa/psModules/src/psmodules.h
===================================================================
--- /branches/rel10_ifa/psModules/src/psmodules.h	(revision 6825)
+++ /branches/rel10_ifa/psModules/src/psmodules.h	(revision 6826)
@@ -8,4 +8,5 @@
 #include <psImageJpeg.h>
 #include <psImageFlip.h>
+#include <psImageUnbin.h>
 #include <psLine.h>
 #include <psPolynomialUtils.h>