拓撲重建
這是http://www.imagepy.org/的作者原創，我只是對其理解之后改進和說明，歡迎大家使用這個小軟件！
第一個用Python寫的小項目，不算自創的，大部分都是參考別人的，由於第一次用python寫opencv遇到了無數的問題，最后也算完成了。
opencv的入門我就不寫了，以前都學過差不多，在此記錄一下基礎！
基礎操作
首先說一下python條件下用opencv，重點記錄，不重點看文檔就好！
1.如何創建一個空白圖片
A.
1 test = np.array([[0,255,0],[0,0,0],[0,0,0]])
2 np.array([[0,255,0],[0,0,0],[0,0,0]],np.uint8)#千萬別寫成這種格式，他把格式也存進去了
B.
1 showImage = np.zeros((inputImage.shape[0],inputImage.shape[1],3),np.uint8)
C.
1 showImage = inputImage.copy()
2 showImage = cv2.cvtColor(showImage,cv2.COLOR_GRAY2BGR)
2.如何讀寫一個圖片
A.
1 img[j,i] = 255#單通道讀寫
B.
1 #三通道讀寫
2 img[j,i,0]= 255
3 img[j,i,1]= 255
4 img[j,i,2]= 255
C.
1 img.itemset((10,10,2),100)#寫入
2 img.item(10,10,2)#讀
3.畫圖形
1 cv2.circle(showImage,(j+1,i+1), 5, (0,0,255), -1)#畫圓注意圓心是X,Y
4.python大部分函數都有返回值
#python的函數貌似都有返回值 img = cv2.cvtColor(img,cv2.COLOR_BGR2GRAY) ret2,img = cv2.threshold(img,0,255,cv2.THRESH_BINARY | cv2.THRESH_OTSU)
5.
1 Shape of image is accessed by img.shape. It returns a tuple of number of rows, columns and channels
1 Total number of pixels is accessed by img.size:
1 Image datatype is obtained by img.dtype:
骨架提取
廢話說了不少開始實戰了：
在此說明：
本程序和算法主要參考：http://www.cnblogs.com/xianglan/archive/2011/01/01/1923779.html寫的非常好，大家初學直接看大佬的博客就行了。
本文后續拓普重建主要是：https://github.com/yxdragon/sknw沒有看這位大佬的程序，但是他的思想幫助我很多。
因為得到別人幫助才有我自己的收獲，所以樂於分享其他初學者！
以下是和大佬的對話，沒經過本人同意所以打了馬，非常樂於助人的大佬！
1 import cv2 2 import numpy as np 3 from matplotlib import pyplot as plt 4 5 array = [0,0,1,1,0,0,1,1,1,1,0,1,1,1,0,1,\ 6 1,1,0,0,1,1,1,1,0,0,0,0,0,0,0,1,\ 7 0,0,1,1,0,0,1,1,1,1,0,1,1,1,0,1,\ 8 1,1,0,0,1,1,1,1,0,0,0,0,0,0,0,1,\ 9 1,1,0,0,1,1,0,0,0,0,0,0,0,0,0,0,\ 10 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,\ 11 1,1,0,0,1,1,0,0,1,1,0,1,1,1,0,1,\ 12 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,\ 13 0,0,1,1,0,0,1,1,1,1,0,1,1,1,0,1,\ 14 1,1,0,0,1,1,1,1,0,0,0,0,0,0,0,1,\ 15 0,0,1,1,0,0,1,1,1,1,0,1,1,1,0,1,\ 16 1,1,0,0,1,1,1,1,0,0,0,0,0,0,0,0,\ 17 1,1,0,0,1,1,0,0,0,0,0,0,0,0,0,0,\ 18 1,1,0,0,1,1,1,1,0,0,0,0,0,0,0,0,\ 19 1,1,0,0,1,1,0,0,1,1,0,1,1,1,0,0,\ 20 1,1,0,0,1,1,1,0,1,1,0,0,1,0,0,0] 21 22 def Thin(image,array): 23 '''未改進算法''' 24 #midImage = image.copy() 25 for i in range(image.shape[0]-2): 26 for j in range(image.shape[1]-2): 27 if image[i+1,j+1] == 0: 28 a = [1]*9 #定義list[9] 29 for k in range(3): 30 for l in range(3): 31 a[k*3+l] = 0 if image[i+k,j+l]==0 else 1 32 sum = a[0]*1+a[1]*2+a[2]*4+a[3]*8+a[5]*16+a[6]*32+a[7]*64+a[8]*128 33 image[i+1,j+1] = array[sum]*255 34 return image 35 36 def HThin(image,array): 37 flag = True #如果該點被刪除，跳過下一個點 38 midImage = image.copy() 39 for i in range(image.shape[0]-2): 40 for j in range(image.shape[1]-2): 41 if flag == False: 42 flag == True 43 else: 44 if image[i+1,j+1] == 0 and (image[i+1,j] != 0 or image[i+1,j+2] != 0):#左右都為黑點不處理 45 a = [1]*9 #定義list[9] 46 for k in range(3): 47 for l in range(3): 48 a[k*3+l] = 0 if midImage[i+k,j+l]==0 else 1 49 sum = a[0]*1+a[1]*2+a[2]*4+a[3]*8+a[5]*16+a[6]*32+a[7]*64+a[8]*128 50 midImage[i+1,j+1] = array[sum]*255 51 return midImage 52 def VThin(image,array): 53 flag = True #如果該點被刪除，跳過下一個點 54 midImage = image.copy() 55 for i in range(image.shape[1]-2): 56 for j in range(image.shape[0]-2): 57 if flag == False: 58 flag == True 59 else: 60 if image[j+1,i+1] == 0 and (image[j,i+1] != 0 or image[j+2,i+1] != 0):#左右都為黑點不處理 61 a = [1]*9 #定義list[9] 62 for k in range(3): 63 for l in range(3): 64 a[k*3+l] = 0 if midImage[j+k,i+l]==0 else 1 65 sum = a[0]*1+a[1]*2+a[2]*4+a[3]*8+a[5]*16+a[6]*32+a[7]*64+a[8]*128 66 midImage[j+1,i+1] = array[sum]*255 67 return midImage 68 69 def wjy_Bone(inputImage,num=100): 70 '''改進算法''' 71 for i in range(num): 72 inputImage = VThin(HThin(inputImage,array),array) 73 return inputImage 74 75 def ThredImage(image,thred): 76 '''二值化圖像''' 77 imageGray = cv2.cvtColor(image,cv2.COLOR_BGR2GRAY) 78 midimage = np.zeros((imageGray.shape[0],imageGray.shape[1]),imageGray.dtype) 79 for i in range(imageGray.shape[0]): 80 for j in range(imageGray.shape[1]): 81 midimage[i,j] = 0 if imageGray[i,j] < int(thred) else 255 82 return midimage 83 84 def drawBone(inputImage): 85 #showImage = np.zeros((inputImage.shape[0],inputImage.shape[1],3),np.uint8) 86 showImage = inputImage.copy() 87 showImage = cv2.cvtColor(showImage,cv2.COLOR_GRAY2BGR) 88 for i in range(inputImage.shape[0]-2): 89 for j in range(inputImage.shape[1]-2): 90 if inputImage[i+1,j+1] == 255 : continue 91 flag = -1 92 for J in range(3): 93 for K in range(3): 94 if inputImage[i+J,j+K] == 0: 95 flag += 1 96 if(flag==1 or flag>=3): 97 showImage = cv2.circle(showImage,(j+1,i+1), 5, (0,0,255), -1) 98 showImage[i+1,j+1,0] = 255 99 showImage[i+1,j+1,1] = 0 100 showImage[i+1,j+1,2] = 0 101 return showImage 102 103 if __name__ == '__main__': 104 img = cv2.imread("5.jpg") 105 #test = np.array([[0,255,0],[0,0,0],[0,0,0]]) 106 #千萬寫寫成np.array([[0,255,0],[0,0,0],[0,0,0]],np.uint8) 107 img = cv2.cvtColor(img,cv2.COLOR_BGR2GRAY) 108 ret2,img = cv2.threshold(img,0,255,cv2.THRESH_BINARY | cv2.THRESH_OTSU) 109 #testImage = Thin(img,array) 110 #testImage = wjy_Bone(img) 111 testImage = drawBone(wjy_Bone(img,10)) 112 plt.imshow(testImage,cmap='gray',interpolation = 'bicubic') 113 plt.show() 114
輪廓追溯
第一步：
對骨架提取之后的圖片變換成（0,1）圖，骨架為0，背景為1
第二步：
提取骨架角點，方法在骨架提取的時候有說明。
第三步：
在提取角點的基礎上標記骨架圖，對角點標記為2，對輪廓標記為1，對背景標記為0
第四步：
對角點進行標記區分，10、11、12、13.......
第五步：（此處是對程序的另行說明，是整個程序核心點）
對標號點進行路徑追溯操作
先說本程序使用的方法：
先從角點進行周圍檢測，找邊緣
以找到的邊界點為中心向四周尋找邊界（角點）
在尋找路徑的過程中把路徑的點記錄下來。
防止尋找重復，就是在尋找之后把路徑清除。
找的結果就是10和11這兩個點
再說說我自己的想法：
我的想法比本程序的稍微簡單一些，思路基本是相同的。就是直接以角點為基礎向四周尋找角點。
第六步：
把得到的角點和路徑存放到一個“圖”中，這樣使用的時候就非常方便了
PS：圖我還沒有學過，這里不過多說明。
上代碼：
1 from cv2 import * 2 import numpy as np 3 import matplotlib.pyplot as plt 4 import networkx as nx 5 from numba import jit 6 7 @jit# get neighbors d index 8 #得到領域的坐標（如果是二維直接算可以，但如果是三維就不一樣了） 9 #這里是作者的改進適合多維圖像，每句話的含義我沒有弄懂 10 def neighbors(shape):#傳進來的是二值圖，所以維度是2 11 '''答案是：[-1921 -1920 -1919 -1 1 1919 1920 1921]''' 12 dim = len(shape) 13 block = np.ones([3]*dim) #維度為dim，數據3X3 14 block[tuple([1]*dim)] = 0 #取[1,1]位置為0 15 idx = np.where(block>0) #非零區域的位置信息 16 idx = np.array(idx, dtype=np.uint8).T #idx轉換成矩陣，並且轉置 17 idx = np.array(idx-[1]*dim) #全部值減一 18 acc = np.cumprod((1,)+shape[::-1][:-1]) #疊乘 19 return np.dot(idx, acc[::-1]) #點乘 20 21 @jit # my mark 22 #標記骨架圖每個點的含義 23 #0：不是骨架 24 #1：骨架路徑上的點 25 #2：端點或者角點 26 def mark(img): # mark the array use (0, 1, 2) 27 nbs = neighbors(img.shape) 28 img = img.ravel() 29 for p in range(len(img)): 30 if img[p]==0:continue 31 s = 0 32 '''判斷中心點領域八個點的情況''' 33 for dp in nbs: 34 if img[p+dp]!=0:s+=1 35 if s==2:img[p]=1 36 else:img[p]=2 37 38 39 @jit # trans index to r, c... 40 # 坐標轉換，圖像被轉換成一行處理，現在得把關鍵點轉換回來（x,y） 41 def idx2rc(idx, acc): 42 rst = np.zeros((len(idx), len(acc)), dtype=np.int16) 43 for i in range(len(idx)): 44 for j in range(len(acc)): 45 rst[i, j] = idx[i] // acc[j] 46 idx[i] -= rst[i, j] * acc[j] 47 return rst 48 49 50 @jit # fill a node (may be two or more points) 51 # 標記節點，把節點打上不同的編碼img.copy()=（2,2,2,2.....）----->>>>（10,11,12,13......） 52 # 並且把行坐標轉換成（x，y）類型的坐標存放 53 def fill(img, p, num, nbs, acc, buf): 54 back = img[p] # 位置點的值存放 55 img[p] = num # 計數值 56 buf[0] = p # 位置點存放 57 cur = 0; 58 s = 1; 59 while True: 60 p = buf[cur] # 61 for dp in nbs: 62 cp = p + dp 63 if img[cp] == back: 64 img[cp] = num 65 buf[s] = cp 66 s += 1 67 cur += 1 68 if cur == s: break 69 return idx2rc(buf[:s], acc) 70 71 72 @jit # trace the edge and use a buffer, then buf.copy, if use [] numba not works 73 #路徑跟蹤找邊緣 74 def trace(img, p, nbs, acc, buf): 75 '''這里的c1和c2是找的一條路徑的兩個端點，存放的端點標記>=10''' 76 c1 = 0 77 c2 = 0 78 newp = 0 #更新之后的位置 79 cur = 0 #計數，用來計算路徑上點的數量 80 81 b = p==97625 #b = (p==97625) 82 #while找的是一條路徑（具體說明請看附圖） 83 while True: 84 buf[cur] = p#位置存儲 85 img[p] = 0 #當前位置置0（搜索過得路徑置0，防止重復搜索） 86 cur += 1 #光標加一（移動一個位置） 87 for dp in nbs: 88 cp = p + dp 89 if img[cp] >= 10:#判斷是否為端點 90 if c1==0: c1=img[cp]#找的第一個端點（） 91 else: c2 = img[cp]# 92 if img[cp] == 1: 93 newp = cp 94 p = newp 95 if c2!=0:break 96 return (c1-10, c2-10, idx2rc(buf[:cur], acc)) 97 98 99 @jit # parse the image then get the nodes and edges 100 # 找節點和邊緣 101 def parse_struc(img): 102 nbs = neighbors(img.shape) 103 acc = np.cumprod((1,) + img.shape[::-1][:-1])[::-1] # （1，img.shape[0]） 104 img = img.ravel() 105 pts = np.array(np.where(img == 2))[0] 106 buf = np.zeros(20, dtype=np.int64) 107 num = 10 108 nodes = [] 109 for p in pts: 110 if img[p] == 2: 111 nds = fill(img, p, num, nbs, acc, buf) 112 num += 1 113 nodes.append(nds) 114 115 buf = np.zeros(10000, dtype=np.int64) 116 edges = [] 117 # 這個for循環是找一個點上對應的多個路徑 118 for p in pts: 119 for dp in nbs: 120 if img[p + dp] == 1: # 找路徑，img路徑點值為1 121 edge = trace(img, p + dp, nbs, acc, buf) 122 edges.append(edge) 123 return nodes, edges 124 125 126 # use nodes and edges build a networkx graph 127 #將建立的節點和路徑存放到NX的圖中 128 def build_graph(nodes, edges): 129 graph = nx.Graph() 130 for i in range(len(nodes)): 131 graph.add_node(i, pts=nodes[i], o=nodes[i].mean(axis=0)) 132 for s,e,pts in edges: 133 l = np.linalg.norm(pts[1:]-pts[:-1], axis=1).sum() 134 graph.add_edge(s,e, pts=pts, weight=l) 135 return graph 136 #建立一個圖 137 def build_sknw(ske): 138 mark(ske) 139 nodes, edges = parse_struc(ske.copy()) 140 return build_graph(nodes, edges) 141 #畫出一個'圖' 142 # draw the graph 143 def draw_graph(img, graph): 144 for idx in graph.nodes(): 145 pts = graph.node[idx]['pts'] 146 img[pts[:,0], pts[:,1]] = 255 147 for (s, e) in graph.edges(): 148 pts = graph[s][e]['pts'] 149 img[pts[:,0], pts[:,1]] = 128 150 151 array = [0,0,1,1,0,0,1,1,1,1,0,1,1,1,0,1, 152 1,1,0,0,1,1,1,1,0,0,0,0,0,0,0,1, 153 0,0,1,1,0,0,1,1,1,1,0,1,1,1,0,1, 154 1,1,0,0,1,1,1,1,0,0,0,0,0,0,0,1, 155 1,1,0,0,1,1,0,0,0,0,0,0,0,0,0,0, 156 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, 157 1,1,0,0,1,1,0,0,1,1,0,1,1,1,0,1, 158 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, 159 0,0,1,1,0,0,1,1,1,1,0,1,1,1,0,1, 160 1,1,0,0,1,1,1,1,0,0,0,0,0,0,0,1, 161 0,0,1,1,0,0,1,1,1,1,0,1,1,1,0,1, 162 1,1,0,0,1,1,1,1,0,0,0,0,0,0,0,0, 163 1,1,0,0,1,1,0,0,0,0,0,0,0,0,0,0, 164 1,1,0,0,1,1,1,1,0,0,0,0,0,0,0,0, 165 1,1,0,0,1,1,0,0,1,1,0,1,1,1,0,0, 166 1,1,0,0,1,1,1,0,1,1,0,0,1,0,0,0] 167 168 @jit#橫向細化 169 def HThin(image,array): 170 flag = True #如果該點被刪除，跳過下一個點 171 midImage = image.copy() 172 for i in range(image.shape[0]-2): 173 for j in range(image.shape[1]-2): 174 if flag == False: 175 flag = True 176 else: 177 if image[i+1,j+1] == 0 and (image[i+1,j] != 0 or image[i+1,j+2] != 0):#左右都為黑點不處理 178 a = [0]*9 #定義list[9] 179 for k in range(3): 180 for l in range(3): 181 a[k*3+l] = 0 if midImage[i+k,j+l]==0 else 1 182 sum = a[0]*1+a[1]*2+a[2]*4+a[3]*8+a[5]*16+a[6]*32+a[7]*64+a[8]*128 183 midImage[i+1,j+1] = array[sum]*255 184 return midImage 185 @jit#縱向細化 186 def VThin(image,array): 187 flag = True #如果該點被刪除，跳過下一個點 188 midImage = image.copy() 189 for i in range(image.shape[1]-2): 190 for j in range(image.shape[0]-2): 191 if flag == False: 192 flag = True 193 else: 194 if image[j+1,i+1] == 0 and (image[j,i+1] != 0 or image[j+2,i+1] != 0):#左右都為黑點不處理 195 a = [0]*9 #定義list[9] 196 for k in range(3): 197 for l in range(3): 198 a[k*3+l] = 0 if midImage[j+k,i+l]==0 else 1 199 sum = a[0]*1+a[1]*2+a[2]*4+a[3]*8+a[5]*16+a[6]*32+a[7]*64+a[8]*128 200 midImage[j+1,i+1] = array[sum]*255 201 return midImage 202 @jit#橫向和縱向合並 203 def wjy_Bone(inputImage,num=100): 204 '''改進算法''' 205 for i in range(num): 206 inputImage = VThin(HThin(inputImage,array),array) 207 return inputImage 208 @jit#骨架提取 209 def drawBone(inputImage): 210 #showImage = np.zeros((inputImage.shape[0],inputImage.shape[1],3),np.uint8) 211 showImage = inputImage.copy() 212 showImage = cvtColor(showImage,COLOR_GRAY2BGR) 213 for i in range(inputImage.shape[0]-2): 214 for j in range(inputImage.shape[1]-2): 215 if inputImage[i+1,j+1] > 0 : continue 216 flag = -1 217 for J in range(3): 218 for K in range(3): 219 if inputImage[i+J,j+K] == 0: 220 flag += 1 221 if(flag==1 or flag>=3): 222 showImage = circle(showImage,(j+1,i+1), 5, (0,0,255), -1) 223 showImage[i+1,j+1,0] = 255 224 showImage[i+1,j+1,1] = 0 225 showImage[i+1,j+1,2] = 0 226 return showImage 227 228 if __name__ == '__main__': 229 img = imread('5.jpg') 230 img = cvtColor(img,COLOR_BGR2GRAY) 231 ret2, img = threshold(img, 0, 255, THRESH_BINARY | THRESH_OTSU) 232 #wjy2 = drawBone(wjy_Bone(img, num = 500)) 233 img = wjy_Bone(img,num=500).astype(np.uint16)#細化之后從uint8轉換到uint16 234 img = np.where(img > 0, np.uint16(img/255),0) 235 img = np.where(img > 0, img - 1 , img + 1) 236 graph = build_sknw(img)#這里傳進來的圖像數值為（0 or 1） ，類型是uint16 237 238 plt.imshow(img, cmap='gray') 239 for (s, e) in graph.edges(): 240 ps = graph[s][e]['pts'] 241 plt.plot(ps[:, 1], ps[:, 0], 'green') 242 243 node, nodes = graph.node, graph.nodes() 244 ps = np.array([node[i]['o'] for i in nodes]) 245 plt.plot(ps[:, 1], ps[:, 0], 'r.') 246 plt.title('Build Graph') 247 plt.show()
不是太完美的圖，應該是細化出問題了
推廣一下大佬的開源項目：http://www.imagepy.org/開源的精神值得我們學習，問問題的時候盡量發個紅包，不在多少在心意。
參考：
http://blog.csdn.net/sunny2038/article/details/9080047 （基礎操作）
http://opencv-python-tutroals.readthedocs.io/en/latest/py_tutorials/py_core/py_basic_ops/py_basic_ops.html#basic-ops （英文文檔）