segmentation / graph_gui

Browse Code »

Commit 9f9f17883ce16db5a1cfa46f0dfa21853c0f4f60

Authored by Pavel Govyadinov 2019-08-05 16:13:56 -0500

0 parents

clead up version of the GUI submitted to vis 2019

Showing 9 changed files with 5125 additions and 0 deletions Show diff stats

Inline Side-by-side

GraphCanvas.py 0 → 100644

Wrap text Show/Hide comments View file @9f9f178

	1	+++ a/GraphCanvas.py
	1	+#!/usr/bin/env python3
	2	+# -- coding: utf-8 --
	3	+"""
	4	+Created on Mon Aug 5 15:43:59 2019
	5	+
	6	+@author: pavel
	7	+ The GraphCanvas class that extends the scene class in vispy in order to draw
	8	+ the graph object. This class is wrapped in a QTcanvas.
	9	+"""
	10	+
	11	+from vispy import gloo, scene
	12	+from vispy.gloo import set_viewport, set_state, clear, set_blend_color
	13	+from vispy.util.transforms import perspective, translate, rotate, scale
	14	+import vispy.gloo.gl as glcore
	15	+
	16	+import numpy as np
	17	+import math
	18	+import network_dep as nwt
	19	+
	20	+#import the graph shaders.
	21	+from graph_shaders import vert, frag, vs, fs
	22	+from subgraph_shaders import vert_s, frag_s, vs_s, fs_s
	23	+
	24	+#The graph canvas class that
	25	+class GraphCanvas(scene.SceneCanvas):
	26	+
	27	+ """
	28	+ Initialization method.
	29	+ Generates the 512x512 canvas, makes it available for drawing the fills all
	30	+ the GLSL shaders with dummy data.
	31	+ """
	32	+ def __init__(self, **kwargs):
	33	+ # Initialize the canvas for real
	34	+ scene.SceneCanvas.__init__(self, size=(512, 512), **kwargs)
	35	+ #Unfreeze the canvas to make dynamic interaction possible
	36	+ self.unfreeze()
	37	+
	38	+ #initialize all the boolean and dictionary variables.
	39	+ ps = self.pixel_scale
	40	+ self.subgraphs = False
	41	+ self.position = 50, 50
	42	+ self.down=False;
	43	+
	44	+ #Dictionaries to store the unique color ID, cluster ID, and edge-to-ID
	45	+ #dictionaries.
	46	+ self.color_dict = {}
	47	+ self.cluster_dict = {}
	48	+ self.edge_dict = {}
	49	+
	50	+ #Booleans the storage for the current "Path", i.e. edges the user selected.
	51	+ self.pathing = False
	52	+ self.path = []
	53	+ self.full_path = []
	54	+
	55	+ #utility variables used for storing the cluster being moved and all the
	56	+ #nodes and edges that belong to that cluster and move along with it.
	57	+ self.moving = False
	58	+ self.moving_cluster = False
	59	+ self.selection = False
	60	+
	61	+ n = 10
	62	+ ne = 10
	63	+ #Init dummy structures
	64	+ self.uniforms = [('u_graph_size', np.float32, 3)]
	65	+ self.data = np.zeros(n, dtype=[('a_position', np.float32, 3),
	66	+ ('a_fg_color', np.float32, 4),
	67	+ ('a_bg_color', np.float32, 4),
	68	+ ('a_size', np.float32, 1),
	69	+ ('a_linewidth', np.float32, 1),
	70	+ ('a_unique_id', np.float32, 4),
	71	+ ])
	72	+
	73	+ self.clusters = np.zeros(n, dtype=[('a_position', np.float32, 3),
	74	+ ('a_bg_color', np.float32, 4),
	75	+ ('a_value', np.float32, 2),
	76	+ ('a_unique_id', np.float32, 4),
	77	+ ])
	78	+
	79	+ self.edges = np.random.randint(size=(ne, 2), low=0,
	80	+ high=n-1).astype(np.uint32)
	81	+ self.edges_s = np.random.randint(size=(ne, 4), low=0,
	82	+ high=n-1).astype(np.uint32)
	83	+ self.data['a_position'] = np.hstack((0.25 * np.random.randn(n, 2),
	84	+ np.zeros((n, 1))))
	85	+ self.data['a_fg_color'] = 0, 0, 0, 1.0
	86	+ color = np.random.uniform(0.5, 1., (n, 3))
	87	+ self.data['a_bg_color'] = np.hstack((color, np.zeros((n, 1))))
	88	+ self.data['a_size'] = np.random.randint(size=n, low=8ps, high=20ps)
	89	+ self.data['a_linewidth'] = 8.*ps
	90	+ self.data['a_unique_id'] = np.hstack((color, np.ones((n, 1))))
	91	+ #self.uniforms['u_graph_size'] = [1.0, 1.0, 1.0]
	92	+ self.translate = [0,0,0]
	93	+ self.scale = [1,1,1]
	94	+ #color = np.random.uniform(0.5, 1., (ne, 3))
	95	+ #self.linecolor = np.hstack((color, np.ones((ne, 1))))
	96	+ #color = np.random.uniform(0.5, 1., (ne, 3))
	97	+ #self.linecolor = np.hstack((color, np.ones((ne, 1))))
	98	+ self.u_antialias = 1
	99	+
	100	+ #init dummy vertex and index buffers.
	101	+ self.vbo = gloo.VertexBuffer(self.data)
	102	+ self.vbo_s = gloo.VertexBuffer(self.clusters)
	103	+
	104	+ #Need to initialize thick lines.
	105	+ self.index = gloo.IndexBuffer(self.edges)
	106	+ self.index_s = gloo.IndexBuffer(self.edges_s)
	107	+
	108	+ #Set the view matrices.
	109	+ self.view = np.eye(4, dtype=np.float32)
	110	+ self.model = np.eye(4, dtype=np.float32)
	111	+ self.projection = np.eye(4, dtype=np.float32)
	112	+
	113	+ #init shaders used for vertices of the full graph.
	114	+ self.program = gloo.Program(vert, frag)
	115	+ self.program.bind(self.vbo)
	116	+ self.program['u_size'] = 1
	117	+ self.program['u_antialias'] = self.u_antialias
	118	+ self.program['u_model'] = self.model
	119	+ self.program['u_view'] = self.view
	120	+ self.program['u_projection'] = self.projection
	121	+ self.program['u_graph_size'] = [1.0, 1.0, 1.0]
	122	+ self.program['u_picking'] = False
	123	+
	124	+ #init shades used for the edges in the graph
	125	+ self.program_e = gloo.Program(vs, fs)
	126	+ self.program_e['u_size'] = 1
	127	+ self.program_e['u_model'] = self.model
	128	+ self.program_e['u_view'] = self.view
	129	+ self.program_e['u_projection'] = self.projection
	130	+ #self.program_e['l_color'] = self.linecolor.astype(np.float32)
	131	+ self.program_e.bind(self.vbo)
	132	+
	133	+ #init shaders used to the vertices in the subgraph graph.
	134	+ self.program_s = gloo.Program(vert_s, frag_s)
	135	+ self.program_s.bind(self.vbo_s)
	136	+ self.program_s['u_model'] = self.model
	137	+ self.program_s['u_view'] = self.view
	138	+ self.program_s['u_projection'] = self.projection
	139	+ self.program_s['u_graph_size'] = [1.0, 1.0, 1.0]
	140	+ self.program_s['u_picking'] = False
	141	+
	142	+ #init shaders used for the subgraph-edges
	143	+ self.program_e_s = gloo.Program(vs_s, fs_s)
	144	+ self.program_e_s['u_size'] = 1
	145	+ self.program_e_s['u_model'] = self.model
	146	+ self.program_e_s['u_view'] = self.view
	147	+ self.program_e_s['u_projection'] = self.projection
	148	+ #self.program_e['l_color'] = self.linecolor.astype(np.float32)
	149	+ self.program_e_s.bind(self.vbo_s)
	150	+
	151	+
	152	+ #set up the viewport and the gl state.
	153	+ set_viewport(0, 0, *self.physical_size)
	154	+
	155	+ set_state(clear_color='white', depth_test=True, blend=True,
	156	+ blend_func=('src_alpha', 'one_minus_src_alpha'), depth_func = ('lequal'))
	157	+
	158	+ """
	159	+ Function that recolors vertices based on the selected statistic
	160	+ Maps the statisic stored in G to a colormap passed to the function
	161	+ Then updates the necessary color array.
	162	+ """
	163	+ def color_vertices(self, G, vertex_property, dtype = False, cm = 'plasma'):
	164	+ #if we are visualing the clusters we should use a discrete colormap
	165	+ #otherwise use the passed colormap
	166	+ if dtype == True:
	167	+ G.vertex_properties["RGBA"] = nwt.Network.map_property_to_color(G, G.vertex_properties["clusters"])
	168	+ else:
	169	+ G.vertex_properties["RGBA"] = nwt.Network.map_property_to_color(G, G.vertex_properties[vertex_property], colormap=cm)
	170	+ #set the color and update the Vertices.
	171	+ self.current_color = vertex_property
	172	+ color = G.vertex_properties["RGBA"].get_2d_array(range(4)).T
	173	+ self.data['a_bg_color'] = color
	174	+ self.vbo = gloo.VertexBuffer(self.data)
	175	+ self.program.bind(self.vbo)
	176	+ #self.program_e.bind(self.vbo)
	177	+ self.update()
	178	+
	179	+ """
	180	+ Maps a statistic of the vertices based on the size of the canvas to size of
	181	+ the drawn object.
	182	+ """
	183	+ def size_vertices(self, G, propertymap):
	184	+ size = nwt.Network.map_vertices_to_range(G, [30self.pixel_scale, 8self.pixel_scale], propertymap).get_array()
	185	+ self.data['a_size'] = size
	186	+ self.vbo = gloo.VertexBuffer(self.data)
	187	+ self.program.bind(self.vbo)
	188	+ #self.program_e.bind(self.vbo)
	189	+ self.update()
	190	+
	191	+
	192	+ """
	193	+ Function to dim all nodes and edges that do not belong to a cluster chosen
	194	+ in the graph view. Returns a copy of the graph with the alpha channel saved.
	195	+ OPTMIZE HERE: could just return an alpha array to reduce memory usage.
	196	+ """
	197	+ def focus_on_cluster(self, G, c_id):
	198	+ G_copy = nwt.gt.Graph(G, directed=False)
	199	+ e_color = G_copy.edge_properties["RGBA"].get_2d_array(range(4)).T
	200	+ vertices = np.argwhere(self.labels != c_id)
	201	+ for v in range(vertices.shape[0]):
	202	+ idx = vertices[v][0]
	203	+ vtx = G_copy.vertex(idx)
	204	+ for e in vtx.all_edges():
	205	+ if (int(e.source()), int(e.target())) in self.edge_dict.keys():
	206	+ index = int(self.edge_dict[int(e.source()), int(e.target())])
	207	+ if vtx == int(e.source()):
	208	+ e_color[index][3] = 0.05
	209	+ elif vtx == int(e.target()):
	210	+ e_color[index][3] = 0.05
	211	+ else:
	212	+ index = int(self.edge_dict[int(e.target()), int(e.source())])
	213	+ if vtx == int(e.target()):
	214	+ e_color[index][3] = 0.05
	215	+ elif vtx == int(e.source()):
	216	+ e_color[index][3] = 0.05
	217	+
	218	+ G_copy.edge_properties["RGBA"] = G_copy.new_edge_property("vector<double>", vals = e_color)
	219	+
	220	+ return G_copy
	221	+
	222	+ """
	223	+ Function that sets the size of the vertices based on the distance from the
	224	+ camera.
	225	+ """
	226	+ def vertexSizeFromDistance(self, G, camera_pos):
	227	+ location = G.vertex_properties["p"].get_2d_array(range(3)).T
	228	+ cam_array = np.zeros(location.shape, dtype=np.float32)
	229	+ len_array = np.zeros(location.shape[0])
	230	+ offset_array = np.zeros(location.shape, dtype=np.float32)
	231	+ cam_array[:][0:3] = camera_pos
	232	+ offset = [(bbu[0]-bbl[0])/2, (bbu[1]-bbl[1])/2, (bbu[2]-bbl[2])/2]
	233	+ location = location - offset
	234	+ location = location - camera_pos
	235	+ for i in range(location.shape[0]):
	236	+ len_array[i] = np.sqrt(np.power(location[i][0],2) + np.power(location[i][1],2) + np.power(location[i][2],2))
	237	+
	238	+ G.vertex_properties['dist_from_camera'] = G.new_vertex_property('float', vals=len_array)
	239	+ self.data['a_size'] = nwt.Network.map_vertices_to_range(G, [1self.pixel_scale, 60self.pixel_scale], 'dist_from_camera').get_array()
	240	+
	241	+
	242	+ size = nwt.Network.map_vertices_to_range(G, [1.0, 0.5], 'dist_from_camera').get_array()
	243	+ edges = G.get_edges()
	244	+ for e in range(edges.shape[0]):
	245	+ idx = int(4*edges[e][2])
	246	+ self.line_data['a_linewidth'][idx] = size[edges[e][0]]
	247	+ self.line_data['a_linewidth'][idx+1] = size[edges[e][1]]
	248	+ self.line_data['a_linewidth'][idx+2] = size[edges[e][0]]
	249	+ self.line_data['a_linewidth'][idx+3] = size[edges[e][1]]
	250	+
	251	+
	252	+ #self.vbo = gloo.VertexBuffer(self.data)
	253	+ #self.vbo_line = gloo.VertexBuffer(self.line_data)
	254	+ #self.program.bind(self.vbo)
	255	+ #self.program_e.bind(self.vbo_line)
	256	+ #self.update()
	257	+
	258	+ """
	259	+ Function that scales the alpha channel of each vertex in the graph based on
	260	+ The distance from the camera.
	261	+ Sometimes needs to be done separetly.
	262	+ """
	263	+ def vertexAlphaFromDistance(self, G, camera_pos):
	264	+ location = G.vertex_properties["p"].get_2d_array(range(3)).T
	265	+ cam_array = np.zeros(location.shape, dtype=np.float32)
	266	+ len_array = np.zeros(location.shape[0])
	267	+ #offset_array = np.zeros(location.shape, dtype=np.float32)
	268	+ cam_array[:][0:3] = camera_pos
	269	+ offset = [(bbu[0]-bbl[0])/2, (bbu[1]-bbl[1])/2, (bbu[2]-bbl[2])/2]
	270	+ location = location - offset
	271	+ location = location - camera_pos
	272	+ for i in range(location.shape[0]):
	273	+ len_array[i] = np.sqrt(np.power(location[i][0],2) + np.power(location[i][1],2) + np.power(location[i][2],2))
	274	+
	275	+
	276	+ test = nwt.Network.map_vertices_to_range(G, [0.0, 1.0], 'dist_from_camera').get_array()
	277	+ color = G.vertex_properties["RGBA"].get_2d_array(range(4)).T
	278	+ for i in range(location.shape[0]):
	279	+ color[i][3] = test[i]
	280	+ G.vertex_properties["RGBA"] = G.new_vertex_property("vector<double>", vals = color)
	281	+ self.data['a_bg_color'] = color
	282	+
	283	+ edges = G.get_edges()
	284	+ for e in range(edges.shape[0]):
	285	+ idx = int(4*edges[e][2])
	286	+ self.line_data['a_fg_color'][idx] = color[edges[e][0]]
	287	+ self.line_data['a_fg_color'][idx+1] = color[edges[e][1]]
	288	+ self.line_data['a_fg_color'][idx+2] = color[edges[e][0]]
	289	+ self.line_data['a_fg_color'][idx+3] = color[edges[e][1]]
	290	+
	291	+
	292	+
	293	+ self.vbo = gloo.VertexBuffer(self.data)
	294	+ self.vbo_line = gloo.VertexBuffer(self.line_data)
	295	+ self.program.bind(self.vbo)
	296	+ self.program_e.bind(self.vbo_line)
	297	+ self.update()
	298	+
	299	+ """
	300	+ Sets the edge color based on the the cluster the vertices belongs to
	301	+ Propertymap is a VERTEXPROPERTYMAP since the color of the edges is based
	302	+ on the clusters the edges belong to.
	303	+ """
	304	+ def color_edges(self, G, propertymap="clusters"):
	305	+ if propertymap == "clusters":
	306	+ for e in G.edges():
	307	+ if G.vertex_properties[propertymap][e.source()] == G.vertex_properties[propertymap][e.target()]:
	308	+ G.edge_properties["RGBA"][e] = G.vertex_properties["RGBA"][e.source()]
	309	+ else:
	310	+ G.edge_properties["RGBA"][e] = [0.0, 0.0, 0.0, 0.8]
	311	+
	312	+
	313	+ """
	314	+ Helper function that generates the framebuffer object that stores the vertices
	315	+ Generates the vertex buffer based on the graph G that is passed to the function
	316	+ Sets the color, generates the graph and subgraph color if necessary.
	317	+ """
	318	+ def gen_vertex_vbo(self, G):
	319	+ color = G.vertex_properties["RGBA"].get_2d_array(range(4)).T
	320	+ size = nwt.Network.map_vertices_to_range(G, [30self.pixel_scale, 8self.pixel_scale], 'degree').get_array()
	321	+
	322	+ position = G.vertex_properties["pos"].get_2d_array(range(3)).T
	323	+ #for p in range(position.shape[0]):
	324	+ # position[p][0] = position[p][0] + self.clusters["a_position"][G.vertex_properties["clusters"][G.vertex(p)]][0]
	325	+ # position[p][1] = position[p][1] + self.clusters["a_position"][G.vertex_properties["clusters"][G.vertex(p)]][1]
	326	+ # position[p][2] = position[p][2] + self.clusters["a_position"][G.vertex_properties["clusters"][G.vertex(p)]][2]
	327	+ #G.vertex_properties["pos"] = G.new_vertex_property("vector<double>", vals = position)
	328	+ edges = G.get_edges();
	329	+ edges = edges[:, 0:2]
	330	+ #width = nwt.Network.map_edges_to_range(G, [1self.pixel_scale, 5self.pixel_scale], 'volume').get_array()
	331	+ #ecolor = G.edge_properties["RGBA"].get_2d_array(range(4)).T
	332	+
	333	+ self.data = np.zeros(G.num_vertices(), dtype=[('a_position', np.float32, 3),
	334	+ ('a_fg_color', np.float32, 4),
	335	+ ('a_bg_color', np.float32, 4),
	336	+ ('a_size', np.float32, 1),
	337	+ ('a_linewidth', np.float32, 1),
	338	+ ('a_unique_id', np.float32, 4),
	339	+ ('a_selection', np.float32, 1),
	340	+ ])
	341	+
	342	+ #self.edges = edges.astype(np.uint32)
	343	+ self.data['a_position'] = position
	344	+ #fg color is the color of the ring
	345	+ self.data['a_fg_color'] = 0, 0, 0, 1
	346	+ self.data['a_bg_color'] = color
	347	+ self.data['a_size'] = size
	348	+ self.data['a_linewidth'] = 4.*self.pixel_scale
	349	+ self.data['a_unique_id'] = self.gen_vertex_id(G)
	350	+ self.data['a_selection'] = G.vertex_properties["selection"].get_array()
	351	+ #self.data['a_graph_size'] = [bbu-bbl]
	352	+
	353	+ self.program['u_graph_size'] = [bbu-bbl]
	354	+
	355	+ self.vbo = gloo.VertexBuffer(self.data)
	356	+ self.gen_line_vbo(G)
	357	+ if(self.subgraphs):
	358	+ self.vbo_s = gloo.VertexBuffer(self.clusters)
	359	+ self.index_s = gloo.IndexBuffer(self.edges_s)
	360	+ #self.index = gloo.IndexBuffer(self.edges)
	361	+ self.program_e.bind(self.vbo_line)
	362	+ self.program.bind(self.vbo)
	363	+ if(self.subgraphs):
	364	+ #self.program_e_s.bind(self.vbo_s)
	365	+ self.program_s.bind(self.vbo_s)
	366	+ print(self.view)
	367	+ self.update()
	368	+
	369	+ """
	370	+ Helper function that creates colored "block" lines based on the edges
	371	+ in the graph. Generates the framebuffer object and fills it with the relavant data.
	372	+ Note that each line segment is saved as a two triangles that share the same
	373	+ two points on the centerline, but are offset according to the normal of the
	374	+ line segmente to control thickness dynamically.
	375	+ """
	376	+ def gen_line_vbo(self, G):
	377	+ #Set the data.
	378	+ self.line_data = np.zeros(G.num_edges()*4, dtype=[('a_position', np.float32, 3),
	379	+ ('a_normal', np.float32, 2),
	380	+ ('a_fg_color', np.float32, 4),
	381	+ ('a_linewidth', np.float32, 1),
	382	+ ])
	383	+ self.edges = np.random.randint(size=(G.num_edges()*2, 3), low=0,
	384	+ high=(G.num_edges()-1)).astype(np.uint32)
	385	+ color = G.edge_properties["RGBA"].get_2d_array(range(4)).T
	386	+ edges = G.get_edges()
	387	+ #size need to be changed to the size based on the current property map
	388	+ size = nwt.Network.map_vertices_to_range(G, [1.0, 0.5], 'degree').get_array()
	389	+ for e in range(edges.shape[0]):
	390	+ idx = int(4*edges[e][2])
	391	+ p0 = G.vertex_properties["pos"][G.vertex(edges[e][0])]
	392	+ p1 = G.vertex_properties["pos"][G.vertex(edges[e][1])]
	393	+ d = np.subtract(p1, p0)
	394	+ #d_norm = np.multiply(d, 1/np.sqrt(np.power(d[0],2) + np.power(d[1],2)))
	395	+ d_norm = d[0:2]
	396	+ d_norm = d_norm / np.sqrt(np.power(d_norm[0],2) + np.power(d_norm[1],2))
	397	+ norm = np.zeros((2,), dtype=np.float32)
	398	+ norm[0] = d_norm[1]
	399	+ norm[1] = d_norm[0]*-1
	400	+ #print(np.sqrt(norm[0]norm[0] + norm[1]norm[1]))
	401	+ #thickness = G.edge_properties["thickness"][e]
	402	+ thickness = 1.0
	403	+ self.edge_dict[int(edges[e][0]), int(edges[e][1])] = int(edges[e][2])
	404	+ self.line_data['a_position'][idx] = p0
	405	+ self.line_data['a_normal'][idx] = norm
	406	+ self.line_data['a_fg_color'][idx] = color[edges[e][2]]
	407	+ #a_linewidth is a vector.
	408	+ self.line_data['a_linewidth'][idx] = size[edges[e][0]]
	409	+
	410	+ self.line_data['a_position'][idx+1] = p1
	411	+ self.line_data['a_normal'][idx+1] = norm
	412	+ self.line_data['a_fg_color'][idx+1] = color[edges[e][2]]
	413	+ self.line_data['a_linewidth'][idx+1] = size[edges[e][1]]
	414	+
	415	+ self.line_data['a_position'][idx+2] = p0
	416	+ self.line_data['a_normal'][idx+2] = -norm
	417	+ self.line_data['a_fg_color'][idx+2] = color[edges[e][2]]
	418	+ self.line_data['a_linewidth'][idx+2] = size[edges[e][0]]
	419	+
	420	+ self.line_data['a_position'][idx+3] = p1
	421	+ self.line_data['a_normal'][idx+3] = -norm
	422	+ self.line_data['a_fg_color'][idx+3] = color[edges[e][2]]
	423	+ self.line_data['a_linewidth'][idx+3] = size[edges[e][1]]
	424	+
	425	+ self.edges[e*2] = [idx, idx+1, idx+3]
	426	+ self.edges[e*2+1] = [idx, idx+2, idx+3]
	427	+
	428	+ #Set the buffer object and update the shader programs.
	429	+ self.program_e = gloo.Program(vs, fs)
	430	+ #self.program_e['l_color'] = self.linecolor.astype(np.float32)
	431	+ self.vbo_line = gloo.VertexBuffer(self.line_data)
	432	+ self.index = gloo.IndexBuffer(self.edges)
	433	+ self.program_e['u_size'] = 1
	434	+ self.program_e['u_model'] = self.model
	435	+ self.program_e['u_view'] = self.view
	436	+ self.program_e['u_projection'] = self.projection
	437	+ self.program_e.bind(self.vbo_line)
	438	+
	439	+
	440	+ """
	441	+ Helper function that generates the edges between the cluster in the layout.
	442	+ Color is based on the cluster source/target color and transitions between the
	443	+ two.
	444	+ """
	445	+ def gen_cluster_line_vbo(self, G):
	446	+ #create a graph that stores the edges of between the clusters
	447	+ self.G_cluster = nwt.gt.Graph(directed=False)
	448	+ self.G_cluster.vertex_properties["pos"] = self.G_cluster.new_vertex_property("vector<double>", val=np.zeros((3,1), dtype=np.float32))
	449	+ self.G_cluster.vertex_properties["RGBA"] = self.G_cluster.new_vertex_property("vector<double>", val=np.zeros((4,1), dtype=np.float32))
	450	+ for v in range(len(self.cluster_pos)):
	451	+ self.G_cluster.add_vertex()
	452	+ self.G_cluster.vertex_properties["pos"][self.G_cluster.vertex(v)] = np.asarray(self.cluster_pos[v], dtype=np.float32)
	453	+ self.G_cluster.edge_properties["weight"] = self.G_cluster.new_edge_property("int", val = 0)
	454	+ #for each edge in the original graph, generate appropriate subgraph edges without repretiions
	455	+ #i.e. controls the thichness of the edges in the subgraph view.
	456	+ for e in G.edges():
	457	+ #if the source and target cluster is not equal to each other
	458	+ #add an inter subgraph edge.
	459	+ if(G.vertex_properties["clusters"][e.source()] != G.vertex_properties["clusters"][e.target()]):
	460	+ #temp_e.append([G.vertex_properties["clusters"][e.source()], G.vertex_properties["clusters"][e.target()]])
	461	+ self.G_cluster.add_edge(self.G_cluster.vertex(G.vertex_properties["clusters"][e.source()]), \
	462	+ self.G_cluster.vertex(G.vertex_properties["clusters"][e.target()]))
	463	+ self.G_cluster.edge_properties["weight"][self.G_cluster.edge(self.G_cluster.vertex(G.vertex_properties["clusters"][e.source()]), \
	464	+ self.G_cluster.vertex(G.vertex_properties["clusters"][e.target()]))] += 1
	465	+ self.G_cluster.vertex_properties["RGBA"][self.G_cluster.vertex(G.vertex_properties["clusters"][e.source()])] \
	466	+ = G.vertex_properties["RGBA"][e.source()]
	467	+ self.G_cluster.vertex_properties["RGBA"][self.G_cluster.vertex(G.vertex_properties["clusters"][e.target()])] \
	468	+ = G.vertex_properties["RGBA"][e.target()]
	469	+
	470	+ self.cluster_line_data = np.zeros(self.G_cluster.num_edges()*4, dtype=[('a_position', np.float32, 3),
	471	+ ('a_normal', np.float32, 2),
	472	+ ('a_fg_color', np.float32, 4),
	473	+ ('a_linewidth', np.float32, 1),
	474	+ ])
	475	+ self.cluster_edges = np.random.randint(size=(self.G_cluster.num_edges()*2, 3), low=0,
	476	+ high=(G.num_edges()-1)).astype(np.uint32)
	477	+
	478	+ edges = self.G_cluster.get_edges()
	479	+ #size need to be changed to the size based on the current property map
	480	+ size = nwt.Network.map_edges_to_range(self.G_cluster, [1.0, 0.5], 'weight').get_array()
	481	+ color = self.G_cluster.vertex_properties["RGBA"].get_2d_array(range(4)).T
	482	+ #generate the vertex buffer and the connections buffer.
	483	+ for e in range(edges.shape[0]):
	484	+ idx = int(4*edges[e][2])
	485	+ p0 = self.G_cluster.vertex_properties["pos"][self.G_cluster.vertex(edges[e][0])]
	486	+ p1 = self.G_cluster.vertex_properties["pos"][self.G_cluster.vertex(edges[e][1])]
	487	+ d = np.subtract(p1, p0)
	488	+ #d_norm = np.multiply(d, 1/np.sqrt(np.power(d[0],2) + np.power(d[1],2)))
	489	+ d_norm = d[0:2]
	490	+ d_norm = d_norm / np.sqrt(np.power(d_norm[0],2) + np.power(d_norm[1],2))
	491	+ norm = np.zeros((2,), dtype=np.float32)
	492	+ norm[0] = d_norm[1]
	493	+ norm[1] = d_norm[0]*-1
	494	+ #print(np.sqrt(norm[0]norm[0] + norm[1]norm[1]))
	495	+ #thickness = G.edge_properties["thickness"][e]
	496	+ self.cluster_dict[int(edges[e][0]), int(edges[e][1])] = int(edges[e][2])
	497	+ self.cluster_line_data['a_position'][idx] = p0
	498	+ self.cluster_line_data['a_normal'][idx] = norm
	499	+ self.cluster_line_data['a_fg_color'][idx] = color[edges[e][0]]
	500	+ self.cluster_line_data['a_linewidth'][idx] = size[e]
	501	+
	502	+ self.cluster_line_data['a_position'][idx+1] = p1
	503	+ self.cluster_line_data['a_normal'][idx+1] = norm
	504	+ self.cluster_line_data['a_fg_color'][idx+1] = color[edges[e][1]]
	505	+ self.cluster_line_data['a_linewidth'][idx+1] = size[e]
	506	+
	507	+ self.cluster_line_data['a_position'][idx+2] = p0
	508	+ self.cluster_line_data['a_normal'][idx+2] = -norm
	509	+ self.cluster_line_data['a_fg_color'][idx+2] = color[edges[e][0]]
	510	+ self.cluster_line_data['a_linewidth'][idx+2] = size[e]
	511	+
	512	+ self.cluster_line_data['a_position'][idx+3] = p1
	513	+ self.cluster_line_data['a_normal'][idx+3] = -norm
	514	+ self.cluster_line_data['a_fg_color'][idx+3] = color[edges[e][1]]
	515	+ self.cluster_line_data['a_linewidth'][idx+3] = size[e]
	516	+
	517	+ self.cluster_edges[e*2] = [idx, idx+1, idx+3]
	518	+ self.cluster_edges[e*2+1] = [idx, idx+2, idx+3]
	519	+
	520	+
	521	+ self.program_e_s = gloo.Program(vs_s, fs_s)
	522	+ self.index_clusters_s = gloo.IndexBuffer(self.cluster_edges)
	523	+ self.vbo_cluster_lines = gloo.VertexBuffer(self.cluster_line_data)
	524	+
	525	+ self.program_e_s['u_size'] = 1
	526	+ self.program_e_s['u_model'] = self.model
	527	+ self.program_e_s['u_view'] = self.view
	528	+ self.program_e_s['u_projection'] = self.projection
	529	+ self.program_e_s.bind(self.vbo_cluster_lines)
	530	+
	531	+
	532	+ """
	533	+ Updates the vertex buffers based on the current position of the cluster.
	534	+ Updates it's position.
	535	+ """
	536	+ def update_cluster_line_vbo(self):
	537	+
	538	+ for v in range(len(self.cluster_pos)):
	539	+ self.G_cluster.vertex_properties["pos"][self.G_cluster.vertex(v)] = np.asarray(self.cluster_pos[v], dtype=np.float32)
	540	+ #OPTIMIZE HERE to update only one cluster at a time.
	541	+ edges = self.G_cluster.get_edges()
	542	+ #size need to be changed to the size based on the current property map
	543	+ size = nwt.Network.map_edges_to_range(self.G_cluster, [1.0, 0.5], 'weight').get_array()
	544	+ color = self.G_cluster.vertex_properties["RGBA"].get_2d_array(range(4)).T
	545	+ for e in range(edges.shape[0]):
	546	+ idx = int(4*edges[e][2])
	547	+ p0 = self.G_cluster.vertex_properties["pos"][self.G_cluster.vertex(edges[e][0])]
	548	+ p1 = self.G_cluster.vertex_properties["pos"][self.G_cluster.vertex(edges[e][1])]
	549	+ d = np.subtract(p1, p0)
	550	+ #d_norm = np.multiply(d, 1/np.sqrt(np.power(d[0],2) + np.power(d[1],2)))
	551	+ d_norm = d[0:2]
	552	+ d_norm = d_norm / np.sqrt(np.power(d_norm[0],2) + np.power(d_norm[1],2))
	553	+ norm = np.zeros((2,), dtype=np.float32)
	554	+ norm[0] = d_norm[1]
	555	+ norm[1] = d_norm[0]*-1
	556	+ #print(np.sqrt(norm[0]norm[0] + norm[1]norm[1]))
	557	+ #thickness = G.edge_properties["thickness"][e]
	558	+ self.cluster_dict[int(edges[e][0]), int(edges[e][1])] = int(edges[e][2])
	559	+ self.cluster_line_data['a_position'][idx] = p0
	560	+ self.cluster_line_data['a_normal'][idx] = norm
	561	+ self.cluster_line_data['a_fg_color'][idx] = color[edges[e][0]]
	562	+ self.cluster_line_data['a_linewidth'][idx] = size[e]
	563	+
	564	+ self.cluster_line_data['a_position'][idx+1] = p1
	565	+ self.cluster_line_data['a_normal'][idx+1] = norm
	566	+ self.cluster_line_data['a_fg_color'][idx+1] = color[edges[e][1]]
	567	+ self.cluster_line_data['a_linewidth'][idx+1] = size[e]
	568	+
	569	+ self.cluster_line_data['a_position'][idx+2] = p0
	570	+ self.cluster_line_data['a_normal'][idx+2] = -norm
	571	+ self.cluster_line_data['a_fg_color'][idx+2] = color[edges[e][0]]
	572	+ self.cluster_line_data['a_linewidth'][idx+2] = size[e]
	573	+
	574	+ self.cluster_line_data['a_position'][idx+3] = p1
	575	+ self.cluster_line_data['a_normal'][idx+3] = -norm
	576	+ self.cluster_line_data['a_fg_color'][idx+3] = color[edges[e][1]]
	577	+ self.cluster_line_data['a_linewidth'][idx+3] = size[e]
	578	+
	579	+ self.program_e_s = gloo.Program(vs_s, fs_s)
	580	+ self.index_clusters_s = gloo.IndexBuffer(self.cluster_edges)
	581	+ self.vbo_cluster_lines = gloo.VertexBuffer(self.cluster_line_data)
	582	+
	583	+ self.program_e_s['u_size'] = 1
	584	+ self.program_e_s['u_model'] = self.model
	585	+ self.program_e_s['u_view'] = self.view
	586	+ self.program_e_s['u_projection'] = self.projection
	587	+ self.program_e_s.bind(self.vbo_cluster_lines)
	588	+
	589	+ """
	590	+ Genererates a unique index for every vertex.
	591	+ """
	592	+ def gen_vertex_id(self, G):
	593	+ self.color_dict = {}
	594	+ base = [0, 0, 0, 255]
	595	+ idx = 0
	596	+ #colors = cm.get_cmap('Wistia', G.num_vertices()*2)
	597	+ v_id = np.zeros((G.num_vertices(), 4), dtype=np.float32)
	598	+ for v in G.vertices():
	599	+ color = np.multiply(base, 1/255.0)
	600	+ v_id[int(v)] = color
	601	+ self.color_dict[tuple(color)] = int(v)
	602	+ idx += 1
	603	+ base = [int(idx/(255*255)), int((idx/255)%255), int(idx%255), 255]
	604	+
	605	+ return(v_id)
	606	+
	607	+ """
	608	+ Generates a unique index for every cluster.
	609	+ """
	610	+ def gen_cluster_id(self, G):
	611	+ self.cluster_dict = {}
	612	+ base = [0, 0, 0, 255]
	613	+ idx = 0
	614	+ #colors = cm.get_cmap('Wistia', G.num_vertices()*2)
	615	+ v_id = np.zeros((self.n_c, 4), dtype=np.float32)
	616	+ for v in range(self.n_c):
	617	+ color = np.multiply(base, 1/255.0)
	618	+ v_id[int(v)] = color
	619	+ self.cluster_dict[tuple(color)] = int(v)
	620	+ idx += 1
	621	+ base = [int(idx/(255*255)), int((idx/255)%255), int(idx%255), 255]
	622	+
	623	+ return(v_id)
	624	+
	625	+
	626	+ """
	627	+ Generates the bounding box of the radial glyph.
	628	+ """
	629	+ def gen_cluster_coords(self, center, diameter):
	630	+ radius = diameter/2.0
	631	+ top = center[1]+radius
	632	+ bottom = center[1]-radius
	633	+ left = center[0]-radius
	634	+ right = center[0]+radius
	635	+
	636	+ positions = [[right, bottom, center[2]],
	637	+ [right, top, center[2]],
	638	+ [left, top, center[2]],
	639	+ [left, bottom, center[2]]]
	640	+
	641	+
	642	+ values = [[1.0, -1.0],
	643	+ [1.0, 1.0,],
	644	+ [-1.0, 1.0],
	645	+ [-1.0, -1.0]]
	646	+ return positions, values
	647	+
	648	+
	649	+ """
	650	+ Layout algorithm that expands the cluster based on the location of the of the clusters
	651	+ """
	652	+ def expand_based_on_clusters(self, G, n):
	653	+ pos = G.vertex_properties["pos"].get_2d_array(range(3)).T
	654	+ for p in range(pos.shape[0]):
	655	+ pos[p][0] = pos[p][0] + self.clusters["a_position"][G.vertex_properties["clusters"][G.vertex(p)]][0]
	656	+ pos[p][1] = pos[p][1] + self.clusters["a_position"][G.vertex_properties["clusters"][G.vertex(p)]][1]
	657	+ pos[p][2] = pos[p][2] + self.clusters["a_position"][G.vertex_properties["clusters"][G.vertex(p)]][2]
	658	+ G.vertex_properties["pos"] = G.new_vertex_property("vector<double>", vals = pos)
	659	+ for i in range(n):
	660	+ index = 4*i
	661	+ #generate the vertex filter for this cluster
	662	+ num_v_in_cluster = len(np.argwhere(self.labels == i))
	663	+ vfilt = np.zeros([G.num_vertices(), 1], dtype="bool")
	664	+ vfilt[np.argwhere(self.labels == i)] = 1
	665	+ vfilt_prop = G.new_vertex_property("bool", vals = vfilt)
	666	+ G.set_vertex_filter(vfilt_prop)
	667	+
	668	+ #get the filtered properties
	669	+ g = nwt.gt.Graph(G, prune=True, directed=False)
	670	+ positions = g.vertex_properties["pos"].get_2d_array(range(3)).T
	671	+ position = np.sum(positions, 0)/num_v_in_cluster
	672	+ p, v = self.gen_cluster_coords(position, np.sum(g.vertex_properties['degree'].get_array()))
	673	+ self.clusters['a_position'][index:index+4] = np.asarray(p, dtype=np.float32)
	674	+ self.clusters['a_value'][index:index+4] = np.asarray(v, dtype=np.float32)
	675	+ G.clear_filters()
	676	+ self.cluster_pos[i] = position
	677	+ color = G.vertex_properties["RGBA"].get_2d_array(range(4)).T
	678	+ size = nwt.Network.map_vertices_to_range(G, [30self.pixel_scale, 8self.pixel_scale], 'degree').get_array()
	679	+
	680	+ position = G.vertex_properties["pos"].get_2d_array(range(3)).T
	681	+ #for p in range(position.shape[0]):
	682	+ # position[p][0] = position[p][0] + self.clusters["a_position"][G.vertex_properties["clusters"][G.vertex(p)]][0]
	683	+ # position[p][1] = position[p][1] + self.clusters["a_position"][G.vertex_properties["clusters"][G.vertex(p)]][1]
	684	+ # position[p][2] = position[p][2] + self.clusters["a_position"][G.vertex_properties["clusters"][G.vertex(p)]][2]
	685	+ #G.vertex_properties["pos"] = G.new_vertex_property("vector<double>", vals = position)
	686	+ edges = G.get_edges();
	687	+ edges = edges[:, 0:2]
	688	+ #width = nwt.Network.map_edges_to_range(G, [1self.pixel_scale, 5self.pixel_scale], 'volume').get_array()
	689	+ #ecolor = G.edge_properties["RGBA"].get_2d_array(range(4)).T
	690	+
	691	+ self.data = np.zeros(G.num_vertices(), dtype=[('a_position', np.float32, 3),
	692	+ ('a_fg_color', np.float32, 4),
	693	+ ('a_bg_color', np.float32, 4),
	694	+ ('a_size', np.float32, 1),
	695	+ ('a_linewidth', np.float32, 1),
	696	+ ('a_unique_id', np.float32, 4),
	697	+ ('a_selection', np.float32, 1),
	698	+ ])
	699	+
	700	+ #self.edges = edges.astype(np.uint32)
	701	+ self.data['a_position'] = position
	702	+ #fg color is the color of the ring
	703	+ self.data['a_fg_color'] = 0, 0, 0, 1
	704	+ self.data['a_bg_color'] = color
	705	+ self.data['a_size'] = size
	706	+ self.data['a_linewidth'] = 4.*self.pixel_scale
	707	+ self.data['a_unique_id'] = self.gen_vertex_id(G)
	708	+ self.data['a_selection'] = G.vertex_properties["selection"].get_array()
	709	+ #self.data['a_graph_size'] = [bbu-bbl]
	710	+
	711	+ self.program['u_graph_size'] = [bbu-bbl]
	712	+ self.vbo = gloo.VertexBuffer(self.data)
	713	+ self.gen_line_vbo(G)
	714	+ if(self.subgraphs):
	715	+ self.vbo_s = gloo.VertexBuffer(self.clusters)
	716	+ self.index_s = gloo.IndexBuffer(self.edges_s)
	717	+ #self.index = gloo.IndexBuffer(self.edges)
	718	+ self.program_e.bind(self.vbo_line)
	719	+ self.program.bind(self.vbo)
	720	+ if(self.subgraphs):
	721	+ #self.program_e_s.bind(self.vbo_s)
	722	+ self.program_s.bind(self.vbo_s)
	723	+ print(self.view)
	724	+ self.update()
	725	+
	726	+
	727	+ """
	728	+ Function that generates the clusters for an unclustered graph
	729	+ These are to be represented by the arcs
	730	+ """
	731	+ def gen_clusters(self, G, bbl, bbu, n_c = None, edge_metric = 'volume', vertex_metric = 'degree'):
	732	+
	733	+ #Generate the clusters
	734	+ self.labels = nwt.Network.spectral_clustering(G,'length', n_clusters = n_c)
	735	+ #self.labels = nwt.Network.spectral_clustering(G,'length')
	736	+
	737	+ #Add clusters as a vertex property
	738	+ G.vertex_properties["clusters"] = G.new_vertex_property("int", vals=self.labels)
	739	+ num_clusters = len(np.unique(self.labels))
	740	+ self.n_c = n_c
	741	+
	742	+ #add colormap
	743	+ G.vertex_properties["RGBA"] = nwt.Network.map_property_to_color(G, G.vertex_properties["clusters"])
	744	+
	745	+ #generate an empty property set for the clusters.
	746	+ self.clusters = np.zeros(num_clusters*4, dtype=[('a_position', np.float32, 3),
	747	+ ('a_value', np.float32, 2),
	748	+ ('a_bg_color', np.float32, 4),
	749	+ ('a_cluster_color', np.float32, 4),
	750	+ ('a_arc_length', np.float32, 1),
	751	+ ('a_outer_arc_length', np.float32, 4),
	752	+ ('a_unique_id', np.float32, 4),
	753	+ ])
	754	+ self.edges_s = np.random.randint(size=(num_clusters*2, 3), low=0,
	755	+ high=4).astype(np.uint32)
	756	+ #fill the foreground color as halo
	757	+ #self.clusters['a_fg_color'] = 1., 1., 1., 0.0
	758	+ #self.clusters['a_linewidth'] = 4.*self.pixel_scale
	759	+
	760	+ G.vertex_properties["pos"] = nwt.gt.sfdp_layout(G, groups = G.vertex_properties["clusters"], pos = G.vertex_properties["pos"])
	761	+ temp = [];
	762	+ temp_pos = [];
	763	+ #Find the global total of the metric.
	764	+ global_metric = np.sum(G.edge_properties[edge_metric].get_array(), 0)
	765	+ unique_color = self.gen_cluster_id(G)
	766	+
	767	+ #generate the property values for every cluster
	768	+ for i in range(num_clusters):
	769	+ idx = 4*i
	770	+ #generate the vertex filter for this cluster
	771	+ num_v_in_cluster = len(np.argwhere(self.labels == i))
	772	+ vfilt = np.zeros([G.num_vertices(), 1], dtype="bool")
	773	+ vfilt[np.argwhere(self.labels == i)] = 1
	774	+ vfilt_prop = G.new_vertex_property("bool", vals = vfilt)
	775	+ G.set_vertex_filter(vfilt_prop)
	776	+
	777	+ #get the filtered properties
	778	+ g = nwt.gt.Graph(G, prune=True, directed=False)
	779	+ positions = g.vertex_properties["pos"].get_2d_array(range(3)).T
	780	+ position = np.sum(positions, 0)/num_v_in_cluster
	781	+
	782	+ #calculate the arclength for the global statistic
	783	+ arc_length = np.sum(g.edge_properties[edge_metric].get_array(), 0)/global_metricnp.pi2
	784	+ arc_length_vertex = np.ones((4,1), dtype = np.float32)
	785	+ array = g.vertex_properties[vertex_metric].get_array()
	786	+
	787	+ #calculate metric distribution and turn it into arc_lengths
	788	+ t_vertex_metric = np.sum(array)
	789	+ arc_length_vertex[0] = np.sum(array < 2)/t_vertex_metric
	790	+ arc_length_vertex[1] = np.sum(array == 2)/t_vertex_metric
	791	+ arc_length_vertex[2] = np.sum(array == 3)/t_vertex_metric
	792	+ arc_length_vertex[3] = np.sum(array > 3)/t_vertex_metric
	793	+
	794	+ #arc_length_vertex = np.asarray(arc_length_vertex, dtype = np.float32)
	795	+ #arc_length_vertex = (max(arc_length_vertex) - min(arc_length_vertex)) \
	796	+ #* (arc_length_vertex- min(arc_length_vertex))
	797	+ for j in range(len(arc_length_vertex)):
	798	+ if j != 0:
	799	+ arc_length_vertex[j] += arc_length_vertex[j-1]
	800	+ print("arc_length before ", arc_length_vertex, " and sum to ", sum(arc_length_vertex))
	801	+ arc_length_vertex = np.asarray(arc_length_vertex, dtype = np.float32)
	802	+ arc_length_vertex = (np.pi - -np.pi)/(max(arc_length_vertex) - min(arc_length_vertex)) \
	803	+ * (arc_length_vertex- min(arc_length_vertex)) + (-np.pi)
	804	+ print(arc_length_vertex)
	805	+ #print(arc_length)
	806	+
	807	+
	808	+ temp_pos.append(position)
	809	+
	810	+ #generate the color for every vertex,
	811	+ #since all vertices belong to the same cluster we can check only
	812	+ #one vertex for the cluster color.
	813	+
	814	+ self.clusters['a_cluster_color'][idx:idx+4] = g.vertex_properties["RGBA"][g.vertex(0)]
	815	+ self.clusters['a_bg_color'][idx:idx+4] = [0.1, 0.1, 0.1, 1.0]
	816	+ self.clusters['a_unique_id'][idx:idx+4] = unique_color[i]
	817	+
	818	+ #The arc-length representing one global metric.
	819	+ self.clusters['a_arc_length'][idx:idx+4] = arc_length
	820	+ self.clusters['a_outer_arc_length'][idx:idx+4] = arc_length_vertex[:].T
	821	+
	822	+ temp.append(np.sum(g.vertex_properties['degree'].get_array()))
	823	+ G.clear_filters()
	824	+ print(self.clusters['a_outer_arc_length'])
	825	+ maximum = max(temp)
	826	+ minimum = min(temp)
	827	+ if len(temp) > 1:
	828	+ temp = ((temp-minimum)/(maximum-minimum)(60self.pixel_scale)+20*self.pixel_scale)
	829	+ else:
	830	+ temp = [60*self.pixel_scale]
	831	+ for i in range(num_clusters):
	832	+ index = i*4
	833	+ index_t = i*2
	834	+ p, v = self.gen_cluster_coords(temp_pos[i], temp[i]*2.0)
	835	+ self.clusters['a_position'][index:index+4] = np.asarray(p, dtype=np.float32)
	836	+ self.clusters['a_value'][index:index+4] = np.asarray(v, dtype=np.float32)
	837	+
	838	+ self.edges_s[index_t] = [index, index+1, index+2]
	839	+ self.edges_s[index_t+1] = [index, index+2, index+3]
	840	+ #self.edges_s[i][0:4] = np.asarray(range(index, index+4), dtype=np.uint32)
	841	+ #self.edges_s[i]
	842	+ self.cluster_pos = temp_pos
	843	+ #self.expand_based_on_clusters(G, self.n_c)
	844	+ G.clear_filters()
	845	+# self.edges_s[1][0:4] = np.asarray(range(0, 0+4), dtype=np.uint32)
	846	+# self.edges_s[1][4] = 0
	847	+# self.edges_s[1][5] = 0+2
	848	+#
	849	+# self.edges_s[0][0:4] = np.asarray(range(index, index+4), dtype=np.uint32)
	850	+# self.edges_s[0][4] = index
	851	+# self.edges_s[0][5] = index+2
	852	+ #self.clusters['a_size'] = temp
	853	+ self.gen_cluster_line_vbo(G)
	854	+ self.program_s['u_graph_size'] = [bbu-bbl]
	855	+ #if len(temp_e) > 0:
	856	+ # self.edges_s = np.unique(np.asarray(temp_e, np.uint32), axis=0)
	857	+ #else:
	858	+ # self.edges_s = []
	859	+ #print(self.edges_s)
	860	+
	861	+
	862	+ """
	863	+ Function that expands that generates the layout and updates the buffer
	864	+
	865	+ """
	866	+ def expand_clusters(self, G, n_c):
	867	+ self.expand_based_on_clusters(G, n_c)
	868	+ self.gen_cluster_line_vbo(G)
	869	+ if(self.subgraphs):
	870	+ self.vbo_s = gloo.VertexBuffer(self.clusters)
	871	+ self.index_s = gloo.IndexBuffer(self.edges_s)
	872	+ self.program_e.bind(self.vbo_line)
	873	+ self.program.bind(self.vbo)
	874	+ if(self.subgraphs):
	875	+ self.program_s.bind(self.vbo_s)
	876	+ print(self.view)
	877	+ self.update()
	878	+
	879	+ """
	880	+ Loads the data G and generates all the buffers necessary as well as performs
	881	+ spectral clustering on the graph passed if the subgraph is set to true.
	882	+ """
	883	+ def set_data(self, G, bbl, bbu, subgraph=True):
	884	+ print("Setting data")
	885	+ self.G = G
	886	+ self.bbl = bbl
	887	+ self.bbu = bbu
	888	+ clear(color=True, depth=True)
	889	+ self.subgraphs = True
	890	+ self.current_color = "clusters"
	891	+ if(subgraph==True):
	892	+ self.gen_clusters(G, bbl, bbu, n_c=19)
	893	+
	894	+ #color based on clusters
	895	+ self.color_edges(G)
	896	+
	897	+ color = G.vertex_properties["RGBA"].get_2d_array(range(4)).T
	898	+ size = nwt.Network.map_vertices_to_range(G, [30self.pixel_scale, 8self.pixel_scale], 'degree').get_array()
	899	+
	900	+ position = G.vertex_properties["pos"].get_2d_array(range(3)).T
	901	+ #for p in range(position.shape[0]):
	902	+ # position[p][0] = position[p][0] + self.clusters["a_position"][G.vertex_properties["clusters"][G.vertex(p)]][0]
	903	+ # position[p][1] = position[p][1] + self.clusters["a_position"][G.vertex_properties["clusters"][G.vertex(p)]][1]
	904	+ # position[p][2] = position[p][2] + self.clusters["a_position"][G.vertex_properties["clusters"][G.vertex(p)]][2]
	905	+ #G.vertex_properties["pos"] = G.new_vertex_property("vector<double>", vals = position)
	906	+ edges = G.get_edges();
	907	+ edges = edges[:, 0:2]
	908	+ #width = nwt.Network.map_edges_to_range(G, [1self.pixel_scale, 5self.pixel_scale], 'volume').get_array()
	909	+ #ecolor = G.edge_properties["RGBA"].get_2d_array(range(4)).T
	910	+
	911	+ self.data = np.zeros(G.num_vertices(), dtype=[('a_position', np.float32, 3),
	912	+ ('a_fg_color', np.float32, 4),
	913	+ ('a_bg_color', np.float32, 4),
	914	+ ('a_size', np.float32, 1),
	915	+ ('a_linewidth', np.float32, 1),
	916	+ ('a_unique_id', np.float32, 4),
	917	+ ('a_selection', np.float32, 1),
	918	+ ])
	919	+
	920	+ #self.edges = edges.astype(np.uint32)
	921	+ self.data['a_position'] = position
	922	+ #fg color is the color of the ring
	923	+ self.data['a_fg_color'] = 0, 0, 0, 1
	924	+ self.data['a_bg_color'] = color
	925	+ self.data['a_size'] = size
	926	+ self.data['a_linewidth'] = 4.*self.pixel_scale
	927	+ self.data['a_unique_id'] = self.gen_vertex_id(G)
	928	+ self.data['a_selection'] = G.vertex_properties["selection"].get_array()
	929	+ #self.data['a_graph_size'] = [bbu-bbl]
	930	+
	931	+ self.program['u_graph_size'] = [bbu-bbl]
	932	+
	933	+ self.vbo = gloo.VertexBuffer(self.data)
	934	+ self.gen_line_vbo(G)
	935	+ #self.gen_cylinder_vbo(G)
	936	+ if(self.subgraphs):
	937	+ self.vbo_s = gloo.VertexBuffer(self.clusters)
	938	+ self.index_s = gloo.IndexBuffer(self.edges_s)
	939	+ #self.index = gloo.IndexBuffer(self.edges)
	940	+ self.program_e.bind(self.vbo_line)
	941	+ self.program.bind(self.vbo)
	942	+ if(self.subgraphs):
	943	+ #self.program_e_s.bind(self.vbo_s)
	944	+ self.program_s.bind(self.vbo_s)
	945	+ print(self.view)
	946	+ self.update()
	947	+
	948	+ """
	949	+ Function that changes and redraws the buffer during a resize event.
	950	+ """
	951	+ def on_resize(self, event):
	952	+ set_viewport(0, 0, *event.physical_size)
	953	+ self.fbo = gloo.FrameBuffer(color=gloo.RenderBuffer(self.size[::-1]), depth=gloo.RenderBuffer(self.size[::-1]))
	954	+
	955	+
	956	+ """
	957	+ Overloaded function that is called during every self.update() call
	958	+ """
	959	+ def on_draw(self, event):
	960	+ clear(color='white', depth=True)
	961	+ self.program_e.draw('triangles', indices=self.index)
	962	+ self.program.draw('points')
	963	+ #self.program_e.draw('lines')
	964	+ if(self.subgraphs):
	965	+ self.program_e_s.draw('triangles', indices=self.index_clusters_s)
	966	+ self.program_s.draw('triangles', indices=self.index_s)
	967	+
	968	+ """
	969	+ Function performed during a mouse click (either right or left)
	970	+ gets the unique id of the object drawn underneath the cursor
	971	+ handles the cases depending on whether the click happened to a cluster
	972	+ or a vertex. Edges are not interactable (yet)
	973	+ """
	974	+ def get_clicked_id(self, event, clusters = False):
	975	+ #Get the framebuffer coordinates of the click
	976	+ coord = self.transforms.get_transform('canvas', 'framebuffer').map(event.pos)
	977	+
	978	+ #get the framebuffer where each element is rendered as a unique color
	979	+ size = self.size;
	980	+ self.fbo = gloo.FrameBuffer(color=gloo.RenderBuffer(size[::-1]), depth=gloo.RenderBuffer(size[::-1]))
	981	+ buff = gloo.read_pixels((0,0,self.physical_size[0], self.physical_size[1]))
	982	+ #imsave("test_ori.png", buff)
	983	+ self.fbo.activate()
	984	+ if clusters == False:
	985	+ self.program['u_picking'] = True
	986	+ clear(color='white', depth=True)
	987	+ self.program.draw('points')
	988	+ buff = gloo.read_pixels((0,0,self.physical_size[0], self.physical_size[1]))
	989	+ #imsave("test.png", buff)
	990	+
	991	+ #return to the original state
	992	+ self.fbo.deactivate()
	993	+ self.program['u_picking'] = False
	994	+ else:
	995	+ self.program_s['u_picking'] = True
	996	+ clear(color='white', depth=True)
	997	+ self.program_s.draw('triangles', indices=self.index_s)
	998	+ buff = gloo.read_pixels((0,0,self.physical_size[0], self.physical_size[1]))
	999	+ #imsave("test.png", buff)
	1000	+
	1001	+ #return to the original state
	1002	+ self.fbo.deactivate()
	1003	+ self.program_s['u_picking'] = False
	1004	+
	1005	+ #print(buff[self.physical_size[1]-int(coord[1]), int(coord[0])])
	1006	+
	1007	+ #Get the color under the click.
	1008	+ #Keep in mind that the buff is y, x
	1009	+ #And 0,0 is in the top RIGHT corner.
	1010	+ #IGNORE THE DOCUMENTATION
	1011	+ color = np.multiply(buff[self.physical_size[1]-int(coord[1]), int(coord[0])], 1/255.0)
	1012	+ #if (tuple(color) not in self.color_dict):
	1013	+ # print("clicked on nothing")
	1014	+ #else:
	1015	+ # print(self.color_dict[tuple(color)])
	1016	+
	1017	+ #reset the original buffer
	1018	+ self.update()
	1019	+
	1020	+ #Return the element under the click.
	1021	+ if clusters == False:
	1022	+ if(tuple(color) not in self.color_dict):
	1023	+ return None
	1024	+ else:
	1025	+ return self.color_dict[tuple(color)]
	1026	+ else:
	1027	+ if(tuple(color) not in self.cluster_dict):
	1028	+ return None
	1029	+ else:
	1030	+ return self.cluster_dict[tuple(color)]
	1031	+
	1032	+
	1033	+ """
	1034	+ Top level handle-mouse presee event for either left or right click
	1035	+ """
	1036	+ def on_mouse_press(self, event):
	1037	+
	1038	+ def update_view():
	1039	+ self.location = event.pos
	1040	+ self.program['u_view'] = self.view
	1041	+ self.program_e['u_view'] = self.view
	1042	+ self.program_s['u_view'] = self.view
	1043	+ self.program_e_s['u_view'] = self.view
	1044	+ self.down = True
	1045	+
	1046	+
	1047	+# if(event.button == 2):
	1048	+## menu = QtWidgets.QMenu(self.parent)
	1049	+## NS = menu.addAction('Node Size')
	1050	+## NC = menu.addAction('Node Color')
	1051	+## action = menu.exec_(self.parent.globalPos())
	1052	+## if action == NS:
	1053	+# print("right_click")
	1054	+# #if menu.exec_(event.globalPos()):
	1055	+# # print(item.text())
	1056	+ if(event.button == 1):
	1057	+ if(self.view[0][0] > 0.0010):
	1058	+ c_id = self.get_clicked_id(event)
	1059	+ if(c_id != None):
	1060	+ self.original_point = G.vertex_properties["pos"][G.vertex(c_id)]
	1061	+ self.location = event.pos
	1062	+ self.moving = True
	1063	+ self.down = True
	1064	+ self.c_id = [c_id]
	1065	+ else:
	1066	+ update_view()
	1067	+ #print("Clicked on:", event.pos)
	1068	+ else:
	1069	+ c_id = self.get_clicked_id(event, True)
	1070	+ print(c_id)
	1071	+ if(c_id != None):
	1072	+ self.original_point = self.cluster_pos[c_id]
	1073	+ self.location = event.pos
	1074	+ self.moving = True
	1075	+ self.down = True
	1076	+ self.c_id = [c_id]
	1077	+ self.moving_cluster = True
	1078	+ else:
	1079	+ update_view()
	1080	+
	1081	+
	1082	+ """
	1083	+ Handles the double click event that it responsible for path selection.
	1084	+ Generates paths our of consecutive paths out of the selected vertices.
	1085	+ """
	1086	+ def on_mouse_double_click(self, event):
	1087	+ def update_vbo(self):
	1088	+ self.vbo = gloo.VertexBuffer(self.data)
	1089	+ self.program.bind(self.vbo)
	1090	+ self.update()
	1091	+
	1092	+ def add_to_path(self, source, target):
	1093	+ vl, el = nwt.gt.graph_tool.topology.shortest_path(G, G.vertex(source), G.vertex(target), weights=G.edge_properties["av_radius"])
	1094	+ for v in vl:
	1095	+ if(int(v) not in self.path):
	1096	+ G.vertex_properties["selection"][v] = 2.0
	1097	+ self.data['a_selection'][int(v)] = 2.0
	1098	+ if(int(v) not in self.full_path):
	1099	+ self.full_path.append(int(v))
	1100	+
	1101	+ if (event.button == 1):
	1102	+ if(self.view[0][0] > 0.0010):
	1103	+ c_id = self.get_clicked_id(event)
	1104	+ if(c_id != None):
	1105	+ #check whether this is the first node to be selected
	1106	+ if(self.pathing == False):
	1107	+ #if it is, select that node and turn the pathing variable on.
	1108	+ G.vertex_properties["selection"][G.vertex(c_id)] = 1.0
	1109	+ self.pathing = True
	1110	+ if(c_id not in self.path):
	1111	+ self.path.append(c_id)
	1112	+ self.data['a_selection'][c_id] = 1.0
	1113	+ update_vbo(self)
	1114	+ print("I turned on the first node")
	1115	+ else:
	1116	+ if(G.vertex_properties["selection"][G.vertex(c_id)] == 1.0):
	1117	+ G.vertex_properties["selection"][G.vertex(c_id)] = 0.0
	1118	+ self.path.remove(c_id)
	1119	+ self.data['a_selection'][c_id] = 0.0
	1120	+ update_vbo(self)
	1121	+ print("I turned off a node")
	1122	+ elif(G.vertex_properties["selection"][G.vertex(c_id)] == 0.0):
	1123	+ G.vertex_properties["selection"][G.vertex(c_id)] = 1.0
	1124	+ if(c_id not in self.path):
	1125	+ self.path.append(c_id)
	1126	+ self.data['a_selection'][c_id] = 1.0
	1127	+ update_vbo(self)
	1128	+ print("I turned on a node")
	1129	+ if(len(self.path) >= 2):
	1130	+ for i in range(len(self.path)-1):
	1131	+ add_to_path(self, self.path[i], self.path[i+1])
	1132	+ update_vbo(self)
	1133	+ #THIS IS WHERE I LEFT IT OFF.
	1134	+ if(np.sum(G.vertex_properties["selection"].get_array()) == 0):
	1135	+ self.pathing = False
	1136	+
	1137	+
	1138	+
	1139	+# elif(np.sum(self.G.vertex_properties["selection"].get_array()) :
	1140	+# self.G.vertex_properties["selection"][self.G.vertex(c_id)] == False
	1141	+ print("clicked on: ", c_id, " ", self.path)
	1142	+
	1143	+ """
	1144	+ Resets the variables that are used during the pressdown and move events
	1145	+ """
	1146	+ def on_mouse_release(self, event):
	1147	+ self.down = False
	1148	+ self.moving = False
	1149	+ self.moving_cluster = False
	1150	+ self.c_id = []
	1151	+ #self.location = event.pos
	1152	+ #print("Clicked off:", event.pos)
	1153	+
	1154	+ """
	1155	+ used during the drag evern to update the position of the clusters
	1156	+ """
	1157	+ def update_cluster_position(self, G, pos, offset, c_id):
	1158	+ v_pos = G.vertex_properties["pos"].get_2d_array(range(3)).T
	1159	+ vertices = np.argwhere(self.labels == c_id)
	1160	+ for v in range(vertices.shape[0]):
	1161	+ idx = vertices[v][0]
	1162	+ v_pos[idx][0] = v_pos[idx][0] + offset[0]
	1163	+ v_pos[idx][1] = v_pos[idx][1] + offset[1]
	1164	+ v_pos[idx][2] = v_pos[idx][2] + offset[2]
	1165	+ self.data['a_position'][idx] = np.asarray([v_pos[idx][0], v_pos[idx][1], v_pos[idx][2]], dtype = np.float32)
	1166	+ #update the edge data by finding all edges connected to the vertex
	1167	+ vtx = self.G.vertex(idx)
	1168	+ for e in vtx.all_edges():
	1169	+ d = np.subtract(G.vertex_properties["pos"][e.source()], G.vertex_properties["pos"][e.target()])
	1170	+ d_norm = d[0:2]
	1171	+ d_norm = d_norm / np.sqrt(np.power(d_norm[0],2) + np.power(d_norm[1],2))
	1172	+ norm = np.zeros((2,), dtype=np.float32)
	1173	+ norm[0] = d_norm[1]
	1174	+ norm[1] = d_norm[0]*-1
	1175	+ if (int(e.source()), int(e.target())) in self.edge_dict.keys():
	1176	+ index = int(self.edge_dict[int(e.source()), int(e.target())])
	1177	+ if vtx == int(e.source()):
	1178	+ self.line_data['a_position'][index*4] = v_pos[idx]
	1179	+ self.line_data['a_position'][index*4+2] = v_pos[idx]
	1180	+ self.line_data['a_normal'][index*4] = norm
	1181	+ self.line_data['a_normal'][index*4+2] = -norm
	1182	+ self.line_data['a_normal'][index*4+1] = norm
	1183	+ self.line_data['a_normal'][index*4+3] = -norm
	1184	+ elif vtx == int(e.target()):
	1185	+ self.line_data['a_position'][index*4+1] = v_pos[idx]
	1186	+ self.line_data['a_position'][index*4+3] = v_pos[idx]
	1187	+ self.line_data['a_normal'][index*4] = norm
	1188	+ self.line_data['a_normal'][index*4+2] = -norm
	1189	+ self.line_data['a_normal'][index*4+1] = norm
	1190	+ self.line_data['a_normal'][index*4+3] = -norm
	1191	+ else:
	1192	+ index = int(self.edge_dict[int(e.target()), int(e.source())])
	1193	+ if vtx == int(e.target()):
	1194	+ self.line_data['a_position'][index*4] = v_pos[idx]
	1195	+ self.line_data['a_position'][index*4+2] = v_pos[idx]
	1196	+ self.line_data['a_normal'][index*4] = norm
	1197	+ self.line_data['a_normal'][index*4+2] = -norm
	1198	+ self.line_data['a_normal'][index*4+1] = norm
	1199	+ self.line_data['a_normal'][index*4+3] = -norm
	1200	+ elif vtx == int(e.source()):
	1201	+ self.line_data['a_position'][index*4+1] = v_pos[idx]
	1202	+ self.line_data['a_position'][index*4+3] = v_pos[idx]
	1203	+ self.line_data['a_normal'][index*4] = norm
	1204	+ self.line_data['a_normal'][index*4+2] = -norm
	1205	+ self.line_data['a_normal'][index*4+1] = norm
	1206	+ self.line_data['a_normal'][index*4+3] = -norm
	1207	+
	1208	+
	1209	+ G.vertex_properties["pos"] = G.new_vertex_property("vector<double>", vals = v_pos)
	1210	+ index = 4*c_id
	1211	+ #generate the vertex filter for this cluster
	1212	+ vfilt = np.zeros([G.num_vertices(), 1], dtype="bool")
	1213	+ vfilt[np.argwhere(self.labels == c_id)] = 1
	1214	+ vfilt_prop = G.new_vertex_property("bool", vals = vfilt)
	1215	+ G.set_vertex_filter(vfilt_prop)
	1216	+
	1217	+ #get the filtered properties
	1218	+ g = nwt.gt.Graph(G, prune=True, directed=False)
	1219	+ p, v = self.gen_cluster_coords(pos, np.sum(g.vertex_properties['degree'].get_array()))
	1220	+ self.clusters['a_position'][index:index+4] = np.asarray(p, dtype=np.float32)
	1221	+ self.clusters['a_value'][index:index+4] = np.asarray(v, dtype=np.float32)
	1222	+ G.clear_filters()
	1223	+ self.cluster_pos[c_id] = pos
	1224	+ self.original_point = pos
	1225	+
	1226	+
	1227	+ """
	1228	+ function that handles the mouse move event in a way that depends on a set
	1229	+ of variables: state of the mouse button, the type of object selected and
	1230	+ the number of objects.
	1231	+ """
	1232	+ def on_mouse_move(self, event):
	1233	+ if(self.down == True):
	1234	+ if(self.moving == True and self.moving_cluster == False):
	1235	+ if(len(self.c_id) < 2):
	1236	+ #Project into GLSpace and get before and after move coordinates
	1237	+ coord = self.transforms.get_transform('canvas', 'render').map(self.location)[:2]
	1238	+ coord2 = self.transforms.get_transform('canvas', 'render').map(event.pos)[:2]
	1239	+ cur_pos = G.vertex_properties["pos"][G.vertex(self.c_id[0])]
	1240	+ #print(cur_pos, " Before")
	1241	+
	1242	+ #Adjust the position of the node based on the current view matrix.
	1243	+ cur_pos[0] = cur_pos[0] - (coord[0]-coord2[0])/self.view[0][0]
	1244	+ cur_pos[1] = cur_pos[1] - (coord[1]-coord2[1])/self.view[0][0]
	1245	+
	1246	+ #print(cur_pos, " After")
	1247	+ #Upload the changed data.
	1248	+ G.vertex_properties["pos"][G.vertex(self.c_id[0])] = cur_pos
	1249	+ self.data['a_position'][self.c_id[0]] = cur_pos
	1250	+
	1251	+ #update the edge data by finding all edges connected to the vertex
	1252	+ v = self.G.vertex(self.c_id[0])
	1253	+ for e in v.all_edges():
	1254	+ d = np.subtract(G.vertex_properties["pos"][e.source()], G.vertex_properties["pos"][e.target()])
	1255	+ d_norm = d[0:2]
	1256	+ d_norm = d_norm / np.sqrt(np.power(d_norm[0],2) + np.power(d_norm[1],2))
	1257	+ norm = np.zeros((2,), dtype=np.float32)
	1258	+ norm[0] = d_norm[1]
	1259	+ norm[1] = d_norm[0]*-1
	1260	+ if (int(e.source()), int(e.target())) in self.edge_dict.keys():
	1261	+ idx = int(self.edge_dict[int(e.source()), int(e.target())])
	1262	+ if self.c_id[0] == int(e.source()):
	1263	+ self.line_data['a_position'][idx*4] = cur_pos
	1264	+ self.line_data['a_position'][idx*4+2] = cur_pos
	1265	+ self.line_data['a_normal'][idx*4] = norm
	1266	+ self.line_data['a_normal'][idx*4+2] = -norm
	1267	+ self.line_data['a_normal'][idx*4+1] = norm
	1268	+ self.line_data['a_normal'][idx*4+3] = -norm
	1269	+ elif self.c_id[0] == int(e.target()):
	1270	+ self.line_data['a_position'][idx*4+1] = cur_pos
	1271	+ self.line_data['a_position'][idx*4+3] = cur_pos
	1272	+ self.line_data['a_normal'][idx*4] = norm
	1273	+ self.line_data['a_normal'][idx*4+2] = -norm
	1274	+ self.line_data['a_normal'][idx*4+1] = norm
	1275	+ self.line_data['a_normal'][idx*4+3] = -norm
	1276	+ else:
	1277	+ idx = int(self.edge_dict[int(e.target()), int(e.source())])
	1278	+ if self.c_id[0] == int(e.target()):
	1279	+ self.line_data['a_position'][idx*4] = cur_pos
	1280	+ self.line_data['a_position'][idx*4+2] = cur_pos
	1281	+ self.line_data['a_normal'][idx*4] = norm
	1282	+ self.line_data['a_normal'][idx*4+2] = -norm
	1283	+ self.line_data['a_normal'][idx*4+1] = norm
	1284	+ self.line_data['a_normal'][idx*4+3] = -norm
	1285	+ elif self.c_id[0] == int(e.source()):
	1286	+ self.line_data['a_position'][idx*4+1] = cur_pos
	1287	+ self.line_data['a_position'][idx*4+3] = cur_pos
	1288	+ self.line_data['a_normal'][idx*4] = norm
	1289	+ self.line_data['a_normal'][idx*4+2] = -norm
	1290	+ self.line_data['a_normal'][idx*4+1] = norm
	1291	+ self.line_data['a_normal'][idx*4+3] = -norm
	1292	+ #self.line_data['a_position'][self.c_id[0]] =
	1293	+ self.vbo = gloo.VertexBuffer(self.data)
	1294	+ self.vbo_line = gloo.VertexBuffer(self.line_data)
	1295	+ #Bind the buffer and redraw.
	1296	+ self.program.bind(self.vbo)
	1297	+ self.program_e.bind(self.vbo_line)
	1298	+ #self.program.draw('points')
	1299	+ self.location = event.pos
	1300	+ self.update()
	1301	+ elif(self.moving == True and self.moving_cluster == True):
	1302	+ if(len(self.c_id) < 2):
	1303	+ #Project into GLSpace and get before and after move coordinates
	1304	+ coord = self.transforms.get_transform('canvas', 'render').map(self.location)[:2]
	1305	+ coord2 = self.transforms.get_transform('canvas', 'render').map(event.pos)[:2]
	1306	+ cur_pos = np.zeros(self.cluster_pos[self.c_id[0]].shape, dtype = np.float32)
	1307	+ offset = np.zeros(self.cluster_pos[self.c_id[0]].shape, dtype = np.float32)
	1308	+ cur_pos[0] = self.cluster_pos[self.c_id[0]][0]
	1309	+ cur_pos[1] = self.cluster_pos[self.c_id[0]][1]
	1310	+ cur_pos[2] = self.cluster_pos[self.c_id[0]][2]
	1311	+ offset[0] = self.cluster_pos[self.c_id[0]][0]
	1312	+ offset[1] = self.cluster_pos[self.c_id[0]][1]
	1313	+ offset[2] = self.cluster_pos[self.c_id[0]][2]
	1314	+# ofset = self.cluster_pos[self.c_id[0]]
	1315	+
	1316	+ #Adjust the position of the node based on the current view matrix.
	1317	+ offset[0] = self.original_point[0] - cur_pos[0] - (coord[0]-coord2[0])/self.view[0][0]
	1318	+ offset[1] = self.original_point[1] - cur_pos[1] - (coord[1]-coord2[1])/self.view[0][0]
	1319	+ cur_pos[0] = cur_pos[0] - (coord[0]-coord2[0])/self.view[0][0]
	1320	+ cur_pos[1] = cur_pos[1] - (coord[1]-coord2[1])/self.view[0][0]
	1321	+
	1322	+ self.update_cluster_position(G, cur_pos, offset, self.c_id[0])
	1323	+ #self.original_point = cur_pos
	1324	+ self.vbo = gloo.VertexBuffer(self.data)
	1325	+ self.vbo_line = gloo.VertexBuffer(self.line_data)
	1326	+ #Bind the buffer and redraw.
	1327	+ self.program.bind(self.vbo)
	1328	+ self.program_e.bind(self.vbo_line)
	1329	+ #self.program.draw('points')
	1330	+ self.location = event.pos
	1331	+ if(self.subgraphs):
	1332	+ self.vbo_s = gloo.VertexBuffer(self.clusters)
	1333	+ self.program_s.bind(self.vbo_s)
	1334	+ self.update_cluster_line_vbo()
	1335	+ self.update()
	1336	+
	1337	+
	1338	+ else:
	1339	+ #print("Mouse at:", event.pos)
	1340	+ #new_model = np.eye(4, dtype=np.float32)
	1341	+ coord = self.transforms.get_transform('canvas', 'render').map(self.location)[:2]
	1342	+ coord2 = self.transforms.get_transform('canvas', 'render').map(event.pos)[:2]
	1343	+ self.translate[0] += (coord[0]-coord2[0])/self.view[0][0]
	1344	+ self.translate[1] += (coord[1]-coord2[1])/self.view[1][1]
	1345	+ #self.view[3][0] = self.view[3][0]-(self.location[0]-event.pos[0])/10000.0
	1346	+ #self.view[3][1] = self.view[3][1]+(self.location[1]-event.pos[1])/10000.0
	1347	+
	1348	+ self.view = np.matmul(translate((self.translate[0], self.translate[1], 0)), scale((self.scale[0], self.scale[1], 0)))
	1349	+
	1350	+ self.program['u_view'] = self.view
	1351	+ self.program_e['u_view'] = self.view
	1352	+ self.program_s['u_view'] = self.view
	1353	+ self.program_e_s['u_view'] = self.view
	1354	+ self.location = event.pos
	1355	+ self.update()
	1356	+
	1357	+ """
	1358	+ Handles the mouse wheel zoom event.
	1359	+ """
	1360	+ def on_mouse_wheel(self, event):
	1361	+
	1362	+ #print(self.view)
	1363	+ #TO_DO IMPLEMENT ZOOM TO CURSOR
	1364	+ #self.view[3][0] = self.view[3][0]-event.pos[0]/10000.0
	1365	+ #self.view[3][1] = self.view[3][1]-event.pos[1]/10000.0
	1366	+ #print(self.scale[0] , self.scale[0]event.delta[1]0.05)
	1367	+ self.scale[0] = self.scale[0] + self.scale[0]event.delta[1]0.05
	1368	+ self.scale[1] = self.scale[1] + self.scale[1]event.delta[1]0.05
	1369	+
	1370	+ self.view = np.matmul(translate((self.translate[0], self.translate[1], 0)),
	1371	+ scale((self.scale[0], self.scale[1], 0)))
	1372	+
	1373	+ #self.view[0][0] = self.view[0][0]+self.view[0][0]event.delta[1]0.05
	1374	+ #self.view[1][1] = self.view[1][1]+self.view[1][1]event.delta[1]0.05
	1375	+ #print(self.view[0][0], " ",self.view[1][1])
	1376	+ #print(self.view)
	1377	+ self.program['u_view'] = self.view
	1378	+ self.program_e['u_view'] = self.view
	1379	+ self.program_s['u_view'] = self.view
	1380	+ self.program_e_s['u_view'] = self.view
	1381	+ #print(event.delta[1])
	1382	+ self.update()
