#pragma once

#include <cmath>
#include <vector>

#include "BoxSizes.h"
#include "Image3D.h"

template <class T>
struct BoxGaus3D
{
    void boxFilter(std::vector<T>& input, unsigned w, unsigned h, unsigned d, unsigned filterSize)
    {
        std::vector<T> buffer(input.size());

        boxFilter(input.data(), buffer.data(), w, h, d, filterSize);

        input.assign(buffer.begin(), buffer.end());
    }

    void boxFilter(T* input, T* output, unsigned w, unsigned h, unsigned d, unsigned filterSize)
    {
        assertMsg((filterSize % 2 == 1), "filterSize must be odd");
        unsigned radius = filterSize / 2;

        ImageView3D<T> in (w, h, d, input);
        ImageView3D<T> out(w, h, d, output);

        boxBlur(in, out, radius, radius, radius);
    }

	void approxGaus(Image3D<T>& input, _Point3<T> sigma, unsigned nFilt)
	{
		approxGaus(input.data(), input.getWidth(), input.getHeight(), input.getDepth(), sigma, nFilt);
	}

    void approxGaus(std::vector<T>& input, unsigned w, unsigned h, unsigned d, _Point3<T> sigma, unsigned nFilt)
    {
        std::vector<T> buffer(input.size());

        approxGaus(input.data(), buffer.data(), w, h, d, sigma, nFilt);

        input.assign(buffer.begin(), buffer.end());
    }

    void approxGaus(T* input, T* output, unsigned w, unsigned h, unsigned d, _Point3<T> sigma, unsigned nFilt)
    {
        ImageView3D<T> in (w, h, d, input);
        ImageView3D<T> out(w, h, d, output);

        approxGaus(in, out, sigma, nFilt);
    }

    void approxGaus(ImageView3D<T>& in, ImageView3D<T>& out, _Point3<T> sigma, unsigned nFilt)
    {
		BoxSizes<T> bsX(sigma.x, nFilt);
		BoxSizes<T> bsY(sigma.y, nFilt);
		BoxSizes<T> bsZ(sigma.z, nFilt);
		
		assertMsg((nFilt % 2 == 1), "nFilt must be odd");

		// if equal, we can save some cond's inside the loop
		if (bsX.m == bsY.m && bsX.m == bsZ.m)
		{
			const size_t m = bsX.m;

			for (size_t i = 1; i <= m; i++)
			{
				if (i % 2) {
					boxBlur(in, out, bsX.wl, bsY.wl, bsZ.wl);
				} else {
					boxBlur(out, in, bsX.wl, bsY.wl, bsZ.wl);
				}
			}

			for (size_t i = 1; i <= nFilt - m; i++)
			{
				if (i % 2) {
					boxBlur(in, out, bsX.wu, bsY.wu, bsZ.wu);
				} else {
					boxBlur(out, in, bsX.wu, bsY.wu, bsZ.wu);
				}
			}
		}
		else
		{
			for (size_t i = 1; i <= nFilt; i++)
			{
				const size_t rX = (i < bsX.m ? bsX.wl : bsX.wu);
				const size_t rY = (i < bsY.m ? bsY.wl : bsY.wu);
				const size_t rZ = (i < bsZ.m ? bsZ.wl : bsZ.wu);

				if (i % 2) {
					boxBlur(in, out, rX, rY, rZ);
				} else{
					boxBlur(out, in, rX, rY, rZ);
				}
			}
		}
    }

private:
    void boxBlur(ImageView3D<T>& in, ImageView3D<T>& out, size_t rX, size_t rY, size_t rZ)
    {
        //box(in, out, r);
        boxBlurX(in, out, rX);
        boxBlurY(out, in, rY);
        boxBlurZ(in, out, rZ);
    }

    void boxBlurX(const ImageView3D<T>& src, ImageView3D<T>& dst, size_t r)
    {
        T iarr = (T)1.0 / (r + r + 1);

        for (size_t z = 0; z < src.getDepth(); z++)
        {
            for (size_t y = 0; y < src.getHeight(); y++)
            {
                auto srcView = src.getViewYZ(y, z);
                auto dstView = dst.getViewYZ(y, z);

                boxBlur1D<SlicePlane::YZ>(srcView, dstView, src.getWidth(), r, iarr);
            }
        }
    }

