图像的平滑滤波、边缘检测、阈值化、腐蚀膨胀等原理介绍

一幅图像可以定义为一个二维函数f(x,y)f(x,y)f(x,y)，其中xxx和yyy是空间/平面坐标，而在任何一对空间坐标(x,y)(x,y)(x,y)处的幅值fff称为图像在该点处的强度或灰度。彩色图像是由三个二维灰度图像 f(x,y)f(x,y)f(x,y)组成，三个通道分别是RGB或HSV。一般来讲，对图像进行处理的主要目的有三个方面：

改善图像的质量。比如，进行图像的亮度、彩色变换，增强、抑制某些成分，对图像进行几何变换等。
提取图像中所包含的某些特征或特殊信息。提取的特征可以包括很多方面，如频域特征、灰度或颜色特征、边界特征、区域特征、纹理特征、形状特征、拓扑特征和关系结构等，这些被提取的特征往往为计算机分析图像提供便利。
便于图像的存储和传输。比如图像数据的变换、编码和压缩，为了解决编码冗余、空间冗余和时间冗余。

基础概念

在图像空间域上的操作是直接对于图像像素的运算，通用表达式为：g(x,y)=T[f(x,y)]g(x,y)=T[f(x,y)]g(x,y)=T[f(x,y)]，g(x,y)g(x,y)g(x,y)为输入图像经过处理得到的输出图像，T[⋅]T[\cdot]T[⋅]是在点(x,y)(x,y)(x,y)的邻域上定义的一种算子。像素邻域有三种类型：四邻域，位于坐标(x,y)(x,y)(x,y)的像素点PPP有水平和垂直的四个相邻像素，表示为N4(P)N_4(P)N4(P)；对角邻域，位于坐标(x,y)(x,y)(x,y)的像素点PPP在对角上的四个相邻像素，表示为ND(P)N_D(P)ND(P)；八邻域，为N4(P)N_4(P)N4(P)与ND(P)N_D(P)ND(P)点的集合，表示为N8(P)N_8(P)N8(P)。如下图所示：

三种距离度量如下：

对像素的基本操作如下图所示，邻域原点依次从一个像素向另一个像素移动，对邻域内的像素应用算子 TTT，并在该位置产生输出。一般从输入图像的左上角开始，以水平扫描的方式逐像素处理，每次一行。当该邻域原点位于图像边界上时，部分邻域将位于图像的外部，此时可以在用 TTT 做指定计算时忽略外侧邻点，或者用0/其他指定的灰度值填充图像边缘。

考虑一种特殊的情况，即最小邻域1×1时，g(x,y)g(x,y)g(x,y)仅取决于每一个(x,y)(x,y)(x,y)处的fff值，这时的 TTT 则成为一个灰度级的映射函数。若(x,y)(x,y)(x,y)处变换前的灰度为rrr，变换后的灰度为s，则根据不同的变换函数，如线性函数、对数函数、幂律函数、分段线性函数等，可以完成图像灰度级的扩展/压缩，进行对比度拉伸、阈值化处理等图像增强技术。

平滑空间滤波

滤波主要通过一个空间滤波器（或称空间掩模、核、模板和窗口）直接作用于图像本身实现的，空间滤波器由一个邻域及对该邻域包围的图像像素执行的预定义操作组成。滤波器的中心访问输入图像的每个像素，产生的新像素坐标等于邻域中心的坐标，新像素的值是滤波操作的结果。平滑滤波器用于模糊处理和降低噪声。模糊处理经常用于预处理任务中，例如在大目标提取之前去除图像中的一些琐碎细节，以及桥接直线或曲线的缝隙。通过线性滤波和非线性滤波模糊处理，可以降低噪声。

