/**********************************************************************************************//**
 * @file	CUDA_Func_Library.cu
 *
 * @brief	Cuda function library class.
 **************************************************************************************************/


#include "CUDA_Func_Library.cuh"

//
// CUDA device info functions
//

/**********************************************************************************************//**
 * @fn	int cudaInfo::cudaDevCount()
 *
 * @brief	Function for getting number of CUDA devices available in system.
 *
 * @author	Tomas Bravenec
 *
 * @return	Number of CUDA devices.
 **************************************************************************************************/
int cudaInfo::cudaDevCount()
{
	int devC;
	if (cudaGetDeviceCount(&devC) != cudaSuccess)
	{
		devC = 0;
	};
	return devC;
}

/**********************************************************************************************//**
 * @fn	int cudaInfo::cudaGetCUDACores(int device)
 *
 * @brief	Function for calculating CUDA core count depending on GPU architecture.
 *
 * @author	Tomas Bravenec
 *
 * @param	device Index of GPU.
 *
 * @return	CUDA core count.
 **************************************************************************************************/
int cudaInfo::cudaGetCUDACores(int device)
{
	int coreCount;

	cudaDeviceProp prop;
	cudaGetDeviceProperties(&prop, device);

	switch (prop.major)
	{
	case 2: // Fermi
		if (prop.minor == 1)
		{
			coreCount = prop.multiProcessorCount * 48;
		}
		else
		{
			coreCount = prop.multiProcessorCount * 32;
		}
		break;
	case 3: // Kepler
		coreCount = prop.multiProcessorCount * 192;
		break;
	case 5: // Maxwell
		coreCount = prop.multiProcessorCount * 128;
		break;
	case 6: // Pascal
		if (prop.minor == 1)
		{
			coreCount = prop.multiProcessorCount * 128;
		}
		else
		{
			coreCount = prop.multiProcessorCount * 64;
		}
		break;
	default:
		coreCount = 0;
		break;
	}
	return coreCount;
}

/**********************************************************************************************//**
 * @fn	char* cudaInfo::cudaGetError(cudaError_t cudaStatus)
 *
 * @brief	Gets readable explanation of cudaError code.
 *
 * @author	Tomas Bravenec
 *
 * @param	cudaStatus	CUDA error code.
 *
 * @return	Pointer to a char array containing meaning of error code.
 **************************************************************************************************/
char* cudaInfo::cudaGetError(cudaError_t cudaStatus)
{
	switch (cudaStatus)
	{
	case cudaSuccess:
		return "No errors";

	case cudaErrorMissingConfiguration:
		return "Missing configuration error";

	case cudaErrorMemoryAllocation:
		return "Memory allocation error";

	case cudaErrorInitializationError:
		return "Initialization error";

	case cudaErrorLaunchFailure:
		return "Launch failure";

	case cudaErrorPriorLaunchFailure:
		return "Prior launch failure";

	case cudaErrorLaunchTimeout:
		return "Launch timeout error";

	case cudaErrorLaunchOutOfResources:
		return "Launch out of resources error";

	case cudaErrorInvalidDeviceFunction:
		return "Invalid device function";

	case cudaErrorInvalidConfiguration:
		return "Invalid configuration";

	case cudaErrorInvalidDevice:
		return "Invalid device";

	case cudaErrorInvalidValue:
		return "Invalid value";

	case cudaErrorInvalidPitchValue:
		return "Invalid pitch value";

	case cudaErrorInvalidSymbol:
		return "Invalid symbol";

	case cudaErrorMapBufferObjectFailed:
		return "Map buffer object failed";

	case cudaErrorUnmapBufferObjectFailed:
		return "Unmap buffer object failed";

	case cudaErrorInvalidHostPointer:
		return "Invalid host pointer";

	case cudaErrorInvalidDevicePointer:
		return "Invalid device pointer";

	case cudaErrorInvalidTexture:
		return "Invalid texture";

	case cudaErrorInvalidTextureBinding:
		return "Invalid texture binding";

	case cudaErrorInvalidChannelDescriptor:
		return "Invalid channel descriptor";

	case cudaErrorInvalidMemcpyDirection:
		return "Invalid memcpy direction";

	case cudaErrorAddressOfConstant:
		return "Address of constant error";

	case cudaErrorTextureFetchFailed:
		return "Texture fetch failed";

	case cudaErrorTextureNotBound:
		return "Texture not bound error";

	case cudaErrorSynchronizationError:
		return "Synchronization error";

	case cudaErrorInvalidFilterSetting:
		return "Invalid filter setting";

	case cudaErrorInvalidNormSetting:
		return "Invalid norm setting";

	case cudaErrorMixedDeviceExecution:
		return "Mixed device execution";

	case cudaErrorCudartUnloading:
		return "CUDA runtime unloading";

	case cudaErrorUnknown:
		return "Unknown error condition";

	case cudaErrorNotYetImplemented:
		return "Function not yet implemented";

	case cudaErrorMemoryValueTooLarge:
		return "Memory value too large";

	case cudaErrorInvalidResourceHandle:
		return "Invalid resource handle";

	case cudaErrorNotReady:
		return "Not ready error";

	case cudaErrorInsufficientDriver:
		return "CUDA runtime is newer than driver";

	case cudaErrorSetOnActiveProcess:
		return "Set on active process error";

	case cudaErrorNoDevice:
		return "No available CUDA device";
	
	case cudaErrorStartupFailure:
		return "Startup failure";

	case cudaErrorApiFailureBase:
		return "API failure base";
	}

	return "<unknown>";
}


/**********************************************************************************************//**
 * 
 * @fn	void cudaCluster::randomizeCentroid()
 * 		
 * @brief	Randomizes centroid.
 *
 * @author	Tomas Bravenec
 **************************************************************************************************/
void cudaCluster::randomizeCentroid()
{
	centroidBlue = rand();
	centroidGreen = rand();
	centroidRed = rand();
}

/**********************************************************************************************//**
 * @fn	void cudaCluster::setCentroid(BYTE Blue, BYTE Green, BYTE Red)
 *
 * @brief	Sets a centroid to manual color.
 *
 * @author	Tomas Bravenec
 *
 * @param	Blue 	The blue.
 * @param	Green	The green.
 * @param	Red  	The red.
 **************************************************************************************************/
void cudaCluster::setCentroid(BYTE Blue, BYTE Green, BYTE Red)
{
	centroidRed = Red;
	centroidGreen = Green;
	centroidBlue = Blue;
}

//
// Helper functions allowing to launch CUDA kernels on GPU
//
/**********************************************************************************************//**
 * @fn	cudaError_t cudaKernelLauncher::cudaRGB2Gray(BYTE *rgb, BYTE *gray, int height, int width, int stride, int threads)
 *
 * @brief	Launcher of CUDA kernel for conversion of RGB image to grayscale image.
 *
 * @author	Tomas Bravenec
 *
 * @param [in]		rgb	   	Pointer to first byte of image in 24bpp BGR format.
 * @param [out]		gray   	Pointer to first byte of grayscale image.
 * @param 		  	height 	Height of the image in pixels.
 * @param 		  	width  	Width of the image in pixels.
 * @param 		  	stride 	Length of one row in bytes, including padding.
 * @param 		  	threads	Square root of maximum allowed threads per block.
 *
 * @return A cudaError_t containing cudaError code, cudaSucces if no error.
 **************************************************************************************************/