...	...

GraphWidget.py 0 → 100644

Wrap text Show/Hide comments View file @9f9f178

	1	+++ a/GraphWidget.py
	1	+#!/usr/bin/env python3
	2	+# -- coding: utf-8 --
	3	+"""
	4	+Created on Mon Aug 5 15:42:19 2019
	5	+
	6	+@author: pavel
	7	+"""
	8	+
	9	+from GraphCanvas import GraphCanvas
	10	+from pyqtgraph.Qt import QtCore, QtGui, QtWidgets
	11	+import network_dep as nwt
	12	+
	13	+
	14	+"""
	15	+Initializes the entire QTlayout and sets the mouse press events.
	16	+These are connected to the slots such that each is processes by this class
	17	+and the class it wraps.
	18	+"""
	19	+class GraphWidget(QtGui.QWidget):
	20	+ def __init__(self):
	21	+ super(GraphWidget, self).__init__()
	22	+ box = QtGui.QVBoxLayout(self)
	23	+ self.resize(500,500)
	24	+ self.setLayout(box)
	25	+ self.canvas = GraphCanvas()
	26	+ #self.canvas.create_native()
	27	+ box.addWidget(self.canvas.native)
	28	+ self.color = True
	29	+ self.use_3D = False
	30	+ self.camera = [0,0,0]
	31	+
	32	+ self.canvas.events.mouse_press.connect(self.on_mouse_press)
	33	+ self.canvas.events.mouse_release.connect(self.on_mouse_release)
	34	+ self.canvas.events.mouse_move.connect(self.on_mouse_move)
	35	+ self.canvas.events.mouse_double_click.connect(self.on_mouse_double_click)
	36	+
	37	+ #self.show()
	38	+
	39	+ """
	40	+ Handles the right click mouse event at the QWidget level.
	41	+ Depending on where the mouse button in clicked and the current zoom level,
	42	+ the function gets the unique id of the node, cluster or background under the
	43	+ cursor and opens a context menu based on the ID.
	44	+
	45	+ The context menu level actions are also coded in this section of the code.
	46	+ """
	47	+ def on_mouse_press(self, event):
	48	+ if event.button == 2:
	49	+ if self.canvas.view[0][0] >= 0.0010:
	50	+ c_id = self.canvas.get_clicked_id(event)
	51	+ else:
	52	+ c_id = self.canvas.get_clicked_id(event, clusters=True)
	53	+ #right click
	54	+ if(c_id == None):
	55	+ menu = QtWidgets.QMenu(self)
	56	+ NS = menu.addAction('Node Size')
	57	+
	58	+ #tmp = menu.addAction('temp')
	59	+ tmp = menu.addMenu('Node Color')
	60	+ vertex_props = self.canvas.G.vertex_properties.keys()
	61	+ vertex_actions = []
	62	+ for i in range(len(vertex_props)):
	63	+ if(vertex_props[i] != "map" or vertex_props[i] != "RGBA"):
	64	+ vertex_actions.append(tmp.addAction(vertex_props[i]))
	65	+ #tmp.addAction("Cluster")
	66	+ #NC = menu.addAction('Node Color')
	67	+ #EXP = menu.addAction('Export VTK')
	68	+ #EXP_adv = menu.addAction('Export VTK Advanced')
	69	+ NL = menu.addAction('New Layout')
	70	+ action = menu.exec_(event.native.globalPos())
	71	+ print(action)
	72	+ for i in range(len(vertex_actions)):
	73	+ if action == vertex_actions[i]:
	74	+ if vertex_props[i] == "clusters":
	75	+ self.canvas.color_vertices(self.canvas.G, vertex_props[i], dtype=True)
	76	+ else:
	77	+ self.canvas.color_vertices(self.canvas.G, vertex_props[i])
	78	+ print(vertex_props[i])
	79	+ if action == NS:
	80	+ if self.use_3D == False:
	81	+ self.use_3D = True
	82	+ else:
	83	+ self.use_3D = False
	84	+ self.canvas.size_vertices(self.canvas.G, 'degree' )
	85	+ self.canvas.color_vertices(self.canvas.G)
	86	+ #if action == NC:
	87	+ # self.canvas.color_vertices(self.canvas.G, not self.color)
	88	+ # self.color = not self.color
	89	+# if action == EXP:
	90	+# nwt.Network.write_vtk(self.canvas.G, "./nwt.vtk")
	91	+# if action == EXP_adv:
	92	+# nwt.Network.write_vtk(self.canvas.G, "./nwt_median_binned.vtk")
	93	+# nwt.Network.write_vtk(self.canvas.G, "./nwt_median_non_binned.vtk", binning=False)
	94	+# nwt.Network.write_vtk(self.canvas.G, "./nwt_points_wise_binned.vtk", camera=self.camera, binning = True)
	95	+# nwt.Network.write_vtk(self.canvas.G, "./nwt_points_wise_non_binned.vtk", camera=self.camera, binning = False)
	96	+ if action == tmp:
	97	+ print("Stuff")
	98	+ #self.cb = QtWidgets.QComboBox()
	99	+ #vertex_props = self.canvas.G.vertex_properties.keys()
	100	+ #for i in range(len(vertex_props)):
	101	+ # self.cb.addItem(vertex_props[i])
	102	+ #self.cb.currentIndexChanged.connect(self.selection_change)
	103	+ #self.cb.show()
	104	+ if action == NL:
	105	+ self.canvas.G = nwt.Network.gen_new_fd_layout(self.canvas.G)
	106	+ self.canvas.gen_vertex_vbo(self.canvas.G)
	107	+ #self.canvas.set_data(self.canvas.G, self.canvas.bbl, self.canvas.bbu)
	108	+ self.canvas.expand_clusters(self.canvas.G, self.canvas.n_c)
	109	+
	110	+ #self.canvas.size_vertices(self.canvas.G, 'degree_volume')
	111	+ else:
	112	+ if self.canvas.view[0][0] >= 0.0010:
	113	+ menu = QtWidgets.QMenu(self)
	114	+ NS = menu.addAction('Display Node Info')
	115	+ action = menu.exec_(event.native.globalPos())
	116	+ if action == NS:
	117	+ print("Display Node Info")
	118	+ else:
	119	+ menu = QtWidgets.QMenu(self)
	120	+ MnOC = menu.addAction('Minimize Other Clusters')
	121	+ RA = menu.addAction('Reset Alpha')
	122	+ action = menu.exec_(event.native.globalPos())
	123	+ if action == MnOC:
	124	+ new_Graph = self.canvas.focus_on_cluster(self.canvas.G, c_id)
	125	+ self.test.draw_new(new_Graph)
	126	+ if action == RA:
	127	+ self.canvas.reset_alpha(c_id[0])
	128	+
	129	+
	130	+# def selection_change(self, i):
	131	+# self.canvas.size_vertices(self.canvas.G, self.cb.currentText())
	132	+
	133	+
	134	+ """
	135	+ Handles the mouse double click event.
	136	+ Pass-through function.
	137	+ """
	138	+ def on_mouse_double_click(self, event):
	139	+ #print("in applet")
	140	+ n = 0
	141	+
	142	+ """
	143	+ Handles the mouse release event.
	144	+ Pass-through function.
	145	+ """
	146	+ def on_mouse_release(self, event):
	147	+ n = 1
	148	+
	149	+ """
	150	+ Handles the mouse move event.
	151	+ Pass-through function.
	152	+ """
	153	+ def on_mouse_move(self, event):
	154	+ n = 2
	155	+
	156	+ """
	157	+ Internal function that interacts with the tube visualization.
	158	+ """
	159	+ def set_g_view(self, test):
	160	+ self.test = test
	161	+
	162	+ """
	163	+ Function to handle the pyqt Slot that receives the interacted-with signal
	164	+ from the Tube visualization. The signal sends the camera position in 3D coordinates
	165	+ every time the camera is moved.
	166	+ """
	167	+ @QtCore.pyqtSlot(float, float, float)
	168	+ def sigUpdate(self, x, y, z):
	169	+ self.camera = [x,y,z]
	170	+ if(self.use_3D == True):
	171	+ self.camera = [x,y,z]
	172	+ self.canvas.vertexSizeFromDistance(self.canvas.G, [x,y,z])
	173	+ self.canvas.vertexAlphaFromDistance(self.canvas.G, [x,y,z])
	174	+ #print("got signal", x, y, z)
	175	+
	176	+ """
	177	+ Initialization method that connects the slot to the function that handles
	178	+ the information being sent.
	179	+ """
	180	+ def connect(self, signal_object):
	181	+ signal_object.sigUpdate.connect(self.sigUpdate)
	182	+
	183	+
	184	+
