中轴线算法
这是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 }