第3章空域图像增强

图像增强的目的是根据主观视觉感受或者满足某些特定分析的需要来改善图像的质量，有选择地强调图像中感兴趣的整体或局部特征或者抑制图像中不需要的特征，使图像更适合特定的应用。当以人类视觉解释为目的时，依据观察者的主观判断作为图像质量评价标准； 当以机器识别为目的时，依据机器识别的准确率作为图像质量评价标准。图像增强是一个主观过程，不存在通用的数学模型。图像增强方法主要分为两类： 空域方法和频域方法。空域指图像空间本身，这类方法是对图像的像素直接进行处理。频域指图像的傅里叶变换域，这类方法通过离散傅里叶变换将图像从空域转换到频域，利用频域滤波器对图像的不同频率成分进行处理。本章讨论空域图像增强方法，第4章和第5章将分别讨论频域变换和频域增强方法。
3.1背景知识
空域图像增强是指直接对图像的像素本身进行操作，大致可分为两类： 点处理和邻域处理。对于点处理，在像素(x,y)处g(x,y)的值仅取决于f(x,y)的值，可表示为


g(x,y)=E［f(x,y)］(31)



式中，E［·］表示增强操作，f(x,y)为输入图像，g(x,y)为增强图像。设r表示输入图像的灰度级，s表示增强图像的灰度级，点处理操作即为灰度级映射，可表示为


s=T(r)(32)



式中，T(r)为变换函数。
式(32)表明，对图像中任意像素的增强仅依赖于该像素的灰度级，因此这类方法称为点处理。点处理又包括灰度级变换和直方图处理。
对于邻域处理，像素(x,y)的邻域定义为中心位于像素(x,y)处的局部区域，例如，图31显示了中心在坐标(x,y)处的3×3正方形邻域，g(x,y)的值不仅取决于f(x,y)的值，而且取决于其邻域像素的值。对于图像中的每一个像素(x,y)，增强操作是以像素(x,y)为中心的邻域内所有像素值的函数： 


g(x,y)=E{f(s,t)},(s,t)∈Sxy(33)



式中，Sxy表示以像素(x,y)为中心的邻域像素集合。
式(33)表明，对图像中任意像素的增强，作用于该像素邻域内的所有像素，因此这类方法称为邻域处理。空域滤波方法属于邻域处理，包括空域图像平滑和空域图像锐化等。



图31邻域处理示意图



3.2直方图概念
直方图是数字图像处理中一个重要的基础工具。在讨论各种空域处理技术之前，首先需要清楚直方图的概念和表示意义。直方图提供了图像的统计信息，为理解多种空域增强技术的内涵提供了铺垫。直方图操作能够直接用于图像增强，见3.4节的直方图图像增强技术。
3.2.1直方图
直方图是数字图像的统计表征量。对于一幅灰度图像，灰度直方图反映了图像中不同灰度像素出现的统计情况，描述了像素灰度分布情况。数字图像的灰度直方图是一维离散函数，定义为


h(rk)=nk,k=0,1,…,L-1(34)



式中，rk表示第k个灰度级，nk表示灰度值为rk的像素在图像中出现的频数，L表示灰度级数。
在一幅图像中，灰度值为rk的像素出现的频数与像素总数的比值，称为概率直方图，可表示为


p(rk)=nkn,k=0,1,…,L-1(35)



式中，n为图像中的像素总数。p(rk)实际上表示灰度级rk,k=0,1,…,L-1的分布律。在概率直方图中，所有灰度级的概率之和等于1。
目前大多数数字照相机都有显示所拍摄照片直方图的功能，直方图可以显示出整张照片的灰度分布情况，根据直方图所示的灰度分布可以判断曝光是否恰当，有助于拍摄前的各种参数设置，如感光度、光圈、快门速度和曝光时间等。图32给出了数字单反照相机液晶显示器（liquid crystal display,LCD）上的直方图示例，其中横轴表示灰度级（对于彩色图像，实际上是亮度级），纵轴表示像素数。



图32数字单反照相机LCD屏上的直方图



图像的直方图具有如下三个主要性质。
(1) 直方图是总体灰度的概念，反映了整幅图像的灰度分布情况。亮图像的直方图倾向于灰度级高的一侧，暗图像的直方图倾向于灰度级低的一侧； 低对比度图像的直方图集中于灰度级中很窄的范围，高对比度图像的直方图覆盖的灰度级很宽而且像素的分布比较均匀。
(2) 直方图表示一幅数字图像中不同灰度像素出现的统计信息，它只能反映该图像中不同灰度值出现的频数（概率），不能表示出像素的位置等其他信息。图像与直方图之间是多对一的映射关系，具体地说，任意一幅图像，都能唯一地确定与之对应的直方图，但是，不同的图像可能具有相同的直方图。如图33所示，对图33(a)作行列置换，随机交换图像中的两行和两列不会改变图像的直方图，图33(b)和图33(c)分别为50次和200次行列置换的结果，三幅图像的直方图如图33(d)所示。



图33行列置换的图像具有相同的直方图



(3) 直方图具有可叠加性，将图像划分为区域，各个区域可分别统计直方图，整幅图像的总直方图为各个区域直方图之和。如图34所示，将图33(a)所示图像划分为四个区域，四个区域的直方图相加等于整幅图像的直方图。



图34直方图的可叠加性



在图像处理中图像的直方图主要有如下四方面的应用。
(1) 直方图可以用于判断一幅图像是否合理地利用了全部可能的灰度级。直观上来说，若一幅图像的像素值充分占有整个灰度级范围并且分布比较均匀，则这样的图像具有高对比度和多灰阶。通过查看直方图，可以确定设备参数调整方向，或者灰度级变换规则，如伽马校正。
(2) 图像的视觉效果与其直方图之间存在对应关系，改变直方图的形状对图像产生对应的影响。因此，处理直方图可以起到图像增强的作用，例如，直方图均衡化、直方图规定化、直方图拉伸处理等。
(3) 通过观察直方图，可在图像分割中确定合适的阈值，尤其适用于双峰模式直方图的全局阈值分割问题，并能够根据直方图对区域进行像素数统计。
(4) 将直方图作为特征，通过度量直方图之间的距离，可用于模板匹配、目标识别等任务。
例31不同灰度范围图像的直方图
一幅图像的直方图可以提供图像的如下信息： ①各个灰度级像素出现的概率； ②图像灰度级的动态范围； ③整幅图像的大致平均亮度； ④图像的整体对比度情况。图35给出了直方图的四种基本类型。从这四幅图中可以直观地看出不同直方图类型对应的图像特征。图35(a)所示的图像整体亮度偏暗，在对应的直方图中像素值集中分布在灰度级的暗端； 相反，图35(b)所示的图像整体亮度偏亮，在对应的直方图中像素值集中分布在灰度级的亮端。图35(c)所示的图像灰度层次少，对比度低，看起来比较暗淡，在对应的直方图中，像素值集中分布在灰度级中部较窄的动态范围。图35(d)所示的图像看起来对比度高，且有层次感、细节丰富，在对应的直方图中像素值覆盖了较宽的灰度动态范围。从直方图也可以看到，该图像的像素值主要集中分布在低灰度级和高灰度级，中间灰度级的像素较少。



图35四种基本类型的图像直方图



3.2.2累积直方图
累积直方图（cumulative histogram）实际上是概率直方图p(rk)关于灰度级rk,k=0,1,…,L-1的累积概率，定义为


c(rk)=∑kj=0p(rj)=∑kj=0njn,0≤rj≤1;k=0,1,…,L-1(36)



式中，rk表示第k个灰度级，nk表示灰度值为rk的像素在图像中出现的频数，n为图像中的像素总数，c(rk)表示灰度值落在区间［0,rk］内的像素在图像中出现的总概率，L表示灰度级数。累积直方图一定是单调递增的（不一定严格单调递增），且第L个灰度级的累积概率值c(rk=L-1)=1。
例32不同灰度范围图像的累积直方图
图36为图35所示的四种直方图类型图像的累积直方图。由于图35(a)~(c)的直方图集中在狭窄的灰度级范围内，因此累积直方图的递增曲线呈现陡坡的上升趋势，如图36(a)~(c)所示。而图35(d)的直方图占用整个灰度级范围，其累积直方图的递增曲线平缓地上升，特别是在像素数少的灰度级区间，如图36(d)所示。



图36图35所示四种直方图类型图像的累积直方图



3.3灰度级变换


图37图像增强中的基本灰度级变换曲线

灰度级变换是一类最简单的图像增强技术，设r和s分别表示输入图像和输出图像的灰度级，如前所述，s=T(r)表示将灰度级r映射为灰度级s的变换。线性函数/反转函数、对数函数/指数函数和幂次函数是图像增强中三类常用的基本灰度级变换曲线。图37显示了这三类函数的典型例子，①和②分别为线性函数和反转函数，③和④是互为反函数的对数函数和指数函数，⑤和⑥是互为反函数的幂次函数。其中，输入灰度级r的范围归一化到区间［0,1］内。

3.3.1对数变换
对数变换的一般表达式为


图38不同底数的对数变换曲线




s=cloga(1+r)(37)


式中，c为常数，a为对数的底数。对数变换的作用是压缩图像的动态范围。图38显示了不同底数的对数变换曲线，图中c=1，底数越大，压缩程度越大。



对数变换的一个典型应用是傅里叶谱第4章将详细讨论傅里叶变换及其傅里叶谱。的显示。傅里叶谱的值为0~106数量级甚至更高，具有很大的动态范围，一般的成像和显示系统都采用均匀量化，因此，在
8位灰度级系统中无法显示傅里叶谱的细节。对数变换可以起到压缩图像动态范围的作用，减小了最大值与最小值之间的反差值。
例33对数变换在傅里叶谱显示中的应用


图39对数变换应用于傅里叶谱的显示





以图像方式显示傅里叶谱是对傅里叶变换可视化的一种手段，在显示之前通过计算傅里叶谱的对数变换来压缩傅里叶谱图像显示的动态范围，增强灰度细节。为了直观地说明对数变换在傅里叶谱图像显示中的作用，对图39(a)所示的灰度图像显示其傅里叶谱以及对数变换的傅里叶谱，分别如图39(b)和图39(c)所示，其中左图为图像显示，右图为对角方向的径向剖面图。傅里叶谱的值域范围为［9.144,7.227×106］，当这样的值在8位灰度级系统中均匀量化而显示时，如图39(b)所示，傅里叶谱中的最大幅度淹没了很大范围内相对低的幅度，对应的高频区域的细节信息几乎完全丢失了，在图像中显示为黑色，只能观察到相当高的幅度表现出的亮点。根据式(37)的对数变换对傅里叶谱进行动态范围压缩，则值域范围缩至［2.317,15.793］，再用8位灰度级来显示对数傅里叶谱，如图39(c)所示，也就是将对数傅里叶谱的范围线性拉伸到8位灰度级的量化范围［0,255］。与图39(b)所示的直接显示相比，傅里叶谱的细节可见度是显而易见的。在本书中显示的傅里叶谱都使用对数变换和重标度进行了处理。


3.3.2指数变换


图310不同底数的指数变换曲线

指数变换为对数变换的反函数，一般表达式为


s=c(ar-1)(38)


式中，c为常数，a为指数的底数。指数变换的作用是拉伸图像的动态范围。图310显示了不同底数的指数变换曲线，图中c=1，底数越大，拉伸程度越大。



指数变换的一个典型应用是对数变换的对消，在利用照度反射模型Land提出的Retinex算法和5.7节中的同态滤波都是以照度反射模型为基础的典型方法，这类方法的目的是从图像中去除场景照度分量的影响，估计实际反射分量。的图像增强方法中，通常首先对图像进行对数变换，在对数域中将乘法运算转换为加法运算，这样能够利用线性滤波进行图像处理，最后再使用指数变换对图像进行反变换。
3.3.3幂次变换
幂次变换的一般表达式为


图311不同α值的幂函数曲线，常数c=1


s=crα(39)


式中，c为常数，α为指数的幂。当常数c=1时，α和1/α的幂函数互为反函数。图311显示了不同α值的幂函数曲线，其中，输入灰度级r的范围归一化到区间［0,1］。α＞1和α＜1的幂函数具有相反的作用，当α＜1时，拉伸直方图灰度级暗端的动态范围，而压缩灰度级亮端的动态范围； 当α＞1时，拉伸直方图灰度级亮端的动态范围，而压缩灰度级暗端的动态范围。
在图像增强中，α＜1的幂函数的作用是提高图像中暗区域的对比度，而降低亮区域的对比度，因此，对于灰度级整体偏暗的图像，使用α＜1的幂函数可以增大灰度动态范围； α＞1的幂函数的作用是提高图像中亮区域的对比度，而降低暗区域的对比度，对于灰度级整体偏亮的图像，使用α＞1的幂函数可以增大灰度动态范围。

