#ifndef STIM_CUDA_GAUSS2_H
#define STIM_CUDA_GAUSS2_H

#include <stim/image/image.h>
#include <stim/math/filters/sepconv2.cuh>
#include <stim/math/constants.h>

namespace stim {

	template<typename T>
	T gauss1d(T x, T mu, T std) {
		return ((T)1 / (T)sqrt(stim::TAU * std * std) * exp(-(x - mu)*(x - mu) / (2 * std*std)));
	}
	///Perform a 2D gaussian convolution on an input image.
	///@param in is a pointer to the input image
	///@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>
	stim::image<T> cpu_gauss2(const stim::image<T>& in, K stdx, K stdy, size_t nstds = 3) {
		size_t kx = stdx * nstds * 2 + 1;					//calculate the kernel sizes
		size_t ky = stdy * nstds * 2 + 1;
		size_t X = in.width() - kx + 1;					//calculate the size of the output image
		size_t Y = in.height() - ky + 1;
		stim::image<T> r(X, Y, in.channels());		//create an output image

		K* gx = (K*)malloc(kx * sizeof(K));			//allocate space for the gaussian kernels
		K* gy = (K*)malloc(ky * sizeof(K));
		K mux = (K)kx / (K)2;						//calculate the mean of the gaussian (center of the kernel)
		K muy = (K)ky / (K)2;
		for (size_t xi = 0; xi < kx; xi++)			//evaluate the kernels
			gx[xi] = gauss1d((K)xi, mux, stdx);
		for (size_t yi = 0; yi < ky; yi++)
			gy[yi] = gauss1d((K)yi, muy, stdy);

		std::vector< stim::image<T> > clist = in.split();	//split the input image into channels
		std::vector< stim::image<T> > R = r.split();		//split the output image into channels
		for (size_t c = 0; c < in.channels(); c++)		//for each channel
			cpu_sepconv2(R[c].data(), clist[c].data(), gx, gy, clist[c].width(), clist[c].height(), kx, ky);

		r.merge(R);										//merge the blurred channels into the final image
		return r;
	}
}


#endif