network_dep.py 113 KB
``````# -*- 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
'''
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]

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):

#print(volume)
return volume

class NWT:

'''
Writes the header given and open file descripion, number of verticies and number of edges.
'''
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('<f',vertex.p[0]))
open_file.write(struct.pack('<f',vertex.p[1]))
open_file.write(struct.pack('<f',vertex.p[2]))
open_file.write(struct.pack('i', len(vertex.o)))
open_file.write(struct.pack('i', len(vertex.i)))
for j in range(len(vertex.o)):
open_file.write(struct.pack('i',vertex.o[j]))

for j in range(len(vertex.i)):
open_file.write(struct.pack('i', vertex.i[j]))

return

'''
Writes a single fiber to a file.
'''
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)))
for j in range(len(edge.points)):
open_file.write(struct.pack('<f', edge.points[j][0]))
open_file.write(struct.pack('<f', edge.points[j][1]))
open_file.write(struct.pack('<f', edge.points[j][2]))

return

'''
Writes the entire network to a file in str given the vertices array and the edges array.
'''
def exportNWT(str, vertices, edges):
with open(str, "wb") as file:
for i in range(len(vertices)):
NWT.writeVertex(file, vertices[i])

for i in range(len(edges)):
NWT.writeFiber(file, edges[i])

return

'''
Reads a single vertex from an open file and returns a node Object.
'''
points = np.tile(0., 3)
points[0] = struct.unpack('f', bytes)[0]
points[1] = struct.unpack('f', bytes)[0]
points[2] = struct.unpack('f', bytes)[0]

numO = int.from_bytes(bytes, byteorder='little')
outgoing = np.tile(0, numO)
numI = int.from_bytes(bts, byteorder='little')
incoming = np.tile(0, numI)
for j in range(numO):
outgoing[j] = int.from_bytes(bytes, byteorder='little')

for j in range(numI):
incoming[j] = int.from_bytes(bytes, byteorder='little')

node = Node(points, outgoing, incoming)
return node

'''
Reads a single fiber from an open file and returns a Fiber object .
'''
vtx0 = int.from_bytes(bytes, byteorder = 'little')
vtx1 = int.from_bytes(bytes, byteorder = 'little')
numVerts = int.from_bytes(bytes, byteorder = 'little')
pts = []

for j in range(numVerts):
point = np.tile(0., 3)
point[0] = struct.unpack('f', bytes)[0]
point[1] = struct.unpack('f', bytes)[0]
point[2] = struct.unpack('f', bytes)[0]
pts.append(point)

F = Fiber(vtx0, vtx1, pts, rads)

return F

'''
Imports a NWT file at location str.
Returns a list of Nodes objects and a list of Fiber objects.
'''

'''
Class defining an aabb around the graph
'''

class AABB():
def __init__(self, G, is_dual=False):
#designate whether we will be generating a bounding box for a dual graph
#   or a normal graph.
self.is_dual = is_dual
#minimum vertex
self.A = np.full(3, 1000000.0, dtype=float)
#maximum vertex
self.B = np.full(3, -1000000.0, dtype=float)
if(is_dual == False):
#find the minumum and the maximum of the graph.
for v in G.vertices():
if G.vertex_properties["p"][v][0] < self.A[0]:
self.A[0] = G.vertex_properties["p"][v][0]
if G.vertex_properties["p"][v][0] > 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=None, clock=False):
if clock:
start_time = time.time()
if filename!=None:
with open(filename, "rb") as file:
numVertex = int.from_bytes(bytes, byteorder='little')
numEdges = int.from_bytes(bytes, byteorder='little')

self.N = []
self.F = []
for i in range(numVertex):
self.N.append(node)

for i in range(numEdges):
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 ##########################
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)):
for i in range(len(self.F)):
G[self.F[i].v0][self.F[i].v1]['pts'] = self.F[i].points

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:
if G.vp["exclude"][v] == False:
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):
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
'''
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)
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)
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)
else:
G1 = self.filterDisconnected(G, return_largest)
for graph in G1:
graph = self.recalculate_metrics(graph, is_dual=dual)
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)

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

'''
Generates a layout based a 3 neighbors spring system.
Vertices within 1 edge, 2 edges and 3 edges act upon one another based on
Hooke's law. 1e springs are the stronges and strength decreases as edges get further away

'''
def gen_spring_onion_layout(G, v_property, k, n_steps, time_step, min_distance, Sum=True):
start = time.time()

