#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;
}

psStats *psImageClippedStats (psImage *image, psImage *mask, psU8 maskValue, double fmin, double fmax)
{
    double value;
    int nx = image->numCols;
    int ny = image->numRows;

    psStats *stats = psStatsAlloc (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_QUARTILE);
    if (nx*ny <= 0)
        return stats;

    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 && 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;

    stats->robustMedian = value;
    stats->robustUQ = values->data.F32[imax];
    stats->robustLQ = values->data.F32[imin];

    psFree (values);
    return stats;
    // XXX correct for selection bias??
}

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;
    }
}