...	...

GuiVisPy_tube.py 0 → 100644

Wrap text Show/Hide comments View file @9f9f178

	1	+++ a/GuiVisPy_tube.py
	1	+#!/usr/bin/env python3
	2	+# -- coding: utf-8 --
	3	+"""
	4	+Created on Thu Jan 31 15:29:40 2019
	5	+
	6	+@author: pavel
	7	+`"""
	8	+
	9	+from vispy import app
	10	+
	11	+
	12	+import network_dep as nwt
	13	+
	14	+from pyqtgraph.Qt import QtCore, QtGui, QtWidgets
	15	+import pyqtgraph as pg
	16	+import numpy as np
	17	+import math
	18	+
	19	+from mpl_toolkits.mplot3d import Axes3D
	20	+import matplotlib
	21	+import matplotlib.pyplot as plt
	22	+matplotlib.use('Qt5Agg')
	23	+
	24	+from GraphWidget import GraphWidget
	25	+from TubeWidget import TubeWidget
	26	+
	27	+
	28	+#set the backend. for different versions of PyQt
	29	+app.use_app(backend_name='PyQt5')
	30	+
	31	+#Define a top level application
	32	+appMain = QtGui.QApplication([])
	33	+
	34	+##Define a toplevel Widget
	35	+top = QtGui.QWidget()
	36	+top.resize(900, 900)
	37	+
	38	+
	39	+hist = pg.PlotWidget()
	40	+#fibers = FiberView()
	41	+fibers = TubeWidget()
	42	+fibers.canvas.create_native()
	43	+#fibers = gl.GLViewWidget()
	44	+
	45	+#plt = hist.addPlot()
	46	+
	47	+layout = QtGui.QGridLayout()
	48	+graph = GraphWidget()
	49	+graph.canvas.create_native()
	50	+
	51	+graph.connect(fibers)
	52	+
	53	+#initialize the layout
	54	+layout.addWidget(graph, 0, 0, 2, 3)
	55	+layout.addWidget(fibers, 0, 2, 2, 3)
	56	+layout.setColumnStretch(0, 2)
	57	+layout.setRowStretch(0, 2)
	58	+layout.setColumnStretch(2, 1)
	59	+layout.setRowStretch(0, 1)
	60	+
	61	+top.setLayout(layout)
	62	+top.show()
	63	+
	64	+
	65	+def draw_histogram(G, value):
	66	+ vals = G.edge_properties[value].get_array()
	67	+ y, x = np.histogram(vals,40)
	68	+ hist.plot(x,y, stepMode=True, fillLevel=0, brush=(0,0,255,150))
	69	+
	70	+def load_nwt(filepath):
	71	+ net = nwt.Network(filepath)
	72	+ G = net.createFullGraph_gt()
	73	+ G = net.filterDisconnected(G)
	74	+ #G = net.filterFullGraph_gt(G, erode=True)
	75	+ #G = net.filterFullGraph_gt(G, erode=True)
	76	+ #G = net.gen_new_fd_layout(G)
	77	+ #G = net.filterBorder(G)
	78	+ #nwt.Network.saveGraph_gt(nwt, G, "FilteredNetwork_JACK_3.nwt")
	79	+ color = np.zeros(4, dtype = np.double)
	80	+ color = [0.0, 1.0, 0.0, 1.0]
	81	+ G.edge_properties["RGBA"] = G.new_edge_property("vector<double>", val=color)
	82	+ color = [1.0, 0.0, 0.0, 0.9]
	83	+ G.vertex_properties["RGBA"] = G.new_vertex_property("vector<double>", val=color)
	84	+ bbl, bbu = net.aabb()
	85	+
	86	+
	87	+ return G, bbl, bbu
	88	+
	89	+G, bbl, bbu = load_nwt("/home/pavel/Documents/Python/GraphGuiQt/network_4.nwt")
	90	+#nwt.Network.write_vtk(G)
	91	+
	92	+#G, bbl, bbu = load_nwt("/home/pavel/Documents/Python/GraphGuiQt/network_test_251_1kx3_short_NEW.nwt")
	93	+#ret = nwt.Network.get_affinity_matrix(G, "length")
	94	+node_image = QtGui.QImage("/home/pavel/Documents/Python/GraphGuiQt/node_tex.jpg")
	95	+node_tex = QtGui.QBitmap.fromImage(node_image)
	96	+print(node_tex.depth())
	97	+center = (bbu-bbl)/2.0
	98	+#fibers.opts['distance'] = 5
	99	+# item = NodeItem(G)
	100	+# graph.addItem(item)
	101	+draw_histogram(G, "length")
	102	+graph.canvas.set_data(G, bbl, bbu)
	103	+fibers.canvas.set_data(G, bbl, bbu, 16)
	104	+#fibers.draw_all(G, center, fibers, graph, node_tex)
	105	+graph.set_g_view(fibers)
	106	+## Start Qt event loop unless running in interactive mode.
	107	+if __name__ == '__main__':
	108	+ import sys
	109	+ if (sys.flags.interactive != 1) or not hasattr(QtCore, 'PYQT_VERSION'):
	110	+ QtGui.QApplication.instance().exec_()
	111	+
	112	+
...	...

TubeCanvas.py 0 → 100644