cudaError_t cudaKernelLauncher::cudaRGB2Gray(BYTE *rgb, BYTE *gray, int height, int width, int stride, int threads)
{
	cudaError_t cudaStatus = cudaSuccess;

	dim3 dimBlock(threads, threads); 
	dim3 dimGrid(width / dimBlock.x + 1, height / dimBlock.y + 1);

	rgbToGray_CUDA<<< dimGrid, dimBlock >>>(rgb, gray, height, width, stride);

	// get last CUDA error and return it
	cudaStatus = cudaGetLastError();
	return cudaStatus;
}

/**********************************************************************************************//**
 * @fn	cudaError_t cudaKernelLauncher::cudaGray2RGB(BYTE *rgb, BYTE *gray, int height, int width, int stride, int threads)
 *
 * @brief	Launcher of CUDA kernel for conversion of grayscale image to RGB image format.
 *
 * @author	Tomas Bravenec
 *
 * @param [out]		rgb	   	Pointer to first byte of image in 24bpp BGR format.
 * @param [in]		gray   	Pointer to first byte of grayscale image.
 * @param 		  	height 	Height of the image in pixels.
 * @param 		  	width  	Width of the image in pixels.
 * @param 		  	stride 	Length of one row in bytes, including padding.
 * @param 		  	threads	Square root of maximum allowed threads per block.
 *
 * @return A cudaError_t containing cudaError code, cudaSucces if no error..
 **************************************************************************************************/
cudaError_t cudaKernelLauncher::cudaGray2RGB(BYTE *rgb, BYTE *gray, int height, int width, int stride, int threads)
{
	cudaError_t cudaStatus;

	dim3 dimBlock(threads, threads);
	dim3 dimGrid(width / dimBlock.x + 1, height / dimBlock.y + 1);

	grayToRGB_CUDA<<< dimGrid, dimBlock >>>(rgb, gray, height, width, stride);

	// get last CUDA error and return it
	cudaStatus = cudaGetLastError();
	return cudaStatus;
}

/**********************************************************************************************//**
 * @fn	cudaError_t cudaKernelLauncher::cudaEdgeDetection(BYTE *gray, BYTE *output, int height, int width, int type, int threads)
 *
 * @brief	Launcher of CUDA kernel for edge detection in 8bpp grayscale image.
 *
 * @author	Tomas Bravenec
 *
 * @param [in]		gray   	Pointer to first byte of grayscale image.
 * @param [out]		output 	If non-null, the output.
 * @param 		  	height 	Height of the image in pixels.
 * @param 		  	width  	Width of the image in pixels.
 * @param 		  	type   	The type.
 * @param 		  	threads	Square root of maximum allowed threads per block.
 *
 * @return A cudaError_t containing cudaError code, cudaSucces if no error..
 **************************************************************************************************/
cudaError_t cudaKernelLauncher::cudaEdgeDetection(BYTE *gray, BYTE *output, int height, int width, int type, int threads)
{
	cudaError_t cudaStatus = cudaSuccess;
	int *d_filterX, *d_filterY;
	int	filterX[9], filterY[9];

	dim3 dimBlock(threads, threads);
	dim3 dimGrid(width / dimBlock.x + 1, height / dimBlock.y + 1);
	
	// filter declaration: default equals to case 0 - Sobel filter, 1 - Prewitt, 2 - Kirsch
	switch (type)
	{
	case 0:		// Sobel filter
		filterX[0] = -1; filterX[1] = 0; filterX[2] = 1; filterX[3] = -2; filterX[4] = 0; filterX[5] = 2; filterX[6] = -1; filterX[7] = 0; filterX[8] = 1; 
		filterY[0] = 1; filterY[1] = 2; filterY[2] = 1; filterY[3] = 0; filterY[4] = 0; filterY[5] = 0; filterY[6] = -1; filterY[7] = -2; filterY[8] = -1;
		break;
	case 1:		// Prewitt filter
		filterX[0] = -1; filterX[1] = 0; filterX[2] = 1; filterX[3] = -1; filterX[4] = 0; filterX[5] = 1; filterX[6] = -1; filterX[7] = 0; filterX[8] = 1;
		filterY[0] = 1; filterY[1] = 1; filterY[2] = 1; filterY[3] = 0; filterY[4] = 0; filterY[5] = 0; filterY[6] = -1; filterY[7] = -1; filterY[8] = -1;
		break;
	case 2:		// Kirsch
		filterX[0] = -5; filterX[1] = 3; filterX[2] = 3; filterX[3] = -5; filterX[4] = 0; filterX[5] = 3; filterX[6] = -5; filterX[7] = 3; filterX[8] = 3;
		filterY[0] = 3; filterY[1] = 3; filterY[2] = 3; filterY[3] = 3; filterY[4] = 0; filterY[5] = 3; filterY[6] = -5; filterY[7] = -5; filterY[8] = -5;
		break;
	default:	// Sobel filter
		filterX[0] = -1; filterX[1] = 0; filterX[2] = 1; filterX[3] = -2; filterX[4] = 0; filterX[5] = 2; filterX[6] = -1; filterX[7] = 0; filterX[8] = 1;
		filterY[0] = 1; filterY[1] = 2; filterY[2] = 1; filterY[3] = 0; filterY[4] = 0; filterY[5] = 0; filterY[6] = -1; filterY[7] = -2; filterY[8] = -1;
		break;
	}

	// allocate memory on GPU for filters
	cudaStatus = cudaMalloc(&d_filterX, 9 * sizeof(int));
	if (cudaStatus != cudaSuccess)
	{
		goto cudaCancel;
	}
	cudaStatus = cudaMalloc(&d_filterY, 9 * sizeof(int));
	if (cudaStatus != cudaSuccess)
	{
		goto cudaCancel;
	}

	// copy filters to GPU memory
	cudaStatus = cudaMemcpy(d_filterX, &filterX, 9 * sizeof(int), cudaMemcpyHostToDevice);
	if (cudaStatus != cudaSuccess)
	{
		goto cudaCancel;
	}
	cudaStatus = cudaMemcpy(d_filterY, &filterY, 9 * sizeof(int), cudaMemcpyHostToDevice);
	if (cudaStatus != cudaSuccess)
	{
		goto cudaCancel;
	}

	// run kernel
	edge_CUDA<<< dimGrid, dimBlock >>>(gray, output, height, width, d_filterX, d_filterY);

	// get last CUDA error
	cudaStatus = cudaGetLastError();

cudaCancel:
	// free memory used by edge detection
	cudaFree(d_filterX);
	cudaFree(d_filterY);
	// return last cuda Error
	return cudaStatus;
}

