数字图像处理学习笔记九形态学处理

形态学处理

（一）腐蚀与膨胀
- 1.1 结构元素
- 1.2 图像膨胀
- 1.3 图像腐蚀
- 1.4 图像膨胀与腐蚀实验
（二）开闭运算
- 2.1 概念介绍
- 2.2 实验
- 2.3 击中与击不中
（三）形态学运算
- 3.1 顶帽运算和黑帽运算
- 3.2 边界提取
- 3.3 检测拐角
- 3.4 孔洞填充
- 3.5 连通区域
- 3.6 删除小块区域
- 3.7 细化（骨架提取）

前言：首先明确一点形态学操作都是针对二值图像，在灰度图像的形态学处理也是先根据一些前置操作转化为二值图像（如高帽运算后再阈值化），这是因为形态学的任何操作都是需要击中这个概念，膨胀开闭其实可以看成包含击中的是一个点，其他操作可以看成击中的是一个模版。

（一）腐蚀与膨胀

1.1 结构元素

形态学处理的核心就是定义结构元素，在OpenCV-Python中，可以使用其自带的getStructuringElement函数，也可以直接使用NumPy的ndarray来定义一个结构元素。如下：
十字架：element = cv2.getStructuringElement(cv2.MORPH_CROSS,(5,5))

还可以定义椭圆/矩形等：
椭圆：cv2.getStructuringElement(cv2.MORPH_ELLIPSE,(5,5))
矩形：cv2.getStructuringElement(cv2.MORPH_RECT,(5,5))

图像的膨胀（Dilation）和腐蚀（Erosion）是两种基本的形态学运算，主要用来寻找图像中的极大区域和极小区域。其中膨胀类似于“领域扩张”，将图像中的高亮区域或白色部分进行扩张，其运行结果图比原图的高亮区域更大；腐蚀类似于“领域被蚕食”，将图像中的高亮区域或白色部分进行缩减细化，其运行结果图比原图的高亮区域更小。

1.2 图像膨胀

膨胀的运算符是“⊕”，其定义如下：

表示用B来对图像A进行膨胀处理，其中B是一个卷积模板或卷积核，其形状可以为正方形或圆形，通过模板B与图像A进行卷积计算，扫描图像中的每一个像素点，用模板元素与二值图像元素做“与”运算，如果都为0，那么目标像素点为0，否则为1。从而计算B覆盖区域的像素点最大值，并用该值替换参考点的像素值实现膨胀。下图是将左边的原始图像A膨胀处理为右边的效果图A⊕B。

1.3 图像腐蚀

腐蚀的运算符是“－”，其定义如下：

该公式表示图像A用卷积模板B来进行腐蚀处理，通过模板B与图像A进行卷积计算，得出B覆盖区域的像素点最小值，并用这个最小值来替代参考点的像素值。如图所示，将左边的原始图像A腐蚀处理为右边的效果图 A-B。

腐蚀和膨胀关于补集和反射操作呈对偶关系：

1.4 图像膨胀与腐蚀实验

（1）膨胀处理
1.基础理论
图像膨胀是腐蚀操作的逆操作，类似于“领域扩张”，将图像中的高亮区域或白色部分进行扩张，其运行结果图比原图的高亮区域更大，线条变粗了，主要用于去噪。
(1) 图像被腐蚀后，去除了噪声，但是会压缩图像。
(2) 对腐蚀过的图像，进行膨胀处理，可以去除噪声，并且保持原有形状。
(2)卷积核
卷积核是腐蚀中的关键数组，采用numpy库可以生成。卷积核的中心点逐个像素扫描原始图像，如下图所示：

