From db331d8a99a1ff623a601d4ef5b40fb12019b890 Mon Sep 17 00:00:00 2001 From: David Mayerich Date: Thu, 4 Jan 2018 19:42:04 -0600 Subject: [PATCH] added adjacency weight calculations and scaling functions to Network, updated spharmonics functions to take tuples --- python/network.py | 170 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++----------------------------------------------- python/spharmonics.py | 48 ++++++++++++++++++++++++++++++------------------ 2 files changed, 153 insertions(+), 65 deletions(-) diff --git a/python/network.py b/python/network.py index a02cd7b..584e3ad 100644 --- a/python/network.py +++ b/python/network.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(' 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 + + diff --git a/python/spharmonics.py b/python/spharmonics.py index 9093511..fb25eef 100644 --- a/python/spharmonics.py +++ b/python/spharmonics.py @@ -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') -- libgit2 0.21.4