/**********************************************************************************************//**
 * @fn	cudaError_t cudaKernelLauncher::cudaThresholdEdge(BYTE *image, int height, int width, double threshold, int threads)
 *
 * @brief	Launcher of tresholding image after edge detection.
 *
 * @author	Tomas Bravenec
 *
 * @param [in,out]	image	 	Pointer to first byte of grayscale image.
 * @param 		  	height   	Height of the image in pixels.
 * @param 		  	width	 	Width of the image in pixels.
 * @param 		  	threshold	The threshold.
 * @param 		  	threads  	Square root of maximum allowed threads per block.
 *
 * @return A cudaError_t containing cudaError code, cudaSucces if no error..
 **************************************************************************************************/
cudaError_t cudaKernelLauncher::cudaThresholdEdge(BYTE *image, int height, int width, double threshold, int threads)
{
	cudaError_t cudaStatus = cudaSuccess;
	
	BYTE *d_preMax, *preMax;
	BYTE Max = 0;

	// dimensions and number of blocks for maximum kernel
	dim3 dimBlockCol(threads * threads);
	dim3 dimGridCol(width / dimBlockCol.x + 1);

	// dimensions and number of blocks for maximum kernel
	dim3 dimBlockRow(1, threads * threads);
	dim3 dimGridRow(1, height / dimBlockRow.y + 1);

	// dimensions and number of blocks for threasholding kernel
	dim3 dimBlockThreshold(threads, threads);
	dim3 dimGridThreshold(width / dimBlockThreshold.x + 1, height / dimBlockThreshold.y + 1);

	// choose more optimal find maximum, if there is more colums, choose find maximum in columns
	if (width >= height)
	{
		// allocate system memory for parallel maximum output
		preMax = (BYTE*)malloc(width * sizeof(BYTE));

		// allocate GPU memory for parallel maximum output
		cudaStatus = cudaMalloc(&d_preMax, width * sizeof(BYTE));
		if (cudaStatus != cudaSuccess)
		{
			goto cudaCancel;
		}

		// run kernel
		maxCol_CUDA<<< dimGridCol, dimBlockCol >>>(image, d_preMax, height, width);
		cudaStatus = cudaGetLastError();
		if (cudaStatus != cudaSuccess)
		{
			goto cudaCancel;
		}
		// copy maximum of all columns to system memory
		cudaStatus = cudaMemcpy(preMax, d_preMax, width * sizeof(BYTE), cudaMemcpyDeviceToHost);
		if (cudaStatus != cudaSuccess)
		{
			goto cudaCancel;
		}

		// find final maximum on fast CPU core
		for (int i = 0; i < width; i++)
		{
			if (preMax[i] > Max)
			{
				Max = preMax[i];
			}
		}
	}
	else
	{
		// allocate system memory for parallel maximum output
		preMax = (BYTE*)malloc(height * sizeof(BYTE));

		// allocate GPU memory for parallel maximum output
		cudaStatus = cudaMalloc(&d_preMax, height * sizeof(BYTE));
		if (cudaStatus != cudaSuccess)
		{
			goto cudaCancel;
		}

		// run kernel
		maxRow_CUDA<<< dimGridRow, dimBlockRow >>>(image, d_preMax, height, width);
		cudaStatus = cudaGetLastError();
		if (cudaStatus != cudaSuccess)
		{
			goto cudaCancel;
		}

		// copy maximum of all columns to system memory
		cudaStatus = cudaMemcpy(preMax, d_preMax, height * sizeof(BYTE), cudaMemcpyDeviceToHost);
		if (cudaStatus != cudaSuccess)
		{
			goto cudaCancel;
		}

		// find final maximum on fast CPU core
		for (int i = 0; i < height; i++)
		{
			if (preMax[i] > Max)
			{
				Max = preMax[i];
			}
		}
	}
	// run kernel
	simpleThresholdEdge_CUDA << <dimGridThreshold, dimBlockThreshold >> >(image, height, width, Max, threshold);
	cudaStatus = cudaGetLastError();

cudaCancel:
	// free memory used by find maximum
	free(preMax);
	cudaFree(d_preMax);
	return cudaStatus;
}

/**********************************************************************************************//**
 * @fn	cudaError_t cudaKernelLauncher::cudaThreshold(BYTE *image, int height, int width, double top, double bottom, int type, int threads)
 *
 * @brief	Launcher of CUDA kernels for tresholding.
 *
 * @author	Tomas Bravenec
 *
 * @param [in,out]	image  	Pointer to first byte of grayscale image.
 * @param 		  	height 	Height of the image in pixels.
 * @param 		  	width  	Width of the image in pixels.
 * @param 		  	top	   	Threshold for simple tresholding or upper treshold for dual tresholding.
 * @param 		  	bottom 	Lower threshold for dual tresholding, set 0 for simple tresholding.
 * @param 		  	type   	Type of tresholding.
 * @param 		  	threads	Square root of maximum allowed threads per block.
 *
 * @return A cudaError_t containing cudaError code, cudaSucces if no error..
 **************************************************************************************************/
cudaError_t cudaKernelLauncher::cudaThreshold(BYTE *image, int height, int width, double top, double bottom, int type, int threads)
{
	cudaError_t cudaStatus = cudaSuccess;

	BYTE *d_preMax, *preMax;
	BYTE Max = 0;

	// dimensions and number of blocks for maximum kernel
	dim3 dimBlockCol(threads * threads);
	dim3 dimGridCol(width / dimBlockCol.x + 1);

	// dimensions and number of blocks for maximum kernel
	dim3 dimBlockRow(1, threads * threads);
	dim3 dimGridRow(1, height / dimBlockRow.y + 1);

	// dimensions and number of blocks for threasholding kernel
	dim3 dimBlockThreshold(threads, threads);
	dim3 dimGridThreshold(width / dimBlockThreshold.x + 1, height / dimBlockThreshold.y + 1);

	// choose more optimal find maximum, if there is more colums, choose find maximum in columns
	if (width >= height)
	{
		// allocate system memory for parallel maximum output
		preMax = (BYTE*)malloc(width * sizeof(BYTE));

		// allocate GPU memory for parallel maximum output
		cudaStatus = cudaMalloc(&d_preMax, width * sizeof(BYTE));
		if (cudaStatus != cudaSuccess)
		{
			goto cudaCancel;
		}

		// run kernel
		maxCol_CUDA<<< dimGridCol, dimBlockCol>>>(image, d_preMax, height, width);
		cudaStatus = cudaGetLastError();
		if (cudaStatus != cudaSuccess)
		{
			goto cudaCancel;
		}
	}
	else
	{
		// allocate system memory for parallel maximum output
		preMax = (BYTE*)malloc(height * sizeof(BYTE));

		// allocate GPU memory for parallel maximum output
		cudaStatus = cudaMalloc(&d_preMax, height * sizeof(BYTE));
		if (cudaStatus != cudaSuccess)
		{
			goto cudaCancel;
		}

		// run kernel
		maxRow_CUDA<<< dimGridRow, dimBlockRow>>>(image, d_preMax, height, width);
		cudaStatus = cudaGetLastError();
		if (cudaStatus != cudaSuccess)
		{
			goto cudaCancel;
		}
	}

	// copy maximum of all columns to system memory
	cudaStatus = cudaMemcpy(preMax, d_preMax, width * sizeof(BYTE), cudaMemcpyDeviceToHost);
	if (cudaStatus != cudaSuccess)
	{
		goto cudaCancel;
	}

	// find final maximum on fast CPU core
	for (int i = 0; i < width; i++)
	{
		if (preMax[i] > Max)
		{
			Max = preMax[i];
		}
	}

	// select and run user selected type of tresholding, return cudaStatus in case of error
	switch (type)
	{
	case 0:
		simpleThreshold_CUDA<<<dimGridThreshold, dimBlockThreshold>>>(image, height, width, Max, top);
		cudaStatus = cudaGetLastError();
		if (cudaStatus != cudaSuccess)
		{
			goto cudaCancel;
		}
		break;
	case 1:
		binaryThreshold_CUDA<<<dimGridThreshold, dimBlockThreshold>>>(image, height, width, Max, bottom, top);
		cudaStatus = cudaGetLastError();
		if (cudaStatus != cudaSuccess)
		{
			goto cudaCancel;
		}
		break;
	case 2:
		grayThreshold_CUDA<<<dimGridThreshold, dimBlockThreshold>>>(image, height, width, Max, bottom, top);
		cudaStatus = cudaGetLastError();
		if (cudaStatus != cudaSuccess)
		{
			goto cudaCancel;
		}
		break;
	default:
		break;
	}
cudaCancel:
	// free memory used by find maximum
	free(preMax);
	cudaFree(d_preMax);
	return cudaStatus;
}

