@@ -35,109 +35,148 @@ def idx2rc(idx, acc):
3535
3636@jit (nopython = True ) # fill a node (may be two or more points)
3737def fill (img , p , num , nbs , acc , buf ):
38- back = img [p ]
3938 img [p ] = num
4039 buf [0 ] = p
41- cur = 0 ; s = 1 ;
40+ cur = 0 ; s = 1 ; iso = True ;
4241
4342 while True :
4443 p = buf [cur ]
4544 for dp in nbs :
4645 cp = p + dp
47- if img [cp ]== back :
46+ if img [cp ]== 2 :
4847 img [cp ] = num
4948 buf [s ] = cp
5049 s += 1
50+ if img [cp ]== 1 : iso = False
5151 cur += 1
5252 if cur == s :break
53- return idx2rc (buf [:s ], acc )
53+ return iso , idx2rc (buf [:s ], acc )
5454
5555@jit (nopython = True ) # trace the edge and use a buffer, then buf.copy, if use [] numba not works
5656def trace (img , p , nbs , acc , buf ):
5757 c1 = 0 ; c2 = 0 ;
5858 newp = 0
59- cur = 0
60-
59+ cur = 1
6160 while True :
6261 buf [cur ] = p
6362 img [p ] = 0
6463 cur += 1
6564 for dp in nbs :
6665 cp = p + dp
6766 if img [cp ] >= 10 :
68- if c1 == 0 :c1 = img [cp ]
69- else : c2 = img [cp ]
67+ if c1 == 0 :
68+ c1 = img [cp ]
69+ buf [0 ] = cp
70+ else :
71+ c2 = img [cp ]
72+ buf [cur ] = cp
7073 if img [cp ] == 1 :
7174 newp = cp
7275 p = newp
7376 if c2 != 0 :break
74- return (c1 - 10 , c2 - 10 , idx2rc (buf [:cur ], acc ))
77+ return (c1 - 10 , c2 - 10 , idx2rc (buf [:cur + 1 ], acc ))
7578
7679@jit (nopython = True ) # parse the image then get the nodes and edges
77- def parse_struc (img , pts , nbs , acc ):
80+ def parse_struc (img , nbs , acc , iso , ring ):
7881 img = img .ravel ()
7982 buf = np .zeros (131072 , dtype = np .int64 )
8083 num = 10
8184 nodes = []
82- for p in pts :
85+ for p in range ( len ( img )) :
8386 if img [p ] == 2 :
84- nds = fill (img , p , num , nbs , acc , buf )
87+ isiso , nds = fill (img , p , num , nbs , acc , buf )
88+ if isiso and not iso : continue
8589 num += 1
8690 nodes .append (nds )
87-
8891 edges = []
89- for p in pts :
92+ for p in range (len (img )):
93+ if img [p ] < 10 : continue
94+ for dp in nbs :
95+ if img [p + dp ]== 1 :
96+ edge = trace (img , p + dp , nbs , acc , buf )
97+ edges .append (edge )
98+ if not ring : return nodes , edges
99+ for p in range (len (img )):
100+ if img [p ]!= 1 : continue
101+ img [p ] = num ; num += 1
102+ nodes .append (idx2rc ([p ], acc ))
90103 for dp in nbs :
91104 if img [p + dp ]== 1 :
92105 edge = trace (img , p + dp , nbs , acc , buf )
93106 edges .append (edge )
94107 return nodes , edges
95108
96109# use nodes and edges build a networkx graph
97- def build_graph (nodes , edges , multi = False ):
110+ def build_graph (nodes , edges , multi = False , full = True ):
111+ os = np .array ([i .mean (axis = 0 ) for i in nodes ])
112+ if full : os = os .round ().astype (np .uint16 )
98113 graph = nx .MultiGraph () if multi else nx .Graph ()
99114 for i in range (len (nodes )):
100- graph .add_node (i , pts = nodes [i ], o = nodes [i ]. mean ( axis = 0 ) )
115+ graph .add_node (i , pts = nodes [i ], o = os [i ])
101116 for s ,e ,pts in edges :
117+ if full : pts [[0 ,- 1 ]] = os [[s ,e ]]
102118 l = np .linalg .norm (pts [1 :]- pts [:- 1 ], axis = 1 ).sum ()
103119 graph .add_edge (s ,e , pts = pts , weight = l )
104120 return graph
105121
106- def buffer (ske ):
107- buf = np .zeros (tuple (np .array (ske .shape )+ 2 ), dtype = np .uint16 )
108- buf [tuple ([slice (1 ,- 1 )]* buf .ndim )] = ske
122+ def mark_node (ske ):
123+ buf = np .pad (ske , (1 ,1 ), mode = 'constant' )
124+ nbs = neighbors (buf .shape )
125+ acc = np .cumprod ((1 ,)+ buf .shape [::- 1 ][:- 1 ])[::- 1 ]
126+ mark (buf , nbs )
109127 return buf
110-
111- def build_sknw (ske , multi = False ):
112- buf = buffer (ske )
128+
129+ def build_sknw (ske , multi = False , iso = True , ring = True , full = True ):
130+ buf = np . pad (ske , ( 1 , 1 ), mode = 'constant' )
113131 nbs = neighbors (buf .shape )
114132 acc = np .cumprod ((1 ,)+ buf .shape [::- 1 ][:- 1 ])[::- 1 ]
115133 mark (buf , nbs )
116- pts = np .array (np .where (buf .ravel ()== 2 ))[0 ]
117- nodes , edges = parse_struc (buf , pts , nbs , acc )
118- return build_graph (nodes , edges , multi )
134+ nodes , edges = parse_struc (buf , nbs , acc , iso , ring )
135+ return build_graph (nodes , edges , multi , full )
119136
120137# draw the graph
121138def draw_graph (img , graph , cn = 255 , ce = 128 ):
122139 acc = np .cumprod ((1 ,)+ img .shape [::- 1 ][:- 1 ])[::- 1 ]
123140 img = img .ravel ()
141+ for (s , e ) in graph .edges ():
142+ eds = graph [s ][e ]
143+ if isinstance (graph , nx .MultiGraph ):
144+ for i in eds :
145+ pts = eds [i ]['pts' ]
146+ img [np .dot (pts , acc )] = ce
147+ else : img [np .dot (eds ['pts' ], acc )] = ce
124148 for idx in graph .nodes ():
125149 pts = graph .nodes [idx ]['pts' ]
126150 img [np .dot (pts , acc )] = cn
127- for (s , e ) in graph .edges ():
128- eds = graph [s ][e ]
129- for i in eds :
130- pts = eds [i ]['pts' ]
131- img [np .dot (pts , acc )] = ce
132151
133152if __name__ == '__main__' :
134- g = nx .MultiGraph ()
135- g .add_nodes_from ([1 ,2 ,3 ,4 ,5 ])
136- g .add_edges_from ([(1 ,2 ),(1 ,3 ),(2 ,3 ),(4 ,5 ),(5 ,4 )])
137- print (g .nodes ())
138- print (g .edges ())
139- a = g .subgraph (1 )
140- print ('d' )
141- print (a )
142- print ('d' )
143-
153+ import matplotlib .pyplot as plt
154+
155+ img = np .array ([
156+ [0 ,0 ,0 ,1 ,0 ,0 ,0 ,0 ,0 ],
157+ [0 ,0 ,0 ,1 ,0 ,0 ,0 ,1 ,0 ],
158+ [0 ,0 ,0 ,1 ,0 ,0 ,0 ,0 ,0 ],
159+ [1 ,1 ,1 ,1 ,0 ,0 ,0 ,0 ,0 ],
160+ [0 ,0 ,0 ,0 ,1 ,0 ,0 ,0 ,0 ],
161+ [0 ,1 ,0 ,0 ,0 ,1 ,0 ,0 ,0 ],
162+ [1 ,0 ,1 ,0 ,0 ,1 ,1 ,1 ,1 ],
163+ [0 ,1 ,0 ,0 ,1 ,0 ,0 ,0 ,0 ],
164+ [0 ,0 ,0 ,1 ,0 ,0 ,0 ,0 ,0 ]])
165+
166+ node_img = mark_node (img )
167+ graph = build_sknw (img , False , iso = True , ring = True )
168+ plt .imshow (node_img [1 :- 1 ,1 :- 1 ], cmap = 'gray' )
169+
170+ # draw edges by pts
171+ for (s ,e ) in graph .edges ():
172+ ps = graph [s ][e ]['pts' ]
173+ plt .plot (ps [:,1 ], ps [:,0 ], 'green' )
174+
175+ # draw node by o
176+ nodes = graph .nodes ()
177+ ps = np .array ([nodes [i ]['o' ] for i in nodes ])
178+ plt .plot (ps [:,1 ], ps [:,0 ], 'r.' )
179+
180+ # title and show
181+ plt .title ('Build Graph' )
182+ plt .show ()
0 commit comments