图像处理学习笔记-03-灰度变换与空间滤波

术语空间域指图像平面本身，这类图像处理方法直接以图像中的像素操作为基础；空间域处理主要分为灰度变换和空间滤波两类，灰度变换在图像的单个像素上操作，主要以对比度和阈值处理为目的；空间滤波涉及改善性能的操作。

背景知识

本章讨论的空间域处理可由下式表示：
$$g(x,y)=T[f(x,y)]$$
其中$f(x,y)$是输入图像，$g(x,y)$是处理后的图像，$T$是在点$(x,y)$的邻域上定义的关于$f$的一种算子

一些基本灰度变换函数

假设$r,s$分别代表处理前后的像素值，$s=T(r)$，其中$T$是把像素值$r$映射到像素值$s$的一种变换；

• 图像反转：灰度级范围为$[0,L-1]$
  $$s=L-1-r$$
  增强一幅图像的暗区域中的白色或灰色细节
• 对数变换：$c$是一个常数
  $$s=c\log (1+r)$$
  该变换将输入中范围较窄的低灰度值映射为输出中较宽范围的灰度值，使用这种变换来扩展图像中的暗像素的值，同时压缩更高灰度级的值，也就是将暗区域的细节显现出来
• 幂律(伽马)变换：$c,\gamma$为正常数
  $$s=c{r}^{\gamma }$$
  当$\gamma =1$为恒等变换，$\gamma >1$增加了亮区域的细节，使原图片变暗，$\gamma <1$增加了暗区域的细节，使原图片变亮，效果和对数变换类似，用于校正幂律响应现象的处理称为伽马校正
• 分段线性变换函数
  对比度拉伸
  
  灰度级分层
  
  比特平面分层
  

直方图处理

灰度级范围为$[0,L-1]$的数字图像的直方图是离散函数$h(r_k)=n_k$，其中$r_k$是第$k$级灰度值，$n_k$是图像中灰度为$r_k$的像素个数，假设$M,N$分别是图像行和列的维数，那么归一化后的直方图为$p(r_k)=n_k/MN，k=0,1,\cdots ,L-1$，简单的说$p(r_k)$是灰度级$r_k$在图像中出现的概率的一个估计

直方图均衡

用$r$表示待处理图像的灰度，通常我们假设$r$的取值区间为$[0,L-1]$，$r=0$表示黑色，$r=L-1$表示白色，对于输入图像中每个具有$r$值的像素值产生一个灰度值$s$：
$$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$

我们知道高对比度的图片的直方图占据整个灰度级而且分布均匀
那么直方图均衡的转换函数为：
$$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,\cdots ,L-1$$
上式中的变换称为直方图均衡或直方图线性变换
证明：
首先变换后的变量$s$的$PDF$概率密度函数可以由下式得到：
$$p_s(s)=p_r(r)\left|\frac{dr}{ds}\right|$$
当变换函数如下时可以证明$s$的概率密度函数是均匀分布函数，该变换公式满足开始的两个条件
$$s=T(r)=(L-1){\int }_0^r{p}_r(w)dw$$
又有：
$$\frac{ds}{dr}=\frac{dT(r)}{dr}=(L-1)\frac{d}{dr}\left[{\int }_0^r{p}_r(w)dw\right]=(L-1)p_r(r)$$
可以推出：
$$p_s(s)=p_r(r)\left|\frac{dr}{ds}\right|=p_r(r)\left|\frac{1}{(L-1)p_r(r)}\right|=\frac{1}{L-1},0\le s\le L-1$$
是一个均匀概率密度函数，也就是经过变换之后的图像的灰度值的分布是均匀的，灰度值是离散取值的：
$$p_r(r_k)=\frac{n_k}{MN}，k=0,1,\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,\cdots ,L-1$$

直方图匹配

目的：有的时候采用均匀直方图的基本增强不是最好的方法，希望处理后的图像具有规定的直方图形状可能更有用
直方图匹配或直方图规定化：产生处理后有特殊直方图的方法