/**********************************************************************************************//**
 * @fn	cudaError_t cudaKernelLauncher::cudaBrightness(BYTE *image, int height, int width, int brightness, int threads)
 *
 * @brief	Launcher of CUDA kernel for brightness adjustment of grayscale image.
 *
 * @author	Tomas Bravenec
 *
 * @param [in,out]	image	  	Pointer to first byte of grayscale image.
 * @param 		  	height	  	Height of the image in pixels.
 * @param 		  	width	  	Width of the image in pixels.
 * @param 		  	brightness	The brightness.
 * @param 		  	threads   	Square root of maximum allowed threads per block.
 *
 * @return A cudaError_t containing cudaError code, cudaSucces if no error..
 **************************************************************************************************/
cudaError_t cudaKernelLauncher::cudaBrightness(BYTE *image, int height, int width, int brightness, int threads)
{
	cudaError_t cudaStatus;

	dim3 dimBlock(threads, threads);
	dim3 dimGrid(width / dimBlock.x + 1, height / dimBlock.y + 1);

	adjustBrightness_CUDA<<<dimGrid, dimBlock>>>(image, height, width, brightness);

	// get last CUDA error and return it
	cudaStatus = cudaGetLastError();
	return cudaStatus;
}

/**********************************************************************************************//**
 * @fn	cudaError_t cudaKernelLauncher::cudaNegative(BYTE *image, int height, int width, int threads)
 *
 * @brief	Inverts color of grayscale image using GPU.
 *
 * @author	Tomas Bravenec
 *
 * @param [in,out]	image  	Image data.
 * @param 		  	height 	Height of the image in pixels.
 * @param 		  	width  	Width of the image in pixels.
 * @param 		  	threads	Square root of maximum allowed threads per block.
 *
 * @return A cudaError_t containing cudaError code, cudaSucces if no error..
 **************************************************************************************************/
cudaError_t cudaKernelLauncher::cudaNegative(BYTE *image, int height, int width, int threads)
{
	cudaError_t cudaStatus;

	dim3 dimBlock(threads, threads);
	dim3 dimGrid(width / dimBlock.x + 1, height / dimBlock.y + 1);

	negative_CUDA<<<dimGrid, dimBlock>>>(image, height, width);

	// get last CUDA error and return it
	cudaStatus = cudaGetLastError();
	return cudaStatus;
}

/**********************************************************************************************//**
 * @fn	cudaError_t cudaKernelLauncher::cudaKmeans(BYTE *rgb, int div, int k, int maxIter, int *iterations, int height, int width, int stride, int threads)
 *
 * @brief	Starts kmeans clustering using CUDA enabled GPU.
 *
 * @warning	It is not recomended to set 0 iterations, as 0 means no limit in amount of iterations processing may take very long amount of time. 
 *
 * @author	Tomas Bravenec
 *
 * @param [in,out]	rgb		  	Pointer to first byte of image in 24bpp BGR format.
 * @param 		  	div		  	The divider, says how many times has to be the image split to fit into GPU memory.
 * @param 		  	k		  	Sets number of clusters.
 * @param 		  	maxIter   	Sets maximum iterations (0 means no limit).
 * @param [out]	  	iterations	Return number of iterations needed for Kmeans clustering.
 * @param 		  	height	  	Height of the image in pixels.
 * @param 		  	width	  	Width of the image in pixels.
 * @param 		  	stride	  	Length of one row in bytes, including padding.
 * @param 		  	threads   	Square root of maximum allowed threads per block.
 *
 * @return A cudaError_t containing cudaError code, cudaSucces if no error..
 **************************************************************************************************/