    void boxBlurY(const ImageView3D<T>& src, ImageView3D<T>& dst, size_t r)
    {
        T iarr = (T)1.0 / (r + r + 1);

        for (size_t x = 0; x < src.getWidth(); x++)
        {
            for (size_t z = 0; z < src.getDepth(); z++)
            {
                auto srcView = src.getViewXZ(x, z);
                auto dstView = dst.getViewXZ(x, z);

                boxBlur1D<SlicePlane::XZ>(srcView, dstView, src.getHeight(), r, iarr);
            }
        }
    }

    void boxBlurZ(const ImageView3D<T>& src, ImageView3D<T>& dst, size_t r)
    {
        T iarr = (T)1.0 / (r + r + 1);

        for (size_t x = 0; x < src.getWidth(); x++)
        {
            for (size_t y = 0; y < src.getHeight(); y++)
            {
                auto srcView = src.getViewXY(x, y);
                auto dstView = dst.getViewXY(x, y);

                boxBlur1D<SlicePlane::XY>(srcView, dstView, src.getDepth(), r, iarr);
            }
        }
    }


    template <SlicePlane S>
    void boxBlur1D(const typename ImageView3D<T>::template ConstLineView<S>& src, typename ImageView3D<T>::template LineView<S>& dst, size_t len, size_t r, T iarr)
    {
        size_t li = 0; // left index
        size_t ri = r; // right index

        T fv = src(0);
        T lv = src(len - 1);
        T val = (r + 1)*fv;

        for (size_t i = 0; i < r; i++)
            val += src(i);

        // Überhangbereich links vom Bild
        for (size_t i = 0; i <= r; i++)
        {
            val += src(ri++) - fv;
            dst(i) = val*iarr;
        }

        // Bildbereich
        for (size_t i = r + 1; i < len - r; i++)
        {
            val += src(ri++) - src(li++);
            dst(i) = val*iarr;
        }

        // Überhangbereich rechts vom Bild
        for (size_t i = len - r; i < len; i++)
        {
            val += lv - src(li++);
            dst(i) = val*iarr;
        }
    }

    void box(const ImageView3D<T>& src, ImageView3D<T>& dst, size_t r)
    {
        T iarr = (T)1.0 / (r + r + 1);

        ImageView3D<T> buffer(src);
        ImageView3D<T> buffer2(dst.getHeight(), dst.getDepth(), dst.getWidth(), dst.val_begin());

        for (size_t i = 0; i < 3; i++)
        {
            for (size_t idx2 = 0; idx2 < buffer.getHeight(); idx2++)
            {
                for (size_t idx3 = 0; idx3 < buffer.getDepth(); idx3++)
                {
                    box1D(buffer, buffer2, buffer.getWidth(), idx2, idx3, r, iarr);
                }
            }

            // reinteprete buffer sizes after rotation and swap data pointers
            buffer  = ImageView3D<T>(buffer.getHeight(), buffer.getDepth(), buffer.getWidth(), buffer.val_begin());
            buffer2 = ImageView3D<T>(buffer2.getHeight(), buffer2.getDepth(), buffer2.getWidth(), buffer2.val_begin());
            buffer.swap(buffer2);
        }

        buffer.swap(dst);
    }

    void box1D(const ImageView3D<T>& src, ImageView3D<T>& dst, size_t len, size_t idx2, size_t idx3, size_t r, T iarr)
    {
        for (size_t i = 0; i < len; i++)
        {
            // dst(idx2, idx3, i) = src(i, idx2, idx3);

            size_t li = 0; // left index
            size_t ri = r; // right index

            T fv = src(0, idx2, idx3);
            T lv = src(len - 1, idx2, idx3);
            T val = (r + 1)*fv;

            for (size_t i = 0; i < r; i++)
                val += src(i, idx2, idx3);

            // Überhangbereich links vom Bild
            for (size_t i = 0; i <= r; i++)
            {
                val += src(ri++, idx2, idx3) - fv;
                dst(idx2, idx3, i) = val*iarr;
            }

            // Bildbereich
            for (size_t i = r + 1; i < len - r; i++)
            {
                val += src(ri++, idx2, idx3) - src(li++, idx2, idx3);
                dst(idx2, idx3, i) = val*iarr;
            }

            // Überhangbereich rechts vom Bild
            for (size_t i = len - r; i < len; i++)
            {
                val += lv - src(li++, idx2, idx3);
                dst(idx2, idx3, i) = val*iarr;
            }
        }
    }
};


template struct BoxGaus3D<float>;
template struct BoxGaus3D<double>;