## 《图像处理实例》 之 中轴线提取

第一步：

距离变换

第二步：

把距离变换的图像进行像素值的排列（排列返回像素的位置信息）

第二步：

从小到大进行像素的查表操作

注释：这里为什么叫中轴线提取？因为提取的过程是绝对的按照对称来的，距离变换的结果就是前景到背景的距离，所以结果是绝对的集合中心。

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
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 }

posted on 2017-08-21 02:07  影醉阏轩窗  阅读(885)  评论(3编辑  收藏

### 导航

/* 线条鼠标集合 */ /* 鼠标点击求赞文字特效 */ //带头像评论