冈萨雷斯——《数字图像处理》（三）灰度变换和空间滤波

3.1 背景

空间域处理基于表达式：
g ( x , y ) = T [ f ( x , y ) ] g(x,y)=T[f(x,y)] g(x,y)=T[f(x,y)]
式中， f ( x , y ) f(x,y) f(x,y)是输入图像， g ( x , y ) g(x,y) g(x,y)是输出图像， T T T是点 ( x , y ) (x,y) (x,y)的一个邻域上定义的针对 f f f的算子

具体对于灰度图像来说就是：
s = T ( r ) s=T(r) s=T(r)
其中，r和s分别代表图像处理前后的像素值，这些值通常与灰度变换函数T有关

3.2一些基本的灰度变换函数

图像反转

灰度级在区间 [ 0. L − 1 ] [0.L-1] [0.L−1] 之间的反转函数形式为：
s = L − 1 − r s=L-1-r s=L−1−r
采用这种方式反转图像灰度级，会得到类似照片底片的结果
也就是说黑变白，白变黑；深色变浅色，浅色变深色

对数变换

对数变换的通式是：
s = c l o g ( 1 + r ) s=clog(1+r) s=clog(1+r)
式中c是一个常数，根据下图的灰度变换函数图可以看出，这个变换将输入中的范围较窄的低灰度值映射成输出中的范围较宽的灰度级；相反的，反对数函数（指数函数）的功能就是将高灰度值映射为输出中范围较窄的灰度级

幂律（伽马）变换

幂律变换的形式为：
s = c r γ s=cr^\gamma s=crγ
式中， c c c和 γ \gamma γ 是正常数。如下图可知改变 γ \gamma γ 可以得到一系列变换曲线， γ > 1 \gamma>1 γ>1 时的曲线和 γ < 1 \gamma<1 γ<1的曲线刚好相反， c = γ = 1 c=\gamma=1 c=γ=1时，称为恒等变换

幂律变换的一个应用就是伽马矫正 ，在输入显示器之前，使用变换来预处理图像

另一个应用就是增强图像的对比度：
对于偏暗的图像，需要扩展灰度级，所以取 γ < 1 \gamma<1 γ<1来增强更多的细节
对于偏灰（苍白）的图像，需要压缩灰度级，所以取 γ > 1 \gamma>1 γ>1来使对比度更高，更清晰

分段线性变换函数

与前面讨论的几个方法互补的是分段线性变换函数，
其优点是其形式可以任意复杂；缺点是要求用户输入很多参数

（1）对比度拉伸
对比度拉伸可以扩展图像中的灰度级范围，使其覆盖记录介质或显示设备的整个理想灰度范围。如下图是对比度拉伸的典型变换，点 ( r 1 , s 1 ) , ( r 2 , s 2 ) (r_1,s_1),(r_2,s_2) (r1​,s1​),(r2​,s2​)的位置控制变换函数的形状

（2）灰度级分层
目的就是凸显出图像中的特定灰度区间，共有两种基本方法，如下图
一种是将范围内的所有灰度值显示为一个颜色，其他的是另外一种颜色，如a
另一种方法是将范围内的灰度变亮（或变暗），其他位置保持不变

（3）比特平面分层
例如一幅256级灰度图像中，图像值是由8比特组成，8比特图像可视为由8个1比特平面组成
包括最低有效位比特和最高有效位比特
再进而分析每一个比特平面在图像中的相对重要性

3.3直方图处理

首先引出直方图
令 r k r_k rk​ 表示一幅L级灰度数字图像 f ( x , y ) f(x,y) f(x,y)的灰度，其中 k = 0 , 1 , 2 , . . . , L − 1 k=0,1,2,...,L-1 k=0,1,2,...,L−1
那么这副图像的直方图定义为
h ( r k ) = n k , k = 0 , 1 , 2 , . . . , L − 1 h(r_k)=n_{k} , k=0,1,2,...,L-1 h(rk​)=nk​,k=0,1,2,...,L−1
式中， n k n_k nk​ 是图像中灰度为 r − k r-k r−k的像素的数量.
进而，我们就有 M × N M\times N M×N的图像 f f f的归一化直方图定义为：
p ( r k ) = h ( r k ) M N = n k M N p(r_k)=\frac{h(r_k)}{MN}=\frac{n_k}{MN} p(rk​)=MNh(rk​)​=MNnk​​
我们将这种直方图简单的称为直方图或者图像直方图。对于所有的 p ( r k ) p(r_k) p(rk​)的和总为1,或者说 p ( r k ) p(r_k) p(rk​)的分量是对图像中出现的灰度级的概率的估计

直方图的形状与图像的外观有关

对比发现，像素占据整个灰度级范围并且均匀分布的图像，将具有高对比度的外观和多种灰色调，这也将是后续学习的处理方式期望达到的效果

直方图均衡化

（1）灰度连续
前面我们提到直方图可以看作是图像中出现灰度级的概率估计，也就是意味着直方图就是概率密度函数（PDF）
并且我们最终希望得到的效果是占据整个灰度级并且均匀分布的直方图

