//
// Source code for the paper
// 
// D. Heck and T. Schloemer and O. Deussen, "Blue Noise Sampling with
// Controlled Aliasing", ACM Trans. Graph., 2013, in press
//
// Copyright (C) 2012,2013 Daniel Heck and Thomas Schloemer
//
#include "common.hh"
#include <fstream>
#include <iostream>
#include <string.h>
#include <cstdlib>

// Precomputed arrays for Bessel function computation
#define PI_4    0.7853981633974483096f
#define DR1     5.7831859629467845211f

static const float JP[5] = {
    -6.068350350393235E-008f,
    6.388945720783375E-006f,
    -3.969646342510940E-004f,
    1.332913422519003E-002f,
    -1.729150680240724E-001f
};
static const float MO[8] = {
    -6.838999669318810E-002f,
    1.864949361379502E-001f,
    -2.145007480346739E-001f,
    1.197549369473540E-001f,
    -3.560281861530129E-003f,
    -4.969382655296620E-002f,
    -3.355424622293709E-006f,
    7.978845717621440E-001f
};
static const float PH[8] = {
    3.242077816988247E+001f,
    -3.630592630518434E+001f,
    1.756221482109099E+001f,
    -4.974978466280903E+000f,
    1.001973420681837E+000f,
    -1.939906941791308E-001f,
    6.490598792654666E-002f,
    -1.249992184872738E-001f
};

// Single precision approximation to j0, about twice as fast as the standard
// implementation. From cephes lib, available at http://www.netlib.org/cephes/
inline float Polynomial(float x, const float *p, int n) {
    float y = p[0];
    for (int i = 1; i < n; ++i)
        y = y * x + p[i];
    return y;
}

inline float j0f(float f) {
    float x = fabsf(f);
    if (x <= 2.f) {
        float z = x*x;
        if (x < 1.0e-3f)
            return (1.f - 0.25f*z);
        return (z - DR1) * Polynomial(z, JP, 5);
	}
    float q = 1.f / x;
    float w = sqrtf(q);
    float p = w * Polynomial(q, MO, 8);
    w = q*q;
    float xn = q * Polynomial(w, PH, 8) - PI_4;
    return p * cosf(xn + x);
}

bool HasSuffix(const std::string &s, const std::string &suffix) {
    int i = s.size() - suffix.size();
    if (i >= 0)
        return s.substr(i) == suffix;
    return false;
}

Endianness SystemEndianness() {
    union { uint32_t i; char c[4]; } e;
    e.i = 1;
    return e.c[0] == 0 ? BigEndian : LittleEndian;
}

void SwapEndian4(uint8_t *buf, int nn) {
    for (int i=0; i<nn; i++) {
        std::swap(buf[4*i], buf[4*i+3]);
        std::swap(buf[4*i+1], buf[4*i+2]);
    }
}

bool WriteFloats(FILE *fp, const float *p, int n, Endianness endian) {
    Endianness sysEndian = SystemEndianness();
    const int nbuf = 256;
    uint8_t buf[nbuf*sizeof(float)];
    while (n > 0) {
        int nn = std::min(n, nbuf);
        memcpy(buf, p, nn * sizeof(float));
        if (endian != sysEndian)
            SwapEndian4(buf, nn);
        if (fwrite(buf, sizeof(float), nn, fp) != (size_t)nn)
            return false;
        p += nn; n -= nn;
    }
    return true;
}

bool ReadFloats(FILE *fp, float *p, int n, Endianness endian) {
    Endianness sysEndian = SystemEndianness();
    const int nbuf = 256;
    uint8_t buf[nbuf*sizeof(float)];
    while (n > 0) {
        int nn = std::min(n, nbuf);
        if (fread(buf, sizeof(float), nn, fp) != (size_t)nn)
            return false;
        if (endian != sysEndian)
            SwapEndian4(buf, nn);
        memcpy(p, buf, nn*sizeof(float));
        p += nn; n -= nn;
    }
    return true;
}

void ReadPFM(int &width, int &height, float *&data, const char *fname) {
    FILE* fp = fopen(fname, "rb");
    char name[100];
    float f;
    fscanf(fp, "%s\n%d %d %f\n", name, &width, &height, &f);
    printf("f is %f\n", f);
    data = new float[width * height];
    fread((char*) data, sizeof(float), width * height, fp);
    fclose(fp);
}

