#ifndef BILATERALFILTERKERNEL_H__
#define BILATERALFILTERKERNEL_H__

#define PixelType float
#define M_PI 3.1415

__global__ void
d_BilateralFilter(PixelType* id, PixelType* od, int w, int h, int KernelSize, float SigmaD, float SigmaR)
{
	int x = blockIdx.x*blockDim.x + threadIdx.x;
	if ( x >= w*h )
		return;

	// Position
	int Row = (int)(x/w);
	int Column = (int)(x%w);

	// Boundary check
	if ( Row-KernelSize < 0 || Row+KernelSize >= h || Column-KernelSize < 0 || Column+KernelSize >= w )
	{
		od[x] = id[x];
		return;
	}
	
	// Convolution
	float accu = 0;
	float norm = 0;
	for ( int i = Row-KernelSize; i <= Row+KernelSize; i++ )
	{
		for ( int j = Column-KernelSize; j <= Column+KernelSize; j++ )
		{
			// 2D-Gaussian operating in the distance domain
			float f1 = 1.f / 2*M_PI*SigmaD*SigmaD;
			float f2 = i-Row;
			float f3 = j-Column;
			float f4 = f2*f2+f3*f3;
			float f5 = expf(-0.5*f4/(SigmaD*SigmaD));
			float w_d = f1*f5;

			// 1D-Gaussian operating in the intensity range domain
			float g1 = 1.f/(sqrtf(2*M_PI)*SigmaR);
			float g2 = (id[Row*w+Column] - id[i*w+j])/SigmaR;
			float g3 = g2*g2;
			float g4 = exp(-0.5*g3);
			float w_r = g1*g4;
			
			// Accumulate
			accu += w_d*w_r*id[i*w+j];
			norm += w_d*w_r;
		}
	}
	od[x] = (PixelType)(accu/norm);
}

#endif // BILATERALFILTERKERNEL_H__