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

拓扑重建

这是http://www.imagepy.org/的作者原创，我只是对其理解之后改进和说明，欢迎大家使用这个小软件！

第一个用Python写的小项目，不算自创的，大部分都是参考别人的，由于第一次用python写opencv遇到了无数的问题，最后也算完成了。

opencv的入门我就不写了，以前都学过差不多，在此记录一下基础！

基础操作

首先说一下python条件下用opencv，重点记录，不重点看文档就好！

1.如何创建一个空白图片

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)#千万别写成这种格式，他把格式也存进去了

1 showImage = np.zeros((inputImage.shape[0],inputImage.shape[1],3),np.uint8)

1 showImage = inputImage.copy()

2 showImage = cv2.cvtColor(showImage,cv2.COLOR_GRAY2BGR)

2.如何读写一个图片

1 img[j,i] = 255#单通道读写

1 #三通道读写

2 img[j,i,0]= 255

3 img[j,i,1]= 255

4 img[j,i,2]= 255

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)

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()