均值滤波：将中心像素的邻域内的每个像素值用平均像素值替换，这种处理结果降低了图像灰度的尖锐变化，由于典型的随机噪声由灰度级的急剧变化组成，因此可以降低噪声。但是，对噪声图有特别大的孤立点时非常敏感，即使有极少数点存在较大差异也会导致平均值的明显波动，而下面的中值滤波可以避免这点影响。如下图所示，每个模板前面的乘积为1除以所有系数和。
中值滤波：将中心像素的邻域内的每个像素值用中间像素值替换，对处理椒盐噪声非常有效。然而除了噪声，图像边缘也是由图像灰度尖锐变化带来的特性，所以平滑滤波也会存在边缘模糊的负面效应，但是同样尺寸的中值滤波器比均值滤波器的模糊程度低。滤波器随着核大小逐渐变大，会让图像变得更加模糊。3×3的中值滤波器可以表示为：g(x,y)=Median(f1,f2,f3,f4,f5,f6,f7,f8,f9)g(x,y) = Median(f_1,f_2,f_3,f_4,f_5,f_6,f_7,f_8,f_9)g(x,y)=Median(f1,f2,f3,f4,f5,f6,f7,f8,f9)。
高斯滤波：图像高斯滤波也是邻域平均的思想对图像进行平滑的一种方法，不过对图像进行平均时，不同位置的像素被赋予了不同的权重。使用高斯滤波的原因是，真实图像在空间内的像素是缓慢变化的，因此邻近点的像素变化不会很明显。所以高斯滤波让临近的像素具有更高的重要度，对周围像素计算加权平均值，较近的像素具有较大的权重值。也正是因为这个原因，高斯滤波会磨平边缘，在接近边缘处的像素不应该与相邻像素相关。g(x,y)=12πσ2e−x2+y22σ2g(x,y) = \frac{1}{2\pi\sigma^2}e^{-\frac{x^2+y^2}{2\sigma^2}}g(x,y)=2πσ21e−2σ2x2+y2，3×3的高斯模板还要确保9个点的和为1。

边缘检测

边缘检测也可以说是图像的锐化处理，主要目的是突出灰度的过渡部分，与平滑滤波的处理相反，并且图像平滑与积分的过程类似，因此锐化处理可由空间微分实现。常用的边缘检测算子可以分为：一阶微分算子Roberts、Sobel、Prewitt，二阶微分算子Laplacian，非微分边缘检测算子Canny。对于一维函数f(x)f(x)f(x)的一阶微分的定义为差值，即∂f(x)∂x=f(x+1)−f(x)\frac{\partial f(x)}{\partial x}=f(x+1)-f(x)∂x∂f(x)=f(x+1)−f(x)二阶微分定义如下差分，即∂2f(x)∂x2=f(x+1)+f(x−1)−2f(x)\frac{\partial^2 f(x)}{\partial x^2}=f(x+1)+f(x-1)-2f(x)∂x2∂2f(x)=f(x+1)+f(x−1)−2f(x)

Roberts算子：与标准一阶差分不同，Robert采用对角线差分，如下图所示。红色为计算的参考点，也就是计算完成后的值会放在它的位置处。
Sobel算子：Sobel算子考虑了水平、垂直和2个对角共计4个方向对的梯度加权求和，是一个3×3各向异性的梯度算子。定义一个给定邻域方向梯度矢量g的幅度=像素灰度差分/相邻像素的距离，Sobel采用的像素距离是一种城市距离，矢量g的方向可以通过中心像素z5相关邻域的单位矢量给出，这里的邻域是对称出现的，即四个方向对 (z1, z9) ，(z2, z8)， (z3, z7)， (z6, z4)。沿着4个方向对上求其梯度矢量和，可以给出当前像素(z5)的平均梯度估计：g=(z3−z7)/4∗[1,1]+(z1−z9)/4∗[−1,1]+(z2−z8)/2∗[0,1]+(z6−z4)/2∗[1,0]=[(z3−z7−z1+z9)/4+(z6−z4)/2,(z3−z7+z1−z9)/4+(z2−z8)/2]g = (z3-z7)/4 * [1, 1]+ (z1-z9)/4 * [-1, 1]+ (z2-z8)/2 * [0, 1]+ (z6-z4)/2 * [1, 0]\\ = [(z3- z7- z1+ z9)/4 + (z6-z4)/2, (z3- z7+ z1- z9)/4 + (z2-z8)/2] g=(z3−z7)/4∗[1,1]+(z1−z9)/4∗[−1,1]+(z2−z8)/2∗[0,1]+(z6−z4)/2∗[1,0]=[(z3−z7−z1+z9)/4+(z6−z4)/2,(z3−z7+z1−z9)/4+(z2−z8)/2]
其中，4个单位向量 [1, 1]，[-1, 1]，[0, 1], [1, 0] 控制差分的方向，系数1/4, 1/2为距离反比权重。
为了方便起见，将上式乘以4，得到G=[z3+2∗z6+z9−z1−2∗z4−z7,z1+2∗z2+z3−z7−2∗z8−z9]G = [z3 + 2*z6 + z9 - z1 - 2*z4-z7, z1 + 2*z2 + z3 - z7 - 2*z8 - z9]G=[z3+2∗z6+z9−z1−2∗z4−z7,z1+2∗z2+z3−z7−2∗z8−z9]
按x-y方向，可分别写成：Gx=(z3+2∗z6+z9)−(z1+2∗z4+z7)Gy=(z1+2∗z2+z3)−(z7+2∗z8+z9)Gx = (z3 + 2*z6 + z9) - (z1 + 2*z4 + z7)\\ Gy = (z1 + 2*z2 + z3) - (z7 + 2*z8 + z9)Gx=(z3+2∗z6+z9)−(z1+2∗z4+z7)Gy=(z1+2∗z2+z3)−(z7+2∗z8+z9)
如下图所示：

