/**********************************************************************************************//**
 * @file	CPU_Func_Library.cpp
 * 			
 * @author	Tomas Bravenec
 *
 * @brief	Implements the CPU function library class.
 **************************************************************************************************/

#include "stdafx.h"
#include "CPU_Func_Library.h"

/**********************************************************************************************//**
 * @fn	void ImageProcessingCPU::RGB2Gray(array<Byte>^ rgb, array<Byte, 2>^ gray, int height, int width, int stride, int threads)
 *
 * @brief	Function for conversion of image in BGR format to grayscale.
 *
 * @author	Tomas Bravenec
 *
 * @param [in]		rgb	   	1 dimensional managed array containing RGB image.
 * @param [out]		gray   	2 dimensional managed array containing 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	Maximum threads allowed for computing.
 **************************************************************************************************/
void ImageProcessingCPU::RGB2Gray(array<Byte>^ rgb, array<Byte, 2>^ gray, int height, int width, int stride, int threads)
{
	// disable dynamic changes in max threads usage
	omp_set_dynamic(0);
	// set maximum of allowed threads
	omp_set_num_threads(threads);
	// code for multithreaded for
#pragma omp parallel for
	// iterate through rows
	for (int i = 0; i < height; i++)
	{
		// iterate through every column in each row
		for (int j = 0; j < width; j++)
		{
			// fill 2D array with 8bpp grayscale image
			gray[i, j] = 0.3 * rgb[i * stride + 3 * j + 2] + 0.59 * rgb[i * stride + 3 * j + 1] + 0.11 * rgb[i * stride + 3 * j + 0];
		}
	}
}

/**********************************************************************************************//**
 * @fn	void ImageProcessingCPU::Gray2RGB(array<Byte>^ rgb, array<Byte, 2>^ output, int height, int width, int stride, int threads)
 *
 * @brief	Function for conversion of grayscale image to BGR format.
 *
 * @author	Tomas Bravenec
 *
 * @param [out]		rgb	   	1 dimensional managed array containing RGB image.
 * @param [in]		output 	2 dimensional managed array containing 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	Maximum threads allowed for computing.
 **************************************************************************************************/
void ImageProcessingCPU::Gray2RGB(array<Byte>^ rgb, array<Byte, 2>^ gray, int height, int width, int stride, int threads)
{
	// disable dynamic changes in max threads usage
	omp_set_dynamic(0);
	// set maximum of allowed threads
	omp_set_num_threads(threads);
	// code for multithreaded for
#pragma omp parallel for
	// iterate through rows
	for (int i = 0; i < height; i++)
	{
		// iterate through every column in each row
		for (int j = 0; j < width; j++)
		{
			// fill every pixels color with same data for 24 bpp grayscale image 
			rgb[i * stride + 3 * j] = gray[i, j];			// blue
			rgb[i * stride + 3 * j + 1] = gray[i, j];		// green
			rgb[i * stride + 3 * j + 2] = gray[i, j];		// red
		}
	}
}

/**********************************************************************************************//**
 * @fn	void ImageProcessingCPU::edgeDetection(array<Byte, 2>^ gray, array<Byte, 2>^ edge, int filter, int threads)
 *
 * @brief	Function for detecting edges in grayscale image.
 *
 * @author	Tomas Bravenec
 *
 * @param [in]		gray   	2 dimensional managed array containing grayscale image.
 * @param [out]		edge   	2 dimensional managed array containing output image with detected edges.
 * @param			filter 	Specifies the filter:
 * 							0: Sobel
 * 							1: Prewitt  
 * 							2: Kirsch
 * @param			threads	Maximum threads allowed for computing.
 **************************************************************************************************/
