击中击不中变换 (Hit-or-miss transform)¶

击中击不中变换是在二值图像中寻找特定形状的强大工具。其有效性源于其独特的方法：它不仅寻找对象的形状，还同时检查对象的前景模式及其紧邻的背景模式。

为此，它使用一个由两个不同且不相交部分 $E_1$ 和 $E_2$ 组成的复合结构元素 $E$ ：

$E_1$ 是 “击中”元素。它定义了我们想要找到的对象（前景）的形状。
$E_2$ 是 “击不中”元素。它定义了必须围绕该对象的背景像素的模式。

只有当“击中”元素 ( $E_1$ ) 完美地覆盖在前景像素上，并且“击不中”元素 ( $E_2$ ) 同时完美地覆盖在背景像素上时，才会在特定位置找到匹配。这种严格的双重条件确保只检测到确切的所需形状，使得该变换非常精确。

击中击不中变换的结果是由以下两个集合的交集给出的：

图像被 $E_1$ 腐蚀的结果： $I \ominus E_1$ ，
图像背景被 $E_2$ 腐蚀的结果： $I^\mathrm{c} \ominus E_2$

由两个结构元素 $E_1$ 和 $E_2$ 进行的击中击不中变换记为 $I \otimes (E_1, E_2)$ :

I \otimes (E_1, E_2) = (I \ominus E_1) \cap (I^\mathrm{c} \ominus E_2) \\ = (I \ominus E_1) \cap (I \oplus E_2)^\mathrm{c}

(1)

应用¶

上面介绍的四种形态学工具是进行更复杂处理的基本操作。本节我们重点介绍一些特定的处理方法，其他方法可参见 [Bloch 2005, chapter 6]。

顶帽变换¶

顶帽变换是一种强大的形态学运算，用于从图像中提取小元素和细节。它对于校正不均匀光照以及在变化的背景中检测小的亮或暗特征特别有效。

如果在顶帽滤波器图像的示意图中绘制内部和环形区域，滤波器名称的由来就显而易见了。内部区域是帽冠，阈值是其高度，而周围的环形区域是帽檐。此操作特别适用于寻找傅里叶变换功率谱中的尖峰。

顶帽变换有两种类型：

白顶帽变换定义为原始图像 $I$ 与其开运算之差。它提取比结构元素小的明亮特征。
$T_w = I - (I \circ E)$
(2)
黑顶帽变换（或底帽变换）定义为图像的闭运算与原始图像之差。它提取比结构元素小的黑暗特征。
$T_b = (I \bullet E) - I$
(3)

使用顶帽变换的关键是选择一个比感兴趣的特征大（因此它们会被开运算或闭运算移除）但比背景变化小的结构元素。

Figure 2 显示了白顶帽变换（只保留比结构元素窄的白色对象）和黑顶帽变换（只保留比结构元素窄的黑色对象）的效果。请注意，原始图像中存在的结构元素没有被保留，因为它不比自身窄。结果还显示了许多非常小的伪影，可以通过丢弃尺寸小于预定义阈值的对象来轻松移除。

Figure 2:使用半径为3像素的圆形作为结构元素的白顶帽和黑顶帽变换示例（后者在原始图像的右下角可见）。

以下示例更清晰地演示了校正不均匀光照，这是顶帽变换的一个主要应用。创建了一个带有渐变背景的合成图像，并使用该变换来分离小的亮或暗特征。

import numpy as np
import matplotlib.pyplot as plt
from skimage.morphology import white_tophat, black_tophat, square

# --- 白顶帽示例 ---
# 创建一个带有渐变背景的合成图像
image_bright_spots = np.zeros((200, 400), dtype=np.uint8)
gradient = np.linspace(50, 150, 400, dtype=np.uint8)
image_bright_spots[:, :] = gradient
# 添加小的亮点（待检测的特征）
image_bright_spots[50:55, 50:55] = 255
image_bright_spots[150:155, 300:305] = 255
image_bright_spots[100:103, 200:203] = 220

# --- 黑顶帽示例 ---
# 创建一个类似的图像，但带有暗点
image_dark_spots = np.copy(image_bright_spots)
image_dark_spots[50:55, 100:105] = 10
image_dark_spots[150:155, 250:255] = 20

