数字图像处理（四）图像增强之普通平滑、高斯平滑、laplacian、sobelprewitt锐化

本文参考https://blog.csdn.net/eastmount/article/details/46378783做的一些笔记，解决扫描行填充不为4的整数字节，出现白色斜线的问题，修补了一些bug，注意，

文字算法部分为直接复制粘贴的，代码部分有所该进。

  图像平滑是一种区域增强的算法，平滑算法有邻域平均法、中指滤波、边界保持类滤波等。在图像产生、传输和复制过程中，常常会因为多方面原因而被噪声干扰或出现数据丢失，降低了图像的质量（某一像素，如果它与周围像素点相比有明显的不同，则该点被噪声所感染）。这就需要对图像进行一定的增强处理以减小这些缺陷带来的影响。

       1.简单平滑-邻域平均法
        图像简单平滑是指通过邻域简单平均对图像进行平滑处理的方法，用这种方法在一定程度上消除原始图像中的噪声、降低原始图像对比度的作用。它利用卷积运算对图像邻域的像素灰度进行平均，从而达到减小图像中噪声影响、降低图像对比度的目的。
        但邻域平均值主要缺点是在降低噪声的同时使图像变得模糊，特别在边缘和细节处，而且邻域越大，在去噪能力增强的同时模糊程度越严重。

        2.高斯平滑
        为了克服简单局部平均法的弊端(图像模糊)，目前已提出许多保持边缘、细节的局部平滑算法。它们的出发点都集中在如何选择邻域的大小、形状和方向、参数加平均及邻域各店的权重系数等。
        图像高斯平滑也是邻域平均的思想对图像进行平滑的一种方法，在图像高斯平滑中，对图像进行平均时，不同位置的像素被赋予了不同的权重。
        在图像简单平滑中，算法利用卷积模板逐一处理图像中每个像素，这一过程可以形象地比作对原始图像的像素一一进行过滤整理，在图像处理中把邻域像素逐一处理的算法过程称为滤波器。平滑线性滤波器的工作原理是利用模板对邻域内像素灰度进行加权平均，也称为均值滤波器。
        高斯平滑与简单平滑不同，它在对邻域内像素进行平均时，给予不同位置的像素不同的权值，下图的所示的3*3和5*5领域的高斯模板。

        模板越靠近邻域中心位置，其权值越高。在图像细节进行模糊时，可以更多的保留图像总体的灰度分布特征。下图是常用的四个模板和matlab代码实现：

        代码如下：

[plain] view plain copy

1. I1 = imread('blood1.tif');  
2. I=imnoise(I1,‘salt & pepper’,0.04);                %对图像加椒盐噪声  
3. imshow(I);  
4. h1= [0.1 0.1 0.1; 0.1 0.2 0.1; 0.1 0.1 0.1];        %定义4种模板  
5. h2=1/16.*[1 2 1;2 4 2;1 2 1];  
6. h3=1/8.*[1 1 1;1 0 1;1 1 1];  
7. h4=1/2.*[0 1/4 0;1/4 1 1/4;0 1/4 0];  
8. I2=filter2(h1,I);                                   %用4种模板进行滤波处理  
9. I3=filter2(h2,I);  
10. I4=filter2(h3,I);  
11. I5=filter2(h4,I);  
12. figure,imshow(I2,[])                                %显示处理结果  
13. figure,imshow(I3,[])  
14. figure,imshow(I4,[])  
15. figure,imshow(I5,[])  

        运行效果如下图所示：

        3.中值滤波
        在使用邻域平均法去噪的同时也使得边界变得模糊。而中值滤波是非线性的图像处理方法，在去噪的同时可以兼顾到边界信息的保留。
        选一个含有奇数点的窗口W，将这个窗口在图像上扫描，把窗口中所含的像素点按灰度级的升或降序排列，取位于中间的灰度值来代替该点的灰度值。

        例如选择滤波的窗口如下图，是一个一维的窗口，待处理像素的灰度取这个模板中灰度的中值，滤波过程如下：

        常用的窗口还有方形、十字形、圆形和环形。不同形状的窗口产生不同的滤波效果，方形和圆形窗口适合外轮廓线较长的物体图像，而十字形窗口对有尖顶角状的图像效果好。

        中值滤波对于消除孤立点和线段的干扰十分有用，尤其是对于二进噪声，但对消除高斯噪声的影响效果不佳。对于一些细节较多的复杂图像，可以多次使用不同的中值滤波。matlab实现参考：http://blog.csdn.net/timidsmile/article/details/6904381

        4.边界保持类滤波
        K近邻均值滤波器(KNNF)是指在m*m的窗口中，属于同一集合类的像素，它们的灰度值将高度相关。被处理的像素(对应于窗口中心的像素)可以用窗口内与中心像素灰度最接近的k个近邻像素的平均灰度来替代。步骤如下：
        (1).作一个m*m的作用模板
        (2).在其中选择K个与待处理像素的灰度差为最小的像素
        (3).用这K个像素的灰度均值替换掉原来的值

        在K近旁均值滤波器(KNNMF)中，不选K个邻近像素的平均灰度来替代，而选K个邻近像素的中值灰度来替代，上图中2,3,3中选择3即可。
        下面介绍具体MFC VC++6.0代码实现过程。

