基于局部均方差的图像去噪声算法

1.基于局部均方差去噪原理

假设一幅M*N的灰度图像的灰度图像，x(i,j)为模板中心的灰度值，那么在（2*n+1）（2*m+1）的窗口内，其局部均值可以表示为：

基于局部均方差的图像去噪声算法

其局部均方差可以表示为：

基于局部均方差的图像去噪声算法

加性去噪后的结果为：

基于局部均方差的图像去噪声算法

其中K为：

基于局部均方差的图像去噪声算法

关于K计算的公式中，基于局部均方差的图像去噪声算法为用户输入的参数。根据我之前的博客《图像比较之模板匹配》中，提到的方差和均值的关系，我们可以推导出，方差可以表示为：

基于局部均方差的图像去噪声算法

我们知道，方差在统计学中表示的是与中心偏离的程度，用来衡量数据的波动性大小。对于图像而言，当上述局部方差比较小，意味着图像中该局部区域属于灰度平坦区，各个像素灰度值相差不大；相反，当上述局部方差比较大的时候，意味着图像中该局部区域属于边缘或者是其他高频部分区域，各个像素的灰度值相差比较大。

当局部属于平坦区时，方差很小，趋近于0。该点滤波之后的像素就是该点的局部平均值。由于该局部各点像素的灰度值相差不大，其局部平均值也与各个像素的灰度值相差不大；当局部属于边缘区域时，方差较大，相对于用户输入的参数可以基本忽略不计，其图像去噪之后，就等于输入的图像灰度值。可以说，这种方法在一定程度上对边缘具有保留，能够保留边缘的同时，进行去噪。

类似边缘保留的去噪算法，诸如双边滤波等。都是对图像中的平坦区域和边缘区域进行相应的数学描述，再基于该数学描述所对应的图像区域，对不同区域赋予不同的权重，边缘区域权重适当较小；平坦区域权重适当增大；做不同程度的去噪。

2.代码实现

8Bit BMP图像数据格式的关键代码片段：

	for (x = SIZE/2;x < height-SIZE/2;x++)
	{
		for (y = SIZE/2;y < width-SIZE/2;y++)
		{
			Sum = 0;
			Mean = 0;
			Sum2 = 0;
			index = x*lineByte+y;

			for(i = -SIZE/2;i<= SIZE/2;i++)
			{
				for(j = -SIZE/2;j<=SIZE/2;j++)
				{
					tmp = img_data[index+i*lineByte+j];
					Sum += tmp;
					Sum2 += img_data[index+i*lineByte+j]*img_data[index+i*lineByte+j];
				}
			}
				Mean = Sum/size;
				Var = (Sum2-(Sum*Sum)/size)/size;
				k = Var/(Var+Sigmma);

				img_data[index] = CLIP255((int)((1-k)*Mean+k*img_data[index]));
		}
	}

24Bit BMP图像数据格式的关键代码片段：

        for(i = SIZE/2;i < height-SIZE/2;i++)
	{
		for(j = SIZE/2;j < width-SIZE/2;j++)
		{
			index = i*lineByte+3*j;
			SumR2 = 0;
			SumG2 = 0;
			SumB2 = 0;
			SumR = 0;
			SumG = 0;
			SumB = 0;
			MeanR = 0;
			VarR = 0;
			MeanG = 0;
			VarG = 0;
			MeanB = 0;
			VarB = 0;
			       
			for(m = -SIZE/2;m<= SIZE/2;m++)
			{
				for(n = -SIZE/2;n<=SIZE/2;n++)
				{
					SumB  += img_data[index+m*lineByte+n*3];
					SumB2 += img_data[index+m*lineByte+n*3]*img_data[index+m*lineByte+n*3];

					SumG  += img_data[index+1+m*lineByte+n*3];
					SumG2 += img_data[index+1+m*lineByte+n*3]*img_data[index+1+m*lineByte+n*3];

					SumR  += img_data[index+2+m*lineByte+n*3];
					SumR2 += img_data[index+2+m*lineByte+n*3]*img_data[index+2+m*lineByte+n*3];
				}
			}		
			               //计算平均值
					MeanR = SumR/size;
					MeanG = SumG/size;
					MeanB = SumB/size;

					//计算方差
					VarR = (SumR2-(SumR*SumR)/size)/size;
					VarG = (SumG2-(SumG*SumG)/size)/size;
					VarB = (SumB2-(SumB*SumB)/size)/size;
	

					//计算系数
					kr = VarR /(VarR+Sigmma);
					kg = VarG /(VarG+Sigmma);
					kb = VarB /(VarB+Sigmma);
					
					
					img_data[index+2] =  CLIP255((int)((1-kr)*MeanR+kr*img_data[index+2]));
					img_data[index+1] =  CLIP255((int)((1-kg)*MeanG+kg*img_data[index+1]));
					img_data[index]   =  CLIP255((int)((1-kb)*MeanB+kb*img_data[index]));

		}
		
	}