中轴线算法
这是http://www.imagepy.org/的作者原创,我只是对其理解之后改进和说明,欢迎大家使用这个小软件!
如果有需要C++版本的朋友,可以博文结尾留邮箱!
首先上效果图:
算法的流程:
第一步:
距离变换
第二步:
把距离变换的图像进行像素值的排列(排列返回像素的位置信息)
第二步:
从小到大进行像素的查表操作
注释:这里为什么叫中轴线提取?因为提取的过程是绝对的按照对称来的,距离变换的结果就是前景到背景的距离,所以结果是绝对的集合中心。
Python代码:
1 import numpy as np
2 from skimage.data import horse, camera
3 import matplotlib.pyplot as plt
4 import scipy.ndimage as ndimg
5 from numba import jit
6 import cv2
7 from scipy.ndimage import label, generate_binary_structure
8
9 strc = np.ones((3, 3), dtype=np.bool)
10
11
12 # check whether this pixcel can be removed
13 def check(n):
14 a = [(n >> i) & 1 for i in range(8)]
15 a.insert(4, 0) # make the 3x3 unit
16 # if up, down, left, right all are 1, you cannot make a hole
17 # if a[1] & a[3] & a[5] & a[7]:return False
18 a = np.array(a).reshape((3, 3))
19 # segments
20 n = label(a, strc)[1]
21 # if sum is 0, it is a isolate point, you cannot remove it.
22 # if number of segments > 2, you cannot split them.
23 return n < 2
24 return a.sum() > 1 and n < 2
25 if a.sum() == 1 or n > 2: return 2
26 if a.sum() > 1 and n < 2: return 1
27 return 0
28
29
30 lut = np.array([check(n) for n in range(256)])
31 lut = np.dot(lut.reshape((-1, 8)), [1, 2, 4, 8, 16, 32, 64, 128]).astype(np.uint8)
32 '''
33 lut = np.array([200, 206, 220, 204, 0, 207, 0, 204, 0, 207, 221, 51, 1, 207, 221, 51,
34 0, 0, 221, 204, 0, 0, 0, 204, 1, 207, 221, 51, 1, 207, 221, 51], dtype=np.int8)
35 '''
36
37
38 @jit
39 def skel2dp(data, idx, lup):
40 h, w = data.shape
41 data = data.ravel()
42 for id in idx:
43
44 if data[id] == 0: continue
45 i2 = id - w
46 i8 = id + w
47 i1 = i2 - 1
48 i3 = i2 + 1
49 i4 = id - 1
50 i6 = id + 1
51 i7 = i8 - 1
52 i9 = i8 + 1
53 c = (data[i1] > 0) << 0 | (data[i2] > 0) << 1 \
54 | (data[i3] > 0) << 2 | (data[i4] > 0) << 3 \
55 | (data[i6] > 0) << 4 | (data[i7] > 0) << 5 \
56 | (data[i8] > 0) << 6 | (data[i9] > 0) << 7
57 if (lup[c // 8] >> c % 8) & 1: data[id] = 0
58 return 0
59
60
61 def mid_axis(img):
62 dis = ndimg.distance_transform_edt(img)
63 idx = np.argsort(dis.flat).astype(np.int32)
64 skel2dp(dis, idx, lut)
65 return dis
66
67
68 from time import time
69 img = ~horse()*255
70 #img = cv2.imread('123.jpg')
71 #img = cv2.cvtColor(img,cv2.COLOR_BGR2GRAY)
72 #ret2, img = cv2.threshold(img, 0, 255, cv2.THRESH_BINARY | cv2.THRESH_OTSU)
73 dis = ndimg.distance_transform_edt(img)
74 plt.imshow(dis)
75 idx = np.argsort(dis.flat).astype(np.int32)
76 a = skel2dp(dis, idx, lut)
77 #mid_axis(img.copy())
78 t1 = time()
79 a = mid_axis(img)
80 t2 = time()
81 print(t2 - t1)
82 plt.imshow(a)
83 plt.show()
C++代码:
1 void center_axis(InputArray _src, Mat& dst)
2 {
3 typedef struct MyStruct
4 {
5 Point position;
6 float num;
7 }MyStruct;
8 int wjy_array[] = { 0, 0, 1, 1, 0, 0, 1, 1, 1, 1, 0, 1, 1, 1, 0, 1,
9 1, 1, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 1,
10 0, 0, 1, 1, 0, 0, 1, 1, 1, 1, 0, 1, 1, 1, 0, 1,
11 1, 1, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 1,
12 1, 1, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
13 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
14 1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 1, 1, 1, 0, 1,
15 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
16 0, 0, 1, 1, 0, 0, 1, 1, 1, 1, 0, 1, 1, 1, 0, 1,
17 1, 1, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 1,
18 0, 0, 1, 1, 0, 0, 1, 1, 1, 1, 0, 1, 1, 1, 0, 1,
19 1, 1, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0,
20 1, 1, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
21 1, 1, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0,
22 1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 1, 1, 1, 0, 0,
23 1, 1, 0, 0, 1, 1, 1, 0, 1, 1, 0, 0, 1, 0, 0, 0 };
24 uchar lut[] = { 200, 206, 220, 204, 0, 207, 0, 204, 0, 207, 221, 51, 1, 207, 221, 51,
25 0, 0, 221, 204, 0, 0, 0, 204, 1, 207, 221, 51, 1, 207, 221, 51 };
26 Mat src = _src.getMat();
27 //Mat dst = _dst.getMat();
28 distanceTransform(src, src, DIST_L2, DIST_MASK_3, 5);
29 normalize(src, src, 0, 255, NORM_MINMAX);
30 Mat img_row = src.reshape(0, 1);
31 vector<MyStruct> my_vector;
32 for (size_t j = 0; j < img_row.cols; j++)
33 {
34 if (img_row.at<float>(0, j) == 0) continue;
35 MyStruct my_struct;
36 my_struct.num = saturate_cast<float>(img_row.at<float>(0, j));
37 my_struct.position = Point(saturate_cast<int>(j % src.cols), saturate_cast<int>(j / src.cols));
38 my_vector.push_back(my_struct);
39 }
40 for (size_t i = 0; i < my_vector.size(); i++)
41 {
42 if (my_vector[i].num == 0) continue;
43 for (size_t j = i; j < my_vector.size(); j++)
44 {
45 MyStruct temp;
46 if (my_vector[i].num >= my_vector[j].num)
47 {
48 if (my_vector[j].num == 0) continue;
49 temp = my_vector[j];
50 my_vector[j] = my_vector[i];
51 my_vector[i] = temp;
52 }
53 }
54 }
55 for (size_t i = 0; i < my_vector.size(); i++)
56 {
57 if (my_vector[i].position.y == 1
58 || my_vector[i].position.x == 1
59 || my_vector[i].position.y == src.rows - 1
60 || my_vector[i].position.x == src.cols - 1
61 || src.at<float>(my_vector[i].position.y, my_vector[i].position.x) == 0) continue;
62 else
63 {
64 char num[] = { 1,1,1,1,1,1,1,1 };
65 num[0] = src.at<float>(my_vector[i].position.y - 1, my_vector[i].position.x - 1)
66 > 0 ? 1 : 0;
67 num[1] = src.at<float>(my_vector[i].position.y - 1, my_vector[i].position.x)
68 > 0 ? 1 : 0;
69 num[2] = src.at<float>(my_vector[i].position.y - 1, my_vector[i].position.x + 1)
70 > 0 ? 1 : 0;
71 num[3] = src.at<float>(my_vector[i].position.y, my_vector[i].position.x - 1)
72 > 0 ? 1 : 0;
73 num[4] = src.at<float>(my_vector[i].position.y, my_vector[i].position.x + 1)
74 > 0 ? 1 : 0;
75 num[5] = src.at<float>(my_vector[i].position.y + 1, my_vector[i].position.x - 1)
76 > 0 ? 1 : 0;
77 num[6] = src.at<float>(my_vector[i].position.y + 1, my_vector[i].position.x)
78 > 0 ? 1 : 0;
79 num[7] = src.at<float>(my_vector[i].position.y + 1, my_vector[i].position.x + 1)
80 > 0 ? 1 : 0;
81 int sum = num[0] + num[1] * 2 + num[2] * 4 + num[3] * 8
82 + num[4] * 16 + num[5] * 32 + num[6] * 64 + num[7] * 128;
83 src.at<float>(my_vector[i].position.y, my_vector[i].position.x) = ((lut[uchar(sum / 8)] >> sum % 8) & 1) != 1 ? 255 : 0;
84 }
85 }
86 dst = src.clone();
87 dst.convertTo(dst, CV_8UC1);
88 }