• 方法：
  首先将原图均匀化：
  $$s=T(r)=(L-1){\int }_0^r{p}_r(w)dw$$
  之后将我们需要的直方图形式所对应的输出的灰度值$z$均匀化，这个$z$也就是我们需要求的：
  $$G(z)=(L-1){\int }_0^z{p}_z(t)dt=s$$
  之后通过反函数求得$z$
  $$z=G^{-1}[T(r)]=G^{-1}(s)$$
  于是就得到了我们需要的输出的灰度值
• 一个例子：
  假设一幅图像的灰度$PDF$为$p_r(r)=2r/(L-1{)}^2,0\le r\le L-1$，对于其他$r$值有$p_r(r)=0$，寻找一个变换函数使得产生的灰度$PDF$是$p_z(z)=3z^2/(L-1{)}^3,0\le z\le L-1$，对于其他的$z$值有$p_z(z)=0$
  首先对原图像进行直方图均衡变换，得到灰度值为$s$的均匀$PDF$图像：
  $$s=T(r)=(L-1){\int }_0^r{p}_r(w)dw=\frac{2}{L-1}{\int }_0^{r}wdw=\frac{r^2}{L-1}$$
  对希望输出的图像进行直方图均衡变换，也得到灰度值为$s$的均匀$PDF$图像：
  $$G(z)=(L-1){\int }_0^z{p}_z(w)dw=\frac{3}{(L-1{)}^2}{\int }_0^z{w}^{2}dw=\frac{z^3}{(L-1{)}^2}$$
  之后求得反函数，得到$s$到$z$的映射：
  $$z=[(L-1{)}^{2}s{]}^{\frac{1}{3}}$$
  有$s=\frac{r^2}{L-1}$，带到上式可得$r$到$z$的映射：
  $$z=[(L-1{)}^{2}s{]}^{\frac{1}{3}}={\left[(L-1{)}^2\frac{r^2}{L-1}\right]}^{\frac{1}{3}}=[(L-1)r^2{]}^{\frac{1}{3}}$$
  而因为灰度值是离散的，所以可以不需要计算$G$的反变换，我们可以首先计算$G(z_q)=(L-1){\sum }_{i=0}^q{p}_z(z_i)$得到$z$和$s$的一个对应形成一张表，之后对于输入的图像的每一个灰度值$r$计算均衡化之后的$s$值，在表中查找对应的$z$值

局部直方图处理

上面两个小节的图像处理都是基于整个图像而言的，还有一种情况是，每一个像素基于邻域来进行局部增强，在每个位置，计算邻域中的点的直方图，例如直方图均衡化或者规定化变换函数，在这个处理的过程中不断移动邻域中心

在图像增强中使用直方图统计

$r$的均值为$m$：
$$m=\sum _{i=0}^{L-1}{r}_{i}p(r_i)$$
$r$关于其均值的$n$阶矩定义为：
$${\mu }_n(r)=\sum _{i=0}^{L-1}(r_i-m{)}^{n}p(r_i)$$
其中二阶矩称为灰度方差，通常使用${\sigma }^2$表示，是图像对比度的度量，还可以直接用下面的公式计算：
$$\begin{gathered} 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 \end{gathered}$$
局部均值和方差是根据图像中每一像素的邻域内的图像特征进行改变的基础，$S_{xy}$表示规定大小的以$(x,y)$为中心的邻域，也就是子图像，$p_{S_{xy}}$是区域$S_{xy}$的直方图，该领域中像素的均值为：
$$m_{S_{xy}}=\sum _{i=0}^{L-1}{r}_i{p}_{s_{xy}}(r_i)$$
邻域中像素的方差由下式给出：
$${\sigma }_{S_{xy}}^2=\sum _{i=0}^{L-1}(r_i-m_{S_{xy}}{)}^2{p}_{S_{xy}}(r_i)$$
下面是一个局部图像增强的例子