例34幂次变换在图像增强中的应用
通过图312和图313来观察α＜1和α＞1幂次变换在图像增强中的作用。图312(a)为一幅曝光不足的图像，右图为对应的直方图，直方图偏向灰度级的暗端。对于这种整体偏暗的情况，α＜1的幂函数能够拉伸灰度直方图的占用范围。图312(b)为α=0.5的幂次变换的结果及其直方图，由图中可以看到，这种灰度级变换增强了图像的对比度，并提高了脸部等暗区域的亮度。图313(a)为一幅曝光过度的图像，右图为对应的直方图，直方图偏向灰度级的亮端。对于这种整体偏亮的情况，α＞1的幂函数能够拉伸灰度直方图的占用范围。图313(b)为α=2的幂次变换的结果及其直方图，从直方图来看，这种灰度级变换使像素值整体移向直方图灰度级的暗端，并拉伸了图像的对比度。


图312α=0.5的幂次变换图像增强示例






图313α=2的幂次变换图像增强示例


3.3.4灰度反转
灰度反转是对图像求反，灰度反转变换的一般表达式为


s=1-r(310)


式中，输入图像和输出图像的灰度级范围为［0,1］。这样的明暗反转的作用是突出在大片黑色区域中的白色或灰色细节。
例35灰度反转在图像增强中的应用
图314给出了一个灰度反转的图例。图314(a)为一幅极短快门时间拍摄的一滴牛奶溅起水花的高速图像，图314(b)为图314(a)的灰度反转图像，在白色或明亮背景中能够更明显地突出黑色或灰暗的水花。普通黑白胶卷照片的负片与景物的明暗就是灰度反转关系。


图314灰度反转图像增强示例



3.3.5分段线性变换
分段线性变换的一般表达式为


s=s1r1r,0≤r＜r1

s2-s1r2-r1(r-r1)+s1,r1≤r≤r2

s2-1r2-1(r-1)+1,r2＜r≤1(311)


图315(a)显示了分段线性函数，点(r1,s1)和(r2,s2)的位置控制了变换函数的形状，其中r1≤r2，s1≤s2。曝光不足、曝光过度和传感器动态范围过窄等都会造成图像表现出低对比度。分段线性变换的作用是提高图像灰度级的动态范围，是一种重要的灰度级变换。
当s1=0且s2=1时，式(311)的分段线性变换可简化为


s=r-r1r2-r1,r1≤r≤r2(312)


如图315(b)所示，将灰度级范围从［r1,r2］线性拉伸到［0,1］，其中，r1和r2分别为线性拉伸的下限和上限。通过截断一定比例的最亮像素和最暗像素，并使中间亮度像素占有整个灰度级，因而能够提高图像的全局对比度，如图315(c)所示，图中Plow和Phigh分别表示截断的最暗和最亮像素的比例，rlow和rhigh分别表示图像所占灰度级范围的最小值和最大值。在这种情况下，通常称为对比度拉伸或直方图剪裁，广泛应用于图像预处理中。


图315分段线性函数及其一种特殊情况



当r1=s1且r2=s2时，式(311)的分段线性函数退化为线性函数。当r1=r2，s1=0且s2=1时，式(311)的分段线性函数退化为阈值函数。


例36对比度拉伸在图像增强中的应用
图316给出了一个对比度拉伸应用于图像对比度增强的图例。如图316(a)所示，水下图像对比度低，由直方图可见，图像的像素值集中分布在较窄的灰度级范围内。MATLAB图像处理工具箱中的stretchlim函数根据截断像素的比例计算对比度拉伸的上下限，imadjust函数根据上下限进行对比度拉伸处理。截断1%的最暗像素和最亮像素，确定对比度拉伸的下限r1=0.4392和上限r2=0.7412，在直方图中用符号“×”标记。将输入图像的灰度级范围从［r1,r2］线性拉伸到［0,1］。图316(b)为对比度拉伸的结果，由图中可见，图像的对比度明显增强。注意，对比度增强包括全局对比度增强和局部对比度增强。对比度拉伸是一种全局对比度增强方法，这类方法能够增强图像中某些局部区域的对比度，同时也会降低其他局部区域的对比度； 而局部对比度增强利用图像的局部特性，同时增强图像的全局和局部对比度。







图316对比度拉伸图像增强示例



3.3.6灰度切片
灰度切片的作用是在整个灰度级范围内将灰度窗口内的灰度值与其他灰度值分开，突出图像中灰度窗口的区域，灰度窗口是指给定的灰度级范围。灰度切片有两种类型： 清除背景和保持背景。前者将灰度窗口内的像素保持不变［图317(a)左图］或赋值为较亮的值［图317(a)右图］，而其他部分赋值为相同的较暗值。显然，图317(a)右图所示的变换函数处理后产生的是一幅二值图像，细节完全丢失。后者将灰度窗口内赋值为相同的较亮值，而其他部分的灰度值保持不变，如图317(b)所示。



图317灰度切片的灰度级变换函数




图318两种灰度切片图像增强示例






例37灰度切片在图像增强中的应用
图318给出了两种类型灰度切片的图例。对于图318(a)所示的图像，从图318(b)所示的直方图可见，该图像灰度级大致可分为三段： 暗部、中部和亮部。设定灰度窗口为中部灰度区间［0.3137,0.6667］，图中符号“×”标识灰度窗口的上限和下限。图318(c)为清除背景型两种方式的灰度切片结果，左图中灰度窗口内的像素值保持不变，而灰度窗口之外暗部和亮部的全部像素置为最小值0； 右图中将灰度窗口内的像素赋值为最大值1，灰度窗口之外的像素赋值为最小值0，显然，形成一幅二值图像，图中只有0和1两种灰度值。图318(d)为保持背景型的灰度切片结果，将灰度窗口内的全部像素赋值为1，而暗部和亮部的灰度区间维持原状。


3.3.7阈值增强
阈值增强的作用是生成一幅高对比度的图像，阈值增强变换的一般表达式为


s=11+(r0/r)γ(313)



式中，γ控制函数的斜率，r0表示明暗跃变点。当r=r0时，s=1/2。图319显示了阈值增强函数曲线（折线），其中符号“×”标识明暗跃变点r0。如图319(a)所示，软阈值增强函数将输入灰度级小于r0的输出值压缩到狭窄的暗灰度级范围内，将输入灰度级大于r0的输出值压缩到狭窄的亮灰度级范围内。当γ→∞时，图319(a)所示的软阈值增强函数将逼近图319(b)所示的硬阈值增强函数。若输入灰度级小于r0，则输出值为0，否则输出值为1，硬阈值增强函数的输出是一幅二值图像。



图319阈值增强函数




图320阈值增强图像增强示例



图320（续）







例38阈值增强在图像增强中的应用
图320(a)为日全食钻石环的拍摄图像，利用式(313)中阈值增强函数来增强钻石环与背景的对比度，这里的r0取值为0.65。图320(b)为该图像的直方图，图中符号“×”标识明暗跃变点r0。图320(c)使用的是图319(a)所示的软阈值增强函数，γ取值为10，明暗过渡区域的像素为［0,1］之间的值，使过渡区域边缘柔化； 而图320(d)使用的是图319(b)所示的硬阈值增强函数，前景区域赋值为1，而背景区域赋值为0，没有过渡区域，产生一幅阈值化的二值图像。


3.4直方图处理
直方图处理也是一种点处理方法，它通过对灰度直方图进行变换，从而扩展灰度级的动态范围、提高图像的对比度、增强局部细节等。直方图处理使直方图变换为均匀分布函数或者其他特定的函数形式，主要有直方图均衡化和直方图规定化两种方法。直方图均衡化是通过灰度级变换将输入图像转换为具有均匀分布直方图的图像。更一般的情况，直方图规定化是指通过灰度级变换将输入图像转换为具有特定分布直方图的图像。
设r表示输入图像的灰度级，s表示增强图像的灰度级，其中，0≤r,s≤1，直方图处理本质上是选取灰度级变换函数s=T(r)，T(r)的选取应满足如下两个条件： 
(1) T(r)在区间0≤r≤1内为严格单调递增函数； 
(2) 对于0≤r≤1，有0≤T(r)≤1。

第一个条件保证了灰度级从暗到亮的次序不变，并且保证了反函数的存在； 第二个条件保证了输出灰度级与输入灰度级具有相同的灰度级范围。满足上述条件的反函数表示为


r=T-1(s)(314)



反函数同样满足上述两个条件。
3.4.1直方图均衡化
当直方图中的像素值集中在狭窄的灰度级范围内或分布极不均匀时，图像呈现较差的对比度，而当一幅图像的直方图具有较宽的灰度级范围，且分布较为平坦时，这样的图像具有高对比度。直方图均衡化的目的是将直方图的灰度级概率分布变换为均匀分布。为了便于描述，通过连续型随机变量的概率密度函数来说明直方图均衡化的机理。以随机变量R=R（r）表示输入图像的灰度级r，pr(r)表示R的概率密度函数，随机变量S=S(s)表示输出图像的灰度级s，ps(s)表示S的概率密度函数。图321直观地说明了直方图均衡化的目的。



图321直方图均衡化的目的



根据概率论的知识，已知S为均匀分布，且S是关于R的函数s=T（r），推导S与R的映射关系。当概率密度函数存在时，累积分布函数（也称为分布函数）是概率密度函数的积分，描述了连续型随机变量的概率分布。概率密度函数则是分布函数的导函数。设Fs(s)表示随机变量S的分布函数，Fr(r)表示随机变量R的分布函数，根据分布函数的定义，则有


Fs(s)①P{S＜s}②P{T(R)＜s}

③P{R＜T-1(s)}④P{R＜r}⑤Fr(r)(315)


式中，P(·)表示随机事件发生的概率； T(r)为严格单调递增函数，即r=T-1(s)存在。①、⑤应用了分布函数的定义； ②应用了随机变量S是R的函数s=T(r)； ③应用了T(r)为严格单调递增函数； ④应用了r=T-1(s)。式(315)建立了输出灰度级s与输入灰度级r之间的映射关系。
随机变量S为均匀分布，其概率密度函数为

ps(s)=1（316）

随机变量S的分布函数可写为概率密度函数ps(s)的变上限积分形式，且根据式(316)
可知 

Fs(s)=∫s0ps(x)dx=s (317)

结合式(315)
和式(317)
可得 

s=Fs(s)=Fr(r)=∫r0pr(x)dx,0≤r≤1(318)