三. 图像平滑代码实现

        第一步：在资源视图的Menu中添加子菜单“图像增强”，然后添加“图像平滑”四个选项如下图所示：

        第二步：打开类向导，在ImageProcessingView类中添加相应的四个实现函数：

        第三步：就是具体的平滑实现函数。
        1.普通平滑 模板一
        该算法采用的模板如下：

        代码如下：

void CImageProcessingView::OnTxzqptph()
{
	// TODO: 在此添加命令处理程序代码
	if (numPicture == 0)
	{
		AfxMessageBox("请输入一张图片", MB_OK, 0);
		return;
	}
	AfxMessageBox("普通平滑，均值平滑", MB_OK, 0);
	float H1[3][3] = {
		{0.1,0.1,0.1},
		{0.1,0.2,0.1},
		{0.1,0.1,0.1}
	};
	int HWS = 3;
	FILE *fpo = fopen(BmpName, "rb");
	FILE *fpw = fopen(BmpNameLin, "wb+");
	fread(&bfh, sizeof(BITMAPFILEHEADER), 1, fpo);
	fread(&bih, sizeof(BITMAPINFOHEADER), 1, fpo);
	fwrite(&bfh, sizeof(BITMAPFILEHEADER), 1, fpw);
	fwrite(&bih, sizeof(BITMAPINFOHEADER), 1, fpw);
	fread(m_pImage, m_nImage, 1, fpo);

	unsigned char *ImageSize;//动态内存的分配和释放
	ImageSize = new unsigned char[m_nImage];
	float red, green, blue;
	int x=0, y=0;//代表宽和高的位置
	int LR, LG, LB;//记录实际的位置
	int xx, yy;

	int num;//记录需要填充的字节，因为BMP扫描行字节必须是4的倍数，比如一个像素是3字节，实际需要补1个字节
	if (m_nWidth * 3 % 4 != 0)
	{
		num = 4 - m_nWidth * 3 % 4;
	}
	else if (m_nWidth * 3 % 4 == 0)
	{
		num = 0;
	}

	for (x = HWS/2; x < m_nWidth-HWS/2; x++)//这边我默认原图像边界一个像素不处理，也就是直接给0了，x代表原图第列
	{
		for (y = HWS / 2; y < m_nHeight- HWS / 2; y++)//y代表第几行
		{
			red = green = blue = 0;

			LR = (x + y*m_nWidth) * 3 + y*num;//之所以要加num*y，还是因为扫描行的问题
			//LG = (x + y*m_nWidth) * 3 + 1;
			//LB = (x + y*m_nWidth) * 3 + 2;

			for (int j = 0; j < HWS; j++)//第几行
			{
				yy = y + j - 1;
				for (int i = 0; i < HWS; i++)//第几列
				{
					xx = x + i - 1;
					red += H1[j][i] * (float)(m_pImage[(xx + yy*m_nWidth) * 3 + yy*num]);
					green += H1[j][i] * (float)(m_pImage[(xx + yy*m_nWidth) * 3 + 1+ yy*num]);
					blue+= H1[j][i] * (float)(m_pImage[(xx + yy*m_nWidth) * 3 + 2 + yy*num]);
				}
			}

			ImageSize[LR] = (unsigned char)(red);
			ImageSize[LR + 1] = (unsigned char)(green);
			ImageSize[LR + 2] = (unsigned char)(blue);
		}
	}
	fwrite(ImageSize, m_nImage, 1, fpw);
	/*fwrite(ImageSize, m_nImage, 1, fpw);*/
	fclose(fpo);
	fclose(fpw);
	numPicture = 2;
	level = 400;
	Invalidate();
}

