第三章 灰度变换与空间滤波
3.1 背景知识
3.1.1 灰度变换和空间滤波基础
本章节所讨论的图像处理技术都是在空间域进行的。可以表示为下式：
$$g(x, y) = T[f(x,y)]$$
其中$f(x,y)$是输入的图像，$g(x,y)$是处理后的图像，$T$是在点$(x,y)$的邻域上定义的关于$f$的一种算子。
对于图像处理由以下几步组成：邻域原点从一个像素向另一个像素移动，对邻域中的像素应用算子T，并在该位置产生输出。对于边界情况，则根据算子的相关定义处理。
这样的过程称之为空间滤波，其中邻域和预定义操作一起称为空间滤波器（空间掩模）。对图像进行处理，像素点仅仅取决于一个点处的灰度方法称为点处理技术，而取决于邻域的方法称为邻域处理技术。
3.1.2 关于本章的例子
图像增强：对图像进行加工，使其结果对于特定的应用比原始图像更合适的一种处理。
3.2 一些基本的灰度变换函数
r、s分别表示处理前后的像素值，对于8比特环境，包含T的值的一个查找表将有256条记录。下图显示了一些基本的灰度变换函数，主要有线性函数（反转和恒等变换）、对数函数（对数和反对数变换）、幂律函数（n次幂和n次根变换）。
3.2.1 图像反转
$$s=L-1-r$$
反转一副图像的灰度级，得到等效的照片底片。适用于增强嵌入在一幅图像的暗区域中的白色或灰色细节。
3.2.2 对数变换
$$s=clog(1+r)$$
改变换将输入中范围较窄的低灰度值映射为输出中较宽范围的灰度值，我们使用这种类型的变换来扩展图像中的暗像素的值，同时压缩更高灰度级的值。对于傅里叶频谱，用对数变换可以达到一个理想效果。
3.2.3 幂律（伽马）变换
$$s=cr^\gamma$$
我们可以看到$\gamma>1$的值和$\gamma<1$的值所产生的曲线效果完全相反，在$c=\gamma=1$的时候，变成了恒等变换。
上图用于校正幂律响应现象的处理称为伽马校正。伽马变换是一种非线性的改变图像亮度的方法，更适合于人的感觉。
3.2.4 分段线性变换函数
对比度拉伸
它是扩展图像灰度级动态范围的处理，因此它可以跨越记录介质和显示装置的全部灰度范围。例如通过阈值处理函数，可以产生一副二值图像。一般函数是单值且单调递增，这一条件保证了灰度级的次序，从而避免了处理后的图像中产生人为的灰度错误。
灰度级分层
它有两种基本变形方法：一种是将感兴趣的范围内所有灰度值置为一个颜色，而其他灰度值置为另一个颜色，该变换产生一个二值图像；另一种方法是做恒等变换，但在里面感兴趣区域的部分进行灰度值进行改变。下图分别说明了这两种情况。
比特平面分层
将灰度值用比特表示，比如194的灰度值用比特二进制表示为1100010。8比特的图像可以看做8个1比特的平面组成，然后从中抽取相应平面进行重建。例如只抽取第8个比特平面，这相当于进行了二值话处理。
使用比特平面分层，可以分析图像中的每个比特的相对重要性，也可以在图像压缩时起到帮助作用。使用4个高比特平面将允许我们以可接受的细节来重建图像。
3.3 直方图处理
直方图水平轴取对应的灰度值，归一化后纵轴取该灰度值出现的次数比上总像素个数$MN$。若一副图像的像素倾向于占据整个可能的灰度级并且均匀分布，则该图像会有高对比度的外观且展示灰色调的较大变化。
3.3.1 直方图均衡
变换公式：
$$s=T(r), 0\le r\le L-1$$
为了使用这个公式以及它的反变换，我们需要限定两个条件：
- $T(r)$在区间$ 0\le r\le L-1$上为严格单调递增函数。
- 当$0\le r\le L-1$时，$ 0\le T(r)\le L-1$
条件1保证了灰度值变换时不会产生人为缺陷，防止二义性；条件2保证了输入输出范围相同。
随机变量的基本描绘子是概率密度函数（PDF），令$p_r(r)$、$p_s(s)$分别表示随机变量r、s的PDF，由数学推导可知，如果$p_r(r)$和$T(r)$已知，且在感兴趣区域$T(r)$上连续可微，则有重要关系式：
$$p_s(s)=p_r(r)|\frac{dr}{ds}|$$
还有变换函数：
$$s=T(r)=(L-1)\int_{0}^{r}p_r(w)dw$$
公式右边是随机变量$r$的累积分布函数（CDF），根据概率的相关知识，我们可以得到$)\int_{0}^{L-1}p_r(w)dw =1$，这是归一化的要求；还有积分学中的莱布尼茨公式，我们有变上限积分导数的相应公式,$[\int_{-\infty}^{x}f(x)dx]^{'}=f(x)$，有
$$\frac{ds}{dr}=\frac{dT(r)}{dr}=(L-1)\frac{d}{dr}[\int_0^r p_r(w)dw]=(L-1)p_r(r)$$
把这个结果带入上面第一个式子，会有
$$p_s(s)=p_r(r)|\frac{dr}{ds}|=p_r(r)|\frac{1}{(L-1)p_r(r)}|=\frac{1}{L-1} \ 0 \le s \le L-1$$
这是一个均匀概率密度函数，我们可以得出结论，无论$p_r(r)$是什么样子，均一化后$p_s(s)$始终是均匀的。
对于离散值，我们处理概率（直方图的值）与求和来代替处理PDF与积分。公式变为：
$$p_r(r)=\frac{n_k}{MN} \ ,\ k=0,1,2,\cdots ,L-1$$
$$s_k=T(r_k)=(L-1)\sum_{j=0}^{k}p_r(r_j)=\frac{L-1}{MN}\sum_{j=0}^{k}n_j\ ,\ k=0,1,2,\cdots ,L-1$$
在这个公式中，变换$T(r_k)$称为直方图均衡或直方图线性变换。因为直方图是PDF的近似，而且在处理中不允许造成新的灰度级，所以在实际的直方图均衡应用中，很少见到完美平坦的直方图。直方图均衡可以作为自适应对比度增强工具使用。
3.3.2 直方图匹配（规定化）
直方图匹配可以看作将直方图变成人为想要的样子。对于连续变量，我们定义一个有如下特性的随机变量$z$：
$$$G(z)=(L-1)\int_0^z p_z(t)dt=s$
于是会发现$G(z)=T(r)$，这就建立r与z的关系，所以z必须满足：
$$z=G^{-1}[T(r)]$$
这就是r到z的一步变换，在实践中的困难是寻找$T(r)$和$G^{-1}$的有意义表达式。
对于离散函数，我们重写为：
$$G(z_q)=(L-1)\sum_{i=0}^{q}p_z(z_i)=s_k$$
过程是：
- 计算给定图像的直方图$p_r(r)$，并寻找直方图均衡化，其中$s_k$四舍五入到$[0,L-1]$内的整数。
- 对$q=0,1,2,\cdots ,L-1$用上式计算变换函数G中的所有值，也四舍五入到$[0,L-1]$内的整数。
- 对$s_k, k=0,1,2,\cdots ,L-1$用已计算得到的G中的值寻找相对应的$z_k$，并储存相应的s到z的映射。
- 借助上面的过程找到r与z的映射关系。
*直方图规定化在大多时候都是试凑过程。
3.3.3 局部直方图处理
用于增强图像中小区域的细节。由于相当于一个掩模平移，每次更新前一列或者前一行的元素并不会产生"棋盘"效应。
3.3.4 在图像增强中使用直方图统计
对于直方图，我们可以获得类似第二章的数据统计，这里给出r关于其均值m的n阶矩定义：
$$\mu_n(r)=\sum_{i=0}^{L-1}(r_i-m)^n p(r_i)$$
其中二阶矩特别重要，称为灰度方差，通常用$\sigma^2$来表示，开方式标准差$\sigma$，用以度量图像对比度。
在仅处理均值和方差时，可以通过直接取样来估计，而不必计算直方图。这些估计称为取样均值和取样方差：
$$m=\frac{1}{MN}\sum_{x=0}^{M-1}\sum_{y=0}^{N-1}f(x,y)$$
$$\sigma^2=\frac{1}{MN}\sum_{x=0}^{M-1}\sum_{y=0}^{N-1}(f(x,y)-m)^2$$
使用局部直方图的统计信息，可以计算局部的均值和方差，达到局部增强的效果。
$$ g(x,y)=\begin{cases}E\cdot f(x,y) & m_{S_{xy}}\le k_0m_G \ \text{and} \ k_1\sigma_G\le \sigma_{S_{xy}}\le k_2\sigma_G \\ f(x,y) & \text{otherwise}\end{cases}$$
其中$ m_{S_{xy}}\le k_0m_G$，表明这个点是需要增强的，而其他参数又限制了一些细节问题。
3.4 空间滤波基础
线性空间滤波与频率域滤波之间存在一一对应关系。
3.4.1 空间滤波机理
空间滤波器分为线性空间滤波器和非线性空间滤波器，一般滤波器选择奇数维的滤波器，比较方便选取中心点，可以表示为：
$$g(x,y)=\sum_{s=-a}^a\sum_{t=-b}^bw(s,t)f(x+s,y+t)$$
上图展示了如果用模板去遍历图像。
3.4.2 空间相关与卷积
相关是滤波器模板移过图像并计算每个位置的乘积之和的处理；卷积的机理相似，但滤波器首先要旋转180°。需要注意的是第一，相关是滤波器位移的函数；第二，滤波器w与包含有全部0和单个1的函数相关，得到的结果是w的一个拷贝，但旋转了180度。我们将包含单个1其余全是0的函数称为离散单位冲激。
卷积可以理解为使被卷积的函数每一点处都像模板，以离散单位冲激来说，模板先反转180°，然后做卷积再被反转180°，即得到的信号形式就是模板。下面列出相关和卷积的公式：
$w(x,y)$ ☆$f(x,y)= \sum_{s=-a}^a\sum_{t=-b}^bw(s,t)f(x+s,y+t)$
$w(x,y)$ ★$f(x,y)= \sum_{s=-a}^a\sum_{t=-b}^bw(s,t)f(x-s,y-t)$
卷积是线性系统理论的基础，使用相关或者卷积执行空间滤波是优先选择的方法。
3.4.3 线性滤波的向量表示
类似于图像的做法，将一个$m\times m$的矩阵模板看做一个$mm\times 1$的向量，那么卷积或者相关的模板响应R可以看作：
$$R=w_1z_1+w_2z_2+ \cdots +w_{mm}z_{mm}=\boldsymbol{w}^T \boldsymbol {z}$$
3.4.4 空间滤波器模板的产生
我们使用线性滤波所能做的所有事情是实现乘积求和的操作。一些实例我们等会在介绍。
3.5 平滑空间滤波器
平滑滤波器用于模糊处理和降低噪声。模糊处理常用语预处理操作，如提取细节，以及桥接直线或者曲线的缝隙。
3.5.1 平滑线性滤波器
有时也称为均值滤波器，也可以归入低通滤波器，其应用主要是降低噪声。
在所有方向都等效，所有系数都相等的空间均值滤波器，如$R=\frac{1}{9} \sum_{i=1}^{9}z_i$，也称为盒状滤波器。其通式可以表示为：
从模板上来看可能更为直观：
空间均值处理的重要应用是为了对感兴趣的物体得到一个粗略的描述而模糊一副图像，之后再利用阈值处理并基于物体亮度来消除某些物体的操作是很典型的。
3.5.2 统计排序（非线性）滤波器
其中最有名的是中值滤波器，它是将像素邻域内灰度的中值代替该像素的值，这种处理方式在处理椒盐噪声时特别有效。按照此思路，也可以取最大值名为最大值滤波器，取最小值称为最小值滤波器。下图展示了用中值滤波器处理的结果（扫描版的PDF上带图非常不清楚…），还是能看到椒盐噪声被有效地消除：
3.6 锐化空间滤波器
锐化处理主要目的是突出灰度的过度部分，上面介绍的均值处理对于与数学中的积分思想，这里介绍的锐化处理对应于数学中的微分（差分）思想。
3.6.1 基础
对于一阶微分的定义：
- 在恒定灰度区域的微分值为零。
- 在灰度台阶或斜坡处微分值非零。
- 沿着斜坡的微分值非零。
对于二阶微分的定义：
- 在恒定灰度区域的微分值为零。
- 在灰度台阶或斜坡处微分值非零。
- 沿着斜坡的微分值非零，若是斜率不变的斜坡，其二阶微分值为0。
下面我们给出对于离散变量，一阶差分和二阶差分的定义式：
$$\frac{\partial f}{\partial x}=f(x+1)-f(x)$$
$$\frac{\partial^2f}{\partial x^2}=f(x+1)+f(x-1)-2f(x)$$
下图表明了上述的具体含义。注意下图中二阶差分的零交叉点（Zero crossing），零交叉点对于边缘的定位非常有用。
二阶微分在增强细节方面比一阶微分好得多，这是一个适合锐化图像的理想特征。
3.6.2 使用二阶微分进行图像锐化--Laplace算子
最简单的各向同性微分算子是Laplace算子，一个二维图像的函数$f(x,y)$的Laplace算子定义为
$$\nabla^2f=\frac{\partial^2f}{\partial x^2}+\frac{\partial^2f}{\partial y^2}$$
则离散标量的Laplace算子是：
$$\nabla^2f(x,y)=f(x+1,y)+f(x-1,y)+f(x,y+1)+f(x,y-1)-4f(x,y)$$
上述模板是以90°为增量进行旋转的一个各向同性结果，还可以将增量改成45°，如下图表示：
在实践中注意中心正负的问题，如果使用定义负的中心系数，那么必须将原图减去经Laplace变换后的图像得到锐化结果。
3.6.3 非锐化掩蔽和高提升滤波
步骤有以下几步：
- 模糊源图像
- 从原图像中减去模糊图像（产生的差值图像称为模板）
- 将模板加到原图像上。
令$\overline{f}(x,y)$表示模糊图像，那么首先我们得到模板：
$$g_{mask}(x,y)=f(x,y)-\overline{f}(x,y)$$
然后在原图像上加上该模板的一个权重部分：
$$g(x,y)=f(x,y)+k\times g_{mask}(x,y) \ \ k\ge 0$$
上式中，当k=1时，我们得到上面定义的非锐化掩蔽，当k>1时，该处理称为高提升滤波，当k<1时，则不强调非锐化模板的贡献。
需要注意的是，如果k足够大的时候，负值将导致边缘周围有暗的晕轮，会产生不好的效果。
3.6.4 使用一阶微分对（非线性）图像锐化—梯度
一个函数$f$在$f(x,y)$的梯度是定义为二维列向量（这是一个矢量，给定了一个方向）：
然后其幅度值为$M(x,y)=mag(\nabla f)=\sqrt{g_x^2+g_y^2} \approx |g_x|+|g_y|$，其中$M(x,y)$是与原图像大小相同的图像，该图像通常称为梯度图像。后面的表达式仍保留了灰度的变换，但丢失掉了其各向同性。
下面介绍两种算子：一阶Roberts交叉梯度算子和二阶Soble算子：
在$3\times 3$的区域图像中，对于Roberts算子，$g_x=(z_9-z_5)$，$g_y=(z_8-z_6)$；对于Soble算子，$g_x=(z_7+2z_8+z_9)-(z_1+2z_2+z_3)$，$g_y=(z_3+2z_6+z_9)-(z_1+2z_4+z_7)$。
使用梯度进行边缘增强，可以用以突出灰度图像中看不见的小斑点，在下面的例子中甚至可以去除灰度不变或变换缓慢的图案阴影。
3.7 混合空间增强法
本书所提供例子的策略是：用Laplace突出图像中的小细节，然后使用梯度法突出边缘。平滑过得梯度图像用于掩蔽Laplace图像，最后使用灰度变换来增大图像的灰度动态范围。其中降低噪声可以使用中值滤波器或者使用原图像梯度操作的平滑形式形成的一个模板。
图像请直接参考Digital Image Processing （3rd Edition）Page 192。
3.8 使用模糊技术进行灰度变换和空间变换
3.8.1 引言
隶属度函数$\mu(z)$：用以描述一个元素是否属于集合的模糊程度，若是一个阶跃函数，可以认为是一个我们所了解的"干脆的"集合，若是一个分段函数，则可以看做一个模糊集合：
3.8.2模糊集合论原理
模糊集合是一个有$z$值和（赋予$z$成员等级的）相应隶属度函数组成的序对，即
$$A=\{z, \mu_A(z)|z\in Z\}$$
其中隶属度函数是关键，表明了元素$z$到集合的一种对应关系。如果隶属度函数仅有0,1两个值，那么模糊集合退化为"干脆的"集合。下面这张图展示了模糊集合间的一些关系：
有一些常用的隶属度函数，如三角形、梯形、$\sum$型、S型、种型、截尾高斯型等，之后用到了会另作介绍。
3.8.3 模糊集合的应用（略）
3.8.4 使用模糊集合进行灰度变换
对于一副灰度图像，描述它的某个区域"暗的"、"灰的"、"亮的"都是利用模糊的概念，那么我们将给予这三种情况定义三种模糊集合：
那么对于任何输入$z_0$，输出$v_0$为
从书中给予的例子中来看，图像的细节部分得到了比较好的保留，但代价是计算量大大增加。
3.8.5 使用模糊集合进行空间滤波（略）
（MATLAB 版）
3.1背景知识
我们在处理中所使用的邻域范围是以$(x,y)$为中心的一点，处理的结果来代替该点的灰度值。我们的邻域中心从左上角开始一个像素一个像素的遍历整幅图像，从而得到整幅图像的相应修改。
3.2 强度变换函数
3.2.1 imadjust 函数
imadjust – Adjust Image intensity values or colormap.
用法：g = imadjust(f, [low_in high_in], [low_out high_out], gamma)
f是输入图像，格式可以为double、Uint8或者Uint16，输出格式跟输入格式一样，其中low_in 、high_in、low_out 、high_out分别是输入输出的阈值标准，范围是0-1，如果格式是Uint8，那么阈值会是其乘以255的结果；如果high_out<low_out，那么输出的结果将会是反色的图片，如下图；其中只是将图像反色的话，可以使用函数imcomplement
imcomplement – Complement image
用法：g = imcomplement(f);
左边的一张图是原图，右边的图示反色后的结果，语句为g=imadjust(f, [0 1], [1 0])；通过对这四个参数进行修改，可以达到拓展灰度级的作用。最后一个gamma值是对映射方式进行的改变。当gamma=1
时，是一种线性的映射；当$gamma\neq 1$时，映射方式更像幂次函数，关系如下：
3.2.2 对数和对比度拉伸变换
对数的使用是借助log函数（以自然常熟e为底数，还有log10， log2），表达式为$g=c*log(1+double(f))$
变换的前提是将一副图像变为uint8的灰度图像，使用语句gs = im2uint8(mat2gray(g))。
下面的例子中将我们关心的灰度级区域的对比度拉伸，其他不感兴趣将其置为0或者255，，下面这个函数也叫阈值函数：
$$s=T(r)=\frac{1}{1+(m/r)^E}$$
其中r代表输入图像的强度，s代表输出图像的强度，E代表函数的斜率，那么在MATLAB中的代码写为：
g = 1./(1 + (m./(double(f) + eps)).^ E)
其中加上eps是防止出现分母为0的情况。
下图展示了对比度拉伸函数和阈值函数在E上的不同（手绘，请见谅）：
3.2.3 一些强度变换的实用函数
下面将介绍两个自己定义的函数，在此之前先做一些准备工作：
控制输入和输出函数的变量数量
nargin- Number of function input arguments
nargout - Number of function output arguments
nargchk - Validate number of input arguments
用法：msg = nargchk(low, high, number)
error - Display message and abort function
可以使用语句来进行输入变量数的控制：
error(nargchk(2, 3, nargin))
下面介绍两个通用使用参数的两个列表
varargin - Variable-length input argument list
varargout - Variable-length output argument list
一个强度变换函数
Function g = intrans(f, varargin)
可以实现翻转、log变换、gamma、变换，和对比度伸缩变换的强度变化，这里面用到了另一个函数changeclass，这个函数可以将图像的类型转变成另外一个，可能的类型为uint8, uint16, double.
函数表形附件提供，我们利用对比度伸缩变换，由下图可以看到区别，命令语句是：
g = intrans(f, 'stretch', mean2(im2double(f)), 0.9);
Imshow(f), Figure, imshow(g);
灰度伸缩变换函数
g = gscale(f, method, low, high)
该函数的作用是将一个图像的灰度级范围扩展到8-bit或者16-bit等想要的范围。
3.3 直方图处理
3.3.1 直方图处理函数
理论如上，主要用到对的函数有一下几个：
Imhist - Display histogram of image data
用法：h = imhist(f, b) 其中b是容器的个数，对于uint8图像，b的默认值是256 。归一化可以将结果比上总像素数（numel(f)）
bar - Bar graph
用法：bar(horz, v, width)，v 和 horz是等长度的行向量，width是0~1的数值。
stem - Plot discrete sequence data
用法：stem(horz, v, 'color_linestyle_marker', 'fill')，color_linestyle_marker由三个部分组成，参见下表
plot - 2-D line plot
用法：plot(horz, v, 'color_linestyle_marker')
3还有一些对坐标轴、图像等处理的辅助函数，详情可看书中介绍。
3.3.2 直方图均衡化
histeq – Enhance contrast using histogram equalization
用法：g = histeq(f, nlev)，其中nlev的默认值是64。
3.3.3 直方图匹配（规定化）
函数依旧是histeq，用法：g = histeq(f, hspec)，其中hspec为规定的直方图，即一个规定值的行向量。
3.4 空间滤波
邻域处理由如下步骤组成：
- 选区中心点；
- 对预先定义的关于中心点的领域内的像素执行操作；
- 定运算结果为该点处的相应；
- 对图像中每一点重复该处理。
3.4.1 线性空间滤波
原理如上介绍，下面介绍通用函数：
Imfilter – N-D filtering of multidimensional images
用法：g = imfilter(f, w, filtering_mode, boundary_options, size_options)
其中f是图像，w是滤波模板，后面三个选项如下：
3.4.2 非线性空间滤波
常使用的函数有nlfilter和colfilt，前者执行二位操作，后者按列的形式组织数据，并且后者速度更快。
Colfilt – columnwise neighborhood operations
用法： g = colfilt (f, [m n], 'sliding', fun)
其中m，n表示滤波区域的维数，'sliding'表示区域内图像f中逐像素华东，fun是一个函数句柄。
在处理边界问题的时候，需要显示的填充图像，可以使用函数padarray
Padarray – Pad array.
用法：fp = padarray(f, [r c], method, direction)
r和c是用于填充f的行和列，method和direction详见下表：
3.5 图像处理工具箱的标准空间滤波器
3.5.1 线性空间滤波器
Fspecial - Create predefined 2-D filter
用法： w = fspecial('type', parameters)
其中type由下表选择，具体的参数请详见文档：
3.5.2 非线性滤波器
Ordfilt2 - 2-D order-statistic filtering
用法：g = ordfilt2(f, order, domain)，这也称为统计排序滤波器，通过邻域的排序组合，用第order个元素代替f中的元素，邻域由domain指定。例如domain为ones(m, n)，则最大值滤波order可以为m*n，最小值滤波order可以为1，中值滤波order为(m*n)/2 。工具箱中还专门提供一个中值滤波器：
medfilt2 - 2-D median filtering
用法：g = medfilt2(f, [m n], padopt) 其中padopt可以为三种边界填充选项之一：zeros、symmetric、indexed。