## 《图像处理实例》 之 拓扑重建

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__':
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.......

对标号点进行路径追溯操作

先说本程序使用的方法：

在寻找路径的过程中把路径的点记录下来。

防止寻找重复，就是在寻找之后把路径清除。

再说说我自己的想法：

我的想法比本程序的稍微简单一些，思路基本是相同的。就是直接以角点为基础向四周寻找角点。

把得到的角点和路径存放到一个“图”中，这样使用的时候就非常方便了

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)):
132     for s,e,pts in edges:
133         l = np.linalg.norm(pts[1:]-pts[:-1], axis=1).sum()
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__':
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    （基础操作）