被扫描到的原始图像中的像素点，当卷积核对应的元素值只要有一个为1时，其值就为1，否则为0。
2.函数原型
图像膨胀主要使用的函数为dilate，其原型如下：
dst = cv2.dilate(src, kernel, iterations)
参数dst表示处理的结果，src表示原图像，kernel表示卷积核，iterations表示迭代次数。下图表示5*5的卷积核，可以采用函数 np.ones((5,5), np.uint8) 构建。
迭代次数默认是1，表示进行一次膨胀，也可以根据需要进行多次迭代，进行多次膨胀。通常进行1次膨胀即可。
3.代码实现

#设置卷积核
kernel = np.ones((5,5), np.uint8)
#图像膨胀处理
erosion = cv2.dilate(src, kernel)
#显示图像
cv2.imshow("src", src)
cv2.imshow("result of expand", erosion)

4.处理效果：
（二）腐蚀处理
1.基础理论
形态学转换主要针对的是二值图像（0或1）。图像腐蚀类似于“领域被蚕食”，将图像中的高亮区域或白色部分进行缩减细化，其运行结果图比原图的高亮区域更小。其主要包括两个输入对象：
(1)二值图像
(2)卷积核
卷积核是腐蚀中的关键数组，采用numpy库可以生成。卷积核的中心点逐个像素扫描原始图像，如下图所示：
被扫描到的原始图像中的像素点，只有当卷积核对应的元素值均为1时，其值才为1，否则其值修改为0。换句话说，遍历到的黄色点位置，其周围全部是白色，保留白色，否则变为黑色，图像腐蚀变小。
2. 函数原型
图像腐蚀主要使用的函数为erode，其原型如下：
dst = cv2.erode(src, kernel, iterations)
参数dst表示处理的结果，src表示原图像，kernel表示卷积核，iterations表示迭代次数。下图表示5*5的卷积核，可以采用函数 np.ones((5,5), np.uint8) 构建。
迭代次数默认是1，表示进行一次腐蚀，也可以根据需要进行多次迭代，进行多次腐蚀。
3.代码实现

#设置卷积核
kernel = np.ones((5,5), np.uint8)#图像腐蚀处理
erosion = cv2.erode(src, kernel)#显示图像
cv2.imshow("src", src)
cv2.imshow("result of corrode", erosion)

4.处理效果：
（三）膨胀与腐蚀处理实例

处理后：
膨胀：

腐蚀：

实验总结：
膨胀操作类似于“领域扩张”，原二值图像中分为前景和背景，膨胀操作就是把前景区域扩大；腐蚀操作类似于“领域压缩”，图像中高亮的前景部分的区域会缩小。

（二）开闭运算

膨胀与腐蚀使边缘形状发生了变化，膨胀发生了扩张，腐蚀发生了收缩，目标物体变形，对识别时的特征提取会造成影响，可以用开闭操作来补偿。

2.1 概念介绍

开运算和闭运算就是将腐蚀和膨胀按照一定的次序进行处理。
但这两者并不可逆的，即先开后闭并不能得到原先的图像。
为了获取图像中的主要对象：对一副二值图连续使用闭运算和开运算，或者消除图像中的噪声，也可以对图像先用开运算后用闭运算，不过这样也会消除一些破碎的对象。

开运算：先腐蚀后膨胀，用于移除由图像噪音形成的斑点；
闭运算：先膨胀后腐蚀，用来连接被误分为许多小块的对象；

2.2 实验

（一）开运算：
定义矩形结构卷积核元素：

kernel = cv2.getStructuringElement(cv2.MORPH_RECT,(3,3))

迭代系数取1或3对应的开运算：

opened1 = cv2.morphologyEx(gray_res, cv2.MORPH_OPEN, kernel,iterations=1)     #开运算1
opened2 = cv2.morphologyEx(gray_res, cv2.MORPH_OPEN, kernel,iterations=3)     #开运算2

处理效果：

迭代系数为1时：
迭代系数为2时：

结果分析： 开运算看上去把细微连在一起的两块目标分开了，能够除去孤立的小点，毛刺和小桥，而总的位置和形状不变，迭代系数越大，处理效果越有卷积核的影子。