one_neighbors = G.new_vertex_property("vector<int>")
o1_n = []
two_neighbors = G.new_vertex_property("vector<int>")
o2_n = []
three_neighbors = G.new_vertex_property("vector<int>")
o3_n = []
#generate the map of 1, 2 and 3 neighborhood vertices.
for v in G.vertices():
one_n = []
two_n = []
three_n = []
#find all 1 neighbors
for i in v.all_neighbors():
one_n.append(int(i))
if(len(one_n) == 0):
print("WTF")
#find all 2 neighbors
for i in one_n:
for j in G.vertex(i).all_neighbors():
#                    if(int(j) not in one_n):
two_n.append(int(j))

#find all 3 neighbors
for i in two_n:
for j in G.vertex(i).all_neighbors():
#                    if(int(j) not in one_n and int(j) not in two_n):
three_n.append(int(j))

#we do not keep ones that have already counted, i.e. a vertex
#appears in all 3 arrays only once.
#it's ok to add self to any array because the distance -k(x=0) = 0
one_neighbors[v] = np.asarray(one_n)
o1_n.append(np.asarray(one_n))
two_neighbors[v] = np.asarray(two_n)
o2_n.append(np.asarray(two_n))
three_neighbors[v] = np.asarray(three_n)
o3_n.append(np.asarray(three_n))

print("generating neighborhoods: ", time.time() - start)
start = time.time()
G.vertex_properties["1_neighbor"] = one_neighbors
G.vertex_properties["2_neighbor"] = two_neighbors
G.vertex_properties["3_neighbor"] = three_neighbors

cooling = 1.0
cooling_step = cooling/n_steps

#        #get points of the centers of every cluster
#        #generate voronoi region and plot it.
#        fig, ax = plt.subplots(2, 1, sharex='col', sharey='row')
#        fig.tight_layout()
#        grid = plt.GridSpec(2,1)
#        grid.update(wspace=0.025, hspace=0.2)
#        ax[0].axis('off')
#        ax[1].axis('off')
#
#
#        pts = G.vertex_properties["pos"].get_2d_array(range(2)).T
#        before.scatter(pts[:,0], pts[:, 1], marker="*")
#
#        #plot the connections of the top level graph
#        for e in G.edges():
#            coord = G.vertex_properties["pos"][e.source()]
#            coord2 = G.vertex_properties["pos"][e.target()]
#            x = [coord[0], coord2[0]]
#            y = [coord[1], coord2[1]]
#            before.plot(x, y, 'go--', linewidth=1, markersize=1)
#
#        #plot the connections of the top level graph
#        for e in G.edges():
#            coord = G.vertex_properties["pos"][e.source()]
#            coord2 = G.vertex_properties["pos"][e.target()]
#            x = [coord[0], coord2[0]]
#            y = [coord[1], coord2[1]]
#            before.plot(x, y, 'go--', linewidth=1, markersize=1)

velocities = np.zeros([G.num_vertices(), 3])
print("setting neighborhoods: ", time.time() - start)
start = time.time()

pos = G.vertex_properties["pos"].get_2d_array(range(3)).T
#x = np.zeros(3, dtype=float)

weight = G.vertex_properties[v_property].get_array().T

for n in range(n_steps+1):
forces = np.zeros([G.num_vertices(), 3])
center = np.sum(pos, 0)/G.num_vertices()
#pos = G.vertex_properties["pos"].get_2d_array(range(2)).T
if Sum:
for v in G.vertices():
for i in o1_n[int(v)]:
x = pos[int(i), :] - pos[int(v), :]
l = np.sqrt(x[0]*x[0] + x[1]*x[1])
forces[int(v), 0] += -k*(l-1.0)*(x[0])
forces[int(v), 1] += -k*(l-1.0)*(x[1])
for i in o2_n[int(v)]:
x = pos[int(i), :] - pos[int(v), :]
l = np.sqrt(x[0]*x[0] + x[1]*x[1])
forces[int(v), 0] += -k*(l-2.0)*(x[0])/2.0
forces[int(v), 1] += -k*(l-2.0)*(x[1])/2.0
for i in o3_n[int(v)]:
x = pos[int(i), :] - pos[int(v), :]
l = np.sqrt(x[0]*x[0] + x[1]*x[1])
forces[int(v), 0] += -k*(l-3.0)*(x[0])/4.0
forces[int(v), 1] += -k*(l-3.0)*(x[1])/4.0
x = -center[:] + pos[int(v), :]
l = np.sqrt(x[0]*x[0] + x[1]*x[1])
forces[int(v), 0] += -k/l*(x[0])
forces[int(v), 1] += -k/l*(x[1])

