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))
        if DEBUG:
            print("num of edges to begin with = ", G.num_edges())
        indices = np.random.randint(0, int(G.num_edges()), size = [num_removed,1])
        #aabb = AABB(G, is_dual)
        tf = np.full((G.num_edges(), 1), True, dtype=bool)
        for j in range(indices.shape[0]):
            tf[indices[j]] = False
        TFe = G.new_edge_property("bool", vals = tf)
        G.set_edge_filter(TFe)
        G1 = gt.Graph(G, directed=False, prune=True)
        G.clear_filters()
        G1.clear_filters()
        G1 = self.filterDisconnected(G1)
        G1.vertex_properties["degree"] = G1.degree_property_map("total")
        if DEBUG:
            print("num of edges left = ", G1.num_edges())
            print("I should have removed, ", num_removed)
            print("Instead I removed, ", G.num_edges() - G1.num_edges())
        
        return G1
    
    '''
        Adds the mst based- sensitivity metrics. N is the number of random removals
        remove is the percentage*100 of the vessels removed in each iteration.
        The number of vessels removed will be floor(num_edges*remove/100)
    '''
    def add_rst_metric(self, G, N=100, remove=10, is_dual=False):
        #rst_ratio = G.new_graph_property("double")
        rst_r = 0
        #G.nu
        num_removed = int(np.floor(G.num_edges()*remove/100))
        indices = np.random.randint(0, int(G.num_edges()), size = [N, num_removed])
        aabb = AABB(G, is_dual)
        
        for i in range(N):
            tf = np.full((G.num_edges(),1), True, dtype=bool)
            for j in range(num_removed):
                tf[indices[i, j]] = False