void ImageProcessingCPU::edgeDetection(array<Byte, 2>^ gray, array<Byte, 2>^ edge, int filter, int threads)
{
	//first index means:
	// 0: Prewitt.x
	// 1: Prewitt.y
	// 2: Sobel.x
	// 3: Sobel.y
	// 4: Kirsch.x
	// 5: Kirsch.y
	filter_type filters = { { -1, 0, 1,
		-1, 0, 1,
		-1, 0, 1 },
		{ 1, 1, 1,
		0, 0, 0,
		-1, -1, -1 },
		{ -1, 0, 1,
		-2, 0, 2,
		-1, 0, 1 },
		{ 1, 2, 1,
		0, 0, 0,
		-1, -2, -1 },
		{ -5, 3, 3,
		-5, 0, 3,
		-5, 3, 3 },
		{ 3, 3, 3,
		3, 0, 3,
		-5, -5, -5 }, };

	array<Int16, 2>^ edgeX = gcnew array<Int16, 2>(edge->GetLength(0), edge->GetLength(1));
	array<Int16, 2>^ edgeY = gcnew array<Int16, 2>(edge->GetLength(0), edge->GetLength(1));
	int fx, fy;

	// filter selection
	switch (filter)
	{
	case 0:
		fx = 2;
		fy = fx++;
		break;
	case 1:
		fx = 0;
		fy = fx++;
		break;
	case 2:
		fx = 4;
		fy = fx++;
		break;
	default:
		fx = 2;
		fy = fx++;
		break;
	}

	// disable dynamic changes in max threads usage
	omp_set_dynamic(0);
	// set maximum of allowed threads
	omp_set_num_threads(threads);
	// code for multithreaded for
#pragma omp parallel for 
	// iterate through rows
	for (int w = 0; w < edge->GetLength(0); w++)
	{
		// iterate through every column in each row
		for (int h = 0; h < edge->GetLength(1); h++)
		{
			// 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 ((w - 1 + m) >= 0 && (w - 1 + m) < edge->GetLength(0) && (h - 1 + n) >= 0 && (h - 1 + n) < edge->GetLength(1))
					{
						// using horizontal filter
						edgeX[w, h] += gray[w - 1 + m, h - 1 + n] * filters[fx][2 - m][2 - n];
						
						// using vertical filter
						edgeY[w, h] += gray[w - 1 + m, h - 1 + n] * filters[fy][2 - m][2 - n];
					}
				}
			}
			// using faster aproximate of square root of sum of squares
			// saturation control
			if ((System::Math::Abs(edgeX[w, h]) + System::Math::Abs(edgeY[w, h])) > 255)
			{
				edge[w, h] = 255;				
			}
			else
			{
				edge[w, h] = System::Math::Abs(edgeX[w, h]) + System::Math::Abs(edgeY[w, h]);
			}
		}
	}
}

/**********************************************************************************************//**
 * @fn	void ImageProcessingCPU::thresholdEdge(array<Byte, 2>^ image, double threshold, int threads)
 *
 * @brief	Function responsible for thresholding of grayscale image after edge detection
 *
 * @author	Tomas Bravenec
 *
 * @param [in, out]	image	 	2 dimensional managed array containing grayscale image to be tresholded.
 * @param			threshold	The threshold.
 * @param			threads  	Maximum threads allowed for computing.
 **************************************************************************************************/
void ImageProcessingCPU::thresholdEdge(array<Byte, 2>^ image, double threshold, int threads)
{
	double max = 0;
	// disable dynamic changes in max threads usage
	omp_set_dynamic(0);
	// set maximum of allowed threads
	omp_set_num_threads(threads);
	// code for multithreaded for
#pragma omp parallel for
	// iterate through rows
	for (int w = 0; w < image->GetLength(0); w++)
	{
		// iterate through every column in each row
		for (int h = 0; h < image->GetLength(1); h++)
		{
			if (image[w, h] > max)
			{
				max = image[w, h];
			}
		}
	}

#pragma omp parallel for
	// iterate through rows
	for (int w = 0; w < image->GetLength(0); w++)
	{
		// iterate through every column in each row
		for (int h = 0; h < image->GetLength(1); h++)
		{
			image[w, h] = (image[w, h] * (255.0 / max));
			if (image[w, h] < (255.0 * threshold))
			{
				image[w, h] = 0;
			}
		}
	}
}