也就是说我们将一个概率分布函数 p r ( r ) p_r(r) pr​(r)映射成另外一个概率分布函数 p s ( s ) p_s(s) ps​(s)即可，其中 p s ( s ) p_s(s) ps​(s)是一个均匀概率密度函数（PDF），公式如下：
s = T ( r ) = ( L − 1 ) ∫ 0 r p r ( w ) d ⁡ w s=T(r)=(L-1)\int_0^rp_r(w)\operatorname{d}w s=T(r)=(L−1)∫0r​pr​(w)dw

（2）灰度离散
这种情况下，结合上述的归一化直方图定义将公式表示为：
s k = T ( r k ) = ( L − 1 ) ∑ j = 0 k p r ( r j ) , k = 0 , 1 , 2 , ⋯   , L − 1 s_k=T(r_k)=(L-1)\sum_{j=0}^kp_r(r_j),\quad k=0,1,2,\cdots,L-1 sk​=T(rk​)=(L−1)j=0∑k​pr​(rj​),k=0,1,2,⋯,L−1
以一个例子表示说明灰度离散的直方图均衡化如何实现：
一幅大小64 * 64的3比特图像（L=8），其灰度级是区间 [ 0 , L − 1 ] = [ 0 , 7 ] [0,L-1]=[0,7] [0,L−1]=[0,7]之间的整数

进一步我们可以计算得到 s k s_k sk​ ，当然我们还要四舍五入取整一下到区间 [ 0 , 7 ] [0,7] [0,7]内，得到：
s 0 = 1.33 → 1 s 4 = 6.23 → 6 s 1 = 3.08 → 3 s 5 = 6.65 → 7 s 2 = 4.55 → 5 s 6 = 6.86 → 7 s 3 = 5.67 → 6 s 7 = 7.00 → 7 \begin{aligned}s_{0}=1.33\rightarrow1\quad s_{4}=6.23\rightarrow6\\s_{1}=3.08\rightarrow3\quad s_{5}=6.65\rightarrow7\\s_{2}=4.55\rightarrow5\quad s_{6}=6.86\rightarrow7\\s_{3}=5.67\rightarrow6\quad s_{7}=7.00\rightarrow7\end{aligned} s0​=1.33→1s4​=6.23→6s1​=3.08→3s5​=6.65→7s2​=4.55→5s6​=6.86→7s3​=5.67→6s7​=7.00→7​
可以看出，这里的最终只产生了5个不同的灰度级;并且这里的直方图是PDF的近似，在整个过程中不会产生新的灰度级，因此也很少会出现完全平坦的直方图，前后结果如下

3.4 空间滤波基础

空间滤波在图像处理中应用广泛，这是一个非常大的主题，这里主要介绍其图像增强相关的内容

线性空间滤波的原理

具体来说线性空间滤波就是图像 f f f 和滤波器核 w w w之间执行乘积之和运算
核是一个阵列，其大小定义了运算的邻域；其系数决定了该滤波器的性质
下图说明了使用 3×3 核进行线性空间滤波的原理。在图像中的任何一点 ( x , y ) (x,y) (x,y)处，滤波器的响应 g ( x , y ) g(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 , 1 ) f ( x + 1 , y + 1 ) g(x,y)=w(-1,-1)f(x-1,y-1)+w(-1,0)f(x-1,y)+\cdots+w(0,0)f(x,y)+\cdots+w(1,1)f(x+1,y+1) 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,1)f(x+1,y+1)

当坐标x和y变化时，核的中心逐个像素移动，这样就有了生成滤波后的图像g(也就是响应)
一般来说，大小为 m × n m\times n m×n的核对大小为 M × N M\times N M×N的图像的线性空间滤波可以表示为：
g ( x , y ) = ∑ s = − a a ∑ t = − b b w ( s , t ) f ( x + s , y + t ) g(x,y)=\sum_{s=-a}^a\sum_{t=-b}^bw(s,t)f(x+s,y+t) g(x,y)=s=−a∑a​t=−b∑b​w(s,t)f(x+s,y+t)

如何构建空间滤波器

（1）基于其数学性质。例如计算邻域像素平均值会模糊图像；计算图像局部导数的滤波器会瑞华图像
（2）对形状具有所需性质的二维空间函数进行取样，例如后面会讲的低通滤波器
（3）构建的方法式涉及具有规定频率响应的空间滤波器

3.5平滑空间滤波器

平滑滤波器用于降低灰度的急剧过渡。由于噪声通常以灰度的急剧过渡组成，因此平滑的一个明显应用就是降噪，当然也包括减少无关细节
（1）盒式滤波器核
特点是其系数的值相同（通常为1），如下是一个 3 × 3 3\times 3 3×3的盒式滤波器

（2）低通高斯滤波器核
上述的盒式滤波器在视觉上会产生平滑效果，但是往往会产沿垂直方向模糊图像（方向性）
因此我们想要的核是圆对称的，而高斯核就具有这样的性质
G ( r ) = K e − r 2 2 σ 2 G(r)=Ke^{-\frac{r^2}{2\sigma^2}} G(r)=Ke−2σ2r2​
r 是距离，通常是二维图像中像素点与中心点的欧几里得距离； σ \sigma σ是标准差，控制高斯函数的扩散程度。也就是说距离中心的的距离相同的点值相同