# 定义一个比斑点大的结构元素
selem = square(21)

# 应用白顶帽变换寻找亮点
tophat_white = white_tophat(image_bright_spots, selem)

# 应用黑顶帽变换寻找暗点
tophat_black = black_tophat(image_dark_spots, selem)

# 绘图
fig, axes = plt.subplots(2, 2, figsize=(12, 8), sharex=True, sharey=True)
ax = axes.ravel()

ax[0].imshow(image_bright_spots, cmap='gray', vmin=0, vmax=255)
ax[0].set_title("带亮点的原始图像")

ax[1].imshow(tophat_white, cmap='gray', vmin=0, vmax=255)
ax[1].set_title("白顶帽输出")

ax[2].imshow(image_dark_spots, cmap='gray', vmin=0, vmax=255)
ax[2].set_title("带暗点的原始图像")

ax[3].imshow(tophat_black, cmap='gray', vmin=0, vmax=255)
ax[3].set_title("黑顶帽输出")

for a in ax:
    a.axis('off')

plt.tight_layout()
plt.show()

颗粒分析¶

颗粒分析是一种通过连续选择二值图像中尺寸递增的颗粒来进行分析的方法。它包括对图像 $I$ 进行一系列尺寸递增的结构元素 $E_k$ 的连续开运算（通常， $E_k = \underbrace{E\oplus \ldots \oplus E} _{k\text{ 次}}$ ，其中 $E$ 是一个圆盘）。

在开运算的每个阶段，更精细的细节被相继消除，图像的面积也随之减小（面积是对象中像素的总数）。

Figure 4 显示了通过使用不同尺寸的圆盘进行开运算，在Figure 3的二值图像上获得的颗粒分析函数。当开运算的尺寸与大多数对象的特征尺寸相对应时，曲线中会出现一个跳跃。

要计算颗粒度，我们对原始图像进行一系列形态学开运算。每次开运算都使用一个半径不断增加的圆盘形结构元素。每次开运算后，我们计算剩余明亮区域的总面积。该面积随圆盘半径变化的曲线图即为颗粒度曲线，它揭示了图像中物体尺寸的分布情况。

import numpy as np
import matplotlib.pyplot as plt
from skimage.morphology import opening, square
from skimage.draw import disk as draw_disk

# 创建一个包含不同尺寸颗粒的图像
image = np.zeros((200, 200), dtype=np.uint8)
rr, cc = draw_disk((50, 50), 10); image[rr, cc] = 1
rr, cc = draw_disk((100, 150), 20); image[rr, cc] = 1
rr, cc = draw_disk((150, 70), 15); image[rr, cc] = 1
rr, cc = draw_disk((30, 160), 8); image[rr, cc] = 1

# 执行颗粒分析
max_size = 25
areas = []
for i in range(1, max_size + 1):
    selem = square(i)
    opened_image = opening(image, selem)
    areas.append(np.sum(opened_image))

# 颗粒分析曲线是面积的导数
gran_curve = -np.diff(areas)

# 绘图
fig, axes = plt.subplots(1, 2, figsize=(12, 5))
ax = axes.ravel()

ax[0].imshow(image, cmap='gray_r', interpolation='nearest')
ax[0].set_title('带颗粒的原始图像')
ax[0].axis('off')

ax[1].plot(range(1, max_size), gran_curve)
ax[1].set_title('颗粒分析曲线（模式谱）')
ax[1].set_xlabel('结构元素大小')
ax[1].set_ylabel('面积变化')

plt.tight_layout()
plt.show()

击中击不中变换¶

击中击不中变换用于检测特定形状的对象。与迄今为止所见的变换不同，击中击不中变换检查的配置中，某些像素验证与对象的关系，而另一些像素则验证与背景（对象补集）的关系。因此，该变换的结构元素 $E$ 由两个元素 $E_1$ 和 $E_2$ 组成（不相交但具有相同的原点）。

击中击不中变换的结果是由以下两个集合的交集给出的：

图像被 $E_1$ 腐蚀的结果： $I \ominus E_1$ ，
图像背景被 $E_2$ 腐蚀的结果： $I^\mathrm{c} \ominus E_2$