（二）闭运算
定义矩形结构卷积核元素：

kernel = cv2.getStructuringElement(cv2.MORPH_RECT,(3,3))

迭代系数取1或3对应的闭运算：

closed1 = cv2.morphologyEx(gray_res, cv2.MORPH_CLOSE, kernel,iterations=1)    #闭运算1
closed2 = cv2.morphologyEx(gray_res, cv2.MORPH_CLOSE, kernel,iterations=3)    #闭运算2

迭代系数为1时：

迭代系数为3时：

结果分析：
闭运算看上去将两个细微连接的图块封闭在一起，能够填平小湖（即小孔），弥合小裂缝，而总的位置和形状不变，迭代系数越大，处理效果越有卷积核的影子。

开运算总结：

（１）开运算能够除去孤立的小点，毛刺和小桥，而总的位置和形状不变。
（２）开运算是一个基于几何运算的滤波器。
（３）结构元素大小的不同将导致滤波效果的不同。
（４）不同的结构元素的选择导致了不同的分割，即提取出不同的特征。

闭运算总结：

（1）闭运算能够填平小湖（即小孔），弥合小裂缝，而总的位置和形状不变。
（2）闭运算是通过填充图像的凹角来滤波图像的。
（3）结构元素大小的不同将导致滤波效果的不同。
（4）不同结构元素的选择导致了不同的分割。

2.3 击中与击不中

击中击不中变换（HMT），HMT变换可以同时探测图像的内部和外部。研究解决目标图像识别和模式识别等领域，在处理目标图像和背景的关系上能够取得更好的效果。
A中对B进行的匹配（击中）表示为：

图示理解：

（三）形态学运算

3.1 顶帽运算和黑帽运算

（1）顶帽（Top Hat）
顶帽运算（Top Hat）又常常被译为”礼帽“运算。为原图像“开运算“的结果图之差，数学表达式如下：

因为开运算带来的结果是放大了裂缝或者局部低亮度的区域，因此，从原图中减去开运算后的图，得到的效果图突出了比原图轮廓周围的区域更明亮的区域，且这一操作和选择的核的大小相关。

顶帽运算往往用来分离比邻近点亮一些的斑块。当一幅图像具有大幅的背景的时候，而微小物品比较有规律的情况下，可以使用顶帽运算进行背景提取。
代码处理：

TOPHAT_img = cv2.morphologyEx(original_img, cv2.MORPH_TOPHAT, kernel)

处理效果：
处理效果分析：
从原图中减去开运算后的图，得到的效果图突出了比原图轮廓周围的区域更明亮的区域。
（2）黑帽（Black Hat）
黑帽（Black Hat）运算为”闭运算“的结果图与原图像之差。数学表达式为：

黑帽运算后的效果图突出了比原图轮廓周围的区域更暗的区域，且这一操作和选择的核的大小相关。
所以，黑帽运算用来分离比邻近点暗一些的斑块。非常完美的轮廓效果图：

BLACKHAT_img = cv2.morphologyEx(original_img, cv2.MORPH_BLACKHAT, kernel)

处理效果：
处理效果分析： 突出了比原图轮廓周围的区域更暗的区域，轮廓效果较好。

3.2 边界提取

形态学检测边缘的原理很简单：在膨胀时，图像中的物体会想周围“扩张”；腐蚀时，图像中的物体会“收缩”。由于这两幅图像其变化的区域只发生在边缘。所以这时将两幅图像相减，得到的就是图像中物体的边缘。
实验代码：
确定卷积核cv2.MORPH_RECT,(3, 3)；
腐蚀erode_img，与膨胀dilate_img；

kernel = cv2.getStructuringElement(cv2.MORPH_RECT,(3, 3))
dilate_img = cv2.dilate(image, kernel)
erode_img = cv2.erode(image, kernel)

图像相减获得边界：

