先举个例子说明一下什么是平滑(smoothing),如下面两幅图所示:可以看到,图3.2比图3.1柔和一些(也模糊一些)。是不是觉得很神奇?其实实现起来很简单。我们将原图中的每一点的灰度和它周围八个点的灰度相加,然后除以9,作为新图中对应点的灰度,就能实现上面的效果。
图3.1 原图
|
图3.2 经过平滑处理后的图
|
这么做并非瞎蒙,而是有其道理的。大概想一想,也很容易明白。举个例子,就象和面一样,先在中间加点水,然后不断把周围的面和进来,搅拌几次,面就均匀了。
用信号处理的理论来解释,这种做法实现的是一种简单的低通滤波器(low pass filter)。哇,好深奥呀!不要紧,这些理论的内容并不多,而且知道一些理论也是很有好处的。在灰度连续变化的图象中,如果出现了与相邻象素的灰度相差很大的点,比如说一片暗区中突然出现了一个亮点,人眼能很容易觉察到。就象看老电影时,由于胶片太旧,屏幕上经常会出现一些亮斑。这种情况被认为是一种噪声。灰度突变在频域中代表了一种高频分量,低通滤波器的作用就是滤掉高频分量,从而达到减少图象噪声的目的。
为了方便地叙述上面所说的“将原图中的每一点的灰度和它周围八个点的灰度相加,然后除以9,作为新图中对应点的灰度”这一操作,我们采用如下的表示方法:
(3.1)
这种表示方法有点象矩阵,我们称其为模板(template)。中间的黑点表示中心元素,即,用哪个元素做为处理后的元素。例如[2. 1]表示将自身的2倍加上右边的元素作为新值,而[2 1.]表示将自身加上左边元素的2倍作为新值。
通常,模板不允许移出边界,所以结果图象会比原图小,例如模板是
,原图是
,经过模板操作后的图象为
;其中数字代表灰度,x表示边界上无法进行模板操作的点,通常的做法是复制原图的灰度,不进行任何处理。
模板操作实现了一种邻域运算(Neighborhood Operation),即某个象素点的结果灰度不仅和该象素灰度有关,而且和其邻域点的值有关。在以后介绍的细化算法中,我们还将接触到邻域运算。模板运算的数学涵义是一种卷积(或互相关)运算,你不需要知道卷积的确切含义,只要有这么一个概念就可以了。
模板运算在图象处理中经常要用到,可以看出,它是一项非常耗时的运算。以
(3.2)
为例,每个象素完成一次模板操作要用9个乘法、8个加法、1个除法。对于一幅n×n(宽度×高度)的图象,就是9n2个乘法,8n2个加法和n2个除法,算法复杂度为O(n2),这对于大图象来说,是非常可怕的。所以,一般常用的模板并不大,如3×3,4×4。有很多专用的图象处理系统,用硬件来完成模板运算,大大提高了速度。另外,可以设法将二维模板运算转换成一维模板运算,对速度的提高也是非常可观的。例如,(3.2)式可以分解成一个水平模板和一个垂直模板,即,
=
×
=
(3.3)
我们来验证一下。
设图象为
,经过(3.2)式处理后变为
,经过(3.3)式处理后变为
,两者完全一样。如果计算时不考虑周围一圈的象素,前者做了4×(9个乘法,8个加法,1个除法),共36个乘法,32个加法,4个除法;后者做了4×(3个乘法,2个加法)+4×(3个乘法,2个加法)+4个除法,共24个乘法,16个加法,4个除法,运算简化了不少,如果是大图,效率的提高将是非常客观的。
平滑模板的思想是通过将一点和周围8个点作平均,从而去除突然变化的点,滤掉噪声,其代价是图象有一定程度的模糊。上面提到的模板(3.1),就是一种平滑模板,称之为Box模板。Box模板虽然考虑了邻域点的作用,但并没有考虑各点位置的影响,对于所有的9个点都一视同仁,所以平滑的效果并不理想。实际上我们可以想象,离某点越近的点对该点的影响应该越大,为此,我们引入了加权系数,将原来的模板改造成
,可以看出,距离越近的点,加权系数越大。
新的模板也是一个常用的平滑模板,称为高斯(Gauss)模板。为什么叫这个名字,这是因为这个模板是通过采样2维高斯函数得到的。
设图象为
,分别用两种平滑模板处理(周围一圈象素直接从原图拷贝)。采用Box模板的结果为
,采用高斯模板的结果为
。
可以看到,原图中出现噪声的区域是第2行第2列和第3行第2列,灰度从2一下子跳到了6,用Box模板处理后,灰度从3.11跳到4.33;用高斯模板处理后,灰度从3.跳到4.56,都缓和了跳变的幅度,从这一点上看,两者都达到了平滑的目的。但是,原图中的第3,第4行总的来说,灰度值是比较高的,经模板1处理后,第3行第2列元素的灰度变成了4.33,与第3,第4行的总体灰度相比偏小,另外,原图中第3行第2列元素的灰度为6,第3行第3列元素的灰度为4,变换后,后者4.56反而比前者4.33大了。而采用高斯模板没有出现这些问题,究其原因,就是因为它考虑了位置的影响。
举个实际的例子:下图中,从左到右分别是原图,用高斯模板处理的图,用Box模板处理的图,可以看出,采用高斯模板,在实现平滑效果的同时,要比Box模板清晰一些。
在学习锐化后,我们将给出一个通用的3×3模板操作的程序。
图3.3 高斯模板和Box模板的对比图
中值滤波也是一种典型的低通滤波器,它的目的是保护图象边缘的同时去除噪声。所谓中值滤波,是指把以某点(x,y)为中心的小窗口内的所有象素的灰度按从大到小的顺序排列,将中间值作为(x,y)处的灰度值(若窗口中有偶数个象素,则取两个中间值的平均)。中值滤波是如何去除噪声的呢?举个例子就很容易明白了。
原图
|
处理后的图
|
|
|
图中数字代表该处的灰度。可以看出原图中间的6和周围的灰度相差很大,是一个噪声点。经过3×1窗口(即水平3个象素取中间值)的中值滤波,得到右边那幅图,可以看出,噪声点被去除了。
下面将中值滤波和上面介绍的两种平滑模板作个比较,看看中值滤波有什么特点。我们以一维模板为例,只考虑水平方向,大小为3×1(宽×高)。Box模板为
,高斯模板为
。
先考察第一幅图:
原图
|
经Box模板处理后
|
经Gauss模板处理后
|
经中值滤波处理后
|
|
|
|
|
从原图中不难看出左边区域灰度值低,右边区域灰度值高,中间有一条明显的边界,这一类图象称之为“step”(就象灰度上了个台阶)。应用平滑模板后,图象平滑了,但是也使边界模糊了。应用中值滤波,就能很好地保持原来的边界。所以说,中值滤波的特点是保护图象边缘的同时去除噪声。
再看第二幅图:
原图
|
经Box模板处理后
|
经Gauss模板处理后
|
经中值滤波处理后
|
|
|
|
|
不难看出,原图中有很多噪声点(灰度为正代表灰度值高的点,灰度为负代表灰度值低的点),而且是杂乱无章,随机分布的。这也是一类很典型的图,称之为高斯噪声。经过Box平滑,噪声的程度有所下降。Gauss模板对付高斯噪声非常有效。而中值滤波对于高斯噪声则无能为力。
最后看第三幅图:
原图
|
经Box模板处理后
|
经Gauss模板处理后
|
经中值滤波处理后
|
|
|
|
|
从原图中不难看出,中间的灰度要比两边高许多。这也是一类很典型的图,称之为脉冲 (impulse)。可见,中值滤波对脉冲噪声非常有效。
综合以上三类图,不难得出下面的结论:中值滤波容易去除孤立点,线的噪声同时保持图象的边缘;它能很好的去除二值噪声,但对高斯噪声无能为力。要注意的是,当窗口内噪声点的个数大于窗口宽度的一半时,中值滤波的效果不好。这是很显然的。
下面的程序实现了中值滤波,参数Hori是一个布尔变量,若为真,做水平中值滤波,否则,做垂直中值滤波。
BOOL MedianFilter(HWND hWnd,BOOL Hori)