void WritePFM(int width, int height, const float *data, const char *fname) {
    FILE* fp = fopen(fname, "wb");
    fprintf(fp, "Pf\n");
    fprintf(fp, "%d %d\n-1.00\n", width, height);
    fwrite((const char*) data, sizeof(float), width * height, fp);
    fclose(fp);
}


void LoadPoints(const std::string &fname, std::vector<Point> &points) {
    int npts = 0;
    if (HasSuffix(fname, ".rps")) {
        FILE *fp= fopen(fname.c_str(), "rb");
        if (!fp) {
            std::cerr << "Cannot load '" << fname << "'.\n";
            exit(1);
        }
        fseek(fp, 0, SEEK_END);
        npts = ftell(fp) / 8;
        points.resize(npts);
        fseek(fp, 0, SEEK_SET);
        ReadFloats(fp, &points[0].x, npts*2, LittleEndian);
        fclose(fp);
    } else {
        std::ifstream fp(fname.c_str());
        if (!fp) {
            std::cerr << "Cannot load '" << fname << "'.\n";
            exit(1);
        }
        fp >> npts;
        points.reserve(npts);
        Point pt;
        while ((int)points.size() < npts && (fp >> pt.x >> pt.y))
            points.push_back(pt);
    }
}

void WritePoints(const std::string &fname, std::vector<Point> &points) {
    std::ofstream fp(fname.c_str());
    if (!fp) {
        std::cerr << "Cannot create '" << fname << "'.\n";
        exit(1);
    }
    fp << points.size() << "\n";
    for (unsigned i = 0; i < points.size(); i++)
        fp << points[i].x << " " << points[i].y << "\n"; 
}


void SavePGM(float *p, int rows, int cols, const char *fname) {
    FILE *fp = fopen(fname, "wb");
    if (!fp) {
        fprintf(stderr, "Cannot open '%s' for writing.\n", fname);
        exit(1);
    }

    int maxval = 255;
    fprintf(fp, "P2\n%d %d\n%d\n", cols, rows, maxval);

    for (int y=0; y<rows; y++) {
        for (int x=0; x<cols; x++) {
            int pp = (int)(p[y*cols + x] * maxval);
            if (pp < 0) pp = 0;
            if (pp > maxval) pp = maxval;
            fprintf(fp, "% 3d ", pp);
        }
        fprintf(fp, "\n");
    }

    fclose(fp);
}


// -------------------- Curve --------------------

Curve::Curve(const Curve &c) {
    x0 = c.x0; 
    x1 = c.x1;
    dx = c.dx;
    y = c.y;
}

Curve::Curve(int size, float x0, float x1) {
    y.resize(size);
    
    this->x0 = x0;
    this->x1 = x1;
    dx = 1.f/size * (x1 - x0);
}

float Curve::At(float x) const {
    float xx = (x-x0)/dx;
    if (xx < 0) xx = 0;
    int i = (int)floor(xx);
    int ii = i+1;
    float a = xx-i;
    if (i >= (int)y.size()) i = y.size()-1;
    if (ii >= (int)y.size()) ii = y.size()-1;
    return (1-a)*y[i] + a * y[ii];
}

void Curve::Write(const std::string &fname) {
    std::ofstream os(fname.c_str());
    for (int j=0; j<size(); j++)
        os << ToX(j) << " " << y[j] << "\n";
}

Curve Curve::Read(const std::string &fname) {
    std::ifstream is(fname.c_str());
    std::vector<float> x, y;
    float xx, yy;
    while (is >> xx >> yy) {
        x.push_back(xx);
        y.push_back(yy);
    }
    if (x.empty()) {
        std::cerr << "Could not load RDF from '" << fname << "'.\n";
        exit(1);
    }

    int n = x.size();
    float dx = n > 1 ? x[1] - x[0] : 1;
    Curve c(n, x[0]-dx/2, x[n-1] + dx/2);
    for (int i=0; i<n; i++)
        c.y[i] = y[i];
    return c;
}



void Curve::FilterGauss(const Curve &source, float sigma) {
    if (sigma == 0){
        y = source.y;
        return;
    }
  
    float dd = (float)source.size() / size();
  
    for (int i=0; i<size(); i++) {
        float a = 0.0f, sumw = 0.0f;

        int jmin = std::max(0, (int)floor(dd * (i - 5 * sigma)));
        int jmax = std::min((int)source.size()-1, (int)ceil(dd * (i + 5*sigma)));
        for (int j=jmin; j<=jmax; j++) {
            float d = j/dd-i;
            float w = exp(-d*d/(2*sigma*sigma));
            a += source[j] * w;
            sumw += w;
        }
        y[i] = a / sumw;
    }
}