absdiff_img = cv2.absdiff(dilate_img,erode_img);
retval, threshold_img = cv2.threshold(absdiff_img, 40, 255, cv2.THRESH_BINARY);
result = cv2.bitwise_not(threshold_img);

处理效果：
例1：
例2

结果分析： 对于线条密集的图像，虽然能有效提取边界，但视觉上出现重叠与一定的粗线条现象，部分区域影响着判断；对于线条明快而清晰的图像，处理效果很好，边界提取的很不错。

3.3 检测拐角

原理:先用十字形的结构元素膨胀像素，这种情况下只会在边缘处“扩张”，角点不发生变化。接着用菱形的结构元素腐蚀原图像，导致只有在拐角处才会“收缩”，而直线边缘都未发生变化。
第二步是用X形膨胀原图像，角点膨胀的比边要多。这样第二次用方块腐蚀时，角点恢复原状，而边要腐蚀的更多。所以当两幅图像相减时，只保留了拐角处。
实验代码：
构造5×5的结构元素，分别为十字形、菱形、方形和X型

cross = cv2.getStructuringElement(cv2.MORPH_CROSS,(5, 5))
diamond = cv2.getStructuringElement(cv2.MORPH_RECT,(5, 5))
diamond[0, 0] = 0
diamond[0, 1] = 0
diamond[1, 0] = 0
diamond[4, 4] = 0
diamond[4, 3] = 0
diamond[3, 4] = 0
diamond[4, 0] = 0
diamond[4, 1] = 0
diamond[3, 0] = 0
diamond[0, 3] = 0
diamond[0, 4] = 0
diamond[1, 4] = 0
square = cv2.getStructuringElement(cv2.MORPH_RECT,(5, 5))  #构造方形结构元素
x = cv2.getStructuringElement(cv2.MORPH_CROSS,(5, 5))dilate_cross_img = cv2.dilate(image,cross)                #使用cross膨胀图像
erode_diamond_img = cv2.erode(dilate_cross_img, diamond)  #使用菱形腐蚀图像dilate_x_img = cv2.dilate(image, x)                       #使用X膨胀原图像
erode_square_img = cv2.erode(dilate_x_img,square)         #使用方形腐蚀图像

将两幅闭运算的图像相减获得角：

result = cv2.absdiff(erode_square_img, erode_diamond_img)          #将两幅闭运算的图像相减获得角
retval, result = cv2.threshold(result, 40, 255, cv2.THRESH_BINARY)

在原图上用半径为5的圆圈将点标出：

for j in range(result.size):y = int(j / result.shape[0])x = int(j % result.shape[0])if result[x, y] == 255:                                        #result[] 只能传入整型cv2.circle(image,(y,x),5,(255,0,0))

处理效果：
例1：
例2：

效果分析： 如图可见，对拐角的检测很敏感，当出现拐角时可以有效检测到，标注出来，效果不错。

3.4 孔洞填充

孔洞定义：被前景（白色）连通域包围的封闭的背景（黑色）区域，不限于圆形。如图所示

原理：
以原图像的补集作为Mask，用来限制膨胀结果；以带有白色边框的黑色图像为初始Marker，用SE对其进行连续膨胀，直至收敛；最后对Marker取补即得到最终图像，与原图相减可得到填充图像。

实验代码：
Marker 图像：包含变换的起点，将被连续膨胀，直至收敛
Mask 图像：用来约束膨胀结果，即 Mask >= Marker
结构单元（Structuring Element，SE）：定义连通性

# 构造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 = markerdilation = cv.dilate(marker, kernel=SE)marker = np.min((dilation, mask), axis=0)if (marker_pre == marker).all():break
dst = 255 - marker
filling = dst - imgray

处理效果：
例1：
例2：
例3：
选用不同尺寸的SE，进行填充对比：ksize=3x3（上），ksize=7x7（下）
结果分析： 如果选择的结构单元过大，膨胀操作会越过边界，膨胀到孔洞中，导致部分孔洞填充失败。

