diff --git a/python/fibernet.py b/python/fibernet.py new file mode 100644 index 0000000..36512dd --- /dev/null +++ b/python/fibernet.py @@ -0,0 +1,142 @@ +# -*- coding: utf-8 -*- +""" +Created on Mon Mar 19 14:45:39 2018 + +@author: david +""" + +import nwt +import matplotlib.pyplot as plt +import numpy + + +class fibernet: + + def __init__(self, filename): + self.nwt(filename) + self.buildgraph() + + #load a network from a NWT file + def nwt(self, filename): + net = nwt.NWT(filename) + + #first insert the points where fibers are connected + self.P = [] + for vi in range(0, len(net.v)): #store each vertex in the point list + self.P.append( (net.v[vi].p[0], net.v[vi].p[1], net.v[vi].p[2]) ) + + #now insert each fiber, connecting them appropriately using previously inserted points (vertices) + self.F = [] + for ei in range(0, len(net.e)): + f = [numpy.int(net.e[ei].v[0])] #create a new fiber to store point indices, initialize with the first vertex ID + for pi in range(1, net.e[ei].p.shape[0]): + 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 + f.append(len(self.P) - 1) #append the current point ID + f.append(numpy.int(net.e[ei].v[1])) #append the last vertex ID + self.F.append(f) + + #build a graph from the available geometry + def buildgraph(self): + p_to_v = numpy.ones((self.npoints()), numpy.int) * (-1) #create an array that maps a point ID to a vertex ID + + self.V = [] #create an empty vertex list + self.E = [] #create an empty edge list + for fi in range(0, self.nfibers()): #for each fiber in the network + if p_to_v[self.F[fi][0]] == -1: #if the first point hasn't been encountered + p_to_v[self.F[fi][0]] = len(self.V) #add it to the mapping + self.V.append( (self.F[fi][0], []) ) #add a new vertex representing this point + if p_to_v[self.F[fi][-1]] == -1: + p_to_v[self.F[fi][-1]] = len(self.V) + self.V.append( (self.F[fi][-1], []) ) + + v0 = p_to_v[self.F[fi][0]] #get the two vertex indices corresponding to this edge + v1 = p_to_v[self.F[fi][-1]] + + self.E.append( (v0, v1) ) #create an edge representing this fiber and add it to the edge list + self.V[v0][1].append(fi) #add this edge ID to the corresponding vertices to update the adjacency list + self.V[v1][1].append(fi) + + #return the length of a fiber specified by fiber ID fid + def length(self, fid): + p = self.points(fid) #get the points corresponding to the fiber + l = 0 #initialize the length to zero + for pi in range(1, len(p)): #for each line segment + v = p[pi] - p[pi-1] #calculate the vector representing the line segment + l = l + numpy.linalg.norm(numpy.array(v)) #add the length of the line segment + return l #return the total length + + #return a list of vertices adjacent to the vertex specified by vid + def adjacent(self, vid): + adj = [] #initialize an empty adjacency list + for e in self.V[vid][1]: #for each edge connected to vid + if self.E[e][0] != vid: + adj.append(self.E[e][0]) + else: + adj.append(self.E[e][1]) + return adj + + #return the fiber network graph as an adjacency list + def adjacencylist(self): + A = [] + for vi in range(0, len(self.V)): #for each vertex in the graph + A.append( self.adjacent(vi) ) #add the adjacency list for the vertex into the main adjacency list + return A + + #return the fiber network graph as an adjacency matrix + # possible values for the matrix include: + # binary - 1 when there is a connection, 0 otherwise + # edgeID - stores the edge ID corresponding to a connection, and -1 otherwise + # length - stores the length of the connecting fiber, and -1 if there is no connection + def adjacencymatrix(self, value="binary"): + if value == "binary": + M = numpy.ones( (len(self.V), len(self.V)), numpy.bool ) + elif value == "edgeID": + M = numpy.ones( (len(self.V), len(self.V)), numpy.int ) * -1 + elif value == "length": + M = numpy.ones( (len(self.V), len(self.V)) ) * -1 + else: + return [] + + for vi in range(0, len(self.V)): #for each vertex in the graph + A = self.adjacent(vi) #get the adjacency matrix + for a in A: + if value == "binary": + M[vi,a] = M[a,vi] = 1 + if value == "edgeID": + M[vi,a] = M[a,vi] = self.findedge(a, vi) + if value == "length": + M[vi,a] = M[a,vi] = self.length(self.findedge(a,vi)) + return M + + #find the edge ID corresponding to two vertices + def findedge(self, v0, v1): + + edges = self.V[v0][1] #get the list of candidate edges connected to v0 + for e in edges: #for each edge + if self.E[e][0] == v1 or self.E[e][1] == v1: + return e + + return -1 #return -1 if the edge doesn't exist + + #get the points associated with a specified fiber ID fid + def points(self, fid): #get the points corresponding to a specified fiber index + f = self.F[fid] #get the list of point indices + return numpy.array(self.P)[f] + + #get the total number of points in the network + def npoints(self): + return len(self.P) + + #get the total number of fibers in the network + def nfibers(self): + return len(self.F) + + #plot the network as a 3D line plot + def plot(self): + fig = plt.figure() + ax = fig.add_subplot(111, projection='3d') + + for fi in range(0, self.nfibers()): + test = fibers.points(fi) + ax.plot(test[:, 0], test[:, 1], test[:, 2]) + return fig \ No newline at end of file diff --git a/python/nwt.py b/python/nwt.py new file mode 100644 index 0000000..0203075 --- /dev/null +++ b/python/nwt.py @@ -0,0 +1,60 @@ +# -*- coding: utf-8 -*- +""" +Created on Mon Mar 19 12:38:30 2018 + +@author: david +""" + +import numpy + +class Vertex: + def __init__(self, f): #load a vertex given a file ID at that vertex location + self.p = numpy.fromfile(f, dtype=numpy.float32, count=3) #load the vertex coordinates + + EO = int.from_bytes(f.read(4), "little") #get the number of edges moving out of and into this vertex + EI = int.from_bytes(f.read(4), "little") + + self.Eout = numpy.fromfile(f, dtype=numpy.uint32, count=EO) #get the incoming and outgoing edge lists + self.Ein = numpy.fromfile(f, dtype=numpy.uint32, count=EI) + + def __str__(self): + sp = "(" + str(self.p[0]) + ", " + str(self.p[1]) + ", " + str(self.p[2]) + ")" + + seo = "[" + for e in self.Eout: + seo = seo + " " + str(e) + seo = seo + " ]" + + sei = "[" + for e in self.Ein: + sei = sei + " " + str(e) + sei = sei + " ]" + + return sp + " out = " + seo + " in = " + sei + +class Edge: + def __init__(self, f): + self.v = numpy.fromfile(f, dtype=numpy.uint32, count=2) + P = int.from_bytes(f.read(4), "little") + p1d = numpy.fromfile(f, dtype=numpy.float32, count=4*P) + self.p = numpy.reshape(p1d, (P, 4)) + + def __str__(self): + sv = "[ " + str(self.v[0]) + "---> " + str(self.v[1]) + " ]" + return sv + str(self.p) + +class NWT: + def __init__(self, fname): + f = open(fname, "rb") #open the NWT file for binary reading + self.filetype = f.read(14) #this should be "nwtfileformat " + self.desc = f.read(58) #NWT file description + V = int.from_bytes(f.read(4), "little") #retrieve the number of vertices (graph) + E = int.from_bytes(f.read(4), "little") + + self.v = [] + for i in range(0, V): + self.v.append(Vertex(f)) + + self.e = [] + for i in range(1, E): + self.e.append(Edge(f)) \ No newline at end of file -- libgit2 0.21.4