/**********************************************************************************************//**
 * @fn	void ImageProcessingCPU::threshold(array<Byte, 2>^ image, double top, double bottom, int type, int threads)
 *
 * @brief	Functions responsible for thresholding of grayscale image
 *
 * @author	Tomas Bravenec
 *
 * @param [in, out]	image  	2 dimensional managed array containing grayscale image to be tresholded.
 * @param			top	   	Upper treshold.
 * @param			bottom 	Lower treshold.
 * @param			type	Specifies tresholding type:
 * 							0: Simple
 * 							1: Dual Binary  
 * 							2: Dual Grayscale
 * @param			threads	Maximum threads allowed for computing.
 **************************************************************************************************/
void ImageProcessingCPU::threshold(array<Byte, 2>^ image, double top, double bottom, int type, int threads)
{
	double max = 0;
	// disable dynamic changes in max threads usage
	omp_set_dynamic(0);
	// set maximum of allowed threads
	omp_set_num_threads(threads);
	// code for multithreaded for
#pragma omp parallel for
	// iterate through rows
	for (int w = 0; w < image->GetLength(0); w++)
	{
		// iterate through every column in each row
		for (int h = 0; h < image->GetLength(1); h++)
		{
			if (image[w, h] > max)
			{
				max = image[w, h];
			}
		}
	}

#pragma omp parallel for
	// iterate through rows
	for (int w = 0; w < image->GetLength(0); w++)
	{
		// iterate through every column in each row
		for (int h = 0; h < image->GetLength(1); h++)
		{
			image[w, h] = (image[w, h] * (255.0 / max));
			switch (type)
			{
			case 1:
				if (image[w, h] <= (255.0 * bottom) || image[w, h] > (255.0 * top))
				{
					image[w, h] = 0;
				}
				else
				{
					image[w, h] = 255;
				}
				break;
			case 2:
				if (image[w, h] <= (255.0 * bottom))
				{
					image[w, h] = 0;
				}
				else if (image[w, h] > (255.0 * top))
				{
					image[w, h] = 255;
				}
				else
				{
					image[w, h] = 128;
				}
				break;
			default:
				if (image[w, h] < (255.0 * top))
				{
					image[w, h] = 0;
				}
				else
				{
					image[w, h] = 255;
				}
				break;
			}
		}
	}
}

/**********************************************************************************************//**
 * @fn	void ImageProcessingCPU::negative(array<Byte, 2>^ image, int threads)
 *
 * @brief	Function for conversion of grayscale image to negative
 *
 * @author	Tomas Bravenec
 *
 * @param [in, out]	image  	The image to be converted to negative.
 * @param			threads	Maximum threads allowed for computing.
 **************************************************************************************************/
void ImageProcessingCPU::negative(array<Byte, 2>^ image, int threads)
{
	double max = 0;
	// disable dynamic changes in max threads usage
	omp_set_dynamic(0);
	// set maximum of allowed threads
	omp_set_num_threads(threads);
	// code for multithreaded for
#pragma omp parallel for
	// iterate through rows
	for (int w = 0; w < image->GetLength(0); w++)
	{
		// iterate through every column in each row
		for (int h = 0; h < image->GetLength(1); h++)
		{
			image[w, h] = 255 - image[w, h];
		}
	}
}

/**********************************************************************************************//**
 * @fn	void ImageProcessingCPU::brightness(array<Byte, 2>^ image, int brightness, int threads)
 *
 * @brief	Function for brightness adjustion
 *
 * @author	Tomas Bravenec
 *
 * @param [in, out]	image	  	2 dimensional managed array containing grayscale image
 * 								for brightness adjusting.
 * @param			brightness	The brightness adjustment step.
 * @param			threads   	Maximum threads allowed for computing.
 **************************************************************************************************/
void ImageProcessingCPU::brightness(array<Byte, 2>^ image, int brightness, int threads)
{
	// disable dynamic changes in max threads usage
	omp_set_dynamic(0);
	// set maximum of allowed threads
	omp_set_num_threads(threads);
	// code for multithreaded for
#pragma omp parallel for
	// iterate through rows
	for (int w = 0; w < image->GetLength(0); w++)
	{
		// iterate through every column in each row
		for (int h = 0; h < image->GetLength(1); h++)
		{
			if (image[w, h] + brightness > 255)
			{
				image[w, h] = 255;
			}
			else if (image[w, h] + brightness < 0)
			{
				image[w, h] = 255;
			}
			else if (image[w, h] > 128)
			{
				image[w, h] = image[w, h] + brightness;
			}
		}
	}
}

