形态学图像处理
形态学主要应用是从图像中提取对于表达和描绘区域形状有意义的图像分量,为后续的目标识别等做准备;同时细化、像素化和修建毛刺等技术也常常应用于图像的预处理和后处理中。
基本思想是用具有一定形态的结构元素去度量和提取图像中的对应形状来对图像分析识别
预备知识
在数字图像处理的形态学运算中,常把一幅图像或者图像中一个我们感兴趣的区域称为集合,用大写字母A、B、C、D等表示;元素常指一个单个像素。像素在图像中的位置坐标用
z
=
(
z
1
,
z
2
)
z = (z_1,z_2)
z=(z1,z2)来表示。
z
∈
Z
2
z\in Z^2
z∈Z2,其中
Z
2
Z^2
Z2为二元整数序偶对的集合。
集合
B
B
B的反射用
B
^
\hat{B}
B^表示:
B
^
=
{
w
∣
w
=
−
b
,
b
∈
B
}
\hat{B} = \{w | w = -b,b\in B\}
B^={w∣w=−b,b∈B}
平移表示为
(
B
)
z
(B)_z
(B)z则有:
(
B
)
z
=
{
c
∣
c
=
b
+
z
,
b
∈
B
}
(B)_z = \{c | c = b+z,b\in B\}
(B)z={c∣c=b+z,b∈B}
假设B是坐标集合,那么有示例图如下:
反射和平移用来表达基于结构元(SE)的操作:研究一幅图像中感兴趣特性所用的小集合或子图像。结构元就类似于空间滤波时的卷积核(模板)。
SE中的值可以是0或1。常见的结构元结构有圆形、矩形、十字形。
当不在意SE中某些点时,该位置用‘x’标记。
结构元的原点:原点的选取与具体问题有关,当SE对称且未显示原点时,则假定原点位于对称中心处。
对图像操作时,要求结构元是矩形阵列。非矩形阵列的结构元,可通过添加最小适合数量的背景元素形成一个矩形阵列来实现。
二值形态学
腐蚀和膨胀
腐蚀
作为
Z
2
Z^2
Z2中的集合A和B,定义B对A的腐蚀如下:
A
⊖
B
=
{
z
∣
(
B
)
z
⊆
A
}
A\ominus B = \{z| (B)_z\subseteq A\}
A⊖B={z∣(B)z⊆A}
说人话就是将结构元B放到A中,所有的B的中心构成的集合就是腐蚀的结果,类似下图:
也就是说,图像本身、结构元的形状和原点位置的选取都会影响腐蚀操作的结果。
腐蚀能够消融物体的边界,而具体的腐蚀结果与图像本身和结构元的形状有关:
- 如果物体整体上大于结构元,腐蚀的结构是使物体变“瘦”一圈,这一圈到底有多大是由结构元决定的
- 如果物体本身小于结构元,则在腐蚀后的图像中物体将完全消失
- 如物体仅有部分区域小于结构元(如细小的连通),则腐蚀后物体会在细连通处断裂,分离为两部分
根据上述特性,在实际应用中,可以利用腐蚀运算去除物体之间的粘连,消除图像中的小颗粒噪声。
一个简单的例子:
import cv2
import numpy as np
def image_erode(image_path):
# 读取图像
image = cv2.imread(image_path)
# 创建腐蚀核(结构元素)
kernel_3 = np.ones((3, 3), np.uint8)
kernel_5 = np.ones((5, 5), np.uint8)
# 执行腐蚀操作
eroded_image_3 = cv2.erode(image, kernel_3, iterations=1)
eroded_image_5 = cv2.erode(image, kernel_5, iterations=1)
# 显示图像
cv2.imshow("Original", image)
cv2.imshow("Eroded 3x3", eroded_image_3)
cv2.imshow("Eroded 5x5", eroded_image_5)
cv2.waitKey(0)
cv2.destroyAllWindows()
# 图像路径
image_path = "images/lena_gray_512.tif"
# 调用函数进行腐蚀
image_erode(image_path)
结果如下:
可以明显看到5*5的图像更加模糊,这说明腐蚀使得原本的图像丢失了更多的细节,也就是说,更大的结构元容易腐蚀更多的细节。
膨胀
膨胀会“增化”或“粗化”图像中的物体。
类似于腐蚀,作为
Z
2
Z^2
Z2中的集合A和B,定义B对A的腐蚀如下:
A
⊕
B
=
{
z
∣
(
B
^
)
z
∩
A
≠
∅
}
A\oplus B = \{z| (\hat{B})_z\cap A \neq \varnothing \}
A⊕B={z∣(B^)z∩A=∅}
说人话就是,在
Z
2
Z^2
Z2上,只要结构元B有一部分跟A重合,那么B的中心构成的集合就是膨胀的结果。
可视化就是下图:
代码示例:
import cv2
import numpy as np
def image_dilate(image_path):
# 读取图像
image = cv2.imread(image_path)
# 创建膨胀核(结构元素)
kernel_100 = np.ones((100, 100), np.uint8)
kernel_300 = np.ones((300, 300), np.uint8)
# 执行膨胀操作
dilated_image_100 = cv2.dilate(image, kernel_100, iterations=1)
dilated_image_300 = cv2.dilate(image, kernel_300, iterations=1)
# 显示原始图像和膨胀后的图像
cv2.imshow("Original Image", image)
cv2.imshow("Dilated Image 10", dilated_image_100)
cv2.imshow("Dilated Image 30", dilated_image_300)
cv2.waitKey(0)
cv2.destroyAllWindows()
# 图像路径
image_path = "ch5/1.tif"
# 调用函数进行膨胀
image_dilate(image_path)
结果示例:
可以看到上图中,设置的膨胀核越大,原图像膨胀的越明显。这是由于原图像与膨胀核更容易有重叠部分导致的。
开操作与闭操作
- 开操作:先腐蚀后膨胀,开操作一般会平滑物体的轮廓、断开较窄的狭颈和消除细突出物(毛刺)。
- 闭操作:先膨胀后腐蚀,闭操作也会平滑轮廓的一部分,但通常会弥合较窄的间断和细长的沟壑,消除较小的孔洞,填补轮廓线中的断裂。
开操作
定义:
A
∘
B
=
(
A
⊖
B
)
⊕
B
A \circ B = (A \ominus B)\oplus B
A∘B=(A⊖B)⊕B
效果模拟:
示例代码:
import cv2
import numpy as np
def image_open(image_path):
# 读取图像
image = cv2.imread(image_path)
# 创建开操作的核(结构元素)
kernel = np.ones((3, 3), np.uint8)
# 执行开操作
opened_image = cv2.morphologyEx(image, cv2.MORPH_OPEN, kernel)
# 显示原始图像和开操作后的图像
cv2.imshow("Original Image", image)
cv2.imshow("Opened Image", opened_image)
cv2.waitKey(0)
cv2.destroyAllWindows()
# 图像路径
image_path = "ch9/1.tif"
# 调用函数进行开操作
image_open(image_path)
效果对比:
可以看到,指纹以外的很多噪点都被消除掉了,这是因为开操作是先腐蚀再膨胀。在腐蚀的过程中,这些噪点由于无法满足腐蚀的要求因而被消除掉了。同时消除后的图像由于又执行了膨胀操作,因此很多细节又被还原了。所以得到了如图所示的效果。
闭操作
定义
A
∙
B
=
(
A
⊕
B
)
⊖
B
A \bullet B = (A \oplus B)\ominus B
A∙B=(A⊕B)⊖B
效果模拟:
示例代码:
import cv2
import numpy as np
def image_open_and_close(image_path):
# 读取图像
image = cv2.imread(image_path)
# 创建开操作的核(结构元素)
kernel = np.ones((3, 3), np.uint8)
# 执行开操作
opened_image = cv2.morphologyEx(image, cv2.MORPH_OPEN, kernel)
# 执行闭操作
opened_closed_image = cv2.morphologyEx(opened_image, cv2.MORPH_CLOSE, kernel)
# 显示图像
cv2.imshow("Original Image", image)
cv2.imshow("Opened and closed Image", opened_closed_image)
cv2.waitKey(0)
cv2.destroyAllWindows()
# 图像路径
image_path = "ch9/1.tif"
# 调用函数进行开操作
image_open_and_close(image_path)
下图是先进行了开操作然后进行了闭操作的图像,可以看到完美消除噪点的同时,也对指纹进行了一定程度的恢复,但是仍旧有一些纹路断开了。这应该是因为开操作对于原图像的指纹断开的太过彻底所导致的。
击中或击不中变换
简单来说击中-击不中运算常用于二值图像,它用于基于结构元素的配置,从图像中寻找具有某种像素排列特征的目标,如单个像素、颗粒中交叉或纵向的特征、直角边缘或其他用户自定义的特征等。计算时,只有当结构元素与其覆盖的图像区域完全相同时,中心像素的值才会被置为1,否则为0。
例如:
用人话讲就是,只有结构与图中的像素排列完全一致时才将结构中心置为1
击中与击不中变换可表示为:
A
⊛
B
=
(
A
⊖
B
1
)
∩
(
A
c
⊖
B
2
)
A\circledast B = (A\ominus B_1)\cap (A^c \ominus B_2)
A⊛B=(A⊖B1)∩(Ac⊖B2)
其中,
B
1
∩
B
2
=
∅
B_1 \cap B_2 = \varnothing
B1∩B2=∅,
B
=
B
1
∪
B
2
B = B_1\cup B_2
B=B1∪B2。实际上
B
1
B_1
B1代表
B
B
B中我们感兴趣的物体(要检测的形状)对应的集合,而
B
2
B_2
B2为
B
B
B中背景部分对应的集合。
其实
B
1
B_1
B1和
B
2
B_2
B2分别代表两个结构元。输出图像由所有在
B
1
B_1
B1中匹配的像素(击中)和未在
B
2
B_2
B2中匹配的像素(击不中)组成。
代码示例:
import cv2
import numpy as np
def hit_or_miss_transform(image_path):
# 读取图像
image = cv2.imread(image_path, cv2.IMREAD_GRAYSCALE)
# 定义击中/击不中变换的核(结构元素)
kernel = np.array([
[1, 1, 1],
[1, 1, 1],
[1, 1, 1]
], dtype=np.uint8)
# 执行击中/击不中变换
transformed_image = cv2.morphologyEx(image, cv2.MORPH_HITMISS, kernel)
# 显示原始图像和变换后的图像
cv2.imshow("Original Image", image)
cv2.imshow("Transformed Image", transformed_image)
cv2.waitKey(0)
cv2.destroyAllWindows()
# 图像路径
image_path = "ch9/1.tif"
# 调用函数进行击中/击不中变换
hit_or_miss_transform(image_path)
运行结果:
可以看到噪点明显消失,且很多纹路消失掉了,这是因为设置的是一个3*3的全为1的核。那么也就意味着必须该区域全部高亮才能保留中心,所以很多地方都消失掉了。
作为对比我们再来看一组只有第一行为1的核。
可以看到,更多的地方都被保留了下来,这是因为只有第一行设置为1的核更容易满足。
基本的形态学算法
边界提取
集合A的边界
β
(
A
)
\beta (A)
β(A),可以通过先用B(结构元)对A腐蚀,然后用A减去腐蚀的结果获得。
公式如下:
β
(
A
)
=
A
−
(
A
⊖
B
)
\beta (A) = A - (A\ominus B)
β(A)=A−(A⊖B)
示例代码:
import cv2
import numpy as np
def boundary_extraction(image_path):
# 读取图像
image = cv2.imread(image_path, cv2.IMREAD_GRAYSCALE)
# 创建膨胀核和腐蚀核(结构元素)
dilate_kernel = np.ones((3, 3), np.uint8)
erode_kernel = np.ones((3, 3), np.uint8)
# 执行膨胀和腐蚀操作
dilated_image = cv2.dilate(image, dilate_kernel, iterations=1)
eroded_image = cv2.erode(image, erode_kernel, iterations=1)
# 计算边界图像
boundary_image = dilated_image - eroded_image
# 显示原始图像和边界图像
cv2.imshow("Original Image", image)
cv2.imshow("Boundary Image", boundary_image)
cv2.waitKey(0)
cv2.destroyAllWindows()
# 图像路径
image_path = "ch5/1.tif"
# 调用函数进行击中/击不中变换
boundary_extraction(image_path)
结果示例:
保留了宽度为1的边界,这是因为腐蚀的过程中,3*3的核会将最外层的给腐蚀掉,然后减去腐蚀的结果,得到的边界就是1了。
孔洞填充
公式如下:
X
k
=
(
X
k
−
1
⊕
B
)
∩
A
c
X_k = (X_{k-1}\oplus B)\cap A^c
Xk=(Xk−1⊕B)∩Ac
解释公式就是说,先找孔洞的一个点,用结构元去膨胀,然后用原图像的补集进行约束(就是求个交集),不断重复膨胀,约束直至图形不改变(即收敛)就停止,与原图求个交集,孔洞就填上了。
更具体的解释可以参见这篇博客:https://blog.csdn.net/weixin_43682519/article/details/110268506
这个更像是寻找连通集的算法,代码实现略去不表,因为这个算法看起来不太聪明的样子。
连通分量提取
公式如下:
X
k
=
(
X
k
−
1
⊕
B
)
∩
A
X_k = (X_{k-1}\oplus B)\cap A
Xk=(Xk−1⊕B)∩A
这个公式与孔洞填充的一致,只不过在该公式中,起始点
X
0
X_0
X0由原本的孔洞中的点变成了图像中的店,背景
A
c
A^c
Ac变成了
A
A
A
凸壳
首先明确几个概念:
- 凸形:集合 A A A内连接任意两点的直线都在 A A A内部,则称集合 A A A为凸形的, A A A是凸集;
- 凸壳:能够包含任意集合 S S S 的最小凸集 H H H , H H H被称为 S S S的凸壳;
- 凸缺: H − S H-S H−S ;
- 集合 A A A 的凸壳: C ( A ) C(A) C(A) 。
说人话,凸壳就是指能够完全包围给定点集的最小凸多边形。凸缺就是指凸多边形或凸多面体内部的一些点没有被包围在凸壳内部的点。
下图中外围的菱形就是里面四角星的凸壳,蓝色部分就是凸缺。
获得凸壳的算法:
X
k
i
=
(
X
k
−
1
⊛
B
i
)
∪
A
,
i
=
1
,
2
,
3
,
4
k
=
1
,
2
,
3...
X_k^i = (X_{k-1} \circledast B^i) \cup A ,\qquad i=1,2,3,4 \qquad k = 1,2,3...
Xki=(Xk−1⊛Bi)∪A,i=1,2,3,4k=1,2,3...
其中
B
i
B^i
Bi如下图所示:
即,该方法反复使用B1对A做击中击不中变换,当不再发生进一步变化时,与A求并集,得到D1,由结构元Bi(i = 2,3,4)和A进行相同的运算可得Di(i = 2,3,4)。最后,4个D的并组成了A的凸壳。
细化
细化是在图像中将二值物体和形状减小为单个像素宽的线。
定义如下:
A
⊗
B
=
A
−
(
A
⊛
B
)
=
A
∩
(
A
⊛
B
)
c
A\otimes B = A-(A\circledast B) = A \cap (A\circledast B)^c
A⊗B=A−(A⊛B)=A∩(A⊛B)c
对称的细化
A
A
A的定义如下:
A
⊗
{
B
}
=
(
(
.
.
.
(
(
A
⊗
B
1
)
⊗
B
2
)
.
.
.
)
⊗
B
n
)
A\otimes \{B\} = ((...((A \otimes B^1)\otimes B^2)...)\otimes B^n)
A⊗{B}=((...((A⊗B1)⊗B2)...)⊗Bn)
其中,
{
B
}
\{B\}
{B} 的定义如下,且
B
i
B^{i}
Bi为
B
i
−
1
B^{i-1}
Bi−1的旋转形式:
{
B
}
=
{
B
1
,
B
2
,
.
.
.
,
B
n
}
\{B\} = \{B^1,B^2,...,B^n\}
{B}={B1,B2,...,Bn}
形象的描述如下:
粗化
粗化是细化的形态学对偶,定义如下:
A
∙
B
=
A
∪
(
A
⊛
B
)
A \bullet B = A \cup (A \circledast B)
A∙B=A∪(A⊛B)
为了粗化A,通常的做法是:先对A求补
C
=
A
c
C = A^c
C=Ac,然后细化C ,然后再对结果求补。
类似于细化,也可以定义一个序列操作:
A
⊙
{
B
}
=
(
(
.
.
.
(
(
A
⊙
B
1
)
⊙
B
2
)
.
.
.
)
⊙
B
n
)
A\odot \{B\} = ((...((A \odot B^1)\odot B^2)...)\odot B^n)
A⊙{B}=((...((A⊙B1)⊙B2)...)⊙Bn)
可视化过程如下:
骨架
如图集合A和骨架S(A)的概念很简单,由该图我们可以推出:
- 如果z是S(A)的一个点,并且(D)z是A内以z为中心的最大圆盘,则不存在包含(D)z且位于A内的更大圆盘(不必以z为中心)。圆盘(D)z称为最大圆盘。
- 圆盘(D)z在两个或多个不同位置与A的边界接触。
示例代码
import cv2
import numpy as np
def image_skeletonization(image_path):
# 读取图像
image = cv2.imread(image_path, cv2.IMREAD_GRAYSCALE)
# 二值化图像
_, binary_image = cv2.threshold(image, 210, 255, cv2.THRESH_BINARY)
# 创建结构元素
kernel = cv2.getStructuringElement(cv2.MORPH_CROSS, (3, 3))
# 创建空白骨架图像
skeleton = np.zeros_like(binary_image)
# 迭代骨架提取
while True:
# 腐蚀操作
eroded = cv2.erode(binary_image, kernel)
# 膨胀后的开操作
temp = cv2.dilate(eroded, kernel)
temp = cv2.subtract(binary_image, temp)
skeleton = cv2.bitwise_or(skeleton, temp)
# 更新二值图像
binary_image = eroded.copy()
# 检查是否结束迭代
if cv2.countNonZero(binary_image) == 0:
break
# 显示原始图像和骨架图像
cv2.imshow("Original Image", image)
cv2.imshow("Skeleton Image", skeleton)
cv2.waitKey(0)
cv2.destroyAllWindows()
# 图像路径
image_path = "ch9/5.tif"
# 调用函数进行图像骨架提取
image_skeletonization(image_path)
运行结果:
如上图,可以看到基本上是提取到了图像的骨架,但是仍旧有些地方细化不够到位,推测是由于在拐角处,图像的宽度不稳定所导致的。
形态学重建
形态学重建涉及两个图像和一个结构元。一幅图像是标记图像F,一幅图像是模板G。形态学重建算法一般凭借约束和标记,迭代的对图像处理,直到收敛,达到处理自动化的目的。
形态学重建的核心是测地膨胀和测地腐蚀。
这种情况下,每一次的膨胀和腐蚀都被模板G约束。即每一次膨胀和腐蚀之后都会和G取并集或者交集。 有限数量图像的测地膨胀和测地腐蚀经过有限次数的迭代步骤之后,由于存在模板的限制,总会收敛。
测地膨胀
令F表示标记图像,令G表示模板图像。在讨论中,我们假设两幅图像都是二值图像,且F包含于G。标记图像相对于模板大小为1的测地膨胀定义为
D
G
(
1
)
=
(
F
⊕
B
)
∩
G
D_G^{(1)} = (F\oplus B)\cap G
DG(1)=(F⊕B)∩G
F相对于G的大小为n的测地膨胀定义为:
D
G
(
n
)
=
D
G
(
1
)
[
D
G
(
n
−
1
)
(
F
)
]
D_G^{(n)} = D_G^{(1)}[D_G^{(n-1)}(F)]
DG(n)=DG(1)[DG(n−1)(F)]
式中,n≥1是整数,D0=F。在这个递归公式中,每一步都是膨胀后取交集。交集运算可以保证模板G限制标记F的生长(膨胀)。下图是大小为1的测地膨胀。
测地腐蚀
类似于测地膨胀,有定义如下:
E
G
(
1
)
=
(
F
⊖
B
)
∪
G
E_G^{(1)} = (F\ominus B)\cup G
EG(1)=(F⊖B)∪G
E
G
(
n
)
=
E
G
(
1
)
[
E
G
(
n
−
1
)
(
F
)
]
E_G^{(n)} = E_G^{(1)}[E_G^{(n-1)}(F)]
EG(n)=EG(1)[EG(n−1)(F)]
在这个递归公式中,每一步都是腐蚀后取并集,以保证图像的测地腐蚀仍然大于模板图像。下图是大小为1的测地腐蚀:
重建开操作
一幅图像F的大小为n的重建开操作定义为来自F大小为n的腐蚀之后的膨胀重建。这个情况下,腐蚀会删除小的物体而之后的膨胀会恢复物体的形状。
示例代码:
import cv2
import numpy as np
src = cv2.imread("ch9/7.tif", 0)
erodeImg = None
dilateImg = None
# 构造标记图像: 采用图像的腐蚀结果作为膨胀重建的标记
element1 = cv2.getStructuringElement(cv2.MORPH_RECT, (1, 15))
erodeImg = cv2.morphologyEx(src, cv2.MORPH_ERODE, element1)
# 形态学重建
dilateImg_pre = erodeImg.copy()
tmp = np.zeros(src.shape, dtype=np.uint8)
cmp = np.zeros(src.shape, dtype=np.uint8)
element2 = cv2.getStructuringElement(cv2.MORPH_RECT, (3, 3))
n = -1
while n != src.shape[0] * src.shape[1]:
dilateImg = cv2.morphologyEx(dilateImg_pre, cv2.MORPH_DILATE, element2)
dilateImg = cv2.bitwise_and(src, dilateImg)
cmp = cv2.compare(dilateImg_pre, dilateImg, cv2.CMP_EQ)
n = cv2.countNonZero(cmp)
dilateImg_pre = dilateImg.copy()
cv2.imshow("src", src)
cv2.imshow("erodeImg", erodeImg)
cv2.imshow("dilateImg", dilateImg)
cv2.waitKey(0)
cv2.destroyAllWindows()
运行结果:
左图为原始图像,同时作为模板,中间的图为腐蚀后的图像,同时作为标记,右图为重建后的图像,可以看到恢复了很多内容,没有被恢复的部分是由于在腐蚀的时候有些字母被完全去除掉了导致该字母的位置上没有标记进而无法重建。
孔洞填充
从边界开始不断执行腐蚀,两张图求最小值操作直至结果不变。从边界开始腐蚀意味着每次只进一点点,两张图求最小值操作意味着原本被包围的区域会被恢复,这样一来能够保证原本的孔洞被完美填充。大致过程如下。
上图是创造marker的过程,对其取反即可得到目标图像
代码如下:
import numpy as np
import cv2 as cv
from matplotlib import pyplot as plt
img = cv.imread("ch9/2.tif")
# 二值化
imgray = cv.cvtColor(img, cv.COLOR_BGR2GRAY)
imgray[imgray < 100] = 0
imgray[imgray >= 100] = 255
# 原图取补得到MASK图像
mask = 255 - imgray
# 构造Marker图像
marker = np.zeros_like(imgray)
marker[0, :] = 255
marker[-1, :] = 255
marker[:, 0] = 255
marker[:, -1] = 255
marker_0 = marker.copy()
# 形态学重建
SE = cv.getStructuringElement(shape=cv.MORPH_CROSS, ksize=(3, 3))
while True:
marker_pre = marker
dilation = cv.dilate(marker, kernel=SE)
marker = np.min((dilation, mask), axis=0)
cv.imshow("marker", marker)
cv.waitKey(0)
if (marker_pre == marker).all():
break
dst = 255 - marker
filling = dst - imgray
# 显示
plt.figure(figsize=(12, 6)) # width * height
plt.subplot(2, 3, 1), plt.imshow(imgray, cmap='gray'), plt.title('src'), plt.axis("off")
plt.subplot(2, 3, 2), plt.imshow(mask, cmap='gray'), plt.title('Mask'), plt.axis("off")
plt.subplot(2, 3, 3), plt.imshow(marker_0, cmap='gray'), plt.title('Marker 0'), plt.axis("off")
plt.subplot(2, 3, 4), plt.imshow(marker, cmap='gray'), plt.title('Marker'), plt.axis("off")
plt.subplot(2, 3, 5), plt.imshow(dst, cmap='gray'), plt.title('dst'), plt.axis("off")
plt.subplot(2, 3, 6), plt.imshow(filling, cmap='gray'), plt.title('Holes'), plt.axis("off")
plt.show()
结果如下:
可以看到孔洞被很好的填充了,如果能够预先知道边界宽度,那么核应该可以设置的更大,然后运算速度就会更快一些。
边界清除
对于后续形状分析而言,从图像中提取目标是自动图像处理的基本任务。检测接触边界的目标的算法是一个很有用的功能,因为(1)它可以遮蔽图像,以便为进一步处理保留完整的目标,或者(2)它可以作为视野中出现的部分目标的一个信号。
代码实现如下:
import cv2
import numpy as np
from matplotlib import pyplot as plt
imgGray = cv2.imread("ch9/2.tif", flags=0) # flags=0 灰度图像
ret, imgBinInv = cv2.threshold(imgGray, 205, 255, cv2.THRESH_BINARY_INV) # 二值化处理 (黑色0/白色1)
imgBin = cv2.bitwise_not(imgBinInv) # 二值图像的补集 (白色背景)
# 构造标记图像:
marker = np.zeros_like(imgBin, dtype=np.uint8)
marker[0, :] = imgBin[0, :]
marker[-1, :] = imgBin[-1, :]
marker[:, 0] = imgBin[:, 0]
marker[:, -1] = imgBin[:, -1]
markerIni = marker.copy() # 标记图像: 边框 f(x,y)=I(x,y),其它为 0
# 形态学重建
element = cv2.getStructuringElement(cv2.MORPH_CROSS, (3, 3))
while True:
marker_pre = marker # 保存 F(n-1)
dilation = cv2.dilate(marker, kernel=element) # 膨胀重建
marker = cv2.bitwise_and(dilation, imgBin) # 原图像作为模板用来约束重建,按位与,有 0 得 0
if (marker_pre == marker).all(): # F(n)=F(n-1)?,判断是否达到稳定收敛状态
break # 收敛的 marker 就是需要清除的边界字符
imgRebuild = cv2.bitwise_not(imgBinInv + marker) # 对收敛的 marker 求补得到边界清除的重建结果
# 显示
plt.figure(figsize=(10, 6))
plt.subplot(221), plt.imshow(imgGray, cmap='gray'), plt.title("origin image"), plt.axis("off")
plt.subplot(222), plt.imshow(imgBinInv, cmap='gray'), plt.title("mask image"), plt.axis("off")
plt.subplot(223), plt.imshow(marker, cmap='gray'), plt.title("final marker"), plt.axis("off")
plt.subplot(224), plt.imshow(imgRebuild, cmap='gray'), plt.title("rebuild image"), plt.axis("off")
plt.tight_layout()
plt.show()
运行结果:
与填充孔洞类似,不同的是1.marker的边缘初始化为与原图像一致 2.腐蚀过后按位与运算保留腐蚀结果。这样子不断重复直至最后能够使得与边缘有连接的所有图像都被提取出来,也就是上图中的final marker。只需用原图像减去final marker即可得到最终结果。
下面是final marker的生成过程:
灰度级形态学
结构元
平坦结构元
平坦结构元的特点是它的所有元素具有相同的灰度值或二值状态。在灰度图像中,平坦结构元的所有像素具有相同的灰度级别。在二值图像中,平坦结构元的所有像素要么是黑色(0),要么是白色(1),没有其他中间值。
例如:
非平坦结构元
非平坦结构元可以是灰度图像或多值图像中的任意形状,像素的灰度级别或二值状态可以不同。
例如:
腐蚀与膨胀
腐蚀
结构元b在(x,y)处对图像f的腐蚀定义为图像f与b的重合区域的最小值,用公式表示就是:
[
f
⊖
b
]
(
x
,
y
)
=
m
i
n
(
s
,
t
)
∈
b
{
f
(
x
+
s
,
y
+
t
)
}
[f\ominus b](x,y) = \underset{(s,t)\in b}{min}\{f(x+s,y+t)\}
[f⊖b](x,y)=(s,t)∈bmin{f(x+s,y+t)}
举个简单的例子:
当结构元b在如图所示位置进行腐蚀时,将从重合位置玄虚最小值,也就是1来填充原点。
代码示例如下:
import numpy as np
import cv2
def grayscale_erosion(image, kernel):
# 获取图像的行数和列数
rows, cols = image.shape[:2]
# 获取结构元素的行数和列数
k_rows, k_cols = kernel.shape[:2]
k_center = (k_rows // 2, k_cols // 2) # 结构元素中心位置
# 创建一个与原图像大小相同的结果图像
result = np.zeros_like(image)
# 遍历图像的每个像素
for i in range(rows):
for j in range(cols):
min_val = 255 # 初始最小值设为最大灰度值
# 遍历结构元素的每个像素
for m in range(k_rows):
for n in range(k_cols):
# 计算结构元素的位置
x = i + m - k_center[0]
y = j + n - k_center[1]
# 检查结构元素是否在图像范围内
if 0 <= x < rows and 0 <= y < cols:
# 更新最小值
min_val = min(min_val, image[x, y])
# 将最小值赋给结果图像对应位置的像素
result[i, j] = min_val
return result
# 读取灰度图像
image = cv2.imread('images/lena_gray_512.tif', cv2.IMREAD_GRAYSCALE)
# 定义平坦结构元素
kernel = np.ones((5, 5), dtype=np.uint8)
# 进行灰度图像的腐蚀操作
eroded_image = grayscale_erosion(image, kernel)
# 显示原图像和腐蚀后的图像
cv2.imshow('Original Image', image)
cv2.imshow('Eroded Image', eroded_image)
cv2.waitKey(0)
cv2.destroyAllWindows()
说白了就是直接卷积一下,拿着一个滤波器去取卷积核范围内的最小值,可能只是由于在形态学上有用才取了个腐蚀的名字。腐蚀显然与滤波具有相同的特征,比如卷积核越大,细节丢失越多(放在这里就是腐蚀得越厉害)。
膨胀
与腐蚀类似,也是直接卷积一下,只不过腐蚀是拿着一个滤波器去取卷积核范围内的最小值,而膨胀是取最大值,细节略去不表。
开运算与闭运算
开运算,先腐蚀再膨胀;闭运算,先膨胀再腐蚀。代码如下:
import cv2
import numpy as np
def grayscale_opening(image, kernel):
# 进行腐蚀操作
eroded = cv2.erode(image, kernel)
# 进行膨胀操作
opened = cv2.dilate(eroded, kernel)
return opened
def grayscale_closing(image, kernel):
# 进行膨胀操作
dilated = cv2.dilate(image, kernel)
# 进行腐蚀操作
closed = cv2.erode(dilated, kernel)
return closed
# 读取灰度图像
image = cv2.imread('images/lena_gray_512.tif', cv2.IMREAD_GRAYSCALE)
# 定义平坦结构元素(3x3的矩形结构元素)
kernel = np.ones((5, 5), dtype=np.uint8)
# 进行开运算
opened_image = grayscale_opening(image, kernel)
# 进行闭运算
closed_image = grayscale_closing(image, kernel)
# 显示原图像、开运算后的图像和闭运算后的图像
cv2.imshow('Original Image', image)
cv2.imshow('Opened Image', opened_image)
cv2.imshow('Closed Image', closed_image)
cv2.waitKey(0)
cv2.destroyAllWindows()
运行结果:
可以明显看出开运算更暗一些,闭运算更亮一些,这是因为开运算先腐蚀再膨胀,腐蚀时丢失的部分细节无法在膨胀时恢复,进而导致很多小细节被去除。闭运算则是先膨胀再腐蚀,膨胀后的一些部分在后续腐蚀时也没有去除完全,进而导致更多的部分被填充。由于我们的腐蚀与膨胀一般是针对亮部而言,所以开运算更暗,闭运算更亮。
部分灰度级形态学算法
形态学平滑
就是结合开运算和闭运算来平滑图像
示例代码:
import cv2
import numpy as np
import matplotlib.pyplot as plt
def grayscale_opening(image, kernel):
# 进行腐蚀操作
eroded = cv2.erode(image, kernel)
# 进行膨胀操作
opened = cv2.dilate(eroded, kernel)
return opened
def grayscale_closing(image, kernel):
# 进行膨胀操作
dilated = cv2.dilate(image, kernel)
# 进行腐蚀操作
closed = cv2.erode(dilated, kernel)
return closed
# 读取灰度图像
image = cv2.imread('ch9/8.tif', cv2.IMREAD_GRAYSCALE)
# 定义平坦结构元素(3x3的矩形结构元素)
kernel = np.ones((5, 5), dtype=np.uint8)
# 进行开运算
opened_image = grayscale_opening(image, kernel)
# 进行闭运算
closed_image = grayscale_closing(image, kernel)
## 先开后闭
opened_closed_image = grayscale_closing(opened_image, kernel)
# 显示原图像、开运算后的图像和闭运算后的图像
titles = ['original image', 'opened image', 'closed image', 'opened_closed image']
images = [image, opened_image, closed_image, opened_closed_image]
for i in range(4):
plt.subplot(2, 2, i + 1)
plt.imshow(images[i], 'gray')
plt.title(titles[i])
plt.xticks([])
plt.yticks([])
plt.show()
可以看到,先开后闭较好的去除了噪声以及保留了原本的细节,并且平滑了原本的图像
形态学梯度
就是用膨胀的图像减去腐蚀的图像
import cv2
import numpy as np
import matplotlib.pyplot as plt
# 读取灰度图像
image = cv2.imread('images/lena_gray_512.tif', cv2.IMREAD_GRAYSCALE)
kernel = np.ones((3, 3), np.uint8)
# 膨胀
dilated_image = cv2.dilate(image, kernel, iterations=1)
# 腐蚀
eroded_image = cv2.erode(image, kernel, iterations=1)
# 形态学梯度
gradient_image = dilated_image - eroded_image
# 显示
titles = ['image', 'dilated_image', 'eroded_image', 'gradient_image']
images = [image, dilated_image, eroded_image, gradient_image]
for i in range(4):
plt.subplot(2, 2, i + 1), plt.imshow(images[i], 'gray')
plt.title(titles[i])
plt.xticks([]), plt.yticks([])
plt.show()
膨胀后的结果减去腐蚀后的结果显然能够得到梯度图像,至于为什么不是减去原图像,我认为这应该是因为梯度的定义。由于前面的章节已经介绍过梯度了,此处略去不表。
顶帽变换和底帽变换
顶帽变换:图像减去其开操作,突出亮部细节
底帽变换:图像闭操作减去其本身,突出暗部细节
import cv2
import numpy as np
import matplotlib.pyplot as plt
def grayscale_opening(image, kernel):
# 进行腐蚀操作
eroded = cv2.erode(image, kernel)
# 进行膨胀操作
opened = cv2.dilate(eroded, kernel)
return opened
def grayscale_closing(image, kernel):
# 进行膨胀操作
dilated = cv2.dilate(image, kernel)
# 进行腐蚀操作
closed = cv2.erode(dilated, kernel)
return closed
# 顶帽运算
def top_hat(image, kernel):
# 原图像与开运算之间的差值
return cv2.subtract(image, grayscale_opening(image, kernel))
# 闭帽运算
def black_hat(image, kernel):
# 闭运算与原图像之间的差值
return cv2.subtract(grayscale_closing(image, kernel), image)
# 读取原始图像
img = cv2.imread('ch9/9.tif')
# 获取自定义核
kernel = cv2.getStructuringElement(cv2.MORPH_RECT, (15, 15))
# 进行顶帽运算
tophat = top_hat(img, kernel)
# 进行闭帽运算
blackhat = black_hat(img, kernel)
# 显示图像
titles = ['Original', 'Top Hat', 'Black Hat']
images = [img, tophat, blackhat]
for i in range(3):
plt.subplot(1, 3, i + 1), plt.imshow(images[i], 'gray')
plt.title(titles[i])
plt.xticks([]), plt.yticks([])
plt.show()
开操作会去除一些细节,闭运算会填充一些细节,因此顶帽变换会突出亮部细节,而底帽变换会突出暗部细节
粒度测定
粒度测定属于判断图像中颗粒的大小分布的领域。形态学可间接用于估计颗粒的大小分布,而不需要识别并测量图像中的每个颗粒。
对于比背景亮且具有规则形状的颗粒,基本方法是使用逐渐增大的结构元对图像执行开运算。对于每一个开运算,计算该开运算中像素值的和,该和有时称为表面区域。因为开运算会降低亮特征的灰度,故表面区域会随着结构元的增大而减小。
灰度级形态学重建
就是用各种算法来对灰度级图像重建。过程与二值的基本一致,只是对象换成了灰度图像。略去不表