#                rst = gt.graph_tool.topology.min_spanning_tree(G, weights = G.)
#                ratio = np.double(np.sum(rst.get_array()))/np.double(G.num_edges())
#                rst_dist.append(ratio)
#                rst_r = rst_r + ratio
            TFe = G.new_edge_property("bool", vals = tf)
            G.set_edge_filter(TFe)
            G1 = gt.Graph(G, directed=False, prune=True)
            G.clear_filters()
            G1.clear_filters()
            while(np.any(G1.degree_property_map("total").get_array() == True)):
                #print(np.any(G1.degree_property_map("total").get_array()), "        ", i)
                G1 = self.filterDisconnected(G1)
                G1 = self.filterBorder(G1, is_dual)
                G1 = self.recalculate_metrics(G1, is_dual)
                
            if(not is_dual and G1.num_edges() > 2):
                G1.edge_properties["mst"] = gt.graph_tool.topology.min_spanning_tree(G1, weights=G1.edge_properties["volume"])
                tvalues = G1.edge_properties["mst"].get_array()
                #print(tvalues)
                values = G1.edge_properties["inverted_volume"].get_array()
                #print(values)
                for k in range(G1.num_edges()):
                    if(tvalues[k] == True):
                        rst_r = rst_r + values[k]
            elif(is_dual and G1.num_edges() > 2):
                #This is non-sensical for now.
                G1.edge_properties["mst"] = gt.graph_tool.topology.min_spanning_tree(G1)
                tvalues = G1.vertex_properties["mst"].get_array()
                values = G1.edge_properties["inverted_volume"].get_array()
                for k in range(G1.num_edges()):
                    if(tvalues[k] == True):
                        rst_r = rst_r + values[k] 
                
        G.graph_properties["sensitivity"] = G.new_graph_property("double", rst_r/N)
        
        return G
    
    '''
        Re-calculates connectivity based metrics and returns the resulting graph
    '''
    def recalculate_metrics(self, G, is_dual=False, normalize=False):
        gt.graph_tool.centrality.betweenness(G, vprop=G.vertex_properties["bc"], eprop=G.edge_properties["bc"], norm=normalize)
        if(is_dual == False):
            #gt.graph_tool.centrality.betweenness(G, G.vertex_properties["bc"], G.edge_properties["bc"], norm=True, is_dual=False)
            #generate minimum spanning tree
            
            G.edge_properties["mst"] = gt.graph_tool.topology.min_spanning_tree(G, weights=G.edge_properties["length"])
            G.graph_properties["mst_ratio"] = np.double(np.sum(G.edge_properties["mst"].get_array()))/np.double(G.num_edges())
            if DEBUG:
                print(np.double(np.sum(G.edge_properties["mst"].get_array()))/np.double(G.num_edges()))
            G.vertex_properties["degree"] = G.degree_property_map("total")
            G.vertex_properties["degree_volume"] = G.degree_property_map("total", weight=G.edge_properties["volume"])
            G.vertex_properties["degree_tortuosity"] = G.degree_property_map("total", G.edge_properties["tortuosity"])
     
        else:
            #gt.graph_tool.centrality.betweenness(G, G.vertex_properties["bc"], G.edge_properties["bc"], norm=True, is_dual=False)
            G.edge_properties["mst"] = gt.graph_tool.topology.min_spanning_tree(G, weights=G.edge_properties["bc"])                                       #Boolean value deligating belongance to the minimum spanning tree.
            G.graph_properties["mst_ratio"] = np.double(np.sum(G.edge_properties["mst"].get_array()))/np.double(G.num_edges())               
            G.vertex_properties["degree"] = G.degree_property_map("total")
            G.vertex_properties["degree_centrality"] = G.degree_property_map("total", weight=G.edge_properties["bc"])

        return G
        
    '''
    Filter the graph to remove the disconnected components of the graph if
        disconnected is set to true. if borders are set to True, the algorithm 
        cleans up the nodes with degree one that are close to the edge of the volume.
        Returns a Graph object with updated internal property maps.
        
        if return_largerst == False, will return an array of graphs.
        disconnected == apply the disconnected componenet filter
        return_largest == if false returns the largest set of graphs (n > 10)
        borders == True, if true applies the border (deg == 1 near border) filter
        add_rst == True, if true adds the rst_metric
    '''
    def filterFullGraph_gt(self, G, disconnected=True, return_largest=True, borders=True, erode = False, add_rst=False, dual=False):
        if(disconnected==True):
            if(return_largest==True and borders==True):
                G1 = self.filterDisconnected(G, return_largest)
                G1 = self.recalculate_metrics(G1, is_dual=dual)
                G1 = self.filterBorder(G1, is_dual=dual)
                if(erode == True):
                    while(np.any(G1.degree_property_map("total").get_array() == True)):
                        G1 = self.filterBorder(G1, is_dual=dual)
                        G1 = self.recalculate_metrics(G1, is_dual=dual)
                else:
                    G1 = self.recalculate_metrics(G1, is_dual=dual)
                if(add_rst == True):
                    G1 = self.add_rst_metric(G1)
            elif(return_largest==True and borders==False):
                G1 = self.filterDisconnected(G, return_largest)
                if(erode == True):
                    while(np.any(G1.degree_property_map("total").get_array() == True)):
                        G1 = self.filterBorder(G1, is_dual=dual)
                        G1 = self.recalculate_metrics(G1, is_dual=dual)
                else:
                    G1 = self.recalculate_metrics(G1, is_dual=dual)
                if(add_rst == True):
                    G1 = self.add_rst_metric(G1)
            elif(return_largest==False and borders == True):
                G1 = self.filterDisconnected(G, return_largest)
                for graph in G1:
                    graph = self.recalculate_metrics(graph, is_dual=dual)
                    graph = self.filterBorder(graph, is_dual=dual)
                    if(erode == True):
                        while(np.any(graph.degree_property_map("total").get_array() == True)):
                            graph = self.filterBorder(graph, is_dual=dual)
                            graph = self.recalculate_metrics(graph, is_dual=dual)
                    else:
                        graph = self.recalculate_metrics(graph, is_dual=dual)
                    if(add_rst == True):
                        graph = self.add_rst_metric(graph)
            else:
                G1 = self.filterDisconnected(G, return_largest)
                for graph in G1:
                    graph = self.recalculate_metrics(graph, is_dual=dual)
                    if(add_rst == True):
                        graph = self.add_rst_metric(graph)
        else:
            G1 = self.filterBorder(G, is_dual=dual);
            if(erode == True):
                while(np.any(G1.degree_property_map("total").get_array() == True)):
                    if DEBUG:
                        print(G1.num_vertices())
                    G1 = self.filterBorder(G1, is_dual=dual)
                    if DEBUG:
                        print(G1.num_vertices())
                    G1 = self.recalculate_metrics(G1, is_dual=dual, )
            else:
                G1 = self.recalculate_metrics(G1, is_dual=dual)
            if(add_rst==True):
                G1 = self.add_rst_metric(G1)
        
        return G1
    
    '''
        Accretes the graph in G using the nodes and edges in G_0
        Assumes that G_0 has all the same properties as G, including "idx"
        idx is an immutable set of indices that do not change when the graph is
        modified
    '''
    def accrete(self, G, G_0):
        
        v_filters = G_0.new_vertex_property("bool", vals = np.full((G_0.num_vertices(),1), False, dtype=bool))
        e_filters = G_0.new_edge_property("bool", vals = np.full((G_0.num_edges(),1), False, dtype=bool))
        for v in G.vertices():
            #print(v)
            G_0_vertex = G_0.vertex(G.vertex_properties["idx"][v])
            v_filters[G_0_vertex] = True
            for e in G_0_vertex.all_edges():
                e_filters[e] = True
            for v_0 in G_0_vertex.all_neighbors():
                v_filters[v_0] = True
                
        G_0.set_filters(e_filters, v_filters)
        
        graph = gt.Graph(G_0, directed=False, prune=True)
        G_0.clear_filters()
        
        return graph
            
            
    
    
    '''
        Creates a graph from a list of nodes and a list of edges
        The graph is the "inverse" of the original graph,
            Vertices are edges and edges are 
        
    '''
    def createDualGraph_gt(self):
        #Generate the undirected graph
        G = gt.Graph()
        G.set_directed(False)
        
        #add all the required vertex properties
        vpos = G.new_vertex_property("vector<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")
        
    
        #add all graph properties
        mst_ratio = G.new_graph_property("double")

        #add all the edges as vertices
        for i in range(len(self.F)):
            if(self.F[i].length() > 0):
                G.add_vertex()
                x,y,z,r = self.F[i].getcoords()
                x_edge[G.vertex(i)] = x
                y_edge[G.vertex(i)] = y
                z_edge[G.vertex(i)] = z
                r_edge[G.vertex(i)] = r
                l_edge[G.vertex(i)] = self.F[i].length()
                gaussian[G.vertex(i)] = np.exp(-np.power(l_edge[G.vertex(i)], 2.) / (2 * np.power(60.0, 2.)))
                i_edge[G.vertex(i)] = 1/l_edge[G.vertex(i)]
                t_edge[G.vertex(i)] = self.F[i].turtuosity()
                v_edge[G.vertex(i)] = self.F[i].volume()
                iv_edge[G.vertex(i)] = 1/v_edge[G.vertex(i)]
                vpos[G.vertex(i)] = np.array([x[int(len(x)/2)], y[int(len(x)/2)], z[int(len(x)/2)]], dtype=float)
                
        #based on the neighbors add edges. If two edges share a vertex, there is a vertex between them.
        for i in range(len(self.N)):
            #in case there are 3 or more edges attached to a vertex
            if(len(self.N[i].o) > 2):
                for j in range(len(self.N[i].o)-1):
                    G.add_edge(G.vertex(self.N[i].o[j]), G.vertex(self.N[i].o[j+1]))
                    l_edge[G.edge(G.vertex(self.N[i].o[j])), G.vertex(self.N[i].o[j+1])] = 1
                    gaussian[G.edge(G.vertex(self.N[i].o[j])), G.vertex(self.N[i].o[j+1])] = \
                        np.exp(-np.power(l_edge[G.edge(G.vertex(self.N[i].o[j]),G.vertex(self.N[i].o[j+1]))], 2.) / (2 * np.power(60.0, 2.)))
                #add the final edge.        
                G.add_edge(G.vertex(self.N[i].o[0]), G.vertex(self.N[i].o[len(self.N[i].o)-1]))
                l_edge[G.edge(G.vertex(self.N[i].o[j])), G.vertex(self.N[i].o[j+1])] = 1
                gaussian[G.edge(G.vertex(self.N[i].o[j])), G.vertex(self.N[i].o[j+1])] = \
                np.exp(-np.power(l_edge[G.edge(G.vertex(self.N[i].o[j]),G.vertex(self.N[i].o[j+1]))], 2.) / (2 * np.power(60.0, 2.)))
            #in case there are 2 edges attached to a vertex
            elif(len(self.N[i].o) == 2):
                G.add_edge(G.vertex(self.N[i].o[0]), G.vertex(self.N[i].o[1]-1))
                l_edge[G.edge(G.vertex(self.N[i].o[j])), G.vertex(self.N[i].o[j+1])] = 1
                gaussian[G.edge(G.vertex(self.N[i].o[j])), G.vertex(self.N[i].o[j+1])] = \
                    np.exp(-np.power(l_edge[G.edge(G.vertex(self.N[i].o[j]),G.vertex(self.N[i].o[j+1]))], 2.) / (2 * np.power(60.0, 2.)))
            #in case there is only a single edge coming into the vertex we do 
            #   not add anything
            
        #generate centrality map
        gt.graph_tool.centrality.betweenness(G, vbetweeness_centrality, ebetweeness_centrality, norm=True)
        #generate minimum spanning tree
        mst = gt.graph_tool.topology.min_spanning_tree(G, weights=ebetweeness_centrality)
        mst_ratio[G] = np.double(np.sum(mst.get_array()))/np.double(len(self.N))
        degree = G.degree_property_map("total")
        degree_centrality = G.degree_property_map("total", weight=ebetweeness_centrality)
        #degree_tortuosity = G.degree_property_map("total", weight=t_edge)
                
        pos = gt.sfdp_layout(G)
        
        rpos = gt.radial_tree_layout(G, root=G.vertex(np.argmax(vbetweeness_centrality.get_array())))
        
        #set property maps for the vertices
        G.vertex_properties["p"] = vpos
        G.vertex_properties["pos"] = pos
        G.vertex_properties["rpos"] = rpos
        G.vertex_properties["bc"] = vbetweeness_centrality
        G.vertex_properties["degree_centrality"] = degree_centrality
        G.vertex_properties["degree"] = degree
        G.vertex_properties["x"] = x_edge
        G.vertex_properties["y"] = y_edge
        G.vertex_properties["z"] = z_edge
        G.vertex_properties["r"] = r_edge
        G.vertex_properties["inverted"] = i_edge
        G.vertex_properties["inverted_volume"] = iv_edge
        G.vertex_properties["length"] = l_edge
        G.vertex_properties["tortuosity"] = t_edge
        G.vertex_properties["volume"] = v_edge
        G.vertex_properties["bc"] = vbetweeness_centrality

        #G.vertex_properties["pos"] = gt.fruchterman_reingold_layout(G, weight=G.edge_properties["length"], circular = True, n_iter= 10000, pos=G.vertex_properties["pos"], r = 2.0)
        G.vertex_properties["pos"] = gt.sfdp_layout(G, C = 0.5, p = 3.0)
        #set property maps for the edges
        G.edge_properties["bc"] = ebetweeness_centrality
        G.edge_properties["mst"] = mst
        
        #set graph properies
        G.graph_properties["mst_ratio"] = mst_ratio

                
        