3.5 连通区域

连通分量概念：
无向图的极大连通子图称为的连通分量( Connected Component)。任何连通图的连通分量只有一个，即是其自身，非连通的无向图有多个连通分量
邻域定义：
坐标为(x,y)的像素p有两个水平和两个垂直相邻像素，它们的坐标分别为(x+1,y)、(x-1,y)、(x,y+1)、(x,y-1)，这4个相邻像素的集合表示为N4（p）N_{4}（p）N4（p）。p的4个对角线相邻像素的坐标分别是(x-1,y-1)、(x+1,y-1)、(x-1,y++1)、(x+1,y+1)，它们称为 ND（p)N_{D}（p)ND（p)。上面两种邻域的并集是p的8相邻像素，表示为N8（p）N_{8}（p）N8（p）。
邻接定义：
如果q∈N4（p）q \in N_{4}（p）q∈N4（p），则q和p就被称为4邻接的；同样，q若属于8邻域、D邻域内，则q和p就分别是8邻接和D邻接。

如果在像素p和q之间存在一条4连接通路，就称这两个前景像素为4连接的；同理，如果在像素p和q之间存在一条8连接通路，就称这两个前景像素为8连接的。连通分量用通路来定义，而通路的定义按照所用邻接的类型而定。通俗来说，连通分量的性质决定我们所选择的邻接方式，常见的为4邻接和8邻接，同时邻接方式还决定图像中连通分量的数目。

8邻接实验代码：
编写一个函数来生成原始二值图像：

def microstructure(l=256):n = 5x, y = np.ogrid[0:l, 0:l]  # 生成网络mask = np.zeros((l, l))generator = np.random.RandomState(1)  # 随机数种子points = l * generator.rand(2, n ** 2)mask[(points[0]).astype(np.int), (points[1]).astype(np.int)] = 1mask = ndi.gaussian_filter(mask, sigma=l / (4. * n))  # 高斯滤波return mask > mask.mean()

8连通区域标记：

labels = measure.label(data, connectivity=2)
dst = color.label2rgb(labels)  # 根据不同的标记显示不同的颜色
print('regions number:', labels.max() + 1)  # 显示连通区域块数(从0开始标记)

处理效果：

3.6 删除小块区域

当只需要一些大块区域，那些零散的、小块的区域，需要删除掉，则可以使用morphology子模块的remove_small_objects（)函数，返回删除了小块区域的二值图像。
函数格式：skimage.morphology.remove_small_objects(ar, min_size=64, connectivity=1, in_place=False)

参数：

ar: 待操作的bool型数组
min_size: 最小连通区域尺寸，小于该尺寸的都将被删除。默认为64
connectivity: 邻接模式，1表示4邻接，2表示8邻接
in_place: bool型值，如果为True,表示直接在输入图像中删除小块区域，否则进行复制后再删除。默认为False.

实验代码：将面积小于300的小块区域删除（由1变为0）

def microstructure(l=256):n = 5x, y = np.ogrid[0:l, 0:l]  # 生成网络mask = np.zeros((l, l))generator = np.random.RandomState(1)  # 随机数种子points = l * generator.rand(2, n ** 2)mask[(points[0]).astype(np.int), (points[1]).astype(np.int)] = 1mask = ndi.gaussian_filter(mask, sigma=l / (4. * n))  # 高斯滤波return mask > mask.mean()
data = microstructure(l=128)  # 生成测试图片
dst = morphology.remove_small_objects(data, min_size=300, connectivity=1)fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(8, 4))
ax1.imshow(data, plt.cm.gray, interpolation='nearest')
ax2.imshow(dst, plt.cm.gray, interpolation='nearest')fig.tight_layout()

处理效果：

3.7 细化（骨架提取）