如下就是典型的高斯核，其中前面的系数为各数之和，为归一化所设计

（3）统计排序（非线性）滤波器
这里尤其指中值滤波器
就是将核内覆盖的所有灰度值的中值（排序后位置在中间的值） 代替中心像素值
这一方法比均值滤波更加适合去除椒盐噪声
当然相应的也会有最大值滤波器、最小值滤波器等

[!NOTE] Title
椒盐噪声：就是指椒和盐，是将黑点和白点叠加到图像上

3.6锐化（高通）空间滤波器

与前面的平滑滤波器相对，锐化空间滤波器就是突出图像细节或者增强模糊的部分

在图像处理中，像素值随位置的变化可以用导数来描述：
在离散的图像中我们可以类比连续函数中导数的定义
d f d x = lim ⁡ h → 0 f ( x + h ) − f ( x ) h \frac{\mathrm{d}f}{\mathrm{d}x}=\lim_{h\to0}\frac{f(x+h)-f(x)}{h} dxdf​=h→0lim​hf(x+h)−f(x)​
定义离散的情况下步长 h = 1 h=1 h=1，那么就有得到差分定义为：
∂ f ∂ x = f ( x + 1 ) − f ( x ) \frac{\partial f}{\partial x}=f(x+1)-f(x) ∂x∂f​=f(x+1)−f(x)
∂ 2 f ∂ x 2 = f ( x + 1 ) + f ( x − 1 ) − 2 f ( x ) \frac{\partial^2f}{\partial x^2}=f(x+1)+f(x-1)-2f(x) ∂x2∂2f​=f(x+1)+f(x−1)−2f(x)

经过验证也可以发现：
一阶导数： 反映图像亮度的变化率（梯度），常用来检测边缘。
二阶导数： 反映亮度变化的变化率，常用来检测边缘位置或增强边界。

二阶导数的拉普拉斯算子

可以证明，最简单的各向同性导数算子（核）是拉普拉斯，对于包含两个变量的图像 f ( x , y ) f(x,y) f(x,y),它定义为：
∇ 2 f = ∂ 2 f ∂ x 2 + ∂ 2 f ∂ y 2 \nabla^2f=\frac{\partial^2f}{\partial x^2}+\frac{\partial^2f}{\partial y^2} ∇2f=∂x2∂2f​+∂y2∂2f​
进而代入上述公式得到：
∇ 2 f ( x , y ) = f ( x + 1 , y ) + f ( x − 1 , y ) + f ( x , y + 1 ) + f ( x , y − 1 ) − 4 f ( x , y ) \nabla^2f(x,y)=f(x+1,y)+f(x-1,y)+f(x,y+1)+f(x,y-1)-4f(x,y) ∇2f(x,y)=f(x+1,y)+f(x−1,y)+f(x,y+1)+f(x,y−1)−4f(x,y)
同理就可以得到对应的拉普拉斯核（上式对应最左边，其他的为不同形式）

一阶导数锐化图像——梯度

在图像处理中，一阶导数是用梯度幅度实现的。图像 f f f在坐标 ( x , y ) (x,y) (x,y)处的梯度定义为二维列向量

∇ f ≡ grad ⁡ ( f ) = [ g x g y ] = [ ∂ f / ∂ x ∂ f / ∂ y ] \nabla f\equiv\operatorname{grad}(f)=\begin{bmatrix}g_x\\g_y\end{bmatrix}=\begin{bmatrix}\partial f/\partial x\\\partial f/\partial y\end{bmatrix} ∇f≡grad(f)=[gx​gy​​]=[∂f/∂x∂f/∂y​]

参考前面的导数的定义以及下面的图像，我们可以得到满足一阶导的导数的近似是：

g x = ( z 9 − z 5 ) , g y = ( z 8 − z 6 ) g_x=(z_9-z_5)\quad,\quad g_y=(z_8-z_6) gx​=(z9​−z5​),gy​=(z8​−z6​)
如此一来我们就可以得到罗伯特（roberts）交叉梯度算子的形式为

同理，对于更具有对称性的 3 × 3 3\times 3 3×3的核，得到一阶导的近似为：
g x = ∂ f / ∂ x = ( z 7 + 2 z 8 + z 9 ) − ( z 1 + 2 z 2 + z 3 ) g_{x}=\partial f/\partial x=(z_{7}+2z_{8}+z_{9})-(z_{1}+2z_{2}+z_{3}) gx​=∂f/∂x=(z7​+2z8​+z9​)−(z1​+2z2​+z3​)
g y = ∂ f / ∂ y = ( z 3 + 2 z 6 + z 9 ) − ( z 1 + 2 z 4 + z 7 ) g_{y}=\partial f/\partial y=(z_{3}+2z_{6}+z_{9})-(z_{1}+2z_{4}+z_{7}) gy​=∂f/∂y=(z3​+2z6​+z9​)−(z1​+2z4​+z7​)
对应的就是Sobel算子：