added adjacency weight calculations and scaling functions to Network, updated sp…

…harmonics functions to take tuples

added adjacency weight calculations and scaling functions to Network, updated sp…
…harmonics functions to take tuples
David Mayerich
1 parent 8d96a467
Showing 2 changed files with 153 additions and 65 deletions Show diff stats
python/network.py
python/spharmonics.py
@@ -11,6 +11,8 @@ import scipy as sp
 import networkx as nx
 import matplotlib.pyplot as plt
 import math
+import time
+import spharmonics
 '''
     Definition of the Node class
@@ -33,18 +35,12 @@ class Node:
 '''
 class Fiber:
-    def __init__ (self):
-        self.v0 = 0
-        self.v1 = 0
-        self.points = []
-        self.radii = []
-        
-        #NOTE: there is no function overloading in Python
-#    def __init__ (self, p1, p2, pois, rads):
-#        self.v0 = p1
-#        self.v1 = p2
-#        self.points = pois
-#        self.radii = rads
+       
+    def __init__ (self, p1, p2, pois, rads):
+        self.v0 = p1
+        self.v1 = p2
+        self.points = pois
+        self.radii = rads
     '''
         return the length of the fiber.
@@ -83,7 +79,7 @@ class NWT:
     '''
         Writes the header given and open file descripion, number of verticies and number of edges.
     '''
-    def writeHeader(self, open_file, numVerts, numEdges):
+    def writeHeader(open_file, numVerts, numEdges):
         txt = "nwtFileFormat fileid(14B), desc(58B), #vertices(4B), #edges(4B): bindata"
         b = bytearray()
         b.extend(txt.encode())
@@ -95,7 +91,7 @@ class NWT:
     '''
         Writes a single vertex to a file.
     '''
-    def writeVertex(self, open_file, vertex):
+    def writeVertex(open_file, vertex):
         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]))
@@ -112,7 +108,7 @@ class NWT:
     '''
         Writes a single fiber to a file.
     '''
-    def writeFiber(self, open_file, edge):
+    def writeFiber(open_file, edge):
         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)))
@@ -127,14 +123,14 @@ class NWT:
     '''
         Writes the entire network to a file in str given the vertices array and the edges array.
     '''
-    def exportNWT(self, str, vertices, edges):
+    def exportNWT(str, vertices, edges):
         with open(str, "wb") as file:
-            self.writeHeader(file, len(vertices), len(edges))
+            NWT.writeHeader(file, len(vertices), len(edges))
             for i in range(len(vertices)):
-                self.writeVertex(file, vertices[i])
+                NWT.writeVertex(file, vertices[i])
             for i in range(len(edges)):
-                self.writeFiber(file, edges[i])
+                NWT.writeFiber(file, edges[i])
         return
@@ -142,7 +138,7 @@ class NWT:
     '''
         Reads a single vertex from an open file and returns a node Object.
     '''
-    def readVertex(self, open_file):
+    def readVertex(open_file):
         points = np.tile(0., 3)
         bytes = open_file.read(4)
         points[0] = struct.unpack('f', bytes)[0]
@@ -172,7 +168,7 @@ class NWT:
     '''
         Reads a single fiber from an open file and returns a Fiber object .   
     '''
-    def readFiber(self, open_file):
+    def readFiber(open_file):
         bytes = open_file.read(4)
         vtx0 = int.from_bytes(bytes, byteorder = 'little')
         bytes = open_file.read(4)
@@ -206,8 +202,11 @@ class NWT:
 class Network:
-    def __init__(self, str):
-        with open(str, "rb") as file:
+    def __init__(self, filename, clock=False):
+        if clock:
+            start_time = time.time()
+            
+        with open(filename, "rb") as file:
             header = file.read(72)
             bytes = file.read(4)
             numVertex = int.from_bytes(bytes, byteorder='little')
@@ -223,6 +222,8 @@ class Network:
             for i in range(numEdges):
                 edge = NWT.readFiber(file)
                 self.F.append(edge)
+        if clock:
+            print("Network initialization: "  + str(time.time() - start_time) + "s")
     '''
     Creates a graph from a list of nodes and a list of edges.
@@ -287,7 +288,7 @@ class Network:
         lower = self.N[0].p.copy()
         upper = lower.copy()
-        for i in self.nodeList:
+        for i in self.N:
             for c in range(len(lower)):
                 if lower[c] > i.p[c]:
                     lower[c] = i.p[c]
@@ -317,7 +318,7 @@ class Network:
         R = (200, 200, 200)
         #generate a meshgrid of the appropriate size and resolution to surround the network
-        lower, upper = self.aabb(self.N, self.F)           #get the space occupied by the network
+        lower, upper = self.aabb(self.N, self.F)    #get the space occupied by the network
         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])
@@ -333,31 +334,106 @@ class Network:
         return D
     #returns the number of points in the network
