骨架算法(skeleton)的实现_skeleton3d函数图示

一、算法原理

具体的算法原理，参考的是文章A Fast Parallel Algorithm for Thinning Digital Patterns，这也是opencv中骨架提取算法的实现方式。

下面开始造轮子。

1.文章整体的思路

可以看到上面的为图像细化的过程，其中前两幅图像为两次子迭代过程， 第三幅图像是最终结果。

• 第一幅图为第一个子迭代过程，目的是通过一定的规则去除图像中左边和下边的多余的点
• 第二幅图为第二个子迭代过程，目的是通过一定的规则去除图像中右边和上边的多余的点
• 第三幅图为最终结果，实际上要经过多次的子迭代才能得到最终的结果

2. 判断像素点是否要被删除的规则

(1)文章中的第三小节展示了如何判断像素点是否被删除的条件

• 第一次子迭代(a)(b)(c )(d)
• 第二次子迭代(a)(b)(c’)(d’)

(2)条件的同等变换（转化成几何的形式更好理解）

文章中的规则为数学表示，而且比较绕，我同等变换为以下三个条件，如果一个像素需要被删除，那么必须满足以下全部条件：

a.第一次子迭代

1. 像素P的八个相邻像素中非零值靠在一起和零值靠在一起，如下图所示：

2. 像素P的八个相邻的像素中非零值像素数量为2 ~ 6
3. 像素P的邻居4号和6号至少有一个为0，或者邻居2号和8号全部为0(邻居像素的序号如下图所示)

b.第二次子迭代

1. 像素P的八个相邻像素中非零值靠在一起和零值靠在一起，如下图所示：
2. 像素P的八个相邻的像素中非零值像素数量为2 ~ 6
3. 像素P的邻居2号和8号至少有一个为0，或者邻居4号和6号全部为0

3. 数据结构的设计

现在要设计一个数据结构来表示一个点位是否满足被删除。

(1) 表示像素P的八个邻居(uchar)surround

很简单可以想到一个字节的长度正好是8个bit，可以完美表示像素P的八个邻居的状态（非零为1或者为零为0），按照bit最高位对应像素序号为1的邻居，最低位对应序号为8的邻居。

(2) 第3个条件的表示

上面提到的条件1和2表示起来比较方便，重点是第3个条件的表示。

• 假设邻居P的邻居表示为surround
• 定义mask_east为0x10、mask_south为0x04、mask_west为0x01、mask_north为0x40
• 在以上定义的情况下2、4、6、8邻居分别为0可表示为：

综上所述，条件3可表示为：(B || C) || (A && D)。
我们使用他的逆否条件（方便编写程序）就可以表示为：!(!B && !C) && !(!A || !D)。

(3)定义像素属性的描述子Descriptor

/**
 * this struct help to describe the pixel attributes
 * @param surround is a unsigned char which has eight bits and are corresponding to a pixel's eight neighbours
 * @param pattern_cnt is the number of pattern where the adjacent two neighbours are 0 and none-zero in clock-order
 * @param neighbour_cnt is the number of none-zero neighbour
 * @param delete_state whether this pixel has been deleted 
 * @param iter_idx if this pixel has been deleted, the variable represent in which iteration this happened
 */
typedef struct _Descriptor_
{
	uchar surround = 0;
	int pattern_cnt = 0;
	int neighbour_cnt = 0;
	bool delete_state = false;
	int iter_idx = -1;
}Descriptor;
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16

• surround表示邻居像素的组成结构
• pattern_cnt对应条件1，满足条件1时的值为1，大于1的时候不满足
• neighbour_cnt对应条件2，记录非零邻居的数量
• delete_state代表在迭代的过程中像素点是否被删除，未被删除值为false
• iter_idx代表在哪一次迭代过程中该像素点被删除，未被删除时值为-1

(4)编写程序中需要注意的点