#        title = "raw_dual_cortical_7_6_before_delete.pdf"
#        title2 = "raw_dual_radial_cordical_7_6_before_delete.pdf"
#
#        gt.graph_draw(G, pos=gt.radial_tree_layout(G, root=G.vertex(np.argmax(vbetweeness_centrality.get_array()))), edge_pen_width = 2.0, vertex_size=degree, edge_color=mst, vertex_fill_color=vbetweeness_centrality, output=title2, bg_color=[0.0, 0.0,0.0,1.0], vertex_text=G.vertex_index, output_size=(1000,1000),vertex_font_size = 6)
#        gt.graph_draw(G, pos=pos, edge_pen_width = 4.0, vertex_size=degree, edge_color=mst, vertex_fill_color=vbetweeness_centrality, output=title, bg_color=[0.0, 0.0,0.0,1.0], vertex_text=G.vertex_index, output_size=(1000,1000),vertex_font_size = 6)      
#        
        title = "./Original_2D_Graph.png"
#        title2 = "raw_dual_radial_cordical_7_6_after_delete.pdf"
#        
        G1 = self.filterFullGraph_gt(G, borders=False, dual=True)

#        #print(clusters.get_array()[:], clusters.get_array()[:])
#        #pos = gt.sfdp_layout(G)
#        
#        #gt.graph_draw(G, pos=,edge_color=mst, vertex_fill_color=vbetweeness_centrality, edge_pen_width=ebetweeness_centrality, output="raw.pdf")
#        
#        gt.graph_draw(G1, pos=gt.radial_tree_layout(G1, root=G1.vertex(np.argmax(G1.vertex_properties["bc"].get_array()))), edge_pen_width = 4.0, vertex_size=G1.vertex_properties["degree"], edge_color=G1.edge_properties["mst"], vertex_fill_color=G1.vertex_properties["bc"], output=title2, bg_color=[0.0, 0.0,0.0,1.0], vertex_text=G1.vertex_index, output_size=(1000,1000),vertex_font_size = 6)
#        gt.graph_draw(G1, pos=G1.vertex_properties["p"], edge_pen_width = 4.0, vertex_size=G1.vertex_properties["degree"], edge_color=G1.edge_properties["mst"], vertex_fill_color=G1.vertex_properties["bc"], output=title, bg_color=[0.0, 0.0,0.0,1.0], vertex_text=G1.vertex_index, output_size=(1000,1000),vertex_font_size = 6)
        gt.graph_draw(G1, pos=G1.vertex_properties["p"], edge_pen_width = 4.0, vertex_size=G1.vertex_properties["degree"], output=title, bg_color=[1.0, 1.0,1.0,1.0], vertex_text=G1.vertex_index, output_size=(1000,1000),vertex_font_size = 6)
        title = "./Original_2D_Layout.png"
        gt.graph_draw(G1, pos=G1.vertex_properties["p"], edge_pen_width = 4.0, vertex_size=G1.vertex_properties["degree"], edge_color=G1.edge_properties["mst"], vertex_fill_color=G1.vertex_properties["bc"], output=title, bg_color=[1.0, 1.0,1.0,1.0], vertex_text=G1.vertex_index, output_size=(1000,1000), vertex_font_size = 6)
