Indoor/math/boxkde/BoxGaus3D.h

#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 = val + src(ri) - fv;

			// Fix potential float error
			if (val < 0) {
				val = 0;
			}

            dst(i) = val*iarr;
			ri++;
        }

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

			// Fix potential float error
			if (val < 0) {
				val = 0;
			}

            dst(i) = val*iarr;
			ri++;
			li++;
        }

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

			// Fix potential float error
			if (val < 0) {
				val = 0;
			}

            dst(i) = val*iarr;
			li++;
        }
    }
};


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