• 由于该算法是迭代运行的，每次迭代后都会使得部分像素点被删除（骨架化），在单次迭代的过程中删除点位会影响该次迭代过程中剩余步骤的计算。简单的解决办法是记录这次迭代中要删除的点位，这次迭代完成之后在集中进行删除，缺点是会浪费储存空间和指针遍历的时间；我在实现的过程中采用的是使用iter_idx 来进行相应的记录，进而消除 影响
• 对于处理过但是没有删除的点位的Descriptor要进行重新初始化，否则会出现错误。
• 在没有发现有新增删除点位时退出程序。
• 实际上两次子迭代过程也要按次序执行，使用iter_idx来进行前后顺序的区分，这影响了算法执行的效率

二、 算法效果

白色为原图，灰色部分为提取出来的骨架。

三、代码实现

函数实现

// Where you can find this article: A Fast Parallel Algorithm for Thinning Digital Patterns 
// <https://dl.acm.org/doi/10.1145/357994.358023>.


/**
 * this struct help to describe the pixel attributes
 * @param surround is a unsigned char which has eight bits and are corresponding to a pixel's eight neighbours
 * @param pattern_cnt is the number of pattern where the adjacent two neighbours are 0 and none-zero in clock-order
 * @param neighbour_cnt is the number of none-zero neighbour
 * @param delete_state whether this pixel has been deleted 
 * @param iter_idx if this pixel has been deleted, the variable represent in which iteration this happened
 */
typedef struct _Descriptor_
{
	uchar surround = 0;
	int pattern_cnt = 0;
	int neighbour_cnt = 0;
	bool delete_state = false;
	int iter_idx = -1;
}Descriptor;

/// @brief this function is used to find the skeleton of patterns in an image
/// @param srcimage IO- input and output image
void skeletonTry(cv::Mat& srcimage)
{
	int dx[8] = { -1, 0, 1, 1, 1, 0, -1, -1 };
	int dy[8] = { -1, -1, -1, 0, 1, 1, 1, 0 };
	// mask: 1 is the highest bit,
	// 1 2 3     mask      0    north   0
	// 8 P 4  --------->  west    p    east
	// 7 6 5               0    south   0
	uchar mask_east = 1 << 4;
	uchar mask_south = 1 << 2;
	uchar mask_west = 1;
	uchar mask_north = 1 << 6;
	size_t cols = srcimage.cols;
	Descriptor* describ = new Descriptor[cols * (size_t)(srcimage.rows)];
	int cnt_break = 2;
	int iter_cnt = 0;
	bool del_found = false;
	while (true) {
		del_found = false;
		for (int i = 1; i < srcimage.rows - 1; i++) { // test point
			uchar* data_a[3];
			data_a[0] = srcimage.ptr<uchar>(i - 1);
			data_a[1] = srcimage.ptr<uchar>(i);
			data_a[2] = srcimage.ptr<uchar>(i + 1);
			for (int j = 1; j < cols - 1; j++) {
				// 1 2 3 order
				// 8 P 4
				// 7 6 5
				uchar p = *(data_a[1] + j);
				Descriptor* tmp_describ = describ + i * cols + j;
				if (0 != p && tmp_describ->delete_state != true) { // the pixel is none zero
					int surround_flag = 0; // auxiliary flag which is used to calculate pattern_cnt
					for (int z = 0; z < 7; z++) {
						Descriptor* surround_describ = describ + (i + dy[z]) * cols + j + dx[z];
						if (0 < *(data_a[dy[z] + 1] + dx[z] + j) || surround_describ->iter_idx == iter_cnt) { // this neighbour count, surround_flag will be setted as 1
							tmp_describ->surround += 1;   // the neighbour composition 
							tmp_describ->neighbour_cnt++; // counting the neighbour
							if (1 != surround_flag) { //  01 pattern found(initial is 0, so if the successor is 1 the counter will self-add by 1)
								tmp_describ->pattern_cnt++;
								if (tmp_describ->pattern_cnt > cnt_break) {
									break;
								}
							}
							surround_flag = 1;
						}
						else {
							surround_flag = 0;
						}
						tmp_describ->surround = tmp_describ->surround << 1;
					}
					Descriptor* surround_describ = describ + (i + dy[7]) * cols + j + dx[7];
					if ((0 < *(data_a[dy[7] + 1] + dx[7] + j) || surround_describ->iter_idx == iter_cnt) && tmp_describ->pattern_cnt <= cnt_break) {
						tmp_describ->surround += 1;
						tmp_describ->neighbour_cnt++; 
						if (1 != surround_flag) { // the second condition is used to avoid repetitive count of the first place
							tmp_describ->pattern_cnt++;
						}
						if (tmp_describ->surround & 0x80) {
							tmp_describ->pattern_cnt--;
						}
					}
					// only if the specific pattern and neighbour number is satisfied
					// then check the (c) and (d) condition in the paper
					if (tmp_describ->pattern_cnt == 1 && tmp_describ->neighbour_cnt > 1 && tmp_describ->neighbour_cnt < 7) {
						// delete east north points and corner points of west sourth 
						if (!(
								(tmp_describ->surround & mask_east && tmp_describ->surround & mask_south) &&
								(tmp_describ->surround & mask_west || tmp_describ->surround & mask_north)
							) && iter_cnt % 2) { // this condition is derived from {!(s & m_e) || !(s & m_s)} || {!(s & m_w) && !(s & m_n)} by applying converse-negative proposition
							tmp_describ->delete_state = true;
							tmp_describ->iter_idx = iter_cnt;
							del_found = true;
							*(data_a[1] + j) = 0;
						}
						// delete west south points and corner points of east north(ok)
						if (!(
							(tmp_describ->surround & mask_west && tmp_describ->surround & mask_north) &&
							(tmp_describ->surround & mask_east || tmp_describ->surround & mask_south)
							) && (iter_cnt + 1) % 2) { // this condition is derived from {!(s & m_w) || !(s & m_n)} || {!(s & m_e) && !(s & m_s)} by applying converse-negative proposition
							tmp_describ->delete_state = true;
							tmp_describ->iter_idx = iter_cnt;
							del_found = true;
							*(data_a[1] + j) = 0;
						}
					}
				}
				if (!tmp_describ->delete_state) { // reset the descriptor
					tmp_describ->neighbour_cnt = 0;
					tmp_describ->pattern_cnt = 0;
					tmp_describ->surround = 0;
				}
			}
		}
		iter_cnt++;
		if (!del_found) {
			break;
		}
	}
}
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122

