第九章 形态学图像处理

1.前言

本章主要从二值图像的形态学入手，了解了形态学的基础知识，不断组合类似于搭积木的样子得到更多的工具形式，并扩展到灰度图像的处理

2.预备知识

形态学有两类像素集合：目标元素和结构元（SE），具体如下：

• 目标元素：即我们常见的网格像素，通常是网格里的前景像素集合
• 结构元（SE）：按照前景像素和背景像素来定义，有时也包含无关元素

其中，与结构元有关的概念有反射和平移，令集合为$B$，则$B$相对于其原点的反射$\hat{B}$的定义为： B ^ = { w ∣ w = − b , b ∈ B } \widehat B=\{w|w=-b,b\in B\} B ={w∣w=−b,b∈B}集合$B$相对于点$z=(z_1.z_2)$定义如下： ( 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}下面是有关集合反射和平移概念图

3.腐蚀和膨胀

腐蚀

腐蚀是一种消除物体边界点的过程。在数字图像中，腐蚀操作通常用于消除噪声和平滑边界。腐蚀操作将物体的边缘向内部收缩，从而消除物体内部的孤立点。令$A$和$B$是$Z^2$的两个集合，$B$对$A$的腐蚀定义为： A ⊖ B = { z ∣ ( B ) z ⊆ A } A\ominus B=\{z|(B)_z\subseteq A\} A⊖B={z∣(B)z​⊆A}表示在平移$z$后的$B$包含于$A$的条件下，$B$对$A$的腐蚀是所有$z$的集合。

腐蚀操作的过程如下：

1. 确定集合$A$和结构元$B$，以$B$的原点从$A$的左上角开始
2. 仅当$B$完全被$A$包含时，保留$B$的原点，其余情况不算
3. 按照$A$的像素格子一步步自上而下自左而右扫描，重复步骤2
4. 将$A$的每个像素格子都扫描完，留下的$B$原点集合就是$A\ominus B$

如上图，自左而右分别是集合$A$、结构元$B$、$A\ominus B$，其中画圈的是结构元$B$的原点

膨胀

膨胀是与腐蚀相反的过程。在数字图像处理中，膨胀操作用于扩大物体边界，以检测物体的交点和断点。令$A$和$B$是$Z^2$的两个集合，$B$对$A$的膨胀定义为： A ⊕ B = { z ∣ ( B ^ ) z ∩ A ≠ ∅ } A\oplus B=\{z|(\widehat B)_z\cap A\ne \varnothing\} A⊕B={z∣(B )z​∩A=∅}表示在$\hat{B}$的前景元素与$A$的至少一个元素重合的条件下，$B$对$A$的膨胀是所有位移$z$的集合。

膨胀操作的过程如下：

1. 确定集合$A$和结构元$B$，以$B$的原点从$A$的左上角开始
2. $B$的原点在$A$的范围内时，保留$B$的所有像素格
3. 按照$A$的像素格子一步步自上而下自左而右扫描，重复步骤2
4. 将$A$的每个像素格子都扫描完，留下的所有$B$集合就是$A\oplus B$

如上图，自左而右分别是集合$A$、结构元$B$、$A\oplus B$，其中画圈的是结构元$B$的原点

4.开运算与闭运算

开运算和闭运算都是基于腐蚀和膨胀的组合运算。

开运算

开运算是一种先膨胀后腐蚀的操作。具体来说，它首先对图像进行膨胀操作，使得所有的连通区域扩大，然后对这些区域进行腐蚀操作，消除那些在膨胀过程中产生的小孔或者细线。因此，开运算可以消除图像中的小噪声，并且能够平滑物体的边界。

按照上述解释，结构元$B$对集合$A$的开运算定义为：$$A\circ B=(A\ominus B)\oplus B$$

以下是实例

# 读取图像
image = cv2.imread('./img9/Fig0905(a)(wirebond-mask).tif')

# 创建形状为正方形的结构元素
open_image = cv2.morphologyEx(image, cv2.MORPH_OPEN, cv2.getStructuringElement(cv2.MORPH_RECT, (30, 30)))#

plt.figure(figsize=(8,6))

plt.subplot(1, 2, 1)
plt.imshow(image)
plt.title('original_image')
plt.axis('off')