3.高斯平滑
         采用的模板如下：

        代码如下图所示：

void CImageProcessingView::OnGsph()
{
	// TODO: 在此添加命令处理程序代码
	if (numPicture == 0)
	{
		AfxMessageBox("请输入一张图片", MB_OK, 0);
		return;
	}
	AfxMessageBox("普通平滑，均值平滑", MB_OK, 0);
	float H1[3][3] = { 
		{ 1.0 / 16,2.0 / 16,1.0 / 16 }, //模板二:系数1/16 
		{ 2.0 / 16,4.0 / 16,2.0 / 16 },
		{ 1.0 / 16,2.0 / 16,1.0 / 16 } };
	int HWS = 3;
	FILE *fpo = fopen(BmpName, "rb");
	FILE *fpw = fopen(BmpNameLin, "wb+");
	fread(&bfh, sizeof(BITMAPFILEHEADER), 1, fpo);
	fread(&bih, sizeof(BITMAPINFOHEADER), 1, fpo);
	fwrite(&bfh, sizeof(BITMAPFILEHEADER), 1, fpw);
	fwrite(&bih, sizeof(BITMAPINFOHEADER), 1, fpw);
	fread(m_pImage, m_nImage, 1, fpo);

	unsigned char *ImageSize;
	ImageSize = new unsigned char[m_nImage];//new可以有效的内存释放
	float red, green, blue;
	int x = 0, y = 0;//代表宽和高的位置
	int LR, LG, LB;//记录实际的位置
	int xx, yy;
	int num;//记录需要填充的字节
	if (m_nWidth * 3 % 4 != 0)
	{
		num = 4 - m_nWidth * 3 % 4;
	}
	else if (m_nWidth * 3 % 4 == 0)
	{
		num = 0;
	}

	for (x = 0, y = 0; x < m_nWidth; x++)
	{
		LR = (x + y*m_nWidth) * 3;
		ImageSize[LR++] = 0;
		ImageSize[LR++] = 0;
		ImageSize[LR++] = 0;
	}

	for (x = 0, y = 0; y < m_nHeight; y++)
	{
		LR = (x + y*m_nWidth) * 3;
		ImageSize[LR++] = 0;
		ImageSize[LR++] = 0;
		ImageSize[LR++] = 0;
	}

	for (x = HWS / 2; x < m_nWidth - HWS / 2; x++)
	{
		for (y = HWS / 2; y < m_nHeight - HWS / 2; y++)
		{
			red = green = blue = 0;

			LR = (x + y*m_nWidth) * 3 + y*num;
			//LG = (x + y*m_nWidth) * 3 + 1;
			//LB = (x + y*m_nWidth) * 3 + 2;

			for (int j = 0; j < HWS; j++)//第几行
			{
				yy = y + j - 1;
				for (int i = 0; i < HWS; i++)//第几列
				{
					xx = x + i - 1;
					red += H1[j][i] * (float)(m_pImage[(xx + yy*m_nWidth) * 3 + yy * num]);
					green += H1[j][i] * (float)(m_pImage[(xx + yy*m_nWidth) * 3 + 1 + yy * num]);
					blue += H1[j][i] * (float)(m_pImage[(xx + yy*m_nWidth) * 3 + 2 + yy * num]);
				}
			}

			ImageSize[LR] = (unsigned char)(red);
			ImageSize[LR + 1] = (unsigned char)(green);
			ImageSize[LR + 2] = (unsigned char)(blue);
		}
	}
	fwrite(ImageSize, m_nImage, 1, fpw);
	/*fwrite(ImageSize, m_nImage, 1, fpw);*/
	fclose(fpo);
	fclose(fpw);
	numPicture = 2;
	level = 400;
	Invalidate();
}

