《opencv实战》 之 中轴线提取

 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:          #idx是距离变换之后排序的位置信息
43 
44         if data[id] == 0: continue  #背景直接返回
45         #周围八个数据
46         i2 = id - w
47         i8 = id + w
48         i1 = i2 - 1
49         i3 = i2 + 1
50         i4 = id - 1
51         i6 = id + 1
52         i7 = i8 - 1
53         i9 = i8 + 1
54         #下面这段话就是查表，和图像细化的查的表格一样
55         c = (data[i1] > 0) << 0 | (data[i2] > 0) << 1 56             | (data[i3] > 0) << 2 | (data[i4] > 0) << 3 57             | (data[i6] > 0) << 4 | (data[i7] > 0) << 5 58             | (data[i8] > 0) << 6 | (data[i9] > 0) << 7
59         if (lup[c // 8] >> c % 8) & 1: data[id] = 0
60     return 0
61 
62 
63 def mid_axis(img):
64     dis = ndimg.distance_transform_edt(img)
65     idx = np.argsort(dis.flat).astype(np.int32)     #flat是降维操作（二维到一维度），argsort是排序之后的位置（不是数值）
66     skel2dp(dis, idx, lut)
67     return dis
68 
69 
70 from time import time
71 img = ~horse()*255
72 #img = cv2.imread(‘123.jpg‘)
73 #img = cv2.cvtColor(img,cv2.COLOR_BGR2GRAY)
74 #ret2, img = cv2.threshold(img, 0, 255, cv2.THRESH_BINARY | cv2.THRESH_OTSU)
75 dis = ndimg.distance_transform_edt(img)
76 plt.imshow(dis)
77 idx = np.argsort(dis.flat).astype(np.int32)
78 a = skel2dp(dis, idx, lut)
79 #mid_axis(img.copy())
80 t1 = time()
81 a = mid_axis(img)
82 t2 = time()
83 print(t2 - t1)
84 plt.imshow(a)
85 plt.show()