{

DWORD OffBits,BufSize;

LPBITMAPINFOHEADER lpImgData;

LPSTR lpPtr;

HLOCAL hTempImgData;

LPBITMAPINFOHEADER lpTempImgData;

LPSTR lpTempPtr;

HDC hDc;

HFILE hf;

LONG x,y;

int g,g1,g2,g3;

//OffBits为BITMAPINFOHEADER结构长度加调色板的大小

OffBits=bf.bfOffBits-sizeof(BITMAPFILEHEADER);

BufSize=OffBits+bi.biHeight*LineBytes;//要开的缓冲区的大小

if((hTempImgData=LocalAlloc(LHND,BufSize))==NULL)



{

MessageBox(hWnd,"Error alloc memory!","Error Message",MB_OK|

MB_ICONEXCLAMATION);

return FALSE;

}

lpImgData=(LPBITMAPINFOHEADER)GlobalLock(hImgData);

lpTempImgData=(LPBITMAPINFOHEADER)LocalLock(hTempImgData);

//拷贝头信息及位图数据

memcpy(lpTempImgData,lpImgData,BufSize);

//注意边界点不处理,所以y从1到高度-2,x类似

for(y=1;y<bi.biHeight-1;y++)


for(x=1;x<bi.biWidth-1;x++)
{

lpPtr=(char *)lpImgData+(BufSize-LineBytes-y*LineBytes)+x;

lpTempPtr=(char*)lpTempImgData+

(BufSize-LineBytes-y*LineBytes)+x;

g2=(unsigned char)*(lpPtr);


if(Hori)
{ //水平方向

g1=(unsigned char)*(lpPtr-1); //左邻点

g3=(unsigned char)*(lpPtr+1); //右邻点

}


else
{ //垂直方向

g1=(unsigned char)*(lpPtr+LineBytes); //上邻点

g3=(unsigned char)*(lpPtr-LineBytes); //下邻点

}

//三者取中


if(g1>g2)
{

if(g2>g3) g=g2;


else
{

if(g1>g3) g=g3;

else g=g1;

}

}


else
{ //g1<=g2

if(g1>g3) g=g1;


else
{

if(g2>g3) g=g3;

else g=g2;

}

}

*lpTempPtr=(BYTE)g; //存入新的缓冲区内

}

hDc=GetDC(hWnd);

if(hBitmap!=NULL)

DeleteObject(hBitmap);

//产生新的位图

hBitmap=CreateDIBitmap(hDc,(LPBITMAPINFOHEADER)lpTempImgData,

(LONG)CBM_INIT,

(LPSTR)lpTempImgData+

sizeof(BITMAPINFOHEADER)+

NumColors*sizeof(RGBQUAD),

(LPBITMAPINFO)lpTempImgData,

DIB_RGB_COLORS);

if(Hori) //取不同的结果文件名

hf=_lcreat("c:\\hmedian.bmp",0);

else

hf=_lcreat("c:\\vmedian.bmp",0);

_lwrite(hf,(LPSTR)&bf,sizeof(BITMAPFILEHEADER));

_lwrite(hf,(LPSTR)lpTempImgData,BufSize);

_lclose(hf);

//释放内存及资源

ReleaseDC(hWnd,hDc);

LocalUnlock(hTempImgData);

LocalFree(hTempImgData);

GlobalUnlock(hImgData);

return TRUE;

}