运行结果如下所示，可以清楚看出斜线没了上面是有斜线的，下面是没有斜线的

 4.中值滤波
        中值滤波我的理解是：它不但可以去除孤点噪声，而且可以保持图像的边缘特性，不会产生显著的模糊；它的方法是把局部区域的像素按灰度等级进行排序，再取该邻域中灰度的中值作为当前像素的灰度值。其步骤如下：
        (1).将滤波模板(含若干个点的滑动窗口)在图像中漫游，并将模板中心与图像中的某个像素位置重合；
        (2).读取模板中各对应像素的灰度值；
        (3).将这些灰度值从小到大排序；
        (4).取这一列数据的中间数据，将其赋值给对应模板中心位置的像素。
        我采用的是3*3的模本，取矩阵中间位置像素替代原像素。代码如下：

void CImageProcessingView::OnZzlb()
{
	// TODO: 在此添加命令处理程序代码
	if (numPicture == 0)
	{
		AfxMessageBox("请输入一张图片", MB_OK, 0);
		return;
	}
	AfxMessageBox("中值滤波", MB_OK, 0);
	FILE *fpo = fopen(BmpName, "rb");
	FILE *fpw = fopen(BmpNameLin, "wb+");
	fread(&bfh, sizeof(BITMAPFILEHEADER), 1, fpo);
	fread(&bih, sizeof(BITMAPINFOHEADER), 1, fpo);
	fwrite(&bfh, sizeof(BITMAPFILEHEADER), 1, fpw);
	fwrite(&bih, sizeof(BITMAPINFOHEADER), 1, fpw);
	fread(m_pImage, m_nImage, 1, fpo);
	unsigned char *ImageSize;
	ImageSize = new unsigned char[m_nImage];
	int HR[9] = { 0 };
	int HG[9] = { 0 };
	int HB[9] = { 0 };
	int x, y,xx,yy;
	int HWS = 3;
	int LR;
	int num;//记录需要填充的字节
	if (m_nWidth * 3 % 4 != 0)
	{
		num = 4 - m_nWidth * 3 % 4;
	}
	else if (m_nWidth * 3 % 4 == 0)
	{
		num = 0;
	}

	for (x = HWS / 2; x <  m_nWidth - HWS / 2; x++)
	{
		for (y = HWS / 2; y < m_nHeight - HWS / 2; y++)
		{
			LR = (x + y*m_nWidth) * 3 + y*num;
			int number = 0;
			for (int j = 0; j < HWS; j++)
			{
				yy = y + j - 1;
				for (int i = 0; i < HWS; i++)
				{
					xx = x + i - 1;
					HR[number] = m_pImage[(xx + yy*m_nWidth)*3 + yy*num];
					HG[number] = m_pImage[(xx + yy*m_nWidth) * 3 + 1 + yy*num];
					HB[number] = m_pImage[(xx + yy*m_nWidth) * 3 + 2 + yy*num];
					number++;
				}
			}

			bubbleSort(HR, 9);//对RGB三个颜色通道排序
			bubbleSort(HG, 9);
			bubbleSort(HB, 9);

			ImageSize[LR] = HR[4];//这里我没有严格按照中值的定义，但是书上写的是H4
			ImageSize[LR+1] = HG[4];
			ImageSize[LR+2] = HB[4];
		}
	}

	fwrite(ImageSize, m_nImage, 1, fpw);
	fclose(fpo);
	fclose(fpw);
	numPicture = 2;
	level = 400;
	Invalidate();
}