Wrap text Show/Hide comments View file @9f9f178

	1	+++ a/TubeCanvas.py
	1	+#!/usr/bin/env python3
	2	+# -- coding: utf-8 --
	3	+"""
	4	+Created on Mon Aug 5 15:56:47 2019
	5	+
	6	+@author: pavel
	7	+"""
	8	+
	9	+"""
	10	+ Class that extends the vispy SceneCanvas to draw 3D tubes
	11	+"""
	12	+
	13	+from vispy import gloo, scene
	14	+from vispy.gloo import set_viewport, set_state, clear, set_blend_color
	15	+from vispy.util.transforms import perspective, translate, rotate, scale
	16	+import vispy.gloo.gl as glcore
	17	+from vispy.util.quaternion import Quaternion
	18	+
	19	+import numpy as np
	20	+import math
	21	+import network_dep as nwt
	22	+
	23	+
	24	+from tube_shaders import FRAG_SHADER, VERT_SHADER
	25	+
	26	+class TubeDraw(scene.SceneCanvas):
	27	+ #sigUpdate = QtCore.pyqtSignal(float, float, float)
	28	+
	29	+ #Initiates the canvas.
	30	+ def __init__(self, **kwargs):
	31	+ #Initialte the class by calling the superclass
	32	+ scene.SceneCanvas.__init__(self, size=(512,512), keys='interactive', **kwargs)
	33	+ #unfreeze the drawing area to allow for dynamic drawing and interaction
	34	+ self.unfreeze()
	35	+
	36	+ #generate dummy buffers for the meshes
	37	+ self.program = gloo.Program(VERT_SHADER, FRAG_SHADER)
	38	+ self.cylinder_data = np.zeros(5*5, dtype=[('a_position', np.float32, 3),
	39	+ ('a_normal', np.float32, 3),
	40	+ ('a_fg_color', np.float32, 4),
	41	+ #('a_linewidth', np.float32, 1),
	42	+ ])
	43	+ self.triangle_data = np.random.randint(size=(5, 3), low=0,
	44	+ high=(4-1)).astype(np.uint32)
	45	+ self.vbo = gloo.VertexBuffer(self.cylinder_data)
	46	+ self.triangles = gloo.IndexBuffer(self.triangle_data)
	47	+ self.program.bind(self.vbo)
	48	+ self.scale = [1,1,1]
	49	+ self.r1 = np.eye(4, dtype=np.float32)
	50	+ self.r2 = np.eye(4, dtype=np.float32)
	51	+ set_viewport(0,0,*self.physical_size)
	52	+ #set_state(clear_color='white', depth_test=True, blend=True,
	53	+ # blend_func=('src_alpha', 'one_minus_src_alpha'), depth_func = ('less'), cull_face='back')
	54	+ set_state(clear_color='white', depth_test=True, blend=True,
	55	+ blend_func=('src_alpha', 'one_minus_src_alpha'), depth_func = ('lequal'))
	56	+ #set_blend_color(color='black')
	57	+ #set_state('translucent')
	58	+ self.program['u_LightPos'] = [0., 0., -1000.]
	59	+ #self.camera = self.central_widget.add_view()
	60	+ #self.camera.camera = 'turntable'
	61	+ self.down = False
	62	+ self.camera = np.asarray([0.0, 0.0, 200.0], dtype=np.float32)
	63	+ self.up = np.asarray([0., 1., 0.], dtype=np.float32)
	64	+ #self.init_camera = [0.,0.,1000.]
	65	+
	66	+ ##### prototype #####
	67	+ #Set the visualization matrices
	68	+ self.program['u_eye'] = self.camera
	69	+ self.program['u_up'] = self.up
	70	+ self.program['u_target'] = np.asarray([0., 0., 0.], dtype=np.float32)
	71	+
	72	+
	73	+
	74	+ #Load the data necessary to draw all of the microvessels
	75	+ def set_data(self, G, bbu, bbl, num_sides):
	76	+ self.G = G
	77	+ self.num_sides = num_sides
	78	+ self.bbu = bbu
	79	+ self.bbl = bbl
	80	+ bb = nwt.AABB(G).resample_sides(3)
	81	+
	82	+
	83	+ #create program
	84	+ self.gen_cylinder_vbo(self.G, self.num_sides)
	85	+ self.vbo = gloo.VertexBuffer(self.cylinder_data)
	86	+ self.triangles = gloo.IndexBuffer(self.triangle_data)
	87	+
	88	+ #self.view = np.eye(4, dtype=np.float32)
	89	+ self.model = np.eye(4, dtype=np.float32)
	90	+ self.projection = np.eye(4, dtype=np.float32)
	91	+ self.projection = perspective(90.0, self.physical_size[0]/self.physical_size[1], 1.0, 1000.0)
	92	+ #self.projection = perspective(90.0, 1.0, -1.0, 1.0)
	93	+ self.program['u_model'] = self.model
	94	+ self.program['u_LightPos'] = [0., 0., 1000.]
	95	+ #self.program['u_view'] = self.view
	96	+ self.program['u_projection'] = self.projection
	97	+ self.program.bind(self.vbo)
	98	+
	99	+ gloo.set_clear_color('white')
	100	+ self.center = (bbu-bbl)/2.0
	101	+ self.translate = [-self.center[0], -self.center[1], -self.center[2]]
	102	+
	103	+ self.bb = np.ones((26, 3), dtype=np.float32)
	104	+ for i in range(26):
	105	+ for j in range(3):
	106	+ self.bb[i,j] = bb[i][j]
	107	+ self.program['u_bb'] = self.bb
	108	+ print('bb is ', self.bb)
	109	+# for i in range(len(self.translate)):
	110	+# self.camera[i] += self.translate[i]
	111	+
	112	+
	113	+ ##### prototype #####
	114	+ self.camera = self.camera - self.translate
	115	+ self.program['u_eye'] = self.camera
	116	+ self.up = np.cross((np.asarray(self.center, dtype=np.float32)-np.asarray(self.camera, dtype=np.float32)), np.asarray(self.up))
	117	+ self.program['u_up'] = self.up
	118	+ self.program['u_target'] = self.translate
	119	+
	120	+
	121	+
	122	+
	123	+ #self.show()
	124	+
	125	+ #Called during resize of the window in order to redraw the same image in the
	126	+ #larger/smaller area.
	127	+ def on_resize(self, event):
	128	+ width, height = event.physical_size
	129	+ gloo.set_viewport(0, 0, width, height)
	130	+ print(self.physical_size)
	131	+
	132	+ #overloaded function called during the self.update() call to update the current
	133	+ #frame using the GLSL frag/vert shaders
	134	+ def on_draw(self, event):
	135	+ clear(color='white', depth=True)
	136	+ gloo.set_clear_color('white')
	137	+ self.program.draw('triangles', self.triangles)
	138	+ self.projection = perspective(90.0, self.physical_size[0]/self.physical_size[1], 1.0, 1000.0)
	139	+ self.program['u_projection'] = self.projection
	140	+
	141	+ #Creates a cylinder around ever segment in the microvascular network.
	142	+ def gen_cylinder_vbo(self, G, num_sides = 32):
	143	+ i = 0
	144	+ num_pts = 0
	145	+ num_tri = 0
	146	+ for e in G.edges():
	147	+ num_pts += len(self.G.edge_properties["x"][e])
	148	+ num_tri += (len(self.G.edge_properties["x"][e])-1)num_sides2
	149	+ self.cylinder_data = np.zeros(num_pts*num_sides, dtype=[('a_position', np.float32, 3),
	150	+ ('a_normal', np.float32, 3),
	151	+ ('a_fg_color', np.float32, 4),
	152	+ #('a_linewidth', np.float32, 1),
	153	+ ])
	154	+ self.triangle_data = np.random.randint(size=(num_tri, 3), low=0,
	155	+ high=(G.num_edges()-1)).astype(np.uint32)
	156	+ index = 0
	157	+ t_index = 0
	158	+ #for each edge generate a cylinder.
	159	+ for e in G.edges():
	160	+ #print("done")
	161	+ #for each fiber get all the points and the radii
	162	+ X = self.G.edge_properties["x"][e]
	163	+ Y = self.G.edge_properties["y"][e]
	164	+ Z = self.G.edge_properties["z"][e]
	165	+ R = self.G.edge_properties["r"][e]
	166	+ color = G.edge_properties["RGBA"][e]
	167	+ pts = np.array([X,Y,Z]).T
	168	+ circle_pts = np.zeros((pts.shape[0], num_sides, 3), dtype = np.float32)
	169	+ step = 2*np.pi/num_sides
	170	+# U = np.zeros(pts.shape, dtype=np.float32)
	171	+# V = np.zeros(pts.shape, dtype=np.float32)
	172	+# direction = np.zeros(pts.shape, dtype=np.float32)
	173	+
	174	+ #for every point in the edge
	175	+ for p in range(pts.shape[0]):
	176	+ #if first point, generate the circles.
	177	+ #In this case we want to generate a cap if we see the first circle or the last.
	178	+ if(p == 0):
	179	+ #get the direction
	180	+ direction = (pts[p+1] - pts[p])
	181	+ #normalize direction
	182	+ direction = direction/np.sqrt(direction[0]direction[0] + direction[1]direction[1] + direction[2]*direction[2])
	183	+ #generate a vector to use as an element of the cross product
	184	+ Y = np.zeros((3,), dtype = np.float32)
	185	+ Y[0] = 1.
	186	+ #if direction and Y are parallel, change Y
	187	+ if(np.cos(np.dot(Y, direction)) < 0.087):
	188	+ Y[0] = 0.0
	189	+ Y[2] = 1.0
	190	+ #generate first plane vector
	191	+ U = np.cross(direction, Y)
	192	+ U = U/np.sqrt(U[0]U[0] + U[1]U[1] + U[2]*U[2])
	193	+ #generate second plane vector
	194	+ V = np.cross(direction, U)
	195	+ V = V/np.sqrt(V[0]V[0] + V[1]V[1] + V[2]*V[2])
	196	+ #print(R[p],pts[p])
	197	+ #generate circle.
	198	+ for s in range(num_sides):
	199	+ circle_pts[p][s][0] = R[p]np.cos(sstep)V[0]0.5 + R[p]np.sin(sstep)U[0]0.5
	200	+ circle_pts[p][s][1] = R[p]np.cos(sstep)V[1]0.5 + R[p]np.sin(sstep)U[1]0.5
	201	+ circle_pts[p][s][2] = R[p]np.cos(sstep)V[2]0.5 + R[p]np.sin(sstep)U[2]0.5
	202	+ #if last point, copy the previous circle.
	203	+ elif(p == pts.shape[0]-1):
	204	+ for s in range(num_sides):
	205	+ circle_pts[p][s] = circle_pts[p-1][s]
	206	+ for v in range(pts.shape[0]):
	207	+ circle_pts[v,:,0] += pts[v][0]
	208	+ circle_pts[v,:,1] += pts[v][1]
	209	+ circle_pts[v,:,2] += pts[v][2]
	210	+ #otherwise, rotate the circle
	211	+ else:
	212	+ #print(R[p], pts[p])
	213	+ #generate a new direction vector.
	214	+ dir_new = (pts[p+1]-pts[p])
	215	+ dir_new = dir_new/np.sqrt(dir_new[0]dir_new[0] + dir_new[1]dir_new[1] + dir_new[2]*dir_new[2])
	216	+ dir_new2 = (pts[p]-pts[p-1])
	217	+ dir_new2 = dir_new2/np.sqrt(dir_new2[0]dir_new2[0] + dir_new2[1]dir_new2[1] + dir_new2[2]*dir_new2[2])
	218	+ dir_new = dir_new+dir_new2
	219	+ dir_new = dir_new/np.sqrt(dir_new[0]dir_new[0] + dir_new[1]dir_new[1] + dir_new[2]*dir_new[2])
	220	+ #print(dir_new, direction)
	221	+ #generate the quaternion rotation vector for the shortest arc
	222	+ k = 1.0 + np.dot(direction, dir_new)
	223	+ s = 1/np.sqrt(k+k)
	224	+ r = s*np.cross(direction, dir_new)
	225	+ theta = k*s
	226	+ #r = np.cross(direction, dir_new)
	227	+ #r = r/np.sqrt(r[0]r[0] + r[1]r[1] + r[2]*r[2])
	228	+ #calculate the degree of quaternion rotation.
	229	+ #cos_theta = np.sqrt(np.sqrt(np.dot(dir_new, dir_new)) * np.sqrt(np.dot(direction, direction))) + np.dot(dir_new, direction)
	230	+ #cos_theta = np.dot(direction, dir_new)
	231	+ #theta = np.arccos(cos_theta)/2.0
	232	+ #print(r, cos_theta, theta)
	233	+ #quat = np.append(theta, r)
	234	+ q = Quaternion(w=theta, x = r[0], y = r[1], z = r[2]).normalize()
	235	+
	236	+ #print(quat)
	237	+ #q = np.quaternion(quat[0], quat[1], quat[2], quat[3])
	238	+ #rot = Rotation.from_quat(quat, normalized=False)
	239	+ #rot.as_quat()
	240	+ for s in range(num_sides):
	241	+ circle_pts[p][s] = q.rotate_point(circle_pts[p-1][s])
	242	+ #circle_pts[p][s] = rot.apply(circle_pts[p-1][s])
	243	+ #circle_pts[p][s] = q.rotate(circle_pts[p-1][s])
	244	+ #circle_pts[p][s] = np.quaternion.rotate_vectors(q, circle_pts[p][s])
	245	+ #circle_pts[p][s] = q.rotate_vectors(q, circle_pts[p][s])
	246	+ #circle_pts[p][s] = circle_pts[p][s]
	247	+ direction = dir_new
	248	+ #generate the triangles
	249	+ triangles = np.random.randint(size=((pts.shape[0]-1)2(num_sides), 3), low=0,
	250	+ high=(G.num_edges()-1)).astype(np.uint32)
	251	+
	252	+ t_index_temp = 0
	253	+ for ring in range(pts.shape[0]-1):
	254	+ for side in range(num_sides):
	255	+ if(side < num_sides-1):
	256	+ idx_current_point = index+ring*num_sides + side
	257	+ idx_next_ring = index + (ring+1) * num_sides + side
	258	+ triangle1 = [idx_current_point, idx_next_ring, idx_next_ring+1]
	259	+ triangle2 = [idx_next_ring+1, idx_current_point+1, idx_current_point]
	260	+ triangles[t_index_temp] = [idx_current_point, idx_next_ring, idx_next_ring+1]
	261	+ triangles[t_index_temp+1] = [idx_next_ring+1, idx_current_point+1, idx_current_point]
	262	+ self.triangle_data[t_index] = triangle1
	263	+ self.triangle_data[t_index+1] = triangle2
	264	+ t_index += 2
	265	+ t_index_temp += 2
	266	+ else:
	267	+ idx_current_point = index + ring*num_sides + side
	268	+ idx_next_ring = index + (ring+1)*num_sides + side
	269	+ triangle1 = [idx_current_point, idx_next_ring, idx_next_ring-num_sides+1]
	270	+ triangle2 = [idx_next_ring-num_sides+1, idx_current_point-num_sides+1, idx_current_point]
	271	+ triangles[t_index_temp] = [idx_current_point, idx_next_ring-num_sides, idx_next_ring-num_sides+1]
	272	+ triangles[t_index_temp+1] = [idx_next_ring-num_sides+1, idx_current_point-num_sides+1, idx_current_point]
	273	+ self.triangle_data[t_index] = triangle1
	274	+ self.triangle_data[t_index+1] = triangle2
	275	+ t_index += 2
	276	+ t_index_temp += 2
	277	+
	278	+ #generate the points data structure
	279	+ circle_pts_data = circle_pts.reshape((pts.shape[0]*num_sides, 3))
	280	+ self.cylinder_data['a_position'][index:(pts.shape[0]*num_sides+index)] = circle_pts_data
	281	+ self.cylinder_data['a_fg_color'][index:(pts.shape[0]*num_sides+index)] = color
	282	+ #generate the normals data structure
	283	+ pts_normals = circle_pts.copy()
	284	+ for p in range(pts.shape[0]):
	285	+ pts_normals[p][:] = pts_normals[p][:] - pts[p]
	286	+ for s in range(num_sides):
	287	+ pts_normals[p][s] = pts_normals[p][s]/np.sqrt(pts_normals[p][s][0]*pts_normals[p][s][0] \
	288	+ + pts_normals[p][s][1]pts_normals[p][s][1] + pts_normals[p][s][2]pts_normals[p][s][2])
	289	+ self.cylinder_data['a_normal'][index:(pts.shape[0]*num_sides+index)] = \
	290	+ pts_normals.reshape((pts.shape[0]*num_sides, 3))
	291	+
	292	+ index += pts.shape[0]*num_sides
	293	+
	294	+ #Add the caps for each of the endpoints.
	295	+
	296	+
	297	+
	298	+ if(i == 2):
	299	+ fig = plt.figure()
	300	+ ax = fig.add_subplot(111, projection='3d')
	301	+ #ax.scatter(circle_pts[:,:,0], circle_pts[:,:,1], circle_pts[:,:,2])
	302	+ ax.plot(pts[:,0], pts[:,1], pts[:,2])
	303	+ for j in range(pts.shape[0]):
	304	+ ax.plot(circle_pts[j,:,0], circle_pts[j,:,1], circle_pts[j,:,2])
	305	+ for j in range(triangles.shape[0]):
	306	+ tri = np.zeros((3,4))
	307	+ tri[:,0] = self.cylinder_data['a_position'][triangles[j][0]]
	308	+ tri[:,1] = self.cylinder_data['a_position'][triangles[j][1]]
	309	+ tri[:,2] = self.cylinder_data['a_position'][triangles[j][2]]
	310	+ tri[:,3] = self.cylinder_data['a_position'][triangles[j][0]]
	311	+ ax.plot(tri[0,:], tri[1,:], tri[2,:], c='b')
	312	+ for j in range(triangles.shape[0]):
	313	+ tri = np.zeros((3,3))
	314	+ tri[:,0] = self.cylinder_data['a_position'][triangles[j][0]]
	315	+ tri[:,1] = self.cylinder_data['a_position'][triangles[j][1]]
	316	+ tri[:,2] = self.cylinder_data['a_position'][triangles[j][2]]
	317	+ norm = np.zeros((3,3))
	318	+ norm[:,0] = self.cylinder_data['a_normal'][triangles[j][0]]
	319	+ norm[:,1] = self.cylinder_data['a_normal'][triangles[j][1]]
	320	+ norm[:,2] = self.cylinder_data['a_normal'][triangles[j][2]]
	321	+ ax.quiver(tri[0,:], tri[1,:], tri[2,:], norm[0,:], norm[1,:], norm[2,:], colors = 'r')
	322	+ plt.show()
	323	+ i+=1
	324	+ #create the data.
	325	+
	326	+ #Handles the mouse wheel event, i.e., zoom
	327	+ def on_mouse_wheel(self, event):
	328	+# self.scale[0] = self.scale[0] + self.scale[0]event.delta[1]0.05
	329	+# self.scale[1] = self.scale[1] + self.scale[1]event.delta[1]0.05
	330	+# self.scale[2] = self.scale[2] + self.scale[2]event.delta[1]0.05
	331	+## self.view[0][0] = self.scale[0]
	332	+## self.view[1][1] = self.scale[1]
	333	+## self.view[2][2] = self.scale[2]
	334	+# print("in mouse wheel ", self.r1, self.r2)
	335	+# self.view = np.matmul(translate((self.translate[0], self.translate[1], self.translate[2])), self.r1)
	336	+# self.view = np.matmul(self.view, self.r2)
	337	+# self.view = np.matmul(self.view, scale((self.scale[0], self.scale[1], self.scale[2])))
	338	+# #self.view = np.matmul(self.view, self.r1)
	339	+# #self.view = np.matmul(self.view, self.r2)
	340	+# #self.view = np.matmul(translate((self.translate[0], self.translate[1], self.translate[2])), scale((self.scale[0], self.scale[1], self.scale[2])))
	341	+# #self.rotate = [0., 0.]
	342	+# #self.camera =
	343	+#
	344	+# #self.view[1][1] = self.view[1][1]+self.view[1][1]event.delta[1]0.05
	345	+# #print(self.view[0][0], " ",self.view[1][1])
	346	+# #print(self.view)
	347	+# self.camera = [0.0, 0.0, -100.0, 1.0]
	348	+## for i in range(len(self.translate)):
	349	+## self.camera[i] += self.translate[i]
	350	+# self.program['u_view'] = self.view
	351	+
	352	+# if(np.asarray(self.camera).all() == np.asarray([0., 0., 0.]).all()):
	353	+# self.camera = np.asarray([0., 0., 0.])
	354	+# else:
	355	+ direction = np.asarray(self.translate) - np.asarray(self.camera)
	356	+ direction = direction/np.sqrt(np.dot(direction, direction))
	357	+ for i in range(3):
	358	+ self.camera[i] = self.camera[i] + direction[i]event.delta[1]2.0
	359	+
	360	+ self.program['u_eye'] = self.camera
	361	+ #print(self.view)
	362	+ #print(event.delta[1])
	363	+ self.update()
	364	+
	365	+
	366	+
	367	+ #Handles the mouse press event to rotate the camera if the left mouse button
	368	+ #if clicked down.
	369	+ def on_mouse_press(self, event):
	370	+ def update_view():
	371	+ self.location = event.pos
	372	+ #self.program['u_view'] = self.view
	373	+ self.down = True
	374	+
	375	+ if(event.button == 1):
	376	+ update_view()
	377	+
	378	+
	379	+ #Handles the rotation of the camera using a quaternion centered around the
	380	+ #focus point.
	381	+ def on_mouse_move(self, event):
	382	+ if(self.down == True):
	383	+ coord = self.transforms.get_transform('canvas', 'render').map(self.location)[:2]
	384	+ coord2 = self.transforms.get_transform('canvas', 'render').map(event.pos)[:2]
	385	+ #coord[1] = 0
	386	+ #coord2[1] = 0
	387	+
	388	+ theta = (coord[0]-coord2[0])*360.0/2.0/np.pi
	389	+ phi = (coord[1]-coord2[1])*360.0/2.0/np.pi
	390	+ print(theta360.0/2.0/np.pi, -phi360.0/2.0/np.pi)
	391	+ self.camera = self.camera - self.translate
	392	+ q1 = Quaternion.create_from_axis_angle(angle=phi, ax=0.0, ay=1.0, az=0.0, degrees=True)
	393	+ q2 = Quaternion.create_from_axis_angle(angle=theta, ax=1.0, ay=0.0, az=0.0, degrees=True)
	394	+ #q1 = Quaternion(w=theta, x = 0, y = 1, z = 0).inverse().normalize()
	395	+ #q2 = Quaternion(w=-phi, x = 1, y = 0, z = 0).inverse().normalize()
	396	+
	397	+ q = q1*q2
	398	+
	399	+# #print("Angle in Degrees is ", theta, " ", phi, coord[0] - coord2[0])
	400	+# self.r1 = rotate(theta, (0, 1, 0))
	401	+# self.r2 = rotate(-phi, (1, 0, 0))
	402	+#
	403	+# print("in on mouse move ", self.r1, self.r2)
	404	+#
	405	+# self.view = np.matmul(self.view, self.r1)
	406	+# #print("r1", self.view)
	407	+# self.view = np.matmul(self.view, self.r2)
	408	+# #print("r2", self.view)
	409	+#
	410	+## self.view = np.matmul(self.view, q1.get_matrix().T)
	411	+## self.view = np.matmul(self.view, q2.get_matrix().T)
	412	+#
	413	+# self.program['u_view'] = self.view
	414	+#
	415	+ self.location = event.pos
	416	+# #print("Angle in Degrees is ", theta, " ", phi)
	417	+# #print(self.camera)
	418	+ self.camera = np.asarray(q.rotate_point(self.camera), dtype=np.float)
	419	+ self.camera = self.camera + self.translate
	420	+ self.up = np.asarray(q.rotate_point(self.up), dtype=np.float)
	421	+ self.up = self.up/np.sqrt(np.dot(self.up, self.up))
	422	+ #self.rotate[0] = self.rotate[0] + theta
	423	+ #self.rotate[1] = self.rotate[1] + phi
	424	+ #print(self.rotate)f
	425	+ #radius = np.sqrt(np.dot(self.center, self.center))*2
	426	+ #test = np.linalg.inv(self.view).T
	427	+ #print(test)
	428	+
	429	+
	430	+
	431	+ #self.camera = sph2cart(self.rotate[0]/360.02.0np.pi+np.pi, self.rotate[1]/360.02.0np.pi+np.pi, 1000.0)
	432	+ #self.camera[0] = self.camera[0] + self.center[0]
	433	+ #self.camera[1] = self.camera[1] + self.center[1]
	434	+ #self.camera[2] = self.camera[2] - self.center[2]
	435	+ print("current position ", self.camera, " and up vector ", self.up)
	436	+ self.program['u_eye'] = self.camera
	437	+ self.program['u_up'] = self.up
	438	+ self.program['u_LightPos'] = [self.camera[0], self.camera[1], self.camera[2]]
	439	+ self.update()
	440	+
	441	+ #reverts the mouse state during release.
	442	+ def on_mouse_release(self, event):
	443	+ self.down = False
...	...

TubeWidget.py 0 → 100644

Wrap text Show/Hide comments View file @9f9f178

	1	+++ a/TubeWidget.py
	1	+#!/usr/bin/env python3
	2	+# -- coding: utf-8 --
	3	+"""
	4	+Created on Mon Aug 5 15:53:16 2019
	5	+
	6	+@author: pavel
	7	+"""
	8	+
	9	+"""
	10	+ Qt wrapper/container class that handles all the top level events and maintains the
	11	+ methods necessary to use the QT signals and slots API.
	12	+"""
	13	+
	14	+from pyqtgraph.Qt import QtCore, QtGui, QtWidgets
	15	+from TubeCanvas import TubeDraw
	16	+
	17	+class TubeWidget(QtGui.QWidget):
	18	+ sigUpdate = QtCore.pyqtSignal(float, float, float)
	19	+ #Initializes the QT wrapper class.
	20	+ def __init__(self):
	21	+ super(TubeWidget, self).__init__()
	22	+ box = QtGui.QVBoxLayout(self)
	23	+ self.resize(500,500)
	24	+ self.setLayout(box)
	25	+ self.canvas = TubeDraw()
	26	+ #self.canvas.create_native()
	27	+ box.addWidget(self.canvas.native)
	28	+ self.camera = [0,0,0]
	29	+ self.down = False
	30	+
	31	+ self.canvas.events.mouse_press.connect(self.on_mouse_press)
	32	+ self.canvas.events.mouse_release.connect(self.on_mouse_release)
	33	+ self.canvas.events.mouse_move.connect(self.on_mouse_move)
	34	+ self.canvas.events.mouse_wheel.connect(self.on_mouse_wheel)
	35	+
	36	+ #self.show()
	37	+
	38	+ #Handles the mouse release event
	39	+ def on_mouse_release(self, event):
	40	+ self.down=False
	41	+ self.sendCameraInfo()
	42	+
	43	+ #Handles the mouse move event
	44	+ def on_mouse_move(self, event):
	45	+ if self.down:
	46	+ self.sendCameraInfo()
	47	+
	48	+ #Handles the mouse press event
	49	+ def on_mouse_press(self, event):
	50	+ self.down = True
	51	+ n = 3
	52	+
	53	+ #Handles the mouse wheel event
	54	+ def on_mouse_wheel(self, event):
	55	+ self.sendCameraInfo()
	56	+
	57	+ #controls the emit function of the QT class to send a signal to the slot
	58	+ #located in the other window.
	59	+ def sendCameraInfo(self):
	60	+ #self.view = self.canvas.view
	61	+ self.camera = self.canvas.camera
	62	+ #print("stuff", self.view[3, 0], self.view[3, 1], self.view[3, 2])
	63	+ print("stuff", self.camera[0], self.camera[1], self.camera[2])
	64	+ self.sigUpdate.emit(self.camera[0], self.camera[1], self.camera[2])
	65	+
...	...

graph_shaders.py 0 → 100644

Wrap text Show/Hide comments View file @9f9f178

	1	+++ a/graph_shaders.py
	1	+#!/usr/bin/env python3
	2	+# -- coding: utf-8 --
	3	+"""
	4	+Created on Mon Aug 5 14:48:00 2019
	5	+Collection of shaders necessary to draw all elements of the graph.
	6	+@author: pavel
	7	+"""
	8	+
	9	+from vispy import gloo, app, scene
	10	+
	11	+#shaders for drawing nodes
	12	+
	13	+vert = """
	14	+#version 120
	15	+// Uniforms
	16	+// ------------------------------------
	17	+uniform mat4 u_model;
	18	+uniform mat4 u_view;
	19	+uniform mat4 u_projection;
	20	+uniform float u_antialias;
	21	+uniform float u_size;
	22	+uniform bool u_picking;
	23	+
	24	+uniform vec3 u_graph_size;
	25	+// Attributes
	26	+// ------------------------------------
	27	+attribute vec3 a_position;
	28	+attribute vec4 a_fg_color;
	29	+attribute vec4 a_bg_color;
	30	+attribute float a_linewidth;
	31	+attribute float a_size;
	32	+attribute vec4 a_unique_id;
	33	+attribute float a_selection;
	34	+// Varyings
	35	+// ------------------------------------
	36	+varying vec4 v_fg_color;
	37	+varying vec4 v_bg_color;
	38	+varying float v_size;
	39	+varying float v_linewidth;
	40	+varying float v_antialias;
	41	+varying vec4 v_unique_id;
	42	+varying float v_selection;
	43	+
	44	+varying float v_zoom_level;
	45	+void main (void) {
	46	+ v_size = a_size * u_size;
	47	+ v_linewidth = a_linewidth;
	48	+ v_antialias = u_antialias;
	49	+ v_fg_color = a_fg_color;
	50	+ v_bg_color = a_bg_color;
	51	+ gl_Position = u_projection * u_view * u_model *
	52	+ vec4(a_position*u_size,1.0);
	53	+ gl_PointSize = v_size + 2(v_linewidth + 1.5v_antialias);
	54	+ v_zoom_level = u_view[0][0];
	55	+ v_unique_id = a_unique_id;
	56	+ v_selection = a_selection;
	57	+}
	58	+"""
	59	+
	60	+frag = """
	61	+#version 120
	62	+// Constants
	63	+// ------------------------------------
	64	+uniform mat4 u_model;
	65	+uniform mat4 u_view;
	66	+uniform mat4 u_projection;
	67	+uniform float u_antialias;
	68	+uniform float u_size;
	69	+uniform bool u_picking;
	70	+
	71	+// Varyings
	72	+// ------------------------------------
	73	+varying vec4 v_fg_color;
	74	+varying vec4 v_bg_color;
	75	+varying float v_size;
	76	+varying float v_linewidth;
	77	+varying float v_antialias;
	78	+varying vec4 v_unique_id;
	79	+
	80	+varying float v_zoom_level;
	81	+varying float v_selection;
	82	+// Functions
	83	+// ------------------------------------
	84	+float marker(vec2 P, float size);
	85	+float new_alpha(float zoom_level);
	86	+// Main
	87	+// ------------------------------------
	88	+void main()
	89	+{
	90	+ if(!u_picking)
	91	+ {
	92	+ float size = v_size + 2(v_linewidth + 1.5v_antialias);
	93	+ float t = v_linewidth/2.0-v_antialias;
	94	+ // The marker function needs to be linked with this shader
	95	+ float r = marker(gl_PointCoord, size);
	96	+ float d = abs(r) - t;
	97	+ if( r > (v_linewidth/2.0+v_antialias))
	98	+ {
	99	+ discard;
	100	+ }
	101	+ else if( d < 0.0 )
	102	+ {
	103	+ if(v_zoom_level < 0.0010)
	104	+ {
	105	+ float alpha = d/v_antialias;
	106	+ alpha = new_alpha(v_zoom_level);
	107	+ gl_FragColor = mix(vec4(1,1,1,0.9), v_bg_color, max(1.0-alpha, 0.3));
	108	+ }
	109	+ else
	110	+ {
	111	+ if(v_selection == 1.0)
	112	+ gl_FragColor = vec4(1.0, 0.0, 0.0, v_fg_color.a);
	113	+ else if(v_selection == 2.0)
	114	+ gl_FragColor = vec4(0.0, 1.0, 0.0, v_fg_color.a);
	115	+ else
	116	+ gl_FragColor = v_fg_color;
	117	+ }
	118	+ }
	119	+ else
	120	+ {
	121	+ float alpha = d/v_antialias;
	122	+ if(v_zoom_level < 0.0010)
	123	+ alpha = new_alpha(v_zoom_level);
	124	+ else
	125	+ alpha = exp(-alpha*alpha);
	126	+
	127	+ if (r > 0 && v_zoom_level >= 0.0010)
	128	+ gl_FragColor = vec4(v_fg_color.rgb, alpha*v_fg_color.a);
	129	+ else
	130	+ if(v_zoom_level < 0.0010)
	131	+ if(vec3(gl_Color.rgb) != vec3(v_bg_color.rgb))
	132	+ gl_FragColor = mix(vec4(1,1,1,0.9), v_bg_color, max(1.0-alpha, 0.3));
	133	+ else
	134	+ gl_FragColor = vec4(gl_Color.rgb, max(1.0-alpha, 0.1));
	135	+ else
	136	+ gl_FragColor = mix(v_bg_color, v_fg_color, alpha);
	137	+ }
	138	+
	139	+
	140	+ } else {
	141	+ float size = v_size +2(v_linewidth + 1.5v_antialias);
	142	+ float t = v_linewidth/2.0-v_antialias;
	143	+ // The marker function needs to be linked with this shader
	144	+ float r = marker(gl_PointCoord, size);
	145	+ float d = abs(r) - t;
	146	+ float alpha = d/v_antialias;
	147	+ if(v_zoom_level < 0.0010)
	148	+ alpha = new_alpha(v_zoom_level);
	149	+ else
	150	+ alpha = exp(-alpha*alpha);
	151	+ if(r > 0)
	152	+ gl_FragColor = v_unique_id;
	153	+ else
	154	+ gl_FragColor = v_unique_id;
	155	+
	156	+ }
	157	+}
	158	+
	159	+float marker(vec2 P, float size)
	160	+{
	161	+ float r = length((P.xy - vec2(0.5,0.5))*size);
	162	+ r -= v_size/2;
	163	+ return r;
	164	+}
	165	+
	166	+float new_alpha(float zoom_level)
	167	+{
	168	+ float val = (zoom_level - 0.0010)/(0.00075-0.0010);
	169	+ if(val < 0)
	170	+ {
	171	+ val = 0;
	172	+ }
	173	+ return val;
	174	+}
	175	+"""
	176	+
	177	+#Shaders for the lines between the nodes
	178	+
	179	+
	180	+vs = """
	181	+#version 120
	182	+
	183	+// Uniforms
	184	+// ------------------------------------
	185	+uniform mat4 u_model;
	186	+uniform mat4 u_view;
	187	+uniform mat4 u_projection;
	188	+
	189	+// Attributes
	190	+// ------------------------------------
	191	+attribute vec3 a_position;
	192	+attribute vec2 a_normal;
	193	+attribute vec4 a_fg_color;
	194	+attribute float a_linewidth;
	195	+//attribute vec4 a_unique_id;
	196	+//attribute vec4 l_color;
	197	+
	198	+// Varyings
	199	+// ------------------------------------
	200	+varying vec4 v_fg_color;
	201	+varying float v_zoom_level;
	202	+varying vec2 v_normal;
	203	+varying float v_linewidth;
	204	+
	205	+void main() {
	206	+ vec3 delta = vec3(a_normal*a_linewidth/(1-u_view[0][0]), 0);
	207	+ //vec4 pos = u_model * u_view * vec4(a_position, 1.0);
	208	+ //gl_Position = u_projection * (pos + vec4(delta, 1.0));
	209	+ gl_Position = u_model * u_view * u_projection * vec4(a_position+delta, 1.0);
	210	+ //gl_Position = u_projection * u_view * u_model *
	211	+ // vec4(a_position, 1.0);
	212	+ v_zoom_level = u_view[0][0];
	213	+ v_fg_color = a_fg_color;
	214	+ v_normal = a_normal;
	215	+ v_linewidth = a_linewidth;
	216	+
	217	+}
	218	+"""
	219	+
	220	+fs = """
	221	+#version 120
	222	+// Varying
	223	+// ------------------------------------
	224	+varying vec4 v_fg_color;
	225	+varying float v_zoom_level;
	226	+varying vec2 v_normal;
	227	+varying float v_linewidth;
	228	+
	229	+float new_alpha(float zoom_level);
	230	+
	231	+void main()
	232	+{
	233	+ float l = length(v_normal);
	234	+ float feather = 0.5;
	235	+ float alpha = 1.0;
	236	+ if(l > v_linewidth/2.0+feather)
	237	+ discard;
	238	+ else
	239	+ alpha = 0.5;
	240	+
	241	+ if(v_zoom_level < 0.0010)
	242	+ alpha = new_alpha(v_zoom_level);
	243	+ gl_FragColor = vec4(v_fg_color.rgb, alpha);
	244	+}
	245	+
	246	+float new_alpha(float zoom_level)
	247	+{
	248	+ float val = (zoom_level-0.00075)/(0.0010-0.00075);
	249	+ if(val < 0.)
	250	+ {
	251	+ val = 0.;
	252	+ }
	253	+ return val;
	254	+
	255	+}
	256	+"""
0	257	\ No newline at end of file
...	...

network_dep.py 0 → 100755

