图像数学形态学的基本原理与代码实现（腐蚀、膨胀、开闭运算）

图像数学形态学的基本原理与代码实现（腐蚀、膨胀、开闭运算）

1. 集合基础

集合是由一个或多个确定元素所构成的整体。

• 空集：确定元素如果不存在，就是空集。常记为 $\varnothing$
• 并集：所有元素合并到一起的集合。常记为 $A\cup B$
• 交集：相同元素的集合。常记为 $A\cap B$
• 补集：不属于该集合元素组成的集合。常记为 $A^{\mathrm{c}}$

• 包含：A 的元素全部在 B 中；$A\cap B=A$；$A\subseteq B$

• 相离：A 与 B 没有共同元素；$A\cap B=\varnothing$；

• 相交：A 与 B 相同元素；$A\cap B\ne \varnothing$

• 映像：所有元素取反组成的集合；记为 A ^ = { x ∣ x = − a , a ∈ A } \hat A = \{x| x=-a,a\in A\} A^={x∣x=−a,a∈A}

• 平移：所有元素平移 x，记为 ( A ) x = { y ∣ y = a + x , a ∈ A } (A)_x = \{y|y = a+x,a\in A\} (A)x​={y∣y=a+x,a∈A}

2. 腐蚀（erode）

2.1 集合描述

使用相关集合知识定义为：

A, B 为 $Z^2$ 的集合，A 被 B腐蚀，记为 $A\ominus B$，则腐蚀定义为

A ⊖ B = { x ∣ ( B ) x ⊆ A } A\ominus B = \{x\mid (B)_x \subseteq A \} A⊖B={x∣(B)x​⊆A}

2.2 通俗理解

过程

B 通常称为核，A 假设为 10x10 像素的二值图中1区域

先腐蚀 A 区域，过程如下：

对于核 B ，中心点坐标即为原点(0,0)，对于 A，左上角第一个元素坐标为 (1,1)，像素值记 $I(1,1)=0$。

不考虑 padding ，那么核 B 在 A 上遍历，即进行 B 的平移操作，那么第一个平移坐标 x 即为 (2,2)，对应图像如下

根据集合定义 A ⊖ B = { x ∣ ( B ) x ⊆ A } A\ominus B = \{x\mid (B)_x \subseteq A \} A⊖B={x∣(B)x​⊆A} ，对于此时 x = (2,2)，$(B{)}_x$即为淡红色区域，明显不包含于 A，因此(2,2)处值仍为后景像素值0,即 $I(2,2)=0$

按照此规则运算，直至到 x=(4,3)，即达到 A 区域内的像素，此时如下

同样按照定义，此时 $(B{)}_x$ 与 A 有交集，但仍旧不完全包含于 A，因此(4,3)处值为后景像素值0，即$I(4,3)$ 由 1 变为 0，就是腐蚀操作。

未被腐蚀的第一个 A 元素坐标为 (5,4)，此时如下

此时满足 $(B{)}_x\subseteq A$ ，$I(5,4)$ 仍为前景像素值 1.

按照此操作运算，可以得到最终结果

总结

遍历过程中，只有 B 完全在 A 内部，中心点所在坐标仍保留为 前景值1，其他全部为后景像素值 0

3. 膨胀（dilate）

3.1 集合描述

使用相关集合知识的定义为：

A, B 为 $Z^2$ 的集合，A 被 B 膨胀，记为 $A\oplus B$，则膨胀定义为

A ⊕ B = { x ∣ ( B ^ ) x ∩ A ⊆ A } A\oplus B = \{x\mid (\hat B)_x\cap A \subseteq A \} A⊕B={x∣(B^)x​∩A⊆A}

3.2 通俗理解

过程

过程可与腐蚀过程类比理解，值得说明的是，因为 B 定义中心点坐标为 (0,0)，$\hat{B}$ 为元素原点对称后的集合，也即为B。

第一个膨胀点坐标是 (3,2)，此时如下