空间滤波基础

• 空间滤波器由一个邻域，对该邻域包围的图像像素执行的预定义操作组成，滤波产生一个新像素，新像素的坐标等于邻域中心的坐标，像素的值是滤波操作的结果；
• 一般来说，假设$m=2a+1,n=2b+1$，使用大小为$m\times n$的滤波器对大小为$M\times N$的图像进行线性空间滤波，可由下式表示：
  $$g(x,y)=\sum _{s=-a}^a\sum _{t=-b}^{b}w(s,t)f(x+s,y+t)$$
• 相关和卷积：卷积需要将滤波器旋转180°，相关不需要旋转，如果滤波器是对称的，那么相关和卷积的结果相同
  相关：
  $$g(x,y)=\sum _{s=-a}^a\sum _{t=-b}^{b}w(s,t)f(x+s,y+t)$$
  卷积：一个函数与单位冲激的卷积，相当于在单位冲激的位置复制该函数
  $$g(x,y)=\sum _{s=-a}^a\sum _{t=-b}^{b}w(s,t)f(x-s,y-t)$$
  
• 线性滤波的向量表示：$w$是一个大小为$m\times n$的滤波器的系数，$z$为滤波器覆盖的相关图像的灰度值
  $$R=w_1{z}_1+w_2{z}_2+\cdots +w_{mn}{z}_{mn}=\sum _{k=1}^{mn}{w}_k{z}_k=w^{T}z$$

平滑空间滤波器

目的：用于模糊处理和降低噪声

平滑线性滤波器

平滑线性空间滤波器的输出是包含在滤波器模板邻域内的像素的简单平均值，有时也称为均值滤波器，可以归入低通滤波器,由于典型的随机噪声由灰度级的急剧变化组成，因此，常见的平滑处理应用就是降低噪声，但是由于图像边缘也是由图像灰度尖锐变化带来的特性，所以均值滤波处理还是有着边缘模糊的负面效应，所有系数都相等的空间均值滤波器称为盒状滤波器
$$R=\frac{1}{9}\sum _{i=1}^9{z}_i$$
减少边缘模糊的一种方法：加权均值滤波器，也就是权重值随着离中心点的距离增大而减小，一幅$M\times N$的图像经过一个大小为$m\times n$的加权均值滤波器的过程可由下式给出：
$$g(x,y)=\frac{{\sum }_{-a}^a{\sum }_{-b}^{b}w(s,t)f(x+s,y+t)}{{\sum }_{-a}^a{\sum }_{-b}^{b}w(s,t)}$$

统计排序(非线性)滤波器

这种滤波器的响应以滤波器包围的图像区域中所包含的像素的排序为基础，例如中值滤波器，中值滤波器对处理椒盐噪声非常有效，这种噪声是以黑白点的形式叠加在图像上面的

锐化空间滤波器

锐化处理的主要目的就是突出灰度的过渡部分，增强了边缘和其他突变(例如噪声)，削弱了灰度变化缓慢的区域，本节将讨论由数字微分来定义和实现锐化算子的各种方法

• 一阶微分的基本定义：
  $$\frac{\partial f}{\partial x}=f(x+1)-f(x)$$
• 二阶微分的基本定义：
  $$\frac{{\partial }^{2}f}{\partial x^2}=f(x+1)+f(x-1)-2f(x)$$

使用二阶微分进行图像锐化-拉普拉斯算子