cudaError_t cudaKernelLauncher::cudaKmeans(BYTE *rgb, int div, int k, int maxIter, int *iterations, int height, int width, int stride, int threads)
{
	cudaError_t cudaStatus;

	// initialize all clusters and set maximum iterations
	kmeansInit(k, maxIter);

	dim3 dimBlock(threads, threads);
	dim3 dimGrid(width / dimBlock.x + 1, height / dimBlock.y + 1);

	// allocation of memory for clustered image
	BYTE *clustered = (BYTE*)malloc(stride * height * sizeof(BYTE));
	BYTE *d_rgb, *d_clustered;
	// allocate memory for k clusters on GPU
	cudaCluster *d_clusters;

	int offset;

	// allocate memory for clusters on GPU
	cudaStatus = cudaMalloc(&d_clusters, clusterCount * sizeof(cudaCluster));
	if (cudaStatus != cudaSuccess)
	{
		goto cudaCancel;
	}

	// memory allocation for input and output image
	cudaStatus = cudaMalloc(&d_rgb, stride * height / div * sizeof(BYTE));
	if (cudaStatus != cudaSuccess)
	{
		goto cudaCancel;
	}
	cudaStatus = cudaMalloc(&d_clustered, stride * height / div * sizeof(BYTE));
	if (cudaStatus != cudaSuccess)
	{
		goto cudaCancel;
	}

	// copy randomized clusters to GPU memory
	cudaStatus = cudaMemcpy(d_clusters, newCluster, clusterCount * sizeof(cudaCluster), cudaMemcpyHostToDevice);
	if (cudaStatus != cudaSuccess)
	{
		goto cudaCancel;
	}

	// select infinite iterations or limite amount of iterations 
	if (maxIterations == 0)
	{
		iterations[0] = 0;
		// iterate through image until convergence, or limit of iterations
		while (true)
		{
			for (int j = 0; j < div; j++)
			{
				offset = stride * height / div * j;

				// copy of input image into GPU memory
				cudaStatus = cudaMemcpy(d_rgb, rgb + offset, stride * height / div * sizeof(BYTE), cudaMemcpyHostToDevice);
				if (cudaStatus != cudaSuccess)
				{
					goto cudaCancel;
				}

				// run the iteration of kmeans
				kmeans_iteration_CUDA<<<dimGrid, dimBlock>>>(d_rgb, d_clustered, d_clusters, clusterCount, height / div, width, stride);
				cudaStatus = cudaGetLastError();
				if (cudaStatus != cudaSuccess)
				{
					goto cudaCancel;
				}

				// copy clustered image back into system memory
				cudaStatus = cudaMemcpy(clustered + offset, d_clustered, stride * height / div * sizeof(BYTE), cudaMemcpyDeviceToHost);
				if (cudaStatus != cudaSuccess)
				{
					goto cudaCancel;
				}
			}

			// recalculate the centroids
			kmeans_recalculate_CUDA<<<1, 256>>>(d_clusters, clusterCount);
			cudaStatus = cudaGetLastError();
			if (cudaStatus != cudaSuccess)
			{
				goto cudaCancel;
			}

			// clear unnecesary cluster data
			kmeans_clear_CUDA<<<1, 256>>>(d_clusters, clusterCount);
			cudaStatus = cudaGetLastError();
			if (cudaStatus != cudaSuccess)
			{
				goto cudaCancel;
			}

			// copy new cluster objects into system memory for convergency check
			cudaMemcpy(newCluster, d_clusters, clusterCount * sizeof(cudaCluster), cudaMemcpyDeviceToHost);

			// run convergency check
			if (kmeansCheckConvergence() == true)
			{
				break;
			}

			// create copy of new centroids for future comparison
			kmeansNewToOld();
			iterations[0]++;
		}
	}
	else
	{
		// iterate through image until convergence, or limit of iterations
		for (int i = 0; i <= maxIterations; i++)
		{
			for (int j = 0; j < div; j++)
			{
				offset = stride * height / div * j;

				// copy of input image into GPU memory
				cudaStatus = cudaMemcpy(d_rgb, rgb + offset, stride * height / div * sizeof(BYTE), cudaMemcpyHostToDevice);
				if (cudaStatus != cudaSuccess)
				{
					goto cudaCancel;
				}

				// run the iteration of kmeans
				kmeans_iteration_CUDA<<<dimGrid, dimBlock>>>(d_rgb, d_clustered, d_clusters, clusterCount, height / div, width, stride);
				cudaStatus = cudaGetLastError();
				if (cudaStatus != cudaSuccess)
				{
					goto cudaCancel;
				}

				// copy clustered image back into system memory
				cudaStatus = cudaMemcpy(clustered + offset, d_clustered, stride * height / div * sizeof(BYTE), cudaMemcpyDeviceToHost);
				if (cudaStatus != cudaSuccess)
				{
					goto cudaCancel;
				}
			}

			// recalculate the centroids
			kmeans_recalculate_CUDA<<<1, 256>>>(d_clusters, clusterCount);
			cudaStatus = cudaGetLastError();
			if (cudaStatus != cudaSuccess)
			{
				goto cudaCancel;
			}

			// clear unnecesary cluster data
			kmeans_clear_CUDA<<<1, 256>>>(d_clusters, clusterCount);
			cudaStatus = cudaGetLastError();
			if (cudaStatus != cudaSuccess)
			{
				goto cudaCancel;
			}

			// copy new cluster objects into system memory for convergency check
			cudaMemcpy(newCluster, d_clusters, clusterCount * sizeof(cudaCluster), cudaMemcpyDeviceToHost);

			// run convergency check
			if (kmeansCheckConvergence() == true)
			{
				break;
			}

			// create copy of new centroids for future comparison
			kmeansNewToOld();
			iterations[0] = i;
		}
	}

	// copy clustered image back into input array
	cudaMemcpy(rgb, clustered, stride * height * sizeof(BYTE), cudaMemcpyHostToHost);

	// free all memory and return cuda error code
cudaCancel:
	free(clustered);
	cudaFree(d_rgb);
	cudaFree(d_clusters);
	cudaFree(d_clustered);
	return cudaStatus;
}

/**********************************************************************************************//**
 * @fn	cudaKernelLauncher::cudaKernelLauncher()
 *
 * @brief	Default constructor.
 *
 * @author	Tomas Bravenec
 **************************************************************************************************/
cudaKernelLauncher::cudaKernelLauncher()
{
}

/**********************************************************************************************//**
 * @fn	cudaKernelLauncher::~cudaKernelLauncher()
 *
 * @brief	Destructor.
 *
 * @author	Tomas Bravenec
 **************************************************************************************************/
cudaKernelLauncher::~cudaKernelLauncher()
{
	if (clusterCount >= 2)
	{
		free(newCluster);
		free(oldCluster);
	}
}

/**********************************************************************************************//**
 * @fn	void cudaKernelLauncher::kmeansInit(int k)
 *
 * @brief	Initializes kmeans method.
 *
 * @author	Tomas Bravenec
 *
 * @param	k	K clusters to allocate and randomize.
 **************************************************************************************************/
void cudaKernelLauncher::kmeansInit(int k)
{
	clusterCount = k;

	newCluster = (cudaCluster*)malloc(k * sizeof(cudaCluster));
	oldCluster = (cudaCluster*)malloc(k * sizeof(cudaCluster));

	for (int i = 0; i < k; i++)
	{
		newCluster[i].randomizeCentroid();
	}
}

/**********************************************************************************************//**
 * @fn	void cudaKernelLauncher::kmeansInit(int k, int iterations)
 *
 * @brief	Initializes kmeans method and sets maximum number of iterations.
 *
 * @warning	It is not recomended to set 0 iterations, as 0 means no limit in amount of iterations processing may take very long amount of time.
 *
 * @author	Tomas Bravenec
 *
 * @param	k		  	K clusters to allocate and randomize.
 * @param	iterations	Maximum iterations, if 0 there is no limit.
 **************************************************************************************************/
void cudaKernelLauncher::kmeansInit(int k, int iterations)
{
	clusterCount = k;
	maxIterations = iterations;

	newCluster = (cudaCluster*)malloc(k * sizeof(cudaCluster));
	oldCluster = (cudaCluster*)malloc(k * sizeof(cudaCluster));

	for (int i = 0; i < k; i++)
	{
		newCluster[i].randomizeCentroid();
	}
}

/**********************************************************************************************//**
 * @fn	void cudaKernelLauncher::kmeansSetMaxIterations(int value)
 *
 * @brief	Sets maximum iterations for kmeans method.
 *
 * @warning	It is not recomended to set 0 iterations, as 0 means no limit in amount of iterations processing may take very long amount of time.
 *
 * @author	Tomas Bravenec
 *
 * @param	iterations	Maximum iterations, if 0 there is no limit.
 **************************************************************************************************/
void cudaKernelLauncher::kmeansSetMaxIterations(int iterations)
{
	maxIterations = iterations;
}