Wrap text Show/Hide comments View file @9f9f178

	1	+++ a/network_dep.py
	1	+# -- coding: utf-8 --
	2	+"""
	3	+Created on Sat Sep 16 16:34:49 2017
	4	+
	5	+@author: pavel
	6	+"""
	7	+
	8	+import struct
	9	+import numpy as np
	10	+from scipy.sparse.linalg import eigsh
	11	+import scipy as sp
	12	+from sklearn.cluster import SpectralClustering
	13	+
	14	+
	15	+#import networkx as nx
	16	+import matplotlib.pyplot as plt
	17	+from matplotlib import cm
	18	+import math
	19	+import time
	20	+#import spharmonics
	21	+import graph_tool.all as gt
	22	+import copy
	23	+#import matplotlib
	24	+
	25	+#for testing
	26	+#import matplotlib.mlab as mlab
	27	+#import matplotlib.pyplot as plt
	28	+#from matplotlib import cm
	29	+
	30	+'''
	31	+ Definition of the Node class
	32	+ Duplicate of the node class in network
	33	+ Stores the physical position, outgoing edges list and incoming edges list.
	34	+'''
	35	+class Node:
	36	+ def __init__(self, point, outgoing, incoming):
	37	+ self.p = point
	38	+ self.o = outgoing
	39	+ self.i = incoming
	40	+
	41	+# def p():
	42	+# return self.p
	43	+
	44	+'''
	45	+ Definition of the Fiber class.
	46	+ Duplicate of the Node class in network
	47	+ Stores the starting vertex, the ending vertex, the points array and the radius array
	48	+'''
	49	+class Fiber:
	50	+
	51	+
	52	+ def __init__ (self, p1, p2, pois, rads):
	53	+ self.v0 = p1
	54	+ self.v1 = p2
	55	+ self.points = pois
	56	+ self.radii = rads
	57	+ '''
	58	+ return the array-likes of the x,y,z,r coordinates of the fiber for the
	59	+ gt representation
	60	+ '''
	61	+
	62	+ def getcoords(self):
	63	+ x = np.zeros(len(self.points), dtype=np.double)
	64	+ y = np.zeros(len(self.points), dtype=np.double)
	65	+ z = np.zeros(len(self.points), dtype=np.double)
	66	+ r = np.zeros(len(self.points), dtype=np.double)
	67	+ for i in range(len(self.points)):
	68	+ x[i] = self.points[i][0]
	69	+ y[i] = self.points[i][1]
	70	+ z[i] = self.points[i][2]
	71	+ r[i] = self.radii[i]
	72	+
	73	+ return x,y,z,r
	74	+
	75	+
	76	+ '''
	77	+ return the length of the fiber.
	78	+ '''
	79	+ def length(self):
	80	+ length = 0
	81	+ for i in range(len(self.points)-1):
	82	+ 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))
	83	+ if(length == 0):
	84	+ print(self.points)
	85	+ print(len(self.points))
	86	+ print(self.v0, " ", self.v1)
	87	+ return length
	88	+
	89	+ '''
	90	+ returns the turtuosity of the fiber.
	91	+ '''
	92	+ def turtuosity(self):
	93	+ turtuosity = 1.0
	94	+ 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))
	95	+ if(distance > 0):
	96	+ turtuosity = self.length()/distance
	97	+ #print(turtuosity)
	98	+
	99	+ return turtuosity
	100	+
	101	+ '''
	102	+ returns the volume of the fiber.
	103	+ '''
	104	+ def volume(self):
	105	+ volume = 0
	106	+ for i in range(len(self.points)-1):
	107	+ 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-6self.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))
	108	+
	109	+ #print(volume)
	110	+ return volume
	111	+
	112	+ def av_radius(self):
	113	+ radius = 0.0
	114	+ for i in range(len(self.radii)):
	115	+ radius = radius + self.radii[i]
	116	+
	117	+ return radius/len(self.radii)
	118	+
	119	+class NWT:
	120	+
	121	+ '''
	122	+ Writes the header given and open file descripion, number of verticies and number of edges.
	123	+ '''
	124	+ def writeHeader(open_file, numVerts, numEdges):
	125	+ txt = "nwtFileFormat fileid(14B), desc(58B), #vertices(4B), #edges(4B): bindata"
	126	+ b = bytearray()
	127	+ b.extend(txt.encode())
	128	+ open_file.write(b)
	129	+ open_file.write(struct.pack('i', numVerts))
	130	+ open_file.write(struct.pack('i', numEdges))
	131	+
	132	+
	133	+ '''
	134	+ Writes a single vertex to a file.
	135	+ '''
	136	+ def writeVertex(open_file, vertex):
	137	+ open_file.write(struct.pack('<f',vertex.p[0]))
	138	+ open_file.write(struct.pack('<f',vertex.p[1]))
	139	+ open_file.write(struct.pack('<f',vertex.p[2]))
	140	+ open_file.write(struct.pack('i', len(vertex.o)))
	141	+ open_file.write(struct.pack('i', len(vertex.i)))
	142	+ for j in range(len(vertex.o)):
	143	+ open_file.write(struct.pack('i',vertex.o[j]))
	144	+
	145	+ for j in range(len(vertex.i)):
	146	+ open_file.write(struct.pack('i', vertex.i[j]))
	147	+
	148	+ return
	149	+
	150	+ '''
	151	+ Writes a single fiber to a file.
	152	+ '''
	153	+ def writeFiber(open_file, edge):
	154	+ open_file.write(struct.pack('i',edge.v0))
	155	+ open_file.write(struct.pack('i',edge.v1))
	156	+ open_file.write(struct.pack('i',len(edge.points)))
	157	+ for j in range(len(edge.points)):
	158	+ open_file.write(struct.pack('<f', edge.points[j][0]))
	159	+ open_file.write(struct.pack('<f', edge.points[j][1]))
	160	+ open_file.write(struct.pack('<f', edge.points[j][2]))
	161	+ open_file.write(struct.pack('<f', edge.radii[j]))
	162	+
	163	+ return
	164	+
	165	+ '''
	166	+ Writes the entire network to a file in str given the vertices array and the edges array.
	167	+ '''
	168	+ def exportNWT(str, vertices, edges):
	169	+ with open(str, "wb") as file:
	170	+ NWT.writeHeader(file, len(vertices), len(edges))
	171	+ for i in range(len(vertices)):
	172	+ NWT.writeVertex(file, vertices[i])
	173	+
	174	+ for i in range(len(edges)):
	175	+ NWT.writeFiber(file, edges[i])
	176	+
	177	+ return
	178	+
	179	+
	180	+ '''
	181	+ Reads a single vertex from an open file and returns a node Object.
	182	+ '''
	183	+ def readVertex(open_file):
	184	+ points = np.tile(0., 3)
	185	+ bytes = open_file.read(4)
	186	+ points[0] = struct.unpack('f', bytes)[0]
	187	+ bytes = open_file.read(4)
	188	+ points[1] = struct.unpack('f', bytes)[0]
	189	+ bytes = open_file.read(4)
	190	+ points[2] = struct.unpack('f', bytes)[0]
	191	+ bytes = open_file.read(4)
	192	+
	193	+ numO = int.from_bytes(bytes, byteorder='little')
	194	+ outgoing = np.tile(0, numO)
	195	+ bts = open_file.read(4)
	196	+ numI = int.from_bytes(bts, byteorder='little')
	197	+ incoming = np.tile(0, numI)
	198	+ for j in range(numO):
	199	+ bytes = open_file.read(4)
	200	+ outgoing[j] = int.from_bytes(bytes, byteorder='little')
	201	+
	202	+ for j in range(numI):
	203	+ bytes = open_file.read(4)
	204	+ incoming[j] = int.from_bytes(bytes, byteorder='little')
	205	+
	206	+ node = Node(points, outgoing, incoming)
	207	+ return node
	208	+
	209	+
	210	+ '''
	211	+ Reads a single fiber from an open file and returns a Fiber object .
	212	+ '''
	213	+ def readFiber(open_file):
	214	+ bytes = open_file.read(4)
	215	+ vtx0 = int.from_bytes(bytes, byteorder = 'little')
	216	+ bytes = open_file.read(4)
	217	+ vtx1 = int.from_bytes(bytes, byteorder = 'little')
	218	+ bytes = open_file.read(4)
	219	+ numVerts = int.from_bytes(bytes, byteorder = 'little')
	220	+ pts = []
	221	+ rads = []
	222	+
	223	+ for j in range(numVerts):
	224	+ point = np.tile(0., 3)
	225	+ bytes = open_file.read(4)
	226	+ point[0] = struct.unpack('f', bytes)[0]
	227	+ bytes = open_file.read(4)
	228	+ point[1] = struct.unpack('f', bytes)[0]
	229	+ bytes = open_file.read(4)
	230	+ point[2] = struct.unpack('f', bytes)[0]
	231	+ bytes = open_file.read(4)
	232	+ radius = struct.unpack('f', bytes)[0]
	233	+ pts.append(point)
	234	+ rads.append(radius)
	235	+
	236	+ F = Fiber(vtx0, vtx1, pts, rads)
	237	+
	238	+ return F
	239	+
	240	+ '''
	241	+ Imports a NWT file at location str.
	242	+ Returns a list of Nodes objects and a list of Fiber objects.
	243	+ '''
	244	+
	245	+'''
	246	+Class defining an aabb around the graph
	247	+'''
	248	+
	249	+class AABB():
	250	+ def __init__(self, G, is_dual=False):
	251	+ #designate whether we will be generating a bounding box for a dual graph
	252	+ # or a normal graph.
	253	+ self.is_dual = is_dual
	254	+ #minimum vertex
	255	+ self.A = np.full((3,1), 1000000.0, dtype=float)
	256	+ #maximum vertex
	257	+ self.B = np.full((3,1), -1000000.0, dtype=float)
	258	+ if(is_dual == False):
	259	+ #find the minumum and the maximum of the graph.
	260	+ for v in G.vertices():
	261	+ if G.vertex_properties["p"][v][0] < self.A[0]:
	262	+ self.A[0] = G.vertex_properties["p"][v][0]
	263	+ if G.vertex_properties["p"][v][0] > self.B[0]:
	264	+ self.B[0] = G.vertex_properties["p"][v][0]
	265	+ if G.vertex_properties["p"][v][1] < self.A[1]:
	266	+ self.A[1] = G.vertex_properties["p"][v][1]
	267	+ if G.vertex_properties["p"][v][1] > self.B[1]:
	268	+ self.B[1] = G.vertex_properties["p"][v][1]
	269	+ if G.vertex_properties["p"][v][2] < self.A[2]:
	270	+ self.A[2] = G.vertex_properties["p"][v][2]
	271	+ if G.vertex_properties["p"][v][2] > self.B[2]:
	272	+ self.B[2] = G.vertex_properties["p"][v][2]
	273	+ #In case of a dual graph we have to scane the first and last points
	274	+ # of every fiber to find the true bounding box
	275	+ else:
	276	+ for v in G.vertices():
	277	+ l = len(G.vertex_properties["x"][v])
	278	+ if G.vertex_properties["x"][v][0] < self.A[0]:
	279	+ self.A[0] = G.vertex_properties["x"][v][0]
	280	+ if G.vertex_properties["x"][v][l-1] < self.A[0]:
	281	+ self.A[0] = G.vertex_properties["x"][v][l-1]
	282	+ if G.vertex_properties["y"][v][0] < self.A[1]:
	283	+ self.A[1] = G.vertex_properties["y"][v][0]
	284	+ if G.vertex_properties["y"][v][l-1] < self.A[1]:
	285	+ self.A[1] = G.vertex_properties["y"][v][l-1]
	286	+ if G.vertex_properties["z"][v][0] < self.A[2]:
	287	+ self.A[2] = G.vertex_properties["z"][v][0]
	288	+ if G.vertex_properties["z"][v][l-1] < self.A[2]:
	289	+ self.A[2] = G.vertex_properties["z"][v][l-1]
	290	+
	291	+ if G.vertex_properties["x"][v][0] > self.B[0]:
	292	+ self.B[0] = G.vertex_properties["x"][v][0]
	293	+ if G.vertex_properties["x"][v][l-1] > self.B[0]:
	294	+ self.B[0] = G.vertex_properties["x"][v][l-1]
	295	+ if G.vertex_properties["y"][v][0] > self.B[1]:
	296	+ self.B[1] = G.vertex_properties["y"][v][0]
	297	+ if G.vertex_properties["y"][v][l-1] > self.B[1]:
	298	+ self.B[1] = G.vertex_properties["y"][v][l-1]
	299	+ if G.vertex_properties["z"][v][0] > self.B[2]:
	300	+ self.B[2] = G.vertex_properties["z"][v][0]
	301	+ if G.vertex_properties["z"][v][l-1] > self.B[2]:
	302	+ self.B[2] = G.vertex_properties["z"][v][l-1]
	303	+ #print(self.A, self.B)
	304	+ self.O = (self.A+self.B)*0.5
	305	+ self.planes = []
	306	+ #append x planes described by a point and a vector
	307	+ self.planes.append((np.array([self.A[0], 0.0, 0.0]), np.array([1.0, 0.0, 0.0])))
	308	+ self.planes.append((np.array([self.B[0], 0.0, 0.0]), np.array([-1.0, 0.0, 0.0])))
	309	+ #append y planes described by a point and a vector
	310	+ self.planes.append((np.array([0.0, self.A[1], 0.0]), np.array([0.0, 1.0, 0.0])))
	311	+ self.planes.append((np.array([0.0, self.B[1], 0.0]), np.array([0.0, -1.0, 0.0])))
	312	+ #append z planes described by a point and a vector
	313	+ self.planes.append((np.array([0.0, 0.0, self.A[2]]), np.array([0.0, 0.0, 1.0])))
	314	+ self.planes.append((np.array([0.0, 0.0, self.B[2]]), np.array([0.0, 0.0, -1.0])))
	315	+
	316	+
	317	+ def distance(self, pt):
	318	+ dist = 10000000
	319	+ for i in self.planes:
	320	+ V = i[0]-pt
	321	+ if np.dot(V, i[1]) < dist:
	322	+ dist = np.dot(V, i[1])
	323	+
	324	+ return dist
	325	+
	326	+ def project_grid(self, n):
	327	+ #r = abs(self.A - self.B)
	328	+ x = np.linspace(self.A[0], self.B[0], (n+2))
	329	+ x = x[1:-1]
	330	+ y = np.linspace(self.A[1], self.B[1], (n+2))
	331	+ y = y[1:-1]
	332	+ z = np.linspace(self.A[2], self.B[2], (n+2))
	333	+ z = z[1:-1]
	334	+
	335	+ return x,y,z
	336	+
	337	+ #returns a resampled bounding both with nxn points on each side.
	338	+ def resample_sides(self, n):
	339	+ #get the size of the cube in every cardinal direction
	340	+ size = abs(self.A - self.B)
	341	+ #set the stepsize for each subdivided point
	342	+ dist_x = size[0]/(n-1)
	343	+ dist_y = size[1]/(n-1)
	344	+ dist_z = size[2]/(n-1)
	345	+ #generate the original 8 corners
	346	+ vertices = []
	347	+ #generate the points for the yz planes.
	348	+ for i in range(n):
	349	+ for j in range(n):
	350	+ #generate the temporary vectors for both the planes
	351	+ temp_p1 = copy.deepcopy(self.A)
	352	+ temp_p2 = copy.deepcopy(self.A)
	353	+ temp_p2[0] = self.B[0]
	354	+
	355	+ temp_p1[1] = temp_p1[1] + dist_y*i
	356	+ temp_p1[2] = temp_p1[2] + dist_z*j
	357	+
	358	+ temp_p2[1] = temp_p2[1] + dist_y*i
	359	+ temp_p2[2] = temp_p2[2] + dist_z*j
	360	+
	361	+ vertices.append(temp_p1)
	362	+ vertices.append(temp_p2)
	363	+ #generate the points for the xz planes.
	364	+ for i in range(n):
	365	+ for j in range(n):
	366	+ #generate the temporary vectors for both the planes
	367	+ temp_p1 = copy.deepcopy(self.A)
	368	+ temp_p2 = copy.deepcopy(self.A)
	369	+ temp_p2[1] = self.B[1]
	370	+
	371	+ temp_p1[0] = temp_p1[0] + dist_x*i
	372	+ temp_p1[2] = temp_p1[2] + dist_z*j
	373	+
	374	+ temp_p2[0] = temp_p2[0] + dist_x*i
	375	+ temp_p2[2] = temp_p2[2] + dist_z*j
	376	+
	377	+ vertices.append(temp_p1)
	378	+ vertices.append(temp_p2)
	379	+
	380	+ #generate the points for the xy planes.
	381	+ for i in range(n):
	382	+ for j in range(n):
	383	+ #generate the temporary vectors for both the planes
	384	+ temp_p1 = copy.deepcopy(self.A)
	385	+ temp_p2 = copy.deepcopy(self.A)
	386	+ temp_p2[2] = self.B[2]
	387	+
	388	+ temp_p1[0] = temp_p1[0] + dist_x*i
	389	+ temp_p1[1] = temp_p1[1] + dist_y*j
	390	+
	391	+ temp_p2[0] = temp_p2[0] + dist_x*i
	392	+ temp_p2[1] = temp_p2[1] + dist_y*j
	393	+
	394	+ vertices.append(temp_p1)
	395	+ vertices.append(temp_p2)
	396	+
	397	+
	398	+ vertices = list(np.unique(np.array(vertices), axis=0))
	399	+ import matplotlib.pyplot as plt
	400	+ fig = plt.figure()
	401	+ ax = plt.axes(projection='3d')
	402	+ ax.scatter3D(np.array(vertices)[:, 0], np.array(vertices)[:, 1], np.array(vertices)[:, 2])
	403	+
	404	+ print("THERE ARE THIS MANY ", len(vertices))
	405	+ return vertices
	406	+
	407	+ def getVolume(self):
	408	+ size = abs(self.A - self.B)
	409	+ return size[0]size[1]size[2]
	410	+
	411	+ def vertices(self):
	412	+ size = abs(self.A - self.B)
	413	+ points = []
	414	+ points.append(self.A)
	415	+
	416	+
	417	+ temp = copy.deepcopy(self.A)
	418	+ temp[0] = temp[0] + size[0]
	419	+ points.append(temp)
	420	+ print('1', temp)
	421	+
	422	+ temp = copy.deepcopy(self.A)
	423	+ temp[1] = temp[1] + size[1]
	424	+ points.append(temp)
	425	+ print('1', temp)
	426	+
	427	+ temp = copy.deepcopy(self.A)
	428	+ temp[2] = temp[2] + size[2]
	429	+ points.append(temp)
	430	+ print('1', temp)
	431	+
	432	+ temp = copy.deepcopy(self.A)
	433	+ temp[1] = temp[1] + size[1]
	434	+ temp[0] = temp[0] + size[0]
	435	+ points.append(temp)
	436	+ print('1', temp)
	437	+
	438	+ temp = copy.deepcopy(self.A)
	439	+ temp[2] = temp[2] + size[2]
	440	+ temp[0] = temp[0] + size[0]
	441	+ points.append(temp)
	442	+ print('1', temp)
	443	+
	444	+ temp = copy.deepcopy(self.A)
	445	+ temp[1] = temp[1] + size[1]
	446	+ temp[2] = temp[2] + size[2]
	447	+ points.append(temp)
	448	+ print('1', temp)
	449	+
	450	+ temp = copy.deepcopy(self.A)
	451	+ temp[0] = temp[0] + size[0]
	452	+ temp[1] = temp[1] + size[1]
	453	+ temp[2] = temp[2] + size[2]
	454	+ points.append(temp)
	455	+ print('1', temp)
	456	+
	457	+ return points
	458	+
	459	+
	460	+'''
	461	+Class defining the BFS traversal of all the vertices in the graph and labeling
	462	+all interconnected vertices.
	463	+'''
	464	+class VisitorClassDisconnected(gt.BFSVisitor):
	465	+ def __init__(self, ecluster, cluster, c):
	466	+ self.ecluster = ecluster
	467	+ self.cluster = cluster
	468	+ self.c = c
	469	+
	470	+ def examine_vertex(self, u):
	471	+ if(self.cluster[u] == -1):
	472	+ self.cluster[u] = self.c
	473	+
	474	+ def examine_edge(self, e):
	475	+ if(self.ecluster[e] == -1):
	476	+ self.ecluster[e] = self.c
	477	+
	478	+
	479	+'''
	480	+Class defining the BFS traversal of all the vertices in the graph and labeling
	481	+all interconnected vertices.
	482	+'''
	483	+class VisitorClassPartition(gt.BFSVisitor):
	484	+ def __init__(self, cluster, dist, c, iterations):
	485	+ self.cluster = cluster
	486	+ self.dist = dist
	487	+ self.c = c
	488	+ self.total_iter = iterations
	489	+
	490	+ def examine_vertex(self, u):
	491	+ if(self.cluster[u] == -1):
	492	+ self.cluster[u] = self.c
	493	+
	494	+ def tree_edge(self, e):
	495	+ d = self.dist[e.source()]
	496	+ #print(d, self.total_iter)
	497	+ if d <= self.total_iter:
	498	+ self.dist[e.target()] = self.dist[e.source()] + 1
	499	+
	500	+class Network:
	501	+ def __init__(self, filename, clock=False):
	502	+ if clock:
	503	+ start_time = time.time()
	504	+
	505	+ with open(filename, "rb") as file:
	506	+ header = file.read(72)
	507	+ bytes = file.read(4)
	508	+ numVertex = int.from_bytes(bytes, byteorder='little')
	509	+ bytes = file.read(4)
	510	+ numEdges = int.from_bytes(bytes, byteorder='little')
	511	+
	512	+ self.N = []
	513	+ self.F = []
	514	+ for i in range(numVertex):
	515	+ node = NWT.readVertex(file)
	516	+ self.N.append(node)
	517	+
	518	+ for i in range(numEdges):
	519	+ edge = NWT.readFiber(file)
	520	+ self.F.append(edge)
	521	+ if clock:
	522	+ print("Network initialization: " + str(time.time() - start_time) + "s")
	523	+ '''
	524	+ Takes in a graph object (networkx Graph) and saves it as a network file
	525	+ Resulting nwt, does not contain incomind or outgoing fibers
	526	+ '''
	527	+ def saveGraph(self, G, path):
	528	+ G_nodes = list(G.nodes())
	529	+ G_edges = list(G.edges())
	530	+ Nodes = []
	531	+ Edges = []
	532	+ points = list(nx.get_node_attributes(G, 'p').values())
	533	+ for i in range(len(G_nodes)):
	534	+ node = Node(points[i], [], [])
	535	+ Nodes.append(node)
	536	+ for e in G_edges:
	537	+ edge = Fiber(e[0], e[1], G[e[0]][e[1]]['pts'], G[e[0]][e[1]]['rads'])
	538	+ Edges.append(edge)
	539	+
	540	+ NWT.exportNWT(path, Nodes, Edges)
	541	+
	542	+ '''
	543	+ Saves the graph-tool graph as a nwt file with all the variables.
	544	+ '''
	545	+
	546	+ def saveGraph_gt(self, G, path):
	547	+ points = G.vertex_properties["p"]
	548	+ Nodes = []
	549	+ Edges = []
	550	+ for i in range(G.num_vertices()):
	551	+ #array = G.get_out_edges(i)
	552	+ outgoing = []
	553	+ for edge in G.get_out_edges(i):
	554	+ outgoing.append(edge[2])
	555	+
	556	+ node = Node(points[i], outgoing, outgoing)
	557	+ Nodes.append(node)
	558	+
	559	+ for e in G.edges():
	560	+ x = G.edge_properties["x"][e].get_array()
	561	+ y = G.edge_properties["y"][e].get_array()
	562	+ z = G.edge_properties["z"][e].get_array()
	563	+ r = G.edge_properties["r"][e].get_array()
	564	+ pts = []
	565	+ for i in range(len(x)):
	566	+ pt = []
	567	+ pt.append(x[i])
	568	+ pt.append(y[i])
	569	+ pt.append(z[i])
	570	+ pts.append(pt)
	571	+ #pts = np.ndarray.tolist(np.dstack([x,y,z]))
	572	+ #print(pts)
	573	+ edge = Fiber(int(e.source()), int(e.target()), pts, np.ndarray.tolist(r))
	574	+ if(edge.length() > 0.0):
	575	+ Edges.append(edge)
	576	+
	577	+ NWT.exportNWT(path, Nodes, Edges)
	578	+
	579	+
	580	+
	581	+ def saveGraphStatistics_gt(self, G, path, prefix, label = "none", norm=False):
	582	+ location = path + "/" + prefix + "_edges.txt"
	583	+ f = open(location, "a+")
	584	+ f.write("inverted\t")
	585	+ f.write("length\t")
	586	+ f.write("tortuosity\t")
	587	+ f.write("volume\t")
	588	+ f.write("inverted_volume\t")
	589	+ f.write("gaussian\t")
	590	+ f.write("bc\t")
	591	+ f.write("bc*length\t")
	592	+ f.write("mst\t\n")
	593	+ for e in G.edges():
	594	+ f.write("%.15f\t" % G.edge_properties["inverted"][e])
	595	+ f.write("%.15f\t" % G.edge_properties["length"][e])
	596	+ f.write("%.15f\t" % G.edge_properties["tortuosity"][e])
	597	+ f.write("%.15f\t" % G.edge_properties["volume"][e])
	598	+ f.write("%.15f\t" % G.edge_properties["inverted_volume"][e])
	599	+ f.write("%.15f\t" % G.edge_properties["gaussian"][e])
	600	+ f.write("%.15f\t" % G.edge_properties["bc"][e])
	601	+ f.write("%.15f\t" % G.edge_properties["bc*length"][e])
	602	+ f.write("%.15f\t" % G.edge_properties["mst"][e])
	603	+ f.write("%s\n" % label)
	604	+
	605	+ f.close()
	606	+
	607	+ location = path+"/" + prefix + "_vertices.txt"
	608	+ f = open(location, "a+")
	609	+ f.write("bc\t")
	610	+ f.write("degree\t")
	611	+ f.write("degree_volume\t")
	612	+ f.write("degree_tortuosity\t\n")
	613	+ for v in G.vertices():
	614	+ f.write("%.15f\t" % G.vertex_properties["bc"][v])
	615	+ f.write("%.15f\t" % G.vertex_properties["degree"][v])
	616	+ f.write("%.15f\t" % G.vertex_properties["degree_volume"][v])
	617	+ f.write("%.15f\t" % G.vertex_properties["degree_tortuosity"][v])
	618	+ f.write("%s\n" % label)
	619	+
	620	+ '''
	621	+ Saves the graph as a sample. Additionally saves the key as a list of string if writeKey is True
	622	+ saves the file at path with prefix_s.txt
	623	+ '''
	624	+ def saveSample_gt(self, G, path, prefix, is_dual, label = "none", norm=False, writeKey=False):
	625	+ location = path + "/" + prefix + "s.txt"
	626	+ f = open(location, "a+")
	627	+
	628	+
	629	+ aabb = AABB(G, is_dual)
	630	+ array = G.vertex_properties["bc"].get_array()
	631	+ av_array = np.sum(array)/G.num_vertices()
	632	+ G.graph_properties["global_vertex_bc"] = G.new_graph_property("double", val = av_array)
	633	+ f.write("%.15f\t" % av_array)
	634	+
	635	+ array = G.vertex_properties["bc"].get_array()
	636	+ av_array = np.sum(array)/aabb.getVolume()
	637	+ G.graph_properties["global_vascular_vertex_bc"] = G.new_graph_property("double", val = av_array)
	638	+ f.write("%.15f\t" % av_array)
	639	+
	640	+ array = G.edge_properties["bc"].get_array()
	641	+ av_array = np.sum(array)/aabb.getVolume()
	642	+ G.graph_properties["global_vascular_edge_bc"] = G.new_graph_property("double", val = av_array)
	643	+ f.write("%.15f\t" % av_array)
	644	+
	645	+ #G.edge_properties["mst"] = gt.graph_tool.topology.min_spanning_tree(G)
	646	+ 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()))
	647	+ f.write("%.15f\t" % G.graph_properties["global_mst_ratio"])
	648	+
	649	+ #G.edge_properties["mst"] = gt.graph_tool.topology.min_spanning_tree(G, weights=G.edge_properties["volume"])
	650	+ 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())
	651	+ f.write("%.15f\t" % G.graph_properties["global_mst_ratio"])
	652	+
	653	+ ########################Local Statistics###############################
	654	+ self.recalculate_metrics(G, is_dual, norm)
	655	+ ########################Biological Statistics##########################
	656	+ # Some biological statistics are reported as a ratio to physical volume of vessels
	657	+ array = G.edge_properties["length"].get_array()
	658	+ av_array = np.sum(array)/G.num_edges()
	659	+ G.graph_properties["average_length"] = G.new_graph_property("double", val= av_array)
	660	+ f.write("%.15f\t" % av_array)
	661	+
	662	+ array = G.edge_properties["volume"].get_array()
	663	+ av_array = np.sum(array)/G.num_edges()
	664	+ G.graph_properties["average_volume"] = G.new_graph_property("double", val= av_array)
	665	+ f.write("%.15f\t" % av_array)
	666	+
	667	+ array = G.edge_properties["tortuosity"].get_array()
	668	+ av_array = np.sum(array)/G.num_edges()
	669	+ G.graph_properties["average_tortuosity"] = G.new_graph_property("double", val= av_array)
	670	+ f.write("%.15f\t" % av_array)
	671	+
	672	+ array = G.edge_properties["length"].get_array()
	673	+ av_array = np.sum(array)/aabb.getVolume()
	674	+ G.graph_properties["vascular length"] = G.new_graph_property("double", val= av_array)
	675	+ f.write("%.15f\t" % av_array)
	676	+
	677	+ array = G.edge_properties["volume"].get_array()
	678	+ av_array = np.sum(array)/aabb.getVolume()
	679	+ G.graph_properties["vascular_volume"] = G.new_graph_property("double", val= av_array)
	680	+ f.write("%.15f\t" % av_array)
	681	+
	682	+ array = G.edge_properties["tortuosity"].get_array()
	683	+ av_array = np.sum(array)/aabb.getVolume()
	684	+ G.graph_properties["vascular_tortuosity"] = G.new_graph_property("double", val= av_array)
	685	+ f.write("%.15f\t" % av_array)
	686	+ ############################Graph Statistics###########################
	687	+
	688	+ # Some Graph Statistics are reported as values representing the graph
	689	+ # Some statistics are reported as a function of volume.
	690	+ array = G.vertex_properties["degree"].get_array()
	691	+ av_array = np.sum(array)/G.num_vertices()
	692	+ G.graph_properties["local_average_degree"] = G.new_graph_property("double", val = av_array)
	693	+ f.write("%.15f\t" % av_array)
	694	+
	695	+ array = G.vertex_properties["degree_volume"].get_array()
	696	+ av_array = np.sum(array)/G.num_vertices()
	697	+ G.graph_properties["local_average_degree_volume"] = G.new_graph_property("double", val = av_array)
	698	+ f.write("%.15f\t" % av_array)
	699	+
	700	+ array = G.vertex_properties["degree_tortuosity"].get_array()
	701	+ av_array = np.sum(array)/G.num_vertices()
	702	+ G.graph_properties["local_average_degree_tortuosity"] = G.new_graph_property("double", val = av_array)
	703	+ f.write("%.15f\t" % av_array)
	704	+
	705	+ array = G.vertex_properties["bc"].get_array()
	706	+ av_array = np.sum(array)/G.num_vertices()
	707	+ G.graph_properties["local_vertex_bc"] = G.new_graph_property("double", val = av_array)
	708	+ f.write("%.15f\t" % av_array)
	709	+
	710	+ array = G.vertex_properties["degree"].get_array()
	711	+ av_array = np.sum(array)/aabb.getVolume()
	712	+ G.graph_properties["local_vascular_degree"] = G.new_graph_property("double", val = av_array)
	713	+ f.write("%.15f\t" % av_array)
	714	+
	715	+ array = G.vertex_properties["degree_volume"].get_array()
	716	+ av_array = np.sum(array)/aabb.getVolume()
	717	+ G.graph_properties["local_vascular_degree_volume"] = G.new_graph_property("double", val = av_array)
	718	+ f.write("%.15f\t" % av_array)
	719	+
	720	+ array = G.vertex_properties["degree_tortuosity"].get_array()
	721	+ av_array = np.sum(array)/aabb.getVolume()
	722	+ G.graph_properties["local_vascular_degree_tortuosity"] = G.new_graph_property("double", val = av_array)
	723	+ f.write("%.15f\t" % av_array)
	724	+
	725	+ array = G.vertex_properties["bc"].get_array()
	726	+ av_array = np.sum(array)/aabb.getVolume()
	727	+ G.graph_properties["local_vascular_vertex_bc"] = G.new_graph_property("double", val = av_array)
	728	+ f.write("%.15f\t" % av_array)
	729	+
	730	+ array = G.edge_properties["bc"].get_array()
	731	+ av_array = np.sum(array)/aabb.getVolume()
	732	+ G.graph_properties["local_vascular_edge_bc"] = G.new_graph_property("double", val = av_array)
	733	+ f.write("%.15f\t" % av_array)
	734	+
	735	+ G.edge_properties["mst"] = gt.graph_tool.topology.min_spanning_tree(G)
	736	+ 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()))
	737	+ f.write("%.15f\t" % G.graph_properties["local_mst_ratio"])
	738	+
	739	+ G.edge_properties["mst"] = gt.graph_tool.topology.min_spanning_tree(G, weights=G.edge_properties["volume"])
	740	+ 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())
	741	+ f.write("%.15f\t" % G.graph_properties["local_mst_ratio_vascular_volume"])
	742	+
	743	+ [lE, lV] = gt.graph_tool.centrality.eigenvector(G)
	744	+ G.graph_properties["Largest Eigenvector"] = G.new_graph_property("double", val = lE)
	745	+ f.write("%.15f\t" % lE)
	746	+ ############################# hybrid metrics ##########################
	747	+ self.add_rst_metric(G, is_dual=is_dual)
	748	+ f.write("%.15f\t" % G.graph_properties["sensitivity"])
	749	+ f.write("%s" % label)
	750	+ f.write("\n")
	751	+
	752	+ f.close()
	753	+ if(writeKey):
	754	+ location = path + "/key.txt"
	755	+ f = open(location, "w+")
	756	+ f.write("%s\t" % "global_vertex_bc")
	757	+ f.write("%s\t" % "global_vascular_vertex_bc")
	758	+ f.write("%s\t" % "global_vascular_edge_bc")
	759	+ f.write("%s\t" % "global_mst_ratio")
	760	+ f.write("%s\t" % "global_mst_ratio_vascular_volume")
	761	+ f.write("%s\t" % "average_length")
	762	+ f.write("%s\t" % "average_volume")
	763	+ f.write("%s\t" % "average_tortuosity")
	764	+ f.write("%s\t" % "vascular_length")
	765	+ f.write("%s\t" % "vascular_volume")
	766	+ f.write("%s\t" % "vascular_tortuosity")
	767	+ ############################Graph Statistics###########################
	768	+ f.write("%s\t" % "local_average_degree")
	769	+ f.write("%s\t" % "local_average_degree_volume")
	770	+ f.write("%s\t" % "local_average_degree_tortuosity")
	771	+ f.write("%s\t" % "local_vertex_bc")
	772	+ f.write("%s\t" % "local_vascular_degree")
	773	+ f.write("%s\t" % "local_vascular_degree_volume")
	774	+ f.write("%s\t" % "local_vascular_degree_tortuosity")
	775	+ f.write("%s\t" % "local_vascular_vertex_bc")
	776	+ f.write("%s\t" % "local_vascular_edge_bc")
	777	+ f.write("%s\t" % "local_mst_ratio")
	778	+ f.write("%s\t" % "local_mst_ratio_vascular_volume")
	779	+ f.write("%s\t" % "Largest Eigenvector")
	780	+ f.write("%s\t" % "label")
	781	+ f.write("%s\t" % "sensitivity")
	782	+ f.close()
	783	+
	784	+
	785	+ '''
	786	+ Creates a graph from a list of nodes and a list of edges.
	787	+ Uses edge length as weight.
	788	+ Returns a NetworkX Object.
	789	+ '''
	790	+ def createLengthGraph(self):
	791	+ G = nx.Graph()
	792	+ for i in range(len(self.N)):
	793	+ G.add_node(i, p=self.N[i].p)
	794	+ for i in range(len(self.F)):
	795	+ G.add_edge(self.F[i].v0, self.F[i].v1, weight = self.F[i].length())
	796	+ G[self.F[i].v0][self.F[i].v1]['pts'] = self.F[i].points
	797	+ G[self.F[i].v0][self.F[i].v1]['rads'] = self.F[i].radii
	798	+
	799	+ return G
	800	+
	801	+ '''
	802	+ filters all the vertices that are close to the border with deg=1
	803	+ As well as affiliated edges. It is recommended thatthe graph have all the
	804	+ graph metrics refreshed after filtering.
	805	+ '''
	806	+
	807	+ def filterBorder(self, G, is_dual=False):
	808	+ bb = AABB(G, is_dual)
	809	+ TFv = G.new_vertex_property("bool", vals = np.full((G.num_vertices(),1), True, dtype=bool))
	810	+ TFe = G.new_edge_property("bool", vals = np.full((G.num_edges(),1), True, dtype=bool))
	811	+ if(is_dual==False):
	812	+ for v in G.vertices():
	813	+ pt = np.array([G.vertex_properties["p"][v][0], G.vertex_properties["p"][v][1], G.vertex_properties["p"][v][2]])
	814	+ if G.vertex_properties["degree"][v] == 1 and bb.distance(pt) < 5.0:
	815	+ TFv[v] = False
	816	+ for e in v.all_edges():
	817	+ TFe[G.edge(e.source(), e.target())] = False
	818	+ else:
	819	+ #print("I have entered this method")
	820	+ for v in G.vertices():
	821	+ l = len(G.vertex_properties["x"][v])
	822	+ pt_0 = np.array([G.vertex_properties["x"][v][0], G.vertex_properties["y"][v][0], G.vertex_properties["z"][v][0]])
	823	+ 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]])
	824	+ if (G.vertex_properties["degree"][v] == 2):
	825	+ if ((bb.distance(pt_0) < 5.0) or (bb.distance(pt_1) < 5.0)):
	826	+ #print("length", l, ": ", G.vertex_properties["x"][v])
	827	+ #print("degree", G.vertex_properties["degree"][v])
	828	+ TFv[v] = False
	829	+ for e in v.all_edges():
	830	+ TFe[G.edge(e.source(), e.target())] = False
	831	+
	832	+ G.set_filters(TFe, TFv)
	833	+ G1 = gt.Graph(G, directed=False, prune=True)
	834	+ G.clear_filters()
	835	+ G1.clear_filters()
	836	+
	837	+ return G1
	838	+
	839	+
	840	+ '''
	841	+ Returns a Graph object that retains all the properties of the original
	842	+ if return_largest is set as false, then all graphs with more than 20 nodes
	843	+ are returned in a list. Metrics might have to be recalculated after
	844	+ It is recommended that the graph have all the metrics necessary prior to filtering
	845	+ In order to avoid internal boost index errors.
	846	+ '''
	847	+ def filterDisconnected(self, G, return_largest=True):
	848	+ #create masks
	849	+ clusters = G.new_vertex_property("int", vals=np.full((G.num_vertices(), 1), -1, dtype=int))
	850	+ eclusters = G.new_edge_property("int", vals=np.full((G.num_edges(), 1), -1, dtype=int))
	851	+ c = 0
	852	+
	853	+ #run a bfs that sets visited nodes on each vertex
	854	+ for i in G.vertices():
	855	+ if clusters[i] == -1:
	856	+ gt.bfs_search(G, i, VisitorClassDisconnected(eclusters, clusters, c))
	857	+ c = c + 1
	858	+
	859	+ #find the number of clusters
	860	+ unique, counts = np.unique(clusters.get_array(), return_counts=True)
	861	+ eunique, ecounts = np.unique(clusters.get_array(), return_counts=True)
	862	+
	863	+ if(return_largest == True):
	864	+ #find the argument of the largest
	865	+ a = counts.argmax()
	866	+ b = ecounts.argmax()
	867	+
	868	+ #turn the largest cluster into a TF graph mask for the vertices and the edges
	869	+ TFv = G.new_vertex_property("bool", vals = np.full((G.num_vertices(), 1), False, dtype=bool))
	870	+ TFe = G.new_edge_property("bool", vals = np.full((G.num_edges(),1), False, dtype=bool))
	871	+ for i in G.vertices():
	872	+ if clusters[i] == unique[a]:
	873	+ TFv[i] = True
	874	+
	875	+ e = G.get_edges();
	876	+ for i in range(len(e)):
	877	+ if eclusters[G.edge(G.vertex(e[i][0]), G.vertex(e[i][1]))] == eunique[b]:
	878	+ TFe[G.edge(G.vertex(e[i][0]), G.vertex(e[i][1]))] = True
	879	+ #set the filters and return the vertices
	880	+
	881	+ G.set_filters(TFe, TFv)
	882	+ G1=gt.Graph(G, directed=False, prune=True)
	883	+ G.clear_filters()
	884	+ G1.clear_filters()
	885	+ return G1
	886	+ else:
	887	+ Graphs = []
	888	+ #repeat the process above for each unique disconncted component.
	889	+ for j in range(len(unique)):
	890	+ if(counts[j] > 20):
	891	+ TFv = G.new_vertex_property("bool", vals = np.full((G.num_vertices(),1), False, dtype=bool))
	892	+ TFe = G.new_edge_property("bool", vals = np.full((G.num_edges(),1), False, dtype=bool))
	893	+ for i in range(G.num_vertices()):
	894	+ if clusters[G.vertex(i)] == unique[j]:
	895	+ TFv[G.vertex(i)] = True
	896	+
	897	+ e = G.get_edges();
	898	+ for i in range(len(e)):
	899	+ if eclusters[G.edge(G.vertex(e[i][0]), G.vertex(e[i][1]))] == eunique[j]:
	900	+ TFe[G.edge(G.vertex(e[i][0]), G.vertex(e[i][1]))] = True
	901	+ G.set_filters(TFe, TFv)
	902	+ G1=gt.Graph(G, directed=False, prune=True)
	903	+ G.clear_filters()
	904	+ G1.clear_filters()
	905	+ Graphs.append(G1)
	906	+ return Graphs
	907	+
	908	+
	909	+ def simulate_fractures(self, G, remove=10):
	910	+ num_removed = int(np.floor(G.num_edges()*remove/100))
	911	+ print("num of edges to begin with = ", G.num_edges())
	912	+ indices = np.random.randint(0, int(G.num_edges()), size = [num_removed,1])
	913	+ #aabb = AABB(G, is_dual)
	914	+ tf = np.full((G.num_edges(), 1), True, dtype=bool)
	915	+ for j in range(indices.shape[0]):
	916	+ tf[indices[j]] = False
	917	+ TFe = G.new_edge_property("bool", vals = tf)
	918	+ G.set_edge_filter(TFe)
	919	+ G1 = gt.Graph(G, directed=False, prune=True)
	920	+ G.clear_filters()
	921	+ G1.clear_filters()
	922	+ G1 = self.filterDisconnected(G1)
	923	+ G1.vertex_properties["degree"] = G1.degree_property_map("total")
	924	+ print("num of edges left = ", G1.num_edges())
	925	+ print("I should have removed, ", num_removed)
	926	+ print("Instead I removed, ", G.num_edges() - G1.num_edges())
	927	+
	928	+ return G1
	929	+
	930	+ '''
	931	+ Adds the mst based- sensitivity metrics. N is the number of random removals
	932	+ remove is the percentage*100 of the vessels removed in each iteration.
	933	+ The number of vessels removed will be floor(num_edges*remove/100)
	934	+ '''
	935	+ def add_rst_metric(self, G, N=100, remove=10, is_dual=False):
	936	+ #rst_ratio = G.new_graph_property("double")
	937	+ rst_r = 0
	938	+ #G.nu
	939	+ num_removed = int(np.floor(G.num_edges()*remove/100))
	940	+ indices = np.random.randint(0, int(G.num_edges()), size = [N, num_removed])
	941	+ aabb = AABB(G, is_dual)
	942	+
	943	+ for i in range(N):
	944	+ tf = np.full((G.num_edges(),1), True, dtype=bool)
	945	+ for j in range(num_removed):
	946	+ tf[indices[i, j]] = False
	947	+# rst = gt.graph_tool.topology.min_spanning_tree(G, weights = G.)
	948	+# ratio = np.double(np.sum(rst.get_array()))/np.double(G.num_edges())
	949	+# rst_dist.append(ratio)
	950	+# rst_r = rst_r + ratio
	951	+ TFe = G.new_edge_property("bool", vals = tf)
	952	+ G.set_edge_filter(TFe)
	953	+ G1 = gt.Graph(G, directed=False, prune=True)
	954	+ G.clear_filters()
	955	+ G1.clear_filters()
	956	+ while(np.any(G1.degree_property_map("total").get_array() == True)):
	957	+ #print(np.any(G1.degree_property_map("total").get_array()), " ", i)
	958	+ G1 = self.filterDisconnected(G1)
	959	+ G1 = self.filterBorder(G1, is_dual)
	960	+ G1 = self.recalculate_metrics(G1, is_dual)
	961	+
	962	+ if(not is_dual and G1.num_edges() > 2):
	963	+ G1.edge_properties["mst"] = gt.graph_tool.topology.min_spanning_tree(G1, weights=G1.edge_properties["volume"])
	964	+ tvalues = G1.edge_properties["mst"].get_array()
	965	+ #print(tvalues)
	966	+ values = G1.edge_properties["inverted_volume"].get_array()
	967	+ #print(values)
	968	+ for k in range(G1.num_edges()):
	969	+ if(tvalues[k] == True):
	970	+ rst_r = rst_r + values[k]
	971	+ elif(is_dual and G1.num_edges() > 2):
	972	+ #This is non-sensical for now.
	973	+ G1.edge_properties["mst"] = gt.graph_tool.topology.min_spanning_tree(G1)
	974	+ tvalues = G1.vertex_properties["mst"].get_array()
	975	+ values = G1.edge_properties["inverted_volume"].get_array()
	976	+ for k in range(G1.num_edges()):
	977	+ if(tvalues[k] == True):
	978	+ rst_r = rst_r + values[k]
	979	+
	980	+ G.graph_properties["sensitivity"] = G.new_graph_property("double", rst_r/N)
	981	+
	982	+ return G
	983	+
	984	+ '''
	985	+ Re-calculates connectivity based metrics and returns the resulting graph
	986	+ '''
	987	+ def recalculate_metrics(self, G, is_dual=False, normalize=False):
	988	+ gt.graph_tool.centrality.betweenness(G, vprop=G.vertex_properties["bc"], eprop=G.edge_properties["bc"], norm=normalize)
	989	+ if(is_dual == False):
	990	+ #gt.graph_tool.centrality.betweenness(G, G.vertex_properties["bc"], G.edge_properties["bc"], norm=True, is_dual=False)
	991	+ #generate minimum spanning tree
	992	+
	993	+ G.edge_properties["mst"] = gt.graph_tool.topology.min_spanning_tree(G, weights=G.edge_properties["length"])
	994	+ G.graph_properties["mst_ratio"] = np.double(np.sum(G.edge_properties["mst"].get_array()))/np.double(G.num_edges())
	995	+ print(np.double(np.sum(G.edge_properties["mst"].get_array()))/np.double(G.num_edges()))
	996	+ G.vertex_properties["degree"] = G.degree_property_map("total")
	997	+ G.vertex_properties["degree_volume"] = G.degree_property_map("total", weight=G.edge_properties["volume"])
	998	+ G.vertex_properties["degree_tortuosity"] = G.degree_property_map("total", G.edge_properties["tortuosity"])
	999	+
	1000	+ else:
	1001	+ #gt.graph_tool.centrality.betweenness(G, G.vertex_properties["bc"], G.edge_properties["bc"], norm=True, is_dual=False)
	1002	+ 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.
	1003	+ G.graph_properties["mst_ratio"] = np.double(np.sum(G.edge_properties["mst"].get_array()))/np.double(G.num_edges())
	1004	+ G.vertex_properties["degree"] = G.degree_property_map("total")
	1005	+ G.vertex_properties["degree_centrality"] = G.degree_property_map("total", weight=G.edge_properties["bc"])
	1006	+
	1007	+ return G
	1008	+
	1009	+ '''
	1010	+ Filter the graph to remove the disconnected components of the graph if
	1011	+ disconnected is set to true. if borders are set to True, the algorithm
	1012	+ cleans up the nodes with degree one that are close to the edge of the volume.
	1013	+ Returns a Graph object with updated internal property maps.
	1014	+
	1015	+ if return_largerst == False, will return an array of graphs.
	1016	+ disconnected == apply the disconnected componenet filter
	1017	+ return_largest == if false returns the largest set of graphs (n > 10)
	1018	+ borders == True, if true applies the border (deg == 1 near border) filter
	1019	+ add_rst == True, if true adds the rst_metric
	1020	+ '''
	1021	+ def filterFullGraph_gt(self, G, disconnected=True, return_largest=True, borders=True, erode = False, add_rst=False, dual=False):
	1022	+ if(disconnected==True):
	1023	+ if(return_largest==True and borders==True):
	1024	+ G1 = self.filterDisconnected(G, return_largest)
	1025	+ G1 = self.recalculate_metrics(G1, is_dual=dual)
	1026	+ G1 = self.filterBorder(G1, is_dual=dual)
	1027	+ if(erode == True):
	1028	+ while(np.any(G1.degree_property_map("total").get_array() == True)):
	1029	+ G1 = self.filterBorder(G1, is_dual=dual)
	1030	+ G1 = self.recalculate_metrics(G1, is_dual=dual)
	1031	+ else:
	1032	+ G1 = self.recalculate_metrics(G1, is_dual=dual)
	1033	+ if(add_rst == True):
	1034	+ G1 = self.add_rst_metric(G1)
	1035	+ elif(return_largest==True and borders==False):
	1036	+ G1 = self.filterDisconnected(G, return_largest)
	1037	+ if(erode == True):
	1038	+ while(np.any(G1.degree_property_map("total").get_array() == True)):
	1039	+ G1 = self.filterBorder(G1, is_dual=dual)
	1040	+ G1 = self.recalculate_metrics(G1, is_dual=dual)
	1041	+ else:
	1042	+ G1 = self.recalculate_metrics(G1, is_dual=dual)
	1043	+ if(add_rst == True):
	1044	+ G1 = self.add_rst_metric(G1)
	1045	+ elif(return_largest==False and borders == True):
	1046	+ G1 = self.filterDisconnected(G, return_largest)
	1047	+ for graph in G1:
	1048	+ graph = self.recalculate_metrics(graph, is_dual=dual)
	1049	+ graph = self.filterBorder(graph, is_dual=dual)
	1050	+ if(erode == True):
	1051	+ while(np.any(graph.degree_property_map("total").get_array() == True)):
	1052	+ graph = self.filterBorder(graph, is_dual=dual)
	1053	+ graph = self.recalculate_metrics(graph, is_dual=dual)
	1054	+ else:
	1055	+ graph = self.recalculate_metrics(graph, is_dual=dual)
	1056	+ if(add_rst == True):
	1057	+ graph = self.add_rst_metric(graph)
	1058	+ else:
	1059	+ G1 = self.filterDisconnected(G, return_largest)
	1060	+ for graph in G1:
	1061	+ graph = self.recalculate_metrics(graph, is_dual=dual)
	1062	+ if(add_rst == True):
	1063	+ graph = self.add_rst_metric(graph)
	1064	+ else:
	1065	+ G1 = self.filterBorder(G, is_dual=dual);
	1066	+ if(erode == True):
	1067	+ while(np.any(G1.degree_property_map("total").get_array() == True)):
	1068	+ print(G1.num_vertices())
	1069	+ G1 = self.filterBorder(G1, is_dual=dual)
	1070	+ print(G1.num_vertices())
	1071	+ G1 = self.recalculate_metrics(G1, is_dual=dual, )
	1072	+ else:
	1073	+ G1 = self.recalculate_metrics(G1, is_dual=dual)
	1074	+ if(add_rst==True):
	1075	+ G1 = self.add_rst_metric(G1)
	1076	+
	1077	+ return G1
	1078	+
	1079	+ '''
	1080	+ Accretes the graph in G using the nodes and edges in G_0
	1081	+ Assumes that G_0 has all the same properties as G, including "idx"
	1082	+ idx is an immutable set of indices that do not change when the graph is
	1083	+ modified
	1084	+ '''
	1085	+ def accrete(self, G, G_0):
	1086	+
	1087	+ v_filters = G_0.new_vertex_property("bool", vals = np.full((G_0.num_vertices(),1), False, dtype=bool))
	1088	+ e_filters = G_0.new_edge_property("bool", vals = np.full((G_0.num_edges(),1), False, dtype=bool))
	1089	+ for v in G.vertices():
	1090	+ #print(v)
	1091	+ G_0_vertex = G_0.vertex(G.vertex_properties["idx"][v])
	1092	+ v_filters[G_0_vertex] = True
	1093	+ for e in G_0_vertex.all_edges():
	1094	+ e_filters[e] = True
	1095	+ for v_0 in G_0_vertex.all_neighbors():
	1096	+ v_filters[v_0] = True
	1097	+
	1098	+ G_0.set_filters(e_filters, v_filters)
	1099	+
	1100	+ graph = gt.Graph(G_0, directed=False, prune=True)
	1101	+ G_0.clear_filters()
	1102	+
	1103	+ return graph
	1104	+
	1105	+
	1106	+
	1107	+
	1108	+ '''
	1109	+ Creates a graph from a list of nodes and a list of edges
	1110	+ The graph is the "inverse" of the original graph,
	1111	+ Vertices are edges and edges are
	1112	+
	1113	+ '''
	1114	+ def createDualGraph_gt(self):
	1115	+ #Generate the undirected graph
	1116	+ G = gt.Graph()
	1117	+ G.set_directed(False)
	1118	+
	1119	+ #add all the required vertex properties
	1120	+ vpos = G.new_vertex_property("vector<double>") #original position of the edge (placed as midpoin of edge)
	1121	+ pos = G.new_vertex_property("vector<double>") #positions according to the fdl
	1122	+ rpos = G.new_vertex_property("vector<double>") #positions according to the radial layout
	1123	+ x_edge = G.new_vertex_property("vector<double>") #x-coordinated of the edge
	1124	+ y_edge = G.new_vertex_property("vector<double>") #y-coordinates of the edge
	1125	+ z_edge = G.new_vertex_property("vector<double>") #z-coordinates of the edge
	1126	+ r_edge = G.new_vertex_property("vector<double>") #r-values at each x,y,z
	1127	+ l_edge = G.new_vertex_property("double") #length of the edge
	1128	+ i_edge = G.new_vertex_property("double") #1/length
	1129	+ iv_edge = G.new_vertex_property("double") #1/volume
	1130	+ t_edge = G.new_vertex_property("double") #tortuosity of the edge
	1131	+ v_edge = G.new_vertex_property("double") #volume of the edge
	1132	+ vbetweeness_centrality = G.new_vertex_property("double") #betweeneness centrality of the vertex
	1133	+ gaussian = G.new_vertex_property("double") #empty property map for gaussian values downstream
	1134	+ degree = G.new_vertex_property("int")
	1135	+ degree_centrality = G.new_vertex_property("int")
	1136	+
	1137	+ #add all the required edge properties
	1138	+ ebetweeness_centrality = G.new_edge_property("double") #betweeness centrality of the edge
	1139	+ mst = G.new_edge_property("bool") #Boolean value deligating belongance to the minimum spanning tree.
	1140	+ l_edge = G.new_edge_property("double") #length of the edge
	1141	+ gaussian = G.new_edge_property("double")
	1142	+
	1143	+
	1144	+ #add all graph properties
	1145	+ mst_ratio = G.new_graph_property("double")
	1146	+
	1147	+ #add all the edges as vertices
	1148	+ for i in range(len(self.F)):
	1149	+ if(self.F[i].length() > 0):
	1150	+ G.add_vertex()
	1151	+ x,y,z,r = self.F[i].getcoords()
	1152	+ x_edge[G.vertex(i)] = x
	1153	+ y_edge[G.vertex(i)] = y
	1154	+ z_edge[G.vertex(i)] = z
	1155	+ r_edge[G.vertex(i)] = r
	1156	+ l_edge[G.vertex(i)] = self.F[i].length()
	1157	+ gaussian[G.vertex(i)] = np.exp(-np.power(l_edge[G.vertex(i)], 2.) / (2 * np.power(60.0, 2.)))
	1158	+ i_edge[G.vertex(i)] = 1/l_edge[G.vertex(i)]
	1159	+ t_edge[G.vertex(i)] = self.F[i].turtuosity()
	1160	+ v_edge[G.vertex(i)] = self.F[i].volume()
	1161	+ iv_edge[G.vertex(i)] = 1/v_edge[G.vertex(i)]
	1162	+ vpos[G.vertex(i)] = np.array([x[int(len(x)/2)], y[int(len(x)/2)], z[int(len(x)/2)]], dtype=float)
	1163	+
	1164	+ #based on the neighbors add edges. If two edges share a vertex, there is a vertex between them.
	1165	+ for i in range(len(self.N)):
	1166	+ #in case there are 3 or more edges attached to a vertex
	1167	+ if(len(self.N[i].o) > 2):
	1168	+ for j in range(len(self.N[i].o)-1):
	1169	+ G.add_edge(G.vertex(self.N[i].o[j]), G.vertex(self.N[i].o[j+1]))
	1170	+ l_edge[G.edge(G.vertex(self.N[i].o[j])), G.vertex(self.N[i].o[j+1])] = 1
	1171	+ gaussian[G.edge(G.vertex(self.N[i].o[j])), G.vertex(self.N[i].o[j+1])] = \
	1172	+ 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.)))
	1173	+ #add the final edge.
	1174	+ G.add_edge(G.vertex(self.N[i].o[0]), G.vertex(self.N[i].o[len(self.N[i].o)-1]))
	1175	+ l_edge[G.edge(G.vertex(self.N[i].o[j])), G.vertex(self.N[i].o[j+1])] = 1
	1176	+ gaussian[G.edge(G.vertex(self.N[i].o[j])), G.vertex(self.N[i].o[j+1])] = \
	1177	+ 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.)))
	1178	+ #in case there are 2 edges attached to a vertex
	1179	+ elif(len(self.N[i].o) == 2):
	1180	+ G.add_edge(G.vertex(self.N[i].o[0]), G.vertex(self.N[i].o[1]-1))
	1181	+ l_edge[G.edge(G.vertex(self.N[i].o[j])), G.vertex(self.N[i].o[j+1])] = 1
	1182	+ gaussian[G.edge(G.vertex(self.N[i].o[j])), G.vertex(self.N[i].o[j+1])] = \
	1183	+ 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.)))
	1184	+ #in case there is only a single edge coming into the vertex we do
	1185	+ # not add anything
	1186	+
	1187	+ #generate centrality map
	1188	+ gt.graph_tool.centrality.betweenness(G, vbetweeness_centrality, ebetweeness_centrality, norm=True)
	1189	+ #generate minimum spanning tree
	1190	+ mst = gt.graph_tool.topology.min_spanning_tree(G, weights=ebetweeness_centrality)
	1191	+ mst_ratio[G] = np.double(np.sum(mst.get_array()))/np.double(len(self.N))
	1192	+ degree = G.degree_property_map("total")
	1193	+ degree_centrality = G.degree_property_map("total", weight=ebetweeness_centrality)
	1194	+ #degree_tortuosity = G.degree_property_map("total", weight=t_edge)
	1195	+
	1196	+ pos = gt.sfdp_layout(G)
	1197	+
	1198	+ rpos = gt.radial_tree_layout(G, root=G.vertex(np.argmax(vbetweeness_centrality.get_array())))
	1199	+
	1200	+ #set property maps for the vertices
	1201	+ G.vertex_properties["p"] = vpos
	1202	+ G.vertex_properties["pos"] = pos
	1203	+ G.vertex_properties["rpos"] = rpos
	1204	+ G.vertex_properties["bc"] = vbetweeness_centrality
	1205	+ G.vertex_properties["degree_centrality"] = degree_centrality
	1206	+ G.vertex_properties["degree"] = degree
	1207	+ G.vertex_properties["x"] = x_edge
	1208	+ G.vertex_properties["y"] = y_edge
	1209	+ G.vertex_properties["z"] = z_edge
	1210	+ G.vertex_properties["r"] = r_edge
	1211	+ G.vertex_properties["inverted"] = i_edge
	1212	+ G.vertex_properties["inverted_volume"] = iv_edge
	1213	+ G.vertex_properties["length"] = l_edge
	1214	+ G.vertex_properties["tortuosity"] = t_edge
	1215	+ G.vertex_properties["volume"] = v_edge
	1216	+ G.vertex_properties["bc"] = vbetweeness_centrality
	1217	+
	1218	+ #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)
	1219	+ G.vertex_properties["pos"] = gt.sfdp_layout(G, C = 0.5, p = 3.0)
	1220	+ #set property maps for the edges
	1221	+ G.edge_properties["bc"] = ebetweeness_centrality
	1222	+ G.edge_properties["mst"] = mst
	1223	+
	1224	+ #set graph properies
	1225	+ G.graph_properties["mst_ratio"] = mst_ratio
	1226	+
	1227	+
	1228	+
	1229	+# title = "raw_dual_cortical_7_6_before_delete.pdf"
	1230	+# title2 = "raw_dual_radial_cordical_7_6_before_delete.pdf"
	1231	+#
	1232	+# 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)
	1233	+# 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)
	1234	+#
	1235	+ title = "./Original_2D_Graph.png"
	1236	+# title2 = "raw_dual_radial_cordical_7_6_after_delete.pdf"
	1237	+#
	1238	+ G1 = self.filterFullGraph_gt(G, borders=False, dual=True)
	1239	+
	1240	+# #print(clusters.get_array()[:], clusters.get_array()[:])
	1241	+# #pos = gt.sfdp_layout(G)
	1242	+#
	1243	+# #gt.graph_draw(G, pos=,edge_color=mst, vertex_fill_color=vbetweeness_centrality, edge_pen_width=ebetweeness_centrality, output="raw.pdf")
	1244	+#
	1245	+# 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)
	1246	+# 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)
	1247	+ 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)
	1248	+ title = "./Original_2D_Layout.png"
	1249	+ 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)
	1250	+#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)
	1251	+#
	1252	+#
	1253	+# G2 = self.filterFullGraph_gt(G1, add_rst=False, dual=True)
	1254	+# title = "raw_dual_cortical_7_6_after_borders.pdf"
	1255	+# title2 = "raw_dual_radial_cortical_7_6_after_borders.pdf"
	1256	+# 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)
	1257	+# 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)
	1258	+ return G
	1259	+
	1260	+
	1261	+ '''
	1262	+ Create a graph from a list of nodes and a list of edges.
	1263	+ Populates the Graph with a bunch of property maps in respect to edges and nodes
	1264	+ each vertex is assigned 3 positional coordiantes:
	1265	+ The original x,y,z location of the graph.
	1266	+ x,y locations according to force directed layout,
	1267	+ and x,y positions according to a radial layout
	1268	+ Returns a graph_tool object
	1269	+ '''
	1270	+ def createFullGraph_gt(self):
	1271	+ #Generate the undirected graph
	1272	+ G = gt.Graph()
	1273	+ G.set_directed(False)
	1274	+
	1275	+ #add all the required vertex properties
	1276	+ vpos = G.new_vertex_property("vector<double>") #original vertex positions
	1277	+ pos = G.new_vertex_property("vector<double>") #positions according to the fdl
	1278	+ rpos = G.new_vertex_property("vector<double>") #positions according to the radial layout
	1279	+ vbetweeness_centrality = G.new_vertex_property("double") #betweeneness centrality of the vertex
	1280	+ degree = G.new_vertex_property("int") #degree
	1281	+ degree_volume = G.new_vertex_property("double") #degree scaled by the volume of all in fibers
	1282	+ degree_tortuosity = G.new_vertex_property("double") #degree scaled by the tortuosity of all in fibers.
	1283	+
	1284	+
	1285	+ #add all the required edge properties
	1286	+ #G.properties[("x,y,z,r")] = G.new_edge_property("vector<double")
	1287	+ x_edge = G.new_edge_property("vector<double>") #x-coordinated of the edge
	1288	+ y_edge = G.new_edge_property("vector<double>") #y-coordinates of the edge
	1289	+ z_edge = G.new_edge_property("vector<double>") #z-coordinates of the edge
	1290	+ r_edge = G.new_edge_property("vector<double>") #r-values at each x,y,z
	1291	+ av_edge = G.new_edge_property("double") #average edge radius
	1292	+ l_edge = G.new_edge_property("double") #length of the edge
	1293	+ i_edge = G.new_edge_property("double") #1/length
	1294	+ iv_edge = G.new_edge_property("double") #1/volume
	1295	+ t_edge = G.new_edge_property("double") #tortuosity of the edge
	1296	+ v_edge = G.new_edge_property("double") #volume of the edge
	1297	+ ebetweeness_centrality = G.new_edge_property("double") #betweeness centrality of the edge
	1298	+ ebc_length = G.new_edge_property("double") #betweeneness centrality scaled by the length
	1299	+ mst = G.new_edge_property("bool") #Boolean value deligating belongance to the minimum spanning tree.
	1300	+
	1301	+
	1302	+ #This map gets updated.
	1303	+ gaussian = G.new_edge_property("double") #empty property map for gaussian values downstream
	1304	+
	1305	+
	1306	+ #Graph properties (once per region)
	1307	+ mst_ratio = G.new_graph_property("double") #Graph based metric of mst edges/total edges.
	1308	+
	1309	+
	1310	+ #add verticies and set their properties
	1311	+ for i in range(len(self.N)):
	1312	+ G.add_vertex()
	1313	+ vpos[G.vertex(i)] = self.N[i].p
	1314	+
	1315	+ #add edges and set their properties.
	1316	+ for i in range(len(self.F)):
	1317	+ v0 = self.F[i].v0
	1318	+ v1 = self.F[i].v1
	1319	+ if self.F[i].length() > 0.0:
	1320	+ x,y,z,r = self.F[i].getcoords()
	1321	+ G.add_edge(G.vertex(v0), G.vertex(v1))
	1322	+ l_edge[G.edge(v0,v1)] = self.F[i].length()
	1323	+ gaussian[G.edge(v0,v1)] = np.exp(-np.power(l_edge[G.edge(v0,v1)], 2.) / (2 * np.power(60.0, 2.)))
	1324	+ i_edge[G.edge(v0,v1)] = 1/l_edge[G.edge(v0,v1)]
	1325	+ t_edge[G.edge(v0,v1)] = self.F[i].turtuosity()
	1326	+ v_edge[G.edge(v0,v1)] = self.F[i].volume()
	1327	+ iv_edge[G.edge(v0,v1)] = 1/v_edge[G.edge(v0,v1)]
	1328	+ x_edge[G.edge(v0,v1)] = x
	1329	+ y_edge[G.edge(v0,v1)] = y
	1330	+ z_edge[G.edge(v0,v1)] = z
	1331	+ r_edge[G.edge(v0,v1)] = r
	1332	+ av_edge[G.edge(v0,v1)] = self.F[i].av_radius()
	1333	+ else:
	1334	+ print("WTF")
	1335	+
	1336	+
	1337	+ #generate centrality map
	1338	+ gt.graph_tool.centrality.betweenness(G, vprop=vbetweeness_centrality, eprop=ebetweeness_centrality, norm=True)
	1339	+ #generate minimum spanning tree
	1340	+ mst = gt.graph_tool.topology.min_spanning_tree(G, weights=l_edge)
	1341	+ mst_ratio[G] = np.double(np.sum(mst.get_array()))/np.double(len(self.N))
	1342	+ degree = G.degree_property_map("total")
	1343	+ degree_volume = G.degree_property_map("total", weight=v_edge)
	1344	+
	1345	+ dg = degree_volume.get_array()
	1346	+ dg = np.exp(-np.power(dg, 2.) / (2 * np.power(0.000005, 2.)))
	1347	+ degree_volume = G.new_vertex_property("double", vals=dg)
	1348	+ degree_tortuosity = G.degree_property_map("total", weight=t_edge)
	1349	+
	1350	+
	1351	+ #print(clusters.get_array()[:], clusters.get_array()[:])
	1352	+ pos = gt.sfdp_layout(G, C = 1.0, K = 10)
	1353	+ rpos = gt.radial_tree_layout(G, root=G.vertex(np.argmax(vbetweeness_centrality.get_array())))
	1354	+
	1355	+
	1356	+
	1357	+ for e in G.edges():
	1358	+ ebc_length[G.edge(e.source(), e.target())] = ebetweeness_centrality[G.edge(e.source(), e.target())]*l_edge[G.edge(e.source(), e.target())]
	1359	+
	1360	+ #set property maps for the vertices
	1361	+ G.vertex_properties["p"] = vpos
	1362	+ G.vertex_properties["pos"] = pos
	1363	+ G.vertex_properties["rpos"] = rpos
	1364	+ G.vertex_properties["bc"] = vbetweeness_centrality
	1365	+ G.vertex_properties["degree"] = degree
	1366	+ G.vertex_properties["degree_volume"] = degree_volume
	1367	+ G.vertex_properties["degree_tortuosity"] = degree_tortuosity
	1368	+ G.vertex_properties["selection"] = G.new_vertex_property("bool", val=False)
	1369	+
	1370	+ #set property maps for the edges
	1371	+ G.edge_properties["x"] = x_edge
	1372	+ G.edge_properties["y"] = y_edge
	1373	+ G.edge_properties["z"] = z_edge
	1374	+ G.edge_properties["r"] = r_edge
	1375	+ G.edge_properties["inverted"] = i_edge
	1376	+ G.edge_properties["length"] = l_edge
	1377	+ G.edge_properties["tortuosity"] = t_edge
	1378	+ G.edge_properties["av_radius"] = av_edge
	1379	+ G.edge_properties["volume"] = v_edge
	1380	+ G.edge_properties["bc"] = ebetweeness_centrality
	1381	+ G.edge_properties["bc*length"] = ebc_length
	1382	+ G.edge_properties["mst"] = mst
	1383	+ G.edge_properties["inverted_volume"] = iv_edge
	1384	+ G.edge_properties["gaussian"] = gaussian
	1385	+ G.edge_properties["selection"] = G.new_edge_property("double", val=0.0)
	1386	+
	1387	+ #set graph properies
	1388	+ G.graph_properties["mst_ratio"] = mst_ratio
	1389	+
	1390	+ N = 0
	1391	+ for e in G.edges():
	1392	+ N += len(G.edge_properties["x"][e])
	1393	+ point_cloud_x = np.zeros((N, 1), dtype=np.float64)
	1394	+ point_cloud_y = np.zeros((N, 1), dtype=np.float64)
	1395	+ point_cloud_z = np.zeros((N, 1), dtype=np.float64)
	1396	+
	1397	+ n = 0
	1398	+ for e in G.edges():
	1399	+ for p in range(len(G.edge_properties["x"][e])):
	1400	+ point_cloud_x[n][0] = G.edge_properties["x"][e][p]
	1401	+ point_cloud_y[n][0] = G.edge_properties["y"][e][p]
	1402	+ point_cloud_z[n][0] = G.edge_properties["z"][e][p]
	1403	+ n += 1
	1404	+
	1405	+ G.graph_properties["point_cloud_x"] = G.new_graph_property("vector<double>")
	1406	+ G.graph_properties["point_cloud_y"] = G.new_graph_property("vector<double>")
	1407	+ G.graph_properties["point_cloud_z"] = G.new_graph_property("vector<double>")
	1408	+ G.graph_properties["point_cloud_x"] = point_cloud_x
	1409	+ G.graph_properties["point_cloud_y"] = point_cloud_y
	1410	+ G.graph_properties["point_cloud_z"] = point_cloud_z
	1411	+ #save the original graph indices for the edges and the vertices.
	1412	+ G.vertex_properties["idx"] = G.vertex_index
	1413	+ G.edge_properties["idx"] = G.edge_index
	1414	+
	1415	+ title = "~/Pictures/Original_2D_Graph_2.png"
	1416	+# title2 = "raw_dual_radial_cordical_7_6_after_delete.pdf"
	1417	+#
	1418	+# G1 = self.filterFullGraph_gt(G, borders=False, dual=False, add_rst = False)
	1419	+##
	1420	+## #print(clusters.get_array()[:], clusters.get_array()[:])
	1421	+## #pos = gt.sfdp_layout(G)
	1422	+##
	1423	+## #gt.graph_draw(G, pos=,edge_color=mst, vertex_fill_color=vbetweeness_centrality, edge_pen_width=ebetweeness_centrality, output="raw.pdf")
	1424	+##
	1425	+## 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)
	1426	+## 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)
	1427	+# 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)
	1428	+# #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)
	1429	+# title = "~/Pictures/Original_2D_Layout.png"
	1430	+# 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)
	1431	+#
	1432	+#
	1433	+# G1 = self.filterFullGraph_gt(G, borders=True, dual=False, erode=True, add_rst = False)
	1434	+# title = "~/Pictures/Original_2D_Graph_Eroded.png"
	1435	+# 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)
	1436	+# title = "~/Pictures/Original_2D_Layout_Eroded.png"
	1437	+# 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)
	1438	+
	1439	+# title = "raw_full_cortical_7_6_before_delete.pdf"
	1440	+# title2 = "raw_full_radial_cordical_7_6_before_delete.pdf"
	1441	+# #G.graph_properties["rst_ratio"] = rst_ratio
	1442	+# #print("")
	1443	+# #print(G.graph_properties["rst_ratio"])
	1444	+# #print(G.edge_properties["bc"].get_array())
	1445	+# #print()
	1446	+#
	1447	+# #G.vertex_properties[("p","betweeness_centrality")] = [vpos, vbetweeness_centrality]
	1448	+# #G.edge_properties[("x", "y", "z", "r", "length", "tortuosity", "volume", "betweeness_centrality")] = []
	1449	+# #G.set_edge_filter(mst)
	1450	+# #G.set_edge_filter(mst)
	1451	+#
	1452	+# #Experimental code for drawing graphs.
	1453	+# #gt.graph_draw(G,pos=vpos,edge_color=mst, vertex_fill_color=vbetweeness_centrality, edge_pen_width=ebetweeness_centrality, output="raw.pdf")
	1454	+# #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])
	1455	+# #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)
	1456	+# #print(np.argmax(vbetweeness_centrality.get_array()))
	1457	+#
	1458	+# #np.savetxt("vbetweeness_centrality", vbetweeness_centrality.get_array())
	1459	+# 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)
	1460	+# 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)
	1461	+#
	1462	+#
	1463	+# #gt.graph_draw(G,pos=vpos, efilt=mst, vertex_fill_color=vbetweeness_centrality, edge_pen_width=ebetweeness_centrality, output="mst.pdf")
	1464	+# #gt.graph_draw(G,pos=pos, efilt=mst, vertex_fill_color=vbetweeness_centrality, output="mst.pdf")
	1465	+# #need to create subgraphs!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	1466	+# #rst_r = 0
	1467	+# #for i in range(1000):
	1468	+# # rst = gt.graph_tool.topology.random_spanning_tree(G)
	1469	+# # rst_r = rst_r + np.double(np.sum(rst.get_array()))/np.double(len(self.N))
	1470	+#
	1471	+# #rst_ratio[G] = rst_r/1000.0
	1472	+#
	1473	+#
	1474	+#
	1475	+# #G1 = gt.Graph(G)
	1476	+#
	1477	+# #G1.purge_edges()
	1478	+# #G1.purge_vertices(in_place = True)
	1479	+#
	1480	+#
	1481	+# title = "raw_full_cortical_7_6_after_delete.pdf"
	1482	+# title2 = "raw_full_radial_cordical_7_6_after_delete.pdf"
	1483	+#
	1484	+# G1 = self.filterFullGraph_gt(G,borders=False)
	1485	+#
	1486	+# #print(clusters.get_array()[:], clusters.get_array()[:])
	1487	+# #pos = gt.sfdp_layout(G)
	1488	+#
	1489	+# #gt.graph_draw(G,pos=vpos,edge_color=mst, vertex_fill_color=vbetweeness_centrality, edge_pen_width=ebetweeness_centrality, output="raw.pdf")
	1490	+#
	1491	+# gt.graph_draw(pos = gt.sfdp_layout(G1), output=title)
	1492	+# 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)
	1493	+# #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)
	1494	+#
	1495	+#
	1496	+# G2 = self.filterFullGraph_gt(G1)
	1497	+# title = "raw_full_cortical_7_6_after_borders.pdf"
	1498	+# title2 = "raw_full_radial_cortical_7_6_after_borders.pdf"
	1499	+# 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)
	1500	+# 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)
	1501	+# #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])
	1502	+#
	1503	+# #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)
	1504	+# #print(np.argmax(vbetweeness_centrality.get_array()))
	1505	+# #np.savetxt("vbetweeness_centrality", vbetweeness_centrality.get_array())
	1506	+# #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)
	1507	+# #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"])
	1508	+# #gt.draw_hierarchy(gt.minimize_nested_blockmodel_dl(G1, deg_corr=False), output = "test.pdf")
	1509	+# #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)
	1510	+# print(gt.graph_tool.openmp_enabled(), gt.graph_tool.openmp_get_num_threads())
	1511	+
	1512	+ #print(ratio, G2.graph_properties["mst_ratio"])
	1513	+ #n, bins, patches = plt.hist(np.asarray(rst_dist), 50, normed=1, facecolor='green', alpha=0.75)
	1514	+
	1515	+ # add a 'best fit' line
	1516	+ #plt.grid(True)
	1517	+
	1518	+ #plt.show()
	1519	+ #self.partition(G, 5, 7, False)
	1520	+ self.export_points(G)
	1521	+ self.export_bb(G)
	1522	+ return G
	1523	+
	1524	+
	1525	+ def gen_new_fd_layout(G):
	1526	+ pos = gt.sfdp_layout(G, C = 1.0, K = 10)
	1527	+ G.vertex_properties["pos"] = pos
	1528	+ return G
	1529	+
	1530	+ def map_edges_to_range(G, rng, propertymap):
	1531	+ def func(maximum, minimum, new_maximum, new_minimum, value):
	1532	+ return ((value-minimum)/(maximum-minimum)*(new_maximum-new_minimum)+new_minimum)
	1533	+
	1534	+ mx = max(G.edge_properties[propertymap].get_array())
	1535	+ mn = min(G.edge_properties[propertymap].get_array())
	1536	+ G.edge_properties["map"] = G.new_edge_property("float")
	1537	+ gt.map_property_values(G.edge_properties[propertymap], G.edge_properties["map"], lambda x: func(mx, mn, rng[0], rng[1], x))
	1538	+
	1539	+ return G.edge_properties["map"]
	1540	+
	1541	+
	1542	+ def map_vertices_to_range(G, rng, propertymap):
	1543	+ def func(maximum, minimum, new_maximum, new_minimum, value):
	1544	+ return ((value-minimum)/(maximum-minimum)*(new_maximum-new_minimum)+new_minimum)
	1545	+
	1546	+ mx = max(G.vertex_properties[propertymap].get_array())
	1547	+ mn = min(G.vertex_properties[propertymap].get_array())
	1548	+ G.vertex_properties["map"] = G.new_vertex_property("float")
	1549	+ gt.map_property_values(G.vertex_properties[propertymap], G.vertex_properties["map"], lambda x: func(mx, mn, rng[0], rng[1], x))
	1550	+
	1551	+ return G.vertex_properties["map"]
	1552	+
	1553	+ '''
	1554	+ G and propertymap in G are passed in.
	1555	+ Maps select property from range to color, works for int's floats and vec3
	1556	+ and vec4 values
	1557	+ Returns a vector<double> property map that maps the property value to color.
	1558	+ If passed a vertex propertymap, returns a vertex property map
	1559	+ If passed an edge propertymap, returns an edge property map
	1560	+ '''
	1561	+ def map_property_to_color(G, propertymap, colormap = 'tab20'):
	1562	+ key_type = ''
	1563	+ if propertymap.key_type()=='e':
	1564	+ key_type = 'e'
	1565	+ elif propertymap.key_type()=='v':
	1566	+ key_type = 'v'
	1567	+ else:
	1568	+ #TO DO: Write graph property type return.
	1569	+ print("Graph_property passed")
	1570	+ if G.properties[(key_type, 'RGBA')] == None:
	1571	+ color = [0.0, 0.0, 0.0, 1.0]
	1572	+ G.properties[(key_type, 'RGBA')] = G.new_property(key_type, "vector<double>", val=color)
	1573	+
	1574	+ #This is when the property map is integers
	1575	+ #int32_t when integer
	1576	+ value_type = propertymap.value_type()
	1577	+ print(value_type)
	1578	+ if(value_type == "int32_t"):
	1579	+ array = propertymap.get_array()
	1580	+ colors = cm.get_cmap(colormap, len(np.unique(array)))
	1581	+ if key_type == 'v':
	1582	+ for v in G.vertices():
	1583	+ G.vertex_properties["RGBA"][v] = colors(propertymap[v])
	1584	+ return G.vertex_properties["RGBA"]
	1585	+ elif(value_type == 'double' or value_type == 'float'):
	1586	+ array = propertymap.get_array()
	1587	+ norm = cm.colors.Normalize(vmin=min(array), vmax=max(array))
	1588	+ #colors = cm.get_cmap(colormap, len(np.unique(array)))
	1589	+ colors = cm.ScalarMappable(norm=norm, cmap=colormap)
	1590	+ if key_type =='v':
	1591	+ for v in G.vertices():
	1592	+ print(colors.to_rgba(propertymap[v]))
	1593	+ G.vertex_properties["RGBA"][v] = colors.to_rgba(propertymap[v])
	1594	+ return G.vertex_properties["RGBA"]
	1595	+ elif(value_type == 'vector<double>' or value_type == 'vector<float>'):
	1596	+ print("detected a vector of doubles or floats")
	1597	+
	1598	+
	1599	+
	1600	+ '''
	1601	+ Attempts to partition the graph G into b_i = [0, n^3] sections, where each section is the connected componenets of
	1602	+ m steps away from the source node (where the source node is the closest node to point n_i)
	1603	+ '''
	1604	+ def partition(self, G, n, m, is_dual, path, prefix, saveKey, lbl = "none"):
	1605	+ bb = AABB(G, is_dual)
	1606	+ #print("FLJKKHDFLKJFDLKJFDLKJ ", m)
	1607	+ x, y, z = bb.project_grid(n)
	1608	+ #cluster_belongance = G.new_vertex_property("vector<boolean>")
	1609	+ #ecluster_belongance = G.new_edge_property("vector<boolean>")
	1610	+ #X, Y, Z = np.meshgrid(x,y,z)
	1611	+ array_cluster_bel = np.zeros([n**3, G.num_vertices()])
	1612	+ array_ecluster_bel = np.zeros([n**3, G.num_edges()])
	1613	+ c = 0
	1614	+ for i in range(len(x)):
	1615	+ for j in range(len(y)):
	1616	+ for k in range(len(z)):
	1617	+ #for i in range(1):
	1618	+ # for j in range(1):
	1619	+ # for k in range(1):
	1620	+ p = G.vertex_properties["p"].get_2d_array(range(G.num_vertices()))
	1621	+ if(p.shape[1] < 3):
	1622	+ break
	1623	+ #print("stuff", p)
	1624	+ idx = -1
	1625	+ dist = 100000000000
	1626	+ for M in range(p.shape[1]):
	1627	+ d = math.sqrt((x[i] - p[0][M])2 + (y[j] - p[1][M])2 + (z[k] - p[2][M])**2)
	1628	+ if d < dist:
	1629	+ idx = M
	1630	+ dist = d
	1631	+ print(idx, " vertex[",p[0][idx], p[1][idx], p[2][idx], "], point [", x[i], y[j], z[k], "]")
	1632	+ clusters = G.new_vertex_property("int", vals=np.full((G.num_vertices(), 1), -1, dtype=int))
	1633	+ #eclusters = G.new_edge_property("int", vals=np.full((G.num_edges(), 1), -1, dtype=int))
	1634	+ dist = G.new_vertex_property("int", val = 100000)
	1635	+ dist[G.vertex(idx)] = 0
	1636	+ #run a bfs on the index node
	1637	+ gt.bfs_search(G, idx, VisitorClassPartition(clusters, dist, c, m))
	1638	+ #print(dist.get_array());
	1639	+ temp_vertices = np.ma.masked_where(dist.get_array() <= m, dist.get_array()).mask
	1640	+ #print(temp_vertices)
	1641	+ temp_edges = np.zeros([G.num_edges()])
	1642	+ for vertex in np.nditer(np.where(temp_vertices == True)):
	1643	+ for edge in G.get_out_edges(vertex):
	1644	+ #print(edge)
	1645	+ temp_edges[edge[2]] = 1
	1646	+ array_cluster_bel[c,:] = temp_vertices[:]
	1647	+ array_ecluster_bel[c,:] = temp_edges[:]
	1648	+ c = c + 1
	1649	+ #print(G.num_vertices())
	1650	+ #np.savetxt("clusters.txt", array_cluster_bel)
	1651	+ G.vertex_properties["partition"] = G.new_vertex_property("vector<boolean>", array_cluster_bel)
	1652	+ G.edge_properties["partition"] = G.new_edge_property("vector<boolean>", array_ecluster_bel)
	1653	+ j = 0
	1654	+ gt.graph_draw(G, pos = gt.sfdp_layout(G), output="Graph_full.pdf")
	1655	+ for i in range(c):
	1656	+ TFv = G.new_vertex_property("bool", vals=array_cluster_bel[i,:])
	1657	+ TFe = G.new_edge_property("bool", vals=array_ecluster_bel[i,:])
	1658	+ G.set_filters(TFe, TFv)
	1659	+ G1 = gt.Graph(G, directed=False, prune=True)
	1660	+ G.clear_filters()
	1661	+ G1.clear_filters()
	1662	+ iteration = 0
	1663	+ title = "graph_" + str(iteration) + ".pdf"
	1664	+ gt.graph_draw(G1, pos = gt.sfdp_layout(G1), output=title)
	1665	+ while(np.any(G1.degree_property_map("total").get_array() == True)):
	1666	+ #print(np.any(G1.degree_property_map("total").get_array()), " ", i)
	1667	+ iteration = iteration + 1
	1668	+ G1 = self.filterBorder(G1, is_dual)
	1669	+ title = "graph_" + str(iteration) + ".pdf"
	1670	+ gt.graph_draw(G1, pos = gt.sfdp_layout(G1), output=title)
	1671	+ #G1 = self.filterDisconnected(G1, is_dual)
	1672	+ G1 = self.recalculate_metrics(G1, is_dual)
	1673	+
	1674	+ if(G1.num_edges() > 2):
	1675	+ ppath = path + "/partial"+str(j)+".nwt"
	1676	+ self.saveGraph_gt(G1, ppath)
	1677	+ if(i == 0):
	1678	+ self.saveSample_gt(G1, path, prefix, is_dual, writeKey=saveKey, label=lbl)
	1679	+ else:
	1680	+ self.saveSample_gt(G1, path, prefix, is_dual, label=lbl)
	1681	+ j = j + 1
	1682	+
	1683	+
	1684	+ #print(array_cluster_bel)
	1685	+
	1686	+
	1687	+ def find_loops(self, G, path, is_dual, largest):
	1688	+ iteration = 0
	1689	+ if(largest):
	1690	+ G = self.filterDisconnected(G)
	1691	+ gt.graph_draw(G, pos = gt.sfdp_layout(G), output="./graph_full.pdf")
	1692	+ while(np.any(G.degree_property_map("total").get_array() == True)):
	1693	+ #print(np.any(G1.degree_property_map("total").get_array()), " ", i)
	1694	+ iteration = iteration + 1
	1695	+ G = self.filterBorder(G, is_dual)
	1696	+ title = path +"graph_" + str(iteration) + ".pdf"
	1697	+ gt.graph_draw(G, pos = gt.sfdp_layout(G), output=title)
	1698	+ #G1 = self.filterDisconnected(G1, is_dual)
	1699	+ G = self.recalculate_metrics(G, is_dual)
	1700	+ ppath = path + "graph_full_mst.pdf"
	1701	+ 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)
	1702	+ paths = G.new_edge_property("int", vals=np.full((G.num_edges(), 1), 1, dtype=int))
	1703	+ #values = np.full((G.num_vertices(),1), 0, dtype=int)
	1704	+ #for i in range(G.num_vertices()):
	1705	+ # values[i,0] = G.vertex_index[G.vertex(i)]
	1706	+ #G.vertex_properties["original_index"] = values
	1707	+
	1708	+ #add index, then use index to search for new vertecies.
	1709	+ G.edge_properties["path_weights"] = paths
	1710	+ cycle_list = []
	1711	+ for e in G.edges():
	1712	+ if G.edge_properties["mst"][e] == 0:
	1713	+ cycle_list.append(e)
	1714	+ # else:
	1715	+ # G.edge_properties["path_weights"][e] = 1
	1716	+
	1717	+
	1718	+ #for e in cycle_list:
	1719	+ #print("Loop Detected between %i, %i", e.source(), e.target())
	1720	+
	1721	+ G.set_edge_filter(G.edge_properties["mst"])
	1722	+ iteration = 0
	1723	+ G1 = gt.Graph(G, directed=False, prune=True)
	1724	+ G.clear_filters()
	1725	+ for e in cycle_list:
	1726	+ edge = G1.add_edge(e.target(), e.source())
	1727	+ #G.edge_properties["mst"][e] = 1
	1728	+ G1.edge_properties["path_weights"][edge] = 1000
	1729	+ #G.set_edge_filter(G.edge_properties["mst"])
	1730	+ #G1 = gt.Graph(G, directed=False, prune=True)
	1731	+ #s = -1
	1732	+ #t = -1
	1733	+ #for edge in G1.edges():
	1734	+ # if(G1.vertex_properties["original_index"][edge.source()] == int(e.source())):
	1735	+ # s = int(edge.source())
	1736	+ # elif(G1.vertex_properties["original_index"][edge.target()] == int(e.target())):
	1737	+ # t = int(edge.target())
	1738	+
	1739	+ [Vl, El] = gt.shortest_path(G1, e.source(), e.target(), weights=G1.edge_properties["path_weights"])
	1740	+ clusters = G1.new_vertex_property("int", vals=np.full((G1.num_vertices(), 1), 5, dtype=int))
	1741	+ eclusters = G1.new_edge_property("int", vals = np.full((G1.num_edges(), 1), 3, dtype=int))
	1742	+ G1.vertex_properties["cycle"] = clusters
	1743	+ G1.edge_properties["cycle"] = eclusters
	1744	+ print("Number of vertices in Path is:", len(Vl))
	1745	+ for v in Vl:
	1746	+ G1.vertex_properties["cycle"][G1.vertex(v)] = 10
	1747	+ print(str(v))
	1748	+ #Create the arrays to be histogrammed later regarding every loop
	1749	+ length_total = 0
	1750	+ volume_total = 0
	1751	+ tortuosity_total = 0
	1752	+ bc_total = 0
	1753	+ bcl_total = 0
	1754	+ num_edges = 1
	1755	+ for v in range(len(Vl)-1):
	1756	+ G1.edge_properties["cycle"][G1.edge(Vl[v], Vl[v+1])] = 10
	1757	+ length_total = length_total + G1.edge_properties["length"][G1.edge(Vl[v], Vl[v+1])]
	1758	+ volume_total = volume_total + G1.edge_properties["volume"][G1.edge(Vl[v], Vl[v+1])]
	1759	+ tortuosity_total = tortuosity_total + G1.edge_properties["tortuosity"][G1.edge(Vl[v], Vl[v+1])]
	1760	+ bc_total = bc_total + G1.edge_properties["bc"][G1.edge(Vl[v], Vl[v+1])]
	1761	+ bcl_total = bcl_total + G1.edge_properties["bc*length"][G1.edge(Vl[v], Vl[v+1])]
	1762	+ num_edges = num_edges + 1
	1763	+ G1.edge_properties["cycle"][G1.edge(e.source(), e.target())] = 10
	1764	+ length_total = length_total + G.edge_properties["length"][G.edge(e.source(), e.target())]
	1765	+ volume_total = volume_total + G.edge_properties["volume"][G.edge(e.source(), e.target())]
	1766	+ tortuosity_total = tortuosity_total + G.edge_properties["tortuosity"][G.edge(e.source(), e.target())]
	1767	+ bc_total = bc_total + G.edge_properties["bc"][G.edge(e.source(), e.target())]
	1768	+ bcl_total = bcl_total + G.edge_properties["bc*length"][G.edge(e.source(), e.target())]
	1769	+ ppath = path + "Loop_Histograms.txt"
	1770	+ f = open(ppath, "a+")
	1771	+ f.write("%.15f\t" % length_total)
	1772	+ f.write("%.15f\t" % volume_total)
	1773	+ f.write("%.15f\t" % tortuosity_total)
	1774	+ f.write("%.15f\t" % bc_total)
	1775	+ f.write("%.15f\t" % bcl_total)
	1776	+ f.write("%d\t\n" % num_edges)
	1777	+ f.close()
	1778	+
	1779	+ title = path+"cycle" + str(iteration)+".png"
	1780	+ 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)
	1781	+ G1.remove_edge(edge)
	1782	+ #G.clear_filters();
	1783	+ #G.edge_properties["mst"][e] = 0
	1784	+ #G.edge_properties["path_weights"][e] = 1000
	1785	+ iteration = iteration + 1
	1786	+ '''
	1787	+ returns an affinity matrix based on the property edge property, given as a string.
	1788	+ '''
	1789	+ def get_affinity_matrix(G, edge_property, gaussian=False, sigma = None):
	1790	+ affinity = gt.shortest_distance(G, weights=G.edge_properties[edge_property])
	1791	+ affinity = affinity.get_2d_array(range(G.num_vertices()))
	1792	+ if gaussian:
	1793	+ if sigma == None:
	1794	+ sig = 60.0
	1795	+ else:
	1796	+ sig = sigma
	1797	+ for i in range(G.num_vertices()):
	1798	+ for j in range(G.num_vertices()):
	1799	+ affinity[i,j] = np.exp(-np.power(affinity[i, j], 2.) / (2 * np.power(sig, 2.)))
	1800	+ return affinity
	1801	+
	1802	+
	1803	+ '''
	1804	+ Attempts to apporixmate the number of clusters in the graph by finding
	1805	+ the approximate number of minimum graph cuts
	1806	+
	1807	+ TO_DO: Program the gaussian and sigma parameters
	1808	+ '''
	1809	+ def approximate_n_cuts(G, affinity_property):
	1810	+ #L = gt.laplacian(G, weight = G.edge_properties[affinity_property])
	1811	+ G1 = gt.Graph(G, directed=False)
	1812	+ G1 = Network.filterDisconnected(G1, G1)
	1813	+ L = Network.get_affinity_matrix(G1, affinity_property, gaussian=True, sigma=157)
	1814	+ L = sp.sparse.csgraph.laplacian(L, normed=True)
	1815	+ n_components = G1.num_vertices()
	1816	+ eigenvalues, eigenvectors = eigsh(L, k=n_components, which = "LM", sigma=1.0, maxiter = 5000)
	1817	+
	1818	+ plt.figure()
	1819	+ plt.scatter(np.arange(len(eigenvalues)), eigenvalues)
	1820	+ plt.grid()
	1821	+ plt.show()
	1822	+
	1823	+ count = sum(eigenvalues > 1.01)
	1824	+ return count
	1825	+ '''
	1826	+ Cluster the vectices in the graph based on a property passed to the clustering algorithm
	1827	+
	1828	+ TO_DO: program the guassian and sigma parameters
	1829	+ '''
	1830	+ def spectral_clustering(G, affinity_property, gaussian = False, sigma = None, n_clusters = None, num_iters = 1000):
	1831	+ if(n_clusters == None):
	1832	+ n_c = Network.approximate_n_cuts(G, affinity_property);
	1833	+ else:
	1834	+ n_c = n_clusters
	1835	+ sc = SpectralClustering(n_c, affinity='precomputed', n_init = num_iters)
	1836	+ sc.fit(Network.get_affinity_matrix(G, affinity_property, gaussian = True, sigma=157))
	1837	+
	1838	+ return sc.labels_
	1839	+
	1840	+
	1841	+ def export_points(self, G):
	1842	+ location = "./vertex_points.txt"
	1843	+ f = open(location, "a+")
	1844	+ for v in G.vertices():
	1845	+ f.write("%.15f\t" % G.vertex_properties["p"][v][0])
	1846	+ f.write("%.15f\t" % G.vertex_properties["p"][v][1])
	1847	+ f.write("%.15f\n" % G.vertex_properties["p"][v][2])
	1848	+
	1849	+ f.close()
	1850	+
	1851	+ def export_bb(self, G):
	1852	+ location = "./bb_points.txt"
	1853	+ aabb = AABB(G)
	1854	+ points = aabb.vertices()
	1855	+ f = open(location, "a+")
	1856	+ for i in points:
	1857	+ f.write("%.15f\t" % i[0])
	1858	+ f.write("%.15f\t" % i[1])
	1859	+ f.write("%.15f\n" % i[2])
	1860	+
	1861	+ f.close()
	1862	+
	1863	+
	1864	+ def write_vtk(G, location, camera = None, binning = True):
	1865	+ #location = "./nwt.vtk"
	1866	+ f = open(location, "w+")
	1867	+ f.write("# vtk DataFile Version 1.0\n")
	1868	+ f.write("Data values\n")
	1869	+ f.write("ASCII\n\n")
	1870	+ f.write("DATASET POLYDATA\n")
	1871	+ num_pts = 0
	1872	+ for e in G.edges():
	1873	+ X = G.edge_properties["x"][e]
	1874	+ num_pts += len(X)
	1875	+ f.write("POINTS %d float\n" % num_pts)
	1876	+ for e in G.edges():
	1877	+ X = G.edge_properties["x"][e]
	1878	+ Y = G.edge_properties["y"][e]
	1879	+ Z = G.edge_properties["z"][e]
	1880	+ #pts = np.array([X,Y,Z]).T
	1881	+ for p in range(len(X)):
	1882	+ f.write("%.15f\t" % X[p])
	1883	+ f.write("%.15f\t" % Y[p])
	1884	+ f.write("%.15f\n" % Z[p])
	1885	+ f.write("LINES " + str(G.num_edges()) + " " + str(G.num_edges()+num_pts) + "\n")
	1886	+ num_pts = 0
	1887	+ for e in G.edges():
	1888	+ X = G.edge_properties["x"][e]
	1889	+ indices = list(range(0, len(X)))
	1890	+ indices = [x+num_pts for x in indices]
	1891	+ f.write(str(len(X)) + " ")
	1892	+ f.write(" ".join(str(x) for x in indices))
	1893	+ f.write("\n")
	1894	+ num_pts += len(X)
	1895	+ f.write("POINT_DATA %d\n" % num_pts)
	1896	+ f.write("SCALARS radius float 1\n")
	1897	+ f.write("LOOKUP_TABLE radius_table\n")
	1898	+ for e in G.edges():
	1899	+ R = G.edge_properties["r"][e]
	1900	+ #pts = np.array([X,Y,Z]).T
	1901	+ for p in range(len(R)):
	1902	+ f.write("%.15f\n" % R[p])
	1903	+
	1904	+ f.write("COLOR_SCALARS color_table %d\n" % 4)
	1905	+ bins = [0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0]
	1906	+
	1907	+ if camera == None:
	1908	+ for e in G.edges():
	1909	+ X = G.edge_properties["x"][e]
	1910	+ #RGBA = G.edge_properties["RGBA"].get_2d_array(range(4)).T
	1911	+ for p in range(len(X)):
	1912	+ f.write("%.15f\t" % G.edge_properties["RGBA"][e][0])
	1913	+ f.write("%.15f\t" % G.edge_properties["RGBA"][e][1])
	1914	+ f.write("%.15f\t" % G.edge_properties["RGBA"][e][2])
	1915	+ if binning:
	1916	+ index = np.digitize(G.edge_properties["RGBA"][e][3], bins, right=True)
	1917	+ if (index >= len(bins) or index < 0):
	1918	+ print(G.edge_properties["RGBA"][e][3])
	1919	+ print(index)
	1920	+ f.write("%.15f\n" % bins[index])
	1921	+ else:
	1922	+ f.write("%.15f\n" % G.edge_properties["RGBA"][e][3])
	1923	+ else:
	1924	+ ptx = G.graph_properties["point_cloud_x"].get_array()
	1925	+ pty = G.graph_properties["point_cloud_y"].get_array()
	1926	+ ptz = G.graph_properties["point_cloud_z"].get_array()
	1927	+ pts = np.array([ptx,pty,ptz]).T
	1928	+ temp = AABB(G)
	1929	+ bbl = temp.A
	1930	+ bbu = temp.B
	1931	+ tp = np.zeros((1,3), dtype = np.float64)
	1932	+ tp[0][0] = (bbu[0]-bbl[0])/2
	1933	+ tp[0][1] = (bbu[1]-bbl[1])/2
	1934	+ tp[0][2] = (bbu[2]-bbl[2])/2
	1935	+ pts[:] = pts[:] - camera - \
	1936	+ tp
	1937	+ distance = np.zeros(pts.shape[0], dtype = np.float32)
	1938	+ for i in range(distance.shape[0]):
	1939	+ distance[i] = np.sqrt(np.power(pts[i][0],2) + np.power(pts[i][1],2) + np.power(pts[i][2],2))
	1940	+ mx = max(distance)
	1941	+ mn = min(distance)
	1942	+ distance = (distance - mx)/(mn-mx)
	1943	+ idx = 0
	1944	+ for e in G.edges():
	1945	+ X = G.edge_properties["x"][e]
	1946	+ for p in range(len(X)):
	1947	+ f.write("%.15f\t" % G.edge_properties["RGBA"][e][0])
	1948	+ f.write("%.15f\t" % G.edge_properties["RGBA"][e][1])
	1949	+ f.write("%.15f\t" % G.edge_properties["RGBA"][e][2])
	1950	+ if binning:
	1951	+ index = np.digitize(distance[idx], bins, right=True)
	1952	+ f.write("%.15f\n" % bins[index])
	1953	+ else:
	1954	+ f.write("%.15f\n" % distance[idx])
	1955	+ idx += 1
	1956	+
	1957	+ f.close()
	1958	+
	1959	+
	1960	+ '''
	1961	+ Creates a graph from a list of nodes and a list of edges.
	1962	+ Uses edge turtuosity as weight.
	1963	+ Returns a NetworkX Object.
	1964	+ '''
	1965	+ def createTortuosityGraph(self):
	1966	+ G = nx.Graph()
	1967	+ for i in range(len(self.N)):
	1968	+ G.add_node(i, p=self.N[i].p)
	1969	+ for i in range(len(self.F)):
	1970	+ G.add_edge(self.F[i].v0, self.F[i].v1, weight = self.F[i].turtuosity())
	1971	+ G[self.F[i].v0][self.F[i].v1]['pts'] = self.F[i].points
	1972	+ G[self.F[i].v0][self.F[i].v1]['rads'] = self.F[i].radii
	1973	+
	1974	+ return G
	1975	+
	1976	+
	1977	+ '''
	1978	+ Creates a graph from a list of nodes and a list of edges.
	1979	+ Uses edge volume as weight.
	1980	+ Returns a NetworkX Object.
	1981	+ '''
	1982	+ def createVolumeGraph(self):
	1983	+ G = nx.Graph()
	1984	+ for i in range(len(self.N)):
	1985	+ G.add_node(i, p=self.N[i].p)
	1986	+ for i in range(len(self.F)):
	1987	+ G.add_edge(self.F[i].v0, self.F[i].v1, weight = self.F[i].volume())
	1988	+ G[self.F[i].v0][self.F[i].v1]['pts'] = self.F[i].points
	1989	+ G[self.F[i].v0][self.F[i].v1]['rads'] = self.F[i].radii
	1990	+
	1991	+ return G
	1992	+#'''
	1993	+#Returns the positions dictionary for the Circular layout.
	1994	+#'''
	1995	+#def getCircularLayout(graph, dim, radius):
	1996	+# return nx.circular_layout(graph, dim, radius)
	1997	+#
	1998	+#'''
	1999	+#Return the positions dictionary for the Spring layout.
	2000	+#'''
	2001	+#def getSpringLayout(graph, pos, iterations, scale):
	2002	+# return nx.spring_layout(graph, 2, None, pos, iterations, 'weight', scale, None)
	2003	+#
	2004	+#'''
	2005	+#Draws the graph.
	2006	+#'''
	2007	+#def drawGraph(graph, pos):
	2008	+# nx.draw(graph, pos)
	2009	+# return
	2010	+
	2011	+ def aabb(self):
	2012	+
	2013	+ lower = self.N[0].p.copy()
	2014	+ upper = lower.copy()
	2015	+ for i in self.N:
	2016	+ for c in range(len(lower)):
	2017	+ if lower[c] > i.p[c]:
	2018	+ lower[c] = i.p[c]
	2019	+ if upper[c] < i.p[c]:
	2020	+ upper[c] = i.p[c]
	2021	+ return lower, upper
	2022	+
	2023	+ #calculate the distance field at a given resolution
	2024	+ # R (triple) resolution along each dimension
	2025	+ def distancefield(self, R=(100, 100, 100)):
	2026	+
	2027	+ #get a list of all node positions in the network
	2028	+ P = []
	2029	+ for e in self.F:
	2030	+ for p in e.points:
	2031	+ P.append(p)
	2032	+
	2033	+ #turn that list into a Numpy array so that we can create a KD tree
	2034	+ P = np.array(P)
	2035	+
	2036	+ #generate a KD-Tree out of the network point array
	2037	+ tree = sp.spatial.cKDTree(P)
	2038	+
	2039	+ plt.scatter(P[:, 0], P[:, 1])
	2040	+
	2041	+ #specify the resolution of the ouput grid
	2042	+ R = (200, 200, 200)
	2043	+
	2044	+ #generate a meshgrid of the appropriate size and resolution to surround the network
	2045	+ lower, upper = self.aabb(self.N, self.F) #get the space occupied by the network
	2046	+ x = np.linspace(lower[0], upper[0], R[0]) #get the grid points for uniform sampling of this space
	2047	+ y = np.linspace(lower[1], upper[1], R[1])
	2048	+ z = np.linspace(lower[2], upper[2], R[2])
	2049	+ X, Y, Z = np.meshgrid(x, y, z)
	2050	+ #Z = 150 * numpy.ones(X.shape)
	2051	+
	2052	+
	2053	+ Q = np.stack((X, Y, Z), 3)
	2054	+
	2055	+
	2056	+ D, I = tree.query(Q)
	2057	+
	2058	+ return D
	2059	+
	2060	+ #returns the number of points in the network
	2061	+ def npoints(self):
	2062	+ n = 0 #initialize the counter to zero
	2063	+ for f in self.F: #for each fiber
	2064	+ n = n + len(f.points) - 2 #count the number of points in the fiber - ignoring the end points
	2065	+ n = n + len(self.N) #add the number of nodes (shared points) to the node count
	2066	+ return n #return the number of nodes
	2067	+
	2068	+ #returns all of the points in the network
	2069	+ def points(self):
	2070	+ k = self.npoints()
	2071	+ P = np.zeros((3, k)) #allocate space for the point list
	2072	+
	2073	+ idx = 0
	2074	+ for f in self.F: #for each fiber in the network
	2075	+ for ip in range(1, len(f.points)-1): #for each point in the network
	2076	+ P[:, idx] = f.points[ip] #store the point in the raw point list
	2077	+ idx = idx + 1
	2078	+ return P #return the point array
	2079	+
	2080	+ #returns the number of linear segments in the network
	2081	+ def nsegments(self):
	2082	+ n = 0 #initialize the segment counter to 0
	2083	+ for f in self.F: #for each fiber
	2084	+ n = n + len(f.points) - 1 #calculate the number of line segments in the fiber (points - 1)
	2085	+ return n #return the number of line segments
	2086	+
	2087	+ #return a list of line segments representing the network
	2088	+ def segments(self, dtype=np.float32):
	2089	+ k = self.nsegments() #get the number of line segments
	2090	+ start = np.zeros((k, 3),dtype=dtype) #start points for the line segments
	2091	+ end = np.zeros((k, 3), dtype=dtype) #end points for the line segments
	2092	+
	2093	+ idx = 0 #initialize the index counter to zero
	2094	+ for f in self.F: #for each fiber in the network
	2095	+ for ip in range(0, len(f.points)-1): #for each point in the network
	2096	+ start[idx, :] = f.points[ip] #store the point in the raw point list
	2097	+ idx = idx + 1
	2098	+
	2099	+ idx = 0
	2100	+ for f in self.F: #for each fiber in the network
	2101	+ for ip in range(1, len(f.points)): #for each point in the network
	2102	+ end[idx, :] = f.points[ip] #store the point in the raw point list
	2103	+ idx = idx + 1
	2104	+
	2105	+ return start, end
	2106	+
	2107	+ #returns the points inside the fiber.
	2108	+ def fiber(self, idx):
	2109	+ p = self.F[idx].points
	2110	+ X = np.zeros(len(p))
	2111	+ Y = np.zeros(len(p))
	2112	+ Z = np.zeros(len(p))
	2113	+ for i in range(len(p)):
	2114	+ X[i] = p[i][0]
	2115	+ Y[i] = p[i][1]
	2116	+ Z[i] = p[i][2]
	2117	+ return X, Y, Z
	2118	+
	2119	+ #function returns the fiber associated with a given 1D line segment index
	2120	+ def segment2fiber(self, idx):
	2121	+ i = 0
	2122	+ for f in range(len(self.F)): #for each fiber in the network
	2123	+ i = i + len(self.F[f].points)-1 #add the number of points in the fiber to i
	2124	+ if i > idx: #if we encounter idx in this fiber
	2125	+ return self.F[f].points, f #return the fiber associated with idx and the index into the fiber array
	2126	+
	2127	+ def vectors(self, clock=False, dtype=np.float32):
	2128	+ if clock:
	2129	+ start_time = time.time()
	2130	+ start, end = self.segments(dtype) #retrieve all of the line segments
	2131	+ v = end - start #calculate the resulting vectors
	2132	+ l = np.sqrt(v[:, 0]2 + v[:,1]2 + v[:,2]**2) #calculate the fiber lengths
	2133	+ z = l==0 #look for any zero values
	2134	+ nz = z.sum()
	2135	+ if nz > 0:
	2136	+ print("WARNING: " + str(nz) + " line segment(s) of length zero were found in the network and will be removed" )
	2137	+
	2138	+ if clock:
	2139	+ print("Network::vectors: " + str(time.time() - start_time) + "s")
	2140	+
	2141	+ return np.delete(v, np.where(z), 0)
	2142	+
	2143	+ #scale all values in the network by tuple S = (sx, sy, sz)
	2144	+ def scale(self, S):
	2145	+ for f in self.F:
	2146	+ for p in f.points:
	2147	+ p[0] = p[0] * S[0]
	2148	+ p[1] = p[1] * S[1]
	2149	+ p[2] = p[2] * S[2]
	2150	+
	2151	+ for n in self.N:
	2152	+ n.p[0] = n.p[0] * S[0]
	2153	+ n.p[1] = n.p[1] * S[1]
	2154	+ n.p[2] = n.p[2] * S[2]
	2155	+
	2156	+
	2157	+ #calculate the adjacency weighting function for the network given a set of vectors X = (x, y, z) and weight exponent k
	2158	+ def adjacencyweight(self, P, k=200, length_threshold = 25, dtype=np.float32):
	2159	+ V = self.vectors(dtype) #get the vectors representing each segment
	2160	+ #V = V[0:n_vectors, :]
	2161	+ L = np.expand_dims(np.sqrt((V**2).sum(1)), 1) #calculate the length of each vector
	2162	+
	2163	+ outliers = L > length_threshold #remove outliers based on the length_threshold
	2164	+ V = np.delete(V, np.where(outliers), 0)
	2165	+ L = np.delete(L, np.where(outliers))
	2166	+ V = V/L[:,None] #normalize the vectors
	2167	+
	2168	+ P = np.stack(spharmonics.sph2cart(1, P[0], P[1]), P[0].ndim)
	2169	+ PV = P[...,None,:] * V
	2170	+ cos_alpha = PV.sum(PV.ndim-1)
	2171	+ W = np.abs(cos_alpha) ** k
	2172	+
	2173	+ return W, L
	2174	+
	2175	+