由两个结构元素 $E_1$ 和 $E_2$ 进行的击中击不中变换记为 $I \otimes (E_1, E_2)$ :

I \otimes (E_1, E_2) = (I \ominus E_1) \cap (I^\mathrm{c} \ominus E_2) \\ = (I \ominus E_1) \cap (I \oplus E_2)^\mathrm{c}

(4)

通过结构元素 E_1 和 E_2 应用于图像 I 的击中击不中变换示例。
结构元素的原点由蓝点标记。
I^\mathrm{c} \ominus E_2 中的一些白色像素来自被认为是白色的 I^\mathrm{c} 外部的背景。 — Figure 5:通过结构元素 $E_1$ 和 $E_2$ 应用于图像 $I$ 的击中击不中变换示例。结构元素的原点由蓝点标记。 $I^\mathrm{c} \ominus E_2$ 中的一些白色像素来自被认为是白色的 $I^\mathrm{c}$ 外部的背景。

import numpy as np
import matplotlib.pyplot as plt
from scipy.ndimage import binary_hit_or_miss

# 创建一个示例图像
image = np.array([
    [0, 0, 0, 0, 0, 0, 0, 0],
    [0, 1, 1, 1, 0, 0, 1, 0],
    [0, 1, 1, 1, 0, 0, 0, 0],
    [0, 1, 1, 1, 0, 1, 0, 0],
    [0, 0, 1, 0, 0, 0, 0, 0],
    [0, 0, 0, 0, 0, 1, 1, 0],
    [0, 1, 0, 0, 0, 1, 1, 0],
    [0, 0, 0, 0, 0, 0, 0, 0]
], dtype=int)

# 定义一个结构元素来寻找右下角
# 1 表示前景，-1 表示背景，0 表示不关心
selem = np.array([
    [1, 0, -1],
    [1, 1, -1],
    [1, 0, -1]
])

# 执行击中击不中变换
output = binary_hit_or_miss(image, structure=selem).astype(np.int)

# 绘图
fig, axes = plt.subplots(1, 2, figsize=(10, 5))
ax = axes.ravel()

ax[0].imshow(image, cmap='gray_r', interpolation='nearest')
ax[0].set_title('原始图像')

ax[1].imshow(output, cmap='gray_r', interpolation='nearest')
ax[1].set_title('击中击不中变换输出 (角点)')

for a in ax:
    a.axis('off')

plt.show()

骨架化¶

骨架化是将二值图像中的前景区域简化为一个骨架残余的过程，该残余在很大程度上保留了原始区域的范围和连通性。生成的骨架是对象的一像素宽表示，类似于中轴。

骨架是形状分析和模式识别的有力工具，因为它提供了一个对象的简化但拓扑等价的表示。它捕捉了基本的几何结构，包括连通分量的数量和孔洞。

计算骨架的一种方法是通过形态学细化过程，该过程迭代地腐蚀掉边界像素而不改变对象的拓扑结构。另一种由Lantuéjoul提出的定义，通过对一系列形态学开运算结果求和来构造骨架。

从一个结构元素 $E$ 出发，我们定义 $\forall k \in \mathbb{N}^*$ :

E_k = \underbrace{E\oplus \ldots \oplus E}_{k\text{ 次}}

(5)

图像 $I$ 的骨架 $S(I)$ 则定义为连续腐蚀与其开运算之差的并集：

S(I) = \bigcup_{k \in \mathbb{N}^*} \left( (I \ominus E_k) \setminus (I \ominus E_k) \circ E \right)

(6)

该公式有效地找到了可以在不同尺度下适应形状内部的最大圆盘的中心。

import numpy as np
import matplotlib.pyplot as plt
from skimage.morphology import skeletonize
from skimage.data import horse

# 使用一个示例二值图像
image = horse()

# 执行骨架化
skeleton = skeletonize(image)

# 绘图
fig, axes = plt.subplots(1, 2, figsize=(10, 5), sharex=True, sharey=True)
ax = axes.ravel()

ax[0].imshow(image, cmap='gray_r', interpolation='nearest')
ax[0].set_title('原始图像')

ax[1].imshow(skeleton, cmap='gray_r', interpolation='nearest')
ax[1].set_title('骨架')

for a in ax:
    a.axis('off')

plt.tight_layout()
plt.show()