-    def npoints(self):
-        n = 0
-        for f in self.F:
-            n = n + len(f.points) - 2
-        n = n + len(self.N)
-        return n
+    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        
     #returns the number of linear segments in the network
     def nsegments(self):
-        n = 0
-        for f in self.F:
-            n = n + len(f.points) - 1
-        return n
+        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
-    def vectorize(self):
-        #initialize three coordinate vectors ()
-        n = self.nsegments(self.N, self.F)
-        X = np.zeros((n))
-        Y = np.zeros((n))
-        Z = np.zeros((n))
+    #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
+        
         idx = 0
-        for i in range(0, len(self.F)):
-            for j in range(0, len(self.F[i].points)-1):
-                X[idx] = self.F[i].points[j][0]-self.F[i].points[j+1][0]
-                Y[idx] = self.F[i].points[j][1]-self.F[i].points[j+1][1]
-                Z[idx] = self.F[i].points[j][2]-self.F[i].points[j+1][2]
+        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
                 idx = idx + 1
-        return X, Y, Z
+                
+        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
+        
+
@@ -11,16 +11,17 @@ import matplotlib.pyplot as plt
 from matplotlib import cm, colors
 from mpl_toolkits.mplot3d import Axes3D
 import math
+import time
-def sph2cart(theta, phi, r):
+def sph2cart(r, theta, phi):
     x = r * numpy.cos(theta) * numpy.sin(phi)
     y = r * numpy.sin(theta) * numpy.sin(phi)
     z = r * numpy.cos(phi)
     return x, y, z
-def cart2sph(x,y,z):
+def cart2sph(x, y, z):
     r = numpy.sqrt(x**2+y**2+z**2)
     theta = numpy.arctan2(y,x)
     phi = numpy.arccos(z/r)
@@ -41,15 +42,15 @@ def sh(theta, phi, l, m):
     else:
         return scipy.special.sph_harm(0, l, theta, phi).real
-#calculate a spherical harmonic value from a set of coefficients and a spherical coordinate        
-def sh_coeff(theta, phi, C):
+#calculate a spherical harmonic value from a set of coefficients and coordinates P = (theta, phi)    
+def sh_coeff(P, C):
-    s = numpy.zeros(theta.shape)
+    s = numpy.zeros(P[0].shape)
     c = range(len(C))
     for c in range(len(C)):
         l, m = i2lm(c)                  #get the 2D indices
-        s = s + C[c] * sh(theta, phi, l, m)
+        s = s + C[c] * sh(P[0], P[1], l, m)
     return s
@@ -86,34 +87,45 @@ def lm2i(l, m):
     return l * (l+1) + m
 #generates a set of spherical harmonic coefficients from samples using linear least squares
-def linfit(theta, phi, s, nc):
-    #allocate space for the matrix
+def linfit(P, s, nc, clock=False):
+    if clock:
+        start_time = time.time()
+        
+    #allocate space for the matrix and RHS values
     A = numpy.zeros((nc, nc))
+    b = numpy.zeros(nc)
     #calculate each of the matrix coefficients
         #(see SH technical report in the vascular_viz repository)
     for i in range(nc):
         li, mi = i2lm(i)
-        yi = sh(theta, phi, li, mi)
+        yi = sh(P[0], P[1], li, mi)
         for j in range(nc):        
             lj, mj = i2lm(j)
-            yj = sh(theta, phi, lj, mj)
+            yj = sh(P[0], P[1], lj, mj)
             A[i, j] = numpy.sum(yi * yj)
+        b[i] = numpy.sum(yi * s)            #calculate the RHS value
     #calculate the RHS values
-    b = numpy.zeros(nc)
-    for j in range(nc):
-        lj, mj = i2lm(j)
-        yj = sh(theta, phi, lj, mj)
-        b[j] = numpy.sum(yj * s)
+    #for j in range(nc):
+    #    lj, mj = i2lm(j)
+    #    yj = sh(theta, phi, lj, mj)
+    #    b[j] = numpy.sum(yj * s)
+    if clock:
+        print("SH::linfit:matrix "+str(time.time() - start_time)+"s")
     #solve the system of linear equations
-    return numpy.linalg.solve(A, b)
+    R = numpy.linalg.solve(A, b)
+    
+    if clock:
+        print("SH::linfit:solution "+str(time.time() - start_time)+"s")
+    return R
 #generate a scatter plot in 3D using spherical coordinates
-def scatterplot3d(theta, phi, r):
+def scatterplot3d(P):
+    r, theta, phi = P
     #convert all of the samples to cartesian coordinates
-    X, Y, Z = sph2cart(theta, phi, r)
+    X, Y, Z = sph2cart(r, theta, phi)
     fig = plt.figure()
     ax = fig.add_subplot(111, projection='3d')