式中, ps(s)为Fs(s)的导函数, pr(r)为Fr(r)的导函数。式(318)
右边表示连续型随机变量R
的分布函数Fr(r)。由上述推导过程可得出结论： 当输出图像的灰度级s是输入灰度级r的累积分布函数s=Fr(r)时，变换后的概率密度函数ps(s)为均`匀分布。


图322直方图均衡化的机理

图322
直观地描述了直方图均衡化的机理。连续型随机变量的取值落在某个区间之内的概率是由概率密度函数在这个区间上的积分给出。由于随机变量S是随机变量R的函数s=T(r), 随机变量R的取值落在区间［r,r+Δr］内的概率等同于随机变量S的取值落在区间［s,s+Δs］内的概率, 即ps(s)ds=pr(r)dr。当函数的定义域和值域都在实数域中时，函数导数的几何意义是曲线上的切线斜率。若当Δr→0时Δs与Δr之比的极限存在，则曲线s=Fr(r)在点P(r,s)处的切线斜率是函数s=Fr(r)在点r处的导数,即

pr(r)=F′r(r)=limΔr→0ΔsΔr

即ds=pr(r)dr。这从几何意义上说明, 直方图均衡化的变换函数s=Fr(r)为灰度级r的累积分布函数时，输出灰度级s服从均匀分布，即pss=1。
在灰度级为离散形式下，对于L个灰度级的数字图像，第k个灰度级为rk, 具有灰度值rk的像素数为nk, 像素总数为n, 设均匀直方图组(bin)数为M, 且L≥M,变换函数sl=Frrk是概率直方图prk的概率累加，k=0,1,…,L-1； l=0,1,…,M-1。直方图均衡化的具体可执行步骤如下： 
(1) 计算输入图像灰度级rk的分布律prrk, 统计各灰度级rk的像素出现的概率为

prrk=nkn,k=0,1,…,L-1(319)


(2) 计算输入图像灰度级rk的累积概率Frrk ： 

ck=Frrk=∑kj=0prrj=∑kj=0njn,k=0,1,…,L-1(320)


(3) 根据灰度映射规则ck→sl, 确定输入灰度级与输出灰度级之间的映射关系rk→sl,将输入图像中灰度级为rk的像素映射到输出图像中灰度级为sl的对应像素。
(4) 统计输出图像各灰度级sl的像素数，并计算输出图像灰度级sl的分布律pssl。
灰度映射规则ck→sl对于各个k寻找l满足： 

minlck-sl=ck-l+1M,k=0,1,…,L-1(321)

式中，sl=l+1M为均匀直方图的累积直方图， M为直方图的组数。由式(321)
可推导出

l=maxround Mck-1,0(322)

其中，round［·］表示四舍五入操作。显然，在映射过程中可能会造成将多个不同的rk值映射到相同的灰度级sl，不同输入灰度级的像素在灰度级映射后归并到同一输出灰度级sl的像素中。表31给出了一个直方图均衡化的具体计算过程，输入图像的灰度级范围为［0,L-1］，均匀直方图组数为M。图323比较了表31中具体示例在直方图均衡化处理前后的概率直方图，图323(a)为输入直方图，图323(b)为输出直方图。由图中可见，直方图中同一灰度级的像素不能拆分，不同灰度级的像素通常会合并。因此，所使用的灰度级数会减少，从而造成细节信息的丢失。对于离散变量而言，直方图均衡化无法实现完全均匀分布的直方图，输出直方图的概率近似为均匀分布。


表31直方图均衡化的计算步骤，像素总数n为4096，灰度级数L为8，均匀直方图组数M为8



步骤

序号

rk0
1
2
3
4
5
6
7
nk508
821
898
892
552
181
159
85



1
               pr(rk)=nk/n            0.124
0.2
0.219
0.218
0.135
0.0442
0.0388
0.0208
2a
               sk=Fr(rk)            0.124
0.324
0.543
0.761
0.896
0.940
0.979
1
2b
               L=max［round(Mck)-1,0］0
2
3
5
6
7
7
7

3
               rk→sl0→0
1→2
2→3
3→5
4→6
5,6,7→7
4
               ps(sl)            
0.124
0
0.2
0.219
0
0.218
0.135
0.104



图323直方图均衡化前后直方图比较



例39直方图均衡化在图像增强中的应用
图324(a)为一幅低对比度图像及其直方图。从直方图中可以看出，该图像所使用灰度级的动态范围较窄。使用MATLAB图像处理工具箱的直方图均衡化histeq函数进行对比度增强处理。图324(b)为直方图均衡化图像及其直方图，均匀直方图具有64个组，即M=64。当M远小于灰度图像中的灰度级数L时，输出图像的直方图更加平坦。从直方图中可以看出，图像的动态范围扩展到整个灰度级上，并且大致平坦，在图像中表现为对比度明显增强。图324(c)为输出图像与输入图像之间的灰度级变换函数，直方图均衡化的变换函数是输入图像的累积概率直方图，从图中可以看出，输入图像具有较窄的灰度级范围，而变换后图像的动态范围扩展到整个灰度级范围。
需要说明的是，直方图均衡化使一幅图像的像素值占有全部可能的灰度级且分布尽可能均匀，尽管能够从视觉效果上提高图像的对比度，但是，由于直方图中概率较小的灰度级合并为更少的几个或一个灰度级内，从而降低了图像的灰度级分辨率，且某些细节信息处于概率较小的灰度级中，这样的灰度级归并到其他灰度级中，从而造成图像细节信息的丢失。从图324(b)可以看出，由于前景与背景过渡区域的灰度值略不同于区域内部像素的灰度值，且这部分像素值的概率较小，因此，在直方图均衡化过程中，合并到概率较大的灰度级中，从而造成前景轮廓部分发生失真，且灰度层次减少。从图324(c)所示的灰度级变换函数也可以分析出这样的失真出现的原因，由于该变换函数并非严格单调递增，因此，它的反函数是奇异的，输入的多个灰度级会映射到输出更少的灰度级上。
已有的研究表明，一种好的图像增强（或复原）算法应使处理后图像的直方图与原图像的直方图在总体形状上保持一致。因此，直方图均衡化的主要问题是造成直方图的形状发生失真，从而使处理后的图像看起来缺乏真实性和自然性。



图324直方图均衡化图像增强示例



3.4.2直方图规定化
直方图均衡化尽管能够增强图像的对比度，然而它无法控制具体的增强效果。更一般的情况，要求图像具有特定形状的直方图，以便于有选择地对图像中某个特定的灰度级范围进行增强，或使其满足后续处理的特定需求。直方图规定化的目的是将具有已知直方图的图像变换为具有某种特定形状直方图或者用户交互指定形状直方图的图像，也称为直方图匹配。直方图均衡化是直方图规定化的特例，它将输入图像的直方图变换为均匀分布的直方图。

为了表述简便，仍从连续型随机变量的概率密度函数入手对直方图规定化进行讨论。设pr(r)表示输入图像灰度级r的概率密度函数，pv(v)表示规定直方图组v的概率密度函数，通过建立输入直方图pr(r)和规定直方图pv(v)之间的关系，求映射函数v=T(r)。
对于连续变量的输入直方图和规定直方图，利用累积直方图均衡化作为中间桥梁，累积直方图分别写为


s=Fr(r)=∫r0pr(x)dx(323)


ν=Fv(v)=∫v0pv(x)dx(324)



由于累积直方图的输出服从均匀分布（直方图均衡化的过程），因此，s和ν的概率密度函数关系可表示为


ps(s)=pν(ν)=1,0≤s,ν≤1(325)



从统计意义上讲，s和ν是相同的，这样，可推导出v关于r的映射函数v=T(r)为


v=F-1v(ν)=F-1v(s)=F-1v(Fr(r))(326)



从而实现将已知灰度图像转换成规定直方图的图像。图325直观地说明了直方图规定化的实现过程，图中，方框中间的符号“=”表示依概率相等。



图325直方图规定化的过程



设数字图像具有L个灰度级，输入图像的灰度级为rk，具有灰度值rk的像素数为nk，规定直方图组数为M，且L≥M，直方图规定化的具体可执行步骤如下： 
(1) 对输入图像的概率直方图pr(rk)计算其累积直方图： 


sk=Fr(rk)=∑ki=0pr(ri),k=0,1,…,L-1(327)



(2) 对规定直方图pv(vj)计算其累积直方图： 


νl=Fv(vl)=∑lj=0pv(vj),l=0,1,…,M-1(328)



(3) 根据灰度映射规则，将输入灰度级映射到输出灰度级； 
(4) 确定输入灰度级与输出灰度级之间的映射关系rk→vl，将输入图像中灰度级为rk的像素映射到输出图像中灰度级为vl的对应像素； 
(5) 统计输出图像各灰度级vl的像素数，并计算输出图像的概率直方图p^v(vl)。
对于离散的数字图像，存在取整误差的影响，灰度映射规则是直方图规定化的关键。不同的映射关系，会产生不同的规定化结果。常用的灰度映射规则包括单映射规则（single mapping law,SML）和组映射规则（group mapping law,GML）。
单映射规则sk→νl对于各个k寻找l满足： 


minl|sk-νl|=∑ki=0pr(ri)-∑lj=0pv(vj),k=0,1,…,L-1(329)



式中，L为输入图像的灰度级数，M为规定直方图的组数，且L≥M。
组映射规则νl→sk确定满足下式的最小整数函数k(l)： 


mink(l)|sk(l)-νl|=∑k(l)i=0pr(ri)-∑lj=0pv(vj),l=0,1,…,M-1(330)



其中，k(l),l=0,1,…,M-1为整数函数，满足


0≤k(0)≤…≤k(l)≤…≤k(M-1)≤L-1



若l=0，则将i=0到i=k(0)的pr(ri)映射到pv(v0)； 若l≥1，则将i=k(l-1)+1到i=k(l)的pr(ri)映射到pv(vl)。
设输入图像的灰度级范围为［0,L-1］，表32给出了一个直方图规定化的具体计算步骤。直方图规定化的列表计算方式不直观，一种简单且直观的方式是绘图计算方式。绘图计算方式用灰度条表示累积直方图，每一段对应直方图中的一个灰度级。图326直观地说明了直方图规定化绘图计算方式的过程，单映射规则对输入累积直方图中各个灰度级的概率累积值依次向规定累积直方图映射［图326(a)］，组映射规则对规定化累积直方图中各个灰度级的概率累积值依次向输入累积直方图映射［图326(b)］，选择最接近的数值，即最短的连线。图中实线和虚线描述了sk和νl之间可能的近邻关系，实线比虚线更短，数值更接近，因此实线确定了输入累积直方图到规定化累积直方图的映射。


表32直方图规定化的计算步骤，像素总数n为4096，灰度级数L为8，规定直方图组数M为8



步骤

序号

               rk            0
1
2
3
4
5
6
7
               nk            508
898
892
821
552
181
159
85

1a
               pr(rk)=nk/n            0.124
0.219
0.218
0.2
0.135
0.0442
0.0388
0.0208
1b
               sk=Fr(rk)            0.124
0.343
0.561
0.761
0.896
0.94
0.979
1
2a
               pv(vl)            0
0
0.1
0.15
0.2
0.25
0.2
0.1
2b
               νl=Fv(vl)            0
0
0.1
0.25
0.45
0.7
0.9
1
单映射

sk→νl            
3
               vl            2
3
4
5
6
6
7
7

4
               rk→vl            0→2
1→3
2→4
3→5
4,5→6
6,7→7
5
               p^v(vl)            0
0
0.124
0.219
0.218
0.2
0.179
0.06

组映射
νl→sk            
3
               vl            2
3
4
5
6
7
7
7
4
               rk→vl            0→2
1→3
2→4
3→5
4→6
5,6,7→7
5
               p^v(vl)            0
0
0.124
0.219
0.218
0.2
0.135
0.104





图326直方图规定化的绘图计算过程



图327比较了表32中具体示例在直方图规定化处理前后的概率直方图，图327(a)和图327(b)分别为输入直方图和规定直方图，图327(c)和图327(d)分别为单映射和组映射计算的输出直方图。由于直方图是离散化的，不同灰度级的像素可以合并，同一灰度级的像素不能拆分。对于离散变量而言，无法做到与规定直方图完全相同，输出直方图与规定直方图近似保持一致。



图327表32直方图规定化前后直方图比较



例310直方图规定化在图像增强中的应用
从图324(b)所示的直方图均衡化结果可以看到，直方图均衡化引起明显的图像失真。对于图324(a)所示的低对比度图像，图328给出了直方图规定化的结果。如前所述，图像增强应保持原图像直方图的大致形状。为了保持了原直方图的大体形状，选择暗部主体的瑞利分布函数构造规定直方图，如图328(a)所示，规定的直方图利用了整个灰度级的动态范围，其组数为64。直方图规定化也使用MATLAB图像处理工具箱的histeq函数，需要指定规定直方图。图328(b)为直方图规定化图像及其直方图，当规定直方图的组数远小于灰度图像中的灰度级数时，输出图像的直方图更好地匹配规定直方图。相比图324(b)所示的直方图均衡化图像，直方图规定化的图像看起来真实、自然。



图328直方图规定化示例



3.5算术运算
图像中的算术运算是在逐像素基础上对两幅或多幅相同尺寸的输入图像之间进行加、减、乘、除运算，从而达到某种增强的目的。在上述四种算术运算中，减法运算和加法运算在图像增强处理中最为有用。本节将讨论减法运算和加法运算在图像增强中的应用。基于像素的算术运算每次处理一个像素，并与其他像素的处理无关，因而适合于并行操作且容易在硬件上实现。
3.5.1图像相减
设f(x,y)和g(x,y)表示两幅相同尺寸的数字图像，两幅图像的减法运算是计算这两幅图像对应像素的灰度差值d(x,y)，可表示为


d(x,y)=f(x,y)-g(x,y)(331)



在图像处理中，图像相减主要有三方面的应用： ①显示两幅图像的差异； ②检测同一场景两帧图像之间的变化，例如，视频中镜头边界的检测； ③消除图像中不需要的加性分量，例如，缓慢变化的背景、阴影、周期性噪声等。
在运动目标检测与跟踪的应用中，背景减除法和帧间差分法是两种常用的运动目标检测方法，可用于视频流中的运动目标检测。背景减除法是比较视频流中当前帧图像与事先存储的或实时获取的背景图像中对应像素的差异，分割出前景运动目标。帧间差分法是通过视频帧中的同一场景相邻帧之间的变化检测出前景运动目标。
减法运算可能产生负值，而显示系统要求像素值为无符号整型，因此，在显示之前需要将差值进行重标度（尺度化）。一种最常用的重标度方法是将差值线性映射到显示允许的灰度级范围内，如图329所示。重标度的差值图像dscale(x,y)可表示为


dscale(x,y)=roundd(x,y)-dmindmax-dmin·fmax(332)



式中，round［·］表示四舍五入符号，dmin和dmax为差值d(x,y)的最小值和最大值，fmax为显示允许的最大灰度值。通用的显示设备使用8位无符号整型256个灰度级显示图像，在这种情况下，fmax为255。式(332)的数学表达式可分解为如下具体的步骤： ①查找差值中的最小值dmin和最大值dmax，对每一个像素值减去dmin，将差值的范围从［dmin,dmax］线性平移到范围［0,dmax-dmin］； ②乘以尺度因子1/(dmax-dmin)进行比例缩放到归一化范围［0,1］； ③乘以显示允许的最大灰度值fmax进行比例缩放到显示范围［0,fmax］； ④四舍五入为正整数。



图329重标度示意图



例311减法运算在运动目标检测中的应用
结合图330中背景建模和图331中背景减除的示例说明减法运算在固定摄像头运动目标检测中的应用。如图330所示，对连续n帧图像中每一个像素的R、G、B颜色分量排序，求取中间值作为背景的像素值，生成背景图像，这是一种简单且常用的固定摄像头背景建模方法。然后，通过计算当前帧图像与背景图像的差值，来检测前景运动目标。如图331所示，图331(a)和图331(b)分别为当前帧图像和背景图像，对它们的差值图像进行阈值化，从而将前景运动目标从背景中分割出来。在图331(c)所示的二值图像中，白色区域为检测出的运动目标。噪声的影响会使一些背景区域被误检测成前景运动区域，也可能使运动目标内的部分区域漏检。另外，由于背景的扰动（如树枝的轻微摇动），这部分背景区域也会误判为运动目标。在第8章中将使用数学形态学算法或连通分量分析进行图像后处理。


图330背景图像的生成示意图




图331背景减除法用于运动目标检测




例312减法运算在帧间变化检测中的应用
图332展示了一个帧差法用于变化检测的示例。图332(a)和图332(b)是QCIF格式CIF（common intermediate format）格式视频的分辨率为352×288像素，QCIF（quarter common intermediate format）格式视频的分辨率为176×144像素，D1格式视频的分辨率为704×576像素。视频中的连续两帧。两帧之间的残差如图332(c)所示，残差有正值也有负值。在帧差图像中，灰色表示差值为0，更亮和更暗的像素分别对应正值和负值。帧间的变化可能是目标运动引起的，包括刚性目标运动（如运动的轿车），非刚性目标运动（如运动的手臂）； 可能是照相机运动引起，如平移、倾斜、缩放、旋转； 可能由区域遮挡引起，如运动目标从部分场景中移开； 也可能由光照变化引起。除了区域遮挡和光照变化以外，其他的变化都对应帧间的像素运动。因此，估计视频中连续两帧之间像素的运动轨迹是可能的，产生的像素轨迹称为光流场。图333显示了连续两帧之间的光流场，完整的光流场中每一个像素对应一个运动向量，但是为了清楚地显示，将两帧图像划分为8×8的块，每一个块对应一个运动向量。从光流场中可以估计从帧1到帧2的目标运动方向和强度。



图332减法运算用于帧间变化检测





图333图332(a)与图332(b)所示连续两帧之间的光流场（Lucas Kanade算法）



例313数字减影血管造影
在医学图像处理中，数字减影技术实际上是减法运算。数字减影血管造影（digital subtraction angiography,DSA）是通过对注入造影剂之后与注入造影剂之前的血管造影影像进行减法操作，消除不需要的组织影像，仅保留血管影像。图334(a)和图334(b)分别为注入造影剂之前和注入造影剂之后的血管造影影像，图334(c)为图334(b)与图334(a)相减的结果，并将灰度级范围重标度以8位灰度级图像显示。由图中可见，数字减影的特点是目标清晰，对血管病变的观察、血管狭窄的定位测量、诊断及介入治疗提供了必备条件。



图334减法运算用于数字减影血管造影



3.5.2图像相加
设f(x,y)和g(x,y)表示两幅相同尺寸的数字图像，两幅图像的加法运算是计算这两幅图像对应像素的灰度值之和s(x,y)，可表示为


s(x,y)=f(x,y)+g(x,y)(333)



在图像处理中，图像相加主要有两个方面的应用： ①对同一场景的多幅序列图像求取平均值，从而降低加性随机噪声的影响； ②将一幅图像叠加到另一幅图像上，以达到二次曝光的效果。其中，第一个应用更为广泛，下面将具体解释多幅图像相加实现降噪的机理。
若成像系统受加性噪声干扰，则所成图像可用如下的降质模型表示为： 


g(x,y)=f(x,y)+η(x,y)(334)



式中，g(x,y)为有噪图像，f(x,y)为原图像，η(x,y)为加性噪声。假设每一个像素(x,y)处的噪声η(x,y)是独立同分布、且均值为零的随机变量。根据概率论的知识，将有噪图像g(x,y)看作随机变量，其期望和方差为


E［g(x,y)］=f(x,y)+E［η(x,y)］=f(x,y)(335)

Var［g(x,y)］=E{［g(x,y)-E(g(x,y))］2}=E{［g(x,y)-f(x,y)］2}

=E{［η(x,y)］2}=Var［η(x,y)］(336)



式中，E(·)表示期望，Var(·)表示方差。式(335)和式(336)中最后的等号成立均是由于每一个像素(x,y)处噪声η(x,y)的期望为零。
假设对同一场景f(x,y)连续K次成像，K幅不同的有噪图像{gi(x,y)}i=1,2,…,K的样本均值g-(x,y)为


g-(x,y)=1K∑Ki=1gi(x,y)=1K∑Ki=1［f(x,y)+ηi(x,y)］(337)



式中，gi(x,y)为第i次采集的有噪图像，ηi(x,y)为第i次成像中的加性噪声。根据期望和方差的性质，均值图像g-(x,y)的期望和方差分别为


E［g-(x,y)］=f(x,y)(338)

Var［g-(x,y)］=1KVar［η(x,y)］(339)



式(338)表明，f(x,y)是g-(x,y)的无偏估计； 式(339)表明，当帧数K增加时，在每一个像素(x,y)处g-(x,y)偏离f(x,y)的方差减小，噪声的包络收窄。由这两式联立可知，随着所采集图像数目的增加，多幅图像相加使g-(x,y)更好地逼近f(x,y)，从而达到降低随机噪声的目的。在实际应用中，由于无法保证成像场景完全相同，因此需要事先对多幅图像gi(x,y)进行图像配准。
例314视频中字幕文本的分割
文本分割是将文本从具有复杂背景的图像或视频中分割出来，转换为二值图像，提供给光学字符识别引擎进行识别。复杂背景下的文本分割问题是视频中文本识别的一个关键环节。文本分割包括场景文本分割和字幕文本分割。对于场景文本，通常为了可读性，将场景文本叠加在纯色的背景上，且非常醒目。因此，场景文本分割较为容易。字幕文本通常直接叠加在视频上，而视频的复杂性增加了字幕文本分割的难度。由于视频中连续多帧中的字幕文本是静止的，而背景是显著变化的，因此，一种有效的方法是通过序列帧相加来削弱背景的影响。本例从某一视频中截取同一字幕的326帧文本图像，图335(a)显示了其中部分帧的文本图像，对这些序列帧逐帧累加求和，然后计算平均值，如图335(b)所示。对图335(b)进行简单的阈值化操作，阈值取值为0.56，图335(c)为最终分割的二值图像。由图中可见，通过对多个不同背景的同一字幕视频帧进行相加，可以显著削弱背景的影响。



图335序列帧文本图像相加以及阈值分割示例


3.6空域滤波基础
空域滤波是一种邻域处理方法，直接在图像空间中对邻域内像素进行处理，达到图像平滑或锐化的作用。空域滤波是图像处理领域中广泛使用的主要工具。“滤波”这一术语来自频域，频域滤波是指允许或者限制一定的频率成分通过，而空域滤波直接在图像空间中增强图像的某些特征或者削弱图像的某些特征。空域滤波的作用域是像素及其邻域，通常使用空域模板对邻域内像素进行处理而产生该像素的输出值。按照数学形态分类，空域滤波可分为线性滤波和非线性滤波。根据信号处理理论可知，线性滤波可写成卷积运算的形式，卷积的概念是线性滤波的基础，第4章介绍的卷积定理是线性系统分析中最有力的工具之一。按照处理效果分类，空域滤波可分为平滑滤波和锐化滤波。3.7节和3.8节分别从平滑滤波和锐化滤波的角度讨论空域滤波方法。
3.6.1卷积与相关
卷积，也称为线性卷积，是对两个函数f(t)和h(t)生成第三个函数的数学运算，其表征f(t)和经过反转、移位的h(t)的乘积函数所围成的曲边梯形的面积。连续形式的卷积定义为


f(t)*h(t)=(f*h)(t)=∫+∞-∞f(τ)h(t-τ)dτ(340)



式中，积分变量τ是一个虚变量，只对计算积分有用，没有具体含义。h(t-τ)表示在τ的坐标系中h(τ)需要进行反转和移位，然后将f(τ)与h(t-τ)的重叠部分相乘作积分。在离散时间系统中，卷积用离散形式表示为


f(n)*h(n)=(f*h)(n)=∑+∞k=-∞f(k)h(n-k)(341)



通常f(n)表示离散信号，h(n)表示系统的单位冲激响应，由于单位冲激响应的长度远小于信号的长度，一般对h(n)进行反转和移位操作。式(341)中求和变量为k，激励信号f(k)和系统的单位冲激响应h(n-k)都是k的函数。本节从两种不同的角度解释卷积运算的公式，分别描述对应的两种卷积运算的过程。
第一种方法是固定时间n，变换自变量，将n改写为k，对h(k)先反转再移位，将信号f(k)看成整体，一次计算某一时间n的响应值，卷积在某一时间n的响应值的计算过程包括如下四个步骤： 
(1) 反转将h(k)关于k=0(时间原点）反转产生h(-k)。
(2) 移位在时间上对反转序列进行移位，如果n是正数（负数），那么将h(-k)右（左）移n个单位产生h(n-k)。 
(3) 相乘将f(k)和h(n-k)相乘构成乘积序列vn(k)≡f(k)h(n-k)。
(4) 求和将乘积序列vn(k)的所有值相加构成在时间n的输出值。
上述步骤只是计算出在某一时间n的系统响应。为了获取系统响应在所有时间的值，-∞＜n＜+∞，在求和时，对所有可能的时间移位-∞＜n＜+∞，都需要重复步骤(2)~(4)。
例315离散序列的卷积(1)
为了更好地理解卷积运算的步骤，用图形方法描述计算过程。图形有助于解释计算卷积的四个步骤。图336(a)和图336(b)分别为9点的输入信号序列f(k)和7点的冲激响应序列h(k)，k为时间变量。计算卷积的第一步是反转h(k)，反转序列h(-k)如图336(c)所示，反转操作仅执行一次。根据式(341)，当位移n=0时，直接利用h(-k)而无须移位。将图336(a)和图336(c)的序列中对应位置的元素相乘产生乘积序列，最后，将乘积序列的所有项相加产生n=0时的输出，在位移n=0时，f(k)和h(-k)仅在时间k=0处相交，很容易得出n=0时的输出值为1（见图336(r)中n=0时的输出值）。继续计算系统在n=1至n=14时的响应，序列h(n-k)是反转序列h(-k)在时间上右移n个单位，如图336(d)~(q)所示，图中红色标出与图336(a)中输入信号序列重合的部分。将图336(a)分别与这些移位的序列对应元素相乘，再将所有的乘积相加产生在时间n处的响应。对于n＜0或者n＞14，由于这两个序列没有重合部分，乘积序列全为零。图336(r)为这两个序列卷积的计算结果，线性卷积结果的序列长度为15。


图336第一种卷积计算过程的释义



第二种方法是固定移位k，一次计算某一移位k上对采样信号f(k)在所有时间的响应值，再将所有位置上的响应值对位相加。利用如下三个步骤可直接实现卷积的计算过程： 
(1) 移位在时间上对序列h(n)进行移位，如果k是正数（负数），那么将h(n)右（左）移k个单位产生h(n-k)。 





(2) 相乘将f(k)和h(n-k)相乘构成乘积序列v(n-k)≡f(k)h(n-k)。
(3) 求和将所有乘积序列v(n-k)对位相加一次计算出所有时间n的响应值。
乘积序列v(n-k)表示对h(n)延迟k的序列h(n-k)赋予对应位置处的信号采样值f(k)的权重。第二种计算卷积的方法是从卷积的物理意义出发。卷积为单位冲激函数诞生，用于描述系统对输入信号的响应。卷积的原理是将信号分解为冲激信号之和，借助系统的冲激响应，从而求解系统对任意激励信号的响应。因此，卷积公式可以理解为冲激响应函数在所有延迟k［即h(n-k)］的加权和，权系数是信号在时间k处的强度f(k)。卷积表征了系统响应与激励信号和单位冲激响应之间的关系线性移不变系统的输出可以表示为输入信号与系统冲激响应的卷积，这部分内容将在第6章讨论。。
例316离散序列的卷积（2)
根据卷积公式的第二种理解，利用对位相乘求和的方法可以更快地求出卷积结果。对于例315中的两个信号f(n)和h(n)，如图337(a)所示，这里将它们写成序列： 


f(n)={111111111}

h（n）={10.60650.36790.22310.13530.08210.0498}



f(n)是由单位冲激信号组成的序列，序列h(n)与单位冲激信号卷积，相当于在单位冲激信号的位置处复制该序列，如图337(b)所示，图中不同颜色的实心圆圈表示不同位置处序列h(n)的复制，空心圆圈为卷积结果。这可以写成如下对位排列的方式，将两个序列以各自n的原点对齐，逐个将f(n)的值与整个序列h(n)对应相乘，最后将同一列上的乘积按对位求和。最后一行是卷积的结果，下画线标出原点的位置。







图337第二种卷积计算过程的释义



卷积运算是可交换的，也就是说卷积结果与两个序列中任一个序列进行反转和移位是无关的。通过变量置换m=n-k，将求和变量从k变成m，则式(341)变为


(f*h)(n)=∑+∞m=-∞f(n-m)h(m)(342)



式(342)中系统的单位冲激响应不改变，对激励信号的序列进行反转和移位。
互相关，也可简称为相关，是两个信号的相似性度量，其定义与卷积相近，由下式给出： 


f(t)h(t)=(fh)(t)=∫+∞-∞f(τ+t)h(τ)dτ=∫+∞-∞f(τ)h(τ-t)dτ(343)



其离散形式可以写为


f(n)h(n)=(fg)(n)=∑+∞k=-∞f(k+n)g(k)=∑+∞k=-∞f(k)g(k-n)(344)


相关运算不满足交换律，即f(n)h(n)≠h(n)f(n)。借助变量置换容易得出卷积与相关的关系： 


f(n)h(n)=f(n)*h(-n)(345)


可见，将h(n)反转（变量取负号）后与f(n)卷积即得f(n)与h(n)的相关函数。卷积和相关运算都包含移位、相乘和求和三个步骤，其差别在于卷积运算开始时需要对h(n)进行反转，而相关运算不需要反转。
3.6.2线性滤波原理
在图像处理中，线性滤波处理利用空域滤波器与图像的空域卷积来实现，在线性滤波中通常使用空域模板来表示滤波器，也称为滤波模板、卷积模板或卷积核。将式(341)中的一维离散卷积推广到二维表示，二维离散卷积的定义为


g(x,y)=h(x,y)*f(x,y)=∑+∞m=-∞∑+∞n=-∞f(m,n)h(x-m,y-n)(346)


式中，g(x,y)为输出图像，f(x,y)为输入图像，h(x,y)为卷积模板。通常卷积模板的尺寸远小于图像的尺寸，因此对卷积模板进行反转和移位操作。卷积具有交换律，为了简化符号表示，将输入图像f(x,y)与卷积模板h(x,y)的卷积写为如下的形式： 


g(x,y)=∑+∞m=-∞∑+∞n=-∞h(m,n)f(x-m,y-n)(347)


通常要求卷积模板的尺寸为奇数，其中心为坐标原点，假设卷积模板的尺寸为(2k+1)×(2l+1)，根据卷积模板的尺寸，限定式(347)中的求和界限，则有


g(x,y)=∑km=-k∑ln=-lh(m,n)f(x-m,y-n)(348)


式中坐标m和n是整数。一般情况下，线性滤波选取尺寸为奇数的正方形模板。根据一维卷积的计算过程，实现式(348)的模板卷积的主要步骤包括如下四个步骤： ①将卷积模板反转，也就是将模板绕中心旋转180°； ②将模板在图像中遍历，将模板中心与各个像素位置重合； ③将模板的各个系数与模板对应像素值相乘； ④将所有乘积相加，并将求和结果赋值于模板中心对应的像素。
将反转的卷积模板定义为新的变量，即w(m,n)=h(-m,-n)，式(348)可改写为


g(x,y)=∑km=-k∑ln=-lw(-m,-n)f(x-m,y-n)

s=-m,t=-n∑ks=-k∑lt=-lw(s,t)f(x+s,y+t)(349)



式(349)表明，线性空域滤波本质上像素的输出值是计算该像素邻域内像素值的线性组合，将系数矩阵称为模板。当模板中心位于图像中像素坐标(x,y)处时，w(s,t)和f(x+s,y+t)分别为模板系数和模板对应图像中的像素值。根据相关的定义，式(349)实际上是计算图像f(x,y)和系数矩阵w(x,y)的相关f(x,y)w(x,y)。图338以尺寸为3×3的模板为例说明线性滤波的基本原理，线性滤波在图像中像素坐标(x,y)处的输出响应g(x,y)为


g(x,y)=w(-1,-1)f(x-1,y-1)+w(-1,0)f(x-1,y)+…+

w(0,0)f(x,y)+…+w(1,0)f(x+1,y)+w(1,1)f(x+1,y+1)(350)



这里s=-1,0,1； t=-1,0,1。当计算乘积与求和时，模板系数w(0,0)与像素值f(x,y)相对应。线性滤波的输出响应为模板系数与模板对应像素的乘积之和，最后将求和结果赋值于模板中心对应像素。在图像处理的出版物中，空域滤波通常未必使用真正的卷积，相关与卷积的区别仅在于模板的反转。根据实际应用的需求设计系数模板，因此直接使用相关运算来实现空域滤波并没有本质上的区别。



图338线性滤波的基本原理



为了便于描述，线性滤波在图像中像素坐标(x,y)处的输出响应R的另一种简化表达式为


R=∑ni=1wizi(351)



式中，wi表示模板系数，zi表示模板系数对应的像素值，i=1,2,…,n，n=(2k+1)×(2l+1)为模板对应的像素总数。以尺寸为3×3的模板为例，图339(a)为3×3模板的系数表示，图339(b)为模板在图像中对应的3×3邻域的像素值。式(351)表明，线性滤波的输出像素值为模板对应邻域内像素值的加权和，权值为模板系数。



图3393×3模板和3×3邻域的一种表示方式



例317卷积神经网络
卷积神经网络（convolutional neural network,CNN）中的卷积层由一组卷积模板（也称为卷积核、滤波器）组成，是图像特征提取的关键。卷积神经网络与本节中图像增强中所用的空域卷积本质上是相同的，其不同之处在于，在图像增强中，根据应用需求设计模板系数，模板系数是给定的或预设的； 而在卷积神经网络中滤波器的系数是通过标注样本对网络进行训练产生的，在训练过程中使这些滤波器组对特定的模式有较大的输出响应，以达到CNN的分类/检测等目的。若识别图像中某种特定模式，则滤波器组对该模式有较大的输出响应，称为神经元的激活。
卷积具有局部性，即它只关注局部特征，局部作用域取决于滤波器的尺寸。例如，Sobel 算子边缘检测本质是检测图像3×3邻域像素中的灰度变化或灰度突变。与全连接神经网络相比，卷积神经网络具有局部连接和权值共享的特殊结构，降低了网络的参数量，加快了学习速率，同时也在一定程度上减少了过拟合的可能。局部连接是指卷积层的节点仅与其前一层的部分节点相连接，仅能够学习局部特征。权值共享是指由于卷积的空间移不变性，滤波器和图像的卷积与位置无关，不会随着位置的改变而发生变化，这里权值是指滤波器系数。如图340所示，卷积操作就是滤波器系数与其对应像素的乘积之和，对于图中的情况，若仅考虑局部连接，共需要9×16=144个权值参数，而加上权值共享后，仅需9个权值参数，进一步减少了参数数量。图中，相同的滤波器系数用相同的颜色标识。



图340局部连接和权值共享



卷积神经网络包含多个卷积层，通过将低阶特征进行线性组合（激活函数引入非线性变换）形成高阶特征。如图341所示，卷积神经网络的第一个卷积层的滤波器组检测低阶特征，如边缘、角点、曲线等，第一层的输出是第二个卷积层的输入，这一层的滤波器组检测低价特征的组合等情况（半圆、四边形等）。随着卷积层的增加，滤波器组检测越来越复杂的特征。构建卷积神经网络的任务在于构建这些滤波器，通过标注样本训练使滤波器系数能够识别特定的模式。




图341卷积神经网络特征提取示意图（图片源自Ian Goodfellow等《深度学习》）



3.6.3可分离卷积
可分离性是多维卷积以交换律为基础的一个数学性质。利用式(347)中二维卷积的定义来解释可分离卷积。若h(x,y)是可分离的卷积核，即可分离为两个一维卷积核相乘的形式： 


h(x,y)=h1(x)h2(y)(352)



则将式(352)中的卷积核代入式(347)的卷积定义中，得出： 


g(x,y)=h(x,y)*f(x,y)=∑+∞m=-∞∑+∞n=-∞h(m,n)f(x-m,y-n)


=∑+∞m=-∞∑+∞n=-∞h1(m)h2(n)f(x-m,y-n)


=∑+∞n=-∞h2(n)∑+∞m=-∞h1(m)f(x-m,y-n)

=∑+∞n=-∞h2(n)(h1*f)(x,y-n)(353)



式(353)中括号内的项为x方向上一维卷积的定义，由下式给出： 



(h1*f)(x,y-n)=∑+∞m=-∞h1(m)f(x-m,y-n)(354)



由式(353)可以证明，对于可分离的卷积核，二维卷积可以通过下式分解为两步一维卷积完成： 


g(x,y)=h2(y)*［h1(x)*f(x,y)］(355)



图342展示了如何分离卷积核，乘号表示矩阵乘法。卷积核分离用矩阵形式表示为h=h1hT2,其中h1和h2为列向量。假设图像的尺寸为M×N，图像块的尺寸为m×n，对于二维卷积，直接卷积的复杂性为O(MNmn)，而可分离卷积的复杂性为O(MN(m+n))。



图342二维卷积核分离为两个一维卷积核相乘的形式


根据矩阵的秩确定卷积核是否可分离，若卷积核是秩1矩阵，则具有可分离性。矩阵的秩等于1也就是核的奇异值分解中有唯一的奇异值。3.6.2.1.1节将介绍的均值和高斯卷积核以及第9章中的一阶差分卷积核均可以分离为两个一维卷积核。例如，Sobel一阶差分模板可以分离为一维高斯平滑卷积核和一维一阶差分卷积核相乘的形式： 

10-1
20-2
10-1=1
2
110-1

3.6.4其他问题
线性空域滤波与频域滤波存在一一对应的关系。空域滤波可以用于非线性滤波，而频域滤波不能用于非线性滤波。非线性空域滤波也是基于邻域的处理，取决于模板对应的邻域内像素，且模板滑过一幅图像的机理与线性空域滤波一致。然而，不能直接利用式(349)计算乘积与求和，因此，非线性滤波不能通过模板卷积实现。

MATLAB中二维离散卷积conv2函数有full、same和valid三种卷积运算的输出模式，这三种模式是对滤波模板移动范围的不同限制。图343以图像平滑为例直观地描述了这三种卷积运算的输出模式，图中正方块表示平滑滤波模板。如图343(a)所示，full是指滤波模板与图像相交开始就进行卷积运算，返回完整的二维卷积。这与3.6.1.1节中介绍的离散卷积运算完全一致。若输入的两个一维离散序列的长度分别为A和B，则线性卷积的离散序列长度等于A+B-1。对于二维离散卷积，输入矩阵的尺寸为A×B和C×D，输出矩阵的尺寸为(A+C-1)×(B+D-1)。如图343(b)所示，same是指当滤波模板的中心位于图像外边界时开始进行卷积运算，对应于线性卷积结果中与输入图像位置相同的中间部分。由于滤波模板的尺寸通常远小于图像的尺寸，这是空域滤波中最常用的模式，这种模式的输出与输入图像的尺寸一致。这种模式的卷积运算涉及图像外边界（border）像素的邻域处理问题。如图343(c)所示，valid是指当滤波模板完全位于图像内部时进行卷积运算，仅输出完全包含滤波模板的图像部分的卷积结果，不考虑模板超出图像外边界的卷积运算。


图343二维卷积的三种输出模式



实现空域滤波的最后一个问题是对于图像外边界像素邻域的处理。对于尺寸为n×n的空域模板，当模板中心距图像外边界为(n-1)/2像素时，该模板正好在图像内部。若模板的中心继续向图像外边界靠近，则模板的行或列将超出图像之外。当模板中心距图像外边界小于(n-1)/2像素时，模板的行或列超出图像之外。图344显示了模板在图像中遍历出现的四种不同的位置，模板中填充的部分超出了图像之外，图中模板的尺寸为5×5。



图344空域滤波的边界问题示意图



MATLAB有四种延拓（padding）方式解决外边界问题，分别为补零、重复（replicate）、对称（symmetric）和循环（circular）方式。补零方式顾名思义，是指通过补零来扩展图像［图345(a)］。重复方式是指通过复制外边界的值来扩展图像［图345(b)］。对称方式是指通过镜像反射外边界的值来扩展图像［图345(c)］。循环方式是指将图像看成二维周期函数的一个周期来扩展［图345(d)］。为了便于观察，图中的白色线框标出了原图像的部分。在空域滤波完成后，从处理后的图像中裁剪出与原图像对应的部分，使处理后的图像与原图像尺寸相等。



图345四种边界延拓方式




3.7空域平滑滤波
图像平滑的作用包括模糊和降噪。图像模糊处理经常用于预处理阶段，例如，为了提取较大的目标，平滑不必要的细节和纹理、桥接直线或曲线的断裂等。图像降噪处理是为了去除或降低图像中的噪声，通常图像具有很强的空域相关性，相邻像素一般具有相同或相近的灰度值，而噪声的特性造成图像灰度的突变，图像平滑处理利用邻域像素的相似性起到降低噪声的作用。空域平滑滤波方法从去除噪声类型的角度可分为加权均值滤波和中值滤波这两类方法，加权均值滤波适用于降低高斯噪声等非脉冲噪声，而中值滤波对滤除脉冲噪声非常有效，本节从这种分类角度讨论这两类方法。空域平滑滤波又可分为线性平滑滤波和非线性平滑滤波，线性平滑滤波能够通过模板卷积实现，除线性平滑滤波以外均属于非线性平滑滤波，非线性平滑滤波由于不能通过模板卷积实现，因此时间开销比线性滤波大。
3.7.1加权均值滤波
本节讨论线性加权均值滤波、边缘保持平滑滤波和非局部均值（nonlocal means）滤波方法。由于邻域像素之间具有高度的相关性，线性加权均值滤波通过对邻域像素进行（加权）平均来平滑图像，边缘保持平滑滤波通过对邻域颜色相似的像素进行加权平均来平滑图像，非局部均值滤波通过对邻域具有相似纹理的像素进行加权平均来平滑图像。
3.7.1.1线性加权均值滤波
如前所述，线性滤波的实现是将空域模板对应邻域内像素值的加权和作为邻域内中心像素的输出响应。线性空域平滑滤波通过（加权）平均模板实现图像滤波。线性平滑模板的权系数全为正值且系数之和等于1，因而不会增加总体灰度程度。在灰度值一致的区域，线性平滑滤波的输出响应不变。线性平滑滤波等效于低通滤波，将在第5章讨论。
最简单的线性平滑滤波所使用的是均值平滑模板，图346(a)所示均值平滑模板的尺寸为3×3，权系数全为1。为了使权系数之和为1，模板系数再乘以归一化因子1/9。通过频率特性分析（见5.7节）可知，均值平滑模板具有旁瓣泄漏效应旁瓣泄漏是指在频域中频谱主瓣内的能量泄漏到旁瓣内。这样，弱信号的主瓣很容易被强信号的旁瓣淹没，造成频谱的模糊和失真。频谱泄漏与旁瓣有关，如果两侧旁瓣的高度趋于零，而使能量相对集中在主瓣，那么可以较为接近真实的频谱。，反映在图像中是振铃效应，且模板尺寸越大，振铃效应越明显。
另一种更重要的平滑模板是加权平均模板，指模板不同位置对应的像素具有不同的权系数。高斯平滑模板是最常用的加权平均模板。对于模板中心对应像素的输出响应，显然中心像素的灰度值比周围像素的灰度值更重要，因而给模板中心对应的像素赋予最大的权系数，随着模板对应的像素到中心位置的距离增大而减小权系数。图346(b)为近似高斯分布的一个简单的3×3高斯平滑模板，中心像素的权系数为4，由于4邻域的像素到中心位置的距离为1，


图346图像平滑模板

而对角邻域的像素到中心位置的距离为2，对角邻域的像素比4邻域的像素到中心位置的距离更远，为此，4邻域像素的权系数为2，而对角邻域像素的权系数为1。如图346(b)所示，加权平滑模板中权系数再乘以归一化因子1/16使模板权系数之和等于1。

更灵活的方式是直接对二维高斯函数关于中心采样来获取任意尺寸的高斯平滑模板。二维高斯函数的表达式为


Gσ(x,y)=12πσ2e-x2+y22σ2(356)



式中，σ为标准差，决定了图像的平滑程度。二维高斯函数具有钟型形状，图347(a)为σ=1.5的二维高斯函数的图形显示。根据与中心像素的距离，通过二维高斯函数来计算该点的权系数，并对所有的权系数进行归一化处理。
二维高斯可以表示为两个一维高斯函数的乘积，即一个关于x的函数和另一个关于y的函数，即


Gσ(x,y)=12πσe-x22σ212πσe-y22σ2(357)



可知高斯平滑滤波具有可分离性，二维卷积过程可以分离为两步一维卷积过程。高斯函数的σ越大，函数越宽，模板尺寸越大，平滑能力越强，如图347(b)所示。
标准差σ是高斯函数的唯一参数，根据选择的高斯函数的标准差σ确定高斯平滑滤波的模板尺寸。如图347(c)所示，服从高斯分布随机变量的取值落在区间(μ-σ,μ+σ)内的概率为68.27%，在区间(μ-2σ,μ+2σ)内的概率为95.45%，在区间(μ-3σ,μ+3σ)内的概率为99.73%。高斯分布取值落在(μ-3σ,μ+3σ)以外的概率小于0.3%，在实际问题中常认为相应的事件是不会发生的，基本上可以将区间(μ-3σ,μ+3σ)看作高斯分布随机变量实际可能的取值区间，称为高斯分布的3σ原则。因此一般选择使用±3σ作为模糊核宽度的参考。MATLAB图像处理工具箱的edge函数中高斯平滑滤波模板的尺寸默认取值为2

3σ

+1，二维高斯滤波imgaussfilt函数默认滤波模板尺寸为2

2σ

+1，这里
·
表示向上取整。



图347高斯函数



例318均值平滑和高斯平滑滤波
图348(a)为一幅尺寸为256×256的灰度图像，分别使用均值平滑模板和高斯平滑模板对该图像进行平滑滤波，通过对高斯函数采样并归一化生成高斯平滑模板，模板尺寸与标准差σ之间的关系设置为2

2σ

+1，图348(b)和图348(c)分别给出了3×3和5×5的高斯模板系数。图348(d)从左到右依次为尺寸为3、5、9和15的均值平滑模板的滤波结果，而图348(e)从左到右依次为尺寸为3、5、9和15的高斯平滑模板的滤波结果，从图中可以看到，平滑模板的尺寸越大，造成边缘越模糊。从图348(d)与图348(e)的比较可以看出，对于小尺寸的模板，均值平滑和高斯平滑的滤波结果没有明显的区别，但是，模板尺寸越大，均值模板的振铃效应越明显。第5章将从频率特性分析振铃效应发生的原因。



图348不同尺寸平滑模板处理示例



3.7.1.2边缘保持平滑滤波
线性均值滤波使用空域卷积对图像进行平滑处理，卷积具有移不变性，因此线性均值滤波在图像降噪同时，图像中边缘和细节的锐度也会受到不同程度的损失。边缘保持平滑滤波考虑不同像素邻域的特性，在梯度较大的边缘处使用各向异性的模糊核，沿着边缘进行平滑，保持图像中的显著边缘，而对图像中纹理、细节以及平坦区域使用各向同性的模糊核，属于非线性滤波。目前具有边缘保持特性的平滑滤波方法有双边滤波（bilateral filtering）、引导滤波（guided filtering）、局部极值（local extrema）平滑、加权最小二乘（weighted least squares）平滑、三边滤波（trilateral filtering）等。
由于双边滤波理论简单且有快速算法，成为目前最为广泛使用的边缘保持滤波方法。Tomasi和Manduchi于1998年提出了双边滤波的理论。双边滤波是一种保持边缘的非迭代平滑滤波方法，它的权系数由空域（spatial domain）S和值域（range domain）R平滑函数的乘积给出。高斯核双边滤波是一种常用的双边滤波方法，即空域和值域平滑函数均是高斯函数。高斯核双边滤波定义为


gb(xp,yp)=1Wp∑(xq,yq)∈SpGσs(‖p-q‖)Gσr(|f(xp,yp)-f(xq,yq)|)f(xq,yq)(358)



其中，gb(xp,yp)为像素(xp,yp)处的输出值； f(xq,yq)为像素(xq,yq)处的输入值，Sp表示以像素(xp,yp)为中心的邻域像素集合； Gσs(·)为空域平滑的高斯函数，σs为空域高斯函数的标准差； ‖p-q‖表示向量p=［xp,yp］T与q=［xq,yq］T之差的范数，即向量表示的像素坐标之间的距离； Gσr(·)为值域平滑的高斯函数，σr为值域高斯函数的标准差； Wp为归一化系数： 


Wp=∑(xq,yq)∈SpGσs(‖p-q‖)Gσr(|f(xp,yp)-f(xq,yq)|)(359)



式(358)表明随着与中心像素的距离和灰度差值的增大，邻域像素的权系数逐渐减小。如图349所示，图349(a)显示了一个存在灰度阶跃的图像块，图349(b)为图349(a)的三维网格图，从图中可以明显看出灰度值的跃变。图349(c)和图349(d)分别为中心像素的高斯滤波函数Gσs(·)的权系数和双边滤波函数Gσs(·)Gσr(·)的权系数。对于与中心像素距离相近，且灰度值相差较小的像素，双边滤波赋予较大的权重； 而对于距离相近，但灰度值相差较大的像素，赋予较小的权重。因此，双边滤波可以很好地保持图像的边缘。由于双边滤波是一种非线性滤波，空域卷积的快速算法不再适用。



图349双边滤波示意图



双边滤波的作用是平滑灰度值变化较小的细节特征，而保持高对比度变化的边缘以及高灰度差的突变。高斯核双边滤波有两个标准差参数： 空域标准差σs和值域标准差σr。较小的空域标准差σs能够滤除小尺度噪声，而保持大尺度特征边缘的锐化程度。当空域标准差σs增大时，增加更远的邻域像素的贡献，滤除大尺度噪声。当σs→∞时，双边滤波退化为值域滤波。较小的值域标准差σr只能平滑方差较小的邻域（如均匀区域），但不会平滑方差较大的邻域（如强边缘）。当值域标准差σr增大时，能够平滑方差较大的邻域，同时会减弱双边滤波的边缘保持能力。当σr→∞时，双边滤波则退化为空域高斯滤波。总而言之，双边滤波需要权衡边缘保持和平滑能力之间的关系。
例319高斯核双边滤波的平滑处理
对于图348(a)所示的图像，使用MATLAB图像处理工具箱的高斯核双边滤波imbilatfilt函数进行边缘保持的平滑处理，空域模板的尺寸与空域高斯函数标准差σs的关系为2

2σs

+1。图350给出了不同参数的双边滤波结果，图350(a)和图350(b)中空域高斯函数的标准差均为σs=1，值域高斯函数的标准差分别为σr=0.01和σr=0.001。通过比较可以看出，值域标准差越小，则边缘保持越好，但是，对于具有高对比度、小尺度特征的细节，不能起到模糊的作用。图350(b)和图350(c)中，值域高斯函数的标准差均为σr=0.001，空域高斯函数的标准差分别为σs=1和σs=3。通过比较可以看出，增大空域高斯函数的标准差，将增大高斯平滑模板的尺寸，使得图像越平滑。



图350双边滤波平滑处理示例



由于噪声的特性也是图像灰度的突变，因此，双边滤波能够有效地应用于图像去噪。对图348(a)所示图像加入均值为0、方差为0.001的加性高斯噪声，生成一幅有噪图像，如图351(a)所示。图351(b)为高斯模板的平滑滤波图像，高斯函数的标准差为σs=1，模板尺寸为5×5。图351(c)为双边滤波去噪结果，空域高斯函数的标准差为σs=1.5，模板尺寸为7×7，值域高斯函数的标准差为σr=0.01。通过比较可以看出，尽管双边滤波中高斯平滑模板的尺寸大于高斯滤波的模板尺寸，然而双边滤波在平滑噪声的同时，能够很好地保持图像的边缘。
 


图351双边滤波去噪处理示例



3.7.1.3非局部均值滤波
前面讨论的线性平滑滤波和边缘保持平滑滤波利用局部相关性，即邻域内像素的灰度值具有相似性和连续性，通过图像邻域像素的加权平均实现图像降噪，称为局部滤波方法。2005年，Buades等提出了非局部均值滤波方法，该方法利用非局部相似性进行降噪滤波处理。
非局部均值滤波方法对搜索窗口中的所有像素计算加权平均，权系数为像素之间的相似度。非局部均值滤波可写为加权和的形式： 


gn(xp,yp)=∑(xq,yq)∈Ωpw(p,q)f(xq,yq)(360)



式中，gn(xp,yp)为像素(xp,yp)处的输出值； f(xq,yq)为像素(xq,yq)处的输入值； Ωp表示像素(xp,yp)的搜索窗口在非局部均值滤波的文献中，这个搜索窗口是整幅图像。； w(p,q)表示像素(xp,yp)与(xq,yq)之间的权系数，权系数之和为1，即


∑(xq,yq)∈Ωpw(p,q)=1(361)



权系数w(p,q)描述图像中像素(xp,yp)与(xq,yq)之间的相似度。由于噪声的存在，使用孤立的像素度量相似性并不可靠，而图像块能够描述复杂的图像结构。式(360)中像素(xp,yp)与(xq,yq)之间的权系数w(p,q)以图像块（像素以及邻域）为单位来度量其相似性，并使用高斯函数作为权系数，则有


w(p,q)=1WpGh/2(‖Vp-Vq‖2,σ)=1Wpe-‖Vp-Vq‖22,σh2(362)



式中，Vp和Vq表示以像素(xp,yp)和(xq,yq)为中心的图像块，并按列表示为向量形式； ‖Vp-Vq‖2,σ表示这两个图像块之间的高斯加权欧氏距离； σ为加权高斯函数的标准差； h称为为衰减因子； Wp为归一化因子，可表示为


Wp=∑(xq,yq)∈Ωpe-‖Vp-Vq‖22,σh2(363)



式(362)表明，权系数w(p,q)度量以像素(xp,yp)和(xq,yq)为中心的图像块之间的相似性。在非局部均值滤波中，对于图像中的各像素(xp,yp)，使用归一化相似度作为权系数，计算其搜索窗口中所有像素(xq,yq)的加权平均值，是像素(xp,yp)的滤波结果。 
高斯加权欧氏距离‖·‖2,σ通过高斯函数对欧氏距离进行加权，赋予不同位置处的像素不同的权重，权重服从高斯分布，与图像块中心越近的像素赋予更大的权重，与图像块中心越远的像素赋予越小的权重。衰减因子h与高斯噪声的标准差ση正相关，即h=κση，其中κ为h与ση的正比例系数。当图像中噪声较大时，增大衰减因子h，相似像素对当前像素的权系数更大，图像更加平滑，细节丢失也更严重； 当图像中噪声较小时，减小衰减因子h，相似像素对当前像素滤波的贡献越小，边缘保持更好，但是保留更多的噪声。


图352非局部均值滤波的解释

非局部均值滤波方法本质上利用图像块的相似性，将相似图像块的加权平均值作为当前图像块的估计，权系数由两个图像块之间的相似度决定。如图352所示，对于以像素p为中心的图像块，在搜索窗口或整幅图像中寻找其相似图像块，对像素p的加权平均贡献实际上主要源于这些相似图像块的中心像素q，图中标出前6个最相似的图像块，分别用q1,q2,…,q6表示，而非相似图像块的权系数很小，对像素p滤波的贡献则很小。非局部处理利用图像块的相似结构不仅能够有效去除图像中的噪声，而且能够有效保持图像的空间细节。

若图像的尺寸为M×N，图像块的尺寸为m×n，搜索窗口的尺寸为S×S，则非局部均值滤波算法复杂度为O(MNmnS2)。当高斯噪声水平级较大时，增大搜索窗口的尺寸。搜索窗口越大，相似图像块越多，去噪的性能就会越好，同时计算量也呈指数增长。实际中需要根据噪声来选取合适的参数。
例320非局部均值滤波去除高斯噪声
对于图351(a)所示的方差为0.001的加性高斯噪声干扰的图像，使用MATLAB图像处理工具箱的非局部均值滤波imnlmfilt函数去除图像中的高斯噪声。在Buades等的实现中，利用欧氏距离和尺寸为m×n的高斯核的卷积实现高斯加权欧氏距离，imnlmfilt 函数省略了这一步以提高计算效率。考虑算法复杂度，图像块的尺寸设置为5×5。图353(a)~(c)给出了不同搜索窗口尺寸S×S和衰减因子h的非局部均值滤波结果。从图353(a)和图353(b)的比较可见，固定衰减因子，将搜索窗口的尺寸从21×21增大到51×51，由于参与加权的像素更多，图像更加平滑，降噪效果更好。从图353(b)和图353(c)的比较可见，固定搜索窗口的尺寸，增大衰减因子h，参与加权的像素权系数更大，图像更加平滑，降噪效果更好。由图中可见，非局部均值滤波有效适用于高斯去噪，且去噪效果优于图351(c)所示的双边滤波的结果。
        


图353高斯噪声图像的非局部均值滤波去噪示例



3.7.2中值滤波相关
均值滤波对于高斯噪声有效，但是对于椒盐噪声的效果一般，本节介绍与中值滤波相关的图像平滑滤波方法，这类方法均属于非线性平滑滤波，主要用于图像降噪处理。统计排序滤波和自适应中值滤波可以有效地去除椒盐噪声，α剪裁均值滤波可用于高斯噪声和椒盐噪声混合的情况。
3.7.2.1统计排序滤波
统计排序滤波是一种简单的非线性滤波方法，对滤除脉冲噪声非常有效。统计排序滤波是将模板对应的邻域内像素的灰度值进行排序，然后将统计排序结果作为模板中心对应像素的输出值。
统计排序滤波中最常用的是中值滤波，中值滤波查找模板对应的邻域内像素值排序的中间值，作为中心像素的输出值。设Sxy表示以像素(x,y)为中心的邻域像素集合，中值滤波在像素(x,y)处的输出g(x,y)为


g(x,y)=median{f(s,t)},(s,t)∈Sxy(364)



式中，median［·］表示中间值查找操作。对一幅图像进行中值滤波的具体过程为： 将模板在图像中遍历，将模板对应的邻域内像素的灰度值排序，查找中间值，将其赋予与模板中心位置对应的图像像素。例如，对于3×3的邻域，其中间值是灰度值排序后的第5个值； 在5×5的邻域中，中间值是第13个值； 在9×9的邻域中，中间值是第41个值。当邻域内具有多个相同灰度值的像素时，可以选取其中任何一个作为中间值。由于中值滤波需要对像素值进行排序，因此它的计算时间一般比线性滤波长，特别是对于较大尺寸的模板。
线性平滑滤波具有低通滤波特性，在降噪的同时也会模糊图像的边缘细节。中值滤波是一种去除噪声的非线性处理方法。一般情况下，中值滤波的结果优于线性滤波。中值滤波具有如下三个特性： ①中值滤波的冲激响应是零，这一性质使其在滤除脉冲噪声方面非常有效，也就是说，当噪声特性未知时，中值滤波对野点具有鲁棒性； ②中值滤波不会改变信号中的阶跃变化，平滑信号中的噪声，又不会模糊信号的边缘，这一性质使其能够很好地适用于图像空域滤波的相关应用； ③中值滤波不会引入新的像素值，因而不会引起灰度上的失真。图354通过一维信号直观描述中值滤波的特性，图354(a)中一维信号在x=5处出现野点，从x=8到x=9时有一个下降沿，图354(b)和图354(c)分别为高斯平滑滤波和中值滤波的结果，由图中可见，对于线性平滑滤波，野点拉高了邻域内的采样值，同时平滑了下降沿，而中值滤波在滤除野点的同时，很好地保持了阶跃信号。



图354一维信号的中值滤波图释



由于脉冲噪声的形式是以孤立的黑白像素叠加在图像上，在图像处理中常称为椒盐噪声。椒盐噪声中，值为0的“黑”像素称为椒噪声，值为255的“白”像素称为盐噪声。中值滤波的直观解释是使模板中心位置对应像素的灰度值更接近它周围像素的灰度值，以此消除孤立的亮点或暗点。为了完全去除椒盐噪声，m×n模板的中值滤波处理要求邻域内孤立亮点或暗点的像素数小于mn/2，也就是要求小于邻域内像素数的一半。
例321中值滤波的去噪处理
图355(a)为一幅由概率为0.1的椒盐噪声污染的电路板图像，椒噪声和盐噪声的概率分别为0.05。使用MATLAB图像处理工具箱的中值滤波medfilt2函数去除图像中的椒盐噪声。图355(b)为3×3均值平滑模板处理的结果，图355(c)为3×3模板的中值滤波处理结果。由于均值平滑模板的响应值是邻域内像素的平均灰度值，具有最大值和最小值的椒盐噪声使像素均值偏向“黑”像素或“白”像素。因此，由图中可见，对于椒盐噪声的情况，中值滤波的处理效果远优于均值平滑滤波。


图355中值滤波和均值滤波对椒盐噪声的平滑滤波比较



中值滤波是统计排序滤波的特例，所谓中值就是邻域内像素值排序中的第

mn/2

个值。统计排序滤波也可以选取排序中其他位置的像素。将模板对应的邻域内像素值按照递增顺序排列，选取第mn个值是最大值滤波。设Sxy表示以像素(x,y)为中心的邻域像素集合，最大值滤波在像素(x,y)处的输出g(x,y)为


g(x,y)=max{f(s,t)},(s,t)∈Sxy(365)



最大值滤波的输出值是邻域内像素的最亮点，因此，最大值滤波能够有效地滤除椒噪声。相反，选取处于第
1个值是最小值滤波，最小值滤波在像素(x,y)处的输出g(x,y)为


g(x,y)=min{f(s,t)},(s,t)∈Sxy(366)



最小值滤波的输出值是邻域内像素的最暗点，因此，最小值滤波能够有效地滤除盐噪声。
例322最大值滤波和最小值滤波的去噪处理
图356(a)为一幅椒噪声污染的电路板图像，椒噪声的概率为0.1，使用最大值滤波去除椒噪声，图357(b)为最大值滤波处理的结果。最大值滤波适用于图像中椒噪声的滤除，但由图像中的黑色金属片可以看出，它同时从暗物体的边缘滤除了一些黑色像素，使得暗物体收缩，而亮物体膨胀。类似地，图357(b)为盐噪声污染的电路板图像，盐噪声的概率为0.1，使用最小值滤波去除盐噪声，图357(d)为最小值滤波处理的结果。最小值滤波适用于图像中盐噪声的滤除，但它同时从亮物体的边缘滤除了一些白色像素，而使亮物体收缩，暗物体膨胀了。第8章数学形态学图像处理中的灰度图像腐蚀和膨胀操作，实际上等效于最小值滤波和最大值滤波。


图356最大值滤波去除椒噪声示例




图357最小值滤波去除盐噪声示例



3.7.2.2自适应中值滤波
中值滤波使用固定尺寸模板对应的邻域内像素的中间灰度值作为模板中心对应像素的输出值，可能引起边缘细节丢失的问题。自适应中值滤波在中值滤波的过程中根据噪声水平自适应地调整模板的尺寸以及输出值，使其能够去除较大概率的脉冲噪声并能够保持图像中的边缘和细节。
设f(x,y)表示输入图像在像素(x,y)处的灰度值，g(x,y)表示输出图像在像素(x,y)处的灰度值，Sxy表示中心在像素(x,y)处的邻域像素集合，zmin、zmax和zmed分别表示邻域像素集合Sxy中的灰度最小值、最大值和中间值，Smax表示允许的最大模板尺寸。
自适应中值滤波包括两个阶段： 阶段A和阶段B，具体的执行过程如下： 
阶段A若zmin<f(x,y)<zmax, 则输出g(x,y)=f(x,y),
否则，转到阶段B。
阶段B若zmin<zmed<zmax, 则输出g(x,y)=zmed,
否则，增大模板的尺寸； 
若模板的尺寸≤Smax，则重复执行阶段B，
否则，输出g(x,y)=zmed。
阶段A判断待处理像素的灰度值f(x,y)是否为脉冲噪声。若zmin<f(x,y)<zmax,则f(x,y)不是脉冲噪声。在这种条件下，输出的灰度值保持不变，即g(x,y)=f(x,y)。通过不改变这些非脉冲噪声的像素来降低边缘和细节的丢失。若f(x,y)=zmin或f(x,y)=zmax，则f(x,y)为脉冲噪声，在这种条件下，转到阶段B。
阶段B计算中值滤波的输出zmed并判断其是否为脉冲噪声。若zmin<zmed<zmax,则zmed不是脉冲噪声。在这种条件下，输出为阶段B中计算的中间值zmed，即g(x,y)=zmed，也就是中值滤波的输出，通过赋予邻域内像素的中间值来消除脉冲噪声。
若阶段B中条件zmin<zmed<zmax不成立，则中值滤波的输出zmed为脉冲噪声，在这种条件下，增大模板的尺寸并重复执行阶段B，继续阶段B的循环直到邻域内像素的中间值zmed并非脉冲噪声，或者达到允许的最大模板尺寸Smax。若达到了最大模板尺寸Smax，则将zmed作为像素(x,y)处的输出g(x,y)。在这种情况下，不能保证该输出值为非脉冲噪声。显然，随着噪声概率的增大，应增大允许的最大模板尺寸Smax。图像受干扰的噪声概率越小，或者允许的最大模板尺寸越大，自适应中值滤波过程发生提前终止的可能性越大。
例323自适应中值滤波的去噪处理
自适应中值滤波建立在中值滤波的基础上，适用于噪声水平较大的椒盐噪声情况。图358(a)为一幅由概率为0.7的椒盐噪声污染的电路板图像，椒噪声和盐噪声的概率分别为0.35，其噪声的概率为图355(a)所示图像中噪声概率的7倍。从图中可以看出，这幅图像具有很高的噪声级。如图358(b)所示，使用尺寸为5×5的模板进行中值滤波，由图中可见，几乎消除了所有的椒盐噪声，仅有少数椒盐噪声仍然残存。增大模板的尺寸，使用尺寸为7×7的模板进行中值滤波，处理结果如图358(c)所示，尽管完全去除了椒盐噪声，然而这样较大尺寸的中值滤波更加明显地造成了图像细节的丢失。注意观察图358(b)和图358(c)中电路板上圆形插孔的边缘失真、黑色连接片的断裂以及说明文字的模糊。
图358(d)为最大模板尺寸Smax=7的自适应中值滤波结果，不仅完全消除了椒盐噪声，而且能够更好地保持图像的锐度和细节，例如，图358(d)中没有发生如上所述的失真、断裂和模糊。可见对于图358(a)这样的高噪声级，自适应中值滤波也能有良好的去噪性能，这证实了自适应中值滤波具有去除较大概率椒盐噪声的优势。另外，由于自适应中值滤波具有提前终止策略，因而它对Smax的取值并不敏感。



图358中值滤波和自适应中值滤波的去噪结果比较



3.7.2.3α剪裁均值滤波
α剪裁均值滤波是一种结合中值滤波和均值平滑滤波的非线性滤波方法，该方法将模板对应的邻域内像素的灰度值进行排序，排除一定数量的处于排序中首尾位置的灰度值，然后计算其余像素的平均灰度值。
设模板尺寸为m×n，Sxy表示中心在像素(x,y)处的邻域像素集合，pmaxi和pmini表示Sxy中第i个最大灰度值和第i个最小灰度值对应的像素，i=1,2,…,d/2，d为［0,mn-1］上的偶数,α剪裁均值滤波在像素(x,y)处的输出g(x,y)为


g(x,y)=1mn-d∑(s,t)∈Sxy＼{pmaxi,pmini}i=1,2,…,d/2f(s,t)(367)



排除邻域像素集合Sxy中d/2个最大灰度值和d/2个最小灰度值剩余mn-d个像素，因此归一化因子为1/(mn-d)。当d=0时，α剪裁均值滤波退化为均值滤波； 当d=mn-1时，α剪裁均值滤波退化为中值滤波。当d取其他值时，α剪裁均值滤波在包括多种噪声的情况下有效，特别是高斯噪声和椒盐噪声混合的情况。
例324α剪裁均值滤波的降噪处理
α剪裁均值滤波建立在均值平滑滤波和中值滤波的基础上，适用于高斯噪声和椒盐噪声混合的情况。图359(a)为一幅由概率为0.2的椒盐噪声和均值为0、方差为0.01的高斯噪声叠加污染的电路板图像。图359(b)、(c)和(d)分别为5×5模板的均值滤波图像、5×5模板的中值滤波和α剪裁均值滤波的处理图像。考虑到椒噪声和盐噪声的概率分别为0.1，25个像素中椒盐噪声的平均数目各为5个，为了充分去除椒盐噪声，设置d=10。从图359(b)可以看到，由于椒盐噪声的存在，均值滤波不能起到平滑椒盐噪声的作用。从图359(c)可以看出，当邻域内像素的中间值是高斯噪声时，中值滤波显然不能去除高斯噪声。从图359(d)可以看出，对于较大的d值，α剪裁均值滤波接近于中值滤波的性能，但是仍然保留了一定的平滑能力。与中值滤波相比，α剪裁均值滤波在降噪方面取得了更好的效果。注意比较在图359(c)和图359(d)所示图像中左上角四个连接片的轮廓，α剪裁均值滤波在排除椒盐噪声之外，能够较好地平滑高斯噪声，使连接片轮廓变得光滑。



图359均值滤波、中值滤波和α剪裁均值滤波的降噪效果比较



3.8空域锐化滤波
在医学成像、遥感成像、手机摄像和视频捕获等成像设备获取图像的过程中，成像机理、成像环境或成像设备可能限制所成图像的清晰度。图像锐化的作用是增强图像中的边缘和细节。然而，图像锐化在增强图像灰度变化的同时，也会放大噪声。线性空域锐化滤波通过差分模板实现图像滤波，二阶差分滤波适用于图像锐化处理。差分算子的响应程度与图像在该点灰度值的突变程度有关，因此，图像锐化使用差分算子。本节将讨论各种锐化模板及其实现算子。
3.8.1微分与差分
函数在某一点的导数是指这个函数在这一点附近的变化率。数学函数导数主要有两点性质： ①在函数值为常数的区域，导数为0； ②函数值变化越快，导数越大。图像可以看成二维函数，在x或y方向的导数即为偏导数。数字图像处理是对离散变量的操作，灰度级的变化发生在两个相邻像素之间。因此，在数字图像处理中，常用有限差分近似偏导数。对于一维离散函数f(x)，一阶导数用一阶差分近似为


Δf(x)=f(x+1)-f(x)(368)



二阶导数用二阶差分近似为


Δ2f(x)=f(x+1)+f(x-1)-2f(x)(369)



对于二维离散函数f(x,y)，将沿两个坐标轴方向计算一阶差分和二阶差分，分别用于近似一阶偏导数和二阶偏导数。
图360(a)为近似阶跃边缘的灰度级剖面曲线，横坐标表示图像中像素的坐标，纵坐标表示像素的归一化灰度值。根据式(368)和式(369)的定义计算水平方向灰度值的一阶差分和二阶差分，如图360(b)所示。在边缘处，一阶差分有一个强的响应，形成一个峰，在二阶差分中边缘的表现是过零点，在过零点也可称为零交叉点。两侧，分别形成一个峰和一个谷，这是二阶差分算子的双响应。在计算差分时会出现负值，通常要归零处理。若对二阶差分响应值取绝对值，则造成双边缘，可以使用二阶差分过零点来定位图像中的边缘。然而，二阶差分一般不直接用于边缘检测。



图360阶跃边缘及其一阶差分和二阶差分示意图



通过图361来分析一阶差分和二阶差分的性质。图361(a)为一幅实际图像的某水平方向灰度级剖面曲线，从①到④依次为实际情况下的斜坡边缘、脉冲边缘、矩形边缘以及阶跃边缘，实际图像的边缘并非理想的情况，由于冲激响应函数的作用，边缘处总存在模糊，且包含噪声，图361(b)为图361(a)的一阶差分和二阶差分曲线。斜坡边缘和阶跃边缘处于图像中不同灰度值的相邻区域之间。对于斜坡边缘，一阶差分实际上是直线的斜率，这条直线为负斜率，一阶差分为负值（约为-2），而常数的斜率为零，由于受噪声干扰，二阶差分在零附近振荡。对于阶跃边缘，一阶差分产生较宽的边缘，而二阶差分更细一些。此外，二阶差分有从正值回到负值的过渡，在边缘图像中表现为双边缘。脉冲边缘主要对应细条状的灰度值突变区域，例如，噪声点、细线，而矩形边缘主要对应较宽的灰度值突变区域。对于矩



图361实际图像的水平方向灰度级剖面曲线及其一阶差分和二阶差分


形边缘，上升沿和下降沿分别可以看成阶跃信号。对于脉冲边缘，二阶差分比一阶差分的响应强，在图像锐化处理中，二阶差分比一阶差分能够更好地增强图像细节。值得注意的是，二阶差分对脉冲边缘的响应强于阶跃边缘。



通过比较一阶差分算子和二阶差分算子的响应，可得出如下结论： ①一阶差分算子对灰度阶跃的响应会产生较宽的边缘，二阶差分算子对灰度阶跃产生双响应，在边缘图像中表现为双线，称为双边效应； ②一阶差分算子一般对灰度阶跃的响应更强，二阶差分算子对孤立点和细线有更强的响应，其中对孤立点比对细线的响应更强。通常情况下，一阶差分算子利用图像梯度突出边缘和细节，主要用于图像的边缘检测； 二阶差分算子是线性算子，通过线性算子提取的边缘和细节叠加在原图像上，主要用于图像的边缘增强。由于本节的内容是图像锐化滤波，后续的内容仅描述二阶差分的拉普拉斯算子，而基于一阶差分的边缘检测内容将在第9章中具体介绍。
3.8.2拉普拉斯算子
拉普拉斯算子是最简单的各向同性二阶差分算子。各向同性滤波器是旋转不变的，即将图像旋转后进行滤波处理的结果与对图像先滤波再旋转的结果相同。二维函数f(x,y)的拉普拉斯变换定义为


2f(x,y)=2f(x,y)x2+2f(x,y)y2(370)



对于二维离散函数，常用二阶差分近似二阶偏导数。因此，数字图像f(x,y)的拉普拉斯变换表示为


2f(x,y)=Δ2xf(x,y)+Δ2yf(x,y)(371)



式中，Δ2xf(x,y)表示x(垂直）方向上的二阶差分，Δ2yf(x,y)表示y(水平）方向上的二阶差分。将式(369)中一维离散函数f(x)的二阶差分定义扩展到二维离散函数f(x,y)，沿x(垂直）方向和y(水平）方向上二阶差分Δ2xf(x,y)和Δ2yf(x,y)的定义为


Δ2xf(x,y)=f(x+1,y)+f(x-1,y)-2f(x,y)(372)

Δ2yf(x,y)=f(x,y+1)+f(x,y-1)-2f(x,y)(373)



将式(372)和式(373)代入式(371)，数字图像f(x,y)的拉普拉斯变换可由下式实现： 


2f(x,y)=［f(x+1,y)+f(x-1,y)+f(x,y+1)+f(x,y-1)］-4f(x,y)(374)



显然式(374)可以用图362(a)所示的4邻域拉普拉斯模板与图像的空域卷积来实现。若考虑8邻域拉普拉斯变换，只需添加两个对角方向的二阶差分： 


Δ2xyf(x,y)=f(x+1,y+1)+f(x-1,y-1)-2f(x,y)(375)

Δ2-xyf(x,y)=f(x-1,y+1)+f(x+1,y-1)-2f(x,y)(376)



新加项的形式与式(372)和式(373)相似，只是其坐标轴的方向沿着对角方向和反对角方向。将两个对角方向加入式(371)拉普拉斯变换的定义中，可得


2f(x,y)=Δ2xf(x,y)+Δ2yf(x,y)+Δ2xyf(x,y)+Δ2-xyf(x,y)


=［f(x+1,y)+f(x-1,y)+f(x,y+1)+f(x,y-1)］+

［f(x+1,y+1)+f(x-1,y-1)+f(x-1,y+1)+f(x+1,y-1)］-

8f(x,y)(377)




图3624邻域和8邻域拉普拉斯模板

这种情况下的拉普拉斯变换用图362(b)所示的8邻域拉普拉斯模板与图像的空域卷积来实现。由于拉普拉斯模板的所有系数之和等于0，在图像中灰度值为常数或者灰度变化平坦的区域，拉普拉斯模板与图像邻域的卷积是0或者几乎为0。因此，拉普拉斯滤波将输出图像的平均灰度值变为0，也就是消除图像频谱中的零频率成分。


由于拉普拉斯算子是差分算子，它将突出图像中灰度的跃变，而清除灰度变化缓慢的区域。拉普拉斯算子是线性算子，将拉普拉斯图像叠加在原图像上，能够同时保持拉普拉斯图像的边缘信息和原图像的灰度信息，它在图像处理中最主要的应用是边缘增强。使用拉普拉斯变换对图像进行锐化滤波可表示为


g(x,y)=f(x,y)-2f(x,y)(378)



式中，2f(x,y)为拉普拉斯图像，g(x,y)为拉普拉斯锐化图像。图363通过一维形式解释拉普拉斯变换用于边缘增强的原理，对于图363(a)所示的图像边缘的水平灰度值剖面图，图363(b)为拉普拉斯变换的结果，也就是二阶差分的结果，对拉普拉斯变换结果的符号取反，如图363(c)所示，将图363(c)加在图363(a)上，即从视觉上增强了边缘的对比度，如图363(d)所示。因此，将原图像减去拉普拉斯图像，可产生锐化滤波图像。



图363拉普拉斯算子边缘增强的一维解释



将式(374)中的2f(x,y)代入式(378)，4邻域拉普拉斯锐化滤波可表示为


g(x,y)=5f(x,y)-［f(x+1,y)+f(x-1,y)+f(x,y+1)+f(x,y-1)］(379)



显然式(379)可以使用图364(a)所示的4邻域拉普拉斯锐化模板与图像的空域卷积来实现。若考虑对角方向的二阶差分，将式(377)中的2f(x,y)代入式(378)，8邻域拉普拉斯锐化滤波可由下式实现： 


g(x,y)=9f(x,y)-［f(x+1,y)+f(x-1,y)+f(x,y+1)+f(x,y-1)］-

［f(x+1,y+1)+f(x-1,y-1)+f(x-1,y+1)+f(x+1,y-1)］(380)




图364拉普拉斯锐化模板

式(380)可以使用图364(b)所示的8邻域拉普拉斯锐化模板与图像的空域卷积来实现。



从另一角度解释图364中的模板。由于拉普拉斯算子是线性算子，可以将f(x,y)看成中心系数为1，其他系数都为0的模板和其自身卷积的结果，而拉普拉斯图像2f(x,y)是图362中拉普拉斯模板与原图像的卷积结果。根据卷积运算具有的线性性质，将这两个模板相加就形成了图364所示的拉普拉斯锐化模板。
例325拉普拉斯图像锐化处理


图3654邻域和8邻域拉普拉斯模板图像锐化示例

对于图365(a)所示的灰度图像，使用拉普拉斯模板对图像进行锐化处理，增强图像的边缘和细节。差分结果中有正值也有负值，通常有三种图像显示差分结果的方式。图365(b)~(d)为这三种方式显示的4邻域拉普拉斯变换的结果。第一种，将差分的数值直接线性映射到8位灰度级显示范围［0,255］，如图365(b)所示，图中大面积平坦的黑色背景区域呈现灰色。第二种，将差分中小于零的数值归零，然后线性映射到［0,255］，如图365(c)所示，仅显示差分为正值的边缘。第三种，对差分的数值作绝对值运算，再线性映射到［0,255］，如图365(d)所示，由于拉普拉斯算子是二阶差分算子，因此图中表现出双边缘。二阶差分算子是线性算子，很容易应用于边缘增强中。图365(e)和图365(f)分别为4邻域和8邻域拉普拉斯锐化图像。由于8邻域拉普拉斯算子比4邻域拉普拉斯算子具有更强的边缘响应，因此图365(f)比图365(e)的边缘更清晰、锐化。

3.9小结
本章从点处理和邻域处理的分类角度，介绍了对图像像素直接处理的空域图像增强方法。在点处理图像增强方面，介绍了灰度级变换、直方图处理以及图像算术运算，其中幂次变换、直方图均衡化和图像相减是主要的图像增强方法。在邻域处理图像增强方面，本章讨论了空域滤波中主要的图像平滑和图像锐化滤波方法，其中线性平滑滤波、中值滤波以及拉普拉斯滤波是重点内容。