最后，G(x,y)=Gx2+Gy2G(x, y) = \sqrt{Gx^2 + Gy^2}G(x,y)=Gx2+Gy2就是 sobel 边缘检测后的图像了。
Prewitt算子：利用像素点上下、左右邻点的灰度差，在边缘处达到极值检测边缘，在图像空间利用两个方向模板与图像进行邻域卷积来完成，这两个方向模板一个检测水平边缘，一个检测垂直边缘。定义如下图：
Laplacian算子：拉普拉斯算子是最简单的各项同性二阶微分算子，一个二维图像f(x,y)f(x,y)f(x,y)的拉普拉斯算子定义为：∇2f=∂2f(x,y)∂x2+∂2f(x,y)∂y2\nabla^2f=\frac{\partial^2 f(x,y)}{\partial x^2}+\frac{\partial^2 f(x,y)}{\partial y^2}∇2f=∂x2∂2f(x,y)+∂y2∂2f(x,y)
那么在两个方向上有：∂2f(x,y)∂x2=f(x+1,y)+f(x−1,y)−2f(x,y)∂2f(x,y)∂y2=f(x,y+1)+f(x,y−1)−2f(x,y)\frac{\partial^2 f(x,y)}{\partial x^2}=f(x+1,y)+f(x-1,y)-2f(x,y)\\ \frac{\partial^2 f(x,y)}{\partial y^2}=f(x,y+1)+f(x,y-1)-2f(x,y)∂x2∂2f(x,y)=f(x+1,y)+f(x−1,y)−2f(x,y)∂y2∂2f(x,y)=f(x,y+1)+f(x,y−1)−2f(x,y)
则，∇2f(x,y)=f(x+1,y)+f(x−1,y)+f(x,y+1)+f(x,y−1)−4f(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)
由于拉普拉斯是二阶微分算子，因此强调的是图像中灰度的突变，并不强调灰度级缓慢变化的区域，简单地将原图像与拉普拉斯变换后的图像叠加到一起，可以复原背景特性并保持拉普拉斯锐化处理的效果，表示为g(x,y)=f(x,y)+c[∇2f(x,y)]g(x,y)=f(x,y)+c[\nabla^2f(x,y)]g(x,y)=f(x,y)+c[∇2f(x,y)]。如上图所示的拉普拉斯算子c=-1。
Canny算子：从表面效果上来讲，Canny算法是对Sobel、Prewitt等算子效果的进一步细化和更加准确的定位。检测步骤如下：
step1：用高斯滤波器平滑图像；
step2：用一阶偏导的有限差分来计算梯度的幅值和方向（可选用的模板Soble算子、Prewitt算子、Roberts模板等）；
step3：对梯度幅值进行非极大值抑制；
在每一点上，邻域中心(x，y)(x，y)(x，y)与沿着其对应的梯度方向的两个像素相比，若中心像素为最大值，则保留；否则中心置0。这样可以抑制非极大值，保留局部梯度最大的点，以得到细化的边缘。
step4：用双阈值算法检测和连接边缘。
选取系数th_high和th_low，取值比率为2:1或3:1，将小于低阈值的点抛弃，赋值为0；将大于高阈值的点标记，确定其为边缘点，赋值为1或255；将在两个阈值之间的点使用8连通区域确定，只有与th_high像素邻接时才会被接受，成为边缘点，赋值为1或255。

阈值化