根据集合定义，$(\hat{B}{)}_x\cap A$ 即为 A 左上角第一个元素，此时满足 $(\hat{B}{)}_x\cap A\subseteq A$，则 $I(3,2)$ 的元素值由后景像素值 0 变为前景像素值 1，也即膨胀操作。

最后结果

总结

遍历移动过程中， B 与 A 只要有交集，则 B 的中心位置像素值设定为 前景值1 。

4. 代码实现

代码

腐蚀和膨胀操作理论上是等同的，只是同一操作下前景与后景的置换（上述例子中，对 1 区域的膨胀就等价与对 0 区域的腐蚀），因此代码可以归到一起。

使用python实现二维下的腐蚀和膨胀操作：

import numpy as np
import cv2

def morphology_2d(img, kenel, value = 1, ispadding = True):
	kenelx = int(kenel.shape[0]/2)
	kenely = int(kenel.shape[1]/2)
	tempImg = img.copy()
	if ispadding:
		imgPadding = np.pad(img, ((kenelx, kenelx), (kenely,kenely)), mode='edge')
	for x in range(img.shape[0]):
		for y in range(img.shape[1]):
			b = 0
			for i in range(-kenelx, kenelx+1):
				if b:
					break
				for j in range(-kenely, kenely+1):
					if kenel[kenelx+i, kenely+j] == 1 and imgPadding[x+i, y+j] == value:
						tempImg[x, y] = value
						b = 1
						break
	return tempImg

def erode_2d(img, kenel):
	return morphology_2d(img, kenel, value = 0)

def dilate_2d(img, kenel):
	return morphology_2d(img, kenel, value = 255)

if __name__ == '__main__':
	img = cv2.imread('./plus.jpg')
	gray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)
	ret, img = cv2.threshold(gray, 0, 255, cv2.THRESH_BINARY+cv2.THRESH_OTSU)
	img = np.uint8(img)
	kernel = np.ones((13,13))
	# 自编程
	erode = erode_2d(img, kernel)
	dilate = dilate_2d(img, kernel)
	# opencv库函数
	erosion = cv2.erode(img,kernel,iterations = 1)
	dilation = cv2.dilate(img, kernel, iterations = 1)
    
	cv2.imshow('img', img)	# 原图
	cv2.imshow('erode', erode)	# 自编程腐蚀
	cv2.imshow('erosion', erosion)	# cv库腐蚀
	cv2.imshow('dilate', dilate)	# 膨胀
	cv2.imshow('dilation', dilation)	# cv库膨胀
	cv2.waitKey()

结果

腐蚀

膨胀

说明

• 注意前景与后景像素值，即是前白后黑，还是前黑后白。

• 使用了四层循环，时间复杂度达到$O(n^4)$ ，计算效率较慢。opencv官方库计算较快。也可以使用卷积操作编程。待探究其原因。

• 四层循环中，break 退出两层循环

5.其他

• 三维数据的腐蚀膨胀实现同二维情况下原理相同，只是核变为三维

参考代码

def morphology_3d(imgMatrix, kenel, value = 1, ispadding = False):
	kenelz = int(kenel.shape[0]/2)
	kenelx = int(kenel.shape[1]/2)
	kenely = int(kenel.shape[2]/2)
	tempImg = imgMatrix.copy()
	for z in range(kenelz, imgMatrix.shape[0]-kenelz):
		for x in range(kenelx, imgMatrix.shape[1]-kenelx):
			for y in range(kenely, imgMatrix.shape[2]-kenely):
				b1 = 0
				b2 = 0
				for k in range(-kenelz, kenelz+1):
					if b1:
						break
					for i in range(-kenelx, kenelx+1):
						if b2:
							b1 = 1
							break
						for j in range(-kenely, kenely+1):
							if kenel[kenelz+k, kenelx+i, kenely+j] == 1 and tempImg[z+k, x+i, y+j] == value:
								tempImg[z, x, y] = value
								b2 = 1
								break
	return tempImg

因为$O(n^6)$ 的时间复杂度，代码运行奇慢

• 开运算是先腐蚀后膨胀，去除孤立区，保持主体形状不变
• 闭运算是先膨胀后腐蚀，去除孔洞，保持主体形状不变