From 9f9f17883ce16db5a1cfa46f0dfa21853c0f4f60 Mon Sep 17 00:00:00 2001 From: Pavel Govyadinov Date: Mon, 5 Aug 2019 16:13:56 -0500 Subject: [PATCH] clead up version of the GUI submitted to vis 2019 --- GraphCanvas.py | 1382 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ GraphWidget.py | 184 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ GuiVisPy_tube.py | 112 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ TubeCanvas.py | 443 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ TubeWidget.py | 65 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ graph_shaders.py | 256 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ network_dep.py | 2175 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ subgraph_shaders.py | 276 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ tube_shaders.py | 232 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 9 files changed, 5125 insertions(+), 0 deletions(-) create mode 100644 GraphCanvas.py create mode 100644 GraphWidget.py create mode 100644 GuiVisPy_tube.py create mode 100644 TubeCanvas.py create mode 100644 TubeWidget.py create mode 100644 graph_shaders.py create mode 100755 network_dep.py create mode 100644 subgraph_shaders.py create mode 100644 tube_shaders.py diff --git a/GraphCanvas.py b/GraphCanvas.py new file mode 100644 index 0000000..cb60c20 --- /dev/null +++ b/GraphCanvas.py @@ -0,0 +1,1382 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Created on Mon Aug 5 15:43:59 2019 + +@author: pavel + The GraphCanvas class that extends the scene class in vispy in order to draw + the graph object. This class is wrapped in a QTcanvas. +""" + +from vispy import gloo, scene +from vispy.gloo import set_viewport, set_state, clear, set_blend_color +from vispy.util.transforms import perspective, translate, rotate, scale +import vispy.gloo.gl as glcore + +import numpy as np +import math +import network_dep as nwt + +#import the graph shaders. +from graph_shaders import vert, frag, vs, fs +from subgraph_shaders import vert_s, frag_s, vs_s, fs_s + +#The graph canvas class that +class GraphCanvas(scene.SceneCanvas): + + """ + Initialization method. + Generates the 512x512 canvas, makes it available for drawing the fills all + the GLSL shaders with dummy data. + """ + def __init__(self, **kwargs): + # Initialize the canvas for real + scene.SceneCanvas.__init__(self, size=(512, 512), **kwargs) + #Unfreeze the canvas to make dynamic interaction possible + self.unfreeze() + + #initialize all the boolean and dictionary variables. + ps = self.pixel_scale + self.subgraphs = False + self.position = 50, 50 + self.down=False; + + #Dictionaries to store the unique color ID, cluster ID, and edge-to-ID + #dictionaries. + self.color_dict = {} + self.cluster_dict = {} + self.edge_dict = {} + + #Booleans the storage for the current "Path", i.e. edges the user selected. + self.pathing = False + self.path = [] + self.full_path = [] + + #utility variables used for storing the cluster being moved and all the + #nodes and edges that belong to that cluster and move along with it. + self.moving = False + self.moving_cluster = False + self.selection = False + + n = 10 + ne = 10 + #Init dummy structures + self.uniforms = [('u_graph_size', np.float32, 3)] + self.data = np.zeros(n, dtype=[('a_position', np.float32, 3), + ('a_fg_color', np.float32, 4), + ('a_bg_color', np.float32, 4), + ('a_size', np.float32, 1), + ('a_linewidth', np.float32, 1), + ('a_unique_id', np.float32, 4), + ]) + + self.clusters = np.zeros(n, dtype=[('a_position', np.float32, 3), + ('a_bg_color', np.float32, 4), + ('a_value', np.float32, 2), + ('a_unique_id', np.float32, 4), + ]) + + self.edges = np.random.randint(size=(ne, 2), low=0, + high=n-1).astype(np.uint32) + self.edges_s = np.random.randint(size=(ne, 4), low=0, + high=n-1).astype(np.uint32) + self.data['a_position'] = np.hstack((0.25 * np.random.randn(n, 2), + np.zeros((n, 1)))) + self.data['a_fg_color'] = 0, 0, 0, 1.0 + color = np.random.uniform(0.5, 1., (n, 3)) + self.data['a_bg_color'] = np.hstack((color, np.zeros((n, 1)))) + self.data['a_size'] = np.random.randint(size=n, low=8*ps, high=20*ps) + self.data['a_linewidth'] = 8.*ps + self.data['a_unique_id'] = np.hstack((color, np.ones((n, 1)))) + #self.uniforms['u_graph_size'] = [1.0, 1.0, 1.0] + self.translate = [0,0,0] + self.scale = [1,1,1] + #color = np.random.uniform(0.5, 1., (ne, 3)) + #self.linecolor = np.hstack((color, np.ones((ne, 1)))) + #color = np.random.uniform(0.5, 1., (ne, 3)) + #self.linecolor = np.hstack((color, np.ones((ne, 1)))) + self.u_antialias = 1 + + #init dummy vertex and index buffers. + self.vbo = gloo.VertexBuffer(self.data) + self.vbo_s = gloo.VertexBuffer(self.clusters) + + #Need to initialize thick lines. + self.index = gloo.IndexBuffer(self.edges) + self.index_s = gloo.IndexBuffer(self.edges_s) + + #Set the view matrices. + self.view = np.eye(4, dtype=np.float32) + self.model = np.eye(4, dtype=np.float32) + self.projection = np.eye(4, dtype=np.float32) + + #init shaders used for vertices of the full graph. + self.program = gloo.Program(vert, frag) + self.program.bind(self.vbo) + self.program['u_size'] = 1 + self.program['u_antialias'] = self.u_antialias + self.program['u_model'] = self.model + self.program['u_view'] = self.view + self.program['u_projection'] = self.projection + self.program['u_graph_size'] = [1.0, 1.0, 1.0] + self.program['u_picking'] = False + + #init shades used for the edges in the graph + self.program_e = gloo.Program(vs, fs) + self.program_e['u_size'] = 1 + self.program_e['u_model'] = self.model + self.program_e['u_view'] = self.view + self.program_e['u_projection'] = self.projection + #self.program_e['l_color'] = self.linecolor.astype(np.float32) + self.program_e.bind(self.vbo) + + #init shaders used to the vertices in the subgraph graph. + self.program_s = gloo.Program(vert_s, frag_s) + self.program_s.bind(self.vbo_s) + self.program_s['u_model'] = self.model + self.program_s['u_view'] = self.view + self.program_s['u_projection'] = self.projection + self.program_s['u_graph_size'] = [1.0, 1.0, 1.0] + self.program_s['u_picking'] = False + + #init shaders used for the subgraph-edges + self.program_e_s = gloo.Program(vs_s, fs_s) + self.program_e_s['u_size'] = 1 + self.program_e_s['u_model'] = self.model + self.program_e_s['u_view'] = self.view + self.program_e_s['u_projection'] = self.projection + #self.program_e['l_color'] = self.linecolor.astype(np.float32) + self.program_e_s.bind(self.vbo_s) + + + #set up the viewport and the gl state. + set_viewport(0, 0, *self.physical_size) + + set_state(clear_color='white', depth_test=True, blend=True, + blend_func=('src_alpha', 'one_minus_src_alpha'), depth_func = ('lequal')) + + """ + Function that recolors vertices based on the selected statistic + Maps the statisic stored in G to a colormap passed to the function + Then updates the necessary color array. + """ + def color_vertices(self, G, vertex_property, dtype = False, cm = 'plasma'): + #if we are visualing the clusters we should use a discrete colormap + #otherwise use the passed colormap + if dtype == True: + G.vertex_properties["RGBA"] = nwt.Network.map_property_to_color(G, G.vertex_properties["clusters"]) + else: + G.vertex_properties["RGBA"] = nwt.Network.map_property_to_color(G, G.vertex_properties[vertex_property], colormap=cm) + #set the color and update the Vertices. + self.current_color = vertex_property + color = G.vertex_properties["RGBA"].get_2d_array(range(4)).T + self.data['a_bg_color'] = color + self.vbo = gloo.VertexBuffer(self.data) + self.program.bind(self.vbo) + #self.program_e.bind(self.vbo) + self.update() + + """ + Maps a statistic of the vertices based on the size of the canvas to size of + the drawn object. + """ + def size_vertices(self, G, propertymap): + size = nwt.Network.map_vertices_to_range(G, [30*self.pixel_scale, 8*self.pixel_scale], propertymap).get_array() + self.data['a_size'] = size + self.vbo = gloo.VertexBuffer(self.data) + self.program.bind(self.vbo) + #self.program_e.bind(self.vbo) + self.update() + + + """ + Function to dim all nodes and edges that do not belong to a cluster chosen + in the graph view. Returns a copy of the graph with the alpha channel saved. + OPTMIZE HERE: could just return an alpha array to reduce memory usage. + """ + def focus_on_cluster(self, G, c_id): + G_copy = nwt.gt.Graph(G, directed=False) + e_color = G_copy.edge_properties["RGBA"].get_2d_array(range(4)).T + vertices = np.argwhere(self.labels != c_id) + for v in range(vertices.shape[0]): + idx = vertices[v][0] + vtx = G_copy.vertex(idx) + for e in vtx.all_edges(): + if (int(e.source()), int(e.target())) in self.edge_dict.keys(): + index = int(self.edge_dict[int(e.source()), int(e.target())]) + if vtx == int(e.source()): + e_color[index][3] = 0.05 + elif vtx == int(e.target()): + e_color[index][3] = 0.05 + else: + index = int(self.edge_dict[int(e.target()), int(e.source())]) + if vtx == int(e.target()): + e_color[index][3] = 0.05 + elif vtx == int(e.source()): + e_color[index][3] = 0.05 + + G_copy.edge_properties["RGBA"] = G_copy.new_edge_property("vector", vals = e_color) + + return G_copy + + """ + Function that sets the size of the vertices based on the distance from the + camera. + """ + def vertexSizeFromDistance(self, G, camera_pos): + location = G.vertex_properties["p"].get_2d_array(range(3)).T + cam_array = np.zeros(location.shape, dtype=np.float32) + len_array = np.zeros(location.shape[0]) + offset_array = np.zeros(location.shape, dtype=np.float32) + cam_array[:][0:3] = camera_pos + offset = [(bbu[0]-bbl[0])/2, (bbu[1]-bbl[1])/2, (bbu[2]-bbl[2])/2] + location = location - offset + location = location - camera_pos + for i in range(location.shape[0]): + len_array[i] = np.sqrt(np.power(location[i][0],2) + np.power(location[i][1],2) + np.power(location[i][2],2)) + + G.vertex_properties['dist_from_camera'] = G.new_vertex_property('float', vals=len_array) + self.data['a_size'] = nwt.Network.map_vertices_to_range(G, [1*self.pixel_scale, 60*self.pixel_scale], 'dist_from_camera').get_array() + + + size = nwt.Network.map_vertices_to_range(G, [1.0, 0.5], 'dist_from_camera').get_array() + edges = G.get_edges() + for e in range(edges.shape[0]): + idx = int(4*edges[e][2]) + self.line_data['a_linewidth'][idx] = size[edges[e][0]] + self.line_data['a_linewidth'][idx+1] = size[edges[e][1]] + self.line_data['a_linewidth'][idx+2] = size[edges[e][0]] + self.line_data['a_linewidth'][idx+3] = size[edges[e][1]] + + + #self.vbo = gloo.VertexBuffer(self.data) + #self.vbo_line = gloo.VertexBuffer(self.line_data) + #self.program.bind(self.vbo) + #self.program_e.bind(self.vbo_line) + #self.update() + + """ + Function that scales the alpha channel of each vertex in the graph based on + The distance from the camera. + Sometimes needs to be done separetly. + """ + def vertexAlphaFromDistance(self, G, camera_pos): + location = G.vertex_properties["p"].get_2d_array(range(3)).T + cam_array = np.zeros(location.shape, dtype=np.float32) + len_array = np.zeros(location.shape[0]) + #offset_array = np.zeros(location.shape, dtype=np.float32) + cam_array[:][0:3] = camera_pos + offset = [(bbu[0]-bbl[0])/2, (bbu[1]-bbl[1])/2, (bbu[2]-bbl[2])/2] + location = location - offset + location = location - camera_pos + for i in range(location.shape[0]): + len_array[i] = np.sqrt(np.power(location[i][0],2) + np.power(location[i][1],2) + np.power(location[i][2],2)) + + + test = nwt.Network.map_vertices_to_range(G, [0.0, 1.0], 'dist_from_camera').get_array() + color = G.vertex_properties["RGBA"].get_2d_array(range(4)).T + for i in range(location.shape[0]): + color[i][3] = test[i] + G.vertex_properties["RGBA"] = G.new_vertex_property("vector", vals = color) + self.data['a_bg_color'] = color + + edges = G.get_edges() + for e in range(edges.shape[0]): + idx = int(4*edges[e][2]) + self.line_data['a_fg_color'][idx] = color[edges[e][0]] + self.line_data['a_fg_color'][idx+1] = color[edges[e][1]] + self.line_data['a_fg_color'][idx+2] = color[edges[e][0]] + self.line_data['a_fg_color'][idx+3] = color[edges[e][1]] + + + + self.vbo = gloo.VertexBuffer(self.data) + self.vbo_line = gloo.VertexBuffer(self.line_data) + self.program.bind(self.vbo) + self.program_e.bind(self.vbo_line) + self.update() + + """ + Sets the edge color based on the the cluster the vertices belongs to + Propertymap is a VERTEXPROPERTYMAP since the color of the edges is based + on the clusters the edges belong to. + """ + def color_edges(self, G, propertymap="clusters"): + if propertymap == "clusters": + for e in G.edges(): + if G.vertex_properties[propertymap][e.source()] == G.vertex_properties[propertymap][e.target()]: + G.edge_properties["RGBA"][e] = G.vertex_properties["RGBA"][e.source()] + else: + G.edge_properties["RGBA"][e] = [0.0, 0.0, 0.0, 0.8] + + + """ + Helper function that generates the framebuffer object that stores the vertices + Generates the vertex buffer based on the graph G that is passed to the function + Sets the color, generates the graph and subgraph color if necessary. + """ + def gen_vertex_vbo(self, G): + color = G.vertex_properties["RGBA"].get_2d_array(range(4)).T + size = nwt.Network.map_vertices_to_range(G, [30*self.pixel_scale, 8*self.pixel_scale], 'degree').get_array() + + position = G.vertex_properties["pos"].get_2d_array(range(3)).T + #for p in range(position.shape[0]): + # position[p][0] = position[p][0] + self.clusters["a_position"][G.vertex_properties["clusters"][G.vertex(p)]][0] + # position[p][1] = position[p][1] + self.clusters["a_position"][G.vertex_properties["clusters"][G.vertex(p)]][1] + # position[p][2] = position[p][2] + self.clusters["a_position"][G.vertex_properties["clusters"][G.vertex(p)]][2] + #G.vertex_properties["pos"] = G.new_vertex_property("vector", vals = position) + edges = G.get_edges(); + edges = edges[:, 0:2] + #width = nwt.Network.map_edges_to_range(G, [1*self.pixel_scale, 5*self.pixel_scale], 'volume').get_array() + #ecolor = G.edge_properties["RGBA"].get_2d_array(range(4)).T + + self.data = np.zeros(G.num_vertices(), dtype=[('a_position', np.float32, 3), + ('a_fg_color', np.float32, 4), + ('a_bg_color', np.float32, 4), + ('a_size', np.float32, 1), + ('a_linewidth', np.float32, 1), + ('a_unique_id', np.float32, 4), + ('a_selection', np.float32, 1), + ]) + + #self.edges = edges.astype(np.uint32) + self.data['a_position'] = position + #fg color is the color of the ring + self.data['a_fg_color'] = 0, 0, 0, 1 + self.data['a_bg_color'] = color + self.data['a_size'] = size + self.data['a_linewidth'] = 4.*self.pixel_scale + self.data['a_unique_id'] = self.gen_vertex_id(G) + self.data['a_selection'] = G.vertex_properties["selection"].get_array() + #self.data['a_graph_size'] = [bbu-bbl] + + self.program['u_graph_size'] = [bbu-bbl] + + self.vbo = gloo.VertexBuffer(self.data) + self.gen_line_vbo(G) + if(self.subgraphs): + self.vbo_s = gloo.VertexBuffer(self.clusters) + self.index_s = gloo.IndexBuffer(self.edges_s) + #self.index = gloo.IndexBuffer(self.edges) + self.program_e.bind(self.vbo_line) + self.program.bind(self.vbo) + if(self.subgraphs): + #self.program_e_s.bind(self.vbo_s) + self.program_s.bind(self.vbo_s) + print(self.view) + self.update() + + """ + Helper function that creates colored "block" lines based on the edges + in the graph. Generates the framebuffer object and fills it with the relavant data. + Note that each line segment is saved as a two triangles that share the same + two points on the centerline, but are offset according to the normal of the + line segmente to control thickness dynamically. + """ + def gen_line_vbo(self, G): + #Set the data. + self.line_data = np.zeros(G.num_edges()*4, dtype=[('a_position', np.float32, 3), + ('a_normal', np.float32, 2), + ('a_fg_color', np.float32, 4), + ('a_linewidth', np.float32, 1), + ]) + self.edges = np.random.randint(size=(G.num_edges()*2, 3), low=0, + high=(G.num_edges()-1)).astype(np.uint32) + color = G.edge_properties["RGBA"].get_2d_array(range(4)).T + edges = G.get_edges() + #size need to be changed to the size based on the current property map + size = nwt.Network.map_vertices_to_range(G, [1.0, 0.5], 'degree').get_array() + for e in range(edges.shape[0]): + idx = int(4*edges[e][2]) + p0 = G.vertex_properties["pos"][G.vertex(edges[e][0])] + p1 = G.vertex_properties["pos"][G.vertex(edges[e][1])] + d = np.subtract(p1, p0) + #d_norm = np.multiply(d, 1/np.sqrt(np.power(d[0],2) + np.power(d[1],2))) + d_norm = d[0:2] + d_norm = d_norm / np.sqrt(np.power(d_norm[0],2) + np.power(d_norm[1],2)) + norm = np.zeros((2,), dtype=np.float32) + norm[0] = d_norm[1] + norm[1] = d_norm[0]*-1 + #print(np.sqrt(norm[0]*norm[0] + norm[1]*norm[1])) + #thickness = G.edge_properties["thickness"][e] + thickness = 1.0 + self.edge_dict[int(edges[e][0]), int(edges[e][1])] = int(edges[e][2]) + self.line_data['a_position'][idx] = p0 + self.line_data['a_normal'][idx] = norm + self.line_data['a_fg_color'][idx] = color[edges[e][2]] + #a_linewidth is a vector. + self.line_data['a_linewidth'][idx] = size[edges[e][0]] + + self.line_data['a_position'][idx+1] = p1 + self.line_data['a_normal'][idx+1] = norm + self.line_data['a_fg_color'][idx+1] = color[edges[e][2]] + self.line_data['a_linewidth'][idx+1] = size[edges[e][1]] + + self.line_data['a_position'][idx+2] = p0 + self.line_data['a_normal'][idx+2] = -norm + self.line_data['a_fg_color'][idx+2] = color[edges[e][2]] + self.line_data['a_linewidth'][idx+2] = size[edges[e][0]] + + self.line_data['a_position'][idx+3] = p1 + self.line_data['a_normal'][idx+3] = -norm + self.line_data['a_fg_color'][idx+3] = color[edges[e][2]] + self.line_data['a_linewidth'][idx+3] = size[edges[e][1]] + + self.edges[e*2] = [idx, idx+1, idx+3] + self.edges[e*2+1] = [idx, idx+2, idx+3] + + #Set the buffer object and update the shader programs. + self.program_e = gloo.Program(vs, fs) + #self.program_e['l_color'] = self.linecolor.astype(np.float32) + self.vbo_line = gloo.VertexBuffer(self.line_data) + self.index = gloo.IndexBuffer(self.edges) + self.program_e['u_size'] = 1 + self.program_e['u_model'] = self.model + self.program_e['u_view'] = self.view + self.program_e['u_projection'] = self.projection + self.program_e.bind(self.vbo_line) + + + """ + Helper function that generates the edges between the cluster in the layout. + Color is based on the cluster source/target color and transitions between the + two. + """ + def gen_cluster_line_vbo(self, G): + #create a graph that stores the edges of between the clusters + self.G_cluster = nwt.gt.Graph(directed=False) + self.G_cluster.vertex_properties["pos"] = self.G_cluster.new_vertex_property("vector", val=np.zeros((3,1), dtype=np.float32)) + self.G_cluster.vertex_properties["RGBA"] = self.G_cluster.new_vertex_property("vector", val=np.zeros((4,1), dtype=np.float32)) + for v in range(len(self.cluster_pos)): + self.G_cluster.add_vertex() + self.G_cluster.vertex_properties["pos"][self.G_cluster.vertex(v)] = np.asarray(self.cluster_pos[v], dtype=np.float32) + self.G_cluster.edge_properties["weight"] = self.G_cluster.new_edge_property("int", val = 0) + #for each edge in the original graph, generate appropriate subgraph edges without repretiions + #i.e. controls the thichness of the edges in the subgraph view. + for e in G.edges(): + #if the source and target cluster is not equal to each other + #add an inter subgraph edge. + if(G.vertex_properties["clusters"][e.source()] != G.vertex_properties["clusters"][e.target()]): + #temp_e.append([G.vertex_properties["clusters"][e.source()], G.vertex_properties["clusters"][e.target()]]) + self.G_cluster.add_edge(self.G_cluster.vertex(G.vertex_properties["clusters"][e.source()]), \ + self.G_cluster.vertex(G.vertex_properties["clusters"][e.target()])) + self.G_cluster.edge_properties["weight"][self.G_cluster.edge(self.G_cluster.vertex(G.vertex_properties["clusters"][e.source()]), \ + self.G_cluster.vertex(G.vertex_properties["clusters"][e.target()]))] += 1 + self.G_cluster.vertex_properties["RGBA"][self.G_cluster.vertex(G.vertex_properties["clusters"][e.source()])] \ + = G.vertex_properties["RGBA"][e.source()] + self.G_cluster.vertex_properties["RGBA"][self.G_cluster.vertex(G.vertex_properties["clusters"][e.target()])] \ + = G.vertex_properties["RGBA"][e.target()] + + self.cluster_line_data = np.zeros(self.G_cluster.num_edges()*4, dtype=[('a_position', np.float32, 3), + ('a_normal', np.float32, 2), + ('a_fg_color', np.float32, 4), + ('a_linewidth', np.float32, 1), + ]) + self.cluster_edges = np.random.randint(size=(self.G_cluster.num_edges()*2, 3), low=0, + high=(G.num_edges()-1)).astype(np.uint32) + + edges = self.G_cluster.get_edges() + #size need to be changed to the size based on the current property map + size = nwt.Network.map_edges_to_range(self.G_cluster, [1.0, 0.5], 'weight').get_array() + color = self.G_cluster.vertex_properties["RGBA"].get_2d_array(range(4)).T + #generate the vertex buffer and the connections buffer. + for e in range(edges.shape[0]): + idx = int(4*edges[e][2]) + p0 = self.G_cluster.vertex_properties["pos"][self.G_cluster.vertex(edges[e][0])] + p1 = self.G_cluster.vertex_properties["pos"][self.G_cluster.vertex(edges[e][1])] + d = np.subtract(p1, p0) + #d_norm = np.multiply(d, 1/np.sqrt(np.power(d[0],2) + np.power(d[1],2))) + d_norm = d[0:2] + d_norm = d_norm / np.sqrt(np.power(d_norm[0],2) + np.power(d_norm[1],2)) + norm = np.zeros((2,), dtype=np.float32) + norm[0] = d_norm[1] + norm[1] = d_norm[0]*-1 + #print(np.sqrt(norm[0]*norm[0] + norm[1]*norm[1])) + #thickness = G.edge_properties["thickness"][e] + self.cluster_dict[int(edges[e][0]), int(edges[e][1])] = int(edges[e][2]) + self.cluster_line_data['a_position'][idx] = p0 + self.cluster_line_data['a_normal'][idx] = norm + self.cluster_line_data['a_fg_color'][idx] = color[edges[e][0]] + self.cluster_line_data['a_linewidth'][idx] = size[e] + + self.cluster_line_data['a_position'][idx+1] = p1 + self.cluster_line_data['a_normal'][idx+1] = norm + self.cluster_line_data['a_fg_color'][idx+1] = color[edges[e][1]] + self.cluster_line_data['a_linewidth'][idx+1] = size[e] + + self.cluster_line_data['a_position'][idx+2] = p0 + self.cluster_line_data['a_normal'][idx+2] = -norm + self.cluster_line_data['a_fg_color'][idx+2] = color[edges[e][0]] + self.cluster_line_data['a_linewidth'][idx+2] = size[e] + + self.cluster_line_data['a_position'][idx+3] = p1 + self.cluster_line_data['a_normal'][idx+3] = -norm + self.cluster_line_data['a_fg_color'][idx+3] = color[edges[e][1]] + self.cluster_line_data['a_linewidth'][idx+3] = size[e] + + self.cluster_edges[e*2] = [idx, idx+1, idx+3] + self.cluster_edges[e*2+1] = [idx, idx+2, idx+3] + + + self.program_e_s = gloo.Program(vs_s, fs_s) + self.index_clusters_s = gloo.IndexBuffer(self.cluster_edges) + self.vbo_cluster_lines = gloo.VertexBuffer(self.cluster_line_data) + + self.program_e_s['u_size'] = 1 + self.program_e_s['u_model'] = self.model + self.program_e_s['u_view'] = self.view + self.program_e_s['u_projection'] = self.projection + self.program_e_s.bind(self.vbo_cluster_lines) + + + """ + Updates the vertex buffers based on the current position of the cluster. + Updates it's position. + """ + def update_cluster_line_vbo(self): + + for v in range(len(self.cluster_pos)): + self.G_cluster.vertex_properties["pos"][self.G_cluster.vertex(v)] = np.asarray(self.cluster_pos[v], dtype=np.float32) + #OPTIMIZE HERE to update only one cluster at a time. + edges = self.G_cluster.get_edges() + #size need to be changed to the size based on the current property map + size = nwt.Network.map_edges_to_range(self.G_cluster, [1.0, 0.5], 'weight').get_array() + color = self.G_cluster.vertex_properties["RGBA"].get_2d_array(range(4)).T + for e in range(edges.shape[0]): + idx = int(4*edges[e][2]) + p0 = self.G_cluster.vertex_properties["pos"][self.G_cluster.vertex(edges[e][0])] + p1 = self.G_cluster.vertex_properties["pos"][self.G_cluster.vertex(edges[e][1])] + d = np.subtract(p1, p0) + #d_norm = np.multiply(d, 1/np.sqrt(np.power(d[0],2) + np.power(d[1],2))) + d_norm = d[0:2] + d_norm = d_norm / np.sqrt(np.power(d_norm[0],2) + np.power(d_norm[1],2)) + norm = np.zeros((2,), dtype=np.float32) + norm[0] = d_norm[1] + norm[1] = d_norm[0]*-1 + #print(np.sqrt(norm[0]*norm[0] + norm[1]*norm[1])) + #thickness = G.edge_properties["thickness"][e] + self.cluster_dict[int(edges[e][0]), int(edges[e][1])] = int(edges[e][2]) + self.cluster_line_data['a_position'][idx] = p0 + self.cluster_line_data['a_normal'][idx] = norm + self.cluster_line_data['a_fg_color'][idx] = color[edges[e][0]] + self.cluster_line_data['a_linewidth'][idx] = size[e] + + self.cluster_line_data['a_position'][idx+1] = p1 + self.cluster_line_data['a_normal'][idx+1] = norm + self.cluster_line_data['a_fg_color'][idx+1] = color[edges[e][1]] + self.cluster_line_data['a_linewidth'][idx+1] = size[e] + + self.cluster_line_data['a_position'][idx+2] = p0 + self.cluster_line_data['a_normal'][idx+2] = -norm + self.cluster_line_data['a_fg_color'][idx+2] = color[edges[e][0]] + self.cluster_line_data['a_linewidth'][idx+2] = size[e] + + self.cluster_line_data['a_position'][idx+3] = p1 + self.cluster_line_data['a_normal'][idx+3] = -norm + self.cluster_line_data['a_fg_color'][idx+3] = color[edges[e][1]] + self.cluster_line_data['a_linewidth'][idx+3] = size[e] + + self.program_e_s = gloo.Program(vs_s, fs_s) + self.index_clusters_s = gloo.IndexBuffer(self.cluster_edges) + self.vbo_cluster_lines = gloo.VertexBuffer(self.cluster_line_data) + + self.program_e_s['u_size'] = 1 + self.program_e_s['u_model'] = self.model + self.program_e_s['u_view'] = self.view + self.program_e_s['u_projection'] = self.projection + self.program_e_s.bind(self.vbo_cluster_lines) + + """ + Genererates a unique index for every vertex. + """ + def gen_vertex_id(self, G): + self.color_dict = {} + base = [0, 0, 0, 255] + idx = 0 + #colors = cm.get_cmap('Wistia', G.num_vertices()*2) + v_id = np.zeros((G.num_vertices(), 4), dtype=np.float32) + for v in G.vertices(): + color = np.multiply(base, 1/255.0) + v_id[int(v)] = color + self.color_dict[tuple(color)] = int(v) + idx += 1 + base = [int(idx/(255*255)), int((idx/255)%255), int(idx%255), 255] + + return(v_id) + + """ + Generates a unique index for every cluster. + """ + def gen_cluster_id(self, G): + self.cluster_dict = {} + base = [0, 0, 0, 255] + idx = 0 + #colors = cm.get_cmap('Wistia', G.num_vertices()*2) + v_id = np.zeros((self.n_c, 4), dtype=np.float32) + for v in range(self.n_c): + color = np.multiply(base, 1/255.0) + v_id[int(v)] = color + self.cluster_dict[tuple(color)] = int(v) + idx += 1 + base = [int(idx/(255*255)), int((idx/255)%255), int(idx%255), 255] + + return(v_id) + + + """ + Generates the bounding box of the radial glyph. + """ + def gen_cluster_coords(self, center, diameter): + radius = diameter/2.0 + top = center[1]+radius + bottom = center[1]-radius + left = center[0]-radius + right = center[0]+radius + + positions = [[right, bottom, center[2]], + [right, top, center[2]], + [left, top, center[2]], + [left, bottom, center[2]]] + + + values = [[1.0, -1.0], + [1.0, 1.0,], + [-1.0, 1.0], + [-1.0, -1.0]] + return positions, values + + + """ + Layout algorithm that expands the cluster based on the location of the of the clusters + """ + def expand_based_on_clusters(self, G, n): + pos = G.vertex_properties["pos"].get_2d_array(range(3)).T + for p in range(pos.shape[0]): + pos[p][0] = pos[p][0] + self.clusters["a_position"][G.vertex_properties["clusters"][G.vertex(p)]][0] + pos[p][1] = pos[p][1] + self.clusters["a_position"][G.vertex_properties["clusters"][G.vertex(p)]][1] + pos[p][2] = pos[p][2] + self.clusters["a_position"][G.vertex_properties["clusters"][G.vertex(p)]][2] + G.vertex_properties["pos"] = G.new_vertex_property("vector", vals = pos) + for i in range(n): + index = 4*i + #generate the vertex filter for this cluster + num_v_in_cluster = len(np.argwhere(self.labels == i)) + vfilt = np.zeros([G.num_vertices(), 1], dtype="bool") + vfilt[np.argwhere(self.labels == i)] = 1 + vfilt_prop = G.new_vertex_property("bool", vals = vfilt) + G.set_vertex_filter(vfilt_prop) + + #get the filtered properties + g = nwt.gt.Graph(G, prune=True, directed=False) + positions = g.vertex_properties["pos"].get_2d_array(range(3)).T + position = np.sum(positions, 0)/num_v_in_cluster + p, v = self.gen_cluster_coords(position, np.sum(g.vertex_properties['degree'].get_array())) + self.clusters['a_position'][index:index+4] = np.asarray(p, dtype=np.float32) + self.clusters['a_value'][index:index+4] = np.asarray(v, dtype=np.float32) + G.clear_filters() + self.cluster_pos[i] = position + color = G.vertex_properties["RGBA"].get_2d_array(range(4)).T + size = nwt.Network.map_vertices_to_range(G, [30*self.pixel_scale, 8*self.pixel_scale], 'degree').get_array() + + position = G.vertex_properties["pos"].get_2d_array(range(3)).T + #for p in range(position.shape[0]): + # position[p][0] = position[p][0] + self.clusters["a_position"][G.vertex_properties["clusters"][G.vertex(p)]][0] + # position[p][1] = position[p][1] + self.clusters["a_position"][G.vertex_properties["clusters"][G.vertex(p)]][1] + # position[p][2] = position[p][2] + self.clusters["a_position"][G.vertex_properties["clusters"][G.vertex(p)]][2] + #G.vertex_properties["pos"] = G.new_vertex_property("vector", vals = position) + edges = G.get_edges(); + edges = edges[:, 0:2] + #width = nwt.Network.map_edges_to_range(G, [1*self.pixel_scale, 5*self.pixel_scale], 'volume').get_array() + #ecolor = G.edge_properties["RGBA"].get_2d_array(range(4)).T + + self.data = np.zeros(G.num_vertices(), dtype=[('a_position', np.float32, 3), + ('a_fg_color', np.float32, 4), + ('a_bg_color', np.float32, 4), + ('a_size', np.float32, 1), + ('a_linewidth', np.float32, 1), + ('a_unique_id', np.float32, 4), + ('a_selection', np.float32, 1), + ]) + + #self.edges = edges.astype(np.uint32) + self.data['a_position'] = position + #fg color is the color of the ring + self.data['a_fg_color'] = 0, 0, 0, 1 + self.data['a_bg_color'] = color + self.data['a_size'] = size + self.data['a_linewidth'] = 4.*self.pixel_scale + self.data['a_unique_id'] = self.gen_vertex_id(G) + self.data['a_selection'] = G.vertex_properties["selection"].get_array() + #self.data['a_graph_size'] = [bbu-bbl] + + self.program['u_graph_size'] = [bbu-bbl] + self.vbo = gloo.VertexBuffer(self.data) + self.gen_line_vbo(G) + if(self.subgraphs): + self.vbo_s = gloo.VertexBuffer(self.clusters) + self.index_s = gloo.IndexBuffer(self.edges_s) + #self.index = gloo.IndexBuffer(self.edges) + self.program_e.bind(self.vbo_line) + self.program.bind(self.vbo) + if(self.subgraphs): + #self.program_e_s.bind(self.vbo_s) + self.program_s.bind(self.vbo_s) + print(self.view) + self.update() + + + """ + Function that generates the clusters for an unclustered graph + These are to be represented by the arcs + """ + def gen_clusters(self, G, bbl, bbu, n_c = None, edge_metric = 'volume', vertex_metric = 'degree'): + + #Generate the clusters + self.labels = nwt.Network.spectral_clustering(G,'length', n_clusters = n_c) + #self.labels = nwt.Network.spectral_clustering(G,'length') + + #Add clusters as a vertex property + G.vertex_properties["clusters"] = G.new_vertex_property("int", vals=self.labels) + num_clusters = len(np.unique(self.labels)) + self.n_c = n_c + + #add colormap + G.vertex_properties["RGBA"] = nwt.Network.map_property_to_color(G, G.vertex_properties["clusters"]) + + #generate an empty property set for the clusters. + self.clusters = np.zeros(num_clusters*4, dtype=[('a_position', np.float32, 3), + ('a_value', np.float32, 2), + ('a_bg_color', np.float32, 4), + ('a_cluster_color', np.float32, 4), + ('a_arc_length', np.float32, 1), + ('a_outer_arc_length', np.float32, 4), + ('a_unique_id', np.float32, 4), + ]) + self.edges_s = np.random.randint(size=(num_clusters*2, 3), low=0, + high=4).astype(np.uint32) + #fill the foreground color as halo + #self.clusters['a_fg_color'] = 1., 1., 1., 0.0 + #self.clusters['a_linewidth'] = 4.*self.pixel_scale + + G.vertex_properties["pos"] = nwt.gt.sfdp_layout(G, groups = G.vertex_properties["clusters"], pos = G.vertex_properties["pos"]) + temp = []; + temp_pos = []; + #Find the global total of the metric. + global_metric = np.sum(G.edge_properties[edge_metric].get_array(), 0) + unique_color = self.gen_cluster_id(G) + + #generate the property values for every cluster + for i in range(num_clusters): + idx = 4*i + #generate the vertex filter for this cluster + num_v_in_cluster = len(np.argwhere(self.labels == i)) + vfilt = np.zeros([G.num_vertices(), 1], dtype="bool") + vfilt[np.argwhere(self.labels == i)] = 1 + vfilt_prop = G.new_vertex_property("bool", vals = vfilt) + G.set_vertex_filter(vfilt_prop) + + #get the filtered properties + g = nwt.gt.Graph(G, prune=True, directed=False) + positions = g.vertex_properties["pos"].get_2d_array(range(3)).T + position = np.sum(positions, 0)/num_v_in_cluster + + #calculate the arclength for the global statistic + arc_length = np.sum(g.edge_properties[edge_metric].get_array(), 0)/global_metric*np.pi*2 + arc_length_vertex = np.ones((4,1), dtype = np.float32) + array = g.vertex_properties[vertex_metric].get_array() + + #calculate metric distribution and turn it into arc_lengths + t_vertex_metric = np.sum(array) + arc_length_vertex[0] = np.sum(array < 2)/t_vertex_metric + arc_length_vertex[1] = np.sum(array == 2)/t_vertex_metric + arc_length_vertex[2] = np.sum(array == 3)/t_vertex_metric + arc_length_vertex[3] = np.sum(array > 3)/t_vertex_metric + + #arc_length_vertex = np.asarray(arc_length_vertex, dtype = np.float32) + #arc_length_vertex = (max(arc_length_vertex) - min(arc_length_vertex)) \ + #* (arc_length_vertex- min(arc_length_vertex)) + for j in range(len(arc_length_vertex)): + if j != 0: + arc_length_vertex[j] += arc_length_vertex[j-1] + print("arc_length before ", arc_length_vertex, " and sum to ", sum(arc_length_vertex)) + arc_length_vertex = np.asarray(arc_length_vertex, dtype = np.float32) + arc_length_vertex = (np.pi - -np.pi)/(max(arc_length_vertex) - min(arc_length_vertex)) \ + * (arc_length_vertex- min(arc_length_vertex)) + (-np.pi) + print(arc_length_vertex) + #print(arc_length) + + + temp_pos.append(position) + + #generate the color for every vertex, + #since all vertices belong to the same cluster we can check only + #one vertex for the cluster color. + + self.clusters['a_cluster_color'][idx:idx+4] = g.vertex_properties["RGBA"][g.vertex(0)] + self.clusters['a_bg_color'][idx:idx+4] = [0.1, 0.1, 0.1, 1.0] + self.clusters['a_unique_id'][idx:idx+4] = unique_color[i] + + #The arc-length representing one global metric. + self.clusters['a_arc_length'][idx:idx+4] = arc_length + self.clusters['a_outer_arc_length'][idx:idx+4] = arc_length_vertex[:].T + + temp.append(np.sum(g.vertex_properties['degree'].get_array())) + G.clear_filters() + print(self.clusters['a_outer_arc_length']) + maximum = max(temp) + minimum = min(temp) + if len(temp) > 1: + temp = ((temp-minimum)/(maximum-minimum)*(60*self.pixel_scale)+20*self.pixel_scale) + else: + temp = [60*self.pixel_scale] + for i in range(num_clusters): + index = i*4 + index_t = i*2 + p, v = self.gen_cluster_coords(temp_pos[i], temp[i]*2.0) + self.clusters['a_position'][index:index+4] = np.asarray(p, dtype=np.float32) + self.clusters['a_value'][index:index+4] = np.asarray(v, dtype=np.float32) + + self.edges_s[index_t] = [index, index+1, index+2] + self.edges_s[index_t+1] = [index, index+2, index+3] + #self.edges_s[i][0:4] = np.asarray(range(index, index+4), dtype=np.uint32) + #self.edges_s[i] + self.cluster_pos = temp_pos + #self.expand_based_on_clusters(G, self.n_c) + G.clear_filters() +# self.edges_s[1][0:4] = np.asarray(range(0, 0+4), dtype=np.uint32) +# self.edges_s[1][4] = 0 +# self.edges_s[1][5] = 0+2 +# +# self.edges_s[0][0:4] = np.asarray(range(index, index+4), dtype=np.uint32) +# self.edges_s[0][4] = index +# self.edges_s[0][5] = index+2 + #self.clusters['a_size'] = temp + self.gen_cluster_line_vbo(G) + self.program_s['u_graph_size'] = [bbu-bbl] + #if len(temp_e) > 0: + # self.edges_s = np.unique(np.asarray(temp_e, np.uint32), axis=0) + #else: + # self.edges_s = [] + #print(self.edges_s) + + + """ + Function that expands that generates the layout and updates the buffer + + """ + def expand_clusters(self, G, n_c): + self.expand_based_on_clusters(G, n_c) + self.gen_cluster_line_vbo(G) + if(self.subgraphs): + self.vbo_s = gloo.VertexBuffer(self.clusters) + self.index_s = gloo.IndexBuffer(self.edges_s) + self.program_e.bind(self.vbo_line) + self.program.bind(self.vbo) + if(self.subgraphs): + self.program_s.bind(self.vbo_s) + print(self.view) + self.update() + + """ + Loads the data G and generates all the buffers necessary as well as performs + spectral clustering on the graph passed if the subgraph is set to true. + """ + def set_data(self, G, bbl, bbu, subgraph=True): + print("Setting data") + self.G = G + self.bbl = bbl + self.bbu = bbu + clear(color=True, depth=True) + self.subgraphs = True + self.current_color = "clusters" + if(subgraph==True): + self.gen_clusters(G, bbl, bbu, n_c=19) + + #color based on clusters + self.color_edges(G) + + color = G.vertex_properties["RGBA"].get_2d_array(range(4)).T + size = nwt.Network.map_vertices_to_range(G, [30*self.pixel_scale, 8*self.pixel_scale], 'degree').get_array() + + position = G.vertex_properties["pos"].get_2d_array(range(3)).T + #for p in range(position.shape[0]): + # position[p][0] = position[p][0] + self.clusters["a_position"][G.vertex_properties["clusters"][G.vertex(p)]][0] + # position[p][1] = position[p][1] + self.clusters["a_position"][G.vertex_properties["clusters"][G.vertex(p)]][1] + # position[p][2] = position[p][2] + self.clusters["a_position"][G.vertex_properties["clusters"][G.vertex(p)]][2] + #G.vertex_properties["pos"] = G.new_vertex_property("vector", vals = position) + edges = G.get_edges(); + edges = edges[:, 0:2] + #width = nwt.Network.map_edges_to_range(G, [1*self.pixel_scale, 5*self.pixel_scale], 'volume').get_array() + #ecolor = G.edge_properties["RGBA"].get_2d_array(range(4)).T + + self.data = np.zeros(G.num_vertices(), dtype=[('a_position', np.float32, 3), + ('a_fg_color', np.float32, 4), + ('a_bg_color', np.float32, 4), + ('a_size', np.float32, 1), + ('a_linewidth', np.float32, 1), + ('a_unique_id', np.float32, 4), + ('a_selection', np.float32, 1), + ]) + + #self.edges = edges.astype(np.uint32) + self.data['a_position'] = position + #fg color is the color of the ring + self.data['a_fg_color'] = 0, 0, 0, 1 + self.data['a_bg_color'] = color + self.data['a_size'] = size + self.data['a_linewidth'] = 4.*self.pixel_scale + self.data['a_unique_id'] = self.gen_vertex_id(G) + self.data['a_selection'] = G.vertex_properties["selection"].get_array() + #self.data['a_graph_size'] = [bbu-bbl] + + self.program['u_graph_size'] = [bbu-bbl] + + self.vbo = gloo.VertexBuffer(self.data) + self.gen_line_vbo(G) + #self.gen_cylinder_vbo(G) + if(self.subgraphs): + self.vbo_s = gloo.VertexBuffer(self.clusters) + self.index_s = gloo.IndexBuffer(self.edges_s) + #self.index = gloo.IndexBuffer(self.edges) + self.program_e.bind(self.vbo_line) + self.program.bind(self.vbo) + if(self.subgraphs): + #self.program_e_s.bind(self.vbo_s) + self.program_s.bind(self.vbo_s) + print(self.view) + self.update() + + """ + Function that changes and redraws the buffer during a resize event. + """ + def on_resize(self, event): + set_viewport(0, 0, *event.physical_size) + self.fbo = gloo.FrameBuffer(color=gloo.RenderBuffer(self.size[::-1]), depth=gloo.RenderBuffer(self.size[::-1])) + + + """ + Overloaded function that is called during every self.update() call + """ + def on_draw(self, event): + clear(color='white', depth=True) + self.program_e.draw('triangles', indices=self.index) + self.program.draw('points') + #self.program_e.draw('lines') + if(self.subgraphs): + self.program_e_s.draw('triangles', indices=self.index_clusters_s) + self.program_s.draw('triangles', indices=self.index_s) + + """ + Function performed during a mouse click (either right or left) + gets the unique id of the object drawn underneath the cursor + handles the cases depending on whether the click happened to a cluster + or a vertex. Edges are not interactable (yet) + """ + def get_clicked_id(self, event, clusters = False): + #Get the framebuffer coordinates of the click + coord = self.transforms.get_transform('canvas', 'framebuffer').map(event.pos) + + #get the framebuffer where each element is rendered as a unique color + size = self.size; + self.fbo = gloo.FrameBuffer(color=gloo.RenderBuffer(size[::-1]), depth=gloo.RenderBuffer(size[::-1])) + buff = gloo.read_pixels((0,0,self.physical_size[0], self.physical_size[1])) + #imsave("test_ori.png", buff) + self.fbo.activate() + if clusters == False: + self.program['u_picking'] = True + clear(color='white', depth=True) + self.program.draw('points') + buff = gloo.read_pixels((0,0,self.physical_size[0], self.physical_size[1])) + #imsave("test.png", buff) + + #return to the original state + self.fbo.deactivate() + self.program['u_picking'] = False + else: + self.program_s['u_picking'] = True + clear(color='white', depth=True) + self.program_s.draw('triangles', indices=self.index_s) + buff = gloo.read_pixels((0,0,self.physical_size[0], self.physical_size[1])) + #imsave("test.png", buff) + + #return to the original state + self.fbo.deactivate() + self.program_s['u_picking'] = False + + #print(buff[self.physical_size[1]-int(coord[1]), int(coord[0])]) + + #Get the color under the click. + #Keep in mind that the buff is y, x + #And 0,0 is in the top RIGHT corner. + #IGNORE THE DOCUMENTATION + color = np.multiply(buff[self.physical_size[1]-int(coord[1]), int(coord[0])], 1/255.0) + #if (tuple(color) not in self.color_dict): + # print("clicked on nothing") + #else: + # print(self.color_dict[tuple(color)]) + + #reset the original buffer + self.update() + + #Return the element under the click. + if clusters == False: + if(tuple(color) not in self.color_dict): + return None + else: + return self.color_dict[tuple(color)] + else: + if(tuple(color) not in self.cluster_dict): + return None + else: + return self.cluster_dict[tuple(color)] + + + """ + Top level handle-mouse presee event for either left or right click + """ + def on_mouse_press(self, event): + + def update_view(): + self.location = event.pos + self.program['u_view'] = self.view + self.program_e['u_view'] = self.view + self.program_s['u_view'] = self.view + self.program_e_s['u_view'] = self.view + self.down = True + + +# if(event.button == 2): +## menu = QtWidgets.QMenu(self.parent) +## NS = menu.addAction('Node Size') +## NC = menu.addAction('Node Color') +## action = menu.exec_(self.parent.globalPos()) +## if action == NS: +# print("right_click") +# #if menu.exec_(event.globalPos()): +# # print(item.text()) + if(event.button == 1): + if(self.view[0][0] > 0.0010): + c_id = self.get_clicked_id(event) + if(c_id != None): + self.original_point = G.vertex_properties["pos"][G.vertex(c_id)] + self.location = event.pos + self.moving = True + self.down = True + self.c_id = [c_id] + else: + update_view() + #print("Clicked on:", event.pos) + else: + c_id = self.get_clicked_id(event, True) + print(c_id) + if(c_id != None): + self.original_point = self.cluster_pos[c_id] + self.location = event.pos + self.moving = True + self.down = True + self.c_id = [c_id] + self.moving_cluster = True + else: + update_view() + + + """ + Handles the double click event that it responsible for path selection. + Generates paths our of consecutive paths out of the selected vertices. + """ + def on_mouse_double_click(self, event): + def update_vbo(self): + self.vbo = gloo.VertexBuffer(self.data) + self.program.bind(self.vbo) + self.update() + + def add_to_path(self, source, target): + vl, el = nwt.gt.graph_tool.topology.shortest_path(G, G.vertex(source), G.vertex(target), weights=G.edge_properties["av_radius"]) + for v in vl: + if(int(v) not in self.path): + G.vertex_properties["selection"][v] = 2.0 + self.data['a_selection'][int(v)] = 2.0 + if(int(v) not in self.full_path): + self.full_path.append(int(v)) + + if (event.button == 1): + if(self.view[0][0] > 0.0010): + c_id = self.get_clicked_id(event) + if(c_id != None): + #check whether this is the first node to be selected + if(self.pathing == False): + #if it is, select that node and turn the pathing variable on. + G.vertex_properties["selection"][G.vertex(c_id)] = 1.0 + self.pathing = True + if(c_id not in self.path): + self.path.append(c_id) + self.data['a_selection'][c_id] = 1.0 + update_vbo(self) + print("I turned on the first node") + else: + if(G.vertex_properties["selection"][G.vertex(c_id)] == 1.0): + G.vertex_properties["selection"][G.vertex(c_id)] = 0.0 + self.path.remove(c_id) + self.data['a_selection'][c_id] = 0.0 + update_vbo(self) + print("I turned off a node") + elif(G.vertex_properties["selection"][G.vertex(c_id)] == 0.0): + G.vertex_properties["selection"][G.vertex(c_id)] = 1.0 + if(c_id not in self.path): + self.path.append(c_id) + self.data['a_selection'][c_id] = 1.0 + update_vbo(self) + print("I turned on a node") + if(len(self.path) >= 2): + for i in range(len(self.path)-1): + add_to_path(self, self.path[i], self.path[i+1]) + update_vbo(self) + #THIS IS WHERE I LEFT IT OFF. + if(np.sum(G.vertex_properties["selection"].get_array()) == 0): + self.pathing = False + + + +# elif(np.sum(self.G.vertex_properties["selection"].get_array()) : +# self.G.vertex_properties["selection"][self.G.vertex(c_id)] == False + print("clicked on: ", c_id, " ", self.path) + + """ + Resets the variables that are used during the pressdown and move events + """ + def on_mouse_release(self, event): + self.down = False + self.moving = False + self.moving_cluster = False + self.c_id = [] + #self.location = event.pos + #print("Clicked off:", event.pos) + + """ + used during the drag evern to update the position of the clusters + """ + def update_cluster_position(self, G, pos, offset, c_id): + v_pos = G.vertex_properties["pos"].get_2d_array(range(3)).T + vertices = np.argwhere(self.labels == c_id) + for v in range(vertices.shape[0]): + idx = vertices[v][0] + v_pos[idx][0] = v_pos[idx][0] + offset[0] + v_pos[idx][1] = v_pos[idx][1] + offset[1] + v_pos[idx][2] = v_pos[idx][2] + offset[2] + self.data['a_position'][idx] = np.asarray([v_pos[idx][0], v_pos[idx][1], v_pos[idx][2]], dtype = np.float32) + #update the edge data by finding all edges connected to the vertex + vtx = self.G.vertex(idx) + for e in vtx.all_edges(): + d = np.subtract(G.vertex_properties["pos"][e.source()], G.vertex_properties["pos"][e.target()]) + d_norm = d[0:2] + d_norm = d_norm / np.sqrt(np.power(d_norm[0],2) + np.power(d_norm[1],2)) + norm = np.zeros((2,), dtype=np.float32) + norm[0] = d_norm[1] + norm[1] = d_norm[0]*-1 + if (int(e.source()), int(e.target())) in self.edge_dict.keys(): + index = int(self.edge_dict[int(e.source()), int(e.target())]) + if vtx == int(e.source()): + self.line_data['a_position'][index*4] = v_pos[idx] + self.line_data['a_position'][index*4+2] = v_pos[idx] + self.line_data['a_normal'][index*4] = norm + self.line_data['a_normal'][index*4+2] = -norm + self.line_data['a_normal'][index*4+1] = norm + self.line_data['a_normal'][index*4+3] = -norm + elif vtx == int(e.target()): + self.line_data['a_position'][index*4+1] = v_pos[idx] + self.line_data['a_position'][index*4+3] = v_pos[idx] + self.line_data['a_normal'][index*4] = norm + self.line_data['a_normal'][index*4+2] = -norm + self.line_data['a_normal'][index*4+1] = norm + self.line_data['a_normal'][index*4+3] = -norm + else: + index = int(self.edge_dict[int(e.target()), int(e.source())]) + if vtx == int(e.target()): + self.line_data['a_position'][index*4] = v_pos[idx] + self.line_data['a_position'][index*4+2] = v_pos[idx] + self.line_data['a_normal'][index*4] = norm + self.line_data['a_normal'][index*4+2] = -norm + self.line_data['a_normal'][index*4+1] = norm + self.line_data['a_normal'][index*4+3] = -norm + elif vtx == int(e.source()): + self.line_data['a_position'][index*4+1] = v_pos[idx] + self.line_data['a_position'][index*4+3] = v_pos[idx] + self.line_data['a_normal'][index*4] = norm + self.line_data['a_normal'][index*4+2] = -norm + self.line_data['a_normal'][index*4+1] = norm + self.line_data['a_normal'][index*4+3] = -norm + + + G.vertex_properties["pos"] = G.new_vertex_property("vector", vals = v_pos) + index = 4*c_id + #generate the vertex filter for this cluster + vfilt = np.zeros([G.num_vertices(), 1], dtype="bool") + vfilt[np.argwhere(self.labels == c_id)] = 1 + vfilt_prop = G.new_vertex_property("bool", vals = vfilt) + G.set_vertex_filter(vfilt_prop) + + #get the filtered properties + g = nwt.gt.Graph(G, prune=True, directed=False) + p, v = self.gen_cluster_coords(pos, np.sum(g.vertex_properties['degree'].get_array())) + self.clusters['a_position'][index:index+4] = np.asarray(p, dtype=np.float32) + self.clusters['a_value'][index:index+4] = np.asarray(v, dtype=np.float32) + G.clear_filters() + self.cluster_pos[c_id] = pos + self.original_point = pos + + + """ + function that handles the mouse move event in a way that depends on a set + of variables: state of the mouse button, the type of object selected and + the number of objects. + """ + def on_mouse_move(self, event): + if(self.down == True): + if(self.moving == True and self.moving_cluster == False): + if(len(self.c_id) < 2): + #Project into GLSpace and get before and after move coordinates + coord = self.transforms.get_transform('canvas', 'render').map(self.location)[:2] + coord2 = self.transforms.get_transform('canvas', 'render').map(event.pos)[:2] + cur_pos = G.vertex_properties["pos"][G.vertex(self.c_id[0])] + #print(cur_pos, " Before") + + #Adjust the position of the node based on the current view matrix. + cur_pos[0] = cur_pos[0] - (coord[0]-coord2[0])/self.view[0][0] + cur_pos[1] = cur_pos[1] - (coord[1]-coord2[1])/self.view[0][0] + + #print(cur_pos, " After") + #Upload the changed data. + G.vertex_properties["pos"][G.vertex(self.c_id[0])] = cur_pos + self.data['a_position'][self.c_id[0]] = cur_pos + + #update the edge data by finding all edges connected to the vertex + v = self.G.vertex(self.c_id[0]) + for e in v.all_edges(): + d = np.subtract(G.vertex_properties["pos"][e.source()], G.vertex_properties["pos"][e.target()]) + d_norm = d[0:2] + d_norm = d_norm / np.sqrt(np.power(d_norm[0],2) + np.power(d_norm[1],2)) + norm = np.zeros((2,), dtype=np.float32) + norm[0] = d_norm[1] + norm[1] = d_norm[0]*-1 + if (int(e.source()), int(e.target())) in self.edge_dict.keys(): + idx = int(self.edge_dict[int(e.source()), int(e.target())]) + if self.c_id[0] == int(e.source()): + self.line_data['a_position'][idx*4] = cur_pos + self.line_data['a_position'][idx*4+2] = cur_pos + self.line_data['a_normal'][idx*4] = norm + self.line_data['a_normal'][idx*4+2] = -norm + self.line_data['a_normal'][idx*4+1] = norm + self.line_data['a_normal'][idx*4+3] = -norm + elif self.c_id[0] == int(e.target()): + self.line_data['a_position'][idx*4+1] = cur_pos + self.line_data['a_position'][idx*4+3] = cur_pos + self.line_data['a_normal'][idx*4] = norm + self.line_data['a_normal'][idx*4+2] = -norm + self.line_data['a_normal'][idx*4+1] = norm + self.line_data['a_normal'][idx*4+3] = -norm + else: + idx = int(self.edge_dict[int(e.target()), int(e.source())]) + if self.c_id[0] == int(e.target()): + self.line_data['a_position'][idx*4] = cur_pos + self.line_data['a_position'][idx*4+2] = cur_pos + self.line_data['a_normal'][idx*4] = norm + self.line_data['a_normal'][idx*4+2] = -norm + self.line_data['a_normal'][idx*4+1] = norm + self.line_data['a_normal'][idx*4+3] = -norm + elif self.c_id[0] == int(e.source()): + self.line_data['a_position'][idx*4+1] = cur_pos + self.line_data['a_position'][idx*4+3] = cur_pos + self.line_data['a_normal'][idx*4] = norm + self.line_data['a_normal'][idx*4+2] = -norm + self.line_data['a_normal'][idx*4+1] = norm + self.line_data['a_normal'][idx*4+3] = -norm + #self.line_data['a_position'][self.c_id[0]] = + self.vbo = gloo.VertexBuffer(self.data) + self.vbo_line = gloo.VertexBuffer(self.line_data) + #Bind the buffer and redraw. + self.program.bind(self.vbo) + self.program_e.bind(self.vbo_line) + #self.program.draw('points') + self.location = event.pos + self.update() + elif(self.moving == True and self.moving_cluster == True): + if(len(self.c_id) < 2): + #Project into GLSpace and get before and after move coordinates + coord = self.transforms.get_transform('canvas', 'render').map(self.location)[:2] + coord2 = self.transforms.get_transform('canvas', 'render').map(event.pos)[:2] + cur_pos = np.zeros(self.cluster_pos[self.c_id[0]].shape, dtype = np.float32) + offset = np.zeros(self.cluster_pos[self.c_id[0]].shape, dtype = np.float32) + cur_pos[0] = self.cluster_pos[self.c_id[0]][0] + cur_pos[1] = self.cluster_pos[self.c_id[0]][1] + cur_pos[2] = self.cluster_pos[self.c_id[0]][2] + offset[0] = self.cluster_pos[self.c_id[0]][0] + offset[1] = self.cluster_pos[self.c_id[0]][1] + offset[2] = self.cluster_pos[self.c_id[0]][2] +# ofset = self.cluster_pos[self.c_id[0]] + + #Adjust the position of the node based on the current view matrix. + offset[0] = self.original_point[0] - cur_pos[0] - (coord[0]-coord2[0])/self.view[0][0] + offset[1] = self.original_point[1] - cur_pos[1] - (coord[1]-coord2[1])/self.view[0][0] + cur_pos[0] = cur_pos[0] - (coord[0]-coord2[0])/self.view[0][0] + cur_pos[1] = cur_pos[1] - (coord[1]-coord2[1])/self.view[0][0] + + self.update_cluster_position(G, cur_pos, offset, self.c_id[0]) + #self.original_point = cur_pos + self.vbo = gloo.VertexBuffer(self.data) + self.vbo_line = gloo.VertexBuffer(self.line_data) + #Bind the buffer and redraw. + self.program.bind(self.vbo) + self.program_e.bind(self.vbo_line) + #self.program.draw('points') + self.location = event.pos + if(self.subgraphs): + self.vbo_s = gloo.VertexBuffer(self.clusters) + self.program_s.bind(self.vbo_s) + self.update_cluster_line_vbo() + self.update() + + + else: + #print("Mouse at:", event.pos) + #new_model = np.eye(4, dtype=np.float32) + coord = self.transforms.get_transform('canvas', 'render').map(self.location)[:2] + coord2 = self.transforms.get_transform('canvas', 'render').map(event.pos)[:2] + self.translate[0] += (coord[0]-coord2[0])/self.view[0][0] + self.translate[1] += (coord[1]-coord2[1])/self.view[1][1] + #self.view[3][0] = self.view[3][0]-(self.location[0]-event.pos[0])/10000.0 + #self.view[3][1] = self.view[3][1]+(self.location[1]-event.pos[1])/10000.0 + + self.view = np.matmul(translate((self.translate[0], self.translate[1], 0)), scale((self.scale[0], self.scale[1], 0))) + + self.program['u_view'] = self.view + self.program_e['u_view'] = self.view + self.program_s['u_view'] = self.view + self.program_e_s['u_view'] = self.view + self.location = event.pos + self.update() + + """ + Handles the mouse wheel zoom event. + """ + def on_mouse_wheel(self, event): + + #print(self.view) + #TO_DO IMPLEMENT ZOOM TO CURSOR + #self.view[3][0] = self.view[3][0]-event.pos[0]/10000.0 + #self.view[3][1] = self.view[3][1]-event.pos[1]/10000.0 + #print(self.scale[0] , self.scale[0]*event.delta[1]*0.05) + self.scale[0] = self.scale[0] + self.scale[0]*event.delta[1]*0.05 + self.scale[1] = self.scale[1] + self.scale[1]*event.delta[1]*0.05 + + self.view = np.matmul(translate((self.translate[0], self.translate[1], 0)), + scale((self.scale[0], self.scale[1], 0))) + + #self.view[0][0] = self.view[0][0]+self.view[0][0]*event.delta[1]*0.05 + #self.view[1][1] = self.view[1][1]+self.view[1][1]*event.delta[1]*0.05 + #print(self.view[0][0], " ",self.view[1][1]) + #print(self.view) + self.program['u_view'] = self.view + self.program_e['u_view'] = self.view + self.program_s['u_view'] = self.view + self.program_e_s['u_view'] = self.view + #print(event.delta[1]) + self.update() diff --git a/GraphWidget.py b/GraphWidget.py new file mode 100644 index 0000000..36ceed3 --- /dev/null +++ b/GraphWidget.py @@ -0,0 +1,184 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Created on Mon Aug 5 15:42:19 2019 + +@author: pavel +""" + +from GraphCanvas import GraphCanvas +from pyqtgraph.Qt import QtCore, QtGui, QtWidgets +import network_dep as nwt + + +""" +Initializes the entire QTlayout and sets the mouse press events. +These are connected to the slots such that each is processes by this class +and the class it wraps. +""" +class GraphWidget(QtGui.QWidget): + def __init__(self): + super(GraphWidget, self).__init__() + box = QtGui.QVBoxLayout(self) + self.resize(500,500) + self.setLayout(box) + self.canvas = GraphCanvas() + #self.canvas.create_native() + box.addWidget(self.canvas.native) + self.color = True + self.use_3D = False + self.camera = [0,0,0] + + self.canvas.events.mouse_press.connect(self.on_mouse_press) + self.canvas.events.mouse_release.connect(self.on_mouse_release) + self.canvas.events.mouse_move.connect(self.on_mouse_move) + self.canvas.events.mouse_double_click.connect(self.on_mouse_double_click) + + #self.show() + + """ + Handles the right click mouse event at the QWidget level. + Depending on where the mouse button in clicked and the current zoom level, + the function gets the unique id of the node, cluster or background under the + cursor and opens a context menu based on the ID. + + The context menu level actions are also coded in this section of the code. + """ + def on_mouse_press(self, event): + if event.button == 2: + if self.canvas.view[0][0] >= 0.0010: + c_id = self.canvas.get_clicked_id(event) + else: + c_id = self.canvas.get_clicked_id(event, clusters=True) + #right click + if(c_id == None): + menu = QtWidgets.QMenu(self) + NS = menu.addAction('Node Size') + + #tmp = menu.addAction('temp') + tmp = menu.addMenu('Node Color') + vertex_props = self.canvas.G.vertex_properties.keys() + vertex_actions = [] + for i in range(len(vertex_props)): + if(vertex_props[i] != "map" or vertex_props[i] != "RGBA"): + vertex_actions.append(tmp.addAction(vertex_props[i])) + #tmp.addAction("Cluster") + #NC = menu.addAction('Node Color') + #EXP = menu.addAction('Export VTK') + #EXP_adv = menu.addAction('Export VTK Advanced') + NL = menu.addAction('New Layout') + action = menu.exec_(event.native.globalPos()) + print(action) + for i in range(len(vertex_actions)): + if action == vertex_actions[i]: + if vertex_props[i] == "clusters": + self.canvas.color_vertices(self.canvas.G, vertex_props[i], dtype=True) + else: + self.canvas.color_vertices(self.canvas.G, vertex_props[i]) + print(vertex_props[i]) + if action == NS: + if self.use_3D == False: + self.use_3D = True + else: + self.use_3D = False + self.canvas.size_vertices(self.canvas.G, 'degree' ) + self.canvas.color_vertices(self.canvas.G) + #if action == NC: + # self.canvas.color_vertices(self.canvas.G, not self.color) + # self.color = not self.color +# if action == EXP: +# nwt.Network.write_vtk(self.canvas.G, "./nwt.vtk") +# if action == EXP_adv: +# nwt.Network.write_vtk(self.canvas.G, "./nwt_median_binned.vtk") +# nwt.Network.write_vtk(self.canvas.G, "./nwt_median_non_binned.vtk", binning=False) +# nwt.Network.write_vtk(self.canvas.G, "./nwt_points_wise_binned.vtk", camera=self.camera, binning = True) +# nwt.Network.write_vtk(self.canvas.G, "./nwt_points_wise_non_binned.vtk", camera=self.camera, binning = False) + if action == tmp: + print("Stuff") + #self.cb = QtWidgets.QComboBox() + #vertex_props = self.canvas.G.vertex_properties.keys() + #for i in range(len(vertex_props)): + # self.cb.addItem(vertex_props[i]) + #self.cb.currentIndexChanged.connect(self.selection_change) + #self.cb.show() + if action == NL: + self.canvas.G = nwt.Network.gen_new_fd_layout(self.canvas.G) + self.canvas.gen_vertex_vbo(self.canvas.G) + #self.canvas.set_data(self.canvas.G, self.canvas.bbl, self.canvas.bbu) + self.canvas.expand_clusters(self.canvas.G, self.canvas.n_c) + + #self.canvas.size_vertices(self.canvas.G, 'degree_volume') + else: + if self.canvas.view[0][0] >= 0.0010: + menu = QtWidgets.QMenu(self) + NS = menu.addAction('Display Node Info') + action = menu.exec_(event.native.globalPos()) + if action == NS: + print("Display Node Info") + else: + menu = QtWidgets.QMenu(self) + MnOC = menu.addAction('Minimize Other Clusters') + RA = menu.addAction('Reset Alpha') + action = menu.exec_(event.native.globalPos()) + if action == MnOC: + new_Graph = self.canvas.focus_on_cluster(self.canvas.G, c_id) + self.test.draw_new(new_Graph) + if action == RA: + self.canvas.reset_alpha(c_id[0]) + + +# def selection_change(self, i): +# self.canvas.size_vertices(self.canvas.G, self.cb.currentText()) + + + """ + Handles the mouse double click event. + Pass-through function. + """ + def on_mouse_double_click(self, event): + #print("in applet") + n = 0 + + """ + Handles the mouse release event. + Pass-through function. + """ + def on_mouse_release(self, event): + n = 1 + + """ + Handles the mouse move event. + Pass-through function. + """ + def on_mouse_move(self, event): + n = 2 + + """ + Internal function that interacts with the tube visualization. + """ + def set_g_view(self, test): + self.test = test + + """ + Function to handle the pyqt Slot that receives the interacted-with signal + from the Tube visualization. The signal sends the camera position in 3D coordinates + every time the camera is moved. + """ + @QtCore.pyqtSlot(float, float, float) + def sigUpdate(self, x, y, z): + self.camera = [x,y,z] + if(self.use_3D == True): + self.camera = [x,y,z] + self.canvas.vertexSizeFromDistance(self.canvas.G, [x,y,z]) + self.canvas.vertexAlphaFromDistance(self.canvas.G, [x,y,z]) + #print("got signal", x, y, z) + + """ + Initialization method that connects the slot to the function that handles + the information being sent. + """ + def connect(self, signal_object): + signal_object.sigUpdate.connect(self.sigUpdate) + + + diff --git a/GuiVisPy_tube.py b/GuiVisPy_tube.py new file mode 100644 index 0000000..96271d7 --- /dev/null +++ b/GuiVisPy_tube.py @@ -0,0 +1,112 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Created on Thu Jan 31 15:29:40 2019 + +@author: pavel +`""" + +from vispy import app + + +import network_dep as nwt + +from pyqtgraph.Qt import QtCore, QtGui, QtWidgets +import pyqtgraph as pg +import numpy as np +import math + +from mpl_toolkits.mplot3d import Axes3D +import matplotlib +import matplotlib.pyplot as plt +matplotlib.use('Qt5Agg') + +from GraphWidget import GraphWidget +from TubeWidget import TubeWidget + + +#set the backend. for different versions of PyQt +app.use_app(backend_name='PyQt5') + +#Define a top level application +appMain = QtGui.QApplication([]) + +##Define a toplevel Widget +top = QtGui.QWidget() +top.resize(900, 900) + + +hist = pg.PlotWidget() +#fibers = FiberView() +fibers = TubeWidget() +fibers.canvas.create_native() +#fibers = gl.GLViewWidget() + +#plt = hist.addPlot() + +layout = QtGui.QGridLayout() +graph = GraphWidget() +graph.canvas.create_native() + +graph.connect(fibers) + +#initialize the layout +layout.addWidget(graph, 0, 0, 2, 3) +layout.addWidget(fibers, 0, 2, 2, 3) +layout.setColumnStretch(0, 2) +layout.setRowStretch(0, 2) +layout.setColumnStretch(2, 1) +layout.setRowStretch(0, 1) + +top.setLayout(layout) +top.show() + + +def draw_histogram(G, value): + vals = G.edge_properties[value].get_array() + y, x = np.histogram(vals,40) + hist.plot(x,y, stepMode=True, fillLevel=0, brush=(0,0,255,150)) + +def load_nwt(filepath): + net = nwt.Network(filepath) + G = net.createFullGraph_gt() + G = net.filterDisconnected(G) + #G = net.filterFullGraph_gt(G, erode=True) + #G = net.filterFullGraph_gt(G, erode=True) + #G = net.gen_new_fd_layout(G) + #G = net.filterBorder(G) + #nwt.Network.saveGraph_gt(nwt, G, "FilteredNetwork_JACK_3.nwt") + color = np.zeros(4, dtype = np.double) + color = [0.0, 1.0, 0.0, 1.0] + G.edge_properties["RGBA"] = G.new_edge_property("vector", val=color) + color = [1.0, 0.0, 0.0, 0.9] + G.vertex_properties["RGBA"] = G.new_vertex_property("vector", val=color) + bbl, bbu = net.aabb() + + + return G, bbl, bbu + +G, bbl, bbu = load_nwt("/home/pavel/Documents/Python/GraphGuiQt/network_4.nwt") +#nwt.Network.write_vtk(G) + +#G, bbl, bbu = load_nwt("/home/pavel/Documents/Python/GraphGuiQt/network_test_251_1kx3_short_NEW.nwt") +#ret = nwt.Network.get_affinity_matrix(G, "length") +node_image = QtGui.QImage("/home/pavel/Documents/Python/GraphGuiQt/node_tex.jpg") +node_tex = QtGui.QBitmap.fromImage(node_image) +print(node_tex.depth()) +center = (bbu-bbl)/2.0 +#fibers.opts['distance'] = 5 +# item = NodeItem(G) +# graph.addItem(item) +draw_histogram(G, "length") +graph.canvas.set_data(G, bbl, bbu) +fibers.canvas.set_data(G, bbl, bbu, 16) +#fibers.draw_all(G, center, fibers, graph, node_tex) +graph.set_g_view(fibers) +## Start Qt event loop unless running in interactive mode. +if __name__ == '__main__': + import sys + if (sys.flags.interactive != 1) or not hasattr(QtCore, 'PYQT_VERSION'): + QtGui.QApplication.instance().exec_() + + diff --git a/TubeCanvas.py b/TubeCanvas.py new file mode 100644 index 0000000..96c476c --- /dev/null +++ b/TubeCanvas.py @@ -0,0 +1,443 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Created on Mon Aug 5 15:56:47 2019 + +@author: pavel +""" + +""" + Class that extends the vispy SceneCanvas to draw 3D tubes +""" + +from vispy import gloo, scene +from vispy.gloo import set_viewport, set_state, clear, set_blend_color +from vispy.util.transforms import perspective, translate, rotate, scale +import vispy.gloo.gl as glcore +from vispy.util.quaternion import Quaternion + +import numpy as np +import math +import network_dep as nwt + + +from tube_shaders import FRAG_SHADER, VERT_SHADER + +class TubeDraw(scene.SceneCanvas): + #sigUpdate = QtCore.pyqtSignal(float, float, float) + + #Initiates the canvas. + def __init__(self, **kwargs): + #Initialte the class by calling the superclass + scene.SceneCanvas.__init__(self, size=(512,512), keys='interactive', **kwargs) + #unfreeze the drawing area to allow for dynamic drawing and interaction + self.unfreeze() + + #generate dummy buffers for the meshes + self.program = gloo.Program(VERT_SHADER, FRAG_SHADER) + self.cylinder_data = np.zeros(5*5, dtype=[('a_position', np.float32, 3), + ('a_normal', np.float32, 3), + ('a_fg_color', np.float32, 4), + #('a_linewidth', np.float32, 1), + ]) + self.triangle_data = np.random.randint(size=(5, 3), low=0, + high=(4-1)).astype(np.uint32) + self.vbo = gloo.VertexBuffer(self.cylinder_data) + self.triangles = gloo.IndexBuffer(self.triangle_data) + self.program.bind(self.vbo) + self.scale = [1,1,1] + self.r1 = np.eye(4, dtype=np.float32) + self.r2 = np.eye(4, dtype=np.float32) + set_viewport(0,0,*self.physical_size) + #set_state(clear_color='white', depth_test=True, blend=True, + # blend_func=('src_alpha', 'one_minus_src_alpha'), depth_func = ('less'), cull_face='back') + set_state(clear_color='white', depth_test=True, blend=True, + blend_func=('src_alpha', 'one_minus_src_alpha'), depth_func = ('lequal')) + #set_blend_color(color='black') + #set_state('translucent') + self.program['u_LightPos'] = [0., 0., -1000.] + #self.camera = self.central_widget.add_view() + #self.camera.camera = 'turntable' + self.down = False + self.camera = np.asarray([0.0, 0.0, 200.0], dtype=np.float32) + self.up = np.asarray([0., 1., 0.], dtype=np.float32) + #self.init_camera = [0.,0.,1000.] + + ##### prototype ##### + #Set the visualization matrices + self.program['u_eye'] = self.camera + self.program['u_up'] = self.up + self.program['u_target'] = np.asarray([0., 0., 0.], dtype=np.float32) + + + + #Load the data necessary to draw all of the microvessels + def set_data(self, G, bbu, bbl, num_sides): + self.G = G + self.num_sides = num_sides + self.bbu = bbu + self.bbl = bbl + bb = nwt.AABB(G).resample_sides(3) + + + #create program + self.gen_cylinder_vbo(self.G, self.num_sides) + self.vbo = gloo.VertexBuffer(self.cylinder_data) + self.triangles = gloo.IndexBuffer(self.triangle_data) + + #self.view = np.eye(4, dtype=np.float32) + self.model = np.eye(4, dtype=np.float32) + self.projection = np.eye(4, dtype=np.float32) + self.projection = perspective(90.0, self.physical_size[0]/self.physical_size[1], 1.0, 1000.0) + #self.projection = perspective(90.0, 1.0, -1.0, 1.0) + self.program['u_model'] = self.model + self.program['u_LightPos'] = [0., 0., 1000.] + #self.program['u_view'] = self.view + self.program['u_projection'] = self.projection + self.program.bind(self.vbo) + + gloo.set_clear_color('white') + self.center = (bbu-bbl)/2.0 + self.translate = [-self.center[0], -self.center[1], -self.center[2]] + + self.bb = np.ones((26, 3), dtype=np.float32) + for i in range(26): + for j in range(3): + self.bb[i,j] = bb[i][j] + self.program['u_bb'] = self.bb + print('bb is ', self.bb) +# for i in range(len(self.translate)): +# self.camera[i] += self.translate[i] + + + ##### prototype ##### + self.camera = self.camera - self.translate + self.program['u_eye'] = self.camera + self.up = np.cross((np.asarray(self.center, dtype=np.float32)-np.asarray(self.camera, dtype=np.float32)), np.asarray(self.up)) + self.program['u_up'] = self.up + self.program['u_target'] = self.translate + + + + + #self.show() + + #Called during resize of the window in order to redraw the same image in the + #larger/smaller area. + def on_resize(self, event): + width, height = event.physical_size + gloo.set_viewport(0, 0, width, height) + print(self.physical_size) + + #overloaded function called during the self.update() call to update the current + #frame using the GLSL frag/vert shaders + def on_draw(self, event): + clear(color='white', depth=True) + gloo.set_clear_color('white') + self.program.draw('triangles', self.triangles) + self.projection = perspective(90.0, self.physical_size[0]/self.physical_size[1], 1.0, 1000.0) + self.program['u_projection'] = self.projection + + #Creates a cylinder around ever segment in the microvascular network. + def gen_cylinder_vbo(self, G, num_sides = 32): + i = 0 + num_pts = 0 + num_tri = 0 + for e in G.edges(): + num_pts += len(self.G.edge_properties["x"][e]) + num_tri += (len(self.G.edge_properties["x"][e])-1)*num_sides*2 + self.cylinder_data = np.zeros(num_pts*num_sides, dtype=[('a_position', np.float32, 3), + ('a_normal', np.float32, 3), + ('a_fg_color', np.float32, 4), + #('a_linewidth', np.float32, 1), + ]) + self.triangle_data = np.random.randint(size=(num_tri, 3), low=0, + high=(G.num_edges()-1)).astype(np.uint32) + index = 0 + t_index = 0 + #for each edge generate a cylinder. + for e in G.edges(): + #print("done") + #for each fiber get all the points and the radii + X = self.G.edge_properties["x"][e] + Y = self.G.edge_properties["y"][e] + Z = self.G.edge_properties["z"][e] + R = self.G.edge_properties["r"][e] + color = G.edge_properties["RGBA"][e] + pts = np.array([X,Y,Z]).T + circle_pts = np.zeros((pts.shape[0], num_sides, 3), dtype = np.float32) + step = 2*np.pi/num_sides +# U = np.zeros(pts.shape, dtype=np.float32) +# V = np.zeros(pts.shape, dtype=np.float32) +# direction = np.zeros(pts.shape, dtype=np.float32) + + #for every point in the edge + for p in range(pts.shape[0]): + #if first point, generate the circles. + #In this case we want to generate a cap if we see the first circle or the last. + if(p == 0): + #get the direction + direction = (pts[p+1] - pts[p]) + #normalize direction + direction = direction/np.sqrt(direction[0]*direction[0] + direction[1]*direction[1] + direction[2]*direction[2]) + #generate a vector to use as an element of the cross product + Y = np.zeros((3,), dtype = np.float32) + Y[0] = 1. + #if direction and Y are parallel, change Y + if(np.cos(np.dot(Y, direction)) < 0.087): + Y[0] = 0.0 + Y[2] = 1.0 + #generate first plane vector + U = np.cross(direction, Y) + U = U/np.sqrt(U[0]*U[0] + U[1]*U[1] + U[2]*U[2]) + #generate second plane vector + V = np.cross(direction, U) + V = V/np.sqrt(V[0]*V[0] + V[1]*V[1] + V[2]*V[2]) + #print(R[p],pts[p]) + #generate circle. + for s in range(num_sides): + circle_pts[p][s][0] = R[p]*np.cos(s*step)*V[0]*0.5 + R[p]*np.sin(s*step)*U[0]*0.5 + circle_pts[p][s][1] = R[p]*np.cos(s*step)*V[1]*0.5 + R[p]*np.sin(s*step)*U[1]*0.5 + circle_pts[p][s][2] = R[p]*np.cos(s*step)*V[2]*0.5 + R[p]*np.sin(s*step)*U[2]*0.5 + #if last point, copy the previous circle. + elif(p == pts.shape[0]-1): + for s in range(num_sides): + circle_pts[p][s] = circle_pts[p-1][s] + for v in range(pts.shape[0]): + circle_pts[v,:,0] += pts[v][0] + circle_pts[v,:,1] += pts[v][1] + circle_pts[v,:,2] += pts[v][2] + #otherwise, rotate the circle + else: + #print(R[p], pts[p]) + #generate a new direction vector. + dir_new = (pts[p+1]-pts[p]) + dir_new = dir_new/np.sqrt(dir_new[0]*dir_new[0] + dir_new[1]*dir_new[1] + dir_new[2]*dir_new[2]) + dir_new2 = (pts[p]-pts[p-1]) + dir_new2 = dir_new2/np.sqrt(dir_new2[0]*dir_new2[0] + dir_new2[1]*dir_new2[1] + dir_new2[2]*dir_new2[2]) + dir_new = dir_new+dir_new2 + dir_new = dir_new/np.sqrt(dir_new[0]*dir_new[0] + dir_new[1]*dir_new[1] + dir_new[2]*dir_new[2]) + #print(dir_new, direction) + #generate the quaternion rotation vector for the shortest arc + k = 1.0 + np.dot(direction, dir_new) + s = 1/np.sqrt(k+k) + r = s*np.cross(direction, dir_new) + theta = k*s + #r = np.cross(direction, dir_new) + #r = r/np.sqrt(r[0]*r[0] + r[1]*r[1] + r[2]*r[2]) + #calculate the degree of quaternion rotation. + #cos_theta = np.sqrt(np.sqrt(np.dot(dir_new, dir_new)) * np.sqrt(np.dot(direction, direction))) + np.dot(dir_new, direction) + #cos_theta = np.dot(direction, dir_new) + #theta = np.arccos(cos_theta)/2.0 + #print(r, cos_theta, theta) + #quat = np.append(theta, r) + q = Quaternion(w=theta, x = r[0], y = r[1], z = r[2]).normalize() + + #print(quat) + #q = np.quaternion(quat[0], quat[1], quat[2], quat[3]) + #rot = Rotation.from_quat(quat, normalized=False) + #rot.as_quat() + for s in range(num_sides): + circle_pts[p][s] = q.rotate_point(circle_pts[p-1][s]) + #circle_pts[p][s] = rot.apply(circle_pts[p-1][s]) + #circle_pts[p][s] = q.rotate(circle_pts[p-1][s]) + #circle_pts[p][s] = np.quaternion.rotate_vectors(q, circle_pts[p][s]) + #circle_pts[p][s] = q.rotate_vectors(q, circle_pts[p][s]) + #circle_pts[p][s] = circle_pts[p][s] + direction = dir_new + #generate the triangles + triangles = np.random.randint(size=((pts.shape[0]-1)*2*(num_sides), 3), low=0, + high=(G.num_edges()-1)).astype(np.uint32) + + t_index_temp = 0 + for ring in range(pts.shape[0]-1): + for side in range(num_sides): + if(side < num_sides-1): + idx_current_point = index+ring*num_sides + side + idx_next_ring = index + (ring+1) * num_sides + side + triangle1 = [idx_current_point, idx_next_ring, idx_next_ring+1] + triangle2 = [idx_next_ring+1, idx_current_point+1, idx_current_point] + triangles[t_index_temp] = [idx_current_point, idx_next_ring, idx_next_ring+1] + triangles[t_index_temp+1] = [idx_next_ring+1, idx_current_point+1, idx_current_point] + self.triangle_data[t_index] = triangle1 + self.triangle_data[t_index+1] = triangle2 + t_index += 2 + t_index_temp += 2 + else: + idx_current_point = index + ring*num_sides + side + idx_next_ring = index + (ring+1)*num_sides + side + triangle1 = [idx_current_point, idx_next_ring, idx_next_ring-num_sides+1] + triangle2 = [idx_next_ring-num_sides+1, idx_current_point-num_sides+1, idx_current_point] + triangles[t_index_temp] = [idx_current_point, idx_next_ring-num_sides, idx_next_ring-num_sides+1] + triangles[t_index_temp+1] = [idx_next_ring-num_sides+1, idx_current_point-num_sides+1, idx_current_point] + self.triangle_data[t_index] = triangle1 + self.triangle_data[t_index+1] = triangle2 + t_index += 2 + t_index_temp += 2 + + #generate the points data structure + circle_pts_data = circle_pts.reshape((pts.shape[0]*num_sides, 3)) + self.cylinder_data['a_position'][index:(pts.shape[0]*num_sides+index)] = circle_pts_data + self.cylinder_data['a_fg_color'][index:(pts.shape[0]*num_sides+index)] = color + #generate the normals data structure + pts_normals = circle_pts.copy() + for p in range(pts.shape[0]): + pts_normals[p][:] = pts_normals[p][:] - pts[p] + for s in range(num_sides): + pts_normals[p][s] = pts_normals[p][s]/np.sqrt(pts_normals[p][s][0]*pts_normals[p][s][0] \ + + pts_normals[p][s][1]*pts_normals[p][s][1] + pts_normals[p][s][2]*pts_normals[p][s][2]) + self.cylinder_data['a_normal'][index:(pts.shape[0]*num_sides+index)] = \ + pts_normals.reshape((pts.shape[0]*num_sides, 3)) + + index += pts.shape[0]*num_sides + + #Add the caps for each of the endpoints. + + + + if(i == 2): + fig = plt.figure() + ax = fig.add_subplot(111, projection='3d') + #ax.scatter(circle_pts[:,:,0], circle_pts[:,:,1], circle_pts[:,:,2]) + ax.plot(pts[:,0], pts[:,1], pts[:,2]) + for j in range(pts.shape[0]): + ax.plot(circle_pts[j,:,0], circle_pts[j,:,1], circle_pts[j,:,2]) + for j in range(triangles.shape[0]): + tri = np.zeros((3,4)) + tri[:,0] = self.cylinder_data['a_position'][triangles[j][0]] + tri[:,1] = self.cylinder_data['a_position'][triangles[j][1]] + tri[:,2] = self.cylinder_data['a_position'][triangles[j][2]] + tri[:,3] = self.cylinder_data['a_position'][triangles[j][0]] + ax.plot(tri[0,:], tri[1,:], tri[2,:], c='b') + for j in range(triangles.shape[0]): + tri = np.zeros((3,3)) + tri[:,0] = self.cylinder_data['a_position'][triangles[j][0]] + tri[:,1] = self.cylinder_data['a_position'][triangles[j][1]] + tri[:,2] = self.cylinder_data['a_position'][triangles[j][2]] + norm = np.zeros((3,3)) + norm[:,0] = self.cylinder_data['a_normal'][triangles[j][0]] + norm[:,1] = self.cylinder_data['a_normal'][triangles[j][1]] + norm[:,2] = self.cylinder_data['a_normal'][triangles[j][2]] + ax.quiver(tri[0,:], tri[1,:], tri[2,:], norm[0,:], norm[1,:], norm[2,:], colors = 'r') + plt.show() + i+=1 + #create the data. + + #Handles the mouse wheel event, i.e., zoom + def on_mouse_wheel(self, event): +# self.scale[0] = self.scale[0] + self.scale[0]*event.delta[1]*0.05 +# self.scale[1] = self.scale[1] + self.scale[1]*event.delta[1]*0.05 +# self.scale[2] = self.scale[2] + self.scale[2]*event.delta[1]*0.05 +## self.view[0][0] = self.scale[0] +## self.view[1][1] = self.scale[1] +## self.view[2][2] = self.scale[2] +# print("in mouse wheel ", self.r1, self.r2) +# self.view = np.matmul(translate((self.translate[0], self.translate[1], self.translate[2])), self.r1) +# self.view = np.matmul(self.view, self.r2) +# self.view = np.matmul(self.view, scale((self.scale[0], self.scale[1], self.scale[2]))) +# #self.view = np.matmul(self.view, self.r1) +# #self.view = np.matmul(self.view, self.r2) +# #self.view = np.matmul(translate((self.translate[0], self.translate[1], self.translate[2])), scale((self.scale[0], self.scale[1], self.scale[2]))) +# #self.rotate = [0., 0.] +# #self.camera = +# +# #self.view[1][1] = self.view[1][1]+self.view[1][1]*event.delta[1]*0.05 +# #print(self.view[0][0], " ",self.view[1][1]) +# #print(self.view) +# self.camera = [0.0, 0.0, -100.0, 1.0] +## for i in range(len(self.translate)): +## self.camera[i] += self.translate[i] +# self.program['u_view'] = self.view + +# if(np.asarray(self.camera).all() == np.asarray([0., 0., 0.]).all()): +# self.camera = np.asarray([0., 0., 0.]) +# else: + direction = np.asarray(self.translate) - np.asarray(self.camera) + direction = direction/np.sqrt(np.dot(direction, direction)) + for i in range(3): + self.camera[i] = self.camera[i] + direction[i]*event.delta[1]*2.0 + + self.program['u_eye'] = self.camera + #print(self.view) + #print(event.delta[1]) + self.update() + + + + #Handles the mouse press event to rotate the camera if the left mouse button + #if clicked down. + def on_mouse_press(self, event): + def update_view(): + self.location = event.pos + #self.program['u_view'] = self.view + self.down = True + + if(event.button == 1): + update_view() + + + #Handles the rotation of the camera using a quaternion centered around the + #focus point. + def on_mouse_move(self, event): + if(self.down == True): + coord = self.transforms.get_transform('canvas', 'render').map(self.location)[:2] + coord2 = self.transforms.get_transform('canvas', 'render').map(event.pos)[:2] + #coord[1] = 0 + #coord2[1] = 0 + + theta = (coord[0]-coord2[0])*360.0/2.0/np.pi + phi = (coord[1]-coord2[1])*360.0/2.0/np.pi + print(theta*360.0/2.0/np.pi, -phi*360.0/2.0/np.pi) + self.camera = self.camera - self.translate + q1 = Quaternion.create_from_axis_angle(angle=phi, ax=0.0, ay=1.0, az=0.0, degrees=True) + q2 = Quaternion.create_from_axis_angle(angle=theta, ax=1.0, ay=0.0, az=0.0, degrees=True) + #q1 = Quaternion(w=theta, x = 0, y = 1, z = 0).inverse().normalize() + #q2 = Quaternion(w=-phi, x = 1, y = 0, z = 0).inverse().normalize() + + q = q1*q2 + +# #print("Angle in Degrees is ", theta, " ", phi, coord[0] - coord2[0]) +# self.r1 = rotate(theta, (0, 1, 0)) +# self.r2 = rotate(-phi, (1, 0, 0)) +# +# print("in on mouse move ", self.r1, self.r2) +# +# self.view = np.matmul(self.view, self.r1) +# #print("r1", self.view) +# self.view = np.matmul(self.view, self.r2) +# #print("r2", self.view) +# +## self.view = np.matmul(self.view, q1.get_matrix().T) +## self.view = np.matmul(self.view, q2.get_matrix().T) +# +# self.program['u_view'] = self.view +# + self.location = event.pos +# #print("Angle in Degrees is ", theta, " ", phi) +# #print(self.camera) + self.camera = np.asarray(q.rotate_point(self.camera), dtype=np.float) + self.camera = self.camera + self.translate + self.up = np.asarray(q.rotate_point(self.up), dtype=np.float) + self.up = self.up/np.sqrt(np.dot(self.up, self.up)) + #self.rotate[0] = self.rotate[0] + theta + #self.rotate[1] = self.rotate[1] + phi + #print(self.rotate)f + #radius = np.sqrt(np.dot(self.center, self.center))*2 + #test = np.linalg.inv(self.view).T + #print(test) + + + + #self.camera = sph2cart(self.rotate[0]/360.0*2.0*np.pi+np.pi, self.rotate[1]/360.0*2.0*np.pi+np.pi, 1000.0) + #self.camera[0] = self.camera[0] + self.center[0] + #self.camera[1] = self.camera[1] + self.center[1] + #self.camera[2] = self.camera[2] - self.center[2] + print("current position ", self.camera, " and up vector ", self.up) + self.program['u_eye'] = self.camera + self.program['u_up'] = self.up + self.program['u_LightPos'] = [self.camera[0], self.camera[1], self.camera[2]] + self.update() + + #reverts the mouse state during release. + def on_mouse_release(self, event): + self.down = False diff --git a/TubeWidget.py b/TubeWidget.py new file mode 100644 index 0000000..2747819 --- /dev/null +++ b/TubeWidget.py @@ -0,0 +1,65 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Created on Mon Aug 5 15:53:16 2019 + +@author: pavel +""" + +""" + Qt wrapper/container class that handles all the top level events and maintains the + methods necessary to use the QT signals and slots API. +""" + +from pyqtgraph.Qt import QtCore, QtGui, QtWidgets +from TubeCanvas import TubeDraw + +class TubeWidget(QtGui.QWidget): + sigUpdate = QtCore.pyqtSignal(float, float, float) + #Initializes the QT wrapper class. + def __init__(self): + super(TubeWidget, self).__init__() + box = QtGui.QVBoxLayout(self) + self.resize(500,500) + self.setLayout(box) + self.canvas = TubeDraw() + #self.canvas.create_native() + box.addWidget(self.canvas.native) + self.camera = [0,0,0] + self.down = False + + self.canvas.events.mouse_press.connect(self.on_mouse_press) + self.canvas.events.mouse_release.connect(self.on_mouse_release) + self.canvas.events.mouse_move.connect(self.on_mouse_move) + self.canvas.events.mouse_wheel.connect(self.on_mouse_wheel) + + #self.show() + + #Handles the mouse release event + def on_mouse_release(self, event): + self.down=False + self.sendCameraInfo() + + #Handles the mouse move event + def on_mouse_move(self, event): + if self.down: + self.sendCameraInfo() + + #Handles the mouse press event + def on_mouse_press(self, event): + self.down = True + n = 3 + + #Handles the mouse wheel event + def on_mouse_wheel(self, event): + self.sendCameraInfo() + + #controls the emit function of the QT class to send a signal to the slot + #located in the other window. + def sendCameraInfo(self): + #self.view = self.canvas.view + self.camera = self.canvas.camera + #print("stuff", self.view[3, 0], self.view[3, 1], self.view[3, 2]) + print("stuff", self.camera[0], self.camera[1], self.camera[2]) + self.sigUpdate.emit(self.camera[0], self.camera[1], self.camera[2]) + diff --git a/graph_shaders.py b/graph_shaders.py new file mode 100644 index 0000000..cc799b2 --- /dev/null +++ b/graph_shaders.py @@ -0,0 +1,256 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Created on Mon Aug 5 14:48:00 2019 +Collection of shaders necessary to draw all elements of the graph. +@author: pavel +""" + +from vispy import gloo, app, scene + +#shaders for drawing nodes + +vert = """ +#version 120 +// Uniforms +// ------------------------------------ +uniform mat4 u_model; +uniform mat4 u_view; +uniform mat4 u_projection; +uniform float u_antialias; +uniform float u_size; +uniform bool u_picking; + +uniform vec3 u_graph_size; +// Attributes +// ------------------------------------ +attribute vec3 a_position; +attribute vec4 a_fg_color; +attribute vec4 a_bg_color; +attribute float a_linewidth; +attribute float a_size; +attribute vec4 a_unique_id; +attribute float a_selection; +// Varyings +// ------------------------------------ +varying vec4 v_fg_color; +varying vec4 v_bg_color; +varying float v_size; +varying float v_linewidth; +varying float v_antialias; +varying vec4 v_unique_id; +varying float v_selection; + +varying float v_zoom_level; +void main (void) { + v_size = a_size * u_size; + v_linewidth = a_linewidth; + v_antialias = u_antialias; + v_fg_color = a_fg_color; + v_bg_color = a_bg_color; + gl_Position = u_projection * u_view * u_model * + vec4(a_position*u_size,1.0); + gl_PointSize = v_size + 2*(v_linewidth + 1.5*v_antialias); + v_zoom_level = u_view[0][0]; + v_unique_id = a_unique_id; + v_selection = a_selection; +} +""" + +frag = """ +#version 120 +// Constants +// ------------------------------------ +uniform mat4 u_model; +uniform mat4 u_view; +uniform mat4 u_projection; +uniform float u_antialias; +uniform float u_size; +uniform bool u_picking; + +// Varyings +// ------------------------------------ +varying vec4 v_fg_color; +varying vec4 v_bg_color; +varying float v_size; +varying float v_linewidth; +varying float v_antialias; +varying vec4 v_unique_id; + +varying float v_zoom_level; +varying float v_selection; +// Functions +// ------------------------------------ +float marker(vec2 P, float size); +float new_alpha(float zoom_level); +// Main +// ------------------------------------ +void main() +{ + if(!u_picking) + { + float size = v_size + 2*(v_linewidth + 1.5*v_antialias); + float t = v_linewidth/2.0-v_antialias; + // The marker function needs to be linked with this shader + float r = marker(gl_PointCoord, size); + float d = abs(r) - t; + if( r > (v_linewidth/2.0+v_antialias)) + { + discard; + } + else if( d < 0.0 ) + { + if(v_zoom_level < 0.0010) + { + float alpha = d/v_antialias; + alpha = new_alpha(v_zoom_level); + gl_FragColor = mix(vec4(1,1,1,0.9), v_bg_color, max(1.0-alpha, 0.3)); + } + else + { + if(v_selection == 1.0) + gl_FragColor = vec4(1.0, 0.0, 0.0, v_fg_color.a); + else if(v_selection == 2.0) + gl_FragColor = vec4(0.0, 1.0, 0.0, v_fg_color.a); + else + gl_FragColor = v_fg_color; + } + } + else + { + float alpha = d/v_antialias; + if(v_zoom_level < 0.0010) + alpha = new_alpha(v_zoom_level); + else + alpha = exp(-alpha*alpha); + + if (r > 0 && v_zoom_level >= 0.0010) + gl_FragColor = vec4(v_fg_color.rgb, alpha*v_fg_color.a); + else + if(v_zoom_level < 0.0010) + if(vec3(gl_Color.rgb) != vec3(v_bg_color.rgb)) + gl_FragColor = mix(vec4(1,1,1,0.9), v_bg_color, max(1.0-alpha, 0.3)); + else + gl_FragColor = vec4(gl_Color.rgb, max(1.0-alpha, 0.1)); + else + gl_FragColor = mix(v_bg_color, v_fg_color, alpha); + } + + + } else { + float size = v_size +2*(v_linewidth + 1.5*v_antialias); + float t = v_linewidth/2.0-v_antialias; + // The marker function needs to be linked with this shader + float r = marker(gl_PointCoord, size); + float d = abs(r) - t; + float alpha = d/v_antialias; + if(v_zoom_level < 0.0010) + alpha = new_alpha(v_zoom_level); + else + alpha = exp(-alpha*alpha); + if(r > 0) + gl_FragColor = v_unique_id; + else + gl_FragColor = v_unique_id; + + } +} + +float marker(vec2 P, float size) +{ + float r = length((P.xy - vec2(0.5,0.5))*size); + r -= v_size/2; + return r; +} + +float new_alpha(float zoom_level) +{ + float val = (zoom_level - 0.0010)/(0.00075-0.0010); + if(val < 0) + { + val = 0; + } + return val; +} +""" + +#Shaders for the lines between the nodes + + +vs = """ +#version 120 + +// Uniforms +// ------------------------------------ +uniform mat4 u_model; +uniform mat4 u_view; +uniform mat4 u_projection; + +// Attributes +// ------------------------------------ +attribute vec3 a_position; +attribute vec2 a_normal; +attribute vec4 a_fg_color; +attribute float a_linewidth; +//attribute vec4 a_unique_id; +//attribute vec4 l_color; + +// Varyings +// ------------------------------------ +varying vec4 v_fg_color; +varying float v_zoom_level; +varying vec2 v_normal; +varying float v_linewidth; + +void main() { + vec3 delta = vec3(a_normal*a_linewidth/(1-u_view[0][0]), 0); + //vec4 pos = u_model * u_view * vec4(a_position, 1.0); + //gl_Position = u_projection * (pos + vec4(delta, 1.0)); + gl_Position = u_model * u_view * u_projection * vec4(a_position+delta, 1.0); + //gl_Position = u_projection * u_view * u_model * + // vec4(a_position, 1.0); + v_zoom_level = u_view[0][0]; + v_fg_color = a_fg_color; + v_normal = a_normal; + v_linewidth = a_linewidth; + +} +""" + +fs = """ +#version 120 +// Varying +// ------------------------------------ +varying vec4 v_fg_color; +varying float v_zoom_level; +varying vec2 v_normal; +varying float v_linewidth; + +float new_alpha(float zoom_level); + +void main() +{ + float l = length(v_normal); + float feather = 0.5; + float alpha = 1.0; + if(l > v_linewidth/2.0+feather) + discard; + else + alpha = 0.5; + + if(v_zoom_level < 0.0010) + alpha = new_alpha(v_zoom_level); + gl_FragColor = vec4(v_fg_color.rgb, alpha); +} + +float new_alpha(float zoom_level) +{ + float val = (zoom_level-0.00075)/(0.0010-0.00075); + if(val < 0.) + { + val = 0.; + } + return val; + +} +""" \ No newline at end of file diff --git a/network_dep.py b/network_dep.py new file mode 100755 index 0000000..4e43b88 --- /dev/null +++ b/network_dep.py @@ -0,0 +1,2175 @@ +# -*- coding: utf-8 -*- +""" +Created on Sat Sep 16 16:34:49 2017 + +@author: pavel +""" + +import struct +import numpy as np +from scipy.sparse.linalg import eigsh +import scipy as sp +from sklearn.cluster import SpectralClustering + + +#import networkx as nx +import matplotlib.pyplot as plt +from matplotlib import cm +import math +import time +#import spharmonics +import graph_tool.all as gt +import copy +#import matplotlib + +#for testing +#import matplotlib.mlab as mlab +#import matplotlib.pyplot as plt +#from matplotlib import cm + +''' + Definition of the Node class + Duplicate of the node class in network + Stores the physical position, outgoing edges list and incoming edges list. +''' +class Node: + def __init__(self, point, outgoing, incoming): + self.p = point + self.o = outgoing + self.i = incoming + +# def p(): +# return self.p + +''' + Definition of the Fiber class. + Duplicate of the Node class in network + Stores the starting vertex, the ending vertex, the points array and the radius array +''' +class Fiber: + + + def __init__ (self, p1, p2, pois, rads): + self.v0 = p1 + self.v1 = p2 + self.points = pois + self.radii = rads + ''' + return the array-likes of the x,y,z,r coordinates of the fiber for the + gt representation + ''' + + def getcoords(self): + x = np.zeros(len(self.points), dtype=np.double) + y = np.zeros(len(self.points), dtype=np.double) + z = np.zeros(len(self.points), dtype=np.double) + r = np.zeros(len(self.points), dtype=np.double) + for i in range(len(self.points)): + x[i] = self.points[i][0] + y[i] = self.points[i][1] + z[i] = self.points[i][2] + r[i] = self.radii[i] + + return x,y,z,r + + + ''' + return the length of the fiber. + ''' + def length(self): + length = 0 + for i in range(len(self.points)-1): + length = length + math.sqrt(pow(self.points[i][0]- self.points[i+1][0],2) + pow(self.points[i][1]- self.points[i+1][1],2) + pow(self.points[i][2]- self.points[i+1][2],2)) + if(length == 0): + print(self.points) + print(len(self.points)) + print(self.v0, " ", self.v1) + return length + + ''' + returns the turtuosity of the fiber. + ''' + def turtuosity(self): + turtuosity = 1.0 + distance = math.sqrt(math.pow(self.points[0][0]- self.points[len(self.points)-1][0],2) + math.pow(self.points[0][1]- self.points[len(self.points)-1][1],2) + math.pow(self.points[0][2]- self.points[len(self.points)-1][2],2)) + if(distance > 0): + turtuosity = self.length()/distance + #print(turtuosity) + + return turtuosity + + ''' + returns the volume of the fiber. + ''' + def volume(self): + volume = 0 + for i in range(len(self.points)-1): + volume = volume + 1.0/3.0 * math.pi * (math.pow(self.radii[i]*10e-6,2) + math.pow(self.radii[i+1]*10e-6,2) + self.radii[i]*10e-6*self.radii[i+1]*10e-6) * math.sqrt(math.pow(self.points[i][0]- self.points[i+1][0],2) + math.pow(self.points[i][1]- self.points[i+1][1],2) + math.pow(self.points[i][2]- self.points[i+1][2],2)) + + #print(volume) + return volume + + def av_radius(self): + radius = 0.0 + for i in range(len(self.radii)): + radius = radius + self.radii[i] + + return radius/len(self.radii) + +class NWT: + + ''' + Writes the header given and open file descripion, number of verticies and number of edges. + ''' + def writeHeader(open_file, numVerts, numEdges): + txt = "nwtFileFormat fileid(14B), desc(58B), #vertices(4B), #edges(4B): bindata" + b = bytearray() + b.extend(txt.encode()) + open_file.write(b) + open_file.write(struct.pack('i', numVerts)) + open_file.write(struct.pack('i', numEdges)) + + + ''' + Writes a single vertex to a file. + ''' + def writeVertex(open_file, vertex): + open_file.write(struct.pack(' self.B[0]: + self.B[0] = G.vertex_properties["p"][v][0] + if G.vertex_properties["p"][v][1] < self.A[1]: + self.A[1] = G.vertex_properties["p"][v][1] + if G.vertex_properties["p"][v][1] > self.B[1]: + self.B[1] = G.vertex_properties["p"][v][1] + if G.vertex_properties["p"][v][2] < self.A[2]: + self.A[2] = G.vertex_properties["p"][v][2] + if G.vertex_properties["p"][v][2] > self.B[2]: + self.B[2] = G.vertex_properties["p"][v][2] + #In case of a dual graph we have to scane the first and last points + # of every fiber to find the true bounding box + else: + for v in G.vertices(): + l = len(G.vertex_properties["x"][v]) + if G.vertex_properties["x"][v][0] < self.A[0]: + self.A[0] = G.vertex_properties["x"][v][0] + if G.vertex_properties["x"][v][l-1] < self.A[0]: + self.A[0] = G.vertex_properties["x"][v][l-1] + if G.vertex_properties["y"][v][0] < self.A[1]: + self.A[1] = G.vertex_properties["y"][v][0] + if G.vertex_properties["y"][v][l-1] < self.A[1]: + self.A[1] = G.vertex_properties["y"][v][l-1] + if G.vertex_properties["z"][v][0] < self.A[2]: + self.A[2] = G.vertex_properties["z"][v][0] + if G.vertex_properties["z"][v][l-1] < self.A[2]: + self.A[2] = G.vertex_properties["z"][v][l-1] + + if G.vertex_properties["x"][v][0] > self.B[0]: + self.B[0] = G.vertex_properties["x"][v][0] + if G.vertex_properties["x"][v][l-1] > self.B[0]: + self.B[0] = G.vertex_properties["x"][v][l-1] + if G.vertex_properties["y"][v][0] > self.B[1]: + self.B[1] = G.vertex_properties["y"][v][0] + if G.vertex_properties["y"][v][l-1] > self.B[1]: + self.B[1] = G.vertex_properties["y"][v][l-1] + if G.vertex_properties["z"][v][0] > self.B[2]: + self.B[2] = G.vertex_properties["z"][v][0] + if G.vertex_properties["z"][v][l-1] > self.B[2]: + self.B[2] = G.vertex_properties["z"][v][l-1] + #print(self.A, self.B) + self.O = (self.A+self.B)*0.5 + self.planes = [] + #append x planes described by a point and a vector + self.planes.append((np.array([self.A[0], 0.0, 0.0]), np.array([1.0, 0.0, 0.0]))) + self.planes.append((np.array([self.B[0], 0.0, 0.0]), np.array([-1.0, 0.0, 0.0]))) + #append y planes described by a point and a vector + self.planes.append((np.array([0.0, self.A[1], 0.0]), np.array([0.0, 1.0, 0.0]))) + self.planes.append((np.array([0.0, self.B[1], 0.0]), np.array([0.0, -1.0, 0.0]))) + #append z planes described by a point and a vector + self.planes.append((np.array([0.0, 0.0, self.A[2]]), np.array([0.0, 0.0, 1.0]))) + self.planes.append((np.array([0.0, 0.0, self.B[2]]), np.array([0.0, 0.0, -1.0]))) + + + def distance(self, pt): + dist = 10000000 + for i in self.planes: + V = i[0]-pt + if np.dot(V, i[1]) < dist: + dist = np.dot(V, i[1]) + + return dist + + def project_grid(self, n): + #r = abs(self.A - self.B) + x = np.linspace(self.A[0], self.B[0], (n+2)) + x = x[1:-1] + y = np.linspace(self.A[1], self.B[1], (n+2)) + y = y[1:-1] + z = np.linspace(self.A[2], self.B[2], (n+2)) + z = z[1:-1] + + return x,y,z + + #returns a resampled bounding both with nxn points on each side. + def resample_sides(self, n): + #get the size of the cube in every cardinal direction + size = abs(self.A - self.B) + #set the stepsize for each subdivided point + dist_x = size[0]/(n-1) + dist_y = size[1]/(n-1) + dist_z = size[2]/(n-1) + #generate the original 8 corners + vertices = [] + #generate the points for the yz planes. + for i in range(n): + for j in range(n): + #generate the temporary vectors for both the planes + temp_p1 = copy.deepcopy(self.A) + temp_p2 = copy.deepcopy(self.A) + temp_p2[0] = self.B[0] + + temp_p1[1] = temp_p1[1] + dist_y*i + temp_p1[2] = temp_p1[2] + dist_z*j + + temp_p2[1] = temp_p2[1] + dist_y*i + temp_p2[2] = temp_p2[2] + dist_z*j + + vertices.append(temp_p1) + vertices.append(temp_p2) + #generate the points for the xz planes. + for i in range(n): + for j in range(n): + #generate the temporary vectors for both the planes + temp_p1 = copy.deepcopy(self.A) + temp_p2 = copy.deepcopy(self.A) + temp_p2[1] = self.B[1] + + temp_p1[0] = temp_p1[0] + dist_x*i + temp_p1[2] = temp_p1[2] + dist_z*j + + temp_p2[0] = temp_p2[0] + dist_x*i + temp_p2[2] = temp_p2[2] + dist_z*j + + vertices.append(temp_p1) + vertices.append(temp_p2) + + #generate the points for the xy planes. + for i in range(n): + for j in range(n): + #generate the temporary vectors for both the planes + temp_p1 = copy.deepcopy(self.A) + temp_p2 = copy.deepcopy(self.A) + temp_p2[2] = self.B[2] + + temp_p1[0] = temp_p1[0] + dist_x*i + temp_p1[1] = temp_p1[1] + dist_y*j + + temp_p2[0] = temp_p2[0] + dist_x*i + temp_p2[1] = temp_p2[1] + dist_y*j + + vertices.append(temp_p1) + vertices.append(temp_p2) + + + vertices = list(np.unique(np.array(vertices), axis=0)) + import matplotlib.pyplot as plt + fig = plt.figure() + ax = plt.axes(projection='3d') + ax.scatter3D(np.array(vertices)[:, 0], np.array(vertices)[:, 1], np.array(vertices)[:, 2]) + + print("THERE ARE THIS MANY ", len(vertices)) + return vertices + + def getVolume(self): + size = abs(self.A - self.B) + return size[0]*size[1]*size[2] + + def vertices(self): + size = abs(self.A - self.B) + points = [] + points.append(self.A) + + + temp = copy.deepcopy(self.A) + temp[0] = temp[0] + size[0] + points.append(temp) + print('1', temp) + + temp = copy.deepcopy(self.A) + temp[1] = temp[1] + size[1] + points.append(temp) + print('1', temp) + + temp = copy.deepcopy(self.A) + temp[2] = temp[2] + size[2] + points.append(temp) + print('1', temp) + + temp = copy.deepcopy(self.A) + temp[1] = temp[1] + size[1] + temp[0] = temp[0] + size[0] + points.append(temp) + print('1', temp) + + temp = copy.deepcopy(self.A) + temp[2] = temp[2] + size[2] + temp[0] = temp[0] + size[0] + points.append(temp) + print('1', temp) + + temp = copy.deepcopy(self.A) + temp[1] = temp[1] + size[1] + temp[2] = temp[2] + size[2] + points.append(temp) + print('1', temp) + + temp = copy.deepcopy(self.A) + temp[0] = temp[0] + size[0] + temp[1] = temp[1] + size[1] + temp[2] = temp[2] + size[2] + points.append(temp) + print('1', temp) + + return points + + +''' +Class defining the BFS traversal of all the vertices in the graph and labeling +all interconnected vertices. +''' +class VisitorClassDisconnected(gt.BFSVisitor): + def __init__(self, ecluster, cluster, c): + self.ecluster = ecluster + self.cluster = cluster + self.c = c + + def examine_vertex(self, u): + if(self.cluster[u] == -1): + self.cluster[u] = self.c + + def examine_edge(self, e): + if(self.ecluster[e] == -1): + self.ecluster[e] = self.c + + +''' +Class defining the BFS traversal of all the vertices in the graph and labeling +all interconnected vertices. +''' +class VisitorClassPartition(gt.BFSVisitor): + def __init__(self, cluster, dist, c, iterations): + self.cluster = cluster + self.dist = dist + self.c = c + self.total_iter = iterations + + def examine_vertex(self, u): + if(self.cluster[u] == -1): + self.cluster[u] = self.c + + def tree_edge(self, e): + d = self.dist[e.source()] + #print(d, self.total_iter) + if d <= self.total_iter: + self.dist[e.target()] = self.dist[e.source()] + 1 + +class Network: + def __init__(self, filename, clock=False): + if clock: + start_time = time.time() + + with open(filename, "rb") as file: + header = file.read(72) + bytes = file.read(4) + numVertex = int.from_bytes(bytes, byteorder='little') + bytes = file.read(4) + numEdges = int.from_bytes(bytes, byteorder='little') + + self.N = [] + self.F = [] + for i in range(numVertex): + node = NWT.readVertex(file) + self.N.append(node) + + for i in range(numEdges): + edge = NWT.readFiber(file) + self.F.append(edge) + if clock: + print("Network initialization: " + str(time.time() - start_time) + "s") + ''' + Takes in a graph object (networkx Graph) and saves it as a network file + Resulting nwt, does not contain incomind or outgoing fibers + ''' + def saveGraph(self, G, path): + G_nodes = list(G.nodes()) + G_edges = list(G.edges()) + Nodes = [] + Edges = [] + points = list(nx.get_node_attributes(G, 'p').values()) + for i in range(len(G_nodes)): + node = Node(points[i], [], []) + Nodes.append(node) + for e in G_edges: + edge = Fiber(e[0], e[1], G[e[0]][e[1]]['pts'], G[e[0]][e[1]]['rads']) + Edges.append(edge) + + NWT.exportNWT(path, Nodes, Edges) + + ''' + Saves the graph-tool graph as a nwt file with all the variables. + ''' + + def saveGraph_gt(self, G, path): + points = G.vertex_properties["p"] + Nodes = [] + Edges = [] + for i in range(G.num_vertices()): + #array = G.get_out_edges(i) + outgoing = [] + for edge in G.get_out_edges(i): + outgoing.append(edge[2]) + + node = Node(points[i], outgoing, outgoing) + Nodes.append(node) + + for e in G.edges(): + x = G.edge_properties["x"][e].get_array() + y = G.edge_properties["y"][e].get_array() + z = G.edge_properties["z"][e].get_array() + r = G.edge_properties["r"][e].get_array() + pts = [] + for i in range(len(x)): + pt = [] + pt.append(x[i]) + pt.append(y[i]) + pt.append(z[i]) + pts.append(pt) + #pts = np.ndarray.tolist(np.dstack([x,y,z])) + #print(pts) + edge = Fiber(int(e.source()), int(e.target()), pts, np.ndarray.tolist(r)) + if(edge.length() > 0.0): + Edges.append(edge) + + NWT.exportNWT(path, Nodes, Edges) + + + + def saveGraphStatistics_gt(self, G, path, prefix, label = "none", norm=False): + location = path + "/" + prefix + "_edges.txt" + f = open(location, "a+") + f.write("inverted\t") + f.write("length\t") + f.write("tortuosity\t") + f.write("volume\t") + f.write("inverted_volume\t") + f.write("gaussian\t") + f.write("bc\t") + f.write("bc*length\t") + f.write("mst\t\n") + for e in G.edges(): + f.write("%.15f\t" % G.edge_properties["inverted"][e]) + f.write("%.15f\t" % G.edge_properties["length"][e]) + f.write("%.15f\t" % G.edge_properties["tortuosity"][e]) + f.write("%.15f\t" % G.edge_properties["volume"][e]) + f.write("%.15f\t" % G.edge_properties["inverted_volume"][e]) + f.write("%.15f\t" % G.edge_properties["gaussian"][e]) + f.write("%.15f\t" % G.edge_properties["bc"][e]) + f.write("%.15f\t" % G.edge_properties["bc*length"][e]) + f.write("%.15f\t" % G.edge_properties["mst"][e]) + f.write("%s\n" % label) + + f.close() + + location = path+"/" + prefix + "_vertices.txt" + f = open(location, "a+") + f.write("bc\t") + f.write("degree\t") + f.write("degree_volume\t") + f.write("degree_tortuosity\t\n") + for v in G.vertices(): + f.write("%.15f\t" % G.vertex_properties["bc"][v]) + f.write("%.15f\t" % G.vertex_properties["degree"][v]) + f.write("%.15f\t" % G.vertex_properties["degree_volume"][v]) + f.write("%.15f\t" % G.vertex_properties["degree_tortuosity"][v]) + f.write("%s\n" % label) + + ''' + Saves the graph as a sample. Additionally saves the key as a list of string if writeKey is True + saves the file at path with prefix_s.txt + ''' + def saveSample_gt(self, G, path, prefix, is_dual, label = "none", norm=False, writeKey=False): + location = path + "/" + prefix + "s.txt" + f = open(location, "a+") + + + aabb = AABB(G, is_dual) + array = G.vertex_properties["bc"].get_array() + av_array = np.sum(array)/G.num_vertices() + G.graph_properties["global_vertex_bc"] = G.new_graph_property("double", val = av_array) + f.write("%.15f\t" % av_array) + + array = G.vertex_properties["bc"].get_array() + av_array = np.sum(array)/aabb.getVolume() + G.graph_properties["global_vascular_vertex_bc"] = G.new_graph_property("double", val = av_array) + f.write("%.15f\t" % av_array) + + array = G.edge_properties["bc"].get_array() + av_array = np.sum(array)/aabb.getVolume() + G.graph_properties["global_vascular_edge_bc"] = G.new_graph_property("double", val = av_array) + f.write("%.15f\t" % av_array) + + #G.edge_properties["mst"] = gt.graph_tool.topology.min_spanning_tree(G) + G.graph_properties["global_mst_ratio"] = G.new_graph_property("double", val = np.double(np.sum(G.edge_properties["mst"].get_array()))/np.double(G.num_edges())) + f.write("%.15f\t" % G.graph_properties["global_mst_ratio"]) + + #G.edge_properties["mst"] = gt.graph_tool.topology.min_spanning_tree(G, weights=G.edge_properties["volume"]) + G.graph_properties["global_mst_ratio_vascular_volume"] = G.new_graph_property("double", val = np.double(np.sum(G.edge_properties["mst"].get_array()))/np.double(G.num_edges())/aabb.getVolume()) + f.write("%.15f\t" % G.graph_properties["global_mst_ratio"]) + + ########################Local Statistics############################### + self.recalculate_metrics(G, is_dual, norm) + ########################Biological Statistics########################## + # Some biological statistics are reported as a ratio to physical volume of vessels + array = G.edge_properties["length"].get_array() + av_array = np.sum(array)/G.num_edges() + G.graph_properties["average_length"] = G.new_graph_property("double", val= av_array) + f.write("%.15f\t" % av_array) + + array = G.edge_properties["volume"].get_array() + av_array = np.sum(array)/G.num_edges() + G.graph_properties["average_volume"] = G.new_graph_property("double", val= av_array) + f.write("%.15f\t" % av_array) + + array = G.edge_properties["tortuosity"].get_array() + av_array = np.sum(array)/G.num_edges() + G.graph_properties["average_tortuosity"] = G.new_graph_property("double", val= av_array) + f.write("%.15f\t" % av_array) + + array = G.edge_properties["length"].get_array() + av_array = np.sum(array)/aabb.getVolume() + G.graph_properties["vascular length"] = G.new_graph_property("double", val= av_array) + f.write("%.15f\t" % av_array) + + array = G.edge_properties["volume"].get_array() + av_array = np.sum(array)/aabb.getVolume() + G.graph_properties["vascular_volume"] = G.new_graph_property("double", val= av_array) + f.write("%.15f\t" % av_array) + + array = G.edge_properties["tortuosity"].get_array() + av_array = np.sum(array)/aabb.getVolume() + G.graph_properties["vascular_tortuosity"] = G.new_graph_property("double", val= av_array) + f.write("%.15f\t" % av_array) + ############################Graph Statistics########################### + + # Some Graph Statistics are reported as values representing the graph + # Some statistics are reported as a function of volume. + array = G.vertex_properties["degree"].get_array() + av_array = np.sum(array)/G.num_vertices() + G.graph_properties["local_average_degree"] = G.new_graph_property("double", val = av_array) + f.write("%.15f\t" % av_array) + + array = G.vertex_properties["degree_volume"].get_array() + av_array = np.sum(array)/G.num_vertices() + G.graph_properties["local_average_degree_volume"] = G.new_graph_property("double", val = av_array) + f.write("%.15f\t" % av_array) + + array = G.vertex_properties["degree_tortuosity"].get_array() + av_array = np.sum(array)/G.num_vertices() + G.graph_properties["local_average_degree_tortuosity"] = G.new_graph_property("double", val = av_array) + f.write("%.15f\t" % av_array) + + array = G.vertex_properties["bc"].get_array() + av_array = np.sum(array)/G.num_vertices() + G.graph_properties["local_vertex_bc"] = G.new_graph_property("double", val = av_array) + f.write("%.15f\t" % av_array) + + array = G.vertex_properties["degree"].get_array() + av_array = np.sum(array)/aabb.getVolume() + G.graph_properties["local_vascular_degree"] = G.new_graph_property("double", val = av_array) + f.write("%.15f\t" % av_array) + + array = G.vertex_properties["degree_volume"].get_array() + av_array = np.sum(array)/aabb.getVolume() + G.graph_properties["local_vascular_degree_volume"] = G.new_graph_property("double", val = av_array) + f.write("%.15f\t" % av_array) + + array = G.vertex_properties["degree_tortuosity"].get_array() + av_array = np.sum(array)/aabb.getVolume() + G.graph_properties["local_vascular_degree_tortuosity"] = G.new_graph_property("double", val = av_array) + f.write("%.15f\t" % av_array) + + array = G.vertex_properties["bc"].get_array() + av_array = np.sum(array)/aabb.getVolume() + G.graph_properties["local_vascular_vertex_bc"] = G.new_graph_property("double", val = av_array) + f.write("%.15f\t" % av_array) + + array = G.edge_properties["bc"].get_array() + av_array = np.sum(array)/aabb.getVolume() + G.graph_properties["local_vascular_edge_bc"] = G.new_graph_property("double", val = av_array) + f.write("%.15f\t" % av_array) + + G.edge_properties["mst"] = gt.graph_tool.topology.min_spanning_tree(G) + G.graph_properties["local_mst_ratio"] = G.new_graph_property("double", val = np.double(np.sum(G.edge_properties["mst"].get_array()))/np.double(G.num_edges())) + f.write("%.15f\t" % G.graph_properties["local_mst_ratio"]) + + G.edge_properties["mst"] = gt.graph_tool.topology.min_spanning_tree(G, weights=G.edge_properties["volume"]) + G.graph_properties["local_mst_ratio_vascular_volume"] = G.new_graph_property("double", val = np.double(np.sum(G.edge_properties["mst"].get_array()))/np.double(G.num_edges())/aabb.getVolume()) + f.write("%.15f\t" % G.graph_properties["local_mst_ratio_vascular_volume"]) + + [lE, lV] = gt.graph_tool.centrality.eigenvector(G) + G.graph_properties["Largest Eigenvector"] = G.new_graph_property("double", val = lE) + f.write("%.15f\t" % lE) + ############################# hybrid metrics ########################## + self.add_rst_metric(G, is_dual=is_dual) + f.write("%.15f\t" % G.graph_properties["sensitivity"]) + f.write("%s" % label) + f.write("\n") + + f.close() + if(writeKey): + location = path + "/key.txt" + f = open(location, "w+") + f.write("%s\t" % "global_vertex_bc") + f.write("%s\t" % "global_vascular_vertex_bc") + f.write("%s\t" % "global_vascular_edge_bc") + f.write("%s\t" % "global_mst_ratio") + f.write("%s\t" % "global_mst_ratio_vascular_volume") + f.write("%s\t" % "average_length") + f.write("%s\t" % "average_volume") + f.write("%s\t" % "average_tortuosity") + f.write("%s\t" % "vascular_length") + f.write("%s\t" % "vascular_volume") + f.write("%s\t" % "vascular_tortuosity") + ############################Graph Statistics########################### + f.write("%s\t" % "local_average_degree") + f.write("%s\t" % "local_average_degree_volume") + f.write("%s\t" % "local_average_degree_tortuosity") + f.write("%s\t" % "local_vertex_bc") + f.write("%s\t" % "local_vascular_degree") + f.write("%s\t" % "local_vascular_degree_volume") + f.write("%s\t" % "local_vascular_degree_tortuosity") + f.write("%s\t" % "local_vascular_vertex_bc") + f.write("%s\t" % "local_vascular_edge_bc") + f.write("%s\t" % "local_mst_ratio") + f.write("%s\t" % "local_mst_ratio_vascular_volume") + f.write("%s\t" % "Largest Eigenvector") + f.write("%s\t" % "label") + f.write("%s\t" % "sensitivity") + f.close() + + + ''' + Creates a graph from a list of nodes and a list of edges. + Uses edge length as weight. + Returns a NetworkX Object. + ''' + def createLengthGraph(self): + G = nx.Graph() + for i in range(len(self.N)): + G.add_node(i, p=self.N[i].p) + for i in range(len(self.F)): + G.add_edge(self.F[i].v0, self.F[i].v1, weight = self.F[i].length()) + G[self.F[i].v0][self.F[i].v1]['pts'] = self.F[i].points + G[self.F[i].v0][self.F[i].v1]['rads'] = self.F[i].radii + + return G + + ''' + filters all the vertices that are close to the border with deg=1 + As well as affiliated edges. It is recommended thatthe graph have all the + graph metrics refreshed after filtering. + ''' + + def filterBorder(self, G, is_dual=False): + bb = AABB(G, is_dual) + TFv = G.new_vertex_property("bool", vals = np.full((G.num_vertices(),1), True, dtype=bool)) + TFe = G.new_edge_property("bool", vals = np.full((G.num_edges(),1), True, dtype=bool)) + if(is_dual==False): + for v in G.vertices(): + pt = np.array([G.vertex_properties["p"][v][0], G.vertex_properties["p"][v][1], G.vertex_properties["p"][v][2]]) + if G.vertex_properties["degree"][v] == 1 and bb.distance(pt) < 5.0: + TFv[v] = False + for e in v.all_edges(): + TFe[G.edge(e.source(), e.target())] = False + else: + #print("I have entered this method") + for v in G.vertices(): + l = len(G.vertex_properties["x"][v]) + pt_0 = np.array([G.vertex_properties["x"][v][0], G.vertex_properties["y"][v][0], G.vertex_properties["z"][v][0]]) + pt_1 = np.array([G.vertex_properties["x"][v][l-1], G.vertex_properties["y"][v][l-1], G.vertex_properties["z"][v][l-1]]) + if (G.vertex_properties["degree"][v] == 2): + if ((bb.distance(pt_0) < 5.0) or (bb.distance(pt_1) < 5.0)): + #print("length", l, ": ", G.vertex_properties["x"][v]) + #print("degree", G.vertex_properties["degree"][v]) + TFv[v] = False + for e in v.all_edges(): + TFe[G.edge(e.source(), e.target())] = False + + G.set_filters(TFe, TFv) + G1 = gt.Graph(G, directed=False, prune=True) + G.clear_filters() + G1.clear_filters() + + return G1 + + + ''' + Returns a Graph object that retains all the properties of the original + if return_largest is set as false, then all graphs with more than 20 nodes + are returned in a list. Metrics might have to be recalculated after + It is recommended that the graph have all the metrics necessary prior to filtering + In order to avoid internal boost index errors. + ''' + def filterDisconnected(self, G, return_largest=True): + #create masks + clusters = G.new_vertex_property("int", vals=np.full((G.num_vertices(), 1), -1, dtype=int)) + eclusters = G.new_edge_property("int", vals=np.full((G.num_edges(), 1), -1, dtype=int)) + c = 0 + + #run a bfs that sets visited nodes on each vertex + for i in G.vertices(): + if clusters[i] == -1: + gt.bfs_search(G, i, VisitorClassDisconnected(eclusters, clusters, c)) + c = c + 1 + + #find the number of clusters + unique, counts = np.unique(clusters.get_array(), return_counts=True) + eunique, ecounts = np.unique(clusters.get_array(), return_counts=True) + + if(return_largest == True): + #find the argument of the largest + a = counts.argmax() + b = ecounts.argmax() + + #turn the largest cluster into a TF graph mask for the vertices and the edges + TFv = G.new_vertex_property("bool", vals = np.full((G.num_vertices(), 1), False, dtype=bool)) + TFe = G.new_edge_property("bool", vals = np.full((G.num_edges(),1), False, dtype=bool)) + for i in G.vertices(): + if clusters[i] == unique[a]: + TFv[i] = True + + e = G.get_edges(); + for i in range(len(e)): + if eclusters[G.edge(G.vertex(e[i][0]), G.vertex(e[i][1]))] == eunique[b]: + TFe[G.edge(G.vertex(e[i][0]), G.vertex(e[i][1]))] = True + #set the filters and return the vertices + + G.set_filters(TFe, TFv) + G1=gt.Graph(G, directed=False, prune=True) + G.clear_filters() + G1.clear_filters() + return G1 + else: + Graphs = [] + #repeat the process above for each unique disconncted component. + for j in range(len(unique)): + if(counts[j] > 20): + TFv = G.new_vertex_property("bool", vals = np.full((G.num_vertices(),1), False, dtype=bool)) + TFe = G.new_edge_property("bool", vals = np.full((G.num_edges(),1), False, dtype=bool)) + for i in range(G.num_vertices()): + if clusters[G.vertex(i)] == unique[j]: + TFv[G.vertex(i)] = True + + e = G.get_edges(); + for i in range(len(e)): + if eclusters[G.edge(G.vertex(e[i][0]), G.vertex(e[i][1]))] == eunique[j]: + TFe[G.edge(G.vertex(e[i][0]), G.vertex(e[i][1]))] = True + G.set_filters(TFe, TFv) + G1=gt.Graph(G, directed=False, prune=True) + G.clear_filters() + G1.clear_filters() + Graphs.append(G1) + return Graphs + + + def simulate_fractures(self, G, remove=10): + num_removed = int(np.floor(G.num_edges()*remove/100)) + print("num of edges to begin with = ", G.num_edges()) + indices = np.random.randint(0, int(G.num_edges()), size = [num_removed,1]) + #aabb = AABB(G, is_dual) + tf = np.full((G.num_edges(), 1), True, dtype=bool) + for j in range(indices.shape[0]): + tf[indices[j]] = False + TFe = G.new_edge_property("bool", vals = tf) + G.set_edge_filter(TFe) + G1 = gt.Graph(G, directed=False, prune=True) + G.clear_filters() + G1.clear_filters() + G1 = self.filterDisconnected(G1) + G1.vertex_properties["degree"] = G1.degree_property_map("total") + print("num of edges left = ", G1.num_edges()) + print("I should have removed, ", num_removed) + print("Instead I removed, ", G.num_edges() - G1.num_edges()) + + return G1 + + ''' + Adds the mst based- sensitivity metrics. N is the number of random removals + remove is the percentage*100 of the vessels removed in each iteration. + The number of vessels removed will be floor(num_edges*remove/100) + ''' + def add_rst_metric(self, G, N=100, remove=10, is_dual=False): + #rst_ratio = G.new_graph_property("double") + rst_r = 0 + #G.nu + num_removed = int(np.floor(G.num_edges()*remove/100)) + indices = np.random.randint(0, int(G.num_edges()), size = [N, num_removed]) + aabb = AABB(G, is_dual) + + for i in range(N): + tf = np.full((G.num_edges(),1), True, dtype=bool) + for j in range(num_removed): + tf[indices[i, j]] = False +# rst = gt.graph_tool.topology.min_spanning_tree(G, weights = G.) +# ratio = np.double(np.sum(rst.get_array()))/np.double(G.num_edges()) +# rst_dist.append(ratio) +# rst_r = rst_r + ratio + TFe = G.new_edge_property("bool", vals = tf) + G.set_edge_filter(TFe) + G1 = gt.Graph(G, directed=False, prune=True) + G.clear_filters() + G1.clear_filters() + while(np.any(G1.degree_property_map("total").get_array() == True)): + #print(np.any(G1.degree_property_map("total").get_array()), " ", i) + G1 = self.filterDisconnected(G1) + G1 = self.filterBorder(G1, is_dual) + G1 = self.recalculate_metrics(G1, is_dual) + + if(not is_dual and G1.num_edges() > 2): + G1.edge_properties["mst"] = gt.graph_tool.topology.min_spanning_tree(G1, weights=G1.edge_properties["volume"]) + tvalues = G1.edge_properties["mst"].get_array() + #print(tvalues) + values = G1.edge_properties["inverted_volume"].get_array() + #print(values) + for k in range(G1.num_edges()): + if(tvalues[k] == True): + rst_r = rst_r + values[k] + elif(is_dual and G1.num_edges() > 2): + #This is non-sensical for now. + G1.edge_properties["mst"] = gt.graph_tool.topology.min_spanning_tree(G1) + tvalues = G1.vertex_properties["mst"].get_array() + values = G1.edge_properties["inverted_volume"].get_array() + for k in range(G1.num_edges()): + if(tvalues[k] == True): + rst_r = rst_r + values[k] + + G.graph_properties["sensitivity"] = G.new_graph_property("double", rst_r/N) + + return G + + ''' + Re-calculates connectivity based metrics and returns the resulting graph + ''' + def recalculate_metrics(self, G, is_dual=False, normalize=False): + gt.graph_tool.centrality.betweenness(G, vprop=G.vertex_properties["bc"], eprop=G.edge_properties["bc"], norm=normalize) + if(is_dual == False): + #gt.graph_tool.centrality.betweenness(G, G.vertex_properties["bc"], G.edge_properties["bc"], norm=True, is_dual=False) + #generate minimum spanning tree + + G.edge_properties["mst"] = gt.graph_tool.topology.min_spanning_tree(G, weights=G.edge_properties["length"]) + G.graph_properties["mst_ratio"] = np.double(np.sum(G.edge_properties["mst"].get_array()))/np.double(G.num_edges()) + print(np.double(np.sum(G.edge_properties["mst"].get_array()))/np.double(G.num_edges())) + G.vertex_properties["degree"] = G.degree_property_map("total") + G.vertex_properties["degree_volume"] = G.degree_property_map("total", weight=G.edge_properties["volume"]) + G.vertex_properties["degree_tortuosity"] = G.degree_property_map("total", G.edge_properties["tortuosity"]) + + else: + #gt.graph_tool.centrality.betweenness(G, G.vertex_properties["bc"], G.edge_properties["bc"], norm=True, is_dual=False) + G.edge_properties["mst"] = gt.graph_tool.topology.min_spanning_tree(G, weights=G.edge_properties["bc"]) #Boolean value deligating belongance to the minimum spanning tree. + G.graph_properties["mst_ratio"] = np.double(np.sum(G.edge_properties["mst"].get_array()))/np.double(G.num_edges()) + G.vertex_properties["degree"] = G.degree_property_map("total") + G.vertex_properties["degree_centrality"] = G.degree_property_map("total", weight=G.edge_properties["bc"]) + + return G + + ''' + Filter the graph to remove the disconnected components of the graph if + disconnected is set to true. if borders are set to True, the algorithm + cleans up the nodes with degree one that are close to the edge of the volume. + Returns a Graph object with updated internal property maps. + + if return_largerst == False, will return an array of graphs. + disconnected == apply the disconnected componenet filter + return_largest == if false returns the largest set of graphs (n > 10) + borders == True, if true applies the border (deg == 1 near border) filter + add_rst == True, if true adds the rst_metric + ''' + def filterFullGraph_gt(self, G, disconnected=True, return_largest=True, borders=True, erode = False, add_rst=False, dual=False): + if(disconnected==True): + if(return_largest==True and borders==True): + G1 = self.filterDisconnected(G, return_largest) + G1 = self.recalculate_metrics(G1, is_dual=dual) + G1 = self.filterBorder(G1, is_dual=dual) + if(erode == True): + while(np.any(G1.degree_property_map("total").get_array() == True)): + G1 = self.filterBorder(G1, is_dual=dual) + G1 = self.recalculate_metrics(G1, is_dual=dual) + else: + G1 = self.recalculate_metrics(G1, is_dual=dual) + if(add_rst == True): + G1 = self.add_rst_metric(G1) + elif(return_largest==True and borders==False): + G1 = self.filterDisconnected(G, return_largest) + if(erode == True): + while(np.any(G1.degree_property_map("total").get_array() == True)): + G1 = self.filterBorder(G1, is_dual=dual) + G1 = self.recalculate_metrics(G1, is_dual=dual) + else: + G1 = self.recalculate_metrics(G1, is_dual=dual) + if(add_rst == True): + G1 = self.add_rst_metric(G1) + elif(return_largest==False and borders == True): + G1 = self.filterDisconnected(G, return_largest) + for graph in G1: + graph = self.recalculate_metrics(graph, is_dual=dual) + graph = self.filterBorder(graph, is_dual=dual) + if(erode == True): + while(np.any(graph.degree_property_map("total").get_array() == True)): + graph = self.filterBorder(graph, is_dual=dual) + graph = self.recalculate_metrics(graph, is_dual=dual) + else: + graph = self.recalculate_metrics(graph, is_dual=dual) + if(add_rst == True): + graph = self.add_rst_metric(graph) + else: + G1 = self.filterDisconnected(G, return_largest) + for graph in G1: + graph = self.recalculate_metrics(graph, is_dual=dual) + if(add_rst == True): + graph = self.add_rst_metric(graph) + else: + G1 = self.filterBorder(G, is_dual=dual); + if(erode == True): + while(np.any(G1.degree_property_map("total").get_array() == True)): + print(G1.num_vertices()) + G1 = self.filterBorder(G1, is_dual=dual) + print(G1.num_vertices()) + G1 = self.recalculate_metrics(G1, is_dual=dual, ) + else: + G1 = self.recalculate_metrics(G1, is_dual=dual) + if(add_rst==True): + G1 = self.add_rst_metric(G1) + + return G1 + + ''' + Accretes the graph in G using the nodes and edges in G_0 + Assumes that G_0 has all the same properties as G, including "idx" + idx is an immutable set of indices that do not change when the graph is + modified + ''' + def accrete(self, G, G_0): + + v_filters = G_0.new_vertex_property("bool", vals = np.full((G_0.num_vertices(),1), False, dtype=bool)) + e_filters = G_0.new_edge_property("bool", vals = np.full((G_0.num_edges(),1), False, dtype=bool)) + for v in G.vertices(): + #print(v) + G_0_vertex = G_0.vertex(G.vertex_properties["idx"][v]) + v_filters[G_0_vertex] = True + for e in G_0_vertex.all_edges(): + e_filters[e] = True + for v_0 in G_0_vertex.all_neighbors(): + v_filters[v_0] = True + + G_0.set_filters(e_filters, v_filters) + + graph = gt.Graph(G_0, directed=False, prune=True) + G_0.clear_filters() + + return graph + + + + + ''' + Creates a graph from a list of nodes and a list of edges + The graph is the "inverse" of the original graph, + Vertices are edges and edges are + + ''' + def createDualGraph_gt(self): + #Generate the undirected graph + G = gt.Graph() + G.set_directed(False) + + #add all the required vertex properties + vpos = G.new_vertex_property("vector") #original position of the edge (placed as midpoin of edge) + pos = G.new_vertex_property("vector") #positions according to the fdl + rpos = G.new_vertex_property("vector") #positions according to the radial layout + x_edge = G.new_vertex_property("vector") #x-coordinated of the edge + y_edge = G.new_vertex_property("vector") #y-coordinates of the edge + z_edge = G.new_vertex_property("vector") #z-coordinates of the edge + r_edge = G.new_vertex_property("vector") #r-values at each x,y,z + l_edge = G.new_vertex_property("double") #length of the edge + i_edge = G.new_vertex_property("double") #1/length + iv_edge = G.new_vertex_property("double") #1/volume + t_edge = G.new_vertex_property("double") #tortuosity of the edge + v_edge = G.new_vertex_property("double") #volume of the edge + vbetweeness_centrality = G.new_vertex_property("double") #betweeneness centrality of the vertex + gaussian = G.new_vertex_property("double") #empty property map for gaussian values downstream + degree = G.new_vertex_property("int") + degree_centrality = G.new_vertex_property("int") + + #add all the required edge properties + ebetweeness_centrality = G.new_edge_property("double") #betweeness centrality of the edge + mst = G.new_edge_property("bool") #Boolean value deligating belongance to the minimum spanning tree. + l_edge = G.new_edge_property("double") #length of the edge + gaussian = G.new_edge_property("double") + + + #add all graph properties + mst_ratio = G.new_graph_property("double") + + #add all the edges as vertices + for i in range(len(self.F)): + if(self.F[i].length() > 0): + G.add_vertex() + x,y,z,r = self.F[i].getcoords() + x_edge[G.vertex(i)] = x + y_edge[G.vertex(i)] = y + z_edge[G.vertex(i)] = z + r_edge[G.vertex(i)] = r + l_edge[G.vertex(i)] = self.F[i].length() + gaussian[G.vertex(i)] = np.exp(-np.power(l_edge[G.vertex(i)], 2.) / (2 * np.power(60.0, 2.))) + i_edge[G.vertex(i)] = 1/l_edge[G.vertex(i)] + t_edge[G.vertex(i)] = self.F[i].turtuosity() + v_edge[G.vertex(i)] = self.F[i].volume() + iv_edge[G.vertex(i)] = 1/v_edge[G.vertex(i)] + vpos[G.vertex(i)] = np.array([x[int(len(x)/2)], y[int(len(x)/2)], z[int(len(x)/2)]], dtype=float) + + #based on the neighbors add edges. If two edges share a vertex, there is a vertex between them. + for i in range(len(self.N)): + #in case there are 3 or more edges attached to a vertex + if(len(self.N[i].o) > 2): + for j in range(len(self.N[i].o)-1): + G.add_edge(G.vertex(self.N[i].o[j]), G.vertex(self.N[i].o[j+1])) + l_edge[G.edge(G.vertex(self.N[i].o[j])), G.vertex(self.N[i].o[j+1])] = 1 + gaussian[G.edge(G.vertex(self.N[i].o[j])), G.vertex(self.N[i].o[j+1])] = \ + np.exp(-np.power(l_edge[G.edge(G.vertex(self.N[i].o[j]),G.vertex(self.N[i].o[j+1]))], 2.) / (2 * np.power(60.0, 2.))) + #add the final edge. + G.add_edge(G.vertex(self.N[i].o[0]), G.vertex(self.N[i].o[len(self.N[i].o)-1])) + l_edge[G.edge(G.vertex(self.N[i].o[j])), G.vertex(self.N[i].o[j+1])] = 1 + gaussian[G.edge(G.vertex(self.N[i].o[j])), G.vertex(self.N[i].o[j+1])] = \ + np.exp(-np.power(l_edge[G.edge(G.vertex(self.N[i].o[j]),G.vertex(self.N[i].o[j+1]))], 2.) / (2 * np.power(60.0, 2.))) + #in case there are 2 edges attached to a vertex + elif(len(self.N[i].o) == 2): + G.add_edge(G.vertex(self.N[i].o[0]), G.vertex(self.N[i].o[1]-1)) + l_edge[G.edge(G.vertex(self.N[i].o[j])), G.vertex(self.N[i].o[j+1])] = 1 + gaussian[G.edge(G.vertex(self.N[i].o[j])), G.vertex(self.N[i].o[j+1])] = \ + np.exp(-np.power(l_edge[G.edge(G.vertex(self.N[i].o[j]),G.vertex(self.N[i].o[j+1]))], 2.) / (2 * np.power(60.0, 2.))) + #in case there is only a single edge coming into the vertex we do + # not add anything + + #generate centrality map + gt.graph_tool.centrality.betweenness(G, vbetweeness_centrality, ebetweeness_centrality, norm=True) + #generate minimum spanning tree + mst = gt.graph_tool.topology.min_spanning_tree(G, weights=ebetweeness_centrality) + mst_ratio[G] = np.double(np.sum(mst.get_array()))/np.double(len(self.N)) + degree = G.degree_property_map("total") + degree_centrality = G.degree_property_map("total", weight=ebetweeness_centrality) + #degree_tortuosity = G.degree_property_map("total", weight=t_edge) + + pos = gt.sfdp_layout(G) + + rpos = gt.radial_tree_layout(G, root=G.vertex(np.argmax(vbetweeness_centrality.get_array()))) + + #set property maps for the vertices + G.vertex_properties["p"] = vpos + G.vertex_properties["pos"] = pos + G.vertex_properties["rpos"] = rpos + G.vertex_properties["bc"] = vbetweeness_centrality + G.vertex_properties["degree_centrality"] = degree_centrality + G.vertex_properties["degree"] = degree + G.vertex_properties["x"] = x_edge + G.vertex_properties["y"] = y_edge + G.vertex_properties["z"] = z_edge + G.vertex_properties["r"] = r_edge + G.vertex_properties["inverted"] = i_edge + G.vertex_properties["inverted_volume"] = iv_edge + G.vertex_properties["length"] = l_edge + G.vertex_properties["tortuosity"] = t_edge + G.vertex_properties["volume"] = v_edge + G.vertex_properties["bc"] = vbetweeness_centrality + + #G.vertex_properties["pos"] = gt.fruchterman_reingold_layout(G, weight=G.edge_properties["length"], circular = True, n_iter= 10000, pos=G.vertex_properties["pos"], r = 2.0) + G.vertex_properties["pos"] = gt.sfdp_layout(G, C = 0.5, p = 3.0) + #set property maps for the edges + G.edge_properties["bc"] = ebetweeness_centrality + G.edge_properties["mst"] = mst + + #set graph properies + G.graph_properties["mst_ratio"] = mst_ratio + + + +# title = "raw_dual_cortical_7_6_before_delete.pdf" +# title2 = "raw_dual_radial_cordical_7_6_before_delete.pdf" +# +# gt.graph_draw(G, pos=gt.radial_tree_layout(G, root=G.vertex(np.argmax(vbetweeness_centrality.get_array()))), edge_pen_width = 2.0, vertex_size=degree, edge_color=mst, vertex_fill_color=vbetweeness_centrality, output=title2, bg_color=[0.0, 0.0,0.0,1.0], vertex_text=G.vertex_index, output_size=(1000,1000),vertex_font_size = 6) +# gt.graph_draw(G, pos=pos, edge_pen_width = 4.0, vertex_size=degree, edge_color=mst, vertex_fill_color=vbetweeness_centrality, output=title, bg_color=[0.0, 0.0,0.0,1.0], vertex_text=G.vertex_index, output_size=(1000,1000),vertex_font_size = 6) +# + title = "./Original_2D_Graph.png" +# title2 = "raw_dual_radial_cordical_7_6_after_delete.pdf" +# + G1 = self.filterFullGraph_gt(G, borders=False, dual=True) + +# #print(clusters.get_array()[:], clusters.get_array()[:]) +# #pos = gt.sfdp_layout(G) +# +# #gt.graph_draw(G, pos=,edge_color=mst, vertex_fill_color=vbetweeness_centrality, edge_pen_width=ebetweeness_centrality, output="raw.pdf") +# +# gt.graph_draw(G1, pos=gt.radial_tree_layout(G1, root=G1.vertex(np.argmax(G1.vertex_properties["bc"].get_array()))), edge_pen_width = 4.0, vertex_size=G1.vertex_properties["degree"], edge_color=G1.edge_properties["mst"], vertex_fill_color=G1.vertex_properties["bc"], output=title2, bg_color=[0.0, 0.0,0.0,1.0], vertex_text=G1.vertex_index, output_size=(1000,1000),vertex_font_size = 6) +# gt.graph_draw(G1, pos=G1.vertex_properties["p"], edge_pen_width = 4.0, vertex_size=G1.vertex_properties["degree"], edge_color=G1.edge_properties["mst"], vertex_fill_color=G1.vertex_properties["bc"], output=title, bg_color=[0.0, 0.0,0.0,1.0], vertex_text=G1.vertex_index, output_size=(1000,1000),vertex_font_size = 6) + gt.graph_draw(G1, pos=G1.vertex_properties["p"], edge_pen_width = 4.0, vertex_size=G1.vertex_properties["degree"], output=title, bg_color=[1.0, 1.0,1.0,1.0], vertex_text=G1.vertex_index, output_size=(1000,1000),vertex_font_size = 6) + title = "./Original_2D_Layout.png" + gt.graph_draw(G1, pos=G1.vertex_properties["p"], edge_pen_width = 4.0, vertex_size=G1.vertex_properties["degree"], edge_color=G1.edge_properties["mst"], vertex_fill_color=G1.vertex_properties["bc"], output=title, bg_color=[1.0, 1.0,1.0,1.0], vertex_text=G1.vertex_index, output_size=(1000,1000), vertex_font_size = 6) +#gt.graph_draw(G1, pos=G1.vertex_properties["pos"], edge_pen_width = 4.0, vertex_size=degree, edge_color=G1.edge_properties["mst"], vertex_fill_color=G1.vertex_properties["bc"], output=title) +# +# +# G2 = self.filterFullGraph_gt(G1, add_rst=False, dual=True) +# title = "raw_dual_cortical_7_6_after_borders.pdf" +# title2 = "raw_dual_radial_cortical_7_6_after_borders.pdf" +# gt.graph_draw(G2, pos=G2.vertex_properties["pos"], edge_pen_width = 4.0, vertex_size=G2.vertex_properties["degree"], edge_color=G2.edge_properties["mst"], vertex_fill_color=G2.vertex_properties["bc"], output=title, bg_color=[0.0, 0.0,0.0,1.0], vertex_text=G2.vertex_index, output_size=(1000,1000),vertex_font_size = 6) +# gt.graph_draw(G2, pos=gt.radial_tree_layout(G2, root=G2.vertex(np.argmax(G2.vertex_properties["bc"].get_array()))), edge_pen_width = 4.0, vertex_size=G2.vertex_properties["degree"], edge_color=G2.edge_properties["mst"], vertex_fill_color=G2.vertex_properties["bc"], output=title2, bg_color=[0.0, 0.0,0.0,1.0], vertex_text=G2.vertex_index, output_size=(1000,1000),vertex_font_size = 6) + return G + + + ''' + Create a graph from a list of nodes and a list of edges. + Populates the Graph with a bunch of property maps in respect to edges and nodes + each vertex is assigned 3 positional coordiantes: + The original x,y,z location of the graph. + x,y locations according to force directed layout, + and x,y positions according to a radial layout + Returns a graph_tool object + ''' + def createFullGraph_gt(self): + #Generate the undirected graph + G = gt.Graph() + G.set_directed(False) + + #add all the required vertex properties + vpos = G.new_vertex_property("vector") #original vertex positions + pos = G.new_vertex_property("vector") #positions according to the fdl + rpos = G.new_vertex_property("vector") #positions according to the radial layout + vbetweeness_centrality = G.new_vertex_property("double") #betweeneness centrality of the vertex + degree = G.new_vertex_property("int") #degree + degree_volume = G.new_vertex_property("double") #degree scaled by the volume of all in fibers + degree_tortuosity = G.new_vertex_property("double") #degree scaled by the tortuosity of all in fibers. + + + #add all the required edge properties + #G.properties[("x,y,z,r")] = G.new_edge_property("vector") #x-coordinated of the edge + y_edge = G.new_edge_property("vector") #y-coordinates of the edge + z_edge = G.new_edge_property("vector") #z-coordinates of the edge + r_edge = G.new_edge_property("vector") #r-values at each x,y,z + av_edge = G.new_edge_property("double") #average edge radius + l_edge = G.new_edge_property("double") #length of the edge + i_edge = G.new_edge_property("double") #1/length + iv_edge = G.new_edge_property("double") #1/volume + t_edge = G.new_edge_property("double") #tortuosity of the edge + v_edge = G.new_edge_property("double") #volume of the edge + ebetweeness_centrality = G.new_edge_property("double") #betweeness centrality of the edge + ebc_length = G.new_edge_property("double") #betweeneness centrality scaled by the length + mst = G.new_edge_property("bool") #Boolean value deligating belongance to the minimum spanning tree. + + + #This map gets updated. + gaussian = G.new_edge_property("double") #empty property map for gaussian values downstream + + + #Graph properties (once per region) + mst_ratio = G.new_graph_property("double") #Graph based metric of mst edges/total edges. + + + #add verticies and set their properties + for i in range(len(self.N)): + G.add_vertex() + vpos[G.vertex(i)] = self.N[i].p + + #add edges and set their properties. + for i in range(len(self.F)): + v0 = self.F[i].v0 + v1 = self.F[i].v1 + if self.F[i].length() > 0.0: + x,y,z,r = self.F[i].getcoords() + G.add_edge(G.vertex(v0), G.vertex(v1)) + l_edge[G.edge(v0,v1)] = self.F[i].length() + gaussian[G.edge(v0,v1)] = np.exp(-np.power(l_edge[G.edge(v0,v1)], 2.) / (2 * np.power(60.0, 2.))) + i_edge[G.edge(v0,v1)] = 1/l_edge[G.edge(v0,v1)] + t_edge[G.edge(v0,v1)] = self.F[i].turtuosity() + v_edge[G.edge(v0,v1)] = self.F[i].volume() + iv_edge[G.edge(v0,v1)] = 1/v_edge[G.edge(v0,v1)] + x_edge[G.edge(v0,v1)] = x + y_edge[G.edge(v0,v1)] = y + z_edge[G.edge(v0,v1)] = z + r_edge[G.edge(v0,v1)] = r + av_edge[G.edge(v0,v1)] = self.F[i].av_radius() + else: + print("WTF") + + + #generate centrality map + gt.graph_tool.centrality.betweenness(G, vprop=vbetweeness_centrality, eprop=ebetweeness_centrality, norm=True) + #generate minimum spanning tree + mst = gt.graph_tool.topology.min_spanning_tree(G, weights=l_edge) + mst_ratio[G] = np.double(np.sum(mst.get_array()))/np.double(len(self.N)) + degree = G.degree_property_map("total") + degree_volume = G.degree_property_map("total", weight=v_edge) + + dg = degree_volume.get_array() + dg = np.exp(-np.power(dg, 2.) / (2 * np.power(0.000005, 2.))) + degree_volume = G.new_vertex_property("double", vals=dg) + degree_tortuosity = G.degree_property_map("total", weight=t_edge) + + + #print(clusters.get_array()[:], clusters.get_array()[:]) + pos = gt.sfdp_layout(G, C = 1.0, K = 10) + rpos = gt.radial_tree_layout(G, root=G.vertex(np.argmax(vbetweeness_centrality.get_array()))) + + + + for e in G.edges(): + ebc_length[G.edge(e.source(), e.target())] = ebetweeness_centrality[G.edge(e.source(), e.target())]*l_edge[G.edge(e.source(), e.target())] + + #set property maps for the vertices + G.vertex_properties["p"] = vpos + G.vertex_properties["pos"] = pos + G.vertex_properties["rpos"] = rpos + G.vertex_properties["bc"] = vbetweeness_centrality + G.vertex_properties["degree"] = degree + G.vertex_properties["degree_volume"] = degree_volume + G.vertex_properties["degree_tortuosity"] = degree_tortuosity + G.vertex_properties["selection"] = G.new_vertex_property("bool", val=False) + + #set property maps for the edges + G.edge_properties["x"] = x_edge + G.edge_properties["y"] = y_edge + G.edge_properties["z"] = z_edge + G.edge_properties["r"] = r_edge + G.edge_properties["inverted"] = i_edge + G.edge_properties["length"] = l_edge + G.edge_properties["tortuosity"] = t_edge + G.edge_properties["av_radius"] = av_edge + G.edge_properties["volume"] = v_edge + G.edge_properties["bc"] = ebetweeness_centrality + G.edge_properties["bc*length"] = ebc_length + G.edge_properties["mst"] = mst + G.edge_properties["inverted_volume"] = iv_edge + G.edge_properties["gaussian"] = gaussian + G.edge_properties["selection"] = G.new_edge_property("double", val=0.0) + + #set graph properies + G.graph_properties["mst_ratio"] = mst_ratio + + N = 0 + for e in G.edges(): + N += len(G.edge_properties["x"][e]) + point_cloud_x = np.zeros((N, 1), dtype=np.float64) + point_cloud_y = np.zeros((N, 1), dtype=np.float64) + point_cloud_z = np.zeros((N, 1), dtype=np.float64) + + n = 0 + for e in G.edges(): + for p in range(len(G.edge_properties["x"][e])): + point_cloud_x[n][0] = G.edge_properties["x"][e][p] + point_cloud_y[n][0] = G.edge_properties["y"][e][p] + point_cloud_z[n][0] = G.edge_properties["z"][e][p] + n += 1 + + G.graph_properties["point_cloud_x"] = G.new_graph_property("vector") + G.graph_properties["point_cloud_y"] = G.new_graph_property("vector") + G.graph_properties["point_cloud_z"] = G.new_graph_property("vector") + G.graph_properties["point_cloud_x"] = point_cloud_x + G.graph_properties["point_cloud_y"] = point_cloud_y + G.graph_properties["point_cloud_z"] = point_cloud_z + #save the original graph indices for the edges and the vertices. + G.vertex_properties["idx"] = G.vertex_index + G.edge_properties["idx"] = G.edge_index + + title = "~/Pictures/Original_2D_Graph_2.png" +# title2 = "raw_dual_radial_cordical_7_6_after_delete.pdf" +# +# G1 = self.filterFullGraph_gt(G, borders=False, dual=False, add_rst = False) +## +## #print(clusters.get_array()[:], clusters.get_array()[:]) +## #pos = gt.sfdp_layout(G) +## +## #gt.graph_draw(G, pos=,edge_color=mst, vertex_fill_color=vbetweeness_centrality, edge_pen_width=ebetweeness_centrality, output="raw.pdf") +## +## gt.graph_draw(G1, pos=gt.radial_tree_layout(G1, root=G1.vertex(np.argmax(G1.vertex_properties["bc"].get_array()))), edge_pen_width = 4.0, vertex_size=G1.vertex_properties["degree"], edge_color=G1.edge_properties["mst"], vertex_fill_color=G1.vertex_properties["bc"], output=title2, bg_color=[0.0, 0.0,0.0,1.0], vertex_text=G1.vertex_index, output_size=(1000,1000),vertex_font_size = 6) +## gt.graph_draw(G1, pos=G1.vertex_properties["p"], edge_pen_width = 4.0, vertex_size=G1.vertex_properties["degree"], edge_color=G1.edge_properties["mst"], vertex_fill_color=G1.vertex_properties["bc"], output=title, bg_color=[0.0, 0.0,0.0,1.0], vertex_text=G1.vertex_index, output_size=(1000,1000),vertex_font_size = 6) +# gt.graphviz_draw(G1, pos = G1.vertex_properties["p"], ratio="fill", size = (100,100), layout = None, pin = True, overlap=True, vsize=(G1.vertex_properties["degree_tortuosity"], 2.), vprops={"index":G1.vertex_index}, penwidth= 8.0, vcolor = G1.vertex_properties["bc"], vcmap = cm.get_cmap("inferno"), output=title) +# #gt.graph_draw(G1, pos=G1.vertex_properties["p"], edge_pen_width = 8.0, output=title, bg_color=[1.0, 1.0,1.0,1.0], vertex_size=60, vertex_fill_color=G1.vertex_properties["bc"], vertex_text=G1.vertex_index, output_size=(3200,3200),vertex_font_size = 32) +# title = "~/Pictures/Original_2D_Layout.png" +# gt.graph_draw(G1, pos=G1.vertex_properties["pos"], edge_pen_width = 8.0, edge_color=G1.edge_properties["mst"], vertex_size=40, vertex_fill_color=G1.vertex_properties["bc"], output=title, bg_color=[1.0, 1.0,1.0,1.0], vertex_text=G1.vertex_index, output_size=(3200,3200), vertex_font_size = 32) +# +# +# G1 = self.filterFullGraph_gt(G, borders=True, dual=False, erode=True, add_rst = False) +# title = "~/Pictures/Original_2D_Graph_Eroded.png" +# gt.graph_draw(G1, pos=G1.vertex_properties["p"], edge_pen_width = 8.0, output=title, bg_color=[1.0, 1.0,1.0,1.0], vertex_size=60, vertex_fill_color=G1.vertex_properties["bc"], vertex_text=G1.vertex_index, output_size=(3200,3200),vertex_font_size = 32) +# title = "~/Pictures/Original_2D_Layout_Eroded.png" +# gt.graph_draw(G1, pos=G1.vertex_properties["pos"], edge_pen_width = 8.0, edge_color=G1.edge_properties["mst"], vertex_size=40, vertex_fill_color=G1.vertex_properties["bc"], output=title, bg_color=[1.0, 1.0,1.0,1.0], vertex_text=G1.vertex_index, output_size=(3200,3200), vertex_font_size = 32) + +# title = "raw_full_cortical_7_6_before_delete.pdf" +# title2 = "raw_full_radial_cordical_7_6_before_delete.pdf" +# #G.graph_properties["rst_ratio"] = rst_ratio +# #print("") +# #print(G.graph_properties["rst_ratio"]) +# #print(G.edge_properties["bc"].get_array()) +# #print() +# +# #G.vertex_properties[("p","betweeness_centrality")] = [vpos, vbetweeness_centrality] +# #G.edge_properties[("x", "y", "z", "r", "length", "tortuosity", "volume", "betweeness_centrality")] = [] +# #G.set_edge_filter(mst) +# #G.set_edge_filter(mst) +# +# #Experimental code for drawing graphs. +# #gt.graph_draw(G,pos=vpos,edge_color=mst, vertex_fill_color=vbetweeness_centrality, edge_pen_width=ebetweeness_centrality, output="raw.pdf") +# #gt.graph_draw(G,pos=pos, vertex_size=degree, edge_color=mst, vertex_fill_color=vbetweeness_centrality, output=title, bg_color=[1.0,1.0,1.0]) +# #gt.graph_draw(G,pos=gt.radial_tree_layout(G, root=G.vertex(np.argmax(vbetweeness_centrality.get_array())), node_weight=vbetweeness_centrality), vertex_size=degree, edge_color=mst, vertex_fill_color=vbetweeness_centrality, output=title2, bg_color=[1.0, 1.0,1.0,1.0], vertex_text=G.vertex_index, output_size=(1000,1000), vertex_font_size = 6) +# #print(np.argmax(vbetweeness_centrality.get_array())) +# +# #np.savetxt("vbetweeness_centrality", vbetweeness_centrality.get_array()) +# gt.graph_draw(G,pos=gt.radial_tree_layout(G, root=G.vertex(np.argmax(vbetweeness_centrality.get_array()))), edge_pen_width = 2.0, vertex_size=degree, edge_color=mst, vertex_fill_color=vbetweeness_centrality, output=title2, bg_color=[0.0, 0.0,0.0,1.0], vertex_text=G.vertex_index, output_size=(1000,1000),vertex_font_size = 6) +# gt.graph_draw(G, pos=pos, edge_pen_width = 4.0, vertex_size=degree, edge_color=mst, vertex_fill_color=vbetweeness_centrality, output=title, bg_color=[0.0, 0.0,0.0,1.0], vertex_text=G.vertex_index, output_size=(1000,1000),vertex_font_size = 6) +# +# +# #gt.graph_draw(G,pos=vpos, efilt=mst, vertex_fill_color=vbetweeness_centrality, edge_pen_width=ebetweeness_centrality, output="mst.pdf") +# #gt.graph_draw(G,pos=pos, efilt=mst, vertex_fill_color=vbetweeness_centrality, output="mst.pdf") +# #need to create subgraphs!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +# #rst_r = 0 +# #for i in range(1000): +# # rst = gt.graph_tool.topology.random_spanning_tree(G) +# # rst_r = rst_r + np.double(np.sum(rst.get_array()))/np.double(len(self.N)) +# +# #rst_ratio[G] = rst_r/1000.0 +# +# +# +# #G1 = gt.Graph(G) +# +# #G1.purge_edges() +# #G1.purge_vertices(in_place = True) +# +# +# title = "raw_full_cortical_7_6_after_delete.pdf" +# title2 = "raw_full_radial_cordical_7_6_after_delete.pdf" +# +# G1 = self.filterFullGraph_gt(G,borders=False) +# +# #print(clusters.get_array()[:], clusters.get_array()[:]) +# #pos = gt.sfdp_layout(G) +# +# #gt.graph_draw(G,pos=vpos,edge_color=mst, vertex_fill_color=vbetweeness_centrality, edge_pen_width=ebetweeness_centrality, output="raw.pdf") +# +# gt.graph_draw(pos = gt.sfdp_layout(G1), output=title) +# gt.graph_draw(G1, pos=G1.vertex_properties["pos"], edge_pen_width = 4.0, vertex_size=G1.vertex_properties["degree"], edge_color=G1.edge_properties["mst"], vertex_fill_color=G1.vertex_properties["bc"], output=title, bg_color=[0.0, 0.0,0.0,1.0], vertex_text=G1.vertex_index, output_size=(1000,1000),vertex_font_size = 6) +# #gt.graph_draw(G1, pos=G1.vertex_properties["pos"], edge_pen_width = 4.0, vertex_size=degree, edge_color=G1.edge_properties["mst"], vertex_fill_color=G1.vertex_properties["bc"], output=title) +# +# +# G2 = self.filterFullGraph_gt(G1) +# title = "raw_full_cortical_7_6_after_borders.pdf" +# title2 = "raw_full_radial_cortical_7_6_after_borders.pdf" +# gt.graph_draw(G2, pos=G2.vertex_properties["pos"], edge_pen_width = 4.0, vertex_size=G2.vertex_properties["degree"], edge_color=G2.edge_properties["mst"], vertex_fill_color=G2.vertex_properties["bc"], output=title, bg_color=[0.0, 0.0,0.0,1.0], vertex_text=G2.vertex_index, output_size=(1000,1000),vertex_font_size = 6) +# gt.graph_draw(G2, pos=gt.radial_tree_layout(G2, root=G2.vertex(np.argmax(G2.vertex_properties["bc"].get_array()))), edge_pen_width = 4.0, vertex_size=G2.vertex_properties["degree"], edge_color=G2.edge_properties["mst"], vertex_fill_color=G2.vertex_properties["bc"], output=title2, bg_color=[0.0, 0.0,0.0,1.0], vertex_text=G2.vertex_index, output_size=(1000,1000),vertex_font_size = 6) +# #gt.graph_draw(G,pos=pos, vertex_size=degree, edge_color=mst, vertex_fill_color=vbetweeness_centrality, output=title, bg_color=[1.0,1.0,1.0]) +# +# #gt.graph_draw(G,pos=gt.radial_tree_layout(G, root=G.vertex(np.argmax(vbetweeness_centrality.get_array())), node_weight=vbetweeness_centrality), vertex_size=degree, edge_color=mst, vertex_fill_color=vbetweeness_centrality, output=title2, bg_color=[1.0, 1.0,1.0,1.0], vertex_text=G.vertex_index, output_size=(1000,1000), vertex_font_size = 6) +# #print(np.argmax(vbetweeness_centrality.get_array())) +# #np.savetxt("vbetweeness_centrality", vbetweeness_centrality.get_array()) +# #gt.graph_draw(G1, pos=gt.radial_tree_layout(G1, root=G1.vertex(np.argmax(G1.vertex_properties["bc"].get_array()))), edge_pen_width = 2.0, vertex_size=G1.vertex_properties["degree"], edge_color=G1.edge_properties["mst"], vertex_fill_color=G1.vertex_properties["bc"], output=title2, bg_color=[0.0, 0.0,0.0,1.0], vertex_text=G1.vertex_index, output_size=(1000,1000),vertex_font_size = 6) +# #gt.graph_draw(G1, pos=gt.radial_tree_layout(G1, root=G1.vertex(np.argmax(G1.vertex_properties["bc"].get_array()))), edge_pen_width = 2.0, vertex_size=degree, edge_color=G1.edge_properties["mst"]) +# #gt.draw_hierarchy(gt.minimize_nested_blockmodel_dl(G1, deg_corr=False), output = "test.pdf") +# #gt.graph_draw(G1,pos=gt.radial_tree_layout(G1, root=G1.vertex(np.argmax(G1.edge_properties["mst"].get_array()))), edge_pen_width = 2.0, vertex_size=degree, edge_color=mst, vertex_fill_color=vbetweeness_centrality, output=title2, bg_color=[0.0, 0.0,0.0,1.0], vertex_text=G.vertex_index, output_size=(1000,1000),vertex_font_size = 6) +# print(gt.graph_tool.openmp_enabled(), gt.graph_tool.openmp_get_num_threads()) + + #print(ratio, G2.graph_properties["mst_ratio"]) + #n, bins, patches = plt.hist(np.asarray(rst_dist), 50, normed=1, facecolor='green', alpha=0.75) + + # add a 'best fit' line + #plt.grid(True) + + #plt.show() + #self.partition(G, 5, 7, False) + self.export_points(G) + self.export_bb(G) + return G + + + def gen_new_fd_layout(G): + pos = gt.sfdp_layout(G, C = 1.0, K = 10) + G.vertex_properties["pos"] = pos + return G + + def map_edges_to_range(G, rng, propertymap): + def func(maximum, minimum, new_maximum, new_minimum, value): + return ((value-minimum)/(maximum-minimum)*(new_maximum-new_minimum)+new_minimum) + + mx = max(G.edge_properties[propertymap].get_array()) + mn = min(G.edge_properties[propertymap].get_array()) + G.edge_properties["map"] = G.new_edge_property("float") + gt.map_property_values(G.edge_properties[propertymap], G.edge_properties["map"], lambda x: func(mx, mn, rng[0], rng[1], x)) + + return G.edge_properties["map"] + + + def map_vertices_to_range(G, rng, propertymap): + def func(maximum, minimum, new_maximum, new_minimum, value): + return ((value-minimum)/(maximum-minimum)*(new_maximum-new_minimum)+new_minimum) + + mx = max(G.vertex_properties[propertymap].get_array()) + mn = min(G.vertex_properties[propertymap].get_array()) + G.vertex_properties["map"] = G.new_vertex_property("float") + gt.map_property_values(G.vertex_properties[propertymap], G.vertex_properties["map"], lambda x: func(mx, mn, rng[0], rng[1], x)) + + return G.vertex_properties["map"] + + ''' + G and propertymap in G are passed in. + Maps select property from range to color, works for int's floats and vec3 + and vec4 values + Returns a vector property map that maps the property value to color. + If passed a vertex propertymap, returns a vertex property map + If passed an edge propertymap, returns an edge property map + ''' + def map_property_to_color(G, propertymap, colormap = 'tab20'): + key_type = '' + if propertymap.key_type()=='e': + key_type = 'e' + elif propertymap.key_type()=='v': + key_type = 'v' + else: + #TO DO: Write graph property type return. + print("Graph_property passed") + if G.properties[(key_type, 'RGBA')] == None: + color = [0.0, 0.0, 0.0, 1.0] + G.properties[(key_type, 'RGBA')] = G.new_property(key_type, "vector", val=color) + + #This is when the property map is integers + #int32_t when integer + value_type = propertymap.value_type() + print(value_type) + if(value_type == "int32_t"): + array = propertymap.get_array() + colors = cm.get_cmap(colormap, len(np.unique(array))) + if key_type == 'v': + for v in G.vertices(): + G.vertex_properties["RGBA"][v] = colors(propertymap[v]) + return G.vertex_properties["RGBA"] + elif(value_type == 'double' or value_type == 'float'): + array = propertymap.get_array() + norm = cm.colors.Normalize(vmin=min(array), vmax=max(array)) + #colors = cm.get_cmap(colormap, len(np.unique(array))) + colors = cm.ScalarMappable(norm=norm, cmap=colormap) + if key_type =='v': + for v in G.vertices(): + print(colors.to_rgba(propertymap[v])) + G.vertex_properties["RGBA"][v] = colors.to_rgba(propertymap[v]) + return G.vertex_properties["RGBA"] + elif(value_type == 'vector' or value_type == 'vector'): + print("detected a vector of doubles or floats") + + + + ''' + Attempts to partition the graph G into b_i = [0, n^3] sections, where each section is the connected componenets of + m steps away from the source node (where the source node is the closest node to point n_i) + ''' + def partition(self, G, n, m, is_dual, path, prefix, saveKey, lbl = "none"): + bb = AABB(G, is_dual) + #print("FLJKKHDFLKJFDLKJFDLKJ ", m) + x, y, z = bb.project_grid(n) + #cluster_belongance = G.new_vertex_property("vector") + #ecluster_belongance = G.new_edge_property("vector") + #X, Y, Z = np.meshgrid(x,y,z) + array_cluster_bel = np.zeros([n**3, G.num_vertices()]) + array_ecluster_bel = np.zeros([n**3, G.num_edges()]) + c = 0 + for i in range(len(x)): + for j in range(len(y)): + for k in range(len(z)): + #for i in range(1): + # for j in range(1): + # for k in range(1): + p = G.vertex_properties["p"].get_2d_array(range(G.num_vertices())) + if(p.shape[1] < 3): + break + #print("stuff", p) + idx = -1 + dist = 100000000000 + for M in range(p.shape[1]): + d = math.sqrt((x[i] - p[0][M])**2 + (y[j] - p[1][M])**2 + (z[k] - p[2][M])**2) + if d < dist: + idx = M + dist = d + print(idx, " vertex[",p[0][idx], p[1][idx], p[2][idx], "], point [", x[i], y[j], z[k], "]") + clusters = G.new_vertex_property("int", vals=np.full((G.num_vertices(), 1), -1, dtype=int)) + #eclusters = G.new_edge_property("int", vals=np.full((G.num_edges(), 1), -1, dtype=int)) + dist = G.new_vertex_property("int", val = 100000) + dist[G.vertex(idx)] = 0 + #run a bfs on the index node + gt.bfs_search(G, idx, VisitorClassPartition(clusters, dist, c, m)) + #print(dist.get_array()); + temp_vertices = np.ma.masked_where(dist.get_array() <= m, dist.get_array()).mask + #print(temp_vertices) + temp_edges = np.zeros([G.num_edges()]) + for vertex in np.nditer(np.where(temp_vertices == True)): + for edge in G.get_out_edges(vertex): + #print(edge) + temp_edges[edge[2]] = 1 + array_cluster_bel[c,:] = temp_vertices[:] + array_ecluster_bel[c,:] = temp_edges[:] + c = c + 1 + #print(G.num_vertices()) + #np.savetxt("clusters.txt", array_cluster_bel) + G.vertex_properties["partition"] = G.new_vertex_property("vector", array_cluster_bel) + G.edge_properties["partition"] = G.new_edge_property("vector", array_ecluster_bel) + j = 0 + gt.graph_draw(G, pos = gt.sfdp_layout(G), output="Graph_full.pdf") + for i in range(c): + TFv = G.new_vertex_property("bool", vals=array_cluster_bel[i,:]) + TFe = G.new_edge_property("bool", vals=array_ecluster_bel[i,:]) + G.set_filters(TFe, TFv) + G1 = gt.Graph(G, directed=False, prune=True) + G.clear_filters() + G1.clear_filters() + iteration = 0 + title = "graph_" + str(iteration) + ".pdf" + gt.graph_draw(G1, pos = gt.sfdp_layout(G1), output=title) + while(np.any(G1.degree_property_map("total").get_array() == True)): + #print(np.any(G1.degree_property_map("total").get_array()), " ", i) + iteration = iteration + 1 + G1 = self.filterBorder(G1, is_dual) + title = "graph_" + str(iteration) + ".pdf" + gt.graph_draw(G1, pos = gt.sfdp_layout(G1), output=title) + #G1 = self.filterDisconnected(G1, is_dual) + G1 = self.recalculate_metrics(G1, is_dual) + + if(G1.num_edges() > 2): + ppath = path + "/partial"+str(j)+".nwt" + self.saveGraph_gt(G1, ppath) + if(i == 0): + self.saveSample_gt(G1, path, prefix, is_dual, writeKey=saveKey, label=lbl) + else: + self.saveSample_gt(G1, path, prefix, is_dual, label=lbl) + j = j + 1 + + + #print(array_cluster_bel) + + + def find_loops(self, G, path, is_dual, largest): + iteration = 0 + if(largest): + G = self.filterDisconnected(G) + gt.graph_draw(G, pos = gt.sfdp_layout(G), output="./graph_full.pdf") + while(np.any(G.degree_property_map("total").get_array() == True)): + #print(np.any(G1.degree_property_map("total").get_array()), " ", i) + iteration = iteration + 1 + G = self.filterBorder(G, is_dual) + title = path +"graph_" + str(iteration) + ".pdf" + gt.graph_draw(G, pos = gt.sfdp_layout(G), output=title) + #G1 = self.filterDisconnected(G1, is_dual) + G = self.recalculate_metrics(G, is_dual) + ppath = path + "graph_full_mst.pdf" + gt.graph_draw(G, pos = gt.sfdp_layout(G, C=0.4), output=ppath, edge_color=G.edge_properties["mst"], vertex_text=G.vertex_index,output_size=(1000,1000), vertex_font_size = 6, vertex_size = 3) + paths = G.new_edge_property("int", vals=np.full((G.num_edges(), 1), 1, dtype=int)) + #values = np.full((G.num_vertices(),1), 0, dtype=int) + #for i in range(G.num_vertices()): + # values[i,0] = G.vertex_index[G.vertex(i)] + #G.vertex_properties["original_index"] = values + + #add index, then use index to search for new vertecies. + G.edge_properties["path_weights"] = paths + cycle_list = [] + for e in G.edges(): + if G.edge_properties["mst"][e] == 0: + cycle_list.append(e) + # else: + # G.edge_properties["path_weights"][e] = 1 + + + #for e in cycle_list: + #print("Loop Detected between %i, %i", e.source(), e.target()) + + G.set_edge_filter(G.edge_properties["mst"]) + iteration = 0 + G1 = gt.Graph(G, directed=False, prune=True) + G.clear_filters() + for e in cycle_list: + edge = G1.add_edge(e.target(), e.source()) + #G.edge_properties["mst"][e] = 1 + G1.edge_properties["path_weights"][edge] = 1000 + #G.set_edge_filter(G.edge_properties["mst"]) + #G1 = gt.Graph(G, directed=False, prune=True) + #s = -1 + #t = -1 + #for edge in G1.edges(): + # if(G1.vertex_properties["original_index"][edge.source()] == int(e.source())): + # s = int(edge.source()) + # elif(G1.vertex_properties["original_index"][edge.target()] == int(e.target())): + # t = int(edge.target()) + + [Vl, El] = gt.shortest_path(G1, e.source(), e.target(), weights=G1.edge_properties["path_weights"]) + clusters = G1.new_vertex_property("int", vals=np.full((G1.num_vertices(), 1), 5, dtype=int)) + eclusters = G1.new_edge_property("int", vals = np.full((G1.num_edges(), 1), 3, dtype=int)) + G1.vertex_properties["cycle"] = clusters + G1.edge_properties["cycle"] = eclusters + print("Number of vertices in Path is:", len(Vl)) + for v in Vl: + G1.vertex_properties["cycle"][G1.vertex(v)] = 10 + print(str(v)) + #Create the arrays to be histogrammed later regarding every loop + length_total = 0 + volume_total = 0 + tortuosity_total = 0 + bc_total = 0 + bcl_total = 0 + num_edges = 1 + for v in range(len(Vl)-1): + G1.edge_properties["cycle"][G1.edge(Vl[v], Vl[v+1])] = 10 + length_total = length_total + G1.edge_properties["length"][G1.edge(Vl[v], Vl[v+1])] + volume_total = volume_total + G1.edge_properties["volume"][G1.edge(Vl[v], Vl[v+1])] + tortuosity_total = tortuosity_total + G1.edge_properties["tortuosity"][G1.edge(Vl[v], Vl[v+1])] + bc_total = bc_total + G1.edge_properties["bc"][G1.edge(Vl[v], Vl[v+1])] + bcl_total = bcl_total + G1.edge_properties["bc*length"][G1.edge(Vl[v], Vl[v+1])] + num_edges = num_edges + 1 + G1.edge_properties["cycle"][G1.edge(e.source(), e.target())] = 10 + length_total = length_total + G.edge_properties["length"][G.edge(e.source(), e.target())] + volume_total = volume_total + G.edge_properties["volume"][G.edge(e.source(), e.target())] + tortuosity_total = tortuosity_total + G.edge_properties["tortuosity"][G.edge(e.source(), e.target())] + bc_total = bc_total + G.edge_properties["bc"][G.edge(e.source(), e.target())] + bcl_total = bcl_total + G.edge_properties["bc*length"][G.edge(e.source(), e.target())] + ppath = path + "Loop_Histograms.txt" + f = open(ppath, "a+") + f.write("%.15f\t" % length_total) + f.write("%.15f\t" % volume_total) + f.write("%.15f\t" % tortuosity_total) + f.write("%.15f\t" % bc_total) + f.write("%.15f\t" % bcl_total) + f.write("%d\t\n" % num_edges) + f.close() + + title = path+"cycle" + str(iteration)+".png" + gt.graph_draw(G1, pos=G1.vertex_properties["pos"], output=title, edge_color='red', edge_pen_width=G1.edge_properties["cycle"], vertex_text=G1.vertex_index, output_size=(1920,1280), vertex_font_size = 4, vertex_fill_color = G1.vertex_properties["cycle"], vertex_size = 3) + G1.remove_edge(edge) + #G.clear_filters(); + #G.edge_properties["mst"][e] = 0 + #G.edge_properties["path_weights"][e] = 1000 + iteration = iteration + 1 + ''' + returns an affinity matrix based on the property edge property, given as a string. + ''' + def get_affinity_matrix(G, edge_property, gaussian=False, sigma = None): + affinity = gt.shortest_distance(G, weights=G.edge_properties[edge_property]) + affinity = affinity.get_2d_array(range(G.num_vertices())) + if gaussian: + if sigma == None: + sig = 60.0 + else: + sig = sigma + for i in range(G.num_vertices()): + for j in range(G.num_vertices()): + affinity[i,j] = np.exp(-np.power(affinity[i, j], 2.) / (2 * np.power(sig, 2.))) + return affinity + + + ''' + Attempts to apporixmate the number of clusters in the graph by finding + the approximate number of minimum graph cuts + + TO_DO: Program the gaussian and sigma parameters + ''' + def approximate_n_cuts(G, affinity_property): + #L = gt.laplacian(G, weight = G.edge_properties[affinity_property]) + G1 = gt.Graph(G, directed=False) + G1 = Network.filterDisconnected(G1, G1) + L = Network.get_affinity_matrix(G1, affinity_property, gaussian=True, sigma=157) + L = sp.sparse.csgraph.laplacian(L, normed=True) + n_components = G1.num_vertices() + eigenvalues, eigenvectors = eigsh(L, k=n_components, which = "LM", sigma=1.0, maxiter = 5000) + + plt.figure() + plt.scatter(np.arange(len(eigenvalues)), eigenvalues) + plt.grid() + plt.show() + + count = sum(eigenvalues > 1.01) + return count + ''' + Cluster the vectices in the graph based on a property passed to the clustering algorithm + + TO_DO: program the guassian and sigma parameters + ''' + def spectral_clustering(G, affinity_property, gaussian = False, sigma = None, n_clusters = None, num_iters = 1000): + if(n_clusters == None): + n_c = Network.approximate_n_cuts(G, affinity_property); + else: + n_c = n_clusters + sc = SpectralClustering(n_c, affinity='precomputed', n_init = num_iters) + sc.fit(Network.get_affinity_matrix(G, affinity_property, gaussian = True, sigma=157)) + + return sc.labels_ + + + def export_points(self, G): + location = "./vertex_points.txt" + f = open(location, "a+") + for v in G.vertices(): + f.write("%.15f\t" % G.vertex_properties["p"][v][0]) + f.write("%.15f\t" % G.vertex_properties["p"][v][1]) + f.write("%.15f\n" % G.vertex_properties["p"][v][2]) + + f.close() + + def export_bb(self, G): + location = "./bb_points.txt" + aabb = AABB(G) + points = aabb.vertices() + f = open(location, "a+") + for i in points: + f.write("%.15f\t" % i[0]) + f.write("%.15f\t" % i[1]) + f.write("%.15f\n" % i[2]) + + f.close() + + + def write_vtk(G, location, camera = None, binning = True): + #location = "./nwt.vtk" + f = open(location, "w+") + f.write("# vtk DataFile Version 1.0\n") + f.write("Data values\n") + f.write("ASCII\n\n") + f.write("DATASET POLYDATA\n") + num_pts = 0 + for e in G.edges(): + X = G.edge_properties["x"][e] + num_pts += len(X) + f.write("POINTS %d float\n" % num_pts) + for e in G.edges(): + X = G.edge_properties["x"][e] + Y = G.edge_properties["y"][e] + Z = G.edge_properties["z"][e] + #pts = np.array([X,Y,Z]).T + for p in range(len(X)): + f.write("%.15f\t" % X[p]) + f.write("%.15f\t" % Y[p]) + f.write("%.15f\n" % Z[p]) + f.write("LINES " + str(G.num_edges()) + " " + str(G.num_edges()+num_pts) + "\n") + num_pts = 0 + for e in G.edges(): + X = G.edge_properties["x"][e] + indices = list(range(0, len(X))) + indices = [x+num_pts for x in indices] + f.write(str(len(X)) + " ") + f.write(" ".join(str(x) for x in indices)) + f.write("\n") + num_pts += len(X) + f.write("POINT_DATA %d\n" % num_pts) + f.write("SCALARS radius float 1\n") + f.write("LOOKUP_TABLE radius_table\n") + for e in G.edges(): + R = G.edge_properties["r"][e] + #pts = np.array([X,Y,Z]).T + for p in range(len(R)): + f.write("%.15f\n" % R[p]) + + f.write("COLOR_SCALARS color_table %d\n" % 4) + bins = [0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0] + + if camera == None: + for e in G.edges(): + X = G.edge_properties["x"][e] + #RGBA = G.edge_properties["RGBA"].get_2d_array(range(4)).T + for p in range(len(X)): + f.write("%.15f\t" % G.edge_properties["RGBA"][e][0]) + f.write("%.15f\t" % G.edge_properties["RGBA"][e][1]) + f.write("%.15f\t" % G.edge_properties["RGBA"][e][2]) + if binning: + index = np.digitize(G.edge_properties["RGBA"][e][3], bins, right=True) + if (index >= len(bins) or index < 0): + print(G.edge_properties["RGBA"][e][3]) + print(index) + f.write("%.15f\n" % bins[index]) + else: + f.write("%.15f\n" % G.edge_properties["RGBA"][e][3]) + else: + ptx = G.graph_properties["point_cloud_x"].get_array() + pty = G.graph_properties["point_cloud_y"].get_array() + ptz = G.graph_properties["point_cloud_z"].get_array() + pts = np.array([ptx,pty,ptz]).T + temp = AABB(G) + bbl = temp.A + bbu = temp.B + tp = np.zeros((1,3), dtype = np.float64) + tp[0][0] = (bbu[0]-bbl[0])/2 + tp[0][1] = (bbu[1]-bbl[1])/2 + tp[0][2] = (bbu[2]-bbl[2])/2 + pts[:] = pts[:] - camera - \ + tp + distance = np.zeros(pts.shape[0], dtype = np.float32) + for i in range(distance.shape[0]): + distance[i] = np.sqrt(np.power(pts[i][0],2) + np.power(pts[i][1],2) + np.power(pts[i][2],2)) + mx = max(distance) + mn = min(distance) + distance = (distance - mx)/(mn-mx) + idx = 0 + for e in G.edges(): + X = G.edge_properties["x"][e] + for p in range(len(X)): + f.write("%.15f\t" % G.edge_properties["RGBA"][e][0]) + f.write("%.15f\t" % G.edge_properties["RGBA"][e][1]) + f.write("%.15f\t" % G.edge_properties["RGBA"][e][2]) + if binning: + index = np.digitize(distance[idx], bins, right=True) + f.write("%.15f\n" % bins[index]) + else: + f.write("%.15f\n" % distance[idx]) + idx += 1 + + f.close() + + + ''' + Creates a graph from a list of nodes and a list of edges. + Uses edge turtuosity as weight. + Returns a NetworkX Object. + ''' + def createTortuosityGraph(self): + G = nx.Graph() + for i in range(len(self.N)): + G.add_node(i, p=self.N[i].p) + for i in range(len(self.F)): + G.add_edge(self.F[i].v0, self.F[i].v1, weight = self.F[i].turtuosity()) + G[self.F[i].v0][self.F[i].v1]['pts'] = self.F[i].points + G[self.F[i].v0][self.F[i].v1]['rads'] = self.F[i].radii + + return G + + + ''' + Creates a graph from a list of nodes and a list of edges. + Uses edge volume as weight. + Returns a NetworkX Object. + ''' + def createVolumeGraph(self): + G = nx.Graph() + for i in range(len(self.N)): + G.add_node(i, p=self.N[i].p) + for i in range(len(self.F)): + G.add_edge(self.F[i].v0, self.F[i].v1, weight = self.F[i].volume()) + G[self.F[i].v0][self.F[i].v1]['pts'] = self.F[i].points + G[self.F[i].v0][self.F[i].v1]['rads'] = self.F[i].radii + + return G +#''' +#Returns the positions dictionary for the Circular layout. +#''' +#def getCircularLayout(graph, dim, radius): +# return nx.circular_layout(graph, dim, radius) +# +#''' +#Return the positions dictionary for the Spring layout. +#''' +#def getSpringLayout(graph, pos, iterations, scale): +# return nx.spring_layout(graph, 2, None, pos, iterations, 'weight', scale, None) +# +#''' +#Draws the graph. +#''' +#def drawGraph(graph, pos): +# nx.draw(graph, pos) +# return + + def aabb(self): + + lower = self.N[0].p.copy() + upper = lower.copy() + for i in self.N: + for c in range(len(lower)): + if lower[c] > i.p[c]: + lower[c] = i.p[c] + if upper[c] < i.p[c]: + upper[c] = i.p[c] + return lower, upper + + #calculate the distance field at a given resolution + # R (triple) resolution along each dimension + def distancefield(self, R=(100, 100, 100)): + + #get a list of all node positions in the network + P = [] + for e in self.F: + for p in e.points: + P.append(p) + + #turn that list into a Numpy array so that we can create a KD tree + P = np.array(P) + + #generate a KD-Tree out of the network point array + tree = sp.spatial.cKDTree(P) + + plt.scatter(P[:, 0], P[:, 1]) + + #specify the resolution of the ouput grid + R = (200, 200, 200) + + #generate a meshgrid of the appropriate size and resolution to surround the network + lower, upper = self.aabb(self.N, self.F) #get the space occupied by the network + x = np.linspace(lower[0], upper[0], R[0]) #get the grid points for uniform sampling of this space + y = np.linspace(lower[1], upper[1], R[1]) + z = np.linspace(lower[2], upper[2], R[2]) + X, Y, Z = np.meshgrid(x, y, z) + #Z = 150 * numpy.ones(X.shape) + + + Q = np.stack((X, Y, Z), 3) + + + D, I = tree.query(Q) + + return D + + #returns the number of points in the network + def npoints(self): + n = 0 #initialize the counter to zero + for f in self.F: #for each fiber + n = n + len(f.points) - 2 #count the number of points in the fiber - ignoring the end points + n = n + len(self.N) #add the number of nodes (shared points) to the node count + return n #return the number of nodes + + #returns all of the points in the network + def points(self): + k = self.npoints() + P = np.zeros((3, k)) #allocate space for the point list + + idx = 0 + for f in self.F: #for each fiber in the network + for ip in range(1, len(f.points)-1): #for each point in the network + P[:, idx] = f.points[ip] #store the point in the raw point list + idx = idx + 1 + return P #return the point array + + #returns the number of linear segments in the network + def nsegments(self): + n = 0 #initialize the segment counter to 0 + for f in self.F: #for each fiber + n = n + len(f.points) - 1 #calculate the number of line segments in the fiber (points - 1) + return n #return the number of line segments + + #return a list of line segments representing the network + def segments(self, dtype=np.float32): + k = self.nsegments() #get the number of line segments + start = np.zeros((k, 3),dtype=dtype) #start points for the line segments + end = np.zeros((k, 3), dtype=dtype) #end points for the line segments + + idx = 0 #initialize the index counter to zero + for f in self.F: #for each fiber in the network + for ip in range(0, len(f.points)-1): #for each point in the network + start[idx, :] = f.points[ip] #store the point in the raw point list + idx = idx + 1 + + idx = 0 + for f in self.F: #for each fiber in the network + for ip in range(1, len(f.points)): #for each point in the network + end[idx, :] = f.points[ip] #store the point in the raw point list + idx = idx + 1 + + return start, end + + #returns the points inside the fiber. + def fiber(self, idx): + p = self.F[idx].points + X = np.zeros(len(p)) + Y = np.zeros(len(p)) + Z = np.zeros(len(p)) + for i in range(len(p)): + X[i] = p[i][0] + Y[i] = p[i][1] + Z[i] = p[i][2] + return X, Y, Z + + #function returns the fiber associated with a given 1D line segment index + def segment2fiber(self, idx): + i = 0 + for f in range(len(self.F)): #for each fiber in the network + i = i + len(self.F[f].points)-1 #add the number of points in the fiber to i + if i > idx: #if we encounter idx in this fiber + return self.F[f].points, f #return the fiber associated with idx and the index into the fiber array + + def vectors(self, clock=False, dtype=np.float32): + if clock: + start_time = time.time() + start, end = self.segments(dtype) #retrieve all of the line segments + v = end - start #calculate the resulting vectors + l = np.sqrt(v[:, 0]**2 + v[:,1]**2 + v[:,2]**2) #calculate the fiber lengths + z = l==0 #look for any zero values + nz = z.sum() + if nz > 0: + print("WARNING: " + str(nz) + " line segment(s) of length zero were found in the network and will be removed" ) + + if clock: + print("Network::vectors: " + str(time.time() - start_time) + "s") + + return np.delete(v, np.where(z), 0) + + #scale all values in the network by tuple S = (sx, sy, sz) + def scale(self, S): + for f in self.F: + for p in f.points: + p[0] = p[0] * S[0] + p[1] = p[1] * S[1] + p[2] = p[2] * S[2] + + for n in self.N: + n.p[0] = n.p[0] * S[0] + n.p[1] = n.p[1] * S[1] + n.p[2] = n.p[2] * S[2] + + + #calculate the adjacency weighting function for the network given a set of vectors X = (x, y, z) and weight exponent k + def adjacencyweight(self, P, k=200, length_threshold = 25, dtype=np.float32): + V = self.vectors(dtype) #get the vectors representing each segment + #V = V[0:n_vectors, :] + L = np.expand_dims(np.sqrt((V**2).sum(1)), 1) #calculate the length of each vector + + outliers = L > length_threshold #remove outliers based on the length_threshold + V = np.delete(V, np.where(outliers), 0) + L = np.delete(L, np.where(outliers)) + V = V/L[:,None] #normalize the vectors + + P = np.stack(spharmonics.sph2cart(1, P[0], P[1]), P[0].ndim) + PV = P[...,None,:] * V + cos_alpha = PV.sum(PV.ndim-1) + W = np.abs(cos_alpha) ** k + + return W, L + + diff --git a/subgraph_shaders.py b/subgraph_shaders.py new file mode 100644 index 0000000..aac8838 --- /dev/null +++ b/subgraph_shaders.py @@ -0,0 +1,276 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Created on Mon Aug 5 15:35:04 2019 + +@author: pavel +""" + +#shaders for drawing supernodes +vert_s = """ +#version 120 +// Uniforms +// ------------------------------------ +uniform mat4 u_model; +uniform mat4 u_view; +uniform mat4 u_projection; +uniform vec3 u_graph_size; +uniform bool u_picking; +// Attributes +// ------------------------------------ +attribute vec3 a_position; +attribute vec4 a_bg_color; +attribute vec4 a_cluster_color; +attribute vec2 a_value; +attribute float a_arc_length; +attribute vec4 a_outer_arc_length; +attribute vec4 a_unique_id; + +// Varyings +// ------------------------------------ +varying vec4 v_bg_color; +varying vec4 v_cluster_color; +varying vec2 v_value; +varying float v_zoom_level; +varying float v_arc_length; +varying vec4 v_outer_arc_length; +varying vec4 v_unique_id; + +void main (void) { + v_value = a_value; + v_bg_color = a_bg_color; + v_cluster_color = a_cluster_color; + gl_Position = u_projection * u_view * u_model * + vec4(a_position, 1.0); + v_zoom_level = u_view[0][0]; + v_arc_length = a_arc_length; + v_outer_arc_length = a_outer_arc_length; + v_unique_id = a_unique_id; +} +""" + +frag_s = """ +#version 120 +const float pi = 3.1415926535897932384626433832795; +// Varyings +// ------------------------------------ +varying vec4 v_bg_color; +varying vec2 v_value; +varying float v_zoom_level; +varying vec4 v_cluster_color; +varying float v_arc_length; +varying vec4 v_outer_arc_length; +varying vec4 v_unique_id; + +uniform bool u_picking; + + +// Attributes +// ------------------------------------ + +// Functions +// ------------------------------------ +float new_alpha(float zoom_level); +float atan2(float y, float x); +// Main +// ------------------------------------ +void main() +{ + if(!u_picking) + { + //This is the color of the outer arc lengths + float R = 0.55; + float R2 = 0.3; + float dist = sqrt(dot(v_value, v_value)); + float sm = smoothstep(R, R-0.0010, dist); + float sm2 = smoothstep(R2, R2+0.0010, dist); + float alpha = sm*sm2; + + float n_alpha = 1.0; + if(v_zoom_level < 0.0010) + { + n_alpha = new_alpha(v_zoom_level); + } + else + { + discard; + } + float angle = atan(v_value[1], v_value[0]); + if(dist < 0.25) + { + //gl_FragColor = v_cluster_color; + gl_FragColor = vec4(v_cluster_color.rgb, n_alpha); + } + if(dist > 0.3 && 0.55 > dist) + { + float angle2 = angle+pi; + if(angle2 > v_arc_length) + { + discard; + } + else + { + gl_FragColor = vec4(v_bg_color.rgb, n_alpha*alpha); + } + } + angle = atan(v_value[1]/dist, v_value[0]/dist); + //Need to add antialiasing to all other circles in the form of smoothstep. + + //THIS NEED TO BE FIXED + if(dist > 0.61 && 0.9 > dist) + { + if(angle < v_outer_arc_length[1]) + { + gl_FragColor = vec4(0.32, 0.61, 0.93, n_alpha); + //gl_FragColor = vec4(1.0, 0.0, 0.0, n_alpha); + } + else if(angle > v_outer_arc_length[1] && angle < v_outer_arc_length[2]) + { + gl_FragColor = vec4(0.337, 0.866, 0.415, n_alpha); + //gl_FragColor = vec4(0.0, 1.0, 0.0, n_alpha); + } + else if(angle > v_outer_arc_length[2] && angle < v_outer_arc_length[3]) + { + gl_FragColor = vec4(0.988, 0.631, 0.058, n_alpha); + //gl_FragColor = vec4(0.0, 0.0, 1.0, n_alpha); + } + else //if(angle > v_outer_arc_length[3] && angle < pi) + { + gl_FragColor = vec4(0.93, 0.949, 0.329, n_alpha); + //gl_FragColor = vec4(1.0, 1.0, 0.0, n_alpha); + } + //else + //{ + // discard; + //} + +// if(angle > -pi && angle < -pi/4.0) +// { +// //gl_FragColor = vec4(0.32, 0.61, 0.93, n_alpha); +// gl_FragColor = vec4(1.0, 0.0, 0.0, n_alpha); +// } +// else if(angle > -pi/4.0 && angle < 0.0) +// { +// gl_FragColor = vec4(0.0, 1.0, 0.0, n_alpha); +// } +// else if(angle > 0.0 && angle < pi/2.0) +// { +// gl_FragColor = vec4(0.0, 0.0, 1.0, n_alpha); +// } +// else if(angle > pi/2.0 && angle < pi) +// { +// gl_FragColor = vec4(1.0, 1.0, 0.0, n_alpha); +// } + } + } + else + { + float dist = sqrt(dot(v_value, v_value)); + if (dist > 0.9) + discard; + else + gl_FragColor = v_unique_id; + } +} + +float atan2(float y, float x) +{ + float s; + if(abs(x) > abs(y)) + s = 1.0; + else + s = 0.0; + return (mix(pi/2.0 - atan(x,y), atan(y,x), s)); +} + +//float atan2(in float y, in float x) +//{ +// return x == 0.0 ? sign(y)*pi/2 : atan(y, x); +//} + +float new_alpha(float zoom_level) +{ + float val = 1 - (zoom_level - 0.00075)/(0.0010-0.00075); + if(val > 1.0) + { + val = 1.0; + } + return val; +} +""" + +vs_s = """ +#version 120 + +// Uniforms +// ------------------------------------ +uniform mat4 u_model; +uniform mat4 u_view; +uniform mat4 u_projection; + +// Attributes +// ------------------------------------ +attribute vec3 a_position; +attribute vec2 a_normal; +attribute vec4 a_fg_color; +attribute float a_linewidth; +//attribute vec4 a_unique_id; +//attribute vec4 l_color; + +// Varyings +// ------------------------------------ +varying vec4 v_fg_color; +varying float v_zoom_level; +varying vec2 v_normal; +varying float v_linewidth; + +void main() { + vec3 delta = vec3(a_normal*a_linewidth/((1-u_view[0][0])*0.1), 0); + //vec3 delta = vec3(a_normal*a_linewidth/(1-u_view[0][0]), 0); + gl_Position = u_model * u_view * u_projection * vec4(a_position+delta, 1.0); + //gl_Position = u_projection * u_view * u_model * + // vec4(a_position, 1.0); + v_zoom_level = u_view[0][0]; + v_fg_color = a_fg_color; + v_normal = a_normal; + v_linewidth = a_linewidth; + +} +""" + +fs_s = """ +#version 120 +// Varying +// ------------------------------------ +varying vec4 v_fg_color; +varying float v_zoom_level; +varying vec2 v_normal; +varying float v_linewidth; + +float new_alpha(float zoom_level); + +void main() +{ + float l = length(v_normal); + float feather = 0.5; + float alpha = 1.0; + if(l > v_linewidth/2.0+feather) + discard; + else + alpha = 0.0; + + if(v_zoom_level < 0.0010) + alpha = new_alpha(v_zoom_level); + gl_FragColor = vec4(v_fg_color.rgb, alpha); +} + +float new_alpha(float zoom_level) +{ + float val = 1 - (zoom_level - 0.00075)/(0.0010-0.00075); + if(val > 1.0) + { + val = 1.0; + } + return val; +} +""" \ No newline at end of file diff --git a/tube_shaders.py b/tube_shaders.py new file mode 100644 index 0000000..642eee3 --- /dev/null +++ b/tube_shaders.py @@ -0,0 +1,232 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Created on Mon Aug 5 15:58:30 2019 + +@author: pavel +""" + +#vertex shader for the tube drawing +VERT_SHADER = """ +// Uniforms +// ------------------------------------ +uniform vec4 u_bb[26]; +uniform vec3 u_LightPost; +uniform mat4 u_model; +//uniform mat4 u_view; +uniform mat4 u_projection; + +uniform vec3 u_eye; +uniform vec3 u_up; +uniform vec3 u_target; + +// Attributes +// ------------------------------------ +attribute vec3 a_position; +attribute vec3 a_normal; +attribute vec4 a_fg_color; + +// Varyings +// ------------------------------------ +varying vec3 v_normal; +varying vec4 v_fg_color; +varying vec3 v_position; +varying mat4 u_view; +varying vec3 vt_position; + + +// Functions +// ------------------------------------ + +mat4 make_view(vec3 eye, vec3 target, vec3 up); + +void main (void) { + // Calculate position + mat4 u_view = make_view(u_eye, u_target, u_up); + + mat4 MVP = u_projection * u_view * u_model; + mat4 MV = u_view * u_model; + v_normal = vec3(MV * vec4(a_normal, 0.0)); + v_position = vec3(MV*vec4(a_position, 1.0)); + vt_position = vec3(a_position); + v_fg_color = a_fg_color; + gl_Position = MVP * vec4(a_position, 1.0); +} + +mat4 make_view(vec3 eye, vec3 target, vec3 up) +{ + vec3 zaxis = normalize(eye - target); + vec3 xaxis = normalize(cross(up, zaxis)); + vec3 yaxis = cross(zaxis, xaxis); + + mat4 viewMatrix = mat4( + vec4( xaxis.x, yaxis.x, zaxis.x, 0 ), + vec4( xaxis.y, yaxis.y, zaxis.y, 0 ), + vec4( xaxis.z, yaxis.z, zaxis.z, 0 ), + vec4(-dot( xaxis, eye ), -dot( yaxis, eye ), -dot( zaxis, eye ), 1 ) + ); + +// mat4 viewMatrix = mat4( +// vec4( xaxis.x, xaxis.y, xaxis.z, -dot( xaxis, eye ) ), +// vec4( yaxis.x, yaxis.y, yaxis.z, -dot( yaxis, eye ) ), +// vec4( zaxis.x, zaxis.y, zaxis.z, -dot( zaxis, eye ) ), +// vec4( 0, 0, 0, 1) +// ); + + + return viewMatrix; +} +""" + +#Fragment shader for the tube drawing. +FRAG_SHADER = """ +// Uniforms +// ------------------------------------ +uniform vec3 u_bb[26]; +uniform vec3 u_LightPos; +uniform mat4 u_model; +//uniform mat4 u_view; +uniform mat4 u_projection; + +uniform vec3 u_eye; +uniform vec3 u_up; +uniform vec3 u_target; + +// Varyings +// ------------------------------------ +varying vec3 v_normal; +varying vec4 v_fg_color; +varying vec3 v_position; +varying mat4 u_view; +varying vec3 vt_position; + +float new_alpha(); +float set_view(); +float bin_alpha(float alpha); + +void main() +{ + //Ambient Lighting + float ambientStrength = 0.5; + vec4 ambient_color = ambientStrength*v_fg_color; + + + mat4 MV = u_projection * u_view * u_model; + //vec3 u_LightP = vec3(0., 0., 0.); + vec3 u_LightP = vec3(u_eye); + //mat4 inv = inverse(u_view); + + //vec3 u_LightP = vec3(u_view[0][3], u_view[1][3], u_view[2][3]); + + //Diffuse Lighting + //float distance = length(u_LightP - v_position); + vec3 norm = normalize(v_normal); + vec3 lightVector = normalize(u_LightP - v_position); + float diffuse = max(dot(norm, lightVector), 0.0); + + //diffuse = diffuse * (1.0 / (1.0 + (0.25 * distance * distance))); + + vec4 color = (ambient_color + diffuse)*v_fg_color; + float alpha = new_alpha(); + //if(alpha == 1.0) + // gl_FragColor = vec4(1.0, 0.0, 0.0, 1.0); + // gl_FragColor = vec4(color.rgb, alpha); + //else + gl_FragColor = vec4(color.rgb, alpha); + // gl_FragColor = vec4(1.0, 0.0, 0.0, 1.0); + //gl_FragColor = v_fg_color; +} + +float new_alpha() +{ + float alpha = 0.0; + + //find the min and the max + float mx = -100000.0; + float mn = 100000.0; + vec3 u_light = vec3(u_eye); + for(int i = 0; i < 26; i++) + { + //vec3 rot_bb = vec3(u_projection * u_view * u_model * vec4(u_bb[i], 1.0)); + vec3 rot_bb = u_bb[i]; + if(length(rot_bb - u_light) < mn) + { + mn = length(rot_bb - u_light); + } + if(length(rot_bb - u_light) > mx) + { + mx = length(rot_bb - u_light); + } + } + + float l = length(u_light - vt_position); + alpha = (l-mn)/(mx-mn); + return bin_alpha(alpha); +// if(alpha > 0.9) +// return 1.0; +// else if(alpha > 0.8) +// return 0.9; +// else if(alpha > 0.7) +// return 0.8; +// else if (alpha > 0.6) +// return 0.7; +// else if (alpha > 0.5) +// return 0.6; +// else if (alpha > 0.4) +// return 0.5; +// else if (alpha > 0.3) +// return 0.4; +// else if (alpha > 0.2) +// return 0.3; +// else if (alpha > 0.1) +// return 0.2; +// else if (alpha > 0.0) +// return 0.0; +// else +// return alpha; +} + +float bin_alpha(float alpha) +{ + + + if(alpha > 0.8) + return 0.0; + else if(alpha > 0.6) + return 0.1; + else if (alpha > 0.4) + return 0.4; + else if (alpha > 0.2) + return 0.8; + else if (alpha > 0.0) + return 1.0; + else + return 1.0; + +// if(alpha > 0.9) +// return 0.0; +// else if(alpha > 0.8) +// return 0.1; +// else if(alpha > 0.7) +// return 0.2; +// else if (alpha > 0.6) +// return 0.3; +// else if (alpha > 0.5) +// return 0.4; +// else if (alpha > 0.4) +// return 0.5; +// else if (alpha > 0.3) +// return 0.6; +// else if (alpha > 0.2) +// return 0.7; +// else if (alpha > 0.1) +// return 0.8; +// else if (alpha > 0.0) +// return 0.9; +// else +// return 1.0; +} + +""" + + -- libgit2 0.21.4