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

拓扑重建

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

不是太完美的图，应该是细化出问题了

推广一下大佬的开源项目：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 （英文文档）

相关阅读:
解决PKIX：unable to find valid certification path to requested target 的问题
 Linux 上的常用文件传输方式介绍与比较
 用VNC远程图形化连接Linux桌面的配置方法
 红帽中出现”This system is not registered with RHN”的解决方案
 linux安装时出现your cpu does not support long mode的解决方法
 CentOS SSH配置
 es6扩展运算符及rest运算符总结
 es6解构赋值总结
 tortoisegit安装、clon、推送
 es6环境搭建
原文地址：https://www.cnblogs.com/wjy-lulu/p/7242788.html

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

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