/**********************************************************************************************//**
 * @fn	void ImageProcessingCPU::histData(array<Byte, 2>^ gray, array<double>^ hist)
 *
 * @brief	Function for getting histogram data from a geayscale image
 *
 * @author	Tomas Bravenec
 *
 * @param [in]		gray	2 dimensional managed array containing grayscale image
 * @param [out]		hist	1 dimensional managed array for histogram data.
 **************************************************************************************************/
void ImageProcessingCPU::histData(array<Byte, 2>^ gray, array<double>^ hist)
{
	for (int w = 0; w < gray->GetLength(0); w++)
	{
		// iterate through every column in each row
		for (int h = 0; h < gray->GetLength(1); h++)
		{
			hist[gray[w, h]]++;
		}
	}
}

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

/**********************************************************************************************//**
 * @fn	void Cluster::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 Cluster::setCentroid(BYTE Blue, BYTE Green, BYTE Red)
{
	centroidRed = Red;
	centroidGreen = Green;
	centroidBlue = Blue;
}

/**********************************************************************************************//**
 * @fn	void Cluster::recalculateCentroids()
 *
 * @brief	Recalculates all cluster centroids
 *
 * @author	Tomas Bravenec
 **************************************************************************************************/
void Cluster::recalculateCentroids()
{
	if (count != 0)
	{
		centroidBlue = sumBlue / count;
		centroidRed = sumRed / count;
		centroidGreen = sumGreen / count;
	}
	clearValues();
}

/**********************************************************************************************//**
 * @fn	void Cluster::clearValues()
 *
 * @brief	Clears all unnecessarz data from clusters.
 *
 * @author	Tomas Bravenec
 **************************************************************************************************/
void Cluster::clearValues()
{
	sumBlue = 0;
	sumGreen = 0;
	sumRed = 0;
	count = 0;
}

/**********************************************************************************************//**
 * @fn	ImageProcessingCPU::ImageProcessingCPU()
 *
 * @brief	Constructor.
 *
 * @author	Tomas Bravenec
 **************************************************************************************************/
ImageProcessingCPU::ImageProcessingCPU()
{

}

/**********************************************************************************************//**
 * @fn	ImageProcessingCPU::~ImageProcessingCPU()
 *
 * @brief	Destructor.
 *
 * @author	Tomas Bravenec
 **************************************************************************************************/
ImageProcessingCPU::~ImageProcessingCPU()
{
	free(newCluster);
	free(oldCluster);
}

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

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

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

/**********************************************************************************************//**
 * @fn	void ImageProcessingCPU::init(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 ImageProcessingCPU::init(int k, int iterations)
{
	clusterCount = k;
	maxIterations = iterations;

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

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

/**********************************************************************************************//**
 * @fn	void ImageProcessingCPU::setMaxIterations(int iterations)
 *
 * @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 ImageProcessingCPU::setMaxIterations(int iterations)
{
	maxIterations = iterations;
}

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

	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;
}

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

/**********************************************************************************************//**
 * @fn	void ImageProcessingCPU::iterate(array<Byte>^ rgb, array<Byte>^ clustered, int height, int width, int stride, int threads)
 *
 * @brief	Runs one iteration of kmeans clustering.
 *
 * @author	Tomas Bravenec
 *
 * @param [in]		rgb		 	1 dimensional managed array containing RGB image.
 * @param [out]		clustered	1 dimensional managed array containing clustered 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  	Maximum threads allowed for computing.
 **************************************************************************************************/
