数字图像处理

  本篇文章主要介绍一些图像处理的基础，包括简单的灰度变换、直方图、卷积的简单介绍以及常见空间滤波器的介绍。   本篇内容主要知识来自冈萨雷斯第三版的《数字图像处理》，如对细节内容感兴趣的读者可以查阅相关书籍进行参考。本篇博文的内容与原书内容的顺序有所出入，如有不足欢迎留言指正，谢谢。

  众所周知，一幅常见的彩色图像通常是由M×N×3的矩阵表示的，其中M是图像的长度，N是图像的宽度，而3则代表三个颜色通道R,G,B(稍微提一下，opencv里读取图像默认的按B,G,R的顺序排序的)。通常，对图像的大部分操作并不需要用到彩色，因此将彩色图像转换为灰度图像是非常常见的。最常用的转换公式是： G r a y = 0.299 R + 0.587 G + 0.114 B Gray = 0.299R + 0.587G + 0.114B Gray=0.299R+0.587G+0.114B。   在图像转换为灰度图像后，就可以开始最基本的图像变换操作。以下子节都是建立在灰度图像的基础上进行操作的。其中， s s s 用来表示变换后的灰度值， r r r 用来表示变换前的灰度值， [ 0 , L − 1 ] [0,L-1] [0,L−1]表示灰度值的取值范围，通常L=256。   另外，需要说明的是，通常在灰度变换完毕后，需要对图像进行标准化，因为灰度变换的过程中，通常会出现小数或超过灰度范围255的数值，而最终图像需要正常显示，就需要将变换后的所有数据标准化到[0,255]的范围内。还需要注意的是，在编程过程中(尤其是python)，要注意图像变量的数据类型，如果数据类型是整数，很可能在图像灰度值计算的过程中对浮点数进行了截断，使得计算结果很不准确。所以，在进行计算之前，要将图像变量从uint8型转换为float32或其他浮点型。

s = L − 1 − r s=L-1-r s=L−1−r   简单来说就是黑白反转，有时候黑底作为背景图片更清晰，有时候白底更清晰，根据情况可进行选择。

s = c log ⁡ ( 1 + r ) s=c\log(1+r) s=clog(1+r)   其中，c 是可人为设置的常数。由于对数函数的性质，图像进行对数变换后，会将有较大动态范围(灰度级范围)的图像，压缩成动态范围较小的图像，且范围都比较集中于高亮度区域。

