#ifndef BILATERALFILTERKERNEL_H__
#define BILATERALFILTERKERNEL_H__

#include <cuda_runtime.h>
#include <driver_types.h>
#include <channel_descriptor.h>
#include <cuda_runtime_api.h>

#ifndef M_PI
	#define M_PI 3.1415f
#endif
#define M_PI2 6.283f
#define VN 2.5066f


// Texture that holds the input image
//----------------------------------------------------------------------------
texture<float, 2, cudaReadModeElementType> InputImageTexture;

// Texture that holds the procomputed gaussian
texture<float, 2, cudaReadModeElementType> PrecomputedGaussianTexture;


// Division. If division remainder is neq zero then the result is ceiled
//----------------------------------------------------------------------------
#define DivUp(a,b) ((a % b != 0) ? (a/b + 1) : (a/b))


#if 0
// Using textures only
//----------------------------------------------------------------------------
template<int KernelRadius, unsigned int BlockSizeX, unsigned int BlockSizeY>
__global__ void
BilateralFilterKernel(unsigned int NX, unsigned int NY, float* Output, int NormalizedOutputPitch, float SigmaD, float SigmaD2, float SigmaR)
{
	// 2D index and linear index within this thread block
	int tu = threadIdx.x;
	int tv = threadIdx.y;
	
	// Global 2D index and linear index.
	int gu = blockIdx.x*BlockSizeX+tu;
	int gv = blockIdx.y*BlockSizeY+tv;
	int gl = gv*NormalizedOutputPitch + gu;

	// These threads do not process a pixel
	if ( gl >= (NY-1)*NormalizedOutputPitch+NX ) // if ( gu >= NX || gv >= NY )
		return;

	// Bilateral filtering for the current thread
	float CenterPixel = tex2D(InputImageTexture,gu,gv);
	float accu = 0.f;
	float norm = 0.f;
	float n1 = 1.f / (M_PI2*SigmaD2);
	float n2 = 1.f/(VN*SigmaR);
	float invSigmaR = 1.f/SigmaR;
	float invSigmaD2Half = -0.5f/SigmaD2;
	for (  int v = gv-KernelRadius; v <= gv+KernelRadius; v++ )
	{
		for ( int u = gu-KernelRadius; u <= gu+KernelRadius; u++ )
		{
			// 2D-Gaussian operating in the distance domain
			float f2 = v-gv;
			float f3 = u-gu;
			float f4 = sqrtf(f2*f2+f3*f3);
			float f5 = expf(f4*invSigmaD2Half);
			float w_d = n1*f5;

			//float w_d = tex2D(PrecomputedGaussianTexture,u-gu+KernelRadius,v-gv+KernelRadius);

			// 1D-Gaussian operating in the intensity range domain
			float sample = tex2D(InputImageTexture, u,v);
			 f2 = (CenterPixel-sample)*invSigmaR;
			 f3 = f2*f2;
			 f4 = expf(-0.5f*f3);
			float w_r = n2*f4;
				
			// Accumulate
			accu += w_d*w_r*sample;
			norm += w_d*w_r;
		}
	}

	// Write the result
	Output[gl] = accu/norm;
}


#else
// Using shared memory
//----------------------------------------------------------------------------
template<int KernelRadius, unsigned int BlockSizeX, unsigned int BlockSizeY>
__global__ void
BilateralFilterKernel(unsigned int NX, unsigned int NY, float* Output, int NormalizedOutputPitch, float SigmaD, float SigmaD2, float SigmaR)
{
	// 2D index and linear index within this thread block
	int tu = threadIdx.x;
	int tv = threadIdx.y;
	
	// Global 2D index and linear index.
	int gu = blockIdx.x*BlockSizeX+tu;
	int gv = blockIdx.y*BlockSizeY+tv;
	int gl = gv*NormalizedOutputPitch + gu;


	// Note that we cannot use the following termination criteria as threads processing these
	// regions might have to pull data
	/*if ( gu >= NU || gv >= NV )
		return;*/

	// Shared memory. We're using matrix and NOT texture notation (i.e. row first).
	// This scheme avoids bank conflicts
	const int SharedMemSizeX = BlockSizeX+2*KernelRadius;
	const int SharedMemSizeY = BlockSizeY+2*KernelRadius;
	__shared__ float SharedMem[SharedMemSizeY][SharedMemSizeX];

	// Successively pull data. First along the v-direction
	int sv = tv; // Index in shared memory
	int v = gv-KernelRadius; // Index in global memory
	while ( sv < SharedMemSizeY )
	{
		// Loop over the u-direction
		int su = tu; // Index in shared memory
		int u = gu-KernelRadius; // Index in global memory
		while ( su < SharedMemSizeX )
		{
			// Pull data for this thread
			SharedMem[sv][su] = tex2D(InputImageTexture,u,v);

			// Update u-coordinate in shared and global memory
			su += BlockSizeX;
			u += BlockSizeX;
		}

		// Update v-coordinate in shared memory
		sv += BlockSizeY;
		v += BlockSizeY;
	}

	// Wait until all threads have pulled their data
	__syncthreads();
	
	// These threads do not process a pixel
	if ( gl >= (NY-1)*NormalizedOutputPitch+NX ) // if ( gu >= NX || gv >= NY )
		return;

	// Bilateral filtering for the current thread. Note the radius offset as we had to pull additional data!
	// This means that threadIdx.x actually processes element SM[threadIdx.x+KernelRadius]
	float CenterPixel = SharedMem[tv+KernelRadius][tu+KernelRadius];//S(tv,tu);
	float accu = 0.f;
	float norm = 0.f;
	//float n1 = 1.f / (M_PI2*SigmaD2);
	float n2 = 1.f /(VN*SigmaR);
	float invSigmaR = 1.f/SigmaR;
	//float invSigmaD2Half = -0.5f/SigmaD2;
	for (  int v = tv-KernelRadius; v <= tv+KernelRadius; v++ )
	{
		for ( int u = tu-KernelRadius; u <= tu+KernelRadius; u++ )
		{
			// 2D-Gaussian operating in the distance domain
			/*float f2 = v-tv;
			float f3 = u-tu;
			float f4 = sqrt(f2*f2+f3*f3);
			float f5 = expf(f4*invSigmaD2Half);
			float w_d = n1*f5;*/

			float w_d = tex2D(PrecomputedGaussianTexture,u-tu+KernelRadius,v-tv+KernelRadius);


			// 1D-Gaussian operating in the intensity range domain
			float sample = SharedMem[v+KernelRadius][u+KernelRadius];
			float f2 = (CenterPixel-sample)*invSigmaR;
			float f3 = f2*f2; // sqrtf missing?
			float f4 = expf(-0.5f*f3);
			float w_r = n2*f4;
				
			// Accumulate
			float w = w_d*w_r;
			accu += w*sample;
			norm += w;
		}
	}

	// Write the result
	Output[gl] = accu/norm;
}
#endif


#endif // BILATERALFILTERKERNEL_H__