图像细化：
图像细化主要是针对二值图而言，所谓骨架，可以理解为图像的中轴，，一个长方形的骨架，是它的长方向上的中轴线，圆的骨架是它的圆心，直线的骨架是它自身，孤立点的骨架也是自身。我们来看看典型的图形的骨架（用粗线表示）
细化的算法有很多种，但比较常用的算法是查表法。细化是从原来的图中去掉一些点，但仍要保持原来的形状。实际上是保持原图的骨架。判断一个点是否能去掉是以8个相邻点（八连通）的情况来作为判据的，具体判据为：
1，内部点不能删除
2，鼓励点不能删除
3，直线端点不能删除
4，如果P是边界点，去掉P后，如果连通分量不增加，则P可删除

处理二值图像，背景为黑色=0，要细化的前景物体像素值=255。对于Hilditch算法来说，它并不是一个完全的并行算法，而是串行并行相结合。当前像素是否是能够删除的骨架点，不仅是由它周围的8邻域决定，而且和前面像素的判定结果有关。一个像素判定为可以删除，我们并不直接删除它，而是在目地图像中设置像素值为GRAY=128，这个信息可能会影响之后其它像素的判定。当图像一次扫描迭代完成后，我们把所有置为GRAY的像素设置为0，从而删除它。
在每行水平扫描的过程中，先判断每一点的左右邻居，如果都是黑点，则该点不做处理。另外，如果某个黑点被删除了，则跳过它的右邻居，处理下一点。对矩形这样做完一遍，水平方向会减少两像素。然后我们再改垂直方向扫描，方法一样。这样做一次水平扫描和垂直扫描，原图会“瘦”一圈。多次重复上面的步骤，直到图形不在变化为止。

实验代码：

def Thin(image, array):h, w = image.shape# height 图像的高度 weight 图像的宽度iThin = imagefor i in range(h):for j in range(w):if image[i, j] == 0:  # 如果图像的这个点是黑色的点a = [1] * 9  # 先默认周围八个点都是白色的for k in range(3):   # range(3):0 1 2for l in range(3):# 如果3*3矩阵的点不在边界# 且是黑色的点# 这一段的难点是用k，l来精确的表示3*3的矩阵if -1 < (i - 1 + k) < h and -1 < (j - 1 + l) < w and iThin[i - 1 + k, j - 1 + l] == 0:a[k + 3 * l] = 0# sum没有a[4]，a[4]刚好就是中心的黑点，不需要加进去计算sum = a[0] * 1 + a[1] * 2 + a[2] * 4 + a[3] * 8 + a[5] * 16 + a[6] * 32 + a[7] * 64 + a[8] * 128# 然后根据array表，对iThin的那一点进行赋值。# 将所有点的总价值映射到[0,255]的表中，返回的结果保存到iThin中iThin[i, j] = array[sum] * 255return iThin

映射表：

array = [0, 0, 1, 1, 0, 0, 1, 1, 1, 1, 0, 1, 1, 1, 0, 1, \1, 1, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 1, \0, 0, 1, 1, 0, 0, 1, 1, 1, 1, 0, 1, 1, 1, 0, 1, \1, 1, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 1, \1, 1, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, \0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, \1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 1, 1, 1, 0, 1, \0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, \0, 0, 1, 1, 0, 0, 1, 1, 1, 1, 0, 1, 1, 1, 0, 1, \1, 1, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 1, \0, 0, 1, 1, 0, 0, 1, 1, 1, 1, 0, 1, 1, 1, 0, 1, \1, 1, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, \1, 1, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, \1, 1, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, \1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 1, 1, 1, 0, 0, \1, 1, 0, 0, 1, 1, 1, 0, 1, 1, 0, 0, 1, 0, 0, 0]

处理效果：
对 效果分析： 细化算法似乎有些问题，对物体骨架的提取不理想，对深色物体骨架的提取效果较好，对浅色物体骨架的提取有待改进。