/**********************************************************************************************//**
 * @fn	void cudaKernelLauncher::kmeansNewToOld()
 *
 * @brief	Creates a copy of clusters for comparison after next iteration.
 *
 * @author	Tomas Bravenec
 **************************************************************************************************/
void cudaKernelLauncher::kmeansNewToOld()
{
	for (int i = 0; i < clusterCount; i++)
	{
		oldCluster[i] = newCluster[i];
	}
}

/**********************************************************************************************//**
 * @fn	bool cudaKernelLauncher::kmeansCheckConvergence()
 *
 * @brief	Determines if kmeans method converged.
 *
 * @author	Tomas Bravenec
 *
 * @return	True if it succeeds, false if it fails.
 **************************************************************************************************/
bool cudaKernelLauncher::kmeansCheckConvergence()
{
	bool convergence = true;

	// check if centroid of every color in all clusters haven't changed from last iteration
	for (int i = 0; i < clusterCount; i++)
	{
		if (newCluster[i].centroidBlue != oldCluster[i].centroidBlue)
		{
			convergence = false;
		}
		if (newCluster[i].centroidGreen != oldCluster[i].centroidGreen)
		{
			convergence = false;
		}
		if (newCluster[i].centroidRed != oldCluster[i].centroidRed)
		{
			convergence = false;
		}
	}
	return convergence;
}

//
// CUDA kernels
//

/**********************************************************************************************//**
 * @fn	__global__ void rgbToGray_CUDA(BYTE *rgb, BYTE *gray, int height, int width, int stride)
 *
 * @brief	CUDA kernel for conversion of image in BGR format to grayscale.
 * 			
 * @warning	Do not use this from code, use cudaKernelLauncher class instead.
 *
 * @author	Tomas Bravenec
 *
 * @param [in]		rgb   	Pointer to first byte of image in 24bpp BGR format.
 * @param [out]		gray  	Pointer to first byte of grayscale image.
 * @param 		  	height	Height of the image in pixels.
 * @param 		  	width 	Width of the image in pixels.
 * @param 		  	stride	Length of one row in bytes, including padding.
 **************************************************************************************************/
__global__ void rgbToGray_CUDA(BYTE *rgb, BYTE *gray, int height, int width, int stride)
{
	int row = blockIdx.y * blockDim.y + threadIdx.y;
	int col = blockIdx.x * blockDim.x + threadIdx.x;

	// every thread is processing exactly one pixel
	if (row < height && col < width)
	{
		// calculate value of grayscale pixel from RGB
		gray[row * width + col] = 0.3 * rgb[row * stride + 3 * col + 2] + 0.59 * rgb[row * stride + 3 * col + 1] + 0.11 * rgb[row * stride + 3 * col + 0];
	}
}

/**********************************************************************************************//**
 * @fn	__global__ void grayToRGB_CUDA(BYTE *rgb, BYTE *gray, int height, int width, int stride)
 *
 * @brief	CUDA kernel for conversion of grayscale image to BGR format.
 * 			
 * @warning	Do not use this from code, use cudaKernelLauncher class instead.
 *
 * @author	Tomas Bravenec
 *
 * @param [out]		rgb   	Pointer to first byte of image in 24bpp BGR format.
 * @param [in]		gray  	Pointer to first byte of grayscale image.
 * @param 		  	height	Height of the image in pixels.
 * @param 		  	width 	Width of the image in pixels.
 * @param 		  	stride	Length of one row in bytes, including padding.
 **************************************************************************************************/
__global__ void grayToRGB_CUDA(BYTE *rgb, BYTE *gray, int height, int width, int stride)
{
	int row = blockIdx.y * blockDim.y + threadIdx.y;
	int col = blockIdx.x * blockDim.x + threadIdx.x;

	// every thread is processing exactly one pixel
	if (row < height && col < width)
	{
		// set all subpixels to same value for creating 24bpp grayscale - 8bpp cannot be displayed on screen
		rgb[row * stride + 3 * col] = gray[row * width + col];
		rgb[row * stride + 3 * col + 1] = gray[row * width + col];
		rgb[row * stride + 3 * col + 2] = gray[row * width + col];
	}
}

/**********************************************************************************************//**
 * @fn	__global__ void edge_CUDA(BYTE *gray, BYTE *output, int height, int width, int *filterX, int *filterY)
 *
 * @brief	CUDA kernel for detecting edges in grayscale image.
 * 			
 * @warning	Do not use this from code, use cudaKernelLauncher class instead.
 *
 * @author	Tomas Bravenec
 *
 * @param [in]		gray   	Pointer to first byte of grayscale image.
 * @param [out]		output 	If non-null, the output.
 * @param 		  	height 	Height of the image in pixels.
 * @param 		  	width  	Width of the image in pixels.
 * @param [in]		filterX	If non-null, the filter x coordinate.
 * @param [in]		filterY	If non-null, the filter y coordinate.
 **************************************************************************************************/
__global__ void edge_CUDA(BYTE *gray, BYTE *output, int height, int width, int *filterX, int *filterY)
{
	int row = blockIdx.y * blockDim.y + threadIdx.y;
	int col = blockIdx.x * blockDim.x + threadIdx.x;
	int edgeX = 0, edgeY = 0;
	// every thread is processing exactly one pixel
	if (row < height && col < width)
	{
		// nested for loops working on 2D convolution of pixel and its surroundings
		for (int m = 0; m < 3; m++)
		{
			for (int n = 0; n < 3; n++)
			{
				// at edges do not try to include outer nonexisting data
				if ((col - 1 + m) >= 0 && (col - 1 + m) < width && (row - 1 + n) >= 0 && (row - 1 + n) < height)
				{
					// using horizontal filter
					edgeX += gray[(row - 1 + n) * width + (col - 1 + m)] * filterX[n * 3 + m];

					// using vertical filter
					edgeY += gray[(row - 1 + n) * width + (col - 1 + m)] * filterY[n * 3 + m];
				}
			}
		}
		// using faster aproximate of square root of sum of squares
		// saturation control
		if ((abs(edgeX) + abs(edgeY)) > 255)
		{
			output[row * width + col] = 255;
		}
		else
		{
			output[row * width + col] = abs(edgeX) + abs(edgeY);
		}
	}
}

/**********************************************************************************************//**
 * @fn	__global__ void maxCol_CUDA(BYTE *input, BYTE *output, int height, int width)
 *
 * @brief	CUDA kernel for finding maximum in columns.
 * 			
 * @warning	Do not use this from code, use cudaKernelLauncher class instead.
 *
 * @author	Tomas Bravenec
 *
 * @param [in]		input 	If non-null, the input.
 * @param [out]		output	If non-null, the output.
 * @param 		  	height	Height of the image in pixels.
 * @param 		  	width 	Width of the image in pixels.
 **************************************************************************************************/
__global__ void maxCol_CUDA(BYTE *input, BYTE *output, int height, int width)
{
	int col = blockIdx.x * blockDim.x + threadIdx.x;

	// every thread is processing one column
	if (col < width)
	{
		// find maximum for every column
		for (int row = 0; row < height; row++)
		{
			if (output[col] < input[row * width + col])
			{
				output[col] = input[row * width + col];
			}
		}
	}
}