plt.subplot(1, 2, 2)
plt.imshow(closed_image)
plt.title('closed_image')
![请添加图片描述](https://img-blog.csdnimg.cn/ec7fba4d7b83414e993ce562e5b3ffb1.png)


plt.show()

可以看到，原图像经过开运算处理后许多细线的地方出现了明显的腐蚀效果，白色细线都没了，圆角部分也变成了平角

闭运算

闭运算是一种先腐蚀后膨胀的操作。具体来说，它首先对图像进行腐蚀操作，使得所有的连通区域缩小，然后对这些区域进行膨胀操作，填充那些在腐蚀过程中产生的空洞或者细线。因此，闭运算可以填充图像中的空洞，并且能够平滑物体的边界。

结构元$B$对集合$A$的闭运算定义为：$$A\bullet B=(A\oplus B)\ominus B$$

import cv2
import numpy as np
import matplotlib.pyplot as plt
# 读取图像
image = cv2.imread('./img9/Fig0905(a)(wirebond-mask).tif')

# 创建形状为正方形的结构元素
closing_image = cv2.morphologyEx(image, cv2.MORPH_CLOSE, cv2.getStructuringElement(cv2.MORPH_RECT, (20, 20)))


plt.figure(figsize=(8,6))

plt.subplot(1, 2, 1)
plt.imshow(image)
plt.title('original_image')
plt.axis('off')

plt.subplot(1, 2, 2)
plt.imshow(closing_image)
plt.title('closing_image')
plt.axis('off')

plt.show()

与开运算相反，闭运算使图像部分细节更加膨胀，比如一些白色细线部分的顶端都有明显的加粗

5.击中-击不中变换

击中-击不中变换（HMT）通常描述图像中的目标检测，“击中”指找到图像中的特定目标，“击不中”则相反。通过腐蚀可以展示HMT的过程，令$I$是一幅二值图像，分别使用检测前景和背景的结构元$B_1$与$B_2$，则图像$I$的HMT定义为： I ⊛ B 1 , 2 = { z ∣ ( B 1 ) z ⊆ A 和 ( B 2 ) z ⊆ A c } = ( A ⊖ B 1 ) ∩ ( A c ⊖ B 2 ) I\circledast B_{1,2}=\{z|(B_1)_z\subseteq A 和 (B_2)_z\subseteq A^c \}=(A\ominus B_1)\cap(A^c\ominus B_2) I⊛B1,2​={z∣(B1​)z​⊆A和(B2​)z​⊆Ac}=(A⊖B1​)∩(Ac⊖B2​)通俗解释下，HMT就是：在$B_1$和$B_2$分别在前景和背景中找到了匹配项的条件下，结构元$B_1$和$B_2$同时平移z的集合，图例说明如下

6.基本形态学算法

在处理二值图像时，形态学的应用主要有提取边界、连通分量、凸壳、区域骨架、区域填充、细化、粗化和修建等算法方法，接下来主要介绍其中几种

边界提取

边界提取算法的原理是：使用合适的结构元$B$腐蚀$A$，然后求两者的差，即$$\beta (A)=A-(A\ominus B)$$

孔洞填充

孔洞填充算法由集合膨胀、补集和交集为基础，在一个图像$I$中集合$A$元素为8连通边界，用前景像素填充所有的孔洞，用$X_k$表示当前的填补结果，则$$X_k=(X_{k-1}\oplus B)\cap I^c,k=1,2,\mathrm{3...}$$其中$B$表示使用的结构元，每一步的$X_k$都在前一步的基础上进行填充，如下图所示

凸壳

在欧式平面中如果一个集合（图像）是凸的，那么连接集合的任意两点的线段都在集合内部。同理数字图像中由于都是像素形成的格子，如果连接边界点形成的欧几里得线段都包含在前景点集中，那么这个图像就是凸的

用以下式子找到某一图像的凸壳$$X_k^i=(X_k^i⊛B^i)\cup X_{k-1}^i,i=1,2,3,4，k=1,2,\mathrm{3....}$$其中$B^i$表示上下左右的四个结构元，$A$的凸壳为使用每个结构元收敛时的$X_K^i$并集，即$$C(A)=\bigcup _{i=1}^4{D}^i$$

def convex_hull(img):
    img_right = img.copy()  # 用来迭代图像
    img_down = img.copy()
    img_left = img.copy()
    img_up = img.copy()

    kernel_right = np.array([[1, 0, 0], [1, -1, 0], [1, 0, 0]])  # 向右找凸壳的kernel
    kernel_down = np.array([[1, 1, 1], [0, -1, 0], [0, 0, 0]])  # 向下找凸壳的kernel
    kernel_left = np.array([[0, 0, 1], [0, -1, 1], [0, 0, 1]])  # 向左找凸壳的kernel
    kernel_up = np.array([[0, 0, 0], [0, -1, 0], [1, 1, 1]])  # 向上找凸壳的kernel

    while True:  # 右扩充
        imgRight = cv2.morphologyEx(img_right, cv2.MORPH_HITMISS, kernel_right,
                                    borderType=cv2.BORDER_DEFAULT)  # 向右做击中-击不中变换
        ret_right = img_right + imgRight  # 原图 + 变换的结果图

        if (ret_right == img_right).all():  # 两次迭代图像一致的话，处理结束
            X_right = ret_right  # 将向右扩充的图像保存
            break
        img_right = ret_right  # 将上步处理的结果返回到下一次

    while True:  # 下扩充
        imgDown = cv2.morphologyEx(img_down, cv2.MORPH_HITMISS, kernel_down, borderType=cv2.BORDER_DEFAULT)
        ret_down = img_down + imgDown

        if (ret_down == img_down).all():
            X_down = ret_down
            break
        img_down = ret_down

    while True:  # 左扩充
        imgLeft = cv2.morphologyEx(img_left, cv2.MORPH_HITMISS, kernel_left, borderType=cv2.BORDER_DEFAULT)
        ret_left = img_left + imgLeft

        if (ret_left == img_left).all():
            X_left = ret_left
            break
        img_left = ret_left

    while True:  # 上扩充
        imgUp = cv2.morphologyEx(img_up, cv2.MORPH_HITMISS, kernel_up, borderType=cv2.BORDER_DEFAULT)
        ret_up = img_up + imgUp

        if (ret_up == img_up).all():
            X_up = ret_up
            break
        img_up = ret_up

    # 将四次扩充的结果合在一块
    dst = cv2.bitwise_or(X_up, X_down)
    dst = cv2.bitwise_or(dst, X_left)
    dst = cv2.bitwise_or(dst, X_right)

    # 对凸壳进行限制，凸壳不超过前景像素的垂直和水平尺寸
    x, y = np.where(img)
    canvas = np.zeros(img.shape, dtype=np.uint8)
    canvas[x.min():x.max() + 1, y.min():y.max() + 1] = 255  # 产生一个水平和垂直限制的图像
    dst = cv2.bitwise_and(canvas, dst)  # 将凸壳限制

    return dst.astype(np.uint8)


img = cv2.imread('./img9/Fig0914(a)(licoln from penny).tif', 0)
_, img_bin = cv2.threshold(img, 50, 255, cv2.THRESH_BINARY)  # 二值化处理

dst = convex_hull(img_bin)  # 凸壳

# 显示图像
img_diff = dst - img
img_diff[img_diff > 0] = 200
img_diff += img_bin

cv2.imshow('img', np.hstack((img_bin, dst, img_diff)))
cv2.waitKey()
cv2.destroyAllWindows()
# 以上代码转载于https://blog.csdn.net/qq_44886601/article/details/127576903

如上述代码实现的效果所示，自左到右依次是原图，凸壳和两者的组合，能看到原图人像都被凸壳包含了

7.形态学重建

之前的形态学只涉及一幅图像，但形态学重建涉及两幅图像和一个结构元，一幅为标记图$F$，另一幅为模版$G$，用来约束重建

测地膨胀和腐蚀

测地膨胀指标记图像$F$与结构元膨胀后和模版$G$的交集，$F$相对于$G$大小为1和n的测地膨胀定义分别为$$D_G^{(1)}(F)=(F\oplus B)\cap G$$$$D_G^{(n)}(F)=D_G^{(1)}[D_G^{(n-1)}(F)]$$

与测地膨胀相反，测地腐蚀指标记图像$F$与结构元腐蚀后和模版$G$的并集，$F$相对于$G$大小为1和n的测地腐蚀定义分别为$$E_G^{(1)}(F)=(F\ominus B)\cup G$$$$E_G^{(n)}(F)=E_G^{(1)}[E_G^{(n-1)}(F)]$$

膨胀和腐蚀形态学重建

$F$相对于$G$的膨胀重建$R_G^D(F)$为$F$相对于$G$的测地膨胀，反复迭代至稳定，即$$R_G^D(F)=D_G^{(k)}(F)$$直到$R_G^{(k)}(F)=D_G^{(k+1)}(F)$如下图所示

$G$相对于$F$的腐蚀重建$R_G^E(F)$则相类似，只是把膨胀换成了腐蚀操作

8.灰度级形态学

首先要知道平坦和非平坦结构元含义，平坦结构元指二值灰度的结构元，非平坦表示多灰度值的结构元。结构元$b$的反射表示为$\hat{b}(x,y)=b(-x,-y)$

灰度腐蚀和膨胀

平坦结构元$b$原点在$(x,y)$，则它对图像$f$的灰度腐蚀和灰度膨胀分别为 [ f ⊖ b ] ( x , y ) = m i n ( s , t ) ∈ b { f ( x + s , y + t ) } [f\ominus b](x,y)=min_{(s,t)\in b}\{f(x+s,y+t)\} [f⊖b](x,y)=min(s,t)∈b​{f(x+s,y+t)} [ f ⊕ b ] ( x , y ) = m a x ( s , t ) ∈ b { f ( x − s , y − t ) } [f\oplus b](x,y)=max_{(s,t)\in b}\{f(x-s,y-t)\} [f⊕b](x,y)=max(s,t)∈b​{f(x−s,y−t)}而对于非平坦结构元$b_N$，它对图像$f$的灰度腐蚀和膨胀定义为 [ f ⊖ b N ] ( x , y ) = m i n ( s , t ) ∈ b ^ N { f ( x + s , y + t ) − b ^ N ( s , t ) } [f\ominus b_N](x,y)=min_{(s,t)\in \widehat b_N}\{f(x+s,y+t)-\widehat b_N(s,t)\} [f⊖bN​](x,y)=min(s,t)∈b N​​{f(x+s,y+t)−b N​(s,t)} [ f ⊕ b N ] ( x , y ) = m a x ( s , t ) ∈ b ^ N { f ( x − s , y − t ) + b ^ N ( s , t ) } [f\oplus b_N](x,y)=max_{(s,t)\in \widehat b_N}\{f(x-s,y-t)+\widehat b_N(s,t)\} [f⊕bN​](x,y)=max(s,t)∈b N​​{f(x−s,y−t)+b N​(s,t)}针对灰度膨胀和腐蚀以下是一个实例，可以看到，经过膨胀的后的图像亮特征变大，暗特征变小，左中右的细黑连线几乎已经看不见了。腐蚀的效果则相反

灰度开闭运算

在形式上，开操作和闭操作的表达式与二值图像相同$$f\circ b=(f\ominus b)\oplus b$$$$f\bullet b=(f\oplus b)\ominus b$$下图就是开闭运算的实例，可以看到在开运算中，所有亮特征都变小了，但与腐蚀结果不同，开运算对暗特征和背景的影响很小，闭运算则呈现相反的效果

灰度级形态学

基本的灰度级形态学主要有平滑、梯度、顶帽变换和底帽变换、粒度测定和纹理分割等，以下选择其中进行操作分析

形态学平滑

import cv2
import numpy as np
import matplotlib.pyplot as plt

f = cv2.imread('./img9/Fig0938(a)(cygnusloop_Xray_original).tif', 0)

plt.subplot(3, 2, 1), plt.imshow(f, cmap='gray'), plt.title('(a)Original Image'),plt.axis('off')

se = cv2.getStructuringElement(cv2.MORPH_ELLIPSE, (11, 11))

fo = cv2.morphologyEx(f, cv2.MORPH_OPEN, se)

plt.subplot(3, 2, 2), plt.imshow(fo, cmap='gray'), plt.title('(b)Open By Disk 5'), plt.axis('off')

foc = cv2.morphologyEx(fo, cv2.MORPH_CLOSE, se)

plt.subplot(3, 2, 3), plt.imshow(foc, cmap='gray'), plt.title('(c)Close B'), plt.axis('off')

focd = cv2.morphologyEx(f, cv2.MORPH_CLOSE, se)

plt.subplot(3, 2, 4), plt.imshow(focd, cmap='gray'), plt.title('(d)Close A'), plt.axis('off')

foce = cv2.morphologyEx(focd, cv2.MORPH_OPEN, se)

plt.subplot(3, 2, 5), plt.imshow(foce, cmap='gray'), plt.title('(e)Open D'), plt.axis('off')

fasf = f.copy()

for i in range(2, 6):
    se = cv2.getStructuringElement(cv2.MORPH_ELLIPSE, (2 * i + 1, 2 * i + 1))
    fasf = cv2.morphologyEx(cv2.morphologyEx(fasf, cv2.MORPH_OPEN, se), cv2.MORPH_CLOSE, se)

plt.subplot(3, 2, 6), plt.imshow(fasf, cmap='gray'), plt.title('(f)Open Close Alternate Filtering'),plt.axis('off')

plt.tight_layout()
plt.axis('off')
plt.show()

上述过程首先对原图进行开运算得到图像b，有点类似于锐化的感觉，中心的细节部分都有较明显的改善，再进行闭操作得到图c，对中心部分进行了模糊，图f是经过开闭操作选择滤波后的最终效果，整体的图像与原图相比产生了较好的平滑效果。

纹理分割

# image = cv2.imread('./img9/FigP0934(blobs_in_circular_arrangement).tif', 0)
import cv2
import numpy as np
import matplotlib.pyplot as plt

def morphological_texture_segmentation(image):
    # 定义结构元素
    kernel_radius_1 = 60
    kernel_radius_2 = 120
    kernel_1 = cv2.getStructuringElement(cv2.MORPH_ELLIPSE, (kernel_radius_1, kernel_radius_1))
    kernel_2 = cv2.getStructuringElement(cv2.MORPH_ELLIPSE, (kernel_radius_2, kernel_radius_2))
    gradient_kernel = cv2.getStructuringElement(cv2.MORPH_RECT, (3, 3))

    # 第一步：执行闭运算以去除小斑点区域
    closed_image = cv2.morphologyEx(image, cv2.MORPH_CLOSE, kernel_1)

    # 第二步：执行开运算以去除大斑点间亮间距图像
    opened_image = cv2.morphologyEx(closed_image, cv2.MORPH_OPEN, kernel_2)

    # 第三步：计算边界并叠加到原图像上
    gradient_image = cv2.morphologyEx(opened_image, cv2.MORPH_GRADIENT, gradient_kernel)
    segmented_image = cv2.add(image, gradient_image)

    return closed_image, opened_image, segmented_image

# 加载原始图像
image = cv2.imread('./img9/Fig0943(a)(dark_blobs_on_light_background).tif', 0)

# 执行形态学纹理分割
closed_image, opened_image, segmented_image = morphological_texture_segmentation(image)

# 显示结果
plt.figure(figsize=(12, 8))

plt.subplot(2, 2, 1)
plt.imshow(cv2.cvtColor(image, cv2.COLOR_BGR2RGB))
plt.title('(a) Original')

plt.subplot(2, 2, 2)
plt.imshow(closed_image, cmap='gray')
plt.title('(b) Closed')

plt.subplot(2, 2, 3)
plt.imshow(opened_image, cmap='gray')
plt.title('(c) Opened')

plt.subplot(2, 2, 4)
plt.imshow(cv2.cvtColor(segmented_image, cv2.COLOR_BGR2RGB))
plt.title('(d) Segmented')

plt.tight_layout()
plt.show()

代码首先使用半径为60像素的圆盘形结构元对原图像进行闭运算，然后使用半径120像素的圆盘形结构元对闭运算结果进行开运算。接下来，使用3x3大小的矩形结构元计算开运算结果和原图像之间的梯度，并将边界叠加到原图像上得到分割图像。最后，使用matplotlib将四张图像显示在一张图上。最后的结果是将大斑点和小斑点分开了

9.章总结

本章学习了形态学的相关知识，从二值图像延伸到灰度图像。腐蚀和膨胀是基础，在腐蚀和膨胀的组合之下形成了开闭运算，并可以组合成更为复杂的工具。在二值图像中，学习到了运用形态学能够执行平滑、边缘检测、凸壳等任务，以及形态学重建，和之前的章节对应了起来，并延伸至灰度图像的相关处理，并简要实现了纹理分割的算法实例