else:
for v in G.vertices():
for i in o3_n[int(v)]:
x = pos[int(i)] - pos[int(v)]
if np.sqrt(np.power(-k*(x[0])/4.0, 2) + np.sqrt(np.power(-k*(x[1])/4.0, 2))) > \
np.sqrt(np.power(forces[int(v), 0], 2) + np.sqrt(np.power(forces[int(v), 0], 2))):

forces[int(v), 0] = -k*(x[0])/4.0
forces[int(v), 1] = -k*(x[1])/4.0

x = pos[int(v), :] - center[:]
forces[int(v), 0] += -k*(x[0])/100.0
forces[int(v), 1] += -k*(x[1])/100.0
#
forces[:, 0] = forces[:, 0]/weight[:]
forces[:, 1] = forces[:, 1]/weight[:]
#            if n < 1:
#                velocities = forces * time_step*G.num_vertices()
#            else:
#                pos += velocities * time_step
#                velocities = forces * time_step * cooling
#                cooling -= cooling_step
#

for v in G.vertices():
if n < 1:
velocities = forces * time_step
#pos[int(v), :] = pos[int(v), :] + velocities[int(v), :] * time_step
#G.vertex_properties["pos"][v] = pos[int(v), :] + forces[int(v), :] * time_step
else:
pos[int(v), :] = pos[int(v), :] + velocities[int(v), :] * time_step
velocities = forces * time_step * cooling
cooling -= cooling_step

G.vertex_properties["pos"] = G.new_vertex_property("vector<double>", vals = pos)
print("calculating new positions: ", time.time() - start)
#        pts = G.vertex_properties["pos"].get_2d_array(range(2)).T
#        after.scatter(pts[:,0], pts[:, 1], marker="*")
#
#        #plot the connections of the top level graph
#        for e in G.edges():
#            coord = G.vertex_properties["pos"][e.source()]
#            coord2 = G.vertex_properties["pos"][e.target()]
#            x = [coord[0], coord2[0]]
#            y = [coord[1], coord2[1]]
#            after.plot(x, y, 'go--', linewidth=1, markersize=1)
#
#        plt.show()
return G

'''
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<double>")                          #original position of the edge (placed as midpoin of edge)
pos = G.new_vertex_property("vector<double>")                           #positions according to the fdl
rpos = G.new_vertex_property("vector<double>")                          #positions according to the radial layout
x_edge = G.new_vertex_property("vector<double>")                          #x-coordinated of the edge
y_edge = G.new_vertex_property("vector<double>")                          #y-coordinates of the edge
z_edge = G.new_vertex_property("vector<double>")                          #z-coordinates of the edge
r_edge = G.new_vertex_property("vector<double>")                          #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")

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):
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):
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.)))
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):
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

#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)

#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"
#
#        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"
#
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"
#        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<double>")                          #original vertex positions
pos = G.new_vertex_property("vector<double>")                           #positions according to the fdl
rpos = G.new_vertex_property("vector<double>")                          #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<double")
x_edge = G.new_edge_property("vector<double>")                          #x-coordinated of the edge
y_edge = G.new_edge_property("vector<double>")                          #y-coordinates of the edge
z_edge = G.new_edge_property("vector<double>")                          #z-coordinates of the edge
r_edge = G.new_edge_property("vector<double>")                          #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.
exclude = G.new_edge_property("bool", val = False)

#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)):
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()
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
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)

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)
G.vertex_properties["exclude"] = 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["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)
G.edge_properties["exclude"] = exclude

#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<double>")
G.graph_properties["point_cloud_y"] = G.new_graph_property("vector<double>")
G.graph_properties["point_cloud_z"] = G.new_graph_property("vector<double>")
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"
#
#        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"
#        #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"
#
#        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"
#        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(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):
G.vertex_properties["pos"] = gt.sfdp_layout(G, groups = G.vertex_properties["clusters"], pos = G.vertex_properties["pos"], C = 1.0, K = 10)
#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())
if mx == mn:
mn = 0
G.vertex_properties["map"] = G.new_vertex_property("float")
print("stuff ", mx, mn, propertymap)
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<double> 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<double>", 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<double>' or value_type == 'vector<float>'):
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<boolean>")
#ecluster_belongance = G.new_edge_property("vector<boolean>")
#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());
#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<boolean>", array_cluster_bel)
G.edge_properties["partition"] = G.new_edge_property("vector<boolean>", 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:
#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)
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)):
for i in range(len(self.F)):
G[self.F[i].v0][self.F[i].v1]['pts'] = self.F[i].points

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)):
for i in range(len(self.F)):
G[self.F[i].v0][self.F[i].v1]['pts'] = self.F[i].points

return G
#'''
#Returns the positions dictionary for the Circular layout.
#'''
#
#'''
#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
``````