...	...

subgraph_shaders.py 0 → 100644

Wrap text Show/Hide comments View file @9f9f178

	1	+++ a/subgraph_shaders.py
	1	+#!/usr/bin/env python3
	2	+# -- coding: utf-8 --
	3	+"""
	4	+Created on Mon Aug 5 15:35:04 2019
	5	+
	6	+@author: pavel
	7	+"""
	8	+
	9	+#shaders for drawing supernodes
	10	+vert_s = """
	11	+#version 120
	12	+// Uniforms
	13	+// ------------------------------------
	14	+uniform mat4 u_model;
	15	+uniform mat4 u_view;
	16	+uniform mat4 u_projection;
	17	+uniform vec3 u_graph_size;
	18	+uniform bool u_picking;
	19	+// Attributes
	20	+// ------------------------------------
	21	+attribute vec3 a_position;
	22	+attribute vec4 a_bg_color;
	23	+attribute vec4 a_cluster_color;
	24	+attribute vec2 a_value;
	25	+attribute float a_arc_length;
	26	+attribute vec4 a_outer_arc_length;
	27	+attribute vec4 a_unique_id;
	28	+
	29	+// Varyings
	30	+// ------------------------------------
	31	+varying vec4 v_bg_color;
	32	+varying vec4 v_cluster_color;
	33	+varying vec2 v_value;
	34	+varying float v_zoom_level;
	35	+varying float v_arc_length;
	36	+varying vec4 v_outer_arc_length;
	37	+varying vec4 v_unique_id;
	38	+
	39	+void main (void) {
	40	+ v_value = a_value;
	41	+ v_bg_color = a_bg_color;
	42	+ v_cluster_color = a_cluster_color;
	43	+ gl_Position = u_projection * u_view * u_model *
	44	+ vec4(a_position, 1.0);
	45	+ v_zoom_level = u_view[0][0];
	46	+ v_arc_length = a_arc_length;
	47	+ v_outer_arc_length = a_outer_arc_length;
	48	+ v_unique_id = a_unique_id;
	49	+}
	50	+"""
	51	+
	52	+frag_s = """
	53	+#version 120
	54	+const float pi = 3.1415926535897932384626433832795;
	55	+// Varyings
	56	+// ------------------------------------
	57	+varying vec4 v_bg_color;
	58	+varying vec2 v_value;
	59	+varying float v_zoom_level;
	60	+varying vec4 v_cluster_color;
	61	+varying float v_arc_length;
	62	+varying vec4 v_outer_arc_length;
	63	+varying vec4 v_unique_id;
	64	+
	65	+uniform bool u_picking;
	66	+
	67	+
	68	+// Attributes
	69	+// ------------------------------------
	70	+
	71	+// Functions
	72	+// ------------------------------------
	73	+float new_alpha(float zoom_level);
	74	+float atan2(float y, float x);
	75	+// Main
	76	+// ------------------------------------
	77	+void main()
	78	+{
	79	+ if(!u_picking)
	80	+ {
	81	+ //This is the color of the outer arc lengths
	82	+ float R = 0.55;
	83	+ float R2 = 0.3;
	84	+ float dist = sqrt(dot(v_value, v_value));
	85	+ float sm = smoothstep(R, R-0.0010, dist);
	86	+ float sm2 = smoothstep(R2, R2+0.0010, dist);
	87	+ float alpha = sm*sm2;
	88	+
	89	+ float n_alpha = 1.0;
	90	+ if(v_zoom_level < 0.0010)
	91	+ {
	92	+ n_alpha = new_alpha(v_zoom_level);
	93	+ }
	94	+ else
	95	+ {
	96	+ discard;
	97	+ }
	98	+ float angle = atan(v_value[1], v_value[0]);
	99	+ if(dist < 0.25)
	100	+ {
	101	+ //gl_FragColor = v_cluster_color;
	102	+ gl_FragColor = vec4(v_cluster_color.rgb, n_alpha);
	103	+ }
	104	+ if(dist > 0.3 && 0.55 > dist)
	105	+ {
	106	+ float angle2 = angle+pi;
	107	+ if(angle2 > v_arc_length)
	108	+ {
	109	+ discard;
	110	+ }
	111	+ else
	112	+ {
	113	+ gl_FragColor = vec4(v_bg_color.rgb, n_alpha*alpha);
	114	+ }
	115	+ }
	116	+ angle = atan(v_value[1]/dist, v_value[0]/dist);
	117	+ //Need to add antialiasing to all other circles in the form of smoothstep.
	118	+
	119	+ //THIS NEED TO BE FIXED
	120	+ if(dist > 0.61 && 0.9 > dist)
	121	+ {
	122	+ if(angle < v_outer_arc_length[1])
	123	+ {
	124	+ gl_FragColor = vec4(0.32, 0.61, 0.93, n_alpha);
	125	+ //gl_FragColor = vec4(1.0, 0.0, 0.0, n_alpha);
	126	+ }
	127	+ else if(angle > v_outer_arc_length[1] && angle < v_outer_arc_length[2])
	128	+ {
	129	+ gl_FragColor = vec4(0.337, 0.866, 0.415, n_alpha);
	130	+ //gl_FragColor = vec4(0.0, 1.0, 0.0, n_alpha);
	131	+ }
	132	+ else if(angle > v_outer_arc_length[2] && angle < v_outer_arc_length[3])
	133	+ {
	134	+ gl_FragColor = vec4(0.988, 0.631, 0.058, n_alpha);
	135	+ //gl_FragColor = vec4(0.0, 0.0, 1.0, n_alpha);
	136	+ }
	137	+ else //if(angle > v_outer_arc_length[3] && angle < pi)
	138	+ {
	139	+ gl_FragColor = vec4(0.93, 0.949, 0.329, n_alpha);
	140	+ //gl_FragColor = vec4(1.0, 1.0, 0.0, n_alpha);
	141	+ }
	142	+ //else
	143	+ //{
	144	+ // discard;
	145	+ //}
	146	+
	147	+// if(angle > -pi && angle < -pi/4.0)
	148	+// {
	149	+// //gl_FragColor = vec4(0.32, 0.61, 0.93, n_alpha);
	150	+// gl_FragColor = vec4(1.0, 0.0, 0.0, n_alpha);
	151	+// }
	152	+// else if(angle > -pi/4.0 && angle < 0.0)
	153	+// {
	154	+// gl_FragColor = vec4(0.0, 1.0, 0.0, n_alpha);
	155	+// }
	156	+// else if(angle > 0.0 && angle < pi/2.0)
	157	+// {
	158	+// gl_FragColor = vec4(0.0, 0.0, 1.0, n_alpha);
	159	+// }
	160	+// else if(angle > pi/2.0 && angle < pi)
	161	+// {
	162	+// gl_FragColor = vec4(1.0, 1.0, 0.0, n_alpha);
	163	+// }
	164	+ }
	165	+ }
	166	+ else
	167	+ {
	168	+ float dist = sqrt(dot(v_value, v_value));
	169	+ if (dist > 0.9)
	170	+ discard;
	171	+ else
	172	+ gl_FragColor = v_unique_id;
	173	+ }
	174	+}
	175	+
	176	+float atan2(float y, float x)
	177	+{
	178	+ float s;
	179	+ if(abs(x) > abs(y))
	180	+ s = 1.0;
	181	+ else
	182	+ s = 0.0;
	183	+ return (mix(pi/2.0 - atan(x,y), atan(y,x), s));
	184	+}
	185	+
	186	+//float atan2(in float y, in float x)
	187	+//{
	188	+// return x == 0.0 ? sign(y)*pi/2 : atan(y, x);
	189	+//}
	190	+
	191	+float new_alpha(float zoom_level)
	192	+{
	193	+ float val = 1 - (zoom_level - 0.00075)/(0.0010-0.00075);
	194	+ if(val > 1.0)
	195	+ {
	196	+ val = 1.0;
	197	+ }
	198	+ return val;
	199	+}
	200	+"""
	201	+
	202	+vs_s = """
	203	+#version 120
	204	+
	205	+// Uniforms
	206	+// ------------------------------------
	207	+uniform mat4 u_model;
	208	+uniform mat4 u_view;
	209	+uniform mat4 u_projection;
	210	+
	211	+// Attributes
	212	+// ------------------------------------
	213	+attribute vec3 a_position;
	214	+attribute vec2 a_normal;
	215	+attribute vec4 a_fg_color;
	216	+attribute float a_linewidth;
	217	+//attribute vec4 a_unique_id;
	218	+//attribute vec4 l_color;
	219	+
	220	+// Varyings
	221	+// ------------------------------------
	222	+varying vec4 v_fg_color;
	223	+varying float v_zoom_level;
	224	+varying vec2 v_normal;
	225	+varying float v_linewidth;
	226	+
	227	+void main() {
	228	+ vec3 delta = vec3(a_normala_linewidth/((1-u_view[0][0])0.1), 0);
	229	+ //vec3 delta = vec3(a_normal*a_linewidth/(1-u_view[0][0]), 0);
	230	+ gl_Position = u_model * u_view * u_projection * vec4(a_position+delta, 1.0);
	231	+ //gl_Position = u_projection * u_view * u_model *
	232	+ // vec4(a_position, 1.0);
	233	+ v_zoom_level = u_view[0][0];
	234	+ v_fg_color = a_fg_color;
	235	+ v_normal = a_normal;
	236	+ v_linewidth = a_linewidth;
	237	+
	238	+}
	239	+"""
	240	+
	241	+fs_s = """
	242	+#version 120
	243	+// Varying
	244	+// ------------------------------------
	245	+varying vec4 v_fg_color;
	246	+varying float v_zoom_level;
	247	+varying vec2 v_normal;
	248	+varying float v_linewidth;
	249	+
	250	+float new_alpha(float zoom_level);
	251	+
	252	+void main()
	253	+{
	254	+ float l = length(v_normal);
	255	+ float feather = 0.5;
	256	+ float alpha = 1.0;
	257	+ if(l > v_linewidth/2.0+feather)
	258	+ discard;
	259	+ else
	260	+ alpha = 0.0;
	261	+
	262	+ if(v_zoom_level < 0.0010)
	263	+ alpha = new_alpha(v_zoom_level);
	264	+ gl_FragColor = vec4(v_fg_color.rgb, alpha);
	265	+}
	266	+
	267	+float new_alpha(float zoom_level)
	268	+{
	269	+ float val = 1 - (zoom_level - 0.00075)/(0.0010-0.00075);
	270	+ if(val > 1.0)
	271	+ {
	272	+ val = 1.0;
	273	+ }
	274	+ return val;
	275	+}
	276	+"""
