# -*- coding: utf-8 -*- """ Created on Sat Sep 16 16:34:49 2017 @author: pavel """ import struct import numpy as np from scipy.sparse.linalg import eigsh import scipy as sp from sklearn.cluster import SpectralClustering #import networkx as nx import matplotlib.pyplot as plt from matplotlib import cm import math import time #import spharmonics import graph_tool.all as gt import copy #import matplotlib #for testing #import matplotlib.mlab as mlab #import matplotlib.pyplot as plt #from matplotlib import cm DEBUG = False ''' 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 # def p(): # return self.p ''' 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: def __init__ (self, p1, p2, pois, rads): self.v0 = p1 self.v1 = p2 self.points = pois self.radii = rads ''' return the array-likes of the x,y,z,r coordinates of the fiber for the gt representation ''' def getcoords(self): x = np.zeros(len(self.points), dtype=np.double) y = np.zeros(len(self.points), dtype=np.double) z = np.zeros(len(self.points), dtype=np.double) r = np.zeros(len(self.points), dtype=np.double) for i in range(len(self.points)): x[i] = self.points[i][0] y[i] = self.points[i][1] z[i] = self.points[i][2] r[i] = self.radii[i] return x,y,z,r ''' 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)) if(length == 0): print("NON-CRITICAL ERROR: edge with length 0 is detected: IDX = "\ , i, ", points = ", self.points, " len(points) = ", len(self.points)\ , " vertices = ", self.v0, " ", self.v1) return length ''' returns the turtuosity of the fiber. ''' def turtuosity(self): turtuosity = 1.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)) if(distance > 0): 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]*10e-6,2) + math.pow(self.radii[i+1]*10e-6,2) + self.radii[i]*10e-6*self.radii[i+1]*10e-6) * 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 def av_radius(self): radius = 0.0 for i in range(len(self.radii)): radius = radius + self.radii[i] return radius/len(self.radii) class NWT: ''' Writes the header given and open file descripion, number of verticies and number of edges. ''' def writeHeader(open_file, numVerts, numEdges): 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)) ''' Writes a single vertex to a file. ''' def writeVertex(open_file, vertex): open_file.write(struct.pack(' self.B[0]: self.B[0] = G.vertex_properties["p"][v][0] if G.vertex_properties["p"][v][1] < self.A[1]: self.A[1] = G.vertex_properties["p"][v][1] if G.vertex_properties["p"][v][1] > self.B[1]: self.B[1] = G.vertex_properties["p"][v][1] if G.vertex_properties["p"][v][2] < self.A[2]: self.A[2] = G.vertex_properties["p"][v][2] if G.vertex_properties["p"][v][2] > self.B[2]: self.B[2] = G.vertex_properties["p"][v][2] #In case of a dual graph we have to scane the first and last points # of every fiber to find the true bounding box else: for v in G.vertices(): l = len(G.vertex_properties["x"][v]) if G.vertex_properties["x"][v][0] < self.A[0]: self.A[0] = G.vertex_properties["x"][v][0] if G.vertex_properties["x"][v][l-1] < self.A[0]: self.A[0] = G.vertex_properties["x"][v][l-1] if G.vertex_properties["y"][v][0] < self.A[1]: self.A[1] = G.vertex_properties["y"][v][0] if G.vertex_properties["y"][v][l-1] < self.A[1]: self.A[1] = G.vertex_properties["y"][v][l-1] if G.vertex_properties["z"][v][0] < self.A[2]: self.A[2] = G.vertex_properties["z"][v][0] if G.vertex_properties["z"][v][l-1] < self.A[2]: self.A[2] = G.vertex_properties["z"][v][l-1] if G.vertex_properties["x"][v][0] > self.B[0]: self.B[0] = G.vertex_properties["x"][v][0] if G.vertex_properties["x"][v][l-1] > self.B[0]: self.B[0] = G.vertex_properties["x"][v][l-1] if G.vertex_properties["y"][v][0] > self.B[1]: self.B[1] = G.vertex_properties["y"][v][0] if G.vertex_properties["y"][v][l-1] > self.B[1]: self.B[1] = G.vertex_properties["y"][v][l-1] if G.vertex_properties["z"][v][0] > self.B[2]: self.B[2] = G.vertex_properties["z"][v][0] if G.vertex_properties["z"][v][l-1] > self.B[2]: self.B[2] = G.vertex_properties["z"][v][l-1] #print(self.A, self.B) self.O = (self.A+self.B)*0.5 self.planes = [] #append x planes described by a point and a vector self.planes.append((np.array([self.A[0], 0.0, 0.0]), np.array([1.0, 0.0, 0.0]))) self.planes.append((np.array([self.B[0], 0.0, 0.0]), np.array([-1.0, 0.0, 0.0]))) #append y planes described by a point and a vector self.planes.append((np.array([0.0, self.A[1], 0.0]), np.array([0.0, 1.0, 0.0]))) self.planes.append((np.array([0.0, self.B[1], 0.0]), np.array([0.0, -1.0, 0.0]))) #append z planes described by a point and a vector self.planes.append((np.array([0.0, 0.0, self.A[2]]), np.array([0.0, 0.0, 1.0]))) self.planes.append((np.array([0.0, 0.0, self.B[2]]), np.array([0.0, 0.0, -1.0]))) def distance(self, pt): dist = 10000000 for i in self.planes: V = i[0]-pt if np.dot(V, i[1]) < dist: dist = np.dot(V, i[1]) return dist def project_grid(self, n): #r = abs(self.A - self.B) x = np.linspace(self.A[0], self.B[0], (n+2)) x = x[1:-1] y = np.linspace(self.A[1], self.B[1], (n+2)) y = y[1:-1] z = np.linspace(self.A[2], self.B[2], (n+2)) z = z[1:-1] return x,y,z #returns a resampled bounding both with nxn points on each side. def resample_sides(self, n): #get the size of the cube in every cardinal direction size = abs(self.A - self.B) #set the stepsize for each subdivided point dist_x = size[0]/(n-1) dist_y = size[1]/(n-1) dist_z = size[2]/(n-1) #generate the original 8 corners vertices = [] #generate the points for the yz planes. for i in range(n): for j in range(n): #generate the temporary vectors for both the planes temp_p1 = copy.deepcopy(self.A) temp_p2 = copy.deepcopy(self.A) temp_p2[0] = self.B[0] temp_p1[1] = temp_p1[1] + dist_y*i temp_p1[2] = temp_p1[2] + dist_z*j temp_p2[1] = temp_p2[1] + dist_y*i temp_p2[2] = temp_p2[2] + dist_z*j vertices.append(temp_p1) vertices.append(temp_p2) #generate the points for the xz planes. for i in range(n): for j in range(n): #generate the temporary vectors for both the planes temp_p1 = copy.deepcopy(self.A) temp_p2 = copy.deepcopy(self.A) temp_p2[1] = self.B[1] temp_p1[0] = temp_p1[0] + dist_x*i temp_p1[2] = temp_p1[2] + dist_z*j temp_p2[0] = temp_p2[0] + dist_x*i temp_p2[2] = temp_p2[2] + dist_z*j vertices.append(temp_p1) vertices.append(temp_p2) #generate the points for the xy planes. for i in range(n): for j in range(n): #generate the temporary vectors for both the planes temp_p1 = copy.deepcopy(self.A) temp_p2 = copy.deepcopy(self.A) temp_p2[2] = self.B[2] temp_p1[0] = temp_p1[0] + dist_x*i temp_p1[1] = temp_p1[1] + dist_y*j temp_p2[0] = temp_p2[0] + dist_x*i temp_p2[1] = temp_p2[1] + dist_y*j vertices.append(temp_p1) vertices.append(temp_p2) vertices = list(np.unique(np.array(vertices), axis=0)) #import matplotlib.pyplot as plt #fig = plt.figure() #ax = plt.axes(projection='3d') #ax.scatter3D(np.array(vertices)[:, 0], np.array(vertices)[:, 1], np.array(vertices)[:, 2]) #print("THERE ARE THIS MANY ", len(vertices)) return vertices def getVolume(self): size = abs(self.A - self.B) return size[0]*size[1]*size[2] def vertices(self): size = abs(self.A - self.B) points = [] points.append(self.A) temp = copy.deepcopy(self.A) temp[0] = temp[0] + size[0] points.append(temp) if DEBUG: print('1', temp) temp = copy.deepcopy(self.A) temp[1] = temp[1] + size[1] points.append(temp) if DEBUG: print('1', temp) temp = copy.deepcopy(self.A) temp[2] = temp[2] + size[2] points.append(temp) if DEBUG: print('1', temp) temp = copy.deepcopy(self.A) temp[1] = temp[1] + size[1] temp[0] = temp[0] + size[0] points.append(temp) if DEBUG: print('1', temp) temp = copy.deepcopy(self.A) temp[2] = temp[2] + size[2] temp[0] = temp[0] + size[0] points.append(temp) if DEBUG: print('1', temp) temp = copy.deepcopy(self.A) temp[1] = temp[1] + size[1] temp[2] = temp[2] + size[2] points.append(temp) if DEBUG: print('1', temp) temp = copy.deepcopy(self.A) temp[0] = temp[0] + size[0] temp[1] = temp[1] + size[1] temp[2] = temp[2] + size[2] points.append(temp) if DEBUG: print('1', temp) return points ''' Class defining the BFS traversal of all the vertices in the graph and labeling all interconnected vertices. ''' class VisitorClassDisconnected(gt.BFSVisitor): def __init__(self, ecluster, cluster, c): self.ecluster = ecluster self.cluster = cluster self.c = c def examine_vertex(self, u): if(self.cluster[u] == -1): self.cluster[u] = self.c def examine_edge(self, e): if(self.ecluster[e] == -1): self.ecluster[e] = self.c ''' Class defining the BFS traversal of all the vertices in the graph and labeling all interconnected vertices. ''' class VisitorClassPartition(gt.BFSVisitor): def __init__(self, cluster, dist, c, iterations): self.cluster = cluster self.dist = dist self.c = c self.total_iter = iterations def examine_vertex(self, u): if(self.cluster[u] == -1): self.cluster[u] = self.c def tree_edge(self, e): d = self.dist[e.source()] #print(d, self.total_iter) if d <= self.total_iter: self.dist[e.target()] = self.dist[e.source()] + 1 class Network: 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') bytes = file.read(4) numEdges = int.from_bytes(bytes, byteorder='little') self.N = [] self.F = [] for i in range(numVertex): node = NWT.readVertex(file) self.N.append(node) for i in range(numEdges): edge = NWT.readFiber(file) self.F.append(edge) if clock: print("Network initialization: " + str(time.time() - start_time) + "s") ''' Takes in a graph object (networkx Graph) and saves it as a network file Resulting nwt, does not contain incomind or outgoing fibers ''' def saveGraph(self, G, path): G_nodes = list(G.nodes()) G_edges = list(G.edges()) Nodes = [] Edges = [] points = list(nx.get_node_attributes(G, 'p').values()) for i in range(len(G_nodes)): node = Node(points[i], [], []) Nodes.append(node) for e in G_edges: edge = Fiber(e[0], e[1], G[e[0]][e[1]]['pts'], G[e[0]][e[1]]['rads']) Edges.append(edge) NWT.exportNWT(path, Nodes, Edges) ''' Saves the graph-tool graph as a nwt file with all the variables. ''' def saveGraph_gt(self, G, path): points = G.vertex_properties["p"] Nodes = [] Edges = [] for i in range(G.num_vertices()): #array = G.get_out_edges(i) outgoing = [] for edge in G.get_out_edges(i): outgoing.append(edge[2]) node = Node(points[i], outgoing, outgoing) Nodes.append(node) for e in G.edges(): x = G.edge_properties["x"][e].get_array() y = G.edge_properties["y"][e].get_array() z = G.edge_properties["z"][e].get_array() r = G.edge_properties["r"][e].get_array() pts = [] for i in range(len(x)): pt = [] pt.append(x[i]) pt.append(y[i]) pt.append(z[i]) pts.append(pt) #pts = np.ndarray.tolist(np.dstack([x,y,z])) #print(pts) edge = Fiber(int(e.source()), int(e.target()), pts, np.ndarray.tolist(r)) if(edge.length() > 0.0): Edges.append(edge) NWT.exportNWT(path, Nodes, Edges) def saveGraphStatistics_gt(self, G, path, prefix, label = "none", norm=False): location = path + "/" + prefix + "_edges.txt" f = open(location, "a+") f.write("inverted\t") f.write("length\t") f.write("tortuosity\t") f.write("volume\t") f.write("inverted_volume\t") f.write("gaussian\t") f.write("bc\t") f.write("bc*length\t") f.write("mst\t\n") for e in G.edges(): f.write("%.15f\t" % G.edge_properties["inverted"][e]) f.write("%.15f\t" % G.edge_properties["length"][e]) f.write("%.15f\t" % G.edge_properties["tortuosity"][e]) f.write("%.15f\t" % G.edge_properties["volume"][e]) f.write("%.15f\t" % G.edge_properties["inverted_volume"][e]) f.write("%.15f\t" % G.edge_properties["gaussian"][e]) f.write("%.15f\t" % G.edge_properties["bc"][e]) f.write("%.15f\t" % G.edge_properties["bc*length"][e]) f.write("%.15f\t" % G.edge_properties["mst"][e]) f.write("%s\n" % label) f.close() location = path+"/" + prefix + "_vertices.txt" f = open(location, "a+") f.write("bc\t") f.write("degree\t") f.write("degree_volume\t") f.write("degree_tortuosity\t\n") for v in G.vertices(): f.write("%.15f\t" % G.vertex_properties["bc"][v]) f.write("%.15f\t" % G.vertex_properties["degree"][v]) f.write("%.15f\t" % G.vertex_properties["degree_volume"][v]) f.write("%.15f\t" % G.vertex_properties["degree_tortuosity"][v]) f.write("%s\n" % label) ''' Saves the graph as a sample. Additionally saves the key as a list of string if writeKey is True saves the file at path with prefix_s.txt ''' def saveSample_gt(self, G, path, prefix, is_dual, label = "none", norm=False, writeKey=False): location = path + "/" + prefix + "s.txt" f = open(location, "a+") aabb = AABB(G, is_dual) array = G.vertex_properties["bc"].get_array() av_array = np.sum(array)/G.num_vertices() G.graph_properties["global_vertex_bc"] = G.new_graph_property("double", val = av_array) f.write("%.15f\t" % av_array) array = G.vertex_properties["bc"].get_array() av_array = np.sum(array)/aabb.getVolume() G.graph_properties["global_vascular_vertex_bc"] = G.new_graph_property("double", val = av_array) f.write("%.15f\t" % av_array) array = G.edge_properties["bc"].get_array() av_array = np.sum(array)/aabb.getVolume() G.graph_properties["global_vascular_edge_bc"] = G.new_graph_property("double", val = av_array) f.write("%.15f\t" % av_array) #G.edge_properties["mst"] = gt.graph_tool.topology.min_spanning_tree(G) G.graph_properties["global_mst_ratio"] = G.new_graph_property("double", val = np.double(np.sum(G.edge_properties["mst"].get_array()))/np.double(G.num_edges())) f.write("%.15f\t" % G.graph_properties["global_mst_ratio"]) #G.edge_properties["mst"] = gt.graph_tool.topology.min_spanning_tree(G, weights=G.edge_properties["volume"]) G.graph_properties["global_mst_ratio_vascular_volume"] = G.new_graph_property("double", val = np.double(np.sum(G.edge_properties["mst"].get_array()))/np.double(G.num_edges())/aabb.getVolume()) f.write("%.15f\t" % G.graph_properties["global_mst_ratio"]) ########################Local Statistics############################### self.recalculate_metrics(G, is_dual, norm) ########################Biological Statistics########################## # Some biological statistics are reported as a ratio to physical volume of vessels array = G.edge_properties["length"].get_array() av_array = np.sum(array)/G.num_edges() G.graph_properties["average_length"] = G.new_graph_property("double", val= av_array) f.write("%.15f\t" % av_array) array = G.edge_properties["volume"].get_array() av_array = np.sum(array)/G.num_edges() G.graph_properties["average_volume"] = G.new_graph_property("double", val= av_array) f.write("%.15f\t" % av_array) array = G.edge_properties["tortuosity"].get_array() av_array = np.sum(array)/G.num_edges() G.graph_properties["average_tortuosity"] = G.new_graph_property("double", val= av_array) f.write("%.15f\t" % av_array) array = G.edge_properties["length"].get_array() av_array = np.sum(array)/aabb.getVolume() G.graph_properties["vascular length"] = G.new_graph_property("double", val= av_array) f.write("%.15f\t" % av_array) array = G.edge_properties["volume"].get_array() av_array = np.sum(array)/aabb.getVolume() G.graph_properties["vascular_volume"] = G.new_graph_property("double", val= av_array) f.write("%.15f\t" % av_array) array = G.edge_properties["tortuosity"].get_array() av_array = np.sum(array)/aabb.getVolume() G.graph_properties["vascular_tortuosity"] = G.new_graph_property("double", val= av_array) f.write("%.15f\t" % av_array) ############################Graph Statistics########################### # Some Graph Statistics are reported as values representing the graph # Some statistics are reported as a function of volume. array = G.vertex_properties["degree"].get_array() av_array = np.sum(array)/G.num_vertices() G.graph_properties["local_average_degree"] = G.new_graph_property("double", val = av_array) f.write("%.15f\t" % av_array) array = G.vertex_properties["degree_volume"].get_array() av_array = np.sum(array)/G.num_vertices() G.graph_properties["local_average_degree_volume"] = G.new_graph_property("double", val = av_array) f.write("%.15f\t" % av_array) array = G.vertex_properties["degree_tortuosity"].get_array() av_array = np.sum(array)/G.num_vertices() G.graph_properties["local_average_degree_tortuosity"] = G.new_graph_property("double", val = av_array) f.write("%.15f\t" % av_array) array = G.vertex_properties["bc"].get_array() av_array = np.sum(array)/G.num_vertices() G.graph_properties["local_vertex_bc"] = G.new_graph_property("double", val = av_array) f.write("%.15f\t" % av_array) array = G.vertex_properties["degree"].get_array() av_array = np.sum(array)/aabb.getVolume() G.graph_properties["local_vascular_degree"] = G.new_graph_property("double", val = av_array) f.write("%.15f\t" % av_array) array = G.vertex_properties["degree_volume"].get_array() av_array = np.sum(array)/aabb.getVolume() G.graph_properties["local_vascular_degree_volume"] = G.new_graph_property("double", val = av_array) f.write("%.15f\t" % av_array) array = G.vertex_properties["degree_tortuosity"].get_array() av_array = np.sum(array)/aabb.getVolume() G.graph_properties["local_vascular_degree_tortuosity"] = G.new_graph_property("double", val = av_array) f.write("%.15f\t" % av_array) array = G.vertex_properties["bc"].get_array() av_array = np.sum(array)/aabb.getVolume() G.graph_properties["local_vascular_vertex_bc"] = G.new_graph_property("double", val = av_array) f.write("%.15f\t" % av_array) array = G.edge_properties["bc"].get_array() av_array = np.sum(array)/aabb.getVolume() G.graph_properties["local_vascular_edge_bc"] = G.new_graph_property("double", val = av_array) f.write("%.15f\t" % av_array) G.edge_properties["mst"] = gt.graph_tool.topology.min_spanning_tree(G) G.graph_properties["local_mst_ratio"] = G.new_graph_property("double", val = np.double(np.sum(G.edge_properties["mst"].get_array()))/np.double(G.num_edges())) f.write("%.15f\t" % G.graph_properties["local_mst_ratio"]) G.edge_properties["mst"] = gt.graph_tool.topology.min_spanning_tree(G, weights=G.edge_properties["volume"]) G.graph_properties["local_mst_ratio_vascular_volume"] = G.new_graph_property("double", val = np.double(np.sum(G.edge_properties["mst"].get_array()))/np.double(G.num_edges())/aabb.getVolume()) f.write("%.15f\t" % G.graph_properties["local_mst_ratio_vascular_volume"]) [lE, lV] = gt.graph_tool.centrality.eigenvector(G) G.graph_properties["Largest Eigenvector"] = G.new_graph_property("double", val = lE) f.write("%.15f\t" % lE) ############################# hybrid metrics ########################## self.add_rst_metric(G, is_dual=is_dual) f.write("%.15f\t" % G.graph_properties["sensitivity"]) f.write("%s" % label) f.write("\n") f.close() if(writeKey): location = path + "/key.txt" f = open(location, "w+") f.write("%s\t" % "global_vertex_bc") f.write("%s\t" % "global_vascular_vertex_bc") f.write("%s\t" % "global_vascular_edge_bc") f.write("%s\t" % "global_mst_ratio") f.write("%s\t" % "global_mst_ratio_vascular_volume") f.write("%s\t" % "average_length") f.write("%s\t" % "average_volume") f.write("%s\t" % "average_tortuosity") f.write("%s\t" % "vascular_length") f.write("%s\t" % "vascular_volume") f.write("%s\t" % "vascular_tortuosity") ############################Graph Statistics########################### f.write("%s\t" % "local_average_degree") f.write("%s\t" % "local_average_degree_volume") f.write("%s\t" % "local_average_degree_tortuosity") f.write("%s\t" % "local_vertex_bc") f.write("%s\t" % "local_vascular_degree") f.write("%s\t" % "local_vascular_degree_volume") f.write("%s\t" % "local_vascular_degree_tortuosity") f.write("%s\t" % "local_vascular_vertex_bc") f.write("%s\t" % "local_vascular_edge_bc") f.write("%s\t" % "local_mst_ratio") f.write("%s\t" % "local_mst_ratio_vascular_volume") f.write("%s\t" % "Largest Eigenvector") f.write("%s\t" % "label") f.write("%s\t" % "sensitivity") f.close() ''' 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.N)): G.add_node(i, p=self.N[i].p) for i in range(len(self.F)): G.add_edge(self.F[i].v0, self.F[i].v1, weight = self.F[i].length()) G[self.F[i].v0][self.F[i].v1]['pts'] = self.F[i].points G[self.F[i].v0][self.F[i].v1]['rads'] = self.F[i].radii return G ''' filters all the vertices that are close to the border with deg=1 As well as affiliated edges. It is recommended thatthe graph have all the graph metrics refreshed after filtering. ''' def filterBorder(self, G, is_dual=False): bb = AABB(G, is_dual) TFv = G.new_vertex_property("bool", vals = np.full((G.num_vertices(),1), True, dtype=bool)) TFe = G.new_edge_property("bool", vals = np.full((G.num_edges(),1), True, dtype=bool)) if(is_dual==False): for v in G.vertices(): pt = np.array([G.vertex_properties["p"][v][0], G.vertex_properties["p"][v][1], G.vertex_properties["p"][v][2]]) if G.vertex_properties["degree"][v] == 1 and bb.distance(pt) < 5.0: TFv[v] = False for e in v.all_edges(): TFe[G.edge(e.source(), e.target())] = False else: #print("I have entered this method") for v in G.vertices(): l = len(G.vertex_properties["x"][v]) pt_0 = np.array([G.vertex_properties["x"][v][0], G.vertex_properties["y"][v][0], G.vertex_properties["z"][v][0]]) pt_1 = np.array([G.vertex_properties["x"][v][l-1], G.vertex_properties["y"][v][l-1], G.vertex_properties["z"][v][l-1]]) if (G.vertex_properties["degree"][v] == 2): if ((bb.distance(pt_0) < 5.0) or (bb.distance(pt_1) < 5.0)): #print("length", l, ": ", G.vertex_properties["x"][v]) #print("degree", G.vertex_properties["degree"][v]) TFv[v] = False for e in v.all_edges(): TFe[G.edge(e.source(), e.target())] = False G.set_filters(TFe, TFv) G1 = gt.Graph(G, directed=False, prune=True) G.clear_filters() G1.clear_filters() return G1 ''' Returns a Graph object that retains all the properties of the original if return_largest is set as false, then all graphs with more than 20 nodes are returned in a list. Metrics might have to be recalculated after It is recommended that the graph have all the metrics necessary prior to filtering In order to avoid internal boost index errors. ''' def filterDisconnected(self, G, return_largest=True): #create masks clusters = G.new_vertex_property("int", vals=np.full((G.num_vertices(), 1), -1, dtype=int)) eclusters = G.new_edge_property("int", vals=np.full((G.num_edges(), 1), -1, dtype=int)) c = 0 #run a bfs that sets visited nodes on each vertex for i in G.vertices(): if clusters[i] == -1: gt.bfs_search(G, i, VisitorClassDisconnected(eclusters, clusters, c)) c = c + 1 #find the number of clusters unique, counts = np.unique(clusters.get_array(), return_counts=True) eunique, ecounts = np.unique(clusters.get_array(), return_counts=True) if(return_largest == True): #find the argument of the largest a = counts.argmax() b = ecounts.argmax() #turn the largest cluster into a TF graph mask for the vertices and the edges TFv = G.new_vertex_property("bool", vals = np.full((G.num_vertices(), 1), False, dtype=bool)) TFe = G.new_edge_property("bool", vals = np.full((G.num_edges(),1), False, dtype=bool)) for i in G.vertices(): if clusters[i] == unique[a]: TFv[i] = True e = G.get_edges(); for i in range(len(e)): if eclusters[G.edge(G.vertex(e[i][0]), G.vertex(e[i][1]))] == eunique[b]: TFe[G.edge(G.vertex(e[i][0]), G.vertex(e[i][1]))] = True #set the filters and return the vertices G.set_filters(TFe, TFv) G1=gt.Graph(G, directed=False, prune=True) G.clear_filters() G1.clear_filters() return G1 else: Graphs = [] #repeat the process above for each unique disconncted component. for j in range(len(unique)): if(counts[j] > 20): TFv = G.new_vertex_property("bool", vals = np.full((G.num_vertices(),1), False, dtype=bool)) TFe = G.new_edge_property("bool", vals = np.full((G.num_edges(),1), False, dtype=bool)) for i in range(G.num_vertices()): if clusters[G.vertex(i)] == unique[j]: TFv[G.vertex(i)] = True e = G.get_edges(); for i in range(len(e)): if eclusters[G.edge(G.vertex(e[i][0]), G.vertex(e[i][1]))] == eunique[j]: TFe[G.edge(G.vertex(e[i][0]), G.vertex(e[i][1]))] = True G.set_filters(TFe, TFv) G1=gt.Graph(G, directed=False, prune=True) G.clear_filters() G1.clear_filters() Graphs.append(G1) return Graphs def simulate_fractures(self, G, remove=10): num_removed = int(np.floor(G.num_edges()*remove/100)) if DEBUG: print("num of edges to begin with = ", G.num_edges()) indices = np.random.randint(0, int(G.num_edges()), size = [num_removed,1]) #aabb = AABB(G, is_dual) tf = np.full((G.num_edges(), 1), True, dtype=bool) for j in range(indices.shape[0]): tf[indices[j]] = False TFe = G.new_edge_property("bool", vals = tf) G.set_edge_filter(TFe) G1 = gt.Graph(G, directed=False, prune=True) G.clear_filters() G1.clear_filters() G1 = self.filterDisconnected(G1) G1.vertex_properties["degree"] = G1.degree_property_map("total") if DEBUG: print("num of edges left = ", G1.num_edges()) print("I should have removed, ", num_removed) print("Instead I removed, ", G.num_edges() - G1.num_edges()) return G1 ''' Adds the mst based- sensitivity metrics. N is the number of random removals remove is the percentage*100 of the vessels removed in each iteration. The number of vessels removed will be floor(num_edges*remove/100) ''' def add_rst_metric(self, G, N=100, remove=10, is_dual=False): #rst_ratio = G.new_graph_property("double") rst_r = 0 #G.nu num_removed = int(np.floor(G.num_edges()*remove/100)) indices = np.random.randint(0, int(G.num_edges()), size = [N, num_removed]) aabb = AABB(G, is_dual) for i in range(N): tf = np.full((G.num_edges(),1), True, dtype=bool) for j in range(num_removed): tf[indices[i, j]] = False # rst = gt.graph_tool.topology.min_spanning_tree(G, weights = G.) # ratio = np.double(np.sum(rst.get_array()))/np.double(G.num_edges()) # rst_dist.append(ratio) # rst_r = rst_r + ratio TFe = G.new_edge_property("bool", vals = tf) G.set_edge_filter(TFe) G1 = gt.Graph(G, directed=False, prune=True) G.clear_filters() G1.clear_filters() while(np.any(G1.degree_property_map("total").get_array() == True)): #print(np.any(G1.degree_property_map("total").get_array()), " ", i) G1 = self.filterDisconnected(G1) G1 = self.filterBorder(G1, is_dual) G1 = self.recalculate_metrics(G1, is_dual) if(not is_dual and G1.num_edges() > 2): G1.edge_properties["mst"] = gt.graph_tool.topology.min_spanning_tree(G1, weights=G1.edge_properties["volume"]) tvalues = G1.edge_properties["mst"].get_array() #print(tvalues) values = G1.edge_properties["inverted_volume"].get_array() #print(values) for k in range(G1.num_edges()): if(tvalues[k] == True): rst_r = rst_r + values[k] elif(is_dual and G1.num_edges() > 2): #This is non-sensical for now. G1.edge_properties["mst"] = gt.graph_tool.topology.min_spanning_tree(G1) tvalues = G1.vertex_properties["mst"].get_array() values = G1.edge_properties["inverted_volume"].get_array() for k in range(G1.num_edges()): if(tvalues[k] == True): rst_r = rst_r + values[k] G.graph_properties["sensitivity"] = G.new_graph_property("double", rst_r/N) return G ''' Re-calculates connectivity based metrics and returns the resulting graph ''' def recalculate_metrics(self, G, is_dual=False, normalize=False): gt.graph_tool.centrality.betweenness(G, vprop=G.vertex_properties["bc"], eprop=G.edge_properties["bc"], norm=normalize) if(is_dual == False): #gt.graph_tool.centrality.betweenness(G, G.vertex_properties["bc"], G.edge_properties["bc"], norm=True, is_dual=False) #generate minimum spanning tree G.edge_properties["mst"] = gt.graph_tool.topology.min_spanning_tree(G, weights=G.edge_properties["length"]) G.graph_properties["mst_ratio"] = np.double(np.sum(G.edge_properties["mst"].get_array()))/np.double(G.num_edges()) if DEBUG: print(np.double(np.sum(G.edge_properties["mst"].get_array()))/np.double(G.num_edges())) G.vertex_properties["degree"] = G.degree_property_map("total") G.vertex_properties["degree_volume"] = G.degree_property_map("total", weight=G.edge_properties["volume"]) G.vertex_properties["degree_tortuosity"] = G.degree_property_map("total", G.edge_properties["tortuosity"]) else: #gt.graph_tool.centrality.betweenness(G, G.vertex_properties["bc"], G.edge_properties["bc"], norm=True, is_dual=False) G.edge_properties["mst"] = gt.graph_tool.topology.min_spanning_tree(G, weights=G.edge_properties["bc"]) #Boolean value deligating belongance to the minimum spanning tree. G.graph_properties["mst_ratio"] = np.double(np.sum(G.edge_properties["mst"].get_array()))/np.double(G.num_edges()) G.vertex_properties["degree"] = G.degree_property_map("total") G.vertex_properties["degree_centrality"] = G.degree_property_map("total", weight=G.edge_properties["bc"]) return G ''' Filter the graph to remove the disconnected components of the graph if disconnected is set to true. if borders are set to True, the algorithm cleans up the nodes with degree one that are close to the edge of the volume. Returns a Graph object with updated internal property maps. if return_largerst == False, will return an array of graphs. disconnected == apply the disconnected componenet filter return_largest == if false returns the largest set of graphs (n > 10) borders == True, if true applies the border (deg == 1 near border) filter add_rst == True, if true adds the rst_metric ''' def filterFullGraph_gt(self, G, disconnected=True, return_largest=True, borders=True, erode = False, add_rst=False, dual=False): if(disconnected==True): if(return_largest==True and borders==True): G1 = self.filterDisconnected(G, return_largest) G1 = self.recalculate_metrics(G1, is_dual=dual) G1 = self.filterBorder(G1, is_dual=dual) if(erode == True): while(np.any(G1.degree_property_map("total").get_array() == True)): G1 = self.filterBorder(G1, is_dual=dual) G1 = self.recalculate_metrics(G1, is_dual=dual) else: G1 = self.recalculate_metrics(G1, is_dual=dual) if(add_rst == True): G1 = self.add_rst_metric(G1) elif(return_largest==True and borders==False): G1 = self.filterDisconnected(G, return_largest) if(erode == True): while(np.any(G1.degree_property_map("total").get_array() == True)): G1 = self.filterBorder(G1, is_dual=dual) G1 = self.recalculate_metrics(G1, is_dual=dual) else: G1 = self.recalculate_metrics(G1, is_dual=dual) if(add_rst == True): G1 = self.add_rst_metric(G1) elif(return_largest==False and borders == True): G1 = self.filterDisconnected(G, return_largest) for graph in G1: graph = self.recalculate_metrics(graph, is_dual=dual) graph = self.filterBorder(graph, is_dual=dual) if(erode == True): while(np.any(graph.degree_property_map("total").get_array() == True)): graph = self.filterBorder(graph, is_dual=dual) graph = self.recalculate_metrics(graph, is_dual=dual) else: graph = self.recalculate_metrics(graph, is_dual=dual) if(add_rst == True): graph = self.add_rst_metric(graph) else: G1 = self.filterDisconnected(G, return_largest) for graph in G1: graph = self.recalculate_metrics(graph, is_dual=dual) if(add_rst == True): graph = self.add_rst_metric(graph) else: G1 = self.filterBorder(G, is_dual=dual); if(erode == True): while(np.any(G1.degree_property_map("total").get_array() == True)): if DEBUG: print(G1.num_vertices()) G1 = self.filterBorder(G1, is_dual=dual) if DEBUG: print(G1.num_vertices()) G1 = self.recalculate_metrics(G1, is_dual=dual, ) else: G1 = self.recalculate_metrics(G1, is_dual=dual) if(add_rst==True): G1 = self.add_rst_metric(G1) return G1 ''' Accretes the graph in G using the nodes and edges in G_0 Assumes that G_0 has all the same properties as G, including "idx" idx is an immutable set of indices that do not change when the graph is modified ''' def accrete(self, G, G_0): v_filters = G_0.new_vertex_property("bool", vals = np.full((G_0.num_vertices(),1), False, dtype=bool)) e_filters = G_0.new_edge_property("bool", vals = np.full((G_0.num_edges(),1), False, dtype=bool)) for v in G.vertices(): #print(v) G_0_vertex = G_0.vertex(G.vertex_properties["idx"][v]) v_filters[G_0_vertex] = True for e in G_0_vertex.all_edges(): e_filters[e] = True for v_0 in G_0_vertex.all_neighbors(): v_filters[v_0] = True G_0.set_filters(e_filters, v_filters) graph = gt.Graph(G_0, directed=False, prune=True) G_0.clear_filters() return graph ''' Creates a graph from a list of nodes and a list of edges The graph is the "inverse" of the original graph, Vertices are edges and edges are ''' def createDualGraph_gt(self): #Generate the undirected graph G = gt.Graph() G.set_directed(False) #add all the required vertex properties vpos = G.new_vertex_property("vector") #original position of the edge (placed as midpoin of edge) pos = G.new_vertex_property("vector") #positions according to the fdl rpos = G.new_vertex_property("vector") #positions according to the radial layout x_edge = G.new_vertex_property("vector") #x-coordinated of the edge y_edge = G.new_vertex_property("vector") #y-coordinates of the edge z_edge = G.new_vertex_property("vector") #z-coordinates of the edge r_edge = G.new_vertex_property("vector") #r-values at each x,y,z l_edge = G.new_vertex_property("double") #length of the edge i_edge = G.new_vertex_property("double") #1/length iv_edge = G.new_vertex_property("double") #1/volume t_edge = G.new_vertex_property("double") #tortuosity of the edge v_edge = G.new_vertex_property("double") #volume of the edge vbetweeness_centrality = G.new_vertex_property("double") #betweeneness centrality of the vertex gaussian = G.new_vertex_property("double") #empty property map for gaussian values downstream degree = G.new_vertex_property("int") degree_centrality = G.new_vertex_property("int") #add all the required edge properties ebetweeness_centrality = G.new_edge_property("double") #betweeness centrality of the edge mst = G.new_edge_property("bool") #Boolean value deligating belongance to the minimum spanning tree. l_edge = G.new_edge_property("double") #length of the edge gaussian = G.new_edge_property("double") #add all graph properties mst_ratio = G.new_graph_property("double") #add all the edges as vertices for i in range(len(self.F)): if(self.F[i].length() > 0): G.add_vertex() x,y,z,r = self.F[i].getcoords() x_edge[G.vertex(i)] = x y_edge[G.vertex(i)] = y z_edge[G.vertex(i)] = z r_edge[G.vertex(i)] = r l_edge[G.vertex(i)] = self.F[i].length() gaussian[G.vertex(i)] = np.exp(-np.power(l_edge[G.vertex(i)], 2.) / (2 * np.power(60.0, 2.))) i_edge[G.vertex(i)] = 1/l_edge[G.vertex(i)] t_edge[G.vertex(i)] = self.F[i].turtuosity() v_edge[G.vertex(i)] = self.F[i].volume() iv_edge[G.vertex(i)] = 1/v_edge[G.vertex(i)] vpos[G.vertex(i)] = np.array([x[int(len(x)/2)], y[int(len(x)/2)], z[int(len(x)/2)]], dtype=float) #based on the neighbors add edges. If two edges share a vertex, there is a vertex between them. for i in range(len(self.N)): #in case there are 3 or more edges attached to a vertex if(len(self.N[i].o) > 2): for j in range(len(self.N[i].o)-1): G.add_edge(G.vertex(self.N[i].o[j]), G.vertex(self.N[i].o[j+1])) l_edge[G.edge(G.vertex(self.N[i].o[j])), G.vertex(self.N[i].o[j+1])] = 1 gaussian[G.edge(G.vertex(self.N[i].o[j])), G.vertex(self.N[i].o[j+1])] = \ np.exp(-np.power(l_edge[G.edge(G.vertex(self.N[i].o[j]),G.vertex(self.N[i].o[j+1]))], 2.) / (2 * np.power(60.0, 2.))) #add the final edge. G.add_edge(G.vertex(self.N[i].o[0]), G.vertex(self.N[i].o[len(self.N[i].o)-1])) l_edge[G.edge(G.vertex(self.N[i].o[j])), G.vertex(self.N[i].o[j+1])] = 1 gaussian[G.edge(G.vertex(self.N[i].o[j])), G.vertex(self.N[i].o[j+1])] = \ np.exp(-np.power(l_edge[G.edge(G.vertex(self.N[i].o[j]),G.vertex(self.N[i].o[j+1]))], 2.) / (2 * np.power(60.0, 2.))) #in case there are 2 edges attached to a vertex elif(len(self.N[i].o) == 2): G.add_edge(G.vertex(self.N[i].o[0]), G.vertex(self.N[i].o[1]-1)) l_edge[G.edge(G.vertex(self.N[i].o[j])), G.vertex(self.N[i].o[j+1])] = 1 gaussian[G.edge(G.vertex(self.N[i].o[j])), G.vertex(self.N[i].o[j+1])] = \ np.exp(-np.power(l_edge[G.edge(G.vertex(self.N[i].o[j]),G.vertex(self.N[i].o[j+1]))], 2.) / (2 * np.power(60.0, 2.))) #in case there is only a single edge coming into the vertex we do # not add anything #generate centrality map gt.graph_tool.centrality.betweenness(G, vbetweeness_centrality, ebetweeness_centrality, norm=True) #generate minimum spanning tree mst = gt.graph_tool.topology.min_spanning_tree(G, weights=ebetweeness_centrality) mst_ratio[G] = np.double(np.sum(mst.get_array()))/np.double(len(self.N)) degree = G.degree_property_map("total") degree_centrality = G.degree_property_map("total", weight=ebetweeness_centrality) #degree_tortuosity = G.degree_property_map("total", weight=t_edge) pos = gt.sfdp_layout(G) rpos = gt.radial_tree_layout(G, root=G.vertex(np.argmax(vbetweeness_centrality.get_array()))) #set property maps for the vertices G.vertex_properties["p"] = vpos G.vertex_properties["pos"] = pos G.vertex_properties["rpos"] = rpos G.vertex_properties["bc"] = vbetweeness_centrality G.vertex_properties["degree_centrality"] = degree_centrality G.vertex_properties["degree"] = degree G.vertex_properties["x"] = x_edge G.vertex_properties["y"] = y_edge G.vertex_properties["z"] = z_edge G.vertex_properties["r"] = r_edge G.vertex_properties["inverted"] = i_edge G.vertex_properties["inverted_volume"] = iv_edge G.vertex_properties["length"] = l_edge G.vertex_properties["tortuosity"] = t_edge G.vertex_properties["volume"] = v_edge G.vertex_properties["bc"] = vbetweeness_centrality #G.vertex_properties["pos"] = gt.fruchterman_reingold_layout(G, weight=G.edge_properties["length"], circular = True, n_iter= 10000, pos=G.vertex_properties["pos"], r = 2.0) G.vertex_properties["pos"] = gt.sfdp_layout(G, C = 0.5, p = 3.0) #set property maps for the edges G.edge_properties["bc"] = ebetweeness_centrality G.edge_properties["mst"] = mst #set graph properies G.graph_properties["mst_ratio"] = mst_ratio # title = "raw_dual_cortical_7_6_before_delete.pdf" # title2 = "raw_dual_radial_cordical_7_6_before_delete.pdf" # # gt.graph_draw(G, pos=gt.radial_tree_layout(G, root=G.vertex(np.argmax(vbetweeness_centrality.get_array()))), edge_pen_width = 2.0, vertex_size=degree, edge_color=mst, vertex_fill_color=vbetweeness_centrality, output=title2, bg_color=[0.0, 0.0,0.0,1.0], vertex_text=G.vertex_index, output_size=(1000,1000),vertex_font_size = 6) # gt.graph_draw(G, pos=pos, edge_pen_width = 4.0, vertex_size=degree, edge_color=mst, vertex_fill_color=vbetweeness_centrality, output=title, bg_color=[0.0, 0.0,0.0,1.0], vertex_text=G.vertex_index, output_size=(1000,1000),vertex_font_size = 6) # title = "./Original_2D_Graph.png" # title2 = "raw_dual_radial_cordical_7_6_after_delete.pdf" # G1 = self.filterFullGraph_gt(G, borders=False, dual=True) # #print(clusters.get_array()[:], clusters.get_array()[:]) # #pos = gt.sfdp_layout(G) # # #gt.graph_draw(G, pos=,edge_color=mst, vertex_fill_color=vbetweeness_centrality, edge_pen_width=ebetweeness_centrality, output="raw.pdf") # # gt.graph_draw(G1, pos=gt.radial_tree_layout(G1, root=G1.vertex(np.argmax(G1.vertex_properties["bc"].get_array()))), edge_pen_width = 4.0, vertex_size=G1.vertex_properties["degree"], edge_color=G1.edge_properties["mst"], vertex_fill_color=G1.vertex_properties["bc"], output=title2, bg_color=[0.0, 0.0,0.0,1.0], vertex_text=G1.vertex_index, output_size=(1000,1000),vertex_font_size = 6) # gt.graph_draw(G1, pos=G1.vertex_properties["p"], edge_pen_width = 4.0, vertex_size=G1.vertex_properties["degree"], edge_color=G1.edge_properties["mst"], vertex_fill_color=G1.vertex_properties["bc"], output=title, bg_color=[0.0, 0.0,0.0,1.0], vertex_text=G1.vertex_index, output_size=(1000,1000),vertex_font_size = 6) gt.graph_draw(G1, pos=G1.vertex_properties["p"], edge_pen_width = 4.0, vertex_size=G1.vertex_properties["degree"], output=title, bg_color=[1.0, 1.0,1.0,1.0], vertex_text=G1.vertex_index, output_size=(1000,1000),vertex_font_size = 6) title = "./Original_2D_Layout.png" gt.graph_draw(G1, pos=G1.vertex_properties["p"], edge_pen_width = 4.0, vertex_size=G1.vertex_properties["degree"], edge_color=G1.edge_properties["mst"], vertex_fill_color=G1.vertex_properties["bc"], output=title, bg_color=[1.0, 1.0,1.0,1.0], vertex_text=G1.vertex_index, output_size=(1000,1000), vertex_font_size = 6) #gt.graph_draw(G1, pos=G1.vertex_properties["pos"], edge_pen_width = 4.0, vertex_size=degree, edge_color=G1.edge_properties["mst"], vertex_fill_color=G1.vertex_properties["bc"], output=title) # # # G2 = self.filterFullGraph_gt(G1, add_rst=False, dual=True) # title = "raw_dual_cortical_7_6_after_borders.pdf" # title2 = "raw_dual_radial_cortical_7_6_after_borders.pdf" # gt.graph_draw(G2, pos=G2.vertex_properties["pos"], edge_pen_width = 4.0, vertex_size=G2.vertex_properties["degree"], edge_color=G2.edge_properties["mst"], vertex_fill_color=G2.vertex_properties["bc"], output=title, bg_color=[0.0, 0.0,0.0,1.0], vertex_text=G2.vertex_index, output_size=(1000,1000),vertex_font_size = 6) # gt.graph_draw(G2, pos=gt.radial_tree_layout(G2, root=G2.vertex(np.argmax(G2.vertex_properties["bc"].get_array()))), edge_pen_width = 4.0, vertex_size=G2.vertex_properties["degree"], edge_color=G2.edge_properties["mst"], vertex_fill_color=G2.vertex_properties["bc"], output=title2, bg_color=[0.0, 0.0,0.0,1.0], vertex_text=G2.vertex_index, output_size=(1000,1000),vertex_font_size = 6) return G ''' Create a graph from a list of nodes and a list of edges. Populates the Graph with a bunch of property maps in respect to edges and nodes each vertex is assigned 3 positional coordiantes: The original x,y,z location of the graph. x,y locations according to force directed layout, and x,y positions according to a radial layout Returns a graph_tool object ''' def createFullGraph_gt(self): #Generate the undirected graph G = gt.Graph() G.set_directed(False) #add all the required vertex properties vpos = G.new_vertex_property("vector") #original vertex positions pos = G.new_vertex_property("vector") #positions according to the fdl rpos = G.new_vertex_property("vector") #positions according to the radial layout vbetweeness_centrality = G.new_vertex_property("double") #betweeneness centrality of the vertex degree = G.new_vertex_property("int") #degree degree_volume = G.new_vertex_property("double") #degree scaled by the volume of all in fibers degree_tortuosity = G.new_vertex_property("double") #degree scaled by the tortuosity of all in fibers. #add all the required edge properties #G.properties[("x,y,z,r")] = G.new_edge_property("vector") #x-coordinated of the edge y_edge = G.new_edge_property("vector") #y-coordinates of the edge z_edge = G.new_edge_property("vector") #z-coordinates of the edge r_edge = G.new_edge_property("vector") #r-values at each x,y,z av_edge = G.new_edge_property("double") #average edge radius l_edge = G.new_edge_property("double") #length of the edge i_edge = G.new_edge_property("double") #1/length iv_edge = G.new_edge_property("double") #1/volume t_edge = G.new_edge_property("double") #tortuosity of the edge v_edge = G.new_edge_property("double") #volume of the edge ebetweeness_centrality = G.new_edge_property("double") #betweeness centrality of the edge ebc_length = G.new_edge_property("double") #betweeneness centrality scaled by the length mst = G.new_edge_property("bool") #Boolean value deligating belongance to the minimum spanning tree. #This map gets updated. gaussian = G.new_edge_property("double") #empty property map for gaussian values downstream #Graph properties (once per region) mst_ratio = G.new_graph_property("double") #Graph based metric of mst edges/total edges. #add verticies and set their properties for i in range(len(self.N)): G.add_vertex() vpos[G.vertex(i)] = self.N[i].p #add edges and set their properties. for i in range(len(self.F)): v0 = self.F[i].v0 v1 = self.F[i].v1 if self.F[i].length() > 0.0: x,y,z,r = self.F[i].getcoords() G.add_edge(G.vertex(v0), G.vertex(v1)) l_edge[G.edge(v0,v1)] = self.F[i].length() gaussian[G.edge(v0,v1)] = np.exp(-np.power(l_edge[G.edge(v0,v1)], 2.) / (2 * np.power(60.0, 2.))) i_edge[G.edge(v0,v1)] = 1/l_edge[G.edge(v0,v1)] t_edge[G.edge(v0,v1)] = self.F[i].turtuosity() v_edge[G.edge(v0,v1)] = self.F[i].volume() iv_edge[G.edge(v0,v1)] = 1/v_edge[G.edge(v0,v1)] x_edge[G.edge(v0,v1)] = x y_edge[G.edge(v0,v1)] = y z_edge[G.edge(v0,v1)] = z r_edge[G.edge(v0,v1)] = r av_edge[G.edge(v0,v1)] = self.F[i].av_radius() else: print("NON-CRITICAL ERROR: edge with length 0 detected--SKIPPED") #generate centrality map gt.graph_tool.centrality.betweenness(G, vprop=vbetweeness_centrality, eprop=ebetweeness_centrality, norm=True) #generate minimum spanning tree mst = gt.graph_tool.topology.min_spanning_tree(G, weights=l_edge) mst_ratio[G] = np.double(np.sum(mst.get_array()))/np.double(len(self.N)) degree = G.degree_property_map("total") degree_volume = G.degree_property_map("total", weight=v_edge) dg = degree_volume.get_array() dg = np.exp(-np.power(dg, 2.) / (2 * np.power(0.000005, 2.))) degree_volume = G.new_vertex_property("double", vals=dg) degree_tortuosity = G.degree_property_map("total", weight=t_edge) #print(clusters.get_array()[:], clusters.get_array()[:]) pos = gt.sfdp_layout(G, C = 1.0, K = 10) rpos = gt.radial_tree_layout(G, root=G.vertex(np.argmax(vbetweeness_centrality.get_array()))) for e in G.edges(): ebc_length[G.edge(e.source(), e.target())] = ebetweeness_centrality[G.edge(e.source(), e.target())]*l_edge[G.edge(e.source(), e.target())] #set property maps for the vertices G.vertex_properties["p"] = vpos G.vertex_properties["pos"] = pos G.vertex_properties["rpos"] = rpos G.vertex_properties["bc"] = vbetweeness_centrality G.vertex_properties["degree"] = degree G.vertex_properties["degree_volume"] = degree_volume G.vertex_properties["degree_tortuosity"] = degree_tortuosity G.vertex_properties["selection"] = G.new_vertex_property("bool", val=False) #set property maps for the edges G.edge_properties["x"] = x_edge G.edge_properties["y"] = y_edge G.edge_properties["z"] = z_edge G.edge_properties["r"] = r_edge G.edge_properties["inverted"] = i_edge G.edge_properties["length"] = l_edge G.edge_properties["tortuosity"] = t_edge G.edge_properties["av_radius"] = av_edge G.edge_properties["volume"] = v_edge G.edge_properties["bc"] = ebetweeness_centrality G.edge_properties["bc*length"] = ebc_length G.edge_properties["mst"] = mst G.edge_properties["inverted_volume"] = iv_edge G.edge_properties["gaussian"] = gaussian G.edge_properties["selection"] = G.new_edge_property("double", val=0.0) #set graph properies G.graph_properties["mst_ratio"] = mst_ratio N = 0 for e in G.edges(): N += len(G.edge_properties["x"][e]) point_cloud_x = np.zeros((N, 1), dtype=np.float64) point_cloud_y = np.zeros((N, 1), dtype=np.float64) point_cloud_z = np.zeros((N, 1), dtype=np.float64) n = 0 for e in G.edges(): for p in range(len(G.edge_properties["x"][e])): point_cloud_x[n][0] = G.edge_properties["x"][e][p] point_cloud_y[n][0] = G.edge_properties["y"][e][p] point_cloud_z[n][0] = G.edge_properties["z"][e][p] n += 1 G.graph_properties["point_cloud_x"] = G.new_graph_property("vector") G.graph_properties["point_cloud_y"] = G.new_graph_property("vector") G.graph_properties["point_cloud_z"] = G.new_graph_property("vector") G.graph_properties["point_cloud_x"] = point_cloud_x G.graph_properties["point_cloud_y"] = point_cloud_y G.graph_properties["point_cloud_z"] = point_cloud_z #save the original graph indices for the edges and the vertices. G.vertex_properties["idx"] = G.vertex_index G.edge_properties["idx"] = G.edge_index title = "~/Pictures/Original_2D_Graph_2.png" # title2 = "raw_dual_radial_cordical_7_6_after_delete.pdf" # # G1 = self.filterFullGraph_gt(G, borders=False, dual=False, add_rst = False) ## ## #print(clusters.get_array()[:], clusters.get_array()[:]) ## #pos = gt.sfdp_layout(G) ## ## #gt.graph_draw(G, pos=,edge_color=mst, vertex_fill_color=vbetweeness_centrality, edge_pen_width=ebetweeness_centrality, output="raw.pdf") ## ## gt.graph_draw(G1, pos=gt.radial_tree_layout(G1, root=G1.vertex(np.argmax(G1.vertex_properties["bc"].get_array()))), edge_pen_width = 4.0, vertex_size=G1.vertex_properties["degree"], edge_color=G1.edge_properties["mst"], vertex_fill_color=G1.vertex_properties["bc"], output=title2, bg_color=[0.0, 0.0,0.0,1.0], vertex_text=G1.vertex_index, output_size=(1000,1000),vertex_font_size = 6) ## gt.graph_draw(G1, pos=G1.vertex_properties["p"], edge_pen_width = 4.0, vertex_size=G1.vertex_properties["degree"], edge_color=G1.edge_properties["mst"], vertex_fill_color=G1.vertex_properties["bc"], output=title, bg_color=[0.0, 0.0,0.0,1.0], vertex_text=G1.vertex_index, output_size=(1000,1000),vertex_font_size = 6) # gt.graphviz_draw(G1, pos = G1.vertex_properties["p"], ratio="fill", size = (100,100), layout = None, pin = True, overlap=True, vsize=(G1.vertex_properties["degree_tortuosity"], 2.), vprops={"index":G1.vertex_index}, penwidth= 8.0, vcolor = G1.vertex_properties["bc"], vcmap = cm.get_cmap("inferno"), output=title) # #gt.graph_draw(G1, pos=G1.vertex_properties["p"], edge_pen_width = 8.0, output=title, bg_color=[1.0, 1.0,1.0,1.0], vertex_size=60, vertex_fill_color=G1.vertex_properties["bc"], vertex_text=G1.vertex_index, output_size=(3200,3200),vertex_font_size = 32) # title = "~/Pictures/Original_2D_Layout.png" # gt.graph_draw(G1, pos=G1.vertex_properties["pos"], edge_pen_width = 8.0, edge_color=G1.edge_properties["mst"], vertex_size=40, vertex_fill_color=G1.vertex_properties["bc"], output=title, bg_color=[1.0, 1.0,1.0,1.0], vertex_text=G1.vertex_index, output_size=(3200,3200), vertex_font_size = 32) # # # G1 = self.filterFullGraph_gt(G, borders=True, dual=False, erode=True, add_rst = False) # title = "~/Pictures/Original_2D_Graph_Eroded.png" # gt.graph_draw(G1, pos=G1.vertex_properties["p"], edge_pen_width = 8.0, output=title, bg_color=[1.0, 1.0,1.0,1.0], vertex_size=60, vertex_fill_color=G1.vertex_properties["bc"], vertex_text=G1.vertex_index, output_size=(3200,3200),vertex_font_size = 32) # title = "~/Pictures/Original_2D_Layout_Eroded.png" # gt.graph_draw(G1, pos=G1.vertex_properties["pos"], edge_pen_width = 8.0, edge_color=G1.edge_properties["mst"], vertex_size=40, vertex_fill_color=G1.vertex_properties["bc"], output=title, bg_color=[1.0, 1.0,1.0,1.0], vertex_text=G1.vertex_index, output_size=(3200,3200), vertex_font_size = 32) # title = "raw_full_cortical_7_6_before_delete.pdf" # title2 = "raw_full_radial_cordical_7_6_before_delete.pdf" # #G.graph_properties["rst_ratio"] = rst_ratio # #print("") # #print(G.graph_properties["rst_ratio"]) # #print(G.edge_properties["bc"].get_array()) # #print() # # #G.vertex_properties[("p","betweeness_centrality")] = [vpos, vbetweeness_centrality] # #G.edge_properties[("x", "y", "z", "r", "length", "tortuosity", "volume", "betweeness_centrality")] = [] # #G.set_edge_filter(mst) # #G.set_edge_filter(mst) # # #Experimental code for drawing graphs. # #gt.graph_draw(G,pos=vpos,edge_color=mst, vertex_fill_color=vbetweeness_centrality, edge_pen_width=ebetweeness_centrality, output="raw.pdf") # #gt.graph_draw(G,pos=pos, vertex_size=degree, edge_color=mst, vertex_fill_color=vbetweeness_centrality, output=title, bg_color=[1.0,1.0,1.0]) # #gt.graph_draw(G,pos=gt.radial_tree_layout(G, root=G.vertex(np.argmax(vbetweeness_centrality.get_array())), node_weight=vbetweeness_centrality), vertex_size=degree, edge_color=mst, vertex_fill_color=vbetweeness_centrality, output=title2, bg_color=[1.0, 1.0,1.0,1.0], vertex_text=G.vertex_index, output_size=(1000,1000), vertex_font_size = 6) # #print(np.argmax(vbetweeness_centrality.get_array())) # # #np.savetxt("vbetweeness_centrality", vbetweeness_centrality.get_array()) # gt.graph_draw(G,pos=gt.radial_tree_layout(G, root=G.vertex(np.argmax(vbetweeness_centrality.get_array()))), edge_pen_width = 2.0, vertex_size=degree, edge_color=mst, vertex_fill_color=vbetweeness_centrality, output=title2, bg_color=[0.0, 0.0,0.0,1.0], vertex_text=G.vertex_index, output_size=(1000,1000),vertex_font_size = 6) # gt.graph_draw(G, pos=pos, edge_pen_width = 4.0, vertex_size=degree, edge_color=mst, vertex_fill_color=vbetweeness_centrality, output=title, bg_color=[0.0, 0.0,0.0,1.0], vertex_text=G.vertex_index, output_size=(1000,1000),vertex_font_size = 6) # # # #gt.graph_draw(G,pos=vpos, efilt=mst, vertex_fill_color=vbetweeness_centrality, edge_pen_width=ebetweeness_centrality, output="mst.pdf") # #gt.graph_draw(G,pos=pos, efilt=mst, vertex_fill_color=vbetweeness_centrality, output="mst.pdf") # #need to create subgraphs!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! # #rst_r = 0 # #for i in range(1000): # # rst = gt.graph_tool.topology.random_spanning_tree(G) # # rst_r = rst_r + np.double(np.sum(rst.get_array()))/np.double(len(self.N)) # # #rst_ratio[G] = rst_r/1000.0 # # # # #G1 = gt.Graph(G) # # #G1.purge_edges() # #G1.purge_vertices(in_place = True) # # # title = "raw_full_cortical_7_6_after_delete.pdf" # title2 = "raw_full_radial_cordical_7_6_after_delete.pdf" # # G1 = self.filterFullGraph_gt(G,borders=False) # # #print(clusters.get_array()[:], clusters.get_array()[:]) # #pos = gt.sfdp_layout(G) # # #gt.graph_draw(G,pos=vpos,edge_color=mst, vertex_fill_color=vbetweeness_centrality, edge_pen_width=ebetweeness_centrality, output="raw.pdf") # # gt.graph_draw(pos = gt.sfdp_layout(G1), output=title) # gt.graph_draw(G1, pos=G1.vertex_properties["pos"], edge_pen_width = 4.0, vertex_size=G1.vertex_properties["degree"], edge_color=G1.edge_properties["mst"], vertex_fill_color=G1.vertex_properties["bc"], output=title, bg_color=[0.0, 0.0,0.0,1.0], vertex_text=G1.vertex_index, output_size=(1000,1000),vertex_font_size = 6) # #gt.graph_draw(G1, pos=G1.vertex_properties["pos"], edge_pen_width = 4.0, vertex_size=degree, edge_color=G1.edge_properties["mst"], vertex_fill_color=G1.vertex_properties["bc"], output=title) # # # G2 = self.filterFullGraph_gt(G1) # title = "raw_full_cortical_7_6_after_borders.pdf" # title2 = "raw_full_radial_cortical_7_6_after_borders.pdf" # gt.graph_draw(G2, pos=G2.vertex_properties["pos"], edge_pen_width = 4.0, vertex_size=G2.vertex_properties["degree"], edge_color=G2.edge_properties["mst"], vertex_fill_color=G2.vertex_properties["bc"], output=title, bg_color=[0.0, 0.0,0.0,1.0], vertex_text=G2.vertex_index, output_size=(1000,1000),vertex_font_size = 6) # gt.graph_draw(G2, pos=gt.radial_tree_layout(G2, root=G2.vertex(np.argmax(G2.vertex_properties["bc"].get_array()))), edge_pen_width = 4.0, vertex_size=G2.vertex_properties["degree"], edge_color=G2.edge_properties["mst"], vertex_fill_color=G2.vertex_properties["bc"], output=title2, bg_color=[0.0, 0.0,0.0,1.0], vertex_text=G2.vertex_index, output_size=(1000,1000),vertex_font_size = 6) # #gt.graph_draw(G,pos=pos, vertex_size=degree, edge_color=mst, vertex_fill_color=vbetweeness_centrality, output=title, bg_color=[1.0,1.0,1.0]) # # #gt.graph_draw(G,pos=gt.radial_tree_layout(G, root=G.vertex(np.argmax(vbetweeness_centrality.get_array())), node_weight=vbetweeness_centrality), vertex_size=degree, edge_color=mst, vertex_fill_color=vbetweeness_centrality, output=title2, bg_color=[1.0, 1.0,1.0,1.0], vertex_text=G.vertex_index, output_size=(1000,1000), vertex_font_size = 6) # #print(np.argmax(vbetweeness_centrality.get_array())) # #np.savetxt("vbetweeness_centrality", vbetweeness_centrality.get_array()) # #gt.graph_draw(G1, pos=gt.radial_tree_layout(G1, root=G1.vertex(np.argmax(G1.vertex_properties["bc"].get_array()))), edge_pen_width = 2.0, vertex_size=G1.vertex_properties["degree"], edge_color=G1.edge_properties["mst"], vertex_fill_color=G1.vertex_properties["bc"], output=title2, bg_color=[0.0, 0.0,0.0,1.0], vertex_text=G1.vertex_index, output_size=(1000,1000),vertex_font_size = 6) # #gt.graph_draw(G1, pos=gt.radial_tree_layout(G1, root=G1.vertex(np.argmax(G1.vertex_properties["bc"].get_array()))), edge_pen_width = 2.0, vertex_size=degree, edge_color=G1.edge_properties["mst"]) # #gt.draw_hierarchy(gt.minimize_nested_blockmodel_dl(G1, deg_corr=False), output = "test.pdf") # #gt.graph_draw(G1,pos=gt.radial_tree_layout(G1, root=G1.vertex(np.argmax(G1.edge_properties["mst"].get_array()))), edge_pen_width = 2.0, vertex_size=degree, edge_color=mst, vertex_fill_color=vbetweeness_centrality, output=title2, bg_color=[0.0, 0.0,0.0,1.0], vertex_text=G.vertex_index, output_size=(1000,1000),vertex_font_size = 6) # print(gt.graph_tool.openmp_enabled(), gt.graph_tool.openmp_get_num_threads()) #print(ratio, G2.graph_properties["mst_ratio"]) #n, bins, patches = plt.hist(np.asarray(rst_dist), 50, normed=1, facecolor='green', alpha=0.75) # add a 'best fit' line #plt.grid(True) #plt.show() #self.partition(G, 5, 7, False) self.export_points(G) self.export_bb(G) return G def gen_new_fd_layout(G): pos = gt.sfdp_layout(G, C = 1.0, K = 10) G.vertex_properties["pos"] = pos return G def map_edges_to_range(G, rng, propertymap): def func(maximum, minimum, new_maximum, new_minimum, value): return ((value-minimum)/(maximum-minimum)*(new_maximum-new_minimum)+new_minimum) mx = max(G.edge_properties[propertymap].get_array()) mn = min(G.edge_properties[propertymap].get_array()) G.edge_properties["map"] = G.new_edge_property("float") gt.map_property_values(G.edge_properties[propertymap], G.edge_properties["map"], lambda x: func(mx, mn, rng[0], rng[1], x)) return G.edge_properties["map"] def map_vertices_to_range(G, rng, propertymap): def func(maximum, minimum, new_maximum, new_minimum, value): return ((value-minimum)/(maximum-minimum)*(new_maximum-new_minimum)+new_minimum) mx = max(G.vertex_properties[propertymap].get_array()) mn = min(G.vertex_properties[propertymap].get_array()) G.vertex_properties["map"] = G.new_vertex_property("float") gt.map_property_values(G.vertex_properties[propertymap], G.vertex_properties["map"], lambda x: func(mx, mn, rng[0], rng[1], x)) return G.vertex_properties["map"] ''' G and propertymap in G are passed in. Maps select property from range to color, works for int's floats and vec3 and vec4 values Returns a vector property map that maps the property value to color. If passed a vertex propertymap, returns a vertex property map If passed an edge propertymap, returns an edge property map ''' def map_property_to_color(G, propertymap, colormap = 'tab20'): key_type = '' if propertymap.key_type()=='e': key_type = 'e' elif propertymap.key_type()=='v': key_type = 'v' else: #TO DO: Write graph property type return. print("Graph_property passed") if G.properties[(key_type, 'RGBA')] == None: color = [0.0, 0.0, 0.0, 1.0] G.properties[(key_type, 'RGBA')] = G.new_property(key_type, "vector", val=color) #This is when the property map is integers #int32_t when integer value_type = propertymap.value_type() if DEBUG: print(value_type) if(value_type == "int32_t"): array = propertymap.get_array() colors = cm.get_cmap(colormap, len(np.unique(array))) if key_type == 'v': for v in G.vertices(): G.vertex_properties["RGBA"][v] = colors(propertymap[v]) return G.vertex_properties["RGBA"] elif(value_type == 'double' or value_type == 'float'): array = propertymap.get_array() norm = cm.colors.Normalize(vmin=min(array), vmax=max(array)) #colors = cm.get_cmap(colormap, len(np.unique(array))) colors = cm.ScalarMappable(norm=norm, cmap=colormap) if key_type =='v': for v in G.vertices(): if DEBUG: print(colors.to_rgba(propertymap[v])) G.vertex_properties["RGBA"][v] = colors.to_rgba(propertymap[v]) return G.vertex_properties["RGBA"] elif(value_type == 'vector' or value_type == 'vector'): print("detected a vector of doubles or floats") ''' Attempts to partition the graph G into b_i = [0, n^3] sections, where each section is the connected componenets of m steps away from the source node (where the source node is the closest node to point n_i) ''' def partition(self, G, n, m, is_dual, path, prefix, saveKey, lbl = "none"): bb = AABB(G, is_dual) #print("FLJKKHDFLKJFDLKJFDLKJ ", m) x, y, z = bb.project_grid(n) #cluster_belongance = G.new_vertex_property("vector") #ecluster_belongance = G.new_edge_property("vector") #X, Y, Z = np.meshgrid(x,y,z) array_cluster_bel = np.zeros([n**3, G.num_vertices()]) array_ecluster_bel = np.zeros([n**3, G.num_edges()]) c = 0 for i in range(len(x)): for j in range(len(y)): for k in range(len(z)): #for i in range(1): # for j in range(1): # for k in range(1): p = G.vertex_properties["p"].get_2d_array(range(G.num_vertices())) if(p.shape[1] < 3): break #print("stuff", p) idx = -1 dist = 100000000000 for M in range(p.shape[1]): d = math.sqrt((x[i] - p[0][M])**2 + (y[j] - p[1][M])**2 + (z[k] - p[2][M])**2) if d < dist: idx = M dist = d #if DEBUG: print(idx, " vertex[",p[0][idx], p[1][idx], p[2][idx], "], point [", x[i], y[j], z[k], "]") clusters = G.new_vertex_property("int", vals=np.full((G.num_vertices(), 1), -1, dtype=int)) #eclusters = G.new_edge_property("int", vals=np.full((G.num_edges(), 1), -1, dtype=int)) dist = G.new_vertex_property("int", val = 100000) dist[G.vertex(idx)] = 0 #run a bfs on the index node gt.bfs_search(G, idx, VisitorClassPartition(clusters, dist, c, m)) #print(dist.get_array()); temp_vertices = np.ma.masked_where(dist.get_array() <= m, dist.get_array()).mask #print(temp_vertices) temp_edges = np.zeros([G.num_edges()]) for vertex in np.nditer(np.where(temp_vertices == True)): for edge in G.get_out_edges(vertex): #print(edge) temp_edges[edge[2]] = 1 array_cluster_bel[c,:] = temp_vertices[:] array_ecluster_bel[c,:] = temp_edges[:] c = c + 1 #print(G.num_vertices()) #np.savetxt("clusters.txt", array_cluster_bel) G.vertex_properties["partition"] = G.new_vertex_property("vector", array_cluster_bel) G.edge_properties["partition"] = G.new_edge_property("vector", array_ecluster_bel) j = 0 gt.graph_draw(G, pos = gt.sfdp_layout(G), output="Graph_full.pdf") for i in range(c): TFv = G.new_vertex_property("bool", vals=array_cluster_bel[i,:]) TFe = G.new_edge_property("bool", vals=array_ecluster_bel[i,:]) G.set_filters(TFe, TFv) G1 = gt.Graph(G, directed=False, prune=True) G.clear_filters() G1.clear_filters() iteration = 0 title = "graph_" + str(iteration) + ".pdf" gt.graph_draw(G1, pos = gt.sfdp_layout(G1), output=title) while(np.any(G1.degree_property_map("total").get_array() == True)): #print(np.any(G1.degree_property_map("total").get_array()), " ", i) iteration = iteration + 1 G1 = self.filterBorder(G1, is_dual) title = "graph_" + str(iteration) + ".pdf" gt.graph_draw(G1, pos = gt.sfdp_layout(G1), output=title) #G1 = self.filterDisconnected(G1, is_dual) G1 = self.recalculate_metrics(G1, is_dual) if(G1.num_edges() > 2): ppath = path + "/partial"+str(j)+".nwt" self.saveGraph_gt(G1, ppath) if(i == 0): self.saveSample_gt(G1, path, prefix, is_dual, writeKey=saveKey, label=lbl) else: self.saveSample_gt(G1, path, prefix, is_dual, label=lbl) j = j + 1 #print(array_cluster_bel) def find_loops(self, G, path, is_dual, largest): iteration = 0 if(largest): G = self.filterDisconnected(G) gt.graph_draw(G, pos = gt.sfdp_layout(G), output="./graph_full.pdf") while(np.any(G.degree_property_map("total").get_array() == True)): #print(np.any(G1.degree_property_map("total").get_array()), " ", i) iteration = iteration + 1 G = self.filterBorder(G, is_dual) title = path +"graph_" + str(iteration) + ".pdf" gt.graph_draw(G, pos = gt.sfdp_layout(G), output=title) #G1 = self.filterDisconnected(G1, is_dual) G = self.recalculate_metrics(G, is_dual) ppath = path + "graph_full_mst.pdf" gt.graph_draw(G, pos = gt.sfdp_layout(G, C=0.4), output=ppath, edge_color=G.edge_properties["mst"], vertex_text=G.vertex_index,output_size=(1000,1000), vertex_font_size = 6, vertex_size = 3) paths = G.new_edge_property("int", vals=np.full((G.num_edges(), 1), 1, dtype=int)) #values = np.full((G.num_vertices(),1), 0, dtype=int) #for i in range(G.num_vertices()): # values[i,0] = G.vertex_index[G.vertex(i)] #G.vertex_properties["original_index"] = values #add index, then use index to search for new vertecies. G.edge_properties["path_weights"] = paths cycle_list = [] for e in G.edges(): if G.edge_properties["mst"][e] == 0: cycle_list.append(e) # else: # G.edge_properties["path_weights"][e] = 1 #for e in cycle_list: #print("Loop Detected between %i, %i", e.source(), e.target()) G.set_edge_filter(G.edge_properties["mst"]) iteration = 0 G1 = gt.Graph(G, directed=False, prune=True) G.clear_filters() for e in cycle_list: edge = G1.add_edge(e.target(), e.source()) #G.edge_properties["mst"][e] = 1 G1.edge_properties["path_weights"][edge] = 1000 #G.set_edge_filter(G.edge_properties["mst"]) #G1 = gt.Graph(G, directed=False, prune=True) #s = -1 #t = -1 #for edge in G1.edges(): # if(G1.vertex_properties["original_index"][edge.source()] == int(e.source())): # s = int(edge.source()) # elif(G1.vertex_properties["original_index"][edge.target()] == int(e.target())): # t = int(edge.target()) [Vl, El] = gt.shortest_path(G1, e.source(), e.target(), weights=G1.edge_properties["path_weights"]) clusters = G1.new_vertex_property("int", vals=np.full((G1.num_vertices(), 1), 5, dtype=int)) eclusters = G1.new_edge_property("int", vals = np.full((G1.num_edges(), 1), 3, dtype=int)) G1.vertex_properties["cycle"] = clusters G1.edge_properties["cycle"] = eclusters if DEBUG: print("Number of vertices in Path is:", len(Vl)) for v in Vl: G1.vertex_properties["cycle"][G1.vertex(v)] = 10 if DEBUG: print(str(v)) #Create the arrays to be histogrammed later regarding every loop length_total = 0 volume_total = 0 tortuosity_total = 0 bc_total = 0 bcl_total = 0 num_edges = 1 for v in range(len(Vl)-1): G1.edge_properties["cycle"][G1.edge(Vl[v], Vl[v+1])] = 10 length_total = length_total + G1.edge_properties["length"][G1.edge(Vl[v], Vl[v+1])] volume_total = volume_total + G1.edge_properties["volume"][G1.edge(Vl[v], Vl[v+1])] tortuosity_total = tortuosity_total + G1.edge_properties["tortuosity"][G1.edge(Vl[v], Vl[v+1])] bc_total = bc_total + G1.edge_properties["bc"][G1.edge(Vl[v], Vl[v+1])] bcl_total = bcl_total + G1.edge_properties["bc*length"][G1.edge(Vl[v], Vl[v+1])] num_edges = num_edges + 1 G1.edge_properties["cycle"][G1.edge(e.source(), e.target())] = 10 length_total = length_total + G.edge_properties["length"][G.edge(e.source(), e.target())] volume_total = volume_total + G.edge_properties["volume"][G.edge(e.source(), e.target())] tortuosity_total = tortuosity_total + G.edge_properties["tortuosity"][G.edge(e.source(), e.target())] bc_total = bc_total + G.edge_properties["bc"][G.edge(e.source(), e.target())] bcl_total = bcl_total + G.edge_properties["bc*length"][G.edge(e.source(), e.target())] ppath = path + "Loop_Histograms.txt" f = open(ppath, "a+") f.write("%.15f\t" % length_total) f.write("%.15f\t" % volume_total) f.write("%.15f\t" % tortuosity_total) f.write("%.15f\t" % bc_total) f.write("%.15f\t" % bcl_total) f.write("%d\t\n" % num_edges) f.close() title = path+"cycle" + str(iteration)+".png" gt.graph_draw(G1, pos=G1.vertex_properties["pos"], output=title, edge_color='red', edge_pen_width=G1.edge_properties["cycle"], vertex_text=G1.vertex_index, output_size=(1920,1280), vertex_font_size = 4, vertex_fill_color = G1.vertex_properties["cycle"], vertex_size = 3) G1.remove_edge(edge) #G.clear_filters(); #G.edge_properties["mst"][e] = 0 #G.edge_properties["path_weights"][e] = 1000 iteration = iteration + 1 ''' returns an affinity matrix based on the property edge property, given as a string. ''' def get_affinity_matrix(G, edge_property, gaussian=False, sigma = None): affinity = gt.shortest_distance(G, weights=G.edge_properties[edge_property]) affinity = affinity.get_2d_array(range(G.num_vertices())) if gaussian: if sigma == None: sig = 60.0 else: sig = sigma for i in range(G.num_vertices()): for j in range(G.num_vertices()): affinity[i,j] = np.exp(-np.power(affinity[i, j], 2.) / (2 * np.power(sig, 2.))) return affinity ''' Attempts to apporixmate the number of clusters in the graph by finding the approximate number of minimum graph cuts TO_DO: Program the gaussian and sigma parameters ''' def approximate_n_cuts(G, affinity_property): #L = gt.laplacian(G, weight = G.edge_properties[affinity_property]) G1 = gt.Graph(G, directed=False) G1 = Network.filterDisconnected(G1, G1) L = Network.get_affinity_matrix(G1, affinity_property, gaussian=True, sigma=157) L = sp.sparse.csgraph.laplacian(L, normed=True) n_components = G1.num_vertices() eigenvalues, eigenvectors = eigsh(L, k=n_components, which = "LM", sigma=1.0, maxiter = 5000) #plt.figure() #plt.scatter(np.arange(len(eigenvalues)), eigenvalues) #plt.grid() #plt.show() count = sum(eigenvalues > 1.01) return count ''' Cluster the vectices in the graph based on a property passed to the clustering algorithm TO_DO: program the guassian and sigma parameters ''' def spectral_clustering(G, affinity_property, gaussian = False, sigma = None, n_clusters = None, num_iters = 1000): if(n_clusters == None): n_c = Network.approximate_n_cuts(G, affinity_property); else: n_c = n_clusters sc = SpectralClustering(n_c, affinity='precomputed', n_init = num_iters) sc.fit(Network.get_affinity_matrix(G, affinity_property, gaussian = True, sigma=157)) return sc.labels_ def export_points(self, G): location = "./vertex_points.txt" f = open(location, "a+") for v in G.vertices(): f.write("%.15f\t" % G.vertex_properties["p"][v][0]) f.write("%.15f\t" % G.vertex_properties["p"][v][1]) f.write("%.15f\n" % G.vertex_properties["p"][v][2]) f.close() def export_bb(self, G): location = "./bb_points.txt" aabb = AABB(G) points = aabb.vertices() f = open(location, "a+") for i in points: f.write("%.15f\t" % i[0]) f.write("%.15f\t" % i[1]) f.write("%.15f\n" % i[2]) f.close() def write_vtk(G, location, camera = None, binning = True): #location = "./nwt.vtk" f = open(location, "w+") f.write("# vtk DataFile Version 1.0\n") f.write("Data values\n") f.write("ASCII\n\n") f.write("DATASET POLYDATA\n") num_pts = 0 for e in G.edges(): X = G.edge_properties["x"][e] num_pts += len(X) f.write("POINTS %d float\n" % num_pts) for e in G.edges(): X = G.edge_properties["x"][e] Y = G.edge_properties["y"][e] Z = G.edge_properties["z"][e] #pts = np.array([X,Y,Z]).T for p in range(len(X)): f.write("%.15f\t" % X[p]) f.write("%.15f\t" % Y[p]) f.write("%.15f\n" % Z[p]) f.write("LINES " + str(G.num_edges()) + " " + str(G.num_edges()+num_pts) + "\n") num_pts = 0 for e in G.edges(): X = G.edge_properties["x"][e] indices = list(range(0, len(X))) indices = [x+num_pts for x in indices] f.write(str(len(X)) + " ") f.write(" ".join(str(x) for x in indices)) f.write("\n") num_pts += len(X) f.write("POINT_DATA %d\n" % num_pts) f.write("SCALARS radius float 1\n") f.write("LOOKUP_TABLE radius_table\n") for e in G.edges(): R = G.edge_properties["r"][e] #pts = np.array([X,Y,Z]).T for p in range(len(R)): f.write("%.15f\n" % R[p]) f.write("COLOR_SCALARS color_table %d\n" % 4) bins = [0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0] if camera == None: for e in G.edges(): X = G.edge_properties["x"][e] #RGBA = G.edge_properties["RGBA"].get_2d_array(range(4)).T for p in range(len(X)): f.write("%.15f\t" % G.edge_properties["RGBA"][e][0]) f.write("%.15f\t" % G.edge_properties["RGBA"][e][1]) f.write("%.15f\t" % G.edge_properties["RGBA"][e][2]) if binning: index = np.digitize(G.edge_properties["RGBA"][e][3], bins, right=True) if (index >= len(bins) or index < 0): if DEBUG: print(G.edge_properties["RGBA"][e][3]) print(index) f.write("%.15f\n" % bins[index]) else: f.write("%.15f\n" % G.edge_properties["RGBA"][e][3]) else: ptx = G.graph_properties["point_cloud_x"].get_array() pty = G.graph_properties["point_cloud_y"].get_array() ptz = G.graph_properties["point_cloud_z"].get_array() pts = np.array([ptx,pty,ptz]).T temp = AABB(G) bbl = temp.A bbu = temp.B tp = np.zeros((1,3), dtype = np.float64) tp[0][0] = (bbu[0]-bbl[0])/2 tp[0][1] = (bbu[1]-bbl[1])/2 tp[0][2] = (bbu[2]-bbl[2])/2 pts[:] = pts[:] - camera - \ tp distance = np.zeros(pts.shape[0], dtype = np.float32) for i in range(distance.shape[0]): distance[i] = np.sqrt(np.power(pts[i][0],2) + np.power(pts[i][1],2) + np.power(pts[i][2],2)) mx = max(distance) mn = min(distance) distance = (distance - mx)/(mn-mx) idx = 0 for e in G.edges(): X = G.edge_properties["x"][e] for p in range(len(X)): f.write("%.15f\t" % G.edge_properties["RGBA"][e][0]) f.write("%.15f\t" % G.edge_properties["RGBA"][e][1]) f.write("%.15f\t" % G.edge_properties["RGBA"][e][2]) if binning: index = np.digitize(distance[idx], bins, right=True) f.write("%.15f\n" % bins[index]) else: f.write("%.15f\n" % distance[idx]) idx += 1 f.close() ''' Creates a graph from a list of nodes and a list of edges. Uses edge turtuosity as weight. Returns a NetworkX Object. ''' def createTortuosityGraph(self): G = nx.Graph() for i in range(len(self.N)): G.add_node(i, p=self.N[i].p) for i in range(len(self.F)): G.add_edge(self.F[i].v0, self.F[i].v1, weight = self.F[i].turtuosity()) G[self.F[i].v0][self.F[i].v1]['pts'] = self.F[i].points G[self.F[i].v0][self.F[i].v1]['rads'] = self.F[i].radii return G ''' Creates a graph from a list of nodes and a list of edges. Uses edge volume as weight. Returns a NetworkX Object. ''' def createVolumeGraph(self): G = nx.Graph() for i in range(len(self.N)): G.add_node(i, p=self.N[i].p) for i in range(len(self.F)): G.add_edge(self.F[i].v0, self.F[i].v1, weight = self.F[i].volume()) G[self.F[i].v0][self.F[i].v1]['pts'] = self.F[i].points G[self.F[i].v0][self.F[i].v1]['rads'] = self.F[i].radii 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): lower = self.N[0].p.copy() upper = lower.copy() for i in self.N: 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 #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 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]) 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 #returns the number of points in the network 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 #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 #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 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 start, end #returns the points inside the fiber. def fiber(self, idx): p = self.F[idx].points X = np.zeros(len(p)) Y = np.zeros(len(p)) Z = np.zeros(len(p)) for i in range(len(p)): X[i] = p[i][0] Y[i] = p[i][1] Z[i] = p[i][2] return X, Y, Z #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