中轴线算法
这是的作者原创，我只是对其理解之后改进和说明，欢迎大家使用这个小软件！
如果有需要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 removed13 def check(n):14 a = [(n >> i) & 1 for i in range(8)]15 a.insert(4, 0) # make the 3x3 unit16 # if up, down, left, right all are 1, you cannot make a hole17 # if a[1] & a[3] & a[5] & a[7]:return False18 a = np.array(a).reshape((3, 3))19 # segments20 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 < 224 return a.sum() > 1 and n < 225 if a.sum() == 1 or n > 2: return 226 if a.sum() > 1 and n < 2: return 127 return 028 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 @jit39 def skel2dp(data, idx, lup):40 h, w = data.shape41 data = data.ravel() 42 for id in idx: 43 44 if data[id] == 0: continue 45 i2 = id - w46 i8 = id + w47 i1 = i2 - 148 i3 = i2 + 149 i4 = id - 150 i6 = id + 151 i7 = i8 - 152 i9 = i8 + 153 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) << 757 if (lup[c // 8] >> c % 8) & 1: data[id] = 058 return 059 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 dis66 67 68 from time import time69 img = ~horse()*25570 #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 my_vector;32 for (size_t j = 0; j < img_row.cols; j++)33 {34 if (img_row.at (0, j) == 0) continue;35 MyStruct my_struct;36 my_struct.num = saturate_cast (img_row.at (0, j));37 my_struct.position = Point(saturate_cast (j % src.cols), saturate_cast (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 == 158 || my_vector[i].position.x == 159 || my_vector[i].position.y == src.rows - 160 || my_vector[i].position.x == src.cols - 161 || src.at (my_vector[i].position.y, my_vector[i].position.x) == 0) continue;62 else63 {64 char num[] = { 1,1,1,1,1,1,1,1 };65 num[0] = src.at (my_vector[i].position.y - 1, my_vector[i].position.x - 1)66 > 0 ? 1 : 0;67 num[1] = src.at (my_vector[i].position.y - 1, my_vector[i].position.x)68 > 0 ? 1 : 0;69 num[2] = src.at (my_vector[i].position.y - 1, my_vector[i].position.x + 1)70 > 0 ? 1 : 0;71 num[3] = src.at (my_vector[i].position.y, my_vector[i].position.x - 1)72 > 0 ? 1 : 0;73 num[4] = src.at (my_vector[i].position.y, my_vector[i].position.x + 1)74 > 0 ? 1 : 0;75 num[5] = src.at (my_vector[i].position.y + 1, my_vector[i].position.x - 1)76 > 0 ? 1 : 0;77 num[6] = src.at (my_vector[i].position.y + 1, my_vector[i].position.x)78 > 0 ? 1 : 0;79 num[7] = src.at (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] * 882 + num[4] * 16 + num[5] * 32 + num[6] * 64 + num[7] * 128;83 src.at (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 }