// get Gaussian Kernel Function
void ycGaussianKernel(float* kernel, int sigma, int slab)
{
	int index;
	float dx2, dy2;
	float sum = 0;
	for(int i=0; i<slab; i++)
	{
		dx2 = pow(i - (slab - 1)/2.0, 2);
		for(int j=0; j<slab; j++)
		{
			dy2 = pow(j - (slab - 1)/2.0, 2);
			index = i*slab + j;
			kernel[index] = exp(-(dx2 + dy2)/2/pow(sigma, 2)) / (2*PI*pow(sigma, 2));
			//printf("%f\n", kernel[index]);
			sum += kernel[index];
		}
	}

	for(int k=0; k< slab*slab; k++)
	{
		kernel[k] = kernel[k] / sum;
		//printf("%f\n", kernel[k]);
	}
}

// Gaussian filter
Mat filterGaussian(cv::Mat img, const float sigma, const int slab)
{
	cvtColor(img, img, CV_BGR2GRAY);
	Mat retMat = Mat::zeros(img.rows, img.cols, CV_8UC1);
	for(int i=0; i<img.rows; i++)
		for(int j=0; j<img.cols; j++)
			retMat.at<uchar>(i, j) = img.at<uchar>(i, j);

	// 一维数组模拟二维数组
	float* kernel = new float[slab * slab];
	int* xLocation = new int[slab];
	int* yLocation = new int[slab];
	ycGaussianKernel(kernel, sigma, slab);
	float sum;
	int index;

	// 对于边缘，这里采取直接舍弃不计算的方法。因此，循环起点是 slab/2
	for(int i= slab/2; i<img.rows - slab/2; i++)
	{
		xLocation[slab/2] = i;
		for(int delta = 0; delta <= slab/2; delta++)
		{
			xLocation[slab/2 - delta] = i - delta;
			xLocation[slab/2 + delta] = i + delta;
		}

		for(int j= slab/2; j<img.cols - slab/2; j++)
		{
			yLocation[slab/2] = j;
			for(int delta = 0; delta <= slab/2; delta++)
			{
				yLocation[slab/2 - delta] = j - delta;
				yLocation[slab/2 + delta] = j + delta;
			}

			sum = 0;
			for(int fi=0; fi < slab; fi++)
			{
				for(int fj=0; fj < slab; fj++)
				{
					index = fi*slab + fj;
					sum += kernel[index] * img.at<uchar>(xLocation[fi], yLocation[fj]);
				}
			}
			retMat.at<uchar>(i ,j) = sum;
		}
	}

	return retMat;
}

// Mean filter: 均值滤波器
Mat filterMean(cv::Mat img, const int slab)
{
	cvtColor(img, img, CV_BGR2GRAY);
	Mat retMat = Mat::zeros(img.rows, img.cols, CV_8UC1);
	for(int i=0; i<img.rows; i++)
		for(int j=0; j<img.cols; j++)
			retMat.at<uchar>(i, j) = img.at<uchar>(i, j);

	int* xLocation = new int[slab];
	int* yLocation = new int[slab];
	float sum;
	int index;

	for(int i= slab/2; i<img.rows - slab/2; i++)
	{
		xLocation[slab/2] = i;
		for(int delta = 0; delta <= slab/2; delta++)
		{
			xLocation[slab/2 - delta] = i - delta;
			xLocation[slab/2 + delta] = i + delta;
		}

		for(int j= slab/2; j<img.cols - slab/2; j++)
		{
			yLocation[slab/2] = j;
			for(int delta = 0; delta <= slab/2; delta++)
			{
				yLocation[slab/2 - delta] = j - delta;
				yLocation[slab/2 + delta] = j + delta;
			}

			sum = 0;
			for(int fi=0; fi < slab; fi++)
			{
				for(int fj=0; fj < slab; fj++)
				{
					index = fi*slab + fj;
					sum += img.at<uchar>(xLocation[fi], yLocation[fj]);
				}
			}
			retMat.at<uchar>(i ,j) = sum/(slab * slab);
		}
	}

	return retMat;
}

中值滤波

// Median filter: 中值滤波器
Mat filterMedian(cv::Mat img, const int slab)
{
	cvtColor(img, img, CV_BGR2GRAY);
	Mat retMat = Mat::zeros(img.rows, img.cols, CV_8UC1);
	for(int i=0; i<img.rows; i++)
		for(int j=0; j<img.cols; j++)
			retMat.at<uchar>(i, j) = img.at<uchar>(i, j);

	int* xLocation	= new int[slab];
	int* yLocation	= new int[slab];
	int* tmpArr		= new int[slab * slab];
	float sum;
	int index;

	for(int i= slab/2; i<img.rows - slab/2; i++)
	{
		xLocation[slab/2] = i;
		for(int delta = 0; delta <= slab/2; delta++)
		{
			xLocation[slab/2 - delta] = i - delta;
			xLocation[slab/2 + delta] = i + delta;
		}

		for(int j= slab/2; j<img.cols - slab/2; j++)
		{
			yLocation[slab/2] = j;
			for(int delta = 0; delta <= slab/2; delta++)
			{
				yLocation[slab/2 - delta] = j - delta;
				yLocation[slab/2 + delta] = j + delta;
			}

			for(int fi = 0; fi<slab; fi++)
				for(int fj =0; fj<slab; fj++)
				{
					index = fi*slab + fj;
					tmpArr[index] = img.at<uchar>(xLocation[fi], yLocation[fj]);
				}

			quickSort(tmpArr, 0, slab*slab - 1);
			retMat.at<uchar>(i ,j) = tmpArr[slab * slab / 2];
		}
	}

	return retMat;
}