锐化(sharpening)和平滑恰恰相反,它是通过增强高频分量来减少图象中的模糊,因此又称为高通滤波(high pass filter)。锐化处理在增强图象边缘的同时增加了图象的噪声。
常用的锐化模板是拉普拉斯(Laplacian)模板(见(3.4)式),又是个数学家的名字,可见学好数学,走遍天下都不怕。
(3.4)
容易看出拉普拉斯模板的作法:先将自身与周围的8个象素相减,表示自身与周围象素的差别;再将这个差别加上自身作为新象素的灰度。可见,如果一片暗区出现了一个亮点,那么锐化处理的结果是这个亮点变得更亮,增加了图象的噪声。
因为图象中的边缘就是那些灰度发生跳变的区域,所以锐化模板在边缘检测中很有用,这一点将在后面详细介绍。
图3.1经过拉普拉斯模板处理后,如图3.4所示
图3.4 锐化
下面给出的程序是一个通用的3×3模板的函数,其中第二参数为模板类型,为如下定义的常量:

#define TEMPLATE_SMOOTH_BOX 1 //Box平滑模板


#define TEMPLATE_SMOOTH_GAUSS 2 //高斯平滑模板


#define TEMPLATE_SHARPEN_LAPLACIAN 3 //拉普拉斯锐化模板


对应的模板数组如下



int Template_Smooth_Box[9]=

{1,1,1,1,1,1,1,1,1};



int Template_Smooth_Gauss[9]=

{1,2,1,2,4,2,1,2,1};



int Template_Sharpen_Laplacian[9]=

{-1,-1,-1,-1,9,-1,-1,-1,-1};


以后我们碰到其它的模板,仍然要用这个函数,所做的操作只是增加一个常量标识,及其对应的模板数组。


要注意的是,运算后如果出现了大于255或者小于0的点,称为溢出,溢出点的处理通常是截断,即大于255时,令其等于255;小于0时,取其绝对值。


这段程序和前几章介绍的代码许多地方是很相似的,所以注释简单一些。程序中并没有用到那种分解成两个一维模板的快速算法,你如果有兴趣,可以自己编着试试。


BOOL TemplateOperation(HWND hWnd, int TemplateType)




