#include <cuda.h> 
#include <math_functions.h>
#include <device_launch_parameters.h> 
#include <texture_fetch_functions.h> 
#include <builtin_types.h> 
#include <vector_functions.h> 
#include <float.h>

#define _SIZE_T_DEFINED 
#ifndef __CUDACC__ 
#define __CUDACC__ 
#endif 
#ifndef __cplusplus 
#define __cplusplus 
#endif

#define USE_TEXTURE 0
#define POWER_OF_TWO 1
#define CUDART_PI_F 3.141592654f

#if(USE_TEXTURE)
texture<float, 1, cudaReadModeElementType> texFloat;
#define   LOAD_FLOAT(i) tex1Dfetch(texFloat, i)
#define  SET_FLOAT_BASE checkCudaErrors( cudaBindTexture(0, texFloat, d_Src) )
#else
#define  LOAD_FLOAT(i) d_Src[i]
#define SET_FLOAT_BASE
#endif

#ifdef __CUDACC__
typedef float2 fComplex;
#else
typedef struct
{
    float x;
    float y;
} fComplex;
#endif

extern "C"
{	
	///////////////////////////////////////////////////////////////////////////////
	// Position convolution kernel center at (0, 0) in the image | VALIDATED
	///////////////////////////////////////////////////////////////////////////////
	__global__ void padKernel_kernel(float *d_Dst, float *d_Src, int fftH, int fftW, int kernelH, int kernelW, int kernelY, int kernelX)
	{
		const int y = blockDim.y * blockIdx.y + threadIdx.y;
		const int x = blockDim.x * blockIdx.x + threadIdx.x;

		if(y < kernelH && x < kernelW) // if we are in the kernel array
		{
			int ky = y - kernelY;
			if (ky < 0)
			{
				ky += fftH;
			}
			int kx = x - kernelX;
			if(kx < 0)
			{
				kx += fftW;
			}

			d_Dst[ky * fftW + kx] = LOAD_FLOAT(y * kernelW + x);
		}
	}

	////////////////////////////////////////////////////////////////////////////////
	// Prepare data for "pad to border" addressing mode | VALIDATED
	////////////////////////////////////////////////////////////////////////////////
	__global__ void padDataClampToBorder_kernel(float *d_Dst, float *d_Src, int fftH, int fftW, int dataH,	int dataW,	int kernelH, int kernelW,  int kernelY,	int kernelX)
	{
		const int y = blockDim.y * blockIdx.y + threadIdx.y;
		const int x = blockDim.x * blockIdx.x + threadIdx.x;
		const int borderH = dataH + kernelY;
		const int borderW = dataW + kernelX;

		if(y < fftH && x < fftW) // If we are in the array
		{
			int dy, dx;

			if (y < dataH) // if y is in the dataSize
			{
				dy = y;
			}

			if (x < dataW) // if x is in the dataSize
			{
				dx = x;
			}

			if ( y >= dataH && y < borderH) // if y is out of the dataSize but in the border
			{
				dy = dataH - 1;
			}

			if ( x >= dataW && x < borderW) // if x is out of the dataSize but in the border
			{
				dx = dataW - 1;
			}

			if(y >= borderH) // if y is out of the border
			{
				dy = 0;
			}

			if(x >= borderW) // if x is out of the border
			{
				dx = 0;
			}

			d_Dst[y * fftW + x] = LOAD_FLOAT(dy * dataW + dx);
		}
	}

	////////////////////////////////////////////////////////////////////////////////
	// Modulate Fourier image of padded data by Fourier image of padded kernel
	// and normalize by FFT size | VALIDATED
	////////////////////////////////////////////////////////////////////////////////
	__device__ void mulAndScale(fComplex &a, const fComplex &b, const float &c)
	{
		fComplex t = {c *(a.x * b.x - a.y * b.y), c *(a.y * b.x + a.x * b.y)};
		a = t;
	}

	__global__ void modulateAndNormalize_kernel( fComplex *d_Dst, fComplex *d_Src, int dataSize, float c)//, bool isX)
	{
		const int i = blockDim.x * blockIdx.x + threadIdx.x; // i (x, y)

		if (i >= dataSize)
		{
			return;
		}

		fComplex a = d_Src[i];
		fComplex b = d_Dst[i];

		// Temporal Convolution = PointWise Product
		mulAndScale(a, b, c);
		d_Dst[i] = a;
	}

	////////////////////////////////////////////////////////////////////////////////
	// Modulate Fourier kernel and derive it
	// and normalize by FFT size
	////////////////////////////////////////////////////////////////////////////////
	__global__ void derivate_kernel( fComplex *d_Dst, int dataSize)
	{
		const int i = blockDim.x * blockIdx.x + threadIdx.x;

		if (i >= dataSize)
		{
			return;
		}

		fComplex a = d_Dst[i];
		fComplex j = {0.0f, 2 * CUDART_PI_F * blockIdx.x / rsqrtf(i)};
		mulAndScale(a, j, 1.0f);

		d_Dst[i] = a;
	}

	///////////////////////////////////////////////////////////////////////////////
	// Position convolution kernel center at (0, 0) in the image | VALIDATED
	///////////////////////////////////////////////////////////////////////////////
	__global__ void resetMap_kernel(float *d_Dst, int dataSize)
	{
		const int i = blockDim.x * blockIdx.x + threadIdx.x;
		if (i >= dataSize)
		{
			return;
		}
		d_Dst[i] = 0;
	}

	// Unused
	__global__ void resetMap2D_kernel(float2 *d_Dst, int dataSize)
	{
		const int i = blockDim.x * blockIdx.x + threadIdx.x;
		if (i >= dataSize)
		{
			return;
		}
		d_Dst[i].x = 0.0f;
		d_Dst[i].y = 0.0f;
	}

	///////////////////////////////////////////////////////////////////////////////
	// Splatting line on maps | VALIDATED
	///////////////////////////////////////////////////////////////////////////////
	__global__ void splatAccMap_kernel(float *d_DstDirXX, float *d_DstDirXY, float *d_DstDirYX, float *d_DstDirYY, float2 *div, float2 *dirMap, uint4 *endLine, int dataH, int dataW, int accRes)
	{
		if (threadIdx.x == 0 || threadIdx.x >= endLine[blockIdx.x].x - 2)
		{
			return;
		}
		// d_Dst[y * datW + x]
		int index = (int)(div[blockDim.x * blockIdx.x + threadIdx.x].y * ((float)accRes - 1.0f)) * dataW + (int)(div[blockDim.x * blockIdx.x + threadIdx.x].x * ((float)accRes - 1.0f));
	/*	d_DstDirXX[index] += dirMap[blockDim.x * blockIdx.x + threadIdx.x].x;
		d_DstDirXY[index] += dirMap[blockDim.x * blockIdx.x + threadIdx.x].y;
		d_DstDirYX[index] += dirMap[blockDim.x * blockIdx.x + threadIdx.x].x;
		d_DstDirYY[index] += dirMap[blockDim.x * blockIdx.x + threadIdx.x].y;*/
		atomicAdd(&d_DstDirXX[index], dirMap[blockDim.x * blockIdx.x + threadIdx.x].x);
		atomicAdd(&d_DstDirXY[index], dirMap[blockDim.x * blockIdx.x + threadIdx.x].y);
		atomicAdd(&d_DstDirYX[index], dirMap[blockDim.x * blockIdx.x + threadIdx.x].x);
		atomicAdd(&d_DstDirYY[index], dirMap[blockDim.x * blockIdx.x + threadIdx.x].y);
	
	}

	// Unused
	__global__ void splatDirMap_kernel(float2 *d_Dst, float2 *div, float2 *dirMap, uint4 *endLine, int dataH, int dataW, int accRes)
	{
		if (threadIdx.x >= endLine[blockIdx.x].x)
		{
			return;
		}
		// d_Dst[y * datW + x]
		atomicAdd(&d_Dst[(int)(div[blockDim.x * blockIdx.x + threadIdx.x].y * ((float)accRes - 1.0f)) * dataW + (int)(div[blockDim.x * blockIdx.x + threadIdx.x].x * ((float)accRes - 1.0f))].x, dirMap[blockDim.x * blockIdx.x + threadIdx.x].x);
		atomicAdd(&d_Dst[(int)(div[blockDim.x * blockIdx.x + threadIdx.x].y * ((float)accRes - 1.0f)) * dataW + (int)(div[blockDim.x * blockIdx.x + threadIdx.x].x * ((float)accRes - 1.0f))].y, dirMap[blockDim.x * blockIdx.x + threadIdx.x].y);
	}

	///////////////////////////////////////////////////////////////////////////////
	// Advecting points kernel | 
	///////////////////////////////////////////////////////////////////////////////	
	__device__ void divise(float2 &a, const float &b)
	{
		float2 c = { a.x / b, a.y / b};
		a = c;
	}

	__global__ void advection_kernel(float2 *d_Dst, float2 *dirData, float* d_SrcDirXXMap, float* d_SrcDirXYMap, float* d_SrcDirYXMap, float* d_SrcDirYYMap, uint4 *endLine, int srcW, int accRes, float attractCoef, int maxLine, int halfkSize)
	{
		if(blockIdx.x >= maxLine)
		{
			return;
		}
		if(threadIdx.x >= endLine[blockIdx.x].x - 2 || threadIdx.x <= 1)
		{
			return;
		}
		
		extern __shared__ float2 DstS[];		
		DstS[threadIdx.x] = d_Dst[blockDim.x * blockIdx.x + threadIdx.x];
		__syncthreads();	
		float2 dir = dirData[blockDim.x * blockIdx.x + threadIdx.x];

		int x = (int)(DstS[threadIdx.x].x  * ((float)accRes - 1.0f));
		int y = (int)(DstS[threadIdx.x].y  * ((float)accRes - 1.0f));

		float2 gradToDirX = {d_SrcDirXXMap[y * srcW + x], d_SrcDirYXMap[y * srcW + x]}; // Take this into account depending on dir.x value (from 0 to 1)
		float2 gradToDirY = {d_SrcDirXYMap[y * srcW + x], d_SrcDirYYMap[y * srcW + x]}; // Take this into account depending on dir.y value (from 0 to 1)
		//// WITH DIRECTION
		float2 finalGradToDir = {gradToDirX.x * dir.x + gradToDirY.x * dir.y, gradToDirX.y * dir.x + gradToDirY.y * dir.y}; // Here using scalar product but all bilinear methods works fine
		// NO DIRECTION
		//float2 finalGradToDir = {gradToDirX.x + gradToDirY.x, gradToDirX.y + gradToDirY.y}; // Here using scalar product but all bilinear methods works fine
		float gradNorm = norm3df(finalGradToDir.x, finalGradToDir.y, 0.0);
		if(gradNorm < 0.00001 || isnan(gradNorm))
		{
			gradNorm = 1.0;
		}
		finalGradToDir.x = finalGradToDir.x / gradNorm;
		finalGradToDir.y = finalGradToDir.y / gradNorm;

		float2 actualNormal = {-(DstS[threadIdx.x + 1].y - DstS[threadIdx.x - 1].y), DstS[threadIdx.x + 1].x - DstS[threadIdx.x - 1].x};
		float normNormalVec = norm3df(actualNormal.x, actualNormal.y, 0.0f);
		if(normNormalVec < 0.00001 || isnan(normNormalVec))
		{
			normNormalVec = 1.0;
		}
		actualNormal.x = actualNormal.x / normNormalVec;
		actualNormal.y = actualNormal.y / normNormalVec;

		float2 finalDelta = {(finalGradToDir.x * actualNormal.x + finalGradToDir.y * actualNormal.y) * actualNormal.x, (finalGradToDir.x * actualNormal.x + finalGradToDir.y * actualNormal.y) * actualNormal.y}; //finalGradToDir;//

		DstS[threadIdx.x].x = max(min(DstS[threadIdx.x].x + attractCoef * finalDelta.x, 1.0), 0.0);
		DstS[threadIdx.x].y = max(min(DstS[threadIdx.x].y + attractCoef * finalDelta.y, 1.0), 0.0);// ;

		__syncthreads();

		d_Dst[blockDim.x * blockIdx.x + threadIdx.x] = DstS[threadIdx.x];
	}

	///////////////////////////////////////////////////////////////////////////////
	// Resampling Lines | VALIDATED but long
	///////////////////////////////////////////////////////////////////////////////	
	__device__ void diff(float2 &a, const float2 &b)
	{
		float2 c = { a.x - b.x, a.y - b.y};
		a = c;
	}
	
	// Possible optimisation with sharedMemoryUsage (put a line in shared mem)
	__global__ void resampling_kernel(float2 *d_Dst, int *endL_Dst, float2 *d_Src, uint4 *endL_Src, int lineSize, int maxLine, float splitDist, float remDist)
	{
		const int i = blockDim.x * blockIdx.x + threadIdx.x;
		if (i >= maxLine)
		{
			return;
		}

		d_Dst[i * lineSize] = d_Src[i * lineSize];
		int newIndex = 2;
		d_Dst[i * lineSize + 1] = d_Src[i * lineSize + 1]; 
		for (int j = 1; j < endL_Src[i].x - 1; j++)
		{
			float2 pos = d_Src[i * lineSize + j];
			float2 dist1 = pos;
			diff(dist1, d_Src[i * lineSize + j + 1]);
			float dist = norm3df(dist1.x, dist1.y, 0.0f);
			if (dist > splitDist)
			{
				divise(dist1, 2.0f);
				diff(pos, dist1);
				d_Dst[i * lineSize + newIndex] = pos;
				newIndex++;
			}
			if (dist > remDist)
			{
				d_Dst[i * lineSize + newIndex] = d_Src[i * lineSize + j + 1];
				newIndex++;
			}
			if (j == endL_Src[i].x - 2)
			{
				d_Dst[i * lineSize + newIndex] = d_Src[i * lineSize + j + 1];
				d_Dst[i * lineSize + newIndex + 1] = d_Src[i * lineSize + j + 1];
				newIndex += 2;
			}
		}
		endL_Dst[i] = newIndex;
		//drawArray[i].x = newIndex;
	}
	
	///////////////////////////////////////////////////////////////////////////////
	// Smoothing Points
	///////////////////////////////////////////////////////////////////////////////
	__device__ void laplacian(float2 &a, const float2 &b, const float2 &d, const float &s)
	{
		a.x = s * b.x + (1.0f - 2.0f * s) * a.x + s * d.x;
		a.y = s * b.y + (1.0f - 2.0f * s) * a.y + s * d.y;
	}

	// Possible optimisation with sharedMemoryUsage (put a line in shared mem)
	__global__ void smoothing_kernel(float2 *d_Dst, uint4 *endL_Dst, float2 *d_Src, int *endL_Src, int maxLine, float smoothFactor)
	{
		if(blockIdx.x >= maxLine || threadIdx.x >= endL_Src[blockIdx.x])
		{
			return;
		}
		//int i = blockDim.x * blockIdx.x + threadIdx.x;

		extern __shared__ float2 temp[];
		temp[threadIdx.x] = d_Src[blockDim.x * blockIdx.x + threadIdx.x];
		__syncthreads();

		//float2 res = d_Src[i];
		if(threadIdx.x > 1 && threadIdx.x < endL_Src[blockIdx.x] - 2)
		{
			laplacian(temp[threadIdx.x], temp[threadIdx.x - 1], temp[threadIdx.x + 1], smoothFactor);
		}
		d_Dst[blockDim.x * blockIdx.x + threadIdx.x] = temp[threadIdx.x];
		endL_Dst[blockIdx.x].x = endL_Src[blockIdx.x];
	}

	///////////////////////////////////////////////////////////////////////////////
	// Compute Line Width
	///////////////////////////////////////////////////////////////////////////////
	__global__ void lineWidth_kernel(float *d_DstDirX, float *d_DstDirY, float2 *d_Dst,int maxLine, float2 *dirMap, uint4 *endL_Dst, float *width_Dst, float compatibilityAngle, int srcW, int accRes)
	{
		if(blockIdx.x >= maxLine || threadIdx.x == 0 || threadIdx.x >= endL_Dst[blockIdx.x].x - 1)
		{
			return;
		}

		float2 origDir = dirMap[blockDim.x * blockIdx.x + threadIdx.x];
		float normOrigDir = norm3df(origDir.x, origDir.y, 0.0f);

		int x = (int)(d_Dst[blockDim.x * blockIdx.x + threadIdx.x].x  * ((float)accRes - 1.0f));
		int y = (int)(d_Dst[blockDim.x * blockIdx.x + threadIdx.x].y  * ((float)accRes - 1.0f));
		float2 accDir = {d_DstDirX[y * srcW + x], d_DstDirY[y * srcW + x]};

		float normAccDir = norm3df(accDir.x, accDir.y, 0.0f);
		float scalProd = accDir.x / normAccDir * origDir.x / normOrigDir + accDir.y / normAccDir * origDir.y / normOrigDir;

		float tempWidth = 0.0;
		if(abs(scalProd) >= compatibilityAngle)
		{
			tempWidth = normAccDir;
		} 

		width_Dst[blockDim.x * blockIdx.x + threadIdx.x] = tempWidth;
	}

	__global__ void normalizeLineWidth_kernel(float *width_Dst, float maxWidth, float minWidth, float ratio, int maxLine)
	{
		if(blockIdx.x >= maxLine)
		{
			return;
		}
		width_Dst[blockDim.x * blockIdx.x + threadIdx.x] = (ratio - 1.0f) * (width_Dst[blockDim.x * blockIdx.x + threadIdx.x] - minWidth) / (maxWidth - minWidth) + 1.0f;
	}
}