s = c r γ s=cr^{\gamma} s=crγ   其中，c 和 γ \gamma γ 是正常数。当 γ < 1 \gamma1 γ>1时，图像的直观效果是由亮变暗，和对数变换相反，此时图像灰度的动态范围被扩展；

  其实可以根据实际情况，设置自己的灰度变换函数，对图像灰度值进行分段函数处理。常见的方式有：   1).对比度拉伸   2).灰度级分层   用来突出某段灰度范围的图像，而其他范围的灰度保持不变。极端一点的情况如上图左，就是阈值化处理，将不需要的部分全部删去，将留下的部分全部变白，从而产生二值图的效果。

  3).比特平面分层   因为常见的灰度级都是 [ 0 , 2 8 − 1 ] [0,2^8-1] [0,28−1]，由8个字节(bit)组成。那么我们将每个灰度值都划分为： G r a y = g 1 × 2 0 + g 2 × 2 1 + ⋯ + g 8 × 2 7 Gray=g_1×2^0+g_2×2^1+\dots+g_8×2^7 Gray=g1​×20+g2​×21+⋯+g8​×27。然后利用所有像素的 g 1 g_1 g1​组成一幅图像， g 2 g_2 g2​组成一幅图像，以此类推。最终可以将原图像划分为8张图像，这8张图像就称为比特平面图像。以下是比特平面分层的效果图(图片源自《数字图像处理》，但是是从https://blog.csdn.net/be_solider/article/details/79114263博客中提供，其中博主将原图像灰度范围设定在[0,127]，所以第8比特平面是全黑的)

  直方图其实就是以灰度等级为横轴，图像中该等级灰度值的个数为纵轴的图像。熟悉概率与统计的读者可以这样理解，是以像素灰度值为随机变量而绘制的离散概率密度函数(PDF)图像。   熟悉摄影的读者可能知道，图像如果要有较好的视觉感受，通常要要较高的对比度。而高对比度在直方图中，通常表现为图像的灰度值均匀地分布在每个灰度级上，即图像直方图趋近于水平线。用PS等软件可以很轻易地做到提高图像对比度，实际上提高图像对比度的操作在图像处理中叫做直方图均衡。

  要完成直方图均衡的操作，只需要进行如下灰度变换即可： s = T ( r ) = ( L − 1 ) ∫ 0 r p r ( w ) d w 离 散 后 得 s k = T ( r k ) = ( L − 1 ) ∑ j = 0 k p r ( r j ) s=T(r)=(L-1)\int_0^r p_r(w)dw\quad离散后得s_k=T(r_k)=(L-1)\sum_{j=0}^k p_r(r_j) s=T(r)=(L−1)∫0r​pr​(w)dw离散后得sk​=T(rk​)=(L−1)j=0∑k​pr​(rj​)   上式最后是将连续函数离散化后的结果。其中， p r ( r j ) = n j M N p_r(r_j)=\frac{n_j}{MN} pr​(rj​)=MNnj​​，表示某个灰度级 r j r_j rj​在图像中出现的频率，而 n j n_j nj​就是其频数。 推导过程？   由于直方图是PDF的近似，实际中直方图是离散的且不允许产生新的灰度级，因此，在计算中对得到的结果通常都是取就近的灰度值。因此，尽管进行了直方图均衡，最后得到的直方图也只是较为平缓，而不会是完全平坦的。

  直方图匹配的意思是，有时候并不希望简单地进行直方图均衡，而是希望图像的直方图尽量地靠近另一张图像的直方图或者是人为设计的直方图。这一过程就叫做直方图匹配。   假设，原图像的PDF是 p r ( r ) p_r(r) pr​(r)，想要匹配的直方图的PDF是 p z ( z ) p_z(z) pz​(z)。那么，想要找到两者的变换，其实可以直接借用直方图均衡的操作作为桥梁，让原图像和匹配模板的图像都作直方图均衡化，理论上两者均衡化后的PDF应该是相同的，利用这一点建立等式关系，从而建立两者的转换关系。 s r = T r ( r ) = ( L − 1 ) ∫ 0 r p r ( w ) d w s_r=T_r(r)=(L-1)\int_0^r p_r(w)dw sr​=Tr​(r)=(L−1)∫0r​pr​(w)dw s z = T z ( z ) = ( L − 1 ) ∫ 0 z p z ( w ) d w s_z=T_z(z)=(L-1)\int_0^z p_z(w)dw sz​=Tz​(z)=(L−1)∫0z​pz​(w)dw ∵ s z = s r ∴ z = T z − 1 [ T r ( r ) ] = T z − 1 ( s r ) ∵s_z=s_r\quad∴z=T_z^{-1}[T_r(r)]=T_z^{-1}(s_r) ∵sz​=sr​∴z=Tz−1​[Tr​(r)]=Tz−1​(sr​)   由于图像是离散的而不是连续的，所以其实上述的函数逆变换可以不用求解，只需要找出两组映射关系得到的相同值进行匹配就可以了。综上，直方图匹配的步骤就是：

  另外，需要说明的是，要用直方图匹配的方式来进行图像增强，通常直方图模板的选择是个试凑的过程，目前并没有很好的理论支持。

  有时候，直方图并不一定需要对全局图像进行处理，可能只有图中一块子区域需要进行直方图处理，或者不同区域需要使用不同的直方图处理方式。这时候就需要对图像区别对待。另外说一句，开始介绍局部直方图处理时，就已经渐渐地有卷积的意识了。   想要进行局部直方图处理，通常会选择一个远远小于原图像尺寸的模板作为当前处理图像区域的邻域(掩模)。这个模板覆盖的图像区域，就是进行局部直方图处理的有效区域，其他区域都暂时不管。我们可以用局部邻域内灰度值的均值来考察该区域的亮度情况，用方差来考察该区域的对比度情况等等。通过一定的统计量分析后，再选择对该区域进行适当的直方图处理。当处理完该区域后，将模板移动至下一区域，通常两个区域需要有重叠部分，否则最终会使图像产生“棋盘效应”。

  放在最前面的话：空间滤波相比于频域滤波，最大的特点是可以用于非线性滤波，而频域滤波是做不到的。

  注意，可能很多人都听过卷积，但是却没听过相关，其实相关和卷积是两种非常类似的操作。其实两者的操作过程都相同，只是在开始操作前，卷积需要将滤波器旋转180°，再进行相关操作。为什么要反转180°？下述例子给了很好的解释：   假设有个一维的离散单位冲激序列[0 0 0 1 0 0 0]，滤波器序列为[1 2 3 4 5]，根据空间滤波的机理，进行一维方向的滤波过程如下表：

  上表中，每一次，原始数列都和卷积核对应的位置相乘，得到的结果放在卷积核中心对应的图像位置上进行覆盖。如果卷积核对应的原始序列没有数字，那么用0填充空白区域。填充这一步很重要，在以后的卷积操作中经常会用到。计算完一次后，滤波器向后移动一位，直到滤波器中心到达原始序列的末端为止。   通过上述过程，可以发现，相关操作最终得到的结果是将滤波器系数反转180°后放在原始序列中有单位冲激的地方；而卷积操作由于一开始就把滤波器旋转180°，所以最终得到的结果是正序的滤波器系数。这就是相关和卷积最大的区别，也是为什么卷积要先把滤波器旋转180°的原因。   另外，读者可以尝试以下在原始数列不同位置多放几个1，再重复一下卷积过程，就会发现，在每个1的位置都会复制一遍滤波器的系数，然后在序列的重叠区域会把每个滤波器系数进行叠加。这个过程有点像是多个函数的叠加过程，读者可以自己体会，有助于进一步理解卷积。   上述过程如果扩展到二维空间也是一样的，只是滤波器不只是向右移动，移动到最右端后会返回最左端并且向下移动，循环往复直到滤波器中心到达二维空间的右下角点。而如果学习了傅里叶变换后可以更深刻地理解卷积操作，也能知道，一个函数与单位冲激的卷积，就是在单位冲激位置处复制该函数。   通常，学术界并不严格区分相关和卷积的术语使用，所以阅读文献的适合需要小心判断。   综上所述，可以把每次卷积的操作过程用以下数学方式表达： w ( x , y ) ⋆ f ( x , y ) = g ( x , y ) = ∑ s = − a a ∑ t = − b b w ( s , t ) f ( x − s , y − t ) w(x,y)\star f(x,y)=g(x,y)=\sum_{s=-a}^a\sum_{t=-b}^b w(s,t)f(x-s,y-t) w(x,y)⋆f(x,y)=g(x,y)=s=−a∑a​t=−b∑b​w(s,t)f(x−s,y−t)

  线性滤波的向量表示   如果在二维卷积操作中，将滤波器系数和图像对应区域的灰度值都拉平，成为一维向量后再进行相乘求和，就是将卷积操作用向量来表示，卷积公式就变为： g ( x , y ) = ∑ k = 1 m n w k z k = w T z g(x,y)=\sum_{k=1}^{mn}w_kz_k=w^Tz g(x,y)=k=1∑mn​wk​zk​=wTz   其中，滤波器的尺寸为m×n，z是滤波器对应覆盖的图像区域的像素值。

  通过上文的介绍后，可以知道，使用空间滤波器对图像的卷积操作，就可以逐一对图像的某个子区域进行适当处理。如果设计好滤波器内部的系数，就会产生各种各样的效果。比如，如果滤波器是3×3的，所有系数设置为 1 9 \frac{1}{9} 91​，那么空间滤波器卷积后，就会对图像的每个3×3子区域进行求均值的操作。可想而知，在整个卷积完成后，图像应当比处理之前变得更加平滑，因为每个像素都被以它为中心的3×3子区域中所有像素的灰度均值替代。

  1).均值滤波器   不带权重的均值滤波器如下： 1 9 × [ 1 1 1 1 1 1 1 1 1 ] \frac{1}{9}×\begin{bmatrix}1&1&1\\1&1&1\\1&1&1\end{bmatrix} 91​×⎣⎡​111​111​111​⎦⎤​   带权重的均值滤波器如下，权重可根据情况自己设定： 1 16 × [ 1 2 1 2 4 2 1 2 1 ] \frac{1}{16}×\begin{bmatrix}1&2&1\\2&4&2\\1&2&1\end{bmatrix} 161​×⎣⎡​121​242​121​⎦⎤​   通常，带权重的均值滤波器可以用如下数学式表达： g ( x , y ) = ∑ s = − a a ∑ t = − b b w ( s , t ) f ( x + s , y + t ) ∑ s = − a a ∑ t = − b b w ( s , t ) g(x,y)=\frac{\sum_{s=-a}^a\sum_{t=-b}^bw(s,t)f(x+s,y+t)}{\sum_{s=-a}^a\sum_{t=-b}^bw(s,t)} g(x,y)=∑s=−aa​∑t=−bb​w(s,t)∑s=−aa​∑t=−bb​w(s,t)f(x+s,y+t)​   其中, ∑ s = − a a ∑ t = − b b w ( s , t ) \sum_{s=-a}^a\sum_{t=-b}^bw(s,t) ∑s=−aa​∑t=−bb​w(s,t)其实是常数，只要计算一次即可。   空间均值处理主要的应用是，模糊图像，以便将小物体融入背景，而将感兴趣的大物体变得模糊而易于检测。

  2).统计排序(非线性)滤波器   常见的有中值滤波器，对椒盐噪声的处理效果非常好，系数不好描述，看名字就知道，中值滤波器就是将滤波器覆盖的图像邻域中，取该区域的中位数作为当前像素点的替代值。另外常见的还有最大值、最小值滤波器。详情可以参见“图像退化与复原”的相关章节。

  3).锐化空间滤波器   顾名思义，锐化空间滤波器的作用就是使图像锐化，也就是让图像的边缘更加显著，更加不平滑。实际上，锐化滤波器的作用是提取图像中的边缘信息(也就是灰度值突变的信息)，然后再把信息加到原图像上，就可以起到突出边缘的作用。   首先介绍拉普拉斯算子实现的锐化滤波器如下： [ 0 1 0 1 − 4 1 0 1 0 ] [ 1 1 1 1 − 8 1 1 1 1 ] [ 0 − 1 0 − 1 4 − 1 0 − 1 0 ] [ − 1 − 1 − 1 − 1 8 − 1 − 1 − 1 − 1 ] \begin{bmatrix}0&1&0\\1&-4&1\\0&1&0\end{bmatrix}\quad \begin{bmatrix}1&1&1\\1&-8&1\\1&1&1\end{bmatrix}\quad \begin{bmatrix}0&-1&0\\-1&4&-1\\0&-1&0\end{bmatrix}\quad \begin{bmatrix}-1&-1&-1\\-1&8&-1\\-1&-1&-1\end{bmatrix} ⎣⎡​010​1−41​010​⎦⎤​⎣⎡​111​1−81​111​⎦⎤​⎣⎡​0−10​−14−1​0−10​⎦⎤​⎣⎡​−1−1−1​−18−1​−1−1−1​⎦⎤​   其中，后两个其实就是前两个乘以-1之后的滤波器。乘以-1之后更常用，因为按上述说法，锐化滤波器只是将边缘信息提取出来，还要再加回原图像，因此有如下公式： g ( x , y ) = f ( x , y ) + c [ ▽ 2 f ( x , y ) ] g(x,y)=f(x,y)+c[\bigtriangledown^2f(x,y)] g(x,y)=f(x,y)+c[▽2f(x,y)]   当取前两个滤波器模板时，常数c应当取-1；而当取后两个模板时，常数c直接取1即可。   那么，上述锐化滤波器是如何构造的呢？和拉普拉斯算子又有什么关系？   首先我们都知道，图像的边缘信息其实就是灰度值突变的信息。能度量这种数值突变的程度，通常使用的就是一阶微分和二阶微分。而二阶微分在增强细节方面要比一阶微分好得多，所以选择图像的二阶微分作为边缘信息的代表。   其次，图像数据是离散数据，可以知道图像的一阶微分是如下表达式： ∂ f ∂ x = f ( x + 1 ) − f ( x ) x + 1 − x = f ( x + 1 ) − f ( x ) \frac{\partial f}{\partial x}=\frac{f(x+1)-f(x)}{x+1-x}=f(x+1)-f(x) ∂x∂f​=x+1−xf(x+1)−f(x)​=f(x+1)−f(x)   那么，其二阶微分就是： ∂ 2 f ∂ x 2 = ∂ f ∂ x ( x ) − ∂ f ∂ x ( x − 1 ) = f ( x + 1 ) − f ( x ) − [ f ( x ) − f ( x − 1 ) ] = f ( x + 1 ) + f ( x − 1 ) − 2 f ( x ) \frac{\partial^2 f}{\partial x^2}=\frac{\partial f}{\partial x}(x)-\frac{\partial f}{\partial x}(x-1)=f(x+1)-f(x)-[f(x)-f(x-1)]=f(x+1)+f(x-1)-2f(x) ∂x2∂2f​=∂x∂f​(x)−∂x∂f​(x−1)=f(x+1)−f(x)−[f(x)−f(x−1)]=f(x+1)+f(x−1)−2f(x)   同理，图像是二维函数，所以有： ∂ 2 f ∂ x 2 = f ( x + 1 , y ) + f ( x − 1 , y ) − 2 f ( x , y ) \frac{\partial^2 f}{\partial x^2}=f(x+1,y)+f(x-1,y)-2f(x,y) ∂x2∂2f​=f(x+1,y)+f(x−1,y)−2f(x,y) ∂ 2 f ∂ y 2 = f ( x , y + 1 ) + f ( x , y − 1 ) − 2 f ( x , y ) \frac{\partial^2 f}{\partial y^2}=f(x,y+1)+f(x,y-1)-2f(x,y) ∂y2∂2f​=f(x,y+1)+f(x,y−1)−2f(x,y)   综上，图像的二维梯度，也就是拉普拉斯算子为： ▽ 2 f = ∂ 2 f ∂ x 2 + ∂ 2 f ∂ y 2 = f ( x + 1 , y ) + f ( x − 1 , y ) + f ( x , y + 1 ) + f ( x , y − 1 ) − 4 f ( x , y ) \bigtriangledown^2f=\frac{\partial^2 f}{\partial x^2}+\frac{\partial^2 f}{\partial y^2}=f(x+1,y)+f(x-1,y)+f(x,y+1)+f(x,y-1)-4f(x,y) ▽2f=∂x2∂2f​+∂y2∂2f​=f(x+1,y)+f(x−1,y)+f(x,y+1)+f(x,y−1)−4f(x,y)   是不是很眼熟，没错，如果以像素点(x,y)为中心，拉普拉斯算子的结果其实就已经说明了滤波器的系数，也就是 [ 0 1 0 1 − 4 1 0 1 0 ] \begin{bmatrix}0&1&0\\1&-4&1\\0&1&0\end{bmatrix} ⎣⎡​010​1−41​010​⎦⎤​。   第二个滤波器系数是怎么得到的呢？刚刚求出的滤波器只对x,y两个方向求了二阶微分，如果对两个斜对角方向也求二阶微分，结果和上述情况差不多，只是滤波器的系数就变成了 [ 1 1 1 1 − 8 1 1 1 1 ] \begin{bmatrix}1&1&1\\1&-8&1\\1&1&1\end{bmatrix} ⎣⎡​111​1−81​111​⎦⎤​。

  4).非锐化掩蔽和高提升滤波   其实要得到图像的边缘信息，还有一个很好的思路，那就是先用平滑滤波器让图像变平滑(模糊)，然后再用原图像减去模糊后的图像，那么边缘信息由于被模糊了和原始图像之间肯定存在差值，这样就能提取边缘信息。再将边缘信息叠加到原图像，增强图像边缘信息。   上述做法就叫做非锐化掩蔽，以 f ˉ ( x , y ) \bar f(x,y) fˉ​(x,y)表示平滑处理后的图像，则非锐化掩蔽得到的边缘模板为： g m a s k ( x , y ) = f ( x , y ) − f ˉ ( x , y ) g_{mask}(x,y)=f(x,y)-\bar f(x,y) gmask​(x,y)=f(x,y)−fˉ​(x,y)   将提取到的 g m a s k g_{mask} gmask​叠加到原图像中，就可以增强图像的边缘信息，得到的增强后图像为： g ( x , y ) = f ( x , y ) + k × g m a s k ( x , y ) g(x,y)=f(x,y)+k×g_{mask}(x,y) g(x,y)=f(x,y)+k×gmask​(x,y)   其中，k为权重系数，通常 k ≥ 0 k\ge0 k≥0。当k=1时，这种操作就叫做非锐化掩蔽。当k>1时，就叫做高提升滤波。注意，k太大会产生不好的效果，取值时需小心。

  5).使用一阶微分锐化图像——Sobel算子   当然，也可以用一阶微分来锐化图像。图像处理中，一阶微分是用梯度幅值来实现的，即： m a g ( ▽ f ) = ( ∂ f / ∂ x ) 2 + ( ∂ f / ∂ y ) 2 mag(\bigtriangledown f)=\sqrt{(\partial f/\partial x)^2+(\partial f/\partial y)^2} mag(▽f)=(∂f/∂x)2+(∂f/∂y)2 ​   在使用某些机器计算上述梯度幅值时，上述表达式不太方便，通常用两个方向上的梯度绝对值来近似： m a g ( ▽ f ) ≈ ∣ g x ∣ + ∣ g y ∣ mag(\bigtriangledown f)≈|g_x|+|g_y| mag(▽f)≈∣gx​∣+∣gy​∣   最终研究得到的一阶微分锐化滤波器模板称为"Sobel算子"，其有： g x = ( z 7 + 2 z 8 + z 9 ) − ( z 1 + 2 z 2 + z 3 ) g_x=(z_7+2z_8+z_9)-(z_1+2z_2+z_3) gx​=(z7​+2z8​+z9​)−(z1​+2z2​+z3​) g y = ( z 3 + 2 z 6 + z 9 ) − ( z 1 + 2 z 4 + z 7 ) g_y=(z_3+2z_6+z_9)-(z_1+2z_4+z_7) gy​=(z3​+2z6​+z9​)−(z1​+2z4​+z7​)   其中， z i , i = 1 , 2 , . . . , 9 z_i,i=1,2,...,9 zi​,i=1,2,...,9表示的是3×3模板的对应图像像素值，即 [ z 1 z 2 z 3 z 4 z 5 z 6 z 7 z 8 z 9 ] \begin{bmatrix}z_1&z_2&z_3\\z_4&z_5&z_6\\z_7&z_8&z_9\end{bmatrix} ⎣⎡​z1​z4​z7​​z2​z5​z8​​z3​z6​z9​​⎦⎤​。   所以，Sobel算子有： m a g ( ▽ f ) ≈ ∣ ( z 7 + 2 z 8 + z 9 ) − ( z 1 + 2 z 2 + z 3 ) ∣ + ∣ ( z 3 + 2 z 6 + z 9 ) − ( z 1 + 2 z 4 + z 7 ) ∣ mag(\bigtriangledown f)≈|(z_7+2z_8+z_9)-(z_1+2z_2+z_3)|+|(z_3+2z_6+z_9)-(z_1+2z_4+z_7)| mag(▽f)≈∣(z7​+2z8​+z9​)−(z1​+2z2​+z3​)∣+∣(z3​+2z6​+z9​)−(z1​+2z4​+z7​)∣ 等 价 于 模 板 [ − 1 − 2 − 1 0 0 0 1 2 1 ] 和 [ − 1 0 1 − 2 0 2 − 1 0 1 ] 等价于模板\begin{bmatrix}-1&-2&-1\\0&0&0\\1&2&1\end{bmatrix}和\begin{bmatrix}-1&0&1\\-2&0&2\\-1&0&1\end{bmatrix} 等价于模板⎣⎡​−101​−202​−101​⎦⎤​和⎣⎡​−1−2−1​000​121​⎦⎤​

  通常，只使用一种方法进行图像增强是不够的，我们需要根据经验，结合图像灰度变换、直方图处理、空间滤波等操作，使用多种技巧对图像进行混合处理，才能达到想要的目标。这时候实践经验就很重要了，可以多学习学习其他人图像处理的方式积累经验，了解每种图像处理方法的使用场合和效果，然后及时记录经验总结。毕竟，实践才是检验真理的唯一标准。