#include <render/rfilter.h>
#include <cmath>
#include <vector>

Float ReconstructionFilter::eval(const Vector2 &p, const Float *misc) const
{
    if (dynamic_cast<const TentFilter *>(this))
        return static_cast<const TentFilter *>(this)->eval(p, misc);
    if (dynamic_cast<const BoxFilter *>(this))
        return static_cast<const BoxFilter *>(this)->eval(p, misc);
    if (dynamic_cast<const AnisotropicGaussianFilter *>(this))
        return static_cast<const AnisotropicGaussianFilter *>(this)->eval(p, misc);
    assert(false);
    return 0.;
}


Float ReconstructionFilter::sample(const Vector2 &rnd, Vector2 &p, const Float *misc) const
{
    if (dynamic_cast<const TentFilter *>(this))
        return static_cast<const TentFilter *>(this)->sample(rnd, p, misc);
    if (dynamic_cast<const BoxFilter *>(this))
        return static_cast<const BoxFilter *>(this)->sample(rnd, p, misc);
    if (dynamic_cast<const AnisotropicGaussianFilter *>(this))
        return static_cast<const AnisotropicGaussianFilter *>(this)->sample(rnd, p, misc);
    assert(false);
    return 0.; 
}


// Tent filter

Float TentFilter::eval1D(Float x) const
{
    Float k = detach(0.5 * m_invPadding);
    Float y = (0.5 + m_padding) * k - k * abs(x - 0.5);
    return std::clamp(y, 0., 1.);
}


Float TentFilter::sample1D(Float rnd, Float &pdf) const
{
    Float s1 = m_padding;
    Float s2 = s1 + 1.0 - 2.0 * m_padding;
    Float ret;
    if (rnd <= s1)
        ret = std::sqrt(4.0 * m_padding * rnd) - m_padding;
    else if (rnd <= s2)
        ret = rnd - s1 + m_padding;
    else
    {
        rnd -= s2;
        ret = -std::sqrt(4.0 * m_padding * rnd) + m_padding + 1.0;
    }
    pdf = eval1D(ret);
    return ret;
}


Float TentFilter::eval(const Vector2 &p, const Float *misc) const
{
    assert( misc == nullptr );
    return eval1D(p[0])*eval1D(p[1]);
}


Float TentFilter::sample(const Vector2 &rnd, Vector2 &p, const Float *misc) const
{
    assert( misc == nullptr );
    Float ret = 1., pdf;
    p[0] = sample1D(rnd[0], pdf); ret *= pdf;
    p[1] = sample1D(rnd[1], pdf); ret *= pdf;
    return ret; 
}


// Box filter

Float BoxFilter::eval1D(Float x) const
{
    if (x > 0 && x < 1)
        return 1.0;
    else
        return 0.0;
}


Float BoxFilter::sample1D(Float rnd, Float &pdf) const
{
    pdf = 1.;
    return rnd;
}


Float BoxFilter::eval(const Vector2 &p, const Float *misc) const
{
    assert( misc == nullptr );
    return eval1D(p[0])*eval1D(p[1]);
}


Float BoxFilter::sample(const Vector2 &rnd, Vector2 &p, const Float *misc) const
{
    assert( misc == nullptr );
    Float ret = 1., pdf;
    p[0] = sample1D(rnd[0], pdf); ret *= pdf;
    p[1] = sample1D(rnd[1], pdf); ret *= pdf;
    return ret; 
}


// Anisotropic Gaussian Filter

Float AnisotropicGaussianFilter::eval(const Vector2 &p, const Float *misc) const
{
    assert( misc );
    Float sX = misc[0], sY = misc[1];

    assert( sX > Epsilon && sY > Epsilon );
    Float dx = (p[0] - .5)/sX, dy = (p[1] - .5)/sY;
    return std::exp(-.5*dx*dx)*std::exp(-.5*dy*dy)/(2.*M_PI*sX*sY);
}


Float AnisotropicGaussianFilter::sample(const Vector2 &rnd, Vector2 &p, const Float *misc) const
{
    assert( misc );
    Float sX = misc[0], sY = misc[1];

    assert( sX > Epsilon && sY > Epsilon );
    Float r = std::sqrt(-2.*std::log(rnd.x())), angle = 2.*M_PI*rnd.y();
    p[0] = sX*r*std::cos(angle) + 0.5;
    p[1] = sY*r*std::sin(angle) + 0.5;

    // Same as eval
    Float dx = (p[0] - .5)/sX, dy = (p[1] - .5)/sY;
    return std::exp(-.5*dx*dx)*std::exp(-.5*dy*dy)/(2.*M_PI*sX*sY);
}