四. 图像锐化

       有时还需要加强图像中景物的边缘和轮廓，边缘和轮廓通常位于图像中灰度突出的地方，因而可以直观的想到用灰度的差分对边缘和轮廓进行提取，通常可以通过梯度算子进行提取。图像锐化的目的是提高图像的对比度，从而使图像更清晰，通过提高邻域内像素的灰度差来提高图像的对比度。
       下面介绍图像锐化的几种算子及效果。
        1.拉普拉斯算子(Laplacian)
        拉普拉斯算子是图像邻域内像素灰度差分计算的基础，通过二阶微分推导出的一种图像邻域增强算法。它的基本思想是当邻域的中心像素灰度低于它所在邻域内的其他像素的平均灰度时，此中心像素的灰度应该被进一步降低；当高于时进一步提高中心像素的灰度，从而实现图像锐化处理。
        在算法实现过程中，通过对邻域中心像素的四方向或八方向求梯度，并将梯度和相加来判断中心像素灰度与邻域内其他像素灰度的关系，并用梯度运算的结果对像素灰度进行调整。
        一个连续的二元函数f(x,y)，其拉普拉斯运算定义为：

        对于数字图像，拉普拉斯算子可以简化为：

        也可以表示为卷积的形式：

        其中K=1，I=1时H(r,s)取下式，四方向模板：

        通过模板可以发现，当邻域内像素灰度相同时，模板的卷积运算结果为0；当中心像素灰度高于邻域内其他像素的平均灰度时，模板的卷积运算结果为正数；当中心像素的灰度低于邻域内其他像素的平均灰度时，模板的卷积为负数。对卷积运算的结果用适当的衰弱因子处理并加在原中心像素上，就可以实现图像的锐化处理。
        其中实现过程步骤如下：
        添加子菜单和类向导添加实现函数

        代码如下：

void CImageProcessingView::OnLp()
{
	// TODO: 在此添加命令处理程序代码
	if (numPicture == 0)
	{
		AfxMessageBox("请输入一张图像", MB_OK, 0);
		return;
	}
	AfxMessageBox("拉布拉斯变换", MB_OK, 0);
	int H1[3][3] = {
		{ 0,-1,0}, //模板二:
		{ -1,4,-1},
		{ 0,-1,0 } };
	int HWS = 3;
	FILE *fpo = fopen(BmpName, "rb");
	FILE *fpw = fopen(BmpNameLin, "wb+");
	fread(&bfh, sizeof(BITMAPFILEHEADER), 1, fpo);
	fread(&bih, sizeof(BITMAPINFOHEADER), 1, fpo);
	fwrite(&bfh, sizeof(BITMAPFILEHEADER), 1, fpw);
	fwrite(&bih, sizeof(BITMAPINFOHEADER), 1, fpw);
	fread(m_pImage, m_nImage, 1, fpo);

	unsigned char *ImageSize;
	ImageSize = new unsigned char[m_nImage];
	int red, green, blue;
	int x = 0, y = 0;//代表宽和高的位置
	int LR, LG, LB;//记录实际的位置
	int xx, yy; 
	int num;//记录需要填充的字节
	if (m_nWidth * 3 % 4 != 0)
	{
		num = 4 - m_nWidth * 3 % 4;
	}
	else if (m_nWidth * 3 % 4 == 0)
	{
		num = 0;

	}

	for (x = HWS / 2; x < m_nWidth - HWS / 2; x++)
	{
		for (y = HWS / 2; y < m_nHeight - HWS / 2; y++)
		{
			red = green = blue = 0;

			LR = (x + y*m_nWidth) * 3 + y*num;
			//LG = (x + y*m_nWidth) * 3 + 1;
			//LB = (x + y*m_nWidth) * 3 + 2;

			for (int j = 0; j < HWS; j++)//第几行这两个循环是求权重后的均值
			{
				yy = y + j - 1;
				for (int i = 0; i < HWS; i++)//第几列
				{
					xx = x + i - 1;
					red += H1[j][i] * (m_pImage[(xx + yy*m_nWidth) * 3 + yy*num]);
					green += H1[j][i] * (m_pImage[(xx + yy*m_nWidth) * 3 + 1 + yy*num]);
					blue += H1[j][i] * (m_pImage[(xx + yy*m_nWidth) * 3 + 2 + yy*num]);
				}
			}
			if (red < 0 || red >= 256) red = 0;
			ImageSize[LR] = (unsigned char)(red);
			if (green < 0 || green >= 256) green= 0;
			ImageSize[LR + 1] = (unsigned char)(green);
			if (blue < 0 || blue >= 256) blue = 0;
			ImageSize[LR + 2] = (unsigned char)(blue);
		}
	}
	fwrite(ImageSize, m_nImage, 1, fpw);
	fclose(fpo);
	fclose(fpw);
	numPicture = 2;
	level = 400;
	Invalidate();
}

  运行效果如下图所示：

        2.高通滤波
        常用的高通模板如下所示，其中H2有的书又称为拉普拉斯八方向的锐化模板。

    3.Sobel算子

        代码如下所示，需要注意X和Y两个方向的模板处理：

