#ifndef STIM_CUDA_GAUSS3_H
#define STIM_CUDA_GAUSS3_H
#include <stim/math/filters/sepconv3.cuh>
#include <stim/math/filters/gauss2.cuh>
#include <stim/math/constants.h>

namespace stim
{
	///Perform a 3D gaussian convolution on an input image.
        ///@param in is a pointer to the input data.
	///@param dimx is the size of in* in the x direction.
	///@param dimx is the size of in* in the y direction.
	///@param dimx is the size of in* in the z direction.
        ///@param stdx is the standard deviation (in pixels) along the x axis.
        ///@param stdy is the standard deviation (in pixels) along the y axis.
        ///@param nstds specifies the number of standard deviations of the Gaussian that will be kept in the kernel.
	template<typename T, typename K>
	void cpu_gauss3(T* in, K dimx, K dimy, K dimz, K stdx, K stdy, K stdz, size_t nstds = 3)
	{
		//Set up the sizes of the gaussian Kernels.
		size_t kx = stdx * nstds * 2;
		size_t ky = stdy * nstds * 2;
		size_t kz = stdz * nstds * 2;
	
		//Set up the sizes of the new output, which will be kx, ky, kz, smaller than the input.
		size_t X = dimx - kx +1; 
		size_t Y = dimy - ky +1; 
		size_t Z = dimz - kz +1; 
		T* out = (T*) malloc(X*Y*Z* sizeof(T));

		///Set up the memory that will store the gaussians
		K* gaussx = (K*)malloc(kx *sizeof(K));
		K* gaussy = (K*)malloc(ky *sizeof(K));
		K* gaussz = (K*)malloc(kz *sizeof(K));

		///Set up the midpoints of the gaussians.
		K midgaussx = (K) kx/ (K)2;
		K midgaussy = (K) ky/ (K)2;
		K midgaussz = (K) kz/ (K)2;

		///Evaluate the kernels in each cardinal direction.
		for(size_t i = 0; i < kx; i++)
			gaussx[i] = gauss1d((K) i, midgaussx, stdx);

		for(size_t i = 0; i < kx; i++)
			gaussy[i] = gauss1d((K) i, midgaussy, stdy);

		for(size_t i = 0; i < kx; i++)
			gaussz[i] = gauss1d((K) i, midgaussz, stdz);

		cpu_sepconv3(out, in, gaussx, gaussy, gaussz, dimx, dimy, dimz, kx, ky, kz);

	}
}
#endif