图像阈值化是一种广泛应用的分割技术，利用图像中要提取的目标区域与其背景在灰度特性上的差异，把图像看作具有不同灰度级的两类区域(目标区域和背景区域)的组合，选取一个比较合理的阈值，以确定图像中每个像素点应该属于目标区域还是背景区域，从而产生相应的二值图像。
1、全局阈值处理：
当物体和背景像素的灰度分布十分明显时，可以用适用于整个图像的单个/全局阈值。如下图所示的二值阈值化，重要的是一幅图像如何自动估计阈值，有如下迭代算法：
step1：为全局阈值T选择一个初始估计值；
step2：用图中的公式将该图像产生两组像素：G1由灰度值大于T的所有像素组成，G2由所有小于等于T的像素组成；
step3：对G1和G2的像素分别计算平均灰度值m1和m2；
step4：计算一个新的阈值：T=12(m1+m2)T=\frac12(m_1+m_2)T=21(m1+m2)
step5：重复步骤2-4，直到T值间的差小于一个预定义的参数ΔT\Delta TΔT为止。

其他阈值化方法：

2、局部阈值处理：
根据像素的邻域块的像素值分布来确定该像素位置上的二值化阈值，又叫自适应阈值。这样做的好处在于每个像素位置处的二值化阈值不是固定不变的，而是由其周围邻域像素的分布来决定的。亮度较高的图像区域的二值化阈值通常会较高，而亮度较低的图像区域的二值化阈值则会相适应地变小。不同亮度、对比度、纹理的局部图像区域将会拥有相对应的局部二值化阈值。常用的局部自适应阈值有：1）局部邻域块的均值；2）局部邻域块的高斯加权和。

腐蚀膨胀

1、腐蚀和膨胀：
基于通过阈值化获得的二进制图像，将图像与核（模板）进行卷积，核可以是任意形状和大小，拥有一个单独定义的参考点。多数情况下，核是一个小的中间带有参考点的实心正方形或圆盘。腐蚀是求局部最小值的操作，膨胀是腐蚀的反操作，取得是局部最大值。如下图所示，腐蚀就是缩小了高亮区域，能够消除较细小的凸起；而膨胀会使图像中高亮区域逐渐增长，可以填补凹洞。核与图像的卷积，即计算核覆盖的区域的像素的最小值/最大值，并把这个值赋值给参考点指定的像素，图中的红色五角星为该模板的参考点。

2、开运算和闭运算：
开运算的情况下，需要先腐蚀后膨胀，一般用来统计二值图像中的区域数，可以将细微连在一起的两块目标分开，能够除去孤立的小点、毛刺和小桥，而总的位置和形状不变。闭运算的情况下，需要先膨胀后腐蚀，来消除噪声引起的区域，能够填平小孔、弥合小裂缝，而总的位置和形状不变。对于连通域分析，通常先采用腐蚀或闭运算来消除纯粹由噪声引起的部分，然后用开运算来连接邻近的区域。

拉伸、收缩、扭曲和旋转时的插值方法

在对图像进行空间变换的过程中，典型的情况是在对图像进行放大、旋转处理的时候，图像会出现失真的现象。这是由于在变换之后的图像中，存在着一些变换之前的图像中没有的像素位置。处理这一问题的方法被称为图像灰度级插值。常用的插值方式有三种：最近邻域插值、双线性插值、双三次插值。理论上来讲，最近邻域插值的效果最差，双三次插值的效果最好，双线性插值的效果介于两者之间。不过对于要求不是非常严格的图像插值而言，使用双线性插值通常就足够了。

最近邻插值：对于通过映射尺寸变换得到的一个浮点坐标，对其进行简单的取整（四舍五入），得到一个整数型坐标，这个整数型坐标对应的像素值就是目标像素的像素值。也就是说，取浮点坐标最邻近的左上角点对应的像素值，也就是在原图上找最邻近的像素进行赋值。最邻近插值简单且直观，但得到的图像质量不高，因为插值后图像灰度不连续，因此在灰度变化的地方可能出现明显的锯齿状，放大后的图像有很严重的马赛克，缩小后的图像有很严重的失真。
双线性插值：双线型插值算法就是一种比较好的图像缩放算法，它充分的利用了原图中虚拟点四周的四个真实存在的最近邻像素值来共同决定目标图中的一个像素值，因此缩放效果比简单的最邻近插值要好很多。双线性内插法的计算比最邻近点法复杂，计算量较大，但没有灰度不连续的缺点，结果基本令人满意。它具有低通滤波性质，使高频分量受损，图像轮廓可能会有一点模糊。计算方式如下，公式中对像素值的加权意义在于距离更近的权重更大：
三次线性插值：复杂度更高，具有16个最近邻点。首先对原图像附近的4×4个邻近像素进行三次样条拟合，然后将目标像素对应的三次样条值作为目标图像对应像素点的值。具体的计算过程可以参考https://zhuanlan.zhihu.com/p/62860859。