Commit 68f91f03488fe3943dc0c1941a570f9f494136ca
1 parent
3cc9b7dd
added fibernet for sanity
Showing
1 changed file
with
222 additions
and
0 deletions
Show diff stats
1 | +# -*- coding: utf-8 -*- | ||
2 | +""" | ||
3 | +Created on Sat Jan 19 2018 | ||
4 | + | ||
5 | +@author: Jiabing | ||
6 | +""" | ||
7 | + | ||
8 | +import struct | ||
9 | +import numpy as np | ||
10 | +import scipy as sp | ||
11 | +import networkx as nx | ||
12 | +import matplotlib.pyplot as plt | ||
13 | +import math | ||
14 | +from mpl_toolkits.mplot3d import Axes3D | ||
15 | +from itertools import chain | ||
16 | + | ||
17 | +class Node: | ||
18 | + def __init__(self, point, outgoing, incoming): | ||
19 | + self.p = point | ||
20 | + self.o = outgoing | ||
21 | + self.i = incoming | ||
22 | + | ||
23 | + | ||
24 | +class Fiber: | ||
25 | + def __init__ (self, indices): | ||
26 | + self.indices = indices | ||
27 | + | ||
28 | +class Point: | ||
29 | + def __init__(self, x, y, z, r): | ||
30 | + self.x = x | ||
31 | + self.y = y | ||
32 | + self.z = z | ||
33 | + self.r = r | ||
34 | + | ||
35 | +class Edge: | ||
36 | + def __init__(self, indices_num, pois, radius): | ||
37 | + self.indices_num = indices_num | ||
38 | + self.points = pois | ||
39 | + self.radius = radius | ||
40 | + | ||
41 | + | ||
42 | +class NWT: | ||
43 | + | ||
44 | + def readVertex(open_file): | ||
45 | + points = np.tile(0., 3) | ||
46 | + bytes = open_file.read(4) | ||
47 | + points[0] = struct.unpack('f', bytes)[0] | ||
48 | + bytes = open_file.read(4) | ||
49 | + points[1] = struct.unpack('f', bytes)[0] | ||
50 | + bytes = open_file.read(4) | ||
51 | + points[2] = struct.unpack('f', bytes)[0] | ||
52 | + bytes = open_file.read(4) | ||
53 | + | ||
54 | + numO = int.from_bytes(bytes, byteorder='little') | ||
55 | + outgoing = np.tile(0, numO) | ||
56 | + bts = open_file.read(4) | ||
57 | + numI = int.from_bytes(bts, byteorder='little') | ||
58 | + incoming = np.tile(0, numI) | ||
59 | + for j in range(numO): | ||
60 | + bytes = open_file.read(4) | ||
61 | + outgoing[j] = int.from_bytes(bytes, byteorder='little') | ||
62 | + | ||
63 | + for j in range(numI): | ||
64 | + bytes = open_file.read(4) | ||
65 | + incoming[j] = int.from_bytes(bytes, byteorder='little') | ||
66 | + | ||
67 | + node = Node(points, outgoing, incoming) | ||
68 | + return node | ||
69 | + | ||
70 | + | ||
71 | + ''' | ||
72 | + Reads a single fiber from an open file and returns a Fiber object . | ||
73 | + ''' | ||
74 | + def readFiber(open_file): | ||
75 | + bytes = open_file.read(4) | ||
76 | + vtx0 = int.from_bytes(bytes, byteorder = 'little') | ||
77 | + bytes = open_file.read(4) | ||
78 | + vtx1 = int.from_bytes(bytes, byteorder = 'little') | ||
79 | + bytes = open_file.read(4) | ||
80 | + numVerts = int.from_bytes(bytes, byteorder = 'little') | ||
81 | + pts = [] | ||
82 | + rads = [] | ||
83 | + | ||
84 | + for j in range(numVerts): | ||
85 | + point = np.tile(0., 3) | ||
86 | + bytes = open_file.read(4) | ||
87 | + point[0] = struct.unpack('f', bytes)[0] | ||
88 | + bytes = open_file.read(4) | ||
89 | + point[1] = struct.unpack('f', bytes)[0] | ||
90 | + bytes = open_file.read(4) | ||
91 | + point[2] = struct.unpack('f', bytes)[0] | ||
92 | + bytes = open_file.read(4) | ||
93 | + radius = struct.unpack('f', bytes)[0] | ||
94 | + pts.append(point) | ||
95 | + rads.append(radius) | ||
96 | + | ||
97 | + # F = Fiber(pts) | ||
98 | + E = Edge(numVerts, pts, radius) | ||
99 | + #E = Edge(pts) | ||
100 | + return E | ||
101 | + | ||
102 | + | ||
103 | +class fibernet: | ||
104 | + def __init__(self, filename): | ||
105 | + | ||
106 | + with open(filename, "rb") as file: | ||
107 | + header = file.read(72) | ||
108 | + bytes = file.read(4) | ||
109 | + numVertex = int.from_bytes(bytes, byteorder='little') | ||
110 | + bytes = file.read(4) | ||
111 | + numEdges = int.from_bytes(bytes, byteorder='little') | ||
112 | + | ||
113 | + #self.P = [] | ||
114 | + self.F = [] | ||
115 | + self.N = [] | ||
116 | + | ||
117 | + for i in range(numVertex): | ||
118 | + node = NWT.readVertex(file) | ||
119 | + self.N.append(node) | ||
120 | + | ||
121 | + for i in range( numEdges): | ||
122 | + edge = NWT.readFiber(file) | ||
123 | + #self.F.append(np.arange(num,num+edge.indices_num,1)) | ||
124 | + #self.P= chain(self.P, edge.points) | ||
125 | + self.F.append(edge.points) | ||
126 | + #num += edge.indices_num | ||
127 | + | ||
128 | + def aabb(self): | ||
129 | + | ||
130 | + lower = self.N[0].p.copy() | ||
131 | + upper = lower.copy() | ||
132 | + for i in self.N: | ||
133 | + for c in range(len(lower)): | ||
134 | + if lower[c] > i.p[c]: | ||
135 | + lower[c] = i.p[c] | ||
136 | + if upper[c] < i.p[c]: | ||
137 | + upper[c] = i.p[c] | ||
138 | + return lower, upper | ||
139 | + | ||
140 | + | ||
141 | + | ||
142 | + def distancefield(self, R=(100, 100, 100)): | ||
143 | + | ||
144 | + #generate a meshgrid of the appropriate size and resolution to surround the network | ||
145 | + lower, upper = self.aabb() #get the space occupied by the network | ||
146 | + x = np.linspace(lower[0], upper[0], R[0]) #get the grid points for uniform sampling of this space | ||
147 | + y = np.linspace(lower[1], upper[1], R[1]) | ||
148 | + z = np.linspace(lower[2], upper[2], R[2]) | ||
149 | + X, Y, Z = np.meshgrid(x, y, z) | ||
150 | + #Z = 150 * numpy.ones(X.shape) | ||
151 | + | ||
152 | + Q = np.stack((X, Y, Z), 3) | ||
153 | + d_x = abs(x[1]-x[0]); | ||
154 | + d_y = abs(y[1]-y[0]); | ||
155 | + d_z = abs(z[1]-z[0]); | ||
156 | + dis1 = math.sqrt(pow(d_x,2)+pow(d_y,2)+pow(d_z ,2)) | ||
157 | + #dx = abs(x[1]-x[0]) | ||
158 | + | ||
159 | + #dy = abs(y[1]-y[0]) | ||
160 | + #dz = abs(z[1]-z[0]) | ||
161 | + #get a list of all node positions in the network | ||
162 | + P = [] | ||
163 | + | ||
164 | + for e in self.F[12:13]: #12-17 | ||
165 | + for p in e: | ||
166 | + P.append(p) | ||
167 | + | ||
168 | + for j in range(len(e)-1): | ||
169 | + d_t = e[j+1]-e[j] | ||
170 | + dis2 = math.sqrt(pow(d_t[0],2)+pow(d_t[1],2)+pow(d_t[2],2)) | ||
171 | + ins = max(int(d_t[0]/d_x), int(d_t[1]/d_y), int(d_t[2]/d_z)) | ||
172 | + if( ins>0 ): | ||
173 | + ins = ins+1; | ||
174 | + for k in range(ins): | ||
175 | + p_ins =e[j]+(k+1)*(e[j+1]-e[j])/ins; | ||
176 | + P.append(p_ins); | ||
177 | + #turn that list into a Numpy array so that we can create a KD tree | ||
178 | + P = np.array(P) | ||
179 | + | ||
180 | + #generate a KD-Tree out of the network point array | ||
181 | + tree = sp.spatial.cKDTree(P) | ||
182 | + | ||
183 | + #specify the resolution of the ouput grid | ||
184 | + # R = (200, 200, 200) | ||
185 | + | ||
186 | + D, I = tree.query(Q) | ||
187 | + | ||
188 | + | ||
189 | + | ||
190 | + return D, Q, dis1 | ||
191 | + | ||
192 | + | ||
193 | +''' | ||
194 | +##read NWT file | ||
195 | +f= fibernet("full_seg.nwt") | ||
196 | +#P = tuple(f.P) | ||
197 | +plist = f.F | ||
198 | + | ||
199 | +fig = plt.figure() | ||
200 | +ax = fig.add_subplot(111, projection='3d') | ||
201 | + | ||
202 | +xs = [] | ||
203 | +ys = [] | ||
204 | +zs = [] | ||
205 | +for i in range(10): | ||
206 | + for j in range(len(F[i])): | ||
207 | + xs.append(P[F[i][j]][0]) | ||
208 | + ys.append(P[F[i][j]][1]) | ||
209 | + zs.append(P[F[i][j]][2]) | ||
210 | + | ||
211 | + | ||
212 | +#ax.scatter(xs, ys, zs) | ||
213 | +#ax = fig.gca(projection='3d') | ||
214 | +ax.set_xlabel('X Label') | ||
215 | +ax.set_ylabel('Y Label') | ||
216 | +ax.set_zlabel('Z Label') | ||
217 | + | ||
218 | +ax.plot(xs, ys, zs, label='center lines') | ||
219 | +ax.legend() | ||
220 | +plt.savefig('p.png', dpi=100) | ||
221 | +plt.show() | ||
222 | +''' | ||
0 | \ No newline at end of file | 223 | \ No newline at end of file |