#gt.graph_draw(G1, pos=G1.vertex_properties["pos"], edge_pen_width = 4.0, vertex_size=degree, edge_color=G1.edge_properties["mst"], vertex_fill_color=G1.vertex_properties["bc"], output=title)
#        
#        
#        G2 = self.filterFullGraph_gt(G1, add_rst=False, dual=True)
#        title = "raw_dual_cortical_7_6_after_borders.pdf"
#        title2 = "raw_dual_radial_cortical_7_6_after_borders.pdf"
#        gt.graph_draw(G2, pos=G2.vertex_properties["pos"], edge_pen_width = 4.0, vertex_size=G2.vertex_properties["degree"], edge_color=G2.edge_properties["mst"], vertex_fill_color=G2.vertex_properties["bc"], output=title, bg_color=[0.0, 0.0,0.0,1.0], vertex_text=G2.vertex_index, output_size=(1000,1000),vertex_font_size = 6)
#        gt.graph_draw(G2, pos=gt.radial_tree_layout(G2, root=G2.vertex(np.argmax(G2.vertex_properties["bc"].get_array()))), edge_pen_width = 4.0, vertex_size=G2.vertex_properties["degree"], edge_color=G2.edge_properties["mst"], vertex_fill_color=G2.vertex_properties["bc"], output=title2, bg_color=[0.0, 0.0,0.0,1.0], vertex_text=G2.vertex_index, output_size=(1000,1000),vertex_font_size = 6)
        return G
            
    
    '''
    Create a graph from a list of nodes and a list of edges.
    Populates the Graph with a bunch of property maps in respect to edges and nodes
    each vertex is assigned 3 positional coordiantes:
        The original x,y,z location of the graph.
        x,y locations according to force directed layout,
        and x,y positions according to a radial layout
    Returns a graph_tool object
    '''
    def createFullGraph_gt(self):
        #Generate the undirected graph
        G = gt.Graph()
        G.set_directed(False)
        
        #add all the required vertex properties
        vpos = G.new_vertex_property("vector<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.

        
        #This map gets updated.
        gaussian = G.new_edge_property("double")                                #empty property map for gaussian values downstream
        
        
        #Graph properties (once per region)
        mst_ratio = G.new_graph_property("double")                              #Graph based metric of mst edges/total edges.                              

        
        #add verticies and set their properties
        for i in range(len(self.N)):
            G.add_vertex()
            vpos[G.vertex(i)] = self.N[i].p
        
        #add edges and set their properties.
        for i in range(len(self.F)):
            v0 = self.F[i].v0
            v1 = self.F[i].v1
            if self.F[i].length() > 0.0:
                x,y,z,r = self.F[i].getcoords()
                G.add_edge(G.vertex(v0), G.vertex(v1))
                l_edge[G.edge(v0,v1)] = self.F[i].length()
                gaussian[G.edge(v0,v1)] = np.exp(-np.power(l_edge[G.edge(v0,v1)], 2.) / (2 * np.power(60.0, 2.)))
                i_edge[G.edge(v0,v1)] = 1/l_edge[G.edge(v0,v1)]
                t_edge[G.edge(v0,v1)] = self.F[i].turtuosity()
                v_edge[G.edge(v0,v1)] = self.F[i].volume()
                iv_edge[G.edge(v0,v1)] = 1/v_edge[G.edge(v0,v1)]
                x_edge[G.edge(v0,v1)] = x
                y_edge[G.edge(v0,v1)] = y
                z_edge[G.edge(v0,v1)] = z
                r_edge[G.edge(v0,v1)] = r
                av_edge[G.edge(v0,v1)] = self.F[i].av_radius()
            else:
                print("NON-CRITICAL ERROR: edge with length 0 detected--SKIPPED")      
        
        
        #generate centrality map
        gt.graph_tool.centrality.betweenness(G, vprop=vbetweeness_centrality, eprop=ebetweeness_centrality, norm=True)
        #generate minimum spanning tree
        mst = gt.graph_tool.topology.min_spanning_tree(G, weights=l_edge)
        mst_ratio[G] = np.double(np.sum(mst.get_array()))/np.double(len(self.N))
        degree = G.degree_property_map("total")
        degree_volume = G.degree_property_map("total", weight=v_edge)
        
        dg = degree_volume.get_array()
        dg = np.exp(-np.power(dg, 2.) / (2 * np.power(0.000005, 2.)))
        degree_volume = G.new_vertex_property("double", vals=dg)
        degree_tortuosity = G.degree_property_map("total", weight=t_edge)
        
        
        #print(clusters.get_array()[:], clusters.get_array()[:])
        pos = gt.sfdp_layout(G, C = 1.0, K = 10)
        rpos = gt.radial_tree_layout(G, root=G.vertex(np.argmax(vbetweeness_centrality.get_array())))
        
        
        
        for e in G.edges():
            ebc_length[G.edge(e.source(), e.target())] = ebetweeness_centrality[G.edge(e.source(), e.target())]*l_edge[G.edge(e.source(), e.target())]
        
        #set property maps for the vertices
        G.vertex_properties["p"] = vpos
        G.vertex_properties["pos"] = pos
        G.vertex_properties["rpos"] = rpos
        G.vertex_properties["bc"] = vbetweeness_centrality
        G.vertex_properties["degree"] = degree
        G.vertex_properties["degree_volume"] = degree_volume
        G.vertex_properties["degree_tortuosity"] = degree_tortuosity
        G.vertex_properties["selection"] = G.new_vertex_property("bool", val=False)
        
        #set property maps for the edges
        G.edge_properties["x"] = x_edge
        G.edge_properties["y"] = y_edge
        G.edge_properties["z"] = z_edge
        G.edge_properties["r"] = r_edge
        G.edge_properties["inverted"] = i_edge
        G.edge_properties["length"] = l_edge
        G.edge_properties["tortuosity"] = t_edge
        G.edge_properties["av_radius"] = av_edge 
        G.edge_properties["volume"] = v_edge
        G.edge_properties["bc"] = ebetweeness_centrality
        G.edge_properties["bc*length"] = ebc_length
        G.edge_properties["mst"] = mst
        G.edge_properties["inverted_volume"] = iv_edge
        G.edge_properties["gaussian"] = gaussian
        G.edge_properties["selection"] = G.new_edge_property("double", val=0.0)
        
        #set graph properies
        G.graph_properties["mst_ratio"] = mst_ratio
        
        N = 0
        for e in G.edges():
            N += len(G.edge_properties["x"][e])
        point_cloud_x = np.zeros((N, 1), dtype=np.float64)
        point_cloud_y = np.zeros((N, 1), dtype=np.float64)
        point_cloud_z = np.zeros((N, 1), dtype=np.float64)
        
        n = 0
        for e in G.edges():
            for p in range(len(G.edge_properties["x"][e])):
                point_cloud_x[n][0] = G.edge_properties["x"][e][p]
                point_cloud_y[n][0] = G.edge_properties["y"][e][p]
                point_cloud_z[n][0] = G.edge_properties["z"][e][p]
                n += 1
        
        G.graph_properties["point_cloud_x"] = G.new_graph_property("vector<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"
#        title2 = "raw_dual_radial_cordical_7_6_after_delete.pdf"
#        
#        G1 = self.filterFullGraph_gt(G, borders=False, dual=False, add_rst = False)
##        
##        #print(clusters.get_array()[:], clusters.get_array()[:])
##        #pos = gt.sfdp_layout(G)
##        
##        #gt.graph_draw(G, pos=,edge_color=mst, vertex_fill_color=vbetweeness_centrality, edge_pen_width=ebetweeness_centrality, output="raw.pdf")
##        
##        gt.graph_draw(G1, pos=gt.radial_tree_layout(G1, root=G1.vertex(np.argmax(G1.vertex_properties["bc"].get_array()))), edge_pen_width = 4.0, vertex_size=G1.vertex_properties["degree"], edge_color=G1.edge_properties["mst"], vertex_fill_color=G1.vertex_properties["bc"], output=title2, bg_color=[0.0, 0.0,0.0,1.0], vertex_text=G1.vertex_index, output_size=(1000,1000),vertex_font_size = 6)
##        gt.graph_draw(G1, pos=G1.vertex_properties["p"], edge_pen_width = 4.0, vertex_size=G1.vertex_properties["degree"], edge_color=G1.edge_properties["mst"], vertex_fill_color=G1.vertex_properties["bc"], output=title, bg_color=[0.0, 0.0,0.0,1.0], vertex_text=G1.vertex_index, output_size=(1000,1000),vertex_font_size = 6)
#        gt.graphviz_draw(G1, pos = G1.vertex_properties["p"], ratio="fill",  size = (100,100), layout = None, pin = True, overlap=True, vsize=(G1.vertex_properties["degree_tortuosity"], 2.), vprops={"index":G1.vertex_index}, penwidth= 8.0, vcolor = G1.vertex_properties["bc"], vcmap = cm.get_cmap("inferno"), output=title)
#        #gt.graph_draw(G1, pos=G1.vertex_properties["p"], edge_pen_width = 8.0, output=title, bg_color=[1.0, 1.0,1.0,1.0], vertex_size=60, vertex_fill_color=G1.vertex_properties["bc"], vertex_text=G1.vertex_index, output_size=(3200,3200),vertex_font_size = 32)
#        title = "~/Pictures/Original_2D_Layout.png"
#        gt.graph_draw(G1, pos=G1.vertex_properties["pos"], edge_pen_width = 8.0, edge_color=G1.edge_properties["mst"], vertex_size=40, vertex_fill_color=G1.vertex_properties["bc"], output=title, bg_color=[1.0, 1.0,1.0,1.0], vertex_text=G1.vertex_index, output_size=(3200,3200), vertex_font_size = 32)
#        
#        
#        G1 = self.filterFullGraph_gt(G, borders=True, dual=False, erode=True, add_rst = False)
#        title = "~/Pictures/Original_2D_Graph_Eroded.png"
#        gt.graph_draw(G1, pos=G1.vertex_properties["p"], edge_pen_width = 8.0, output=title, bg_color=[1.0, 1.0,1.0,1.0], vertex_size=60, vertex_fill_color=G1.vertex_properties["bc"], vertex_text=G1.vertex_index, output_size=(3200,3200),vertex_font_size = 32)
#        title = "~/Pictures/Original_2D_Layout_Eroded.png"
#        gt.graph_draw(G1, pos=G1.vertex_properties["pos"], edge_pen_width = 8.0, edge_color=G1.edge_properties["mst"], vertex_size=40, vertex_fill_color=G1.vertex_properties["bc"], output=title, bg_color=[1.0, 1.0,1.0,1.0], vertex_text=G1.vertex_index, output_size=(3200,3200), vertex_font_size = 32)
              
#        title = "raw_full_cortical_7_6_before_delete.pdf"
#        title2 = "raw_full_radial_cordical_7_6_before_delete.pdf"
#        #G.graph_properties["rst_ratio"] = rst_ratio
#        #print("")
#        #print(G.graph_properties["rst_ratio"])
#        #print(G.edge_properties["bc"].get_array())
#        #print()
#        
#        #G.vertex_properties[("p","betweeness_centrality")] = [vpos, vbetweeness_centrality]
#        #G.edge_properties[("x", "y", "z", "r", "length", "tortuosity", "volume", "betweeness_centrality")] = []
#        #G.set_edge_filter(mst)
#        #G.set_edge_filter(mst)
#        
#        #Experimental code for drawing graphs.
#        #gt.graph_draw(G,pos=vpos,edge_color=mst, vertex_fill_color=vbetweeness_centrality, edge_pen_width=ebetweeness_centrality, output="raw.pdf")
#        #gt.graph_draw(G,pos=pos, vertex_size=degree, edge_color=mst, vertex_fill_color=vbetweeness_centrality, output=title, bg_color=[1.0,1.0,1.0])       
#        #gt.graph_draw(G,pos=gt.radial_tree_layout(G, root=G.vertex(np.argmax(vbetweeness_centrality.get_array())), node_weight=vbetweeness_centrality), vertex_size=degree, edge_color=mst, vertex_fill_color=vbetweeness_centrality, output=title2, bg_color=[1.0, 1.0,1.0,1.0], vertex_text=G.vertex_index, output_size=(1000,1000), vertex_font_size = 6)
#        #print(np.argmax(vbetweeness_centrality.get_array()))
#        
#        #np.savetxt("vbetweeness_centrality", vbetweeness_centrality.get_array())
#        gt.graph_draw(G,pos=gt.radial_tree_layout(G, root=G.vertex(np.argmax(vbetweeness_centrality.get_array()))), edge_pen_width = 2.0, vertex_size=degree, edge_color=mst, vertex_fill_color=vbetweeness_centrality, output=title2, bg_color=[0.0, 0.0,0.0,1.0], vertex_text=G.vertex_index, output_size=(1000,1000),vertex_font_size = 6)
#        gt.graph_draw(G, pos=pos, edge_pen_width = 4.0, vertex_size=degree, edge_color=mst, vertex_fill_color=vbetweeness_centrality, output=title, bg_color=[0.0, 0.0,0.0,1.0], vertex_text=G.vertex_index, output_size=(1000,1000),vertex_font_size = 6)      
#        
#
#        #gt.graph_draw(G,pos=vpos, efilt=mst, vertex_fill_color=vbetweeness_centrality, edge_pen_width=ebetweeness_centrality, output="mst.pdf")
#        #gt.graph_draw(G,pos=pos, efilt=mst, vertex_fill_color=vbetweeness_centrality, output="mst.pdf")
#        #need to create subgraphs!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
#        #rst_r = 0
#        #for i in range(1000):
#        #    rst = gt.graph_tool.topology.random_spanning_tree(G)
#        #    rst_r = rst_r + np.double(np.sum(rst.get_array()))/np.double(len(self.N))
#            
#        #rst_ratio[G] = rst_r/1000.0
#        
#        
#        
#        #G1 = gt.Graph(G)
#
#        #G1.purge_edges()
#        #G1.purge_vertices(in_place = True)
#        
#        
#        title = "raw_full_cortical_7_6_after_delete.pdf"
#        title2 = "raw_full_radial_cordical_7_6_after_delete.pdf"
#        
#        G1 = self.filterFullGraph_gt(G,borders=False)
#        
#        #print(clusters.get_array()[:], clusters.get_array()[:])
#        #pos = gt.sfdp_layout(G)
#        
#        #gt.graph_draw(G,pos=vpos,edge_color=mst, vertex_fill_color=vbetweeness_centrality, edge_pen_width=ebetweeness_centrality, output="raw.pdf")
#        
#        gt.graph_draw(pos = gt.sfdp_layout(G1), output=title)
#        gt.graph_draw(G1, pos=G1.vertex_properties["pos"], edge_pen_width = 4.0, vertex_size=G1.vertex_properties["degree"], edge_color=G1.edge_properties["mst"], vertex_fill_color=G1.vertex_properties["bc"], output=title, bg_color=[0.0, 0.0,0.0,1.0], vertex_text=G1.vertex_index, output_size=(1000,1000),vertex_font_size = 6)
#        #gt.graph_draw(G1, pos=G1.vertex_properties["pos"], edge_pen_width = 4.0, vertex_size=degree, edge_color=G1.edge_properties["mst"], vertex_fill_color=G1.vertex_properties["bc"], output=title)
#        
#        
#        G2 = self.filterFullGraph_gt(G1)
#        title = "raw_full_cortical_7_6_after_borders.pdf"
#        title2 = "raw_full_radial_cortical_7_6_after_borders.pdf"
#        gt.graph_draw(G2, pos=G2.vertex_properties["pos"], edge_pen_width = 4.0, vertex_size=G2.vertex_properties["degree"], edge_color=G2.edge_properties["mst"], vertex_fill_color=G2.vertex_properties["bc"], output=title, bg_color=[0.0, 0.0,0.0,1.0], vertex_text=G2.vertex_index, output_size=(1000,1000),vertex_font_size = 6)
#        gt.graph_draw(G2, pos=gt.radial_tree_layout(G2, root=G2.vertex(np.argmax(G2.vertex_properties["bc"].get_array()))), edge_pen_width = 4.0, vertex_size=G2.vertex_properties["degree"], edge_color=G2.edge_properties["mst"], vertex_fill_color=G2.vertex_properties["bc"], output=title2, bg_color=[0.0, 0.0,0.0,1.0], vertex_text=G2.vertex_index, output_size=(1000,1000),vertex_font_size = 6)
#        #gt.graph_draw(G,pos=pos, vertex_size=degree, edge_color=mst, vertex_fill_color=vbetweeness_centrality, output=title, bg_color=[1.0,1.0,1.0])
#        
#        #gt.graph_draw(G,pos=gt.radial_tree_layout(G, root=G.vertex(np.argmax(vbetweeness_centrality.get_array())), node_weight=vbetweeness_centrality), vertex_size=degree, edge_color=mst, vertex_fill_color=vbetweeness_centrality, output=title2, bg_color=[1.0, 1.0,1.0,1.0], vertex_text=G.vertex_index, output_size=(1000,1000), vertex_font_size = 6)
#        #print(np.argmax(vbetweeness_centrality.get_array()))
#        #np.savetxt("vbetweeness_centrality", vbetweeness_centrality.get_array())
#        #gt.graph_draw(G1, pos=gt.radial_tree_layout(G1, root=G1.vertex(np.argmax(G1.vertex_properties["bc"].get_array()))), edge_pen_width = 2.0, vertex_size=G1.vertex_properties["degree"], edge_color=G1.edge_properties["mst"], vertex_fill_color=G1.vertex_properties["bc"], output=title2, bg_color=[0.0, 0.0,0.0,1.0], vertex_text=G1.vertex_index, output_size=(1000,1000),vertex_font_size = 6)
#        #gt.graph_draw(G1, pos=gt.radial_tree_layout(G1, root=G1.vertex(np.argmax(G1.vertex_properties["bc"].get_array()))), edge_pen_width = 2.0, vertex_size=degree, edge_color=G1.edge_properties["mst"])
#        #gt.draw_hierarchy(gt.minimize_nested_blockmodel_dl(G1, deg_corr=False), output = "test.pdf")
#        #gt.graph_draw(G1,pos=gt.radial_tree_layout(G1, root=G1.vertex(np.argmax(G1.edge_properties["mst"].get_array()))), edge_pen_width = 2.0, vertex_size=degree, edge_color=mst, vertex_fill_color=vbetweeness_centrality, output=title2, bg_color=[0.0, 0.0,0.0,1.0], vertex_text=G.vertex_index, output_size=(1000,1000),vertex_font_size = 6)
#        print(gt.graph_tool.openmp_enabled(), gt.graph_tool.openmp_get_num_threads())
        
            #print(ratio, G2.graph_properties["mst_ratio"]) 
        #n, bins, patches = plt.hist(np.asarray(rst_dist), 50, normed=1, facecolor='green', alpha=0.75)

        # add a 'best fit' line
        #plt.grid(True)

        #plt.show()
        #self.partition(G, 5, 7, False)
        self.export_points(G)
        self.export_bb(G)
        return G
            
    
    def gen_new_fd_layout(G):
        pos = gt.sfdp_layout(G, C = 1.0, K = 10)
        G.vertex_properties["pos"] = pos
        return G
    
    def map_edges_to_range(G, rng, propertymap):
        def func(maximum, minimum, new_maximum, new_minimum, value):
            return ((value-minimum)/(maximum-minimum)*(new_maximum-new_minimum)+new_minimum)
        
        mx = max(G.edge_properties[propertymap].get_array())
        mn = min(G.edge_properties[propertymap].get_array())
        G.edge_properties["map"] = G.new_edge_property("float")
        gt.map_property_values(G.edge_properties[propertymap], G.edge_properties["map"], lambda x: func(mx, mn, rng[0], rng[1], x))
        
        return G.edge_properties["map"]


    def map_vertices_to_range(G, rng, propertymap):
        def func(maximum, minimum, new_maximum, new_minimum, value):
            return ((value-minimum)/(maximum-minimum)*(new_maximum-new_minimum)+new_minimum)
        
        mx = max(G.vertex_properties[propertymap].get_array())
        mn = min(G.vertex_properties[propertymap].get_array())
        G.vertex_properties["map"] = G.new_vertex_property("float")
        gt.map_property_values(G.vertex_properties[propertymap], G.vertex_properties["map"], lambda x: func(mx, mn, rng[0], rng[1], x))
        
        return G.vertex_properties["map"]    

    '''
    G and propertymap in G are passed in.
    Maps select property from range to color, works for int's floats and vec3
    and vec4 values
    Returns a vector<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());
                    temp_vertices = np.ma.masked_where(dist.get_array() <= m, dist.get_array()).mask
                    #print(temp_vertices)
                    temp_edges = np.zeros([G.num_edges()])
                    for vertex in np.nditer(np.where(temp_vertices == True)):
                        for edge in G.get_out_edges(vertex):
                            #print(edge)
                            temp_edges[edge[2]] = 1
                    array_cluster_bel[c,:] = temp_vertices[:]
                    array_ecluster_bel[c,:] = temp_edges[:]
                    c = c + 1
                    #print(G.num_vertices())
        #np.savetxt("clusters.txt", array_cluster_bel)
        G.vertex_properties["partition"] = G.new_vertex_property("vector<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:
            edge = G1.add_edge(e.target(), e.source())
            #G.edge_properties["mst"][e] = 1
            G1.edge_properties["path_weights"][edge] = 1000
            #G.set_edge_filter(G.edge_properties["mst"])
            #G1 = gt.Graph(G, directed=False, prune=True)
            #s = -1
            #t = -1
            #for edge in G1.edges():
            #    if(G1.vertex_properties["original_index"][edge.source()] == int(e.source())):
            #        s = int(edge.source())
            #    elif(G1.vertex_properties["original_index"][edge.target()] == int(e.target())):
            #        t = int(edge.target())
                    
            [Vl, El] = gt.shortest_path(G1, e.source(), e.target(), weights=G1.edge_properties["path_weights"])
            clusters = G1.new_vertex_property("int", vals=np.full((G1.num_vertices(), 1), 5, dtype=int))
            eclusters = G1.new_edge_property("int", vals = np.full((G1.num_edges(), 1), 3, dtype=int))
            G1.vertex_properties["cycle"] = clusters
            G1.edge_properties["cycle"] = eclusters
            if DEBUG:
                print("Number of vertices in Path is:", len(Vl))
            for v in Vl:
                G1.vertex_properties["cycle"][G1.vertex(v)] = 10
                if DEBUG:
                    print(str(v))
            #Create the arrays to be histogrammed later regarding every loop    
            length_total = 0
            volume_total = 0
            tortuosity_total = 0
            bc_total = 0
            bcl_total = 0
            num_edges = 1
            for v in range(len(Vl)-1):
                G1.edge_properties["cycle"][G1.edge(Vl[v], Vl[v+1])] = 10
                length_total = length_total + G1.edge_properties["length"][G1.edge(Vl[v], Vl[v+1])]
                volume_total = volume_total + G1.edge_properties["volume"][G1.edge(Vl[v], Vl[v+1])]
                tortuosity_total = tortuosity_total + G1.edge_properties["tortuosity"][G1.edge(Vl[v], Vl[v+1])]
                bc_total = bc_total + G1.edge_properties["bc"][G1.edge(Vl[v], Vl[v+1])]
                bcl_total = bcl_total + G1.edge_properties["bc*length"][G1.edge(Vl[v], Vl[v+1])]
                num_edges = num_edges + 1
            G1.edge_properties["cycle"][G1.edge(e.source(), e.target())] = 10
            length_total = length_total + G.edge_properties["length"][G.edge(e.source(), e.target())]
            volume_total = volume_total + G.edge_properties["volume"][G.edge(e.source(), e.target())]
            tortuosity_total = tortuosity_total + G.edge_properties["tortuosity"][G.edge(e.source(), e.target())]
            bc_total = bc_total + G.edge_properties["bc"][G.edge(e.source(), e.target())]
            bcl_total = bcl_total + G.edge_properties["bc*length"][G.edge(e.source(), e.target())]
            ppath = path + "Loop_Histograms.txt"
            f = open(ppath, "a+")
            f.write("%.15f\t" % length_total)
            f.write("%.15f\t" % volume_total)
            f.write("%.15f\t" % tortuosity_total)
            f.write("%.15f\t" % bc_total)
            f.write("%.15f\t" % bcl_total)
            f.write("%d\t\n" % num_edges)
            f.close()
            
            title = path+"cycle" + str(iteration)+".png"
            gt.graph_draw(G1, pos=G1.vertex_properties["pos"], output=title, edge_color='red', edge_pen_width=G1.edge_properties["cycle"], vertex_text=G1.vertex_index, output_size=(1920,1280), vertex_font_size = 4, vertex_fill_color = G1.vertex_properties["cycle"], vertex_size = 3)
            G1.remove_edge(edge)
            #G.clear_filters();
            #G.edge_properties["mst"][e] = 0
            #G.edge_properties["path_weights"][e] = 1000
            iteration = iteration + 1
    '''        
    returns an affinity matrix based on the property edge property, given as a string.
    '''       
    def get_affinity_matrix(G, edge_property, gaussian=False, sigma = None):
        affinity = gt.shortest_distance(G, weights=G.edge_properties[edge_property])
        affinity = affinity.get_2d_array(range(G.num_vertices()))
        if gaussian:
            if sigma == None:
                sig = 60.0
            else:
                sig = sigma
            for i in range(G.num_vertices()):
                for j in range(G.num_vertices()):
                    affinity[i,j] = np.exp(-np.power(affinity[i, j], 2.) / (2 * np.power(sig, 2.)))
        return affinity
                                
            
    '''
        Attempts to apporixmate the number of clusters in the graph by finding
        the approximate number of minimum graph cuts
        
        TO_DO: Program the gaussian and sigma parameters
    '''    
    def approximate_n_cuts(G, affinity_property):
        #L = gt.laplacian(G, weight = G.edge_properties[affinity_property])
        G1 = gt.Graph(G, directed=False)
        G1 = Network.filterDisconnected(G1, G1)
        L = Network.get_affinity_matrix(G1, affinity_property, gaussian=True, sigma=157)
        L = sp.sparse.csgraph.laplacian(L, normed=True)
        n_components = G1.num_vertices()
        eigenvalues, eigenvectors = eigsh(L, k=n_components, which = "LM", sigma=1.0, maxiter = 5000)
        
        #plt.figure()
        #plt.scatter(np.arange(len(eigenvalues)), eigenvalues)
        #plt.grid()
        #plt.show()
        
        count = sum(eigenvalues > 1.01)
        return count       
    '''
    Cluster the vectices in the graph based on a property passed to the clustering algorithm
    
    TO_DO: program the guassian and sigma parameters
    '''
    def spectral_clustering(G, affinity_property, gaussian = False, sigma = None, n_clusters = None, num_iters = 1000):
        if(n_clusters == None):
            n_c = Network.approximate_n_cuts(G, affinity_property);
        else:
            n_c = n_clusters
        sc = SpectralClustering(n_c, affinity='precomputed', n_init = num_iters)
        sc.fit(Network.get_affinity_matrix(G, affinity_property, gaussian = True, sigma=157))
        
        return sc.labels_ 
        
    
    def export_points(self, G):
        location = "./vertex_points.txt"
        f = open(location, "a+")
        for v in G.vertices():
            f.write("%.15f\t" % G.vertex_properties["p"][v][0])
            f.write("%.15f\t" % G.vertex_properties["p"][v][1])
            f.write("%.15f\n" % G.vertex_properties["p"][v][2])
            
        f.close()
            
    def export_bb(self, G):
        location = "./bb_points.txt"
        aabb = AABB(G)
        points = aabb.vertices()
        f = open(location, "a+")
        for i in points:
            f.write("%.15f\t" % i[0])
            f.write("%.15f\t" % i[1])
            f.write("%.15f\n" % i[2])
            
        f.close()
        
    
    def write_vtk(G, location, camera = None, binning = True):
        #location = "./nwt.vtk"
        f = open(location, "w+")
        f.write("# vtk DataFile Version 1.0\n")
        f.write("Data values\n")
        f.write("ASCII\n\n")
        f.write("DATASET POLYDATA\n")
        num_pts = 0
        for e in G.edges():
            X = G.edge_properties["x"][e]
            num_pts += len(X)
        f.write("POINTS %d float\n" % num_pts)
        for e in G.edges():
            X = G.edge_properties["x"][e]
            Y = G.edge_properties["y"][e]
            Z = G.edge_properties["z"][e]
            #pts = np.array([X,Y,Z]).T
            for p in range(len(X)):
                f.write("%.15f\t" % X[p])
                f.write("%.15f\t" % Y[p])
                f.write("%.15f\n" % Z[p])
        f.write("LINES " + str(G.num_edges()) + " " + str(G.num_edges()+num_pts) + "\n")
        num_pts = 0
        for e in G.edges():
            X = G.edge_properties["x"][e]
            indices = list(range(0, len(X)))
            indices = [x+num_pts for x in indices]
            f.write(str(len(X)) + " ")
            f.write(" ".join(str(x) for x in indices))
            f.write("\n")
            num_pts += len(X)
        f.write("POINT_DATA %d\n" % num_pts)
        f.write("SCALARS radius float 1\n")
        f.write("LOOKUP_TABLE radius_table\n")
        for e in G.edges():
            R = G.edge_properties["r"][e]
            #pts = np.array([X,Y,Z]).T
            for p in range(len(R)):
                f.write("%.15f\n" % R[p])

        f.write("COLOR_SCALARS color_table %d\n" % 4)
        bins = [0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0]
        
        if camera == None:
            for e in G.edges():
                X = G.edge_properties["x"][e]
                #RGBA = G.edge_properties["RGBA"].get_2d_array(range(4)).T
                for p in range(len(X)):
                    f.write("%.15f\t" % G.edge_properties["RGBA"][e][0])
                    f.write("%.15f\t" % G.edge_properties["RGBA"][e][1])
                    f.write("%.15f\t" % G.edge_properties["RGBA"][e][2])
                    if binning:
                        index = np.digitize(G.edge_properties["RGBA"][e][3], bins, right=True)
                        if (index >= len(bins) or index < 0):
                            if DEBUG:
                                print(G.edge_properties["RGBA"][e][3])
                                print(index)
                        f.write("%.15f\n" % bins[index])
                    else:
                        f.write("%.15f\n" % G.edge_properties["RGBA"][e][3])
        else:
            ptx = G.graph_properties["point_cloud_x"].get_array()
            pty = G.graph_properties["point_cloud_y"].get_array()
            ptz = G.graph_properties["point_cloud_z"].get_array()
            pts = np.array([ptx,pty,ptz]).T
            temp = AABB(G)
            bbl = temp.A
            bbu = temp.B
            tp = np.zeros((1,3), dtype = np.float64)
            tp[0][0] = (bbu[0]-bbl[0])/2
            tp[0][1] = (bbu[1]-bbl[1])/2
            tp[0][2] = (bbu[2]-bbl[2])/2
            pts[:] = pts[:] - camera - \
            tp
            distance = np.zeros(pts.shape[0], dtype = np.float32)
            for i in range(distance.shape[0]):
                distance[i] = np.sqrt(np.power(pts[i][0],2) + np.power(pts[i][1],2) + np.power(pts[i][2],2))
            mx = max(distance)
            mn = min(distance)
            distance = (distance - mx)/(mn-mx)
            idx = 0
            for e in G.edges():
                X = G.edge_properties["x"][e]
                for p in range(len(X)):
                    f.write("%.15f\t" % G.edge_properties["RGBA"][e][0])
                    f.write("%.15f\t" % G.edge_properties["RGBA"][e][1])
                    f.write("%.15f\t" % G.edge_properties["RGBA"][e][2])
                    if binning:
                        index = np.digitize(distance[idx], bins, right=True)
                        f.write("%.15f\n" % bins[index])
                    else:
                        f.write("%.15f\n" % distance[idx])
                    idx += 1
                
        f.close()
        
        
    '''
    Creates a graph from a list of nodes and a list of edges.
    Uses edge turtuosity as weight.
    Returns a NetworkX Object.
    '''    
    def createTortuosityGraph(self):
        G = nx.Graph()
        for i in range(len(self.N)):
            G.add_node(i, p=self.N[i].p)
        for i in range(len(self.F)):
            G.add_edge(self.F[i].v0, self.F[i].v1, weight = self.F[i].turtuosity())
            G[self.F[i].v0][self.F[i].v1]['pts'] = self.F[i].points
            G[self.F[i].v0][self.F[i].v1]['rads'] = self.F[i].radii
            
        return G

    
    '''
    Creates a graph from a list of nodes and a list of edges.
    Uses edge volume as weight.
    Returns a NetworkX Object.
    '''    
    def createVolumeGraph(self):
        G = nx.Graph()
        for i in range(len(self.N)):
            G.add_node(i, p=self.N[i].p)
        for i in range(len(self.F)):
            G.add_edge(self.F[i].v0, self.F[i].v1, weight = self.F[i].volume())
            G[self.F[i].v0][self.F[i].v1]['pts'] = self.F[i].points
            G[self.F[i].v0][self.F[i].v1]['rads'] = self.F[i].radii
            
        return G
#'''
#Returns the positions dictionary for the Circular layout.
#'''    
#def getCircularLayout(graph, dim, radius):
#    return nx.circular_layout(graph, dim, radius)
#
#'''
#Return the positions dictionary for the Spring layout.
#'''    
#def getSpringLayout(graph, pos, iterations, scale):
#    return nx.spring_layout(graph, 2, None, pos, iterations, 'weight', scale, None)
#        
#'''
#Draws the graph.
#'''        
#def drawGraph(graph, pos):
#    nx.draw(graph, pos)
#    return

    def aabb(self):
    
        lower = self.N[0].p.copy()
        upper = lower.copy()
        for i in self.N:
            for c in range(len(lower)):
                if lower[c] > i.p[c]:
                    lower[c] = i.p[c]
                if upper[c] < i.p[c]:
                    upper[c] = i.p[c]
        return lower, upper
    
    #calculate the distance field at a given resolution
    #   R (triple) resolution along each dimension
    def distancefield(self, R=(100, 100, 100)):
        
        #get a list of all node positions in the network
        P = []
        for e in self.F:
            for p in e.points:
                P.append(p)
                
        #turn that list into a Numpy array so that we can create a KD tree
        P = np.array(P)
        
        #generate a KD-Tree out of the network point array
        tree = sp.spatial.cKDTree(P)
        
        #plt.scatter(P[:, 0], P[:, 1])
        
        #specify the resolution of the ouput grid
        R = (200, 200, 200)
        
        #generate a meshgrid of the appropriate size and resolution to surround the network
        lower, upper = self.aabb(self.N, self.F)    #get the space occupied by the network
        x = np.linspace(lower[0], upper[0], R[0])   #get the grid points for uniform sampling of this space
        y = np.linspace(lower[1], upper[1], R[1])
        z = np.linspace(lower[2], upper[2], R[2])
        X, Y, Z = np.meshgrid(x, y, z)
        #Z = 150 * numpy.ones(X.shape)
        
        
        Q = np.stack((X, Y, Z), 3)
        
        
        D, I = tree.query(Q)
        
        return D
    
    #returns the number of points in the network
    def npoints(self):                              
        n = 0                                       #initialize the counter to zero
        for f in self.F:                            #for each fiber
            n = n + len(f.points) - 2               #count the number of points in the fiber - ignoring the end points
        n = n + len(self.N)                         #add the number of nodes (shared points) to the node count
        return n                                    #return the number of nodes
    
    #returns all of the points in the network
    def points(self):
        k = self.npoints()
        P = np.zeros((3, k))                        #allocate space for the point list
        
        idx = 0
        for f in self.F:                            #for each fiber in the network
            for ip in range(1, len(f.points)-1):    #for each point in the network
                P[:, idx] = f.points[ip]            #store the point in the raw point list
                idx = idx + 1
        return P                                    #return the point array        
    
    #returns the number of linear segments in the network
    def nsegments(self):
        n = 0                                       #initialize the segment counter to 0
        for f in self.F:                            #for each fiber
            n = n + len(f.points) - 1               #calculate the number of line segments in the fiber (points - 1)
        return n                                    #return the number of line segments
    
    #return a list of line segments representing the network
    def segments(self, dtype=np.float32):
        k = self.nsegments()                        #get the number of line segments
        start = np.zeros((k, 3),dtype=dtype)                    #start points for the line segments
        end = np.zeros((k, 3), dtype=dtype)                      #end points for the line segments
        
        idx = 0                                     #initialize the index counter to zero
        for f in self.F:                            #for each fiber in the network
            for ip in range(0, len(f.points)-1):    #for each point in the network
                start[idx, :] = f.points[ip]            #store the point in the raw point list
                idx = idx + 1
        
        idx = 0
        for f in self.F:                            #for each fiber in the network
            for ip in range(1, len(f.points)):      #for each point in the network
                end[idx, :] = f.points[ip]            #store the point in the raw point list
                idx = idx + 1
                
        return start, end
    
    #returns the points inside the fiber.
    def fiber(self, idx):
        p = self.F[idx].points
        X = np.zeros(len(p))
        Y = np.zeros(len(p))
        Z = np.zeros(len(p))
        for i in range(len(p)):
            X[i] = p[i][0]
            Y[i] = p[i][1]
            Z[i] = p[i][2]
        return X, Y, Z
    
    #function returns the fiber associated with a given 1D line segment index
    def segment2fiber(self, idx):        
        i = 0
        for f in range(len(self.F)):                #for each fiber in the network
            i = i + len(self.F[f].points)-1         #add the number of points in the fiber to i
            if i > idx:                             #if we encounter idx in this fiber
                return self.F[f].points, f          #return the fiber associated with idx and the index into the fiber array
        
    def vectors(self, clock=False, dtype=np.float32):
        if clock:
            start_time = time.time()
        start, end = self.segments(dtype)                #retrieve all of the line segments
        v = end - start                             #calculate the resulting vectors
        l = np.sqrt(v[:, 0]**2 + v[:,1]**2 + v[:,2]**2) #calculate the fiber lengths
        z = l==0                                    #look for any zero values
        nz = z.sum()
        if nz > 0:
            print("WARNING: " + str(nz) + " line segment(s) of length zero were found in the network and will be removed" )
            
        if clock:
            print("Network::vectors: " + str(time.time() - start_time) + "s")
            
        return np.delete(v, np.where(z), 0)
    
    #scale all values in the network by tuple S = (sx, sy, sz)
    def scale(self, S):
        for f in self.F:
            for p in f.points:
                p[0] = p[0] * S[0]
                p[1] = p[1] * S[1]
                p[2] = p[2] * S[2]
                
        for n in self.N:
            n.p[0] = n.p[0] * S[0]
            n.p[1] = n.p[1] * S[1]
            n.p[2] = n.p[2] * S[2]
        
    
    #calculate the adjacency weighting function for the network given a set of vectors X = (x, y, z) and weight exponent k
    def adjacencyweight(self, P, k=200, length_threshold = 25, dtype=np.float32):
        V = self.vectors(dtype)                                                 #get the vectors representing each segment
        #V = V[0:n_vectors, :]
        L = np.expand_dims(np.sqrt((V**2).sum(1)), 1)                           #calculate the length of each vector
        
        outliers = L > length_threshold                                         #remove outliers based on the length_threshold
        V = np.delete(V, np.where(outliers), 0)
        L = np.delete(L, np.where(outliers))
        V = V/L[:,None]                                                         #normalize the vectors
        
        P = np.stack(spharmonics.sph2cart(1, P[0], P[1]), P[0].ndim)        
        PV = P[...,None,:] * V
        cos_alpha = PV.sum(PV.ndim-1)
        W = np.abs(cos_alpha) ** k

        return W, L