Commit e9d3fb64f1ab7e9d4ca122ef2280530509d65ce1

Authored by David Mayerich
1 parent 0eed0e54

added new fiber network class: fibernet

Showing 2 changed files with 202 additions and 0 deletions   Show diff stats
python/fibernet.py 0 → 100644
  1 +# -*- coding: utf-8 -*-
  2 +"""
  3 +Created on Mon Mar 19 14:45:39 2018
  4 +
  5 +@author: david
  6 +"""
  7 +
  8 +import nwt
  9 +import matplotlib.pyplot as plt
  10 +import numpy
  11 +
  12 +
  13 +class fibernet:
  14 +
  15 + def __init__(self, filename):
  16 + self.nwt(filename)
  17 + self.buildgraph()
  18 +
  19 + #load a network from a NWT file
  20 + def nwt(self, filename):
  21 + net = nwt.NWT(filename)
  22 +
  23 + #first insert the points where fibers are connected
  24 + self.P = []
  25 + for vi in range(0, len(net.v)): #store each vertex in the point list
  26 + self.P.append( (net.v[vi].p[0], net.v[vi].p[1], net.v[vi].p[2]) )
  27 +
  28 + #now insert each fiber, connecting them appropriately using previously inserted points (vertices)
  29 + self.F = []
  30 + for ei in range(0, len(net.e)):
  31 + f = [numpy.int(net.e[ei].v[0])] #create a new fiber to store point indices, initialize with the first vertex ID
  32 + for pi in range(1, net.e[ei].p.shape[0]):
  33 + self.P.append( (net.e[ei].p[pi,0], net.e[ei].p[pi,1], net.e[ei].p[pi,2])) #append the current point to the global point list
  34 + f.append(len(self.P) - 1) #append the current point ID
  35 + f.append(numpy.int(net.e[ei].v[1])) #append the last vertex ID
  36 + self.F.append(f)
  37 +
  38 + #build a graph from the available geometry
  39 + def buildgraph(self):
  40 + p_to_v = numpy.ones((self.npoints()), numpy.int) * (-1) #create an array that maps a point ID to a vertex ID
  41 +
  42 + self.V = [] #create an empty vertex list
  43 + self.E = [] #create an empty edge list
  44 + for fi in range(0, self.nfibers()): #for each fiber in the network
  45 + if p_to_v[self.F[fi][0]] == -1: #if the first point hasn't been encountered
  46 + p_to_v[self.F[fi][0]] = len(self.V) #add it to the mapping
  47 + self.V.append( (self.F[fi][0], []) ) #add a new vertex representing this point
  48 + if p_to_v[self.F[fi][-1]] == -1:
  49 + p_to_v[self.F[fi][-1]] = len(self.V)
  50 + self.V.append( (self.F[fi][-1], []) )
  51 +
  52 + v0 = p_to_v[self.F[fi][0]] #get the two vertex indices corresponding to this edge
  53 + v1 = p_to_v[self.F[fi][-1]]
  54 +
  55 + self.E.append( (v0, v1) ) #create an edge representing this fiber and add it to the edge list
  56 + self.V[v0][1].append(fi) #add this edge ID to the corresponding vertices to update the adjacency list
  57 + self.V[v1][1].append(fi)
  58 +
  59 + #return the length of a fiber specified by fiber ID fid
  60 + def length(self, fid):
  61 + p = self.points(fid) #get the points corresponding to the fiber
  62 + l = 0 #initialize the length to zero
  63 + for pi in range(1, len(p)): #for each line segment
  64 + v = p[pi] - p[pi-1] #calculate the vector representing the line segment
  65 + l = l + numpy.linalg.norm(numpy.array(v)) #add the length of the line segment
  66 + return l #return the total length
  67 +
  68 + #return a list of vertices adjacent to the vertex specified by vid
  69 + def adjacent(self, vid):
  70 + adj = [] #initialize an empty adjacency list
  71 + for e in self.V[vid][1]: #for each edge connected to vid
  72 + if self.E[e][0] != vid:
  73 + adj.append(self.E[e][0])
  74 + else:
  75 + adj.append(self.E[e][1])
  76 + return adj
  77 +
  78 + #return the fiber network graph as an adjacency list
  79 + def adjacencylist(self):
  80 + A = []
  81 + for vi in range(0, len(self.V)): #for each vertex in the graph
  82 + A.append( self.adjacent(vi) ) #add the adjacency list for the vertex into the main adjacency list
  83 + return A
  84 +
  85 + #return the fiber network graph as an adjacency matrix
  86 + # possible values for the matrix include:
  87 + # binary - 1 when there is a connection, 0 otherwise
  88 + # edgeID - stores the edge ID corresponding to a connection, and -1 otherwise
  89 + # length - stores the length of the connecting fiber, and -1 if there is no connection
  90 + def adjacencymatrix(self, value="binary"):
  91 + if value == "binary":
  92 + M = numpy.ones( (len(self.V), len(self.V)), numpy.bool )
  93 + elif value == "edgeID":
  94 + M = numpy.ones( (len(self.V), len(self.V)), numpy.int ) * -1
  95 + elif value == "length":
  96 + M = numpy.ones( (len(self.V), len(self.V)) ) * -1
  97 + else:
  98 + return []
  99 +
  100 + for vi in range(0, len(self.V)): #for each vertex in the graph
  101 + A = self.adjacent(vi) #get the adjacency matrix
  102 + for a in A:
  103 + if value == "binary":
  104 + M[vi,a] = M[a,vi] = 1
  105 + if value == "edgeID":
  106 + M[vi,a] = M[a,vi] = self.findedge(a, vi)
  107 + if value == "length":
  108 + M[vi,a] = M[a,vi] = self.length(self.findedge(a,vi))
  109 + return M
  110 +
  111 + #find the edge ID corresponding to two vertices
  112 + def findedge(self, v0, v1):
  113 +
  114 + edges = self.V[v0][1] #get the list of candidate edges connected to v0
  115 + for e in edges: #for each edge
  116 + if self.E[e][0] == v1 or self.E[e][1] == v1:
  117 + return e
  118 +
  119 + return -1 #return -1 if the edge doesn't exist
  120 +
  121 + #get the points associated with a specified fiber ID fid
  122 + def points(self, fid): #get the points corresponding to a specified fiber index
  123 + f = self.F[fid] #get the list of point indices
  124 + return numpy.array(self.P)[f]
  125 +
  126 + #get the total number of points in the network
  127 + def npoints(self):
  128 + return len(self.P)
  129 +
  130 + #get the total number of fibers in the network
  131 + def nfibers(self):
  132 + return len(self.F)
  133 +
  134 + #plot the network as a 3D line plot
  135 + def plot(self):
  136 + fig = plt.figure()
  137 + ax = fig.add_subplot(111, projection='3d')
  138 +
  139 + for fi in range(0, self.nfibers()):
  140 + test = fibers.points(fi)
  141 + ax.plot(test[:, 0], test[:, 1], test[:, 2])
  142 + return fig
