#pragma once

#include "DataStructures.h"

template<typename T>
struct BoxSizes
{
    static_assert(std::is_floating_point<T>::value, "This class only works with floats.");

    T sigma, sigmaActual;
    T wIdeal, mIdeal;

    unsigned n, m, wl, wu;

    BoxSizes(T sigma, unsigned n)
        : sigma(sigma), n(n)
    {
        assertMsg(sigma >= 0.8, "Sigma values below about 0.8 cannot be represented");

        wIdeal = sqrt(12 * sigma*sigma / n + 1); // Ideal averaging filter width

        // wl is first odd valued integer less than wIdeal
        wl = (unsigned)floor(wIdeal);
        if (wl % 2 == 0)
            wl = wl - 1;

        // wu is the next odd value > wl
        wu = wl + 2;

        // Compute m.Refer to the tech note for derivation of this formula
        mIdeal = (12 * sigma*sigma - n*wl*wl - 4 * n*wl - 3 * n) / (-4 * wl - 4);
        m = (unsigned)round(mIdeal);

        assertMsg(!(m > n || m < 0), "calculation of m has failed");

        // Compute actual sigma that will be achieved
        sigmaActual = sqrt((m*wl*wl + (n - m)*wu*wu - n) / (T)12.0);
    }
};


template<typename T>
struct ExBoxSizes
{
    static_assert(std::is_floating_point<T>::value, "This class only works with floats.");

    T sigma, sigma_actual;
    unsigned n, r;
    T alpha, c, c1, c2;

    T r_f;

    // Special case for h == 1
    ExBoxSizes(T sigma, unsigned n)
        : sigma(sigma), n(n)
    {
        T var = sigma*sigma;
        r_f = 0.5*sqrt((12 * var) / n + 1) - T(0.5);
        r = (unsigned)std::floor(r_f);

        alpha = (2 * r + 1) * ( (r*r + r - 3*var/n) / (6 * (var/n - (r+1)*(r+1))) );
        
        c1 = alpha / (2*alpha + 2*r + 1);
        c2 = (1 - alpha) / (2 * alpha + 2 * r + 1);
        c = c1 + c2;
    }

    ExBoxSizes(T sigma, unsigned d, T h)
        : sigma(sigma), n(d)
    {
        T v = sigma*sigma;
        r_f = sqrt((12 * v) / n + 1)/(2*h) - T(0.5); // (7)
        r = (unsigned)std::floor(std::max(T(0), r_f));

        alpha = (2 * r + 1) * (((r*r + r) - (v*3*h)/(d*h*h*h)) / (6 * ( v/(d*h*h) - (r+1)*(r+1)) )); // (8) (14)
        c1 = alpha / (h*(2 * r + 2 * alpha + 1));       // (8) (13)
        c2 = (1 - alpha) / (h*(2 * r + 2 * alpha + 1)); // (8) (13)

        c = c1 + c2;

        T lambda = h*(2*r + 1 + 2*alpha); // (8)
        T variance_actual = (d*h*h*h) / (3 * lambda) * (2*r*r*r + 3*r*r + r + 6*alpha*(r+1)*(r+1)); // (14)
        sigma_actual = sqrt(variance_actual);
    }
};