segmentation / graph_gui

  # -*- 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('<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]))
              open_file.write(struct.pack('<f', edge.radii[j]))
              
          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:
              NWT.writeHeader(file, len(vertices), len(edges))
              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.
      '''
      def readVertex(open_file):
          points = np.tile(0., 3)
          bytes = open_file.read(4)
          points[0] = struct.unpack('f', bytes)[0]
          bytes = open_file.read(4)
          points[1] = struct.unpack('f', bytes)[0]
          bytes = open_file.read(4)
          points[2] = struct.unpack('f', bytes)[0]
          bytes = open_file.read(4)
          
          numO = int.from_bytes(bytes, byteorder='little')
          outgoing = np.tile(0, numO)
          bts = open_file.read(4)
          numI = int.from_bytes(bts, byteorder='little')
          incoming = np.tile(0, numI)
          for j in range(numO):
              bytes = open_file.read(4)
              outgoing[j] = int.from_bytes(bytes, byteorder='little')
              
          for j in range(numI):
              bytes = open_file.read(4)
              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 .   
      '''
      def readFiber(open_file):
          bytes = open_file.read(4)
          vtx0 = int.from_bytes(bytes, byteorder = 'little')
          bytes = open_file.read(4)
          vtx1 = int.from_bytes(bytes, byteorder = 'little')
          bytes = open_file.read(4)
          numVerts = int.from_bytes(bytes, byteorder = 'little')
          pts = []
          rads = []
          
          for j in range(numVerts):
              point = np.tile(0., 3)
              bytes = open_file.read(4)
              point[0] = struct.unpack('f', bytes)[0]
              bytes = open_file.read(4)
              point[1] = struct.unpack('f', bytes)[0]
              bytes = open_file.read(4)
              point[2] = struct.unpack('f', bytes)[0]
              bytes = open_file.read(4)
              radius = struct.unpack('f', bytes)[0]
              pts.append(point)
              rads.append(radius)
              
          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,1), 1000000.0, dtype=float)
          #maximum vertex
          self.B = np.full((3,1), -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, 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))