0 143 \ No newline at end of file
... ...
python/nwt.py 0 → 100644
  1 +# -*- coding: utf-8 -*-
  2 +"""
  3 +Created on Mon Mar 19 12:38:30 2018
  4 +
  5 +@author: david
  6 +"""
  7 +
  8 +import numpy
  9 +
  10 +class Vertex:
  11 + def __init__(self, f): #load a vertex given a file ID at that vertex location
  12 + self.p = numpy.fromfile(f, dtype=numpy.float32, count=3) #load the vertex coordinates
  13 +
  14 + EO = int.from_bytes(f.read(4), "little") #get the number of edges moving out of and into this vertex
  15 + EI = int.from_bytes(f.read(4), "little")
  16 +
  17 + self.Eout = numpy.fromfile(f, dtype=numpy.uint32, count=EO) #get the incoming and outgoing edge lists
  18 + self.Ein = numpy.fromfile(f, dtype=numpy.uint32, count=EI)
  19 +
  20 + def __str__(self):
  21 + sp = "(" + str(self.p[0]) + ", " + str(self.p[1]) + ", " + str(self.p[2]) + ")"
  22 +
  23 + seo = "["
  24 + for e in self.Eout:
  25 + seo = seo + " " + str(e)
  26 + seo = seo + " ]"
  27 +
  28 + sei = "["
  29 + for e in self.Ein:
  30 + sei = sei + " " + str(e)
  31 + sei = sei + " ]"
  32 +
  33 + return sp + " out = " + seo + " in = " + sei
  34 +
  35 +class Edge:
  36 + def __init__(self, f):
  37 + self.v = numpy.fromfile(f, dtype=numpy.uint32, count=2)
  38 + P = int.from_bytes(f.read(4), "little")
  39 + p1d = numpy.fromfile(f, dtype=numpy.float32, count=4*P)
  40 + self.p = numpy.reshape(p1d, (P, 4))
  41 +
  42 + def __str__(self):
  43 + sv = "[ " + str(self.v[0]) + "---> " + str(self.v[1]) + " ]"
  44 + return sv + str(self.p)
  45 +
  46 +class NWT:
  47 + def __init__(self, fname):
  48 + f = open(fname, "rb") #open the NWT file for binary reading
  49 + self.filetype = f.read(14) #this should be "nwtfileformat "
  50 + self.desc = f.read(58) #NWT file description
  51 + V = int.from_bytes(f.read(4), "little") #retrieve the number of vertices (graph)
  52 + E = int.from_bytes(f.read(4), "little")
  53 +
  54 + self.v = []
  55 + for i in range(0, V):
  56 + self.v.append(Vertex(f))
  57 +
  58 + self.e = []
  59 + for i in range(1, E):
  60 + self.e.append(Edge(f))
0 61 \ No newline at end of file
... ...