bda1bd09
Pavel Govyadinov
Added a network c...
|
1
2
3
4
5
6
7
8
9
|
# -*- coding: utf-8 -*-
"""
Created on Sat Sep 16 16:34:49 2017
@author: pavel
"""
import struct
import numpy as np
|
3df117ca
David Mayerich
edited spherical ...
|
10
|
import scipy as sp
|
bda1bd09
Pavel Govyadinov
Added a network c...
|
11
12
13
|
import networkx as nx
import matplotlib.pyplot as plt
import math
|
db331d8a
David Mayerich
added adjacency w...
|
14
15
|
import time
import spharmonics
|
bda1bd09
Pavel Govyadinov
Added a network c...
|
16
17
18
19
20
21
22
23
24
25
26
27
|
'''
Definition of the Node class
Duplicate of the node class in network
Stores the physical position, outgoing edges list and incoming edges list.
'''
class Node:
def __init__(self, point, outgoing, incoming):
self.p = point
self.o = outgoing
self.i = incoming
|
1ae6fbe0
David Mayerich
error fixes in ne...
|
28
29
30
|
# def p():
# return self.p
|
bda1bd09
Pavel Govyadinov
Added a network c...
|
31
32
33
34
35
36
37
|
'''
Definition of the Fiber class.
Duplicate of the Node class in network
Stores the starting vertex, the ending vertex, the points array and the radius array
'''
class Fiber:
|
db331d8a
David Mayerich
added adjacency w...
|
38
39
40
41
42
43
|
def __init__ (self, p1, p2, pois, rads):
self.v0 = p1
self.v1 = p2
self.points = pois
self.radii = rads
|
bda1bd09
Pavel Govyadinov
Added a network c...
|
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
|
'''
return the length of the fiber.
'''
def length(self):
length = 0
for i in range(len(self.points)-1):
length = length + math.sqrt(pow(self.points[i][0]- self.points[i+1][0],2) + pow(self.points[i][1]- self.points[i+1][1],2) + pow(self.points[i][2]- self.points[i+1][2],2))
return length
'''
returns the turtuosity of the fiber.
'''
def turtuosity(self):
turtuosity = 0
distance = math.sqrt(math.pow(self.points[0][0]- self.points[len(self.points)-1][0],2) + math.pow(self.points[0][1]- self.points[len(self.points)-1][1],2) + math.pow(self.points[0][2]- self.points[len(self.points)-1][2],2))
turtuosity = self.length()/distance
#print(turtuosity)
return turtuosity
'''
returns the volume of the fiber.
'''
def volume(self):
volume = 0
for i in range(len(self.points)-1):
volume = volume + 1.0/3.0 * math.pi * (math.pow(self.radii[i],2) + math.pow(self.radii[i+1],2) + self.radii[i]*self.radii[i+1]) * math.sqrt(math.pow(self.points[i][0]- self.points[i+1][0],2) + math.pow(self.points[i][1]- self.points[i+1][1],2) + math.pow(self.points[i][2]- self.points[i+1][2],2))
#print(volume)
return volume
|
bda1bd09
Pavel Govyadinov
Added a network c...
|
76
|
|
8d96a467
David Mayerich
started convertin...
|
77
78
79
80
81
|
class NWT:
'''
Writes the header given and open file descripion, number of verticies and number of edges.
'''
|
db331d8a
David Mayerich
added adjacency w...
|
82
|
def writeHeader(open_file, numVerts, numEdges):
|
8d96a467
David Mayerich
started convertin...
|
83
84
85
86
87
88
|
txt = "nwtFileFormat fileid(14B), desc(58B), #vertices(4B), #edges(4B): bindata"
b = bytearray()
b.extend(txt.encode())
open_file.write(b)
open_file.write(struct.pack('i', numVerts))
open_file.write(struct.pack('i', numEdges))
|
bda1bd09
Pavel Govyadinov
Added a network c...
|
89
|
|
8d96a467
David Mayerich
started convertin...
|
90
91
92
93
|
'''
Writes a single vertex to a file.
'''
|
db331d8a
David Mayerich
added adjacency w...
|
94
|
def writeVertex(open_file, vertex):
|
8d96a467
David Mayerich
started convertin...
|
95
96
97
98
99
100
101
|
open_file.write(struct.pack('<f',vertex.p[0]))
open_file.write(struct.pack('<f',vertex.p[1]))
open_file.write(struct.pack('<f',vertex.p[2]))
open_file.write(struct.pack('i', len(vertex.o)))
open_file.write(struct.pack('i', len(vertex.i)))
for j in range(len(vertex.o)):
open_file.write(struct.pack('i',vertex.o[j]))
|
bda1bd09
Pavel Govyadinov
Added a network c...
|
102
|
|
8d96a467
David Mayerich
started convertin...
|
103
104
|
for j in range(len(vertex.i)):
open_file.write(struct.pack('i', vertex.i[j]))
|
bda1bd09
Pavel Govyadinov
Added a network c...
|
105
|
|
8d96a467
David Mayerich
started convertin...
|
106
|
return
|
bda1bd09
Pavel Govyadinov
Added a network c...
|
107
|
|
8d96a467
David Mayerich
started convertin...
|
108
109
110
|
'''
Writes a single fiber to a file.
'''
|
db331d8a
David Mayerich
added adjacency w...
|
111
|
def writeFiber(open_file, edge):
|
8d96a467
David Mayerich
started convertin...
|
112
113
114
115
116
117
118
119
120
121
|
open_file.write(struct.pack('i',edge.v0))
open_file.write(struct.pack('i',edge.v1))
open_file.write(struct.pack('i',len(edge.points)))
for j in range(len(edge.points)):
open_file.write(struct.pack('<f', edge.points[j][0]))
open_file.write(struct.pack('<f', edge.points[j][1]))
open_file.write(struct.pack('<f', edge.points[j][2]))
open_file.write(struct.pack('<f', edge.radii[j]))
return
|
bda1bd09
Pavel Govyadinov
Added a network c...
|
122
|
|
8d96a467
David Mayerich
started convertin...
|
123
124
125
|
'''
Writes the entire network to a file in str given the vertices array and the edges array.
'''
|
db331d8a
David Mayerich
added adjacency w...
|
126
|
def exportNWT(str, vertices, edges):
|
8d96a467
David Mayerich
started convertin...
|
127
|
with open(str, "wb") as file:
|
db331d8a
David Mayerich
added adjacency w...
|
128
|
NWT.writeHeader(file, len(vertices), len(edges))
|
8d96a467
David Mayerich
started convertin...
|
129
|
for i in range(len(vertices)):
|
db331d8a
David Mayerich
added adjacency w...
|
130
|
NWT.writeVertex(file, vertices[i])
|
8d96a467
David Mayerich
started convertin...
|
131
132
|
for i in range(len(edges)):
|
db331d8a
David Mayerich
added adjacency w...
|
133
|
NWT.writeFiber(file, edges[i])
|
8d96a467
David Mayerich
started convertin...
|
134
135
|
return
|
bda1bd09
Pavel Govyadinov
Added a network c...
|
136
|
|
bda1bd09
Pavel Govyadinov
Added a network c...
|
137
|
|
8d96a467
David Mayerich
started convertin...
|
138
139
140
|
'''
Reads a single vertex from an open file and returns a node Object.
'''
|
db331d8a
David Mayerich
added adjacency w...
|
141
|
def readVertex(open_file):
|
8d96a467
David Mayerich
started convertin...
|
142
|
points = np.tile(0., 3)
|
bda1bd09
Pavel Govyadinov
Added a network c...
|
143
|
bytes = open_file.read(4)
|
8d96a467
David Mayerich
started convertin...
|
144
|
points[0] = struct.unpack('f', bytes)[0]
|
bda1bd09
Pavel Govyadinov
Added a network c...
|
145
|
bytes = open_file.read(4)
|
8d96a467
David Mayerich
started convertin...
|
146
|
points[1] = struct.unpack('f', bytes)[0]
|
bda1bd09
Pavel Govyadinov
Added a network c...
|
147
|
bytes = open_file.read(4)
|
8d96a467
David Mayerich
started convertin...
|
148
|
points[2] = struct.unpack('f', bytes)[0]
|
bda1bd09
Pavel Govyadinov
Added a network c...
|
149
|
bytes = open_file.read(4)
|
bda1bd09
Pavel Govyadinov
Added a network c...
|
150
|
|
8d96a467
David Mayerich
started convertin...
|
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
|
numO = int.from_bytes(bytes, byteorder='little')
outgoing = np.tile(0, numO)
bts = open_file.read(4)
numI = int.from_bytes(bts, byteorder='little')
incoming = np.tile(0, numI)
for j in range(numO):
bytes = open_file.read(4)
outgoing[j] = int.from_bytes(bytes, byteorder='little')
for j in range(numI):
bytes = open_file.read(4)
incoming[j] = int.from_bytes(bytes, byteorder='little')
node = Node(points, outgoing, incoming)
return node
|
bda1bd09
Pavel Govyadinov
Added a network c...
|
166
|
|
bda1bd09
Pavel Govyadinov
Added a network c...
|
167
|
|
8d96a467
David Mayerich
started convertin...
|
168
169
170
|
'''
Reads a single fiber from an open file and returns a Fiber object .
'''
|
db331d8a
David Mayerich
added adjacency w...
|
171
|
def readFiber(open_file):
|
8d96a467
David Mayerich
started convertin...
|
172
173
174
175
176
177
178
179
|
bytes = open_file.read(4)
vtx0 = int.from_bytes(bytes, byteorder = 'little')
bytes = open_file.read(4)
vtx1 = int.from_bytes(bytes, byteorder = 'little')
bytes = open_file.read(4)
numVerts = int.from_bytes(bytes, byteorder = 'little')
pts = []
rads = []
|
bda1bd09
Pavel Govyadinov
Added a network c...
|
180
|
|
8d96a467
David Mayerich
started convertin...
|
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
|
for j in range(numVerts):
point = np.tile(0., 3)
bytes = open_file.read(4)
point[0] = struct.unpack('f', bytes)[0]
bytes = open_file.read(4)
point[1] = struct.unpack('f', bytes)[0]
bytes = open_file.read(4)
point[2] = struct.unpack('f', bytes)[0]
bytes = open_file.read(4)
radius = struct.unpack('f', bytes)[0]
pts.append(point)
rads.append(radius)
F = Fiber(vtx0, vtx1, pts, rads)
return F
|
bda1bd09
Pavel Govyadinov
Added a network c...
|
197
|
|
8d96a467
David Mayerich
started convertin...
|
198
199
200
201
|
'''
Imports a NWT file at location str.
Returns a list of Nodes objects and a list of Fiber objects.
'''
|
3df117ca
David Mayerich
edited spherical ...
|
202
|
|
8d96a467
David Mayerich
started convertin...
|
203
|
class Network:
|
3df117ca
David Mayerich
edited spherical ...
|
204
|
|
db331d8a
David Mayerich
added adjacency w...
|
205
206
207
208
209
|
def __init__(self, filename, clock=False):
if clock:
start_time = time.time()
with open(filename, "rb") as file:
|
8d96a467
David Mayerich
started convertin...
|
210
211
212
213
214
|
header = file.read(72)
bytes = file.read(4)
numVertex = int.from_bytes(bytes, byteorder='little')
bytes = file.read(4)
numEdges = int.from_bytes(bytes, byteorder='little')
|
3df117ca
David Mayerich
edited spherical ...
|
215
|
|
8d96a467
David Mayerich
started convertin...
|
216
217
218
219
220
|
self.N = []
self.F = []
for i in range(numVertex):
node = NWT.readVertex(file)
self.N.append(node)
|
3df117ca
David Mayerich
edited spherical ...
|
221
|
|
8d96a467
David Mayerich
started convertin...
|
222
223
224
|
for i in range(numEdges):
edge = NWT.readFiber(file)
self.F.append(edge)
|
db331d8a
David Mayerich
added adjacency w...
|
225
226
|
if clock:
print("Network initialization: " + str(time.time() - start_time) + "s")
|
3df117ca
David Mayerich
edited spherical ...
|
227
|
|
8d96a467
David Mayerich
started convertin...
|
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
|
'''
Creates a graph from a list of nodes and a list of edges.
Uses edge length as weight.
Returns a NetworkX Object.
'''
# def createLengthGraph(self):
# G = nx.Graph()
# for i in range(len(self.nodeList)):
# G.add_node(i, p=V[i].p)
# for i in range(len(self.edgeList)):
# G.add_edge(self.edgeList[i].v0, self.edgeList[i].v1, weight = E[i].length())
#
# return G
# '''
# Creates a graph from a list of nodes and a list of edges.
# Uses edge turtuosity as weight.
# Returns a NetworkX Object.
# '''
# def createTortuosityGraph(nodeList, edgeList):
# G = nx.Graph()
# for i in range(len(nodeList)):
# G.add_node(i, p=V[i].p)
# for i in range(len(edgeList)):
# G.add_edge(edgeList[i].v0, edgeList[i].v1, weight = E[i].turtuosity())
#
# return G
|
3df117ca
David Mayerich
edited spherical ...
|
254
|
|
8d96a467
David Mayerich
started convertin...
|
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
|
# '''
# Creates a graph from a list of nodes and a list of edges.
# Uses edge volume as weight.
# Returns a NetworkX Object.
# '''
# def createVolumeGraph(nodeList, edgeList):
# G = nx.Graph()
# for i in range(len(nodeList)):
# G.add_node(i, p=V[i].p)
# for i in range(len(edgeList)):
# G.add_edge(edgeList[i].v0, edgeList[i].v1, weight = E[i].volume())
#
# return G
#'''
#Returns the positions dictionary for the Circular layout.
#'''
#def getCircularLayout(graph, dim, radius):
# return nx.circular_layout(graph, dim, radius)
#
#'''
#Return the positions dictionary for the Spring layout.
#'''
#def getSpringLayout(graph, pos, iterations, scale):
# return nx.spring_layout(graph, 2, None, pos, iterations, 'weight', scale, None)
#
#'''
#Draws the graph.
#'''
#def drawGraph(graph, pos):
# nx.draw(graph, pos)
# return
def aabb(self):
|
3df117ca
David Mayerich
edited spherical ...
|
288
|
|
8d96a467
David Mayerich
started convertin...
|
289
290
|
lower = self.N[0].p.copy()
upper = lower.copy()
|
db331d8a
David Mayerich
added adjacency w...
|
291
|
for i in self.N:
|
8d96a467
David Mayerich
started convertin...
|
292
293
294
295
296
297
|
for c in range(len(lower)):
if lower[c] > i.p[c]:
lower[c] = i.p[c]
if upper[c] < i.p[c]:
upper[c] = i.p[c]
return lower, upper
|
3df117ca
David Mayerich
edited spherical ...
|
298
|
|
8d96a467
David Mayerich
started convertin...
|
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
|
#calculate the distance field at a given resolution
# R (triple) resolution along each dimension
def distancefield(self, R=(100, 100, 100)):
#get a list of all node positions in the network
P = []
for e in self.F:
for p in e.points:
P.append(p)
#turn that list into a Numpy array so that we can create a KD tree
P = np.array(P)
#generate a KD-Tree out of the network point array
tree = sp.spatial.cKDTree(P)
plt.scatter(P[:, 0], P[:, 1])
#specify the resolution of the ouput grid
R = (200, 200, 200)
#generate a meshgrid of the appropriate size and resolution to surround the network
|
db331d8a
David Mayerich
added adjacency w...
|
321
|
lower, upper = self.aabb(self.N, self.F) #get the space occupied by the network
|
8d96a467
David Mayerich
started convertin...
|
322
323
324
325
326
327
328
329
330
331
332
333
334
|
x = np.linspace(lower[0], upper[0], R[0]) #get the grid points for uniform sampling of this space
y = np.linspace(lower[1], upper[1], R[1])
z = np.linspace(lower[2], upper[2], R[2])
X, Y, Z = np.meshgrid(x, y, z)
#Z = 150 * numpy.ones(X.shape)
Q = np.stack((X, Y, Z), 3)
D, I = tree.query(Q)
return D
|
3df117ca
David Mayerich
edited spherical ...
|
335
|
|
8d96a467
David Mayerich
started convertin...
|
336
|
#returns the number of points in the network
|
db331d8a
David Mayerich
added adjacency w...
|
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
|
def npoints(self):
n = 0 #initialize the counter to zero
for f in self.F: #for each fiber
n = n + len(f.points) - 2 #count the number of points in the fiber - ignoring the end points
n = n + len(self.N) #add the number of nodes (shared points) to the node count
return n #return the number of nodes
#returns all of the points in the network
def points(self):
k = self.npoints()
P = np.zeros((3, k)) #allocate space for the point list
idx = 0
for f in self.F: #for each fiber in the network
for ip in range(1, len(f.points)-1): #for each point in the network
P[:, idx] = f.points[ip] #store the point in the raw point list
idx = idx + 1
return P #return the point array
|
3df117ca
David Mayerich
edited spherical ...
|
355
|
|
8d96a467
David Mayerich
started convertin...
|
356
357
|
#returns the number of linear segments in the network
def nsegments(self):
|
db331d8a
David Mayerich
added adjacency w...
|
358
359
360
361
|
n = 0 #initialize the segment counter to 0
for f in self.F: #for each fiber
n = n + len(f.points) - 1 #calculate the number of line segments in the fiber (points - 1)
return n #return the number of line segments
|
3df117ca
David Mayerich
edited spherical ...
|
362
|
|
db331d8a
David Mayerich
added adjacency w...
|
363
364
365
366
367
368
369
370
371
372
373
374
|
#return a list of line segments representing the network
def segments(self, dtype=np.float32):
k = self.nsegments() #get the number of line segments
start = np.zeros((k, 3),dtype=dtype) #start points for the line segments
end = np.zeros((k, 3), dtype=dtype) #end points for the line segments
idx = 0 #initialize the index counter to zero
for f in self.F: #for each fiber in the network
for ip in range(0, len(f.points)-1): #for each point in the network
start[idx, :] = f.points[ip] #store the point in the raw point list
idx = idx + 1
|
8d96a467
David Mayerich
started convertin...
|
375
|
idx = 0
|
db331d8a
David Mayerich
added adjacency w...
|
376
377
378
|
for f in self.F: #for each fiber in the network
for ip in range(1, len(f.points)): #for each point in the network
end[idx, :] = f.points[ip] #store the point in the raw point list
|
8d96a467
David Mayerich
started convertin...
|
379
|
idx = idx + 1
|
db331d8a
David Mayerich
added adjacency w...
|
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
|
return start, end
#function returns the fiber associated with a given 1D line segment index
def segment2fiber(self, idx):
i = 0
for f in range(len(self.F)): #for each fiber in the network
i = i + len(self.F[f].points)-1 #add the number of points in the fiber to i
if i > idx: #if we encounter idx in this fiber
return self.F[f].points, f #return the fiber associated with idx and the index into the fiber array
def vectors(self, clock=False, dtype=np.float32):
if clock:
start_time = time.time()
start, end = self.segments(dtype) #retrieve all of the line segments
v = end - start #calculate the resulting vectors
l = np.sqrt(v[:, 0]**2 + v[:,1]**2 + v[:,2]**2) #calculate the fiber lengths
z = l==0 #look for any zero values
nz = z.sum()
if nz > 0:
print("WARNING: " + str(nz) + " line segment(s) of length zero were found in the network and will be removed" )
if clock:
print("Network::vectors: " + str(time.time() - start_time) + "s")
return np.delete(v, np.where(z), 0)
#scale all values in the network by tuple S = (sx, sy, sz)
def scale(self, S):
for f in self.F:
for p in f.points:
p[0] = p[0] * S[0]
p[1] = p[1] * S[1]
p[2] = p[2] * S[2]
for n in self.N:
n.p[0] = n.p[0] * S[0]
n.p[1] = n.p[1] * S[1]
n.p[2] = n.p[2] * S[2]
#calculate the adjacency weighting function for the network given a set of vectors X = (x, y, z) and weight exponent k
def adjacencyweight(self, P, k=200, length_threshold = 25, dtype=np.float32):
V = self.vectors(dtype) #get the vectors representing each segment
#V = V[0:n_vectors, :]
L = np.expand_dims(np.sqrt((V**2).sum(1)), 1) #calculate the length of each vector
outliers = L > length_threshold #remove outliers based on the length_threshold
V = np.delete(V, np.where(outliers), 0)
L = np.delete(L, np.where(outliers))
V = V/L[:,None] #normalize the vectors
P = np.stack(spharmonics.sph2cart(1, P[0], P[1]), P[0].ndim)
PV = P[...,None,:] * V
cos_alpha = PV.sum(PV.ndim-1)
W = np.abs(cos_alpha) ** k
return W, L
|