/**********************************************************************************************//**
 * @fn	__global__ void maxRow_CUDA(BYTE *input, BYTE *output, int height, int width)
 *
 * @brief	CUDA kernel for finding maximum in rows.
 * 			
 * @warning	Do not use this from code, use cudaKernelLauncher class instead.
 *
 * @author	Tomas Bravenec
 *
 * @param [in]		input 	If non-null, the input.
 * @param [out]		output	If non-null, the output.
 * @param 		  	height	Height of the image in pixels.
 * @param 		  	width 	Width of the image in pixels.
 **************************************************************************************************/
__global__ void maxRow_CUDA(BYTE *input, BYTE *output, int height, int width)
{
	int row = blockIdx.y * blockDim.y + threadIdx.y;

	// every thread is processing one column
	if (row < height)
	{
		// find maximum for every column
		for (int col = 0; col < width; col++)
		{
			if (output[row] < input[row * width + col])
			{
				output[row] = input[row * width + col];
			}
		}
	}
}

/**********************************************************************************************//**
 * @fn	__global__ void simpleThresholdEdge_CUDA(BYTE *image, int height, int width, BYTE max, double threshold)
 *
 * @brief	CUDA kernel for normalization to interval and simple tresholding after Edge detection.
 * 			
 * @warning	Do not use this from code, use cudaKernelLauncher class instead.
 *
 * @author	Tomas Bravenec
 *
 * @param [in,out]	image	 	Pointer to first byte of grayscale image.
 * @param 		  	height   	Height of the image in pixels.
 * @param 		  	width	 	Width of the image in pixels.
 * @param 		  	max		 	The maximum.
 * @param 		  	threshold	The threshold.
 **************************************************************************************************/
__global__ void simpleThresholdEdge_CUDA(BYTE *image, int height, int width, BYTE max, double threshold)
{
	int row = blockIdx.y * blockDim.y + threadIdx.y;
	int col = blockIdx.x * blockDim.x + threadIdx.x;

	// every thread is processing exactly one pixel
	if (row < height && col < width)
	{
		// normalize image pixels
		image[row * width + col] = image[row * width + col] * (255.0 / max);

		// treshold every pixel
		if (image[row * width + col] < (255.0 * threshold))
		{
			image[row * width + col] = 0;
		}
	}
}

/**********************************************************************************************//**
 * @fn	__global__ void simpleThreshold_CUDA(BYTE *image, int height, int width, BYTE max, double threshold)
 *
 * @brief	CUDA kernel for normalization to interval and simple tresholding.
 * 			
 * @warning	Do not use this from code, use cudaKernelLauncher class instead.
 *
 * @author	Tomas Bravenec
 *
 * @param [in,out]	image	 	Pointer to first byte of grayscale image.
 * @param 		  	height   	Height of the image in pixels.
 * @param 		  	width	 	Width of the image in pixels.
 * @param 		  	max		 	The maximum.
 * @param 		  	threshold	The threshold.
 **************************************************************************************************/
__global__ void simpleThreshold_CUDA(BYTE *image, int height, int width, BYTE max, double threshold)
{
	int row = blockIdx.y * blockDim.y + threadIdx.y;
	int col = blockIdx.x * blockDim.x + threadIdx.x;

	// every thread is processing exactly one pixel
	if (row < height && col < width)
	{
		// normalize image pixels
		image[row * width + col] = image[row * width + col] * (255.0 / max);

		// treshold every pixel
		if (image[row * width + col] < (255.0 * threshold))
		{
			image[row * width + col] = 0;
		}
		else
		{
			image[row * width + col] = 255;
		}
	}
}

/**********************************************************************************************//**
 * @fn	__global__ void binaryThreshold_CUDA(BYTE *image, int height, int width, BYTE max, double bottom, double top)
 *
 * @brief	CUDA kernel for normalization to interval and dual binary tresholding.
 * 			
 * @warning	Do not use this from code, use cudaKernelLauncher class instead.
 * 			
 * @author	Tomas Bravenec
 *
 * @param [in,out]	image 	Pointer to first byte of grayscale image.
 * @param 		  	height	Height of the image in pixels.
 * @param 		  	width 	Width of the image in pixels.
 * @param 		  	max   	The maximum.
 * @param 		  	bottom	The bottom.
 * @param 		  	top   	The top.
 **************************************************************************************************/
__global__ void binaryThreshold_CUDA(BYTE *image, int height, int width, BYTE max, double bottom, double top)
{
	int row = blockIdx.y * blockDim.y + threadIdx.y;
	int col = blockIdx.x * blockDim.x + threadIdx.x;

	// every thread is processing exactly one pixel
	if (row < height && col < width)
	{
		// normalize image pixels
		image[row * width + col] = image[row * width + col] * (255.0 / max);

		// treshold every pixel
		if (image[row * width + col] < (255.0 * bottom) || image[row * width + col] > (255.0 * top))
		{
			image[row * width + col] = 0;
		}
		else
		{
			image[row * width + col] = 255;
		}
	}
}

/**********************************************************************************************//**
 * @fn	__global__ void grayThreshold_CUDA(BYTE *image, int height, int width, BYTE max, double bottom, double top)
 *
 * @brief	CUDA kernel for normalization to interval and dual gray tresholding.
 * 			
 * @warning	Do not use this from code, use cudaKernelLauncher class instead.
 *
 * @author	Tomas Bravenec
 *
 * @param [in,out]	image 	Pointer to first byte of grayscale image.
 * @param 		  	height	Height of the image in pixels.
 * @param 		  	width 	Width of the image in pixels.
 * @param 		  	max   	The maximum.
 * @param 		  	bottom	The bottom.
 * @param 		  	top   	The top.
 **************************************************************************************************/
__global__ void grayThreshold_CUDA(BYTE *image, int height, int width, BYTE max, double bottom, double top)
{
	int row = blockIdx.y * blockDim.y + threadIdx.y;
	int col = blockIdx.x * blockDim.x + threadIdx.x;

	// every thread is processing exactly one pixel
	if (row < height && col < width)
	{
		// normalize image pixels
		image[row * width + col] = image[row * width + col] * (255.0 / max);

		// treshold every pixel
		if (image[row * width + col] < (255.0 * bottom))
		{
			image[row * width + col] = 0;
		}
		else if (image[row * width + col] > (255.0 * top))
		{
			image[row * width + col] = 255;
		}
		else
		{
			image[row * width + col] = 128;
		}
	}
}

/**********************************************************************************************//**
 * @fn	__global__ void negative_CUDA(BYTE *image, int height, int width)
 *
 * @brief	CUDA kernel for color inversion of grayscale image.
 * 			
 * @warning	Do not use this from code, use cudaKernelLauncher class instead.
 *
 * @author	Tomas Bravenec
 *
 * @param [in,out]	image 	Pointer to first byte of grayscale image.
 * @param 		  	height	Height of the image in pixels.
 * @param 		  	width 	Width of the image in pixels.
 **************************************************************************************************/