void ImageProcessingCPU::iterate(array<Byte>^ rgb, array<Byte>^ clustered, int height, int width, int stride, int threads)
{
	// disable dynamic changes in max threads usage
	omp_set_dynamic(0);
	// set maximum of allowed threads
	omp_set_num_threads(threads);
	// code for multithreaded for
#pragma omp parallel for
	// iterate through rows

	for (int i = 0; i < height; i++)
	{
		// iterate through every column in each row
		for (int j = 0; j < width * 3; j += 3)
		{
			int shortestDist = 0, shortestCluster = 0;
			for (int k = 0; k < clusterCount; k++)
			{
				if (k == 0)
				{
					shortestDist = (System::Math::Pow(rgb[i * stride + j + 0] - newCluster[k].centroidBlue, 2) +
						System::Math::Pow(rgb[i * stride + j + 1] - newCluster[k].centroidGreen, 2) +
						System::Math::Pow(rgb[i * stride + j + 2] - newCluster[k].centroidRed, 2));
					shortestCluster = k;
				}
				if ((System::Math::Pow(rgb[i * stride + j + 0] - newCluster[k].centroidBlue, 2) +
					System::Math::Pow(rgb[i * stride + j + 1] - newCluster[k].centroidGreen, 2) +
					System::Math::Pow(rgb[i * stride + j + 2] - newCluster[k].centroidRed, 2)) < shortestDist)
				{
					shortestDist = (System::Math::Pow(rgb[i * stride + j + 0] - newCluster[k].centroidBlue, 2) +
						System::Math::Pow(rgb[i * stride + j + 1] - newCluster[k].centroidGreen, 2) +
						System::Math::Pow(rgb[i * stride + j + 2] - newCluster[k].centroidRed, 2));
					shortestCluster = k;
				}
			}
			newCluster[shortestCluster].addBlueValue(rgb[i * stride + j + 0]);
			newCluster[shortestCluster].addGreenValue(rgb[i * stride + j + 1]);
			newCluster[shortestCluster].addRedValue(rgb[i * stride + j + 2]);
			newCluster[shortestCluster].increaseCount();

			clustered[i * stride + j + 0] = newCluster[shortestCluster].centroidBlue;
			clustered[i * stride + j + 1] = newCluster[shortestCluster].centroidGreen;
			clustered[i * stride + j + 2] = newCluster[shortestCluster].centroidRed;
		}
	}
}

/**********************************************************************************************//**
 * @fn	void ImageProcessingCPU::recalculateCentroids()
 *
 * @brief	Recalculates the centroids.
 *
 * @author	Tomas Bravenec
 **************************************************************************************************/
void ImageProcessingCPU::recalculateCentroids()
{
	for (int i = 0; i < clusterCount; i++)
	{
		newCluster[i].recalculateCentroids();
	}
}

/**********************************************************************************************//**
 * @fn	int ImageProcessingCPU::kmeans(array<Byte>^ rgb, array<Byte>^ clustered, BYTE k, int maxIter, int height, int width, int stride, int threads)
 *
 * @brief	Starts kmeans clustering using CPU.
 *
 * @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]		rgb		 	1 dimensional managed array containing RGB image.
 * @param [out]		clustered	1 dimensional managed array containing clustered image.
 * @param			k		 	Sets number of clusters.
 * @param			maxIter  	Sets maximum iterations (0 means no limit).
 * @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  	Maximum threads allowed for computing.
 *
 * @return	Return number of iterations needed for Kmeans clustering.
 **************************************************************************************************/
int ImageProcessingCPU::kmeans(array<Byte>^ rgb, array<Byte>^ clustered, BYTE k, int maxIter, int height, int width, int stride, int threads)
{
	init(k);
	setMaxIterations(maxIter);

	int numberOfIterations = 0;

	if (maxIterations == 0)
	{
		while (true)
		{
			iterate(rgb, clustered, height, width, stride, threads);
			recalculateCentroids();
			if (checkConvergence() == true)
			{
				break;
			}
			newToOld();
			numberOfIterations++;
		}
	}
	else
	{
		for (numberOfIterations = 0; numberOfIterations < maxIterations; numberOfIterations++)
		{
			iterate(rgb, clustered, height, width, stride, threads);
			recalculateCentroids();
			if (checkConvergence() == true)
			{
				break;
			}
			newToOld();
		}
	}

	return numberOfIterations;
}