{


DWORD OffBits,BufSize;


LPBITMAPINFOHEADER lpImgData;


LPSTR lpPtr;


HLOCAL hTempImgData;


LPBITMAPINFOHEADER lpTempImgData;


LPSTR lpTempPtr;


HDC hDc;


HFILE hf;


LONG x,y;


float coef; //模板前面所乘的系数


int CoefArray[9]; //模板数组


float TempNum;


char filename[80];



switch(TemplateType)

{ //判断模板类型


case TEMPLATE_SMOOTH_BOX: //Box平滑模板


coef=(float)(1.0/9.0);


memcpy(CoefArray,Template_Smooth_Box,9*sizeof(int));


strcpy(filename,"c:\\smbox.bmp");


break;


case TEMPLATE_SMOOTH_GAUSS: //高斯平滑模板


coef=(float)(1.0/16.0);


memcpy(CoefArray,Template_Smooth_Gauss,9*sizeof(int));


strcpy(filename,"c:\\smgauss.bmp");


break;


case TEMPLATE_SHARPEN_LAPLACIAN: //拉普拉斯锐化模板


coef=(float)1.0;


memcpy(CoefArray,Template_Sharpen_Laplacian,9*sizeof(int));


strcpy(filename,"c:\\shlaplac.bmp");


break;


}


OffBits=bf.bfOffBits-sizeof(BITMAPFILEHEADER);


BufSize=OffBits+bi.biHeight*LineBytes;


if((hTempImgData=LocalAlloc(LHND,BufSize))==NULL)




{


MessageBox(hWnd,"Error alloc memory!","Error Message",MB_OK|


MB_ICONEXCLAMATION);


return FALSE;


}


lpImgData=(LPBITMAPINFOHEADER)GlobalLock(hImgData);


lpTempImgData=(LPBITMAPINFOHEADER)LocalLock(hTempImgData);


lpPtr=(char *)lpImgData;


lpTempPtr=(char *)lpTempImgData;


//先将原图直接拷贝过来,其实主要是拷贝周围一圈的象素


memcpy(lpTempPtr,lpPtr,BufSize);


for(y=1;y<bi.biHeight-1;y++) //注意y的范围是从1到bi.biHeight-2



for(x=1;x<bi.biWidth-1;x++)

{ //注意x的范围是从1到bi.biWidth-2


lpPtr=(char *)lpImgData+(BufSize-LineBytes-y*LineBytes)+x;


lpTempPtr=(char*)lpTempImgData+


(BufSize-LineBytes-y*LineBytes)+x;


TempNum=(float)((unsigned char)*(lpPtr+LineBytes-1))*


CoefArray[0];


TempNum+=(float)((unsigned char)*(lpPtr+LineBytes))*


CoefArray[1];


TempNum+=(float)((unsigned char)*(lpPtr+LineBytes+1))*


CoefArray[2];


TempNum+=(float)((unsigned char)*(lpPtr-1))*CoefArray[3];


TempNum+=(float)((unsigned char)*lpPtr)*CoefArray[4];


TempNum+=(float)((unsigned char)*(lpPtr+1))*CoefArray[5];


TempNum+=(float)((unsigned char)*(lpPtr-LineBytes-1))*


CoefArray[6];


TempNum+=(float)((unsigned char)*(lpPtr-LineBytes))*


CoefArray[7];


TempNum+=(float)((unsigned char)*(lpPtr-LineBytes+1))*


CoefArray[8];


//最后乘以系数


TempNum*=coef;


//注意对溢出点的处理


if(TempNum>255.0) *lpTempPtr=(BYTE)255;


else if(TempNum<0.0)


*lpTempPtr=(unsigned char)fabs(TempNum);


else *lpTempPtr=(BYTE)TempNum;


}


hDc=GetDC(hWnd);


if(hBitmap!=NULL)


DeleteObject(hBitmap);


hBitmap=CreateDIBitmap(hDc,(LPBITMAPINFOHEADER)lpTempImgData,


(LONG)CBM_INIT,


(LPSTR)lpTempImgData+


sizeof(BITMAPINFOHEADER)+


NumColors*sizeof(RGBQUAD),


(LPBITMAPINFO)lpTempImgData,


DIB_RGB_COLORS);


hf=_lcreat(filename,0);


_lwrite(hf,(LPSTR)&bf,sizeof(BITMAPFILEHEADER));


_lwrite(hf,(LPSTR)lpTempImgData,BufSize);


_lclose(hf);


ReleaseDC(hWnd,hDc);


LocalUnlock(hTempImgData);


LocalFree(hTempImgData);


GlobalUnlock(hImgData);


return TRUE;


}