float Integrate(const Curve &c, float x0, float x1, float *E) {
    bool negate = false;
    if (x0 > x1) {
        std::swap(x0, x1);
        negate = true;
    }

    float T = 0.0f, M = 0.0f;
    int i0 = c.ToIndex(x0);
    int i1 = c.ToIndex(x1+c.dx);
    if (i0 < 0) i0 = 0;
    if (i1 < 0) i1 = 0;
    if (i0 >= c.size()) i0 = c.size()-1;
    if (i1 >= c.size()) i1 = c.size()-1;

    T = 0.5 * (c[i0] + c[i1]);
    for (int i=i0+2; i<i1; i+=2) {
        T += c[i];
        M += c[i-1];
    }
    if (E != NULL)
        *E = c.dx * 2.0f * fabs(T-M);
    float a = c.dx*(T+M);
    return negate ? -a : a;
}


// -------------------- RDF Calculation --------------------

void CalcRDF(Curve &c, size_t npts, const float *pts, float (*distfunc)(const float *, const float *)) {
    std::vector<unsigned long> bins(c.size());
    for (unsigned j=0; j<npts; j++) {
        for (unsigned k=j+1; k<npts; k++) {
            float d = distfunc(&pts[2*j], &pts[2*k]);
            int idx = c.ToIndex(d);
            if (0 <= idx && idx < c.size())
                bins[idx]++;
//            if (d < c.x1)
//                bins[c.ToIndex(d)]++;
        }
    }
    const float scale = npts * (npts-1)/2 * M_PI * c.dx * c.dx;
    for (int j=0; j<c.size(); j++)
        c[j] = bins[j] / (scale * (2*j+1));
}

Curve CalcRDF(int numbins, size_t npts, const float *pts, 
        float sigma, float (*distfunc)(const float *, const float *))
{
    Curve c(numbins, 0, 0.5f);
    CalcRDF(c, npts, pts, distfunc);
    if (sigma == 0)
        return c;
    else
        return FilterGauss(numbins, c, sigma);
}

float HannWindow(float x, float xlim) {
    return x > xlim ? 0.0f : Sqr(cosf(M_PI * x / (2*xlim)));
}

float BlackmanWindow(float x, float xlim) {
    return x > xlim ? 0.0f : 0.43 + 0.5*cosf(M_PI*x/xlim) + 0.08*cosf(2*M_PI*x/xlim);
}

float HammingWindow(float x, float xlim) {
    return x > xlim ? 0.0f : 0.54f + 0.46f * cosf(M_PI * x/xlim);
}

float WelchWindow(float x, float xlim) {
    return x > xlim ? 0.0f : 1 - x*x/(xlim*xlim);
}

Curve RDF2Power(int npts, const Curve &rdf, int nbins, float x0, float x1) {
    Curve power(nbins, x0, x1);
    Curve tmp(rdf);
    for (int i=0; i<nbins; i++) {
        const float u0 = power.ToX(i);
        const float u = 2*M_PI*u0;
        const float wndsize = rdf.x1 * std::min(0.5f, std::max(0.2f, 4*u0/sqrtf(npts)));
        for (int j=0; j<tmp.size(); j++) {
            float x = rdf.ToX(j);
            float wnd = BlackmanWindow(x, wndsize);
//            float wnd = HannWindow(x, wndsize);
            tmp[j] = (rdf[j]-1) * j0f(u*x) * x * wnd;
        }
        power[i] = fabs(1 + 2 * M_PI * Integrate(tmp) * npts);
    }
    return power;
}

Curve Power2RDF(int npts, const Curve &power, int nbins, float x0, float x1, float smoothing) {
    Curve rdf(nbins, x0, x1);
    Curve tmp(power);
    for (int i=0; i<rdf.size(); i++) {
        float u = rdf.ToX(i);
        for (int j=0; j<tmp.size(); j++) {
            float x = tmp.ToX(j);
            float wnd = BlackmanWindow(x, power.x1);
            tmp[j] = (power[j]-1) * jn(0, 2*M_PI * u * x) * x * wnd;
        }
        rdf[i] = 1 + 2 * M_PI * Integrate(tmp) / npts;
    }
    for (int i=0; i<rdf.size(); i++) {
        if (rdf[i] < 0) {
            printf("Warning: RDF has negative components\n");
            break;
        }
    }
    return FilterGauss(nbins, rdf, smoothing);
}