C++ 实现matlab高斯滤波函数imgaussfilt

 最近需要用C++来实现matlab中的空间高斯滤波函数imgaussfilt，做下笔记。

/*
* 高斯滤波
*/
void SeparateGaussianFilt(const std::vector<std::vector<double>>& src, std::vector<std::vector<double>>& dst,const double& sigma)
{
	//创建一维高斯内核
	std::vector<double> kernel;
	double filterSize = 2 * ceil(2 * sigma) + 1; //高斯内核宽度和高度，与matlab保持一致
	double filterRadius = (filterSize - 1) / 2;  //高斯内核半径
	double sum = 0.0;
	
	for (int i = 0; i < filterSize; ++i)
	{
		double g = exp(-(pow((i - filterRadius), 2)) / (2 * sigma * sigma));
		kernel.push_back(g);
		sum += g;
	}
	//高斯内核归一化
	for (int i = 0; i < kernel.size(); i++)
	{
		kernel[i] = double(kernel[i] / sum);
	}

	//数据边界复制填充
	int a = src.size() + (2 * filterRadius);
	int b = src[0].size() + (2 * filterRadius);
	std::vector<std::vector<double>> newSrc(a);
	for (int i = 0; i < a; i++)
	{
		newSrc[i].resize(b);
	}
	
	// 边界行填充
	for (int i = 0; i < filterRadius; i++)
	{
		for (int j = 0; j < (b - 2 * filterRadius); j++)
		{
			newSrc[i][filterRadius + j] = src[0][j];
		}
	}
	for (int i = (a - filterRadius); i < a; i++)
	{
		for (int j = 0; j < (b - 2 * filterRadius); j++)
		{
			newSrc[i][filterRadius + j] = src[a - (2 * filterRadius) - 1][j];
		}
	}

	// 中间部分填充
	for (int i = filterRadius; i < a - filterRadius; i++)
	{
		for (int j = filterRadius; j < b - filterRadius; j++)
		{
			newSrc[i][j] = src[i - filterRadius][j - filterRadius];
		}
	}

	// 边界列填充
	for (int i = 0; i < a; i++)
	{
		for (int j = 0; j < filterRadius; j++)
		{
			newSrc[i][j] = newSrc[i][filterRadius];
		}
	}
	for (int i = 0; i < a; i++)
	{
		for (int j = (b - filterRadius); j < b; j++)
		{
			newSrc[i][j] = newSrc[i][b - filterRadius - 1];
		}
	}

	// 高斯滤波水平方向
	for (int i = filterRadius; i < src.size() + filterRadius; i++)
	{
		for (int j = filterRadius; j < src[0].size() + filterRadius; j++)
		{
			double sum = 0.0;
			for (int r = (-filterRadius); r <= filterRadius; r++)
			{
				sum = sum + newSrc[i][j + r] * kernel[r + filterRadius];
			}
			/*if (sum < 0)
			{
				sum = 0;
			}
			else if (sum > 255)
			{ 
				sum = 255;
			}*/
			dst[i - filterRadius][j - filterRadius] = sum;
		}
	}

	//水平方向处理后的dst边界复制填充，边界行填充
	for (int i = 0; i < filterRadius; i++)
	{
		for (int j = 0; j < (b - 2 * filterRadius); j++)
		{
			newSrc[i][filterRadius + j] = dst[0][j];
		}
	}
	for (int i = (a - filterRadius); i < a; i++)
	{
		for (int j = 0; j < (b - 2 * filterRadius); j++)
		{
			newSrc[i][filterRadius + j] = dst[a - (2 * filterRadius) - 1][j];
		}
	}

	// 中间部分填充
	for (int i = filterRadius; i < a - filterRadius; i++)
	{
		for (int j = filterRadius; j < b - filterRadius; j++)
		{
			newSrc[i][j] = dst[i - filterRadius][j - filterRadius];
		}
	}

	// 边界列填充
	for (int i = 0; i < a; i++)
	{
		for (int j = 0; j < filterRadius; j++)
		{
			newSrc[i][j] = newSrc[i][filterRadius];
		}
	}
	for (int i = 0; i < a; i++)
	{
		for (int j = (b - filterRadius); j < b; j++)
		{
			newSrc[i][j] = newSrc[i][b - filterRadius - 1];
		}
	}

	//高斯滤波垂直方向
	for (int i = filterRadius; i < src.size() + filterRadius; i++)
	{
		for (int j = filterRadius; j < src[0].size() + filterRadius; j++)
		{
			double sum = 0.0;
			for (int r = (-filterRadius); r <= filterRadius; r++)
			{
				sum = sum + newSrc[i + r][j] * kernel[r + filterRadius];
			}
			/*if (sum < 0)
			{
				sum = 0;
			}
			else if (sum > 255)
			{
				sum = 255;
			}*/
			dst[i - filterRadius][j - filterRadius] = sum;
		}
	}
}