#ifdef HAVE_CONFIG_H
# include "config.h"
#endif

#include <stdio.h>

#include "psError.h"
#include "psImage.h"
#include "psAssert.h"
#include "psImageUnbin.h"

// 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)
{
    PS_ASSERT_IMAGE_NON_NULL(in, NULL);
    PS_ASSERT_IMAGE_NON_NULL(out, NULL);
    PS_ASSERT_INT_POSITIVE(DX, NULL);
    PS_ASSERT_INT_POSITIVE(DY, NULL);
    PS_ASSERT_INT_LESS_THAN_OR_EQUAL(dx, DX, NULL);
    PS_ASSERT_INT_LESS_THAN_OR_EQUAL(dy, DY, NULL);

    long nx = in->numCols;
    long ny = in->numRows;
    long Nx = out->numCols;
    long Ny = out->numRows;

    psF32 **vIn = in->data.F32;
    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 ++) {

            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 = 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++) {
                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;
                    //assert (fabs(V - psImageUnbinPixel(ix, iy, in, DX, DY, dx, dy)) < 1e-5*fabs(V));
                    V += dV;
                }
            }
        }
    }

    // side pixels
    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 = PS_MAX (0, PS_MIN (Ny, (Iy + 0.5)*DY - dy));
        int Ye = PS_MAX (0, PS_MIN (Ny, 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 = 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 = PS_MAX (0, PS_MIN (Nx, (Ix + 0.5)*DX - dx));
        int Xe = PS_MAX (0, PS_MIN (Nx, 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[ny-1][nx-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;
}

/************************************************************************************************************/
/*
 * 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_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 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;
    }
    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]
}