• 各向同性滤波器：这种滤波器和滤波器作用的图像的突变方向无关，也就是说各向同性滤波器是旋转不变的，即将原图像旋转后进行滤波处理得到的结果和先对图像滤波然后再旋转的结果相同
• 最简单的各向同性微分算子是拉普拉斯算子，一个二维图像函数$f(x,y)$的拉普拉斯算子定义为
  $${\nabla }^{2}f=\frac{{\partial }^{2}f}{\partial x^2}+\frac{{\partial }^{2}f}{\partial y^2}$$
  以离散形式描述，在$x$方向上：
  $$\frac{{\partial }^{2}f}{\partial x^2}=f(x+1,y)+f(x-1,y)-2f(x,y)$$
  在$y$方向上：
  $$\frac{{\partial }^{2}f}{\partial y^2}=f(x,y+1)+f(x,y-1)-2f(x,y)$$
  由上式可得，拉普拉斯算子是：
  $${\nabla }^{2}f(x,y)=f(x+1,y)+f(x-1,y)+f(x,y+1)+f(x,y-1)-4f(x,y)$$
  上式可以用滤波模板来实现：第一幅是以90°为增量进行旋转的一个各向同性结果，第二幅是以45°为增量的
  
  使用拉普拉斯对图像进行增强的基本方法如下：如果使用$a,b$中的滤波器，$c=-1$，另外两幅的话$c=1$
  $$g(x,y)=f(x,y)+c[{\nabla }^{2}f(x,y)]$$

非锐化掩蔽和高提升滤波

非锐化掩蔽的处理过程：

• 模糊原图像
• 从原图像中减去模糊图像(产生的差值图像称为模板)
• 将模板加到原图像上

令$\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=1$时，得到上面定义的非锐化掩蔽，$k>1$时，该处理称为高提升滤波，$k<1$则不强调非锐化模板的贡献

使用一阶微分对(非线性)图像锐化-梯度

图像处理中的一阶微分是用梯度幅值来实现的，对于函数$f(x,y)$，$f$在坐标$(x,y)$处的梯度定义为二维向量
$${\nabla }_f\equiv grad(f)=\left[\begin{array}{c}{g}_x\\ g_y\end{array}\right]=\left[\begin{array}{c}\frac{\partial f}{\partial x}\\ \frac{\partial f}{\partial y}\end{array}\right]$$
该向量具有重要的几何特性，即它指出了在位置$(x,y)$处$f$的最大变化率的方向，向量${\nabla }_f$的幅度值表示为$M(x,y)$:
$$M(x,y)=mag({\nabla }_f)=\sqrt{g_x^2+g_y^2}$$
$M(x,y)$是与原图像大小相同的图像，该图像称为梯度图像，梯度向量的分量是微分，所以他们是线性算子，然而该向量的幅度不是线性算子，偏微分不是旋转不变的，梯度向量的幅度是旋转不变的，在某些实现中，用绝对值的近似平方和平方根操作更适合计算，近似式在旋转90°的倍数的时候是各向同性的，近似式也是非线性的：
$$M(x,y)\approx |g_x|+|g_y|$$
假设使用交叉差分：
$$\begin{gathered} g_x=(z_9-z_5) \\ {g}_y=(z_8-z_6) \end{gathered}$$
结合上式得到梯度图像：
$$M(x,y)={\left[(z_9-z_5{)}^2+(z_8-z_6{)}^2\right]}^{\frac{1}{2}}$$
如使用近似平方根：
$$M(x,y)=|z_9-z_5|+|z_8-z_6|$$
能实现上述所需偏微分项的滤波器称为罗伯特交叉梯度算子，见图$b,c$，但是偶数尺寸的模板难以实现，使用$3\times 3$的模板：
$$\begin{gathered} g_x=\frac{\partial f}{\partial x}=(z_7+2z_8+z_9)-(z_1+2z_2+z_3) \\ {g}_y=\frac{\partial f}{\partial y}=(z_3+2z_6+z_9)-(z_1+2z_4+z_7) \end{gathered}$$
计算梯度幅值：
$$M(x,y)\approx |(z_7+2z_8+z_9)-(z_1+2z_2+z_3)|+|(z_3+2z_6+z_9)-(z_1+2z_4+z_7)|$$
对应的模板为图$d,e$，称为Soble算子