__global__ void negative_CUDA(BYTE *image, int height, int width)
{
	int row = blockIdx.y * blockDim.y + threadIdx.y;
	int col = blockIdx.x * blockDim.x + threadIdx.x;

	// every thread is processing exactly one pixel
	if (row < height && col < width)
	{
		image[row * width + col] = 255 - image[row * width + col];
	}
}

/**********************************************************************************************//**
 * @fn	__global__ void adjustBrightness_CUDA(BYTE *image ,int height, int width, int bright)
 *
 * @brief	CUDA kernel for brightness adjustment in grayscale image.
 * 			
 * @warning	Do not use this from code, use cudaKernelLauncher class instead.
 *
 * @author	Tomas Bravenec
 *
 * @param [in,out]	image 	Pointer to first byte of grayscale image.
 * @param 		  	height	Height of the image in pixels.
 * @param 		  	width 	Width of the image in pixels.
 * @param 		  	bright	The bright.
 **************************************************************************************************/
__global__ void adjustBrightness_CUDA(BYTE *image ,int height, int width, int bright)
{
	int row = blockIdx.y * blockDim.y + threadIdx.y;
	int col = blockIdx.x * blockDim.x + threadIdx.x;

	// every thread is processing exactly one pixel
	if (row < height && col < width)
	{
		// check for saturation, else add bright to pixel
		if (image[row * width + col] + bright > 255)
		{
			image[row * width + col] = 255;
		}
		else if (image[row * width + col] + bright < 0)
		{
			image[row * width + col] = 255;
		}
		else if (image[row * width + col] > 128)
		{
			image[row * width + col] = image[row * width + col] + bright;
		}
	}
}

/**********************************************************************************************//**
 * @fn	__global__ void kmeans_iteration_CUDA(BYTE *rgb, BYTE *clustered, cudaCluster *clusters, BYTE clusterCount, int height, int width, int stride)
 *
 * @brief	CUDA kernel running a labeling phase of kmeans method.
 * 			
 * @warning	Do not use this from code, use method start from kmeansCUDA class instead.
 *
 * @author	Tomas Bravenec
 *
 * @param [in]		rgb				Pointer to first byte of image in 24bpp BGR format.
 * @param [out]		clustered   	If non-null, the clustered.
 * @param [in]		clusters		If non-null, the clusters.
 * @param 		  	clusterCount	Number of clusters.
 * @param 		  	height			Height of the image in pixels.
 * @param 		  	width			Width of the image in pixels.
 * @param 		  	stride			Length of one row in bytes, including padding.
 **************************************************************************************************/
__global__ void kmeans_iteration_CUDA(BYTE *rgb, BYTE *clustered, cudaCluster *clusters, BYTE clusterCount, int height, int width, int stride)
{
	int row = blockIdx.y * blockDim.y + threadIdx.y;
	int col = blockIdx.x * blockDim.x + threadIdx.x;

	// helper variables for kmeans clustering
	int shortestDist = 0;
	BYTE shortestCluster = 0;

	// every thread is processing exactly one pixel
	if (row < height && col < width)
	{
		// calculate approximation of distance for every cluster and pick the closest
		for (int k = 0; k < clusterCount; k++)
		{
			// first cluster set as closest and then proceed with comparation of distance with others
			if (k == 0)
			{
				shortestDist = powf(rgb[row * stride + 3 * col + 0] - clusters[k].centroidBlue, 2) +
					powf(rgb[row * stride + 3 * col + 1] - clusters[k].centroidGreen, 2) +
					powf(rgb[row * stride + 3 * col + 2] - clusters[k].centroidRed, 2);
				shortestCluster = k;
			}
			else if ((powf(rgb[row * stride + 3 * col + 0] - clusters[k].centroidBlue, 2) +
				powf(rgb[row * stride + 3 * col + 1] - clusters[k].centroidGreen, 2) +
				powf(rgb[row * stride + 3 * col + 2] - clusters[k].centroidRed, 2)) < shortestDist)
			{
				shortestDist = powf(rgb[row * stride + 3 * col + 0] - clusters[k].centroidBlue, 2) +
					powf(rgb[row * stride + 3 * col + 1] - clusters[k].centroidGreen, 2) +
					powf(rgb[row * stride + 3 * col + 2] - clusters[k].centroidRed, 2);
				shortestCluster = k;
			}
		}
		
		// use atomic operations to serialize process of counting data
		atomicAdd(&clusters[shortestCluster].sumBlue, rgb[row * stride + 3 * col + 0]);
		atomicAdd(&clusters[shortestCluster].sumGreen, rgb[row * stride + 3 * col + 1]);
		atomicAdd(&clusters[shortestCluster].sumRed, rgb[row * stride + 3 * col + 2]);
		atomicInc(&clusters[shortestCluster].count, -1);

		// fill output array with segmented image
		clustered[row * stride + 3 * col + 0] = clusters[shortestCluster].centroidBlue;
		clustered[row * stride + 3 * col + 1] = clusters[shortestCluster].centroidGreen;
		clustered[row * stride + 3 * col + 2] = clusters[shortestCluster].centroidRed;
	}
}

/**********************************************************************************************//**
 * @fn	__global__ void kmeans_recalculate_CUDA(cudaCluster *clusters, BYTE clusterCount)
 *
 * @brief	CUDA kernel for recalculation of centroids
 * 			
 * @warning	Do not use this from code, use method start from kmeansCUDA class instead.
 *
 * @author	Tomas Bravenec
 *
 * @param [in,out]	clusters		If non-null, the clusters.
 * @param 		  	clusterCount	Number of clusters.
 **************************************************************************************************/
__global__ void kmeans_recalculate_CUDA(cudaCluster *clusters, BYTE clusterCount)
{
	int k = blockIdx.x * blockDim.x + threadIdx.x;

	// every thread is processing single cluster
	if (k < clusterCount)
	{
		clusters[k].centroidBlue = clusters[k].sumBlue / clusters[k].count;
		clusters[k].centroidGreen = clusters[k].sumGreen / clusters[k].count;
		clusters[k].centroidRed = clusters[k].sumRed / clusters[k].count;
	}
}

/**********************************************************************************************//**
 * @fn	__global__ void kmeans_clear_CUDA(cudaCluster *clusters, BYTE clusterCount)
 *
 * @brief	CUDA kernel for clearing all unnecessary data from clusters.
 * 			
 * @warning	Do not use this from code, use method start from kmeansCUDA class instead.
 * 
 * @author	Tomas Bravenec
 *
 * @param [in,out]	clusters		If non-null, the clusters.
 * @param 		  	clusterCount	Number of clusters.
 **************************************************************************************************/
__global__ void kmeans_clear_CUDA(cudaCluster *clusters, BYTE clusterCount)
{
	int k = blockIdx.x * blockDim.x + threadIdx.x;

	// every thread is processing single cluster
	if (k < clusterCount)
	{
		clusters[k].sumBlue = 0;
		clusters[k].sumGreen = 0;
		clusters[k].sumRed = 0;
		clusters[k].count = 0;
	}
}