0	277	\ No newline at end of file
...	...

tube_shaders.py 0 → 100644

Wrap text Show/Hide comments View file @9f9f178

	1	+++ a/tube_shaders.py
	1	+#!/usr/bin/env python3
	2	+# -- coding: utf-8 --
	3	+"""
	4	+Created on Mon Aug 5 15:58:30 2019
	5	+
	6	+@author: pavel
	7	+"""
	8	+
	9	+#vertex shader for the tube drawing
	10	+VERT_SHADER = """
	11	+// Uniforms
	12	+// ------------------------------------
	13	+uniform vec4 u_bb[26];
	14	+uniform vec3 u_LightPost;
	15	+uniform mat4 u_model;
	16	+//uniform mat4 u_view;
	17	+uniform mat4 u_projection;
	18	+
	19	+uniform vec3 u_eye;
	20	+uniform vec3 u_up;
	21	+uniform vec3 u_target;
	22	+
	23	+// Attributes
	24	+// ------------------------------------
	25	+attribute vec3 a_position;
	26	+attribute vec3 a_normal;
	27	+attribute vec4 a_fg_color;
	28	+
	29	+// Varyings
	30	+// ------------------------------------
	31	+varying vec3 v_normal;
	32	+varying vec4 v_fg_color;
	33	+varying vec3 v_position;
	34	+varying mat4 u_view;
	35	+varying vec3 vt_position;
	36	+
	37	+
	38	+// Functions
	39	+// ------------------------------------
	40	+
	41	+mat4 make_view(vec3 eye, vec3 target, vec3 up);
	42	+
	43	+void main (void) {
	44	+ // Calculate position
	45	+ mat4 u_view = make_view(u_eye, u_target, u_up);
	46	+
	47	+ mat4 MVP = u_projection * u_view * u_model;
	48	+ mat4 MV = u_view * u_model;
	49	+ v_normal = vec3(MV * vec4(a_normal, 0.0));
	50	+ v_position = vec3(MV*vec4(a_position, 1.0));
	51	+ vt_position = vec3(a_position);
	52	+ v_fg_color = a_fg_color;
	53	+ gl_Position = MVP * vec4(a_position, 1.0);
	54	+}
	55	+
	56	+mat4 make_view(vec3 eye, vec3 target, vec3 up)
	57	+{
	58	+ vec3 zaxis = normalize(eye - target);
	59	+ vec3 xaxis = normalize(cross(up, zaxis));
	60	+ vec3 yaxis = cross(zaxis, xaxis);
	61	+
	62	+ mat4 viewMatrix = mat4(
	63	+ vec4( xaxis.x, yaxis.x, zaxis.x, 0 ),
	64	+ vec4( xaxis.y, yaxis.y, zaxis.y, 0 ),
	65	+ vec4( xaxis.z, yaxis.z, zaxis.z, 0 ),
	66	+ vec4(-dot( xaxis, eye ), -dot( yaxis, eye ), -dot( zaxis, eye ), 1 )
	67	+ );
	68	+
	69	+// mat4 viewMatrix = mat4(
	70	+// vec4( xaxis.x, xaxis.y, xaxis.z, -dot( xaxis, eye ) ),
	71	+// vec4( yaxis.x, yaxis.y, yaxis.z, -dot( yaxis, eye ) ),
	72	+// vec4( zaxis.x, zaxis.y, zaxis.z, -dot( zaxis, eye ) ),
	73	+// vec4( 0, 0, 0, 1)
	74	+// );
	75	+
	76	+
	77	+ return viewMatrix;
	78	+}
	79	+"""
	80	+
	81	+#Fragment shader for the tube drawing.
	82	+FRAG_SHADER = """
	83	+// Uniforms
	84	+// ------------------------------------
	85	+uniform vec3 u_bb[26];
	86	+uniform vec3 u_LightPos;
	87	+uniform mat4 u_model;
	88	+//uniform mat4 u_view;
	89	+uniform mat4 u_projection;
	90	+
	91	+uniform vec3 u_eye;
	92	+uniform vec3 u_up;
	93	+uniform vec3 u_target;
	94	+
	95	+// Varyings
	96	+// ------------------------------------
	97	+varying vec3 v_normal;
	98	+varying vec4 v_fg_color;
	99	+varying vec3 v_position;
	100	+varying mat4 u_view;
	101	+varying vec3 vt_position;
	102	+
	103	+float new_alpha();
	104	+float set_view();
	105	+float bin_alpha(float alpha);
	106	+
	107	+void main()
	108	+{
	109	+ //Ambient Lighting
	110	+ float ambientStrength = 0.5;
	111	+ vec4 ambient_color = ambientStrength*v_fg_color;
	112	+
	113	+
	114	+ mat4 MV = u_projection * u_view * u_model;
	115	+ //vec3 u_LightP = vec3(0., 0., 0.);
	116	+ vec3 u_LightP = vec3(u_eye);
	117	+ //mat4 inv = inverse(u_view);
	118	+
	119	+ //vec3 u_LightP = vec3(u_view[0][3], u_view[1][3], u_view[2][3]);
	120	+
	121	+ //Diffuse Lighting
	122	+ //float distance = length(u_LightP - v_position);
	123	+ vec3 norm = normalize(v_normal);
	124	+ vec3 lightVector = normalize(u_LightP - v_position);
	125	+ float diffuse = max(dot(norm, lightVector), 0.0);
	126	+
	127	+ //diffuse = diffuse * (1.0 / (1.0 + (0.25 * distance * distance)));
	128	+
	129	+ vec4 color = (ambient_color + diffuse)*v_fg_color;
	130	+ float alpha = new_alpha();
	131	+ //if(alpha == 1.0)
	132	+ // gl_FragColor = vec4(1.0, 0.0, 0.0, 1.0);
	133	+ // gl_FragColor = vec4(color.rgb, alpha);
	134	+ //else
	135	+ gl_FragColor = vec4(color.rgb, alpha);
	136	+ // gl_FragColor = vec4(1.0, 0.0, 0.0, 1.0);
	137	+ //gl_FragColor = v_fg_color;
	138	+}
	139	+
	140	+float new_alpha()
	141	+{
	142	+ float alpha = 0.0;
	143	+
	144	+ //find the min and the max
	145	+ float mx = -100000.0;
	146	+ float mn = 100000.0;
	147	+ vec3 u_light = vec3(u_eye);
	148	+ for(int i = 0; i < 26; i++)
	149	+ {
	150	+ //vec3 rot_bb = vec3(u_projection * u_view * u_model * vec4(u_bb[i], 1.0));
	151	+ vec3 rot_bb = u_bb[i];
	152	+ if(length(rot_bb - u_light) < mn)
	153	+ {
	154	+ mn = length(rot_bb - u_light);
	155	+ }
	156	+ if(length(rot_bb - u_light) > mx)
	157	+ {
	158	+ mx = length(rot_bb - u_light);
	159	+ }
	160	+ }
	161	+
	162	+ float l = length(u_light - vt_position);
	163	+ alpha = (l-mn)/(mx-mn);
	164	+ return bin_alpha(alpha);
	165	+// if(alpha > 0.9)
	166	+// return 1.0;
	167	+// else if(alpha > 0.8)
	168	+// return 0.9;
	169	+// else if(alpha > 0.7)
	170	+// return 0.8;
	171	+// else if (alpha > 0.6)
	172	+// return 0.7;
	173	+// else if (alpha > 0.5)
	174	+// return 0.6;
	175	+// else if (alpha > 0.4)
	176	+// return 0.5;
	177	+// else if (alpha > 0.3)
	178	+// return 0.4;
	179	+// else if (alpha > 0.2)
	180	+// return 0.3;
	181	+// else if (alpha > 0.1)
	182	+// return 0.2;
	183	+// else if (alpha > 0.0)
	184	+// return 0.0;
	185	+// else
	186	+// return alpha;
	187	+}
	188	+
	189	+float bin_alpha(float alpha)
	190	+{
	191	+
	192	+
	193	+ if(alpha > 0.8)
	194	+ return 0.0;
	195	+ else if(alpha > 0.6)
	196	+ return 0.1;
	197	+ else if (alpha > 0.4)
	198	+ return 0.4;
	199	+ else if (alpha > 0.2)
	200	+ return 0.8;
	201	+ else if (alpha > 0.0)
	202	+ return 1.0;
	203	+ else
	204	+ return 1.0;
	205	+
	206	+// if(alpha > 0.9)
	207	+// return 0.0;
	208	+// else if(alpha > 0.8)
	209	+// return 0.1;
	210	+// else if(alpha > 0.7)
	211	+// return 0.2;
	212	+// else if (alpha > 0.6)
	213	+// return 0.3;
	214	+// else if (alpha > 0.5)
	215	+// return 0.4;
	216	+// else if (alpha > 0.4)
	217	+// return 0.5;
	218	+// else if (alpha > 0.3)
	219	+// return 0.6;
	220	+// else if (alpha > 0.2)
	221	+// return 0.7;
	222	+// else if (alpha > 0.1)
	223	+// return 0.8;
	224	+// else if (alpha > 0.0)
	225	+// return 0.9;
	226	+// else
	227	+// return 1.0;
	228	+}
	229	+
	230	+"""
	231	+
	232	+
...	...