调用

int main(int argc, char* argv)
{
	//callEdgeClose();
	cv::Mat test(cv::Size(112, 112), CV_8UC1, cv::Scalar(0));
	cv::Mat region = test(cv::Rect(1, 1, 100, 100));
	cv::Mat region2 = test(cv::Rect(10, 10, 100, 100));
	cv::Mat region1 = test(cv::Rect(10, 10, 80, 80));
	cv::Mat region3 = test(cv::Rect(90, 90, 5, 5));
	cv::Mat region4 = test(cv::Rect(40, 40, 15, 66));
	cv::Mat region5 = test(cv::Rect(40, 40, 5, 25));
	region = cv::Scalar(255);
	region2 = cv::Scalar(255);
	region1 = cv::Scalar(0);
	region3 = cv::Scalar(0);
	region4 = cv::Scalar(0);
	region5 = cv::Scalar(255);
	cv::circle(test, cv::Point(55, 55), 40, cv::Scalar(255), 2);
	//cv::rectangle(test, cv::Rect(1, 1, 100, 100), cv::Scalar(255), 10);
	std::string im_name = "origin";
	cv::namedWindow(im_name, 0);
	cv::imshow(im_name, test);
	cv::waitKey();
	cv::Mat skeleton_im = test.clone();
	skeletonTry(skeleton_im);
	im_name = "skeleton";
	cv::namedWindow(im_name, 0);
	cv::imshow(im_name, test - skeleton_im / 2);
	cv::waitKey();
	return 0;
}
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30

四、代码改进（TODO）

• 算法的每次迭代都是全图搜索，其实这是不必要的：
  因为第n+1次迭代过程中所经过的闭合轮廓必定和第n次的迭代有所交集，因此我们只需要在第一次迭代过程中记录闭合轮廓到一个轮廓队列中即可，在后面的迭代过程中每次按照记录的轮廓进行搜索，而非全图搜索，并且在新的迭代过程中更新轮廓队列。
• 涉及的两次子迭代过程也存在着无效的全图迭代，完成第一条优化的基础上可以将两次子迭代的过程合并成一个过程。