void CImageProcessingView::OnSobel()
{
	// TODO: 在此添加命令处理程序代码
	if (numPicture == 0)
	{
		AfxMessageBox("Plaece Enter a picture!", MB_OK, 0);
		return;
	}
	AfxMessageBox("Sobel in sessing", MB_OK, 0);
	int HX[3][3] = {{-1,-2,-1},{0,0,0,},{1,2,1}};
	int HY[3][3] = { {-1,0,1},{-2,0,2},{-1,0,1} };

	int HWS = 3;

	int num;//记录需要填充的字节
	if (m_nWidth * 3 % 4 != 0)
	{
		num = 4 - m_nWidth * 3 % 4;
	}
	else if (m_nWidth * 3 % 4 == 0)
	{
		num = 0;

	}

	FILE *fpo = fopen(BmpName, "rb");
	FILE *fpw = fopen(BmpNameLin, "wb+");
	fread(&bfh, sizeof(BITMAPFILEHEADER), 1, fpo);
	fread(&bih, sizeof(BITMAPINFOHEADER), 1, fpo);
	fwrite(&bfh, sizeof(BITMAPFILEHEADER), 1, fpw);
	fwrite(&bih, sizeof(BITMAPINFOHEADER), 1, fpw);
	fread(m_pImage, m_nImage, 1, fpo);


	unsigned char *ImageSize;
	ImageSize = new unsigned char[m_nImage];
	int x, y,xx,yy;
	int redX, greenX, blueX;
	int redY, greenY, blueY;
	int R, G, B;
	int LR;
	unsigned char other = 0;        //补码00H='\0'  



	for (x = HWS / 2; x < m_nWidth - HWS / 2; x++)
	{
		for (y = HWS / 2; y < m_nHeight - HWS / 2; y++)
		{
			redX = greenX = blueX = 0;
			redY = greenY = blueY = 0;
			LR = (y*m_nWidth + x)*3 + y*num;

			for (int j = 0; j < HWS; j++)//第几行这两个循环是求权重后的均值
			{
				yy = y + j - 1;
				for (int i = 0; i < HWS; i++)//第几列
				{
					xx = x + i - 1;
					redX += HX[j][i] * (m_pImage[(xx + yy*m_nWidth) * 3 + yy*num]);
					greenX += HX[j][i] * (m_pImage[(xx + yy*m_nWidth) * 3 + 1 + yy*num]);
					blueX += HX[j][i] * (m_pImage[(xx + yy*m_nWidth) * 3 + 2 + yy*num]);

					redY += HY[j][i] * (m_pImage[(xx + yy*m_nWidth) * 3 + yy*num]);
					greenY += HY[j][i] * (m_pImage[(xx + yy*m_nWidth) * 3 + 1 + yy*num]);
					blueY += HY[j][i] * (m_pImage[(xx + yy*m_nWidth) * 3 + 2 + yy*num]);

				}
			}


			

			R = (int)sqrt(redX*redX*1.0 + redY*redY*1.0);
			G = (int)sqrt(greenX*greenX*1.0 + greenY*greenY*1.0);
			B = (int)sqrt(blueX*blueX*1.0 + blueY*blueY*1.0);
			if (R > 255) R = 255;
			ImageSize[LR] = unsigned char(R);
			if (G > 255)G = 255;
			ImageSize[LR+1] = unsigned char(G);
			if (B > 255)B = 255;
			ImageSize[LR+2] = unsigned char(B);
		}
	}

	fwrite(ImageSize, m_nImage, 1, fpw);
	fclose(fpo);
	fclose(fpw);
	numPicture = 2;
	level = 401;
	Invalidate();

}

 4.Isotropic算子

  5.Prewitt算子