Commit 193cb4c6c21380f415e8a5ef6c16c4b42980e090

Authored by Pavel Govyadinov
1 parent 2282da38

need this to test

@@ -828,8 +828,8 @@ class GraphCanvas(scene.SceneCanvas): @@ -828,8 +828,8 @@ class GraphCanvas(scene.SceneCanvas):
828 #self.clusters['a_linewidth'] = 4.*self.pixel_scale 828 #self.clusters['a_linewidth'] = 4.*self.pixel_scale
829 829
830 G.vertex_properties["pos"] = nwt.gt.sfdp_layout(G, groups = G.vertex_properties["clusters"], pos = G.vertex_properties["pos"]) 830 G.vertex_properties["pos"] = nwt.gt.sfdp_layout(G, groups = G.vertex_properties["clusters"], pos = G.vertex_properties["pos"])
831 - temp = [];  
832 - temp_pos = []; 831 + temp = []
  832 + temp_pos = []
833 #Find the global total of the metric. 833 #Find the global total of the metric.
834 global_metric = np.sum(G.edge_properties[edge_metric].get_array(), 0) 834 global_metric = np.sum(G.edge_properties[edge_metric].get_array(), 0)
835 unique_color = self.gen_cluster_id(G) 835 unique_color = self.gen_cluster_id(G)
@@ -93,7 +93,7 @@ def load_nwt(filepath): @@ -93,7 +93,7 @@ def load_nwt(filepath):
93 G, bbl, bbu = load_nwt("/home/pavel/Documents/Python/GraphGuiQt/network_4.nwt") 93 G, bbl, bbu = load_nwt("/home/pavel/Documents/Python/GraphGuiQt/network_4.nwt")
94 #nwt.Network.write_vtk(G) 94 #nwt.Network.write_vtk(G)
95 95
96 -#G, bbl, bbu = load_nwt("/home/pavel/Documents/Python/GraphGuiQt/network_test_251_1kx3_short_NEW.nwt") 96 +#G, bbl, bbu = load_nwt("/home/pavel/Documents/Python/GraphGuiQt/net`````````work_test_251_1kx3_short_NEW.nwt")
97 #ret = nwt.Network.get_affinity_matrix(G, "length") 97 #ret = nwt.Network.get_affinity_matrix(G, "length")
98 node_image = QtGui.QImage("/home/pavel/Documents/Python/GraphGuiQt/node_tex.jpg") 98 node_image = QtGui.QImage("/home/pavel/Documents/Python/GraphGuiQt/node_tex.jpg")
99 node_tex = QtGui.QBitmap.fromImage(node_image) 99 node_tex = QtGui.QBitmap.fromImage(node_image)
@@ -53,6 +53,10 @@ class TubeDraw(scene.SceneCanvas): @@ -53,6 +53,10 @@ class TubeDraw(scene.SceneCanvas):
53 high=(4-1)).astype(np.uint32) 53 high=(4-1)).astype(np.uint32)
54 self.vbo = gloo.VertexBuffer(self.cylinder_data) 54 self.vbo = gloo.VertexBuffer(self.cylinder_data)
55 self.triangles = gloo.IndexBuffer(self.triangle_data) 55 self.triangles = gloo.IndexBuffer(self.triangle_data)
  56 +
  57 + #generate an index buffer for the caps.
  58 + self.cap_data = np.random.randint(size=(5, 3), low=0, high=(4-1)).astype(np.uint32)
  59 + self.caps = gloo.IndexBuffer(self.cap_data)
56 self.program.bind(self.vbo) 60 self.program.bind(self.vbo)
57 self.scale = [1,1,1] 61 self.scale = [1,1,1]
58 self.r1 = np.eye(4, dtype=np.float32) 62 self.r1 = np.eye(4, dtype=np.float32)
@@ -61,7 +65,7 @@ class TubeDraw(scene.SceneCanvas): @@ -61,7 +65,7 @@ class TubeDraw(scene.SceneCanvas):
61 #set_state(clear_color='white', depth_test=True, blend=True, 65 #set_state(clear_color='white', depth_test=True, blend=True,
62 # blend_func=('src_alpha', 'one_minus_src_alpha'), depth_func = ('less'), cull_face='back') 66 # blend_func=('src_alpha', 'one_minus_src_alpha'), depth_func = ('less'), cull_face='back')
63 set_state(clear_color='white', depth_test=True, blend=True, 67 set_state(clear_color='white', depth_test=True, blend=True,
64 - blend_func=('src_alpha', 'one_minus_src_alpha'), depth_func = ('lequal'), cull_face='back') 68 + blend_func=('src_alpha', 'one_minus_src_alpha'), depth_func = ('lequal'))
65 #set_blend_color(color='black') 69 #set_blend_color(color='black')
66 #set_state('translucent') 70 #set_state('translucent')
67 self.program['u_LightPos'] = [0., 0., -1000.] 71 self.program['u_LightPos'] = [0., 0., -1000.]
@@ -93,7 +97,9 @@ class TubeDraw(scene.SceneCanvas): @@ -93,7 +97,9 @@ class TubeDraw(scene.SceneCanvas):
93 #create program 97 #create program
94 self.gen_cylinder_vbo(self.G, self.num_sides) 98 self.gen_cylinder_vbo(self.G, self.num_sides)
95 self.vbo = gloo.VertexBuffer(self.cylinder_data) 99 self.vbo = gloo.VertexBuffer(self.cylinder_data)
  100 + #self.triangle_data = np.append(self.triangle_data, self.cap_data, axis=0)
96 self.triangles = gloo.IndexBuffer(self.triangle_data) 101 self.triangles = gloo.IndexBuffer(self.triangle_data)
  102 + self.caps = gloo.IndexBuffer(self.cap_data)
97 103
98 #self.view = np.eye(4, dtype=np.float32) 104 #self.view = np.eye(4, dtype=np.float32)
99 self.model = np.eye(4, dtype=np.float32) 105 self.model = np.eye(4, dtype=np.float32)
@@ -153,6 +159,14 @@ class TubeDraw(scene.SceneCanvas): @@ -153,6 +159,14 @@ class TubeDraw(scene.SceneCanvas):
153 set_state(clear_color='white', depth_test=True, blend=True, 159 set_state(clear_color='white', depth_test=True, blend=True,
154 blend_func=('src_alpha', 'one_minus_src_alpha'), depth_func = ('lequal'), cull_face='back') 160 blend_func=('src_alpha', 'one_minus_src_alpha'), depth_func = ('lequal'), cull_face='back')
155 self.program.draw('triangles', self.triangles) 161 self.program.draw('triangles', self.triangles)
  162 + if self.select:
  163 + set_state(clear_color='white', depth_test=True, blend=True,
  164 + blend_func=('src_alpha', 'one_minus_src_alpha'), depth_func = ('always'), cull_face=False)
  165 + else:
  166 + set_state(clear_color='white', depth_test=True, blend=True,
  167 + blend_func=('src_alpha', 'one_minus_src_alpha'), depth_func = ('lequal'), cull_face=False)
  168 + self.program.draw('triangles', self.caps)
  169 +
156 self.projection = perspective(90.0, self.physical_size[0]/self.physical_size[1], 1.0, 1000.0) 170 self.projection = perspective(90.0, self.physical_size[0]/self.physical_size[1], 1.0, 1000.0)
157 self.program['u_projection'] = self.projection 171 self.program['u_projection'] = self.projection
158 172
@@ -161,9 +175,12 @@ class TubeDraw(scene.SceneCanvas): @@ -161,9 +175,12 @@ class TubeDraw(scene.SceneCanvas):
161 i = 0 175 i = 0
162 num_pts = 0 176 num_pts = 0
163 num_tri = 0 177 num_tri = 0
  178 + num_tri_caps = 0
164 for e in G.edges(): 179 for e in G.edges():
165 num_pts += len(self.G.edge_properties["x"][e]) 180 num_pts += len(self.G.edge_properties["x"][e])
166 num_tri += (len(self.G.edge_properties["x"][e])-1)*num_sides*2 181 num_tri += (len(self.G.edge_properties["x"][e])-1)*num_sides*2
  182 + num_tri_caps += (2*(num_sides-2))
  183 +
167 self.cylinder_data = np.zeros(num_pts*num_sides, dtype=[('a_position', np.float32, 3), 184 self.cylinder_data = np.zeros(num_pts*num_sides, dtype=[('a_position', np.float32, 3),
168 ('a_normal', np.float32, 3), 185 ('a_normal', np.float32, 3),
169 ('a_fg_color', np.float32, 4), 186 ('a_fg_color', np.float32, 4),
@@ -171,8 +188,12 @@ class TubeDraw(scene.SceneCanvas): @@ -171,8 +188,12 @@ class TubeDraw(scene.SceneCanvas):
171 ]) 188 ])
172 self.triangle_data = np.random.randint(size=(num_tri, 3), low=0, 189 self.triangle_data = np.random.randint(size=(num_tri, 3), low=0,
173 high=(G.num_edges()-1)).astype(np.uint32) 190 high=(G.num_edges()-1)).astype(np.uint32)
  191 +
  192 + self.cap_data = np.random.randint(size = (num_tri_caps, 3),
  193 + low = 0, high=(G.num_edges()-2)).astype(np.uint32)
174 index = 0 194 index = 0
175 t_index = 0 195 t_index = 0
  196 + c_index = 0
176 #for each edge generate a cylinder. 197 #for each edge generate a cylinder.
177 for e in G.edges(): 198 for e in G.edges():
178 #print("done") 199 #print("done")
@@ -266,9 +287,26 @@ class TubeDraw(scene.SceneCanvas): @@ -266,9 +287,26 @@ class TubeDraw(scene.SceneCanvas):
266 #generate the triangles 287 #generate the triangles
267 triangles = np.random.randint(size=((pts.shape[0]-1)*2*(num_sides), 3), low=0, 288 triangles = np.random.randint(size=((pts.shape[0]-1)*2*(num_sides), 3), low=0,
268 high=(G.num_edges()-1)).astype(np.uint32) 289 high=(G.num_edges()-1)).astype(np.uint32)
269 - 290 + cap_triangles = np.random.randint(size=((num_sides-2)*2 , 3), low=0, high=(num_sides-2)).astype(np.uint32)
270 t_index_temp = 0 291 t_index_temp = 0
  292 + c_index_temp = 0
  293 + #for every ring in the fiber
  294 + for ring in range(0, pts.shape[0], pts.shape[0]-1):
  295 + #if we have the first ring or the last ring add a cap.
  296 + idx_ori = index+ring*num_sides
  297 + for side in range(0, num_sides-2):
  298 + idx_current_point = index+ring*num_sides + side
  299 + if(side < num_sides-2):
  300 + cap_tri = [idx_ori, idx_current_point+1, idx_current_point+2]
  301 + else:
  302 + cap_tri = [idx_ori, idx_current_point+1, idx_ori]
  303 + cap_triangles[c_index_temp] = cap_tri
  304 + self.cap_data[c_index] = cap_tri
  305 + c_index_temp += 1
  306 + c_index += 1
  307 +
271 for ring in range(pts.shape[0]-1): 308 for ring in range(pts.shape[0]-1):
  309 + #otherwise generate the sides.
272 for side in range(num_sides): 310 for side in range(num_sides):
273 if(side < num_sides-1): 311 if(side < num_sides-1):
274 idx_current_point = index+ring*num_sides + side 312 idx_current_point = index+ring*num_sides + side
@@ -322,6 +360,13 @@ class TubeDraw(scene.SceneCanvas): @@ -322,6 +360,13 @@ class TubeDraw(scene.SceneCanvas):
322 ax.plot(pts[:,0], pts[:,1], pts[:,2]) 360 ax.plot(pts[:,0], pts[:,1], pts[:,2])
323 for j in range(pts.shape[0]): 361 for j in range(pts.shape[0]):
324 ax.plot(circle_pts[j,:,0], circle_pts[j,:,1], circle_pts[j,:,2]) 362 ax.plot(circle_pts[j,:,0], circle_pts[j,:,1], circle_pts[j,:,2])
  363 + for j in range(cap_triangles.shape[0]):
  364 + tri = np.zeros((3,4))
  365 + tri[:,0] = self.cylinder_data['a_position'][cap_triangles[j][0]]
  366 + tri[:,1] = self.cylinder_data['a_position'][cap_triangles[j][1]]
  367 + tri[:,2] = self.cylinder_data['a_position'][cap_triangles[j][2]]
  368 + tri[:,3] = self.cylinder_data['a_position'][cap_triangles[j][0]]
  369 + ax.plot(tri[0,:], tri[1,:], tri[2,:], c='b')
325 for j in range(triangles.shape[0]): 370 for j in range(triangles.shape[0]):
326 tri = np.zeros((3,4)) 371 tri = np.zeros((3,4))
327 tri[:,0] = self.cylinder_data['a_position'][triangles[j][0]] 372 tri[:,0] = self.cylinder_data['a_position'][triangles[j][0]]
@@ -13,6 +13,8 @@ Created on Mon Aug 5 15:53:16 2019 @@ -13,6 +13,8 @@ Created on Mon Aug 5 15:53:16 2019
13 13
14 from pyqtgraph.Qt import QtCore, QtGui, QtWidgets 14 from pyqtgraph.Qt import QtCore, QtGui, QtWidgets
15 from TubeCanvas import TubeDraw 15 from TubeCanvas import TubeDraw
  16 +import numpy as np
  17 +import trimesh as tm
16 18
17 DEBUG = False 19 DEBUG = False
18 20
@@ -50,7 +52,13 @@ class TubeWidget(QtGui.QWidget): @@ -50,7 +52,13 @@ class TubeWidget(QtGui.QWidget):
50 #Handles the mouse press events 52 #Handles the mouse press events
51 def on_mouse_press(self, event): 53 def on_mouse_press(self, event):
52 self.down = True 54 self.down = True
53 - n = 3 55 + if event.button == 2:
  56 + menu = QtWidgets.QMenu(self)
  57 + NS = menu.addAction('Export Model')
  58 + action = menu.exec_(event.native.globalPos())
  59 + if action == NS:
  60 + self.save_stl_mesh()
  61 +
54 62
55 #Handles the mouse wheel event 63 #Handles the mouse wheel event
56 def on_mouse_wheel(self, event): 64 def on_mouse_wheel(self, event):
@@ -86,3 +94,21 @@ class TubeWidget(QtGui.QWidget): @@ -86,3 +94,21 @@ class TubeWidget(QtGui.QWidget):
86 """ 94 """
87 def connect(self, signal_object): 95 def connect(self, signal_object):
88 signal_object.select.connect(self.select) 96 signal_object.select.connect(self.select)
  97 +
  98 + def save_stl_mesh(self):
  99 + triangles = self.canvas.triangle_data
  100 + caps = self.canvas.cap_data
  101 + points = self.canvas.cylinder_data['a_position']
  102 + colors = self.canvas.cylinder_data['a_fg_color']
  103 + normals = self.canvas.cylinder_data['a_normal']
  104 +
  105 + mesh = tm.Trimesh(vertices=points, faces=np.append(triangles, caps, axis=0),
  106 + vertex_colors=colors, vertex_normals=normals)
  107 + mesh.export('with_caps.obj')
  108 + mesh.export('with_caps.ply')
  109 +
  110 +
  111 +
  112 +
  113 +
  114 +
@@ -94,7 +94,7 @@ void main() @@ -94,7 +94,7 @@ void main()
94 // The marker function needs to be linked with this shader 94 // The marker function needs to be linked with this shader
95 float r = marker(gl_PointCoord, size); 95 float r = marker(gl_PointCoord, size);
96 float d = abs(r) - t; 96 float d = abs(r) - t;
97 - if( r > (v_linewidth/2.0+v_antialias)) 97 + if( r > (v_linewidth/3.0+v_antialias))
98 { 98 {
99 discard; 99 discard;
100 } 100 }
@@ -113,7 +113,7 @@ void main() @@ -113,7 +113,7 @@ void main()
113 else if(v_selection == 2.0) 113 else if(v_selection == 2.0)
114 gl_FragColor = vec4(0.0, 1.0, 0.0, v_fg_color.a); 114 gl_FragColor = vec4(0.0, 1.0, 0.0, v_fg_color.a);
115 else 115 else
116 - gl_FragColor = v_fg_color; 116 + gl_FragColor = vec4(v_fg_color.rgb, v_bg_color.a);
117 } 117 }
118 } 118 }
119 else 119 else
voronoi_test.py 0 → 100644
  1 +#!/usr/bin/env python3
  2 +# -*- coding: utf-8 -*-
  3 +"""
  4 +Created on Tue Sep 3 13:15:54 2019
  5 +
  6 +@author: pavel
  7 +"""
  8 +
  9 +from scipy.spatial import Voronoi, voronoi_plot_2d
  10 +from scipy.interpolate import interp1d
  11 +#import matplotlib._cntr as cntr
  12 +from shapely.geometry import Point
  13 +from shapely.geometry import Polygon
  14 +import numpy as np
  15 +import scipy as sp
  16 +import math
  17 +import matplotlib.pyplot as plt
  18 +import sys
  19 +import copy
  20 +
  21 +from skimage import measure
  22 +
  23 +from collections import defaultdict
  24 +
  25 +import network_dep as nwt
  26 +
  27 +class Polygon_mass:
  28 + def __init__(self, G):
  29 + self.G = G
  30 + self.get_aabb()
  31 + self.gen_polygon()
  32 + self.torque = []
  33 +
  34 + def add_torque(self, p, f):
  35 + d = self.CoM - p
  36 + r = np.linalg.norm(self.CoM - p)
  37 + theta = math.acos(np.dot(d, f)/np.dot(d, d)/np.dot(f, f))
  38 + torque = r * math.sin(theta) * f
  39 + self.torque.append(torque)
  40 +
  41 + def rotate(self, phi, direction = "counterclock"):
  42 + if("counterclock"):
  43 + for v in self.G.vertices:
  44 + p = self.G.vertex_properties["pos"][v]
  45 + p[0] = self.CoM[0] + math.cos(phi) * (self.CoM[0] - p[0]) - math.sin(phi) * (p[1] - self.CoM[1])
  46 + p[1] = self.CoM[1] + math.sin(phi) * (self.CoM[0] - p[0]) + math.cos(phi) * (p[1] - self.CoM[1])
  47 + self.G.vertex_properties["pos"][v] = p
  48 + else:
  49 + for v in self.G.vertices:
  50 + p = self.G.vertex_properties["pos"][v]
  51 + p[0] = self.CoM[0] + math.cos(phi) * (self.CoM[0] + p[0]) + math.sin(phi) * (p[1] - self.CoM[1])
  52 + p[1] = self.CoM[1] + math.sin(phi) * (self.CoM[0] + p[0]) - math.cos(phi) * (p[1] - self.CoM[1])
  53 + self.G.vertex_properties["pos"][v] = p
  54 +
  55 + def plot_graph(self, D, x, y):
  56 + plt.figure()
  57 + ext = [self.a[0], self.b[0], self.a[1], self.b[1]]
  58 + plt.imshow(D, origin = 'lower', extent=ext)
  59 + p = self.G.vertex_properties["pos"].get_2d_array(range(2)).T
  60 + plt.scatter(p[:,0], p[:,1], color='r')
  61 + plt.scatter(self.CoM[0], self.CoM[1], marker='*')
  62 +
  63 +# for n, contour in enumerate(self.cn):
  64 +# X = interp1d(np.arange(0, x.shape[0]), x)
  65 +# Y = interp1d(np.arange(0, y.shape[0]), y)
  66 +# contour[:, 1] = X(contour[:, 1])
  67 +# contour[: ,0] = Y(contour[:, 0])
  68 +# plt.plot(contour[:, 1], contour[:, 0])
  69 +# mx = np.amax(D)
  70 +# mn = np.amin(D)
  71 +# level = (mx-mn)/5.5
  72 +# cn = plt.contour(x, y, D, levels = [level])
  73 +
  74 + for e in self.G.edges():
  75 + coord = self.G.vertex_properties["pos"][e.source()]
  76 + coord2 = self.G.vertex_properties["pos"][e.target()]
  77 + X = [coord[0], coord2[0]]
  78 + Y = [coord[1], coord2[1]]
  79 + #all_plots.plot(x, y, 'go--', linewidth=1, markersize=1)
  80 + plt.plot(X, Y, 'go--', linewidth=1, markersize=1)
  81 +
  82 + plt.plot(*self.polygon.exterior.xy, color = 'r')
  83 + plt.show()
  84 +
  85 + def get_aabb(self):
  86 + pts = self.G.vertex_properties["pos"].get_2d_array(range(2)).T
  87 + a = np.asarray([100000.0, 100000.0])
  88 + b = np.asarray([-100000.0, -100000.0])
  89 +
  90 + #Find the bounding box based on the vertices.
  91 + for i in pts:
  92 + if(i[0] < a[0]):
  93 + a[0] = i[0]
  94 + if(i[1] < a[1]):
  95 + a[1] = i[1]
  96 + if(i[0] > b[0]):
  97 + b[0] = i[0]
  98 + if(i[1] > b[1]):
  99 + b[1] = i[1]
  100 +
  101 + #add 50% of the bounding box as padding on each side
  102 + d = 0.5*abs(a-b)
  103 + self.a = a - d
  104 + self.b = b + d
  105 +
  106 + def line(self, p1, p2, step1, step2):
  107 + return list(np.asarray(a) for a in zip(np.linspace(p1[0], p2[0], step1+1), np.linspace(p1[1], p2[1], step2+1)))
  108 +
  109 + def gen_polygon(self):
  110 + D, x, y = self.distancefield()
  111 + mx = np.amax(D)
  112 + mn = np.amin(D)
  113 + level = (mx-mn)/5.5
  114 + cn = measure.find_contours(D, level)
  115 + contour = copy.deepcopy(cn[0])
  116 + X = interp1d(np.arange(0, x.shape[0]), x)
  117 + Y = interp1d(np.arange(0, y.shape[0]), y)
  118 + contour[:, 0] = X(cn[0][:, 1])
  119 + contour[: ,1] = Y(cn[0][:, 0])
  120 + self.polygon = Polygon(contour)
  121 + self.CoM = self.centroid_com(contour)
  122 +
  123 + #cn = plt.contour(x, y, D, levels = [level])
  124 + #cn = plt.contour(x, y, D, levels = [level])
  125 +# plt.close()
  126 + #p = cn.collections[0].get_paths()[0]
  127 +# for i in range(len(cn.allsegs[0])):
  128 +#
  129 +# self.p = p
  130 +# v = p.vertices
  131 +# x = v[:, 0]
  132 +# y = v[:, 1]
  133 +# pts = np.array(zip(x, y))
  134 +# #nlist = c.trace(level, level, 0)
  135 +# #segs = nlist[:len(nlist)//2]
  136 + #self.polygon = Polygon(pts)
  137 + self.plot_graph(D, x, y)
  138 +
  139 +
  140 +
  141 + def distancefield(self):
  142 +
  143 + #generate a meshgrid of the appropriate size and resolution to surround the network
  144 + #get the space occupied by the network
  145 + lower = self.a
  146 + upper = self.b
  147 + R = np.asarray(np.floor(abs(lower-upper)), dtype=np.int)
  148 + if(R[0] < 10):
  149 + R[0] = 10
  150 + if(R[1] < 10):
  151 + R[1] = 10
  152 + x = np.linspace(lower[0], upper[0], R[0]) #get the grid points for uniform sampling of this space
  153 + y = np.linspace(lower[1], upper[1], R[1])
  154 + X, Y = np.meshgrid(x, y)
  155 + #Z = 150 * numpy.ones(X.shape)
  156 +
  157 + Q = np.stack((X, Y), 2)
  158 + d_x = abs(x[1]-x[0]);
  159 + d_y = abs(y[1]-y[0]);
  160 + dis1 = math.sqrt(pow(d_x,2)+pow(d_y,2))
  161 + #dx = abs(x[1]-x[0])
  162 +
  163 + #dy = abs(y[1]-y[0])
  164 + #dz = abs(z[1]-z[0])
  165 + #get a list of all node positions in the network
  166 + P = []
  167 +
  168 + for e in self.G.edges(): #12-17
  169 + start = self.G.vertex_properties["pos"][e.source()]
  170 + end = self.G.vertex_properties["pos"][e.target()]
  171 + l = self.line(start, end, 10, 10)
  172 + P = P + l
  173 +
  174 + for j in range(len(l)-1):
  175 + d_t = l[j+1]-l[j]
  176 + dis2 = math.sqrt(pow(d_t[0],2)+pow(d_t[1],2))
  177 + ins = max(int(d_t[0]/d_x), int(d_t[1]/d_y))
  178 + if(ins > 0):
  179 + ins = ins+1
  180 + for k in range(ins):
  181 + p_ins =l[j]+(k+1)*(l[j+1]-l[j])/ins
  182 + P.append(p_ins)
  183 + #turn that list into a Numpy array so that we can create a KD tree
  184 + P = np.array(P)
  185 +
  186 + #generate a KD-Tree out of the network point array
  187 + tree = sp.spatial.cKDTree(P)
  188 +
  189 + #specify the resolution of the ouput grid
  190 + # R = (200, 200, 200)
  191 +
  192 + D, I = tree.query(Q)
  193 +
  194 + return D, x, y
  195 +
  196 + def centroid_com(self, vertices):
  197 + # Polygon's signed area
  198 + A = 0
  199 + # Centroid's x
  200 + C_x = 0
  201 + # Centroid's y
  202 + C_y = 0
  203 + for i in range(0, len(vertices) - 1):
  204 + s = (vertices[i, 0] * vertices[i + 1, 1] - vertices[i + 1, 0] * vertices[i, 1])
  205 + A = A + s
  206 + C_x = C_x + (vertices[i, 0] + vertices[i + 1, 0]) * s
  207 + C_y = C_y + (vertices[i, 1] + vertices[i + 1, 1]) * s
  208 + A = 0.5 * A
  209 + C_x = (1.0 / (6.0 * A)) * C_x
  210 + C_y = (1.0 / (6.0 * A)) * C_y
  211 +
  212 + return np.array([C_x, C_y])
  213 +
  214 +
  215 +
  216 +def voronoi_polygons(voronoi, diameter):
  217 + """Generate shapely.geometry.Polygon objects corresponding to the
  218 + regions of a scipy.spatial.Voronoi object, in the order of the
  219 + input points. The polygons for the infinite regions are large
  220 + enough that all points within a distance 'diameter' of a Voronoi
  221 + vertex are contained in one of the infinite polygons.
  222 +
  223 + """
  224 + centroid = voronoi.points.mean(axis=0)
  225 +
  226 + # Mapping from (input point index, Voronoi point index) to list of
  227 + # unit vectors in the directions of the infinite ridges starting
  228 + # at the Voronoi point and neighbouring the input point.
  229 + ridge_direction = defaultdict(list)
  230 + for (p, q), rv in zip(voronoi.ridge_points, voronoi.ridge_vertices):
  231 + u, v = sorted(rv)
  232 + if u == -1:
  233 + # Infinite ridge starting at ridge point with index v,
  234 + # equidistant from input points with indexes p and q.
  235 + t = voronoi.points[q] - voronoi.points[p] # tangent
  236 + n = np.array([-t[1], t[0]]) / np.linalg.norm(t) # normal
  237 + midpoint = voronoi.points[[p, q]].mean(axis=0)
  238 + direction = np.sign(np.dot(midpoint - centroid, n)) * n
  239 + ridge_direction[p, v].append(direction)
  240 + ridge_direction[q, v].append(direction)
  241 +
  242 + for i, r in enumerate(voronoi.point_region):
  243 + region = voronoi.regions[r]
  244 + if -1 not in region:
  245 + # Finite region.
  246 + yield Polygon(voronoi.vertices[region])
  247 + continue
  248 + # Infinite region.
  249 + inf = region.index(-1) # Index of vertex at infinity.
  250 + j = region[(inf - 1) % len(region)] # Index of previous vertex.
  251 + k = region[(inf + 1) % len(region)] # Index of next vertex.
  252 + if j == k:
  253 + # Region has one Voronoi vertex with two ridges.
  254 + dir_j, dir_k = ridge_direction[i, j]
  255 + else:
  256 + # Region has two Voronoi vertices, each with one ridge.
  257 + dir_j, = ridge_direction[i, j]
  258 + dir_k, = ridge_direction[i, k]
  259 +
  260 + # Length of ridges needed for the extra edge to lie at least
  261 + # 'diameter' away from all Voronoi vertices.
  262 + length = 2 * diameter / np.linalg.norm(dir_j + dir_k)
  263 +
  264 + # Polygon consists of finite part plus an extra edge.
  265 + finite_part = voronoi.vertices[region[inf + 1:] + region[:inf]]
  266 + extra_edge = [voronoi.vertices[j] + dir_j * length,
  267 + voronoi.vertices[k] + dir_k * length]
  268 + yield Polygon(np.concatenate((finite_part, extra_edge)))
  269 +
  270 +
  271 +def load_nwt(filepath):
  272 + net = nwt.Network(filepath)
  273 + G = net.createFullGraph_gt()
  274 + G = net.filterDisconnected(G)
  275 + color = np.zeros(4, dtype = np.double)
  276 + color = [0.0, 1.0, 0.0, 1.0]
  277 + G.edge_properties["RGBA"] = G.new_edge_property("vector<double>", val=color)
  278 + color = [1.0, 0.0, 0.0, 0.9]
  279 + G.vertex_properties["RGBA"] = G.new_vertex_property("vector<double>", val=color)
  280 + bbl, bbu = net.aabb()
  281 +
  282 + return G, bbl, bbu
  283 +
  284 +def gen_cluster_graph(G, num_clusters, cluster_pos):
  285 + #create a graph that stores the edges of between the clusters
  286 + G_cluster = nwt.gt.Graph(directed=False)
  287 + G_cluster.vertex_properties["pos"] = G_cluster.new_vertex_property("vector<double>", val=np.zeros((3,1), dtype=np.float32))
  288 + G_cluster.vertex_properties["RGBA"] = G_cluster.new_vertex_property("vector<double>", val=np.zeros((4,1), dtype=np.float32))
  289 + for v in range(num_clusters):
  290 + G_cluster.add_vertex()
  291 + G_cluster.vertex_properties["pos"][G_cluster.vertex(v)] = np.asarray(cluster_pos[v], dtype=np.float32)
  292 + G_cluster.edge_properties["weight"] = G_cluster.new_edge_property("int", val = 0)
  293 + G_cluster.edge_properties["volume"] = G_cluster.new_edge_property("float", val = 0.0)
  294 + #for each edge in the original graph, generate appropriate subgraph edges without repretiions
  295 + #i.e. controls the thichness of the edges in the subgraph view.
  296 + for e in G.edges():
  297 + #if the source and target cluster is not equal to each other
  298 + #add an inter subgraph edge.
  299 + if(G.vertex_properties["clusters"][e.source()] != G.vertex_properties["clusters"][e.target()]):
  300 + t0 = e.source()
  301 + t1 = e.target()
  302 + ct0 = G_cluster.vertex(G.vertex_properties["clusters"][t0])
  303 + ct1 = G_cluster.vertex(G.vertex_properties["clusters"][t1])
  304 + if(G_cluster.edge(ct0, ct1) == None):
  305 + if(G_cluster.edge(ct1, ct0) == None):
  306 + #temp_e.append([G.vertex_properties["clusters"][e.source()], G.vertex_properties["clusters"][e.target()]])
  307 + G_cluster.add_edge(G_cluster.vertex(G.vertex_properties["clusters"][t0]), \
  308 + G_cluster.vertex(G.vertex_properties["clusters"][t1]))
  309 + G_cluster.edge_properties["weight"][G_cluster.edge(G_cluster.vertex(G.vertex_properties["clusters"][t0]), \
  310 + G_cluster.vertex(G.vertex_properties["clusters"][t1]))] += 1
  311 + G_cluster.edge_properties["volume"][G_cluster.edge(G_cluster.vertex(G.vertex_properties["clusters"][t0]), \
  312 + G_cluster.vertex(G.vertex_properties["clusters"][t1]))] += G.edge_properties["volume"][e]
  313 + G_cluster.vertex_properties["RGBA"][G_cluster.vertex(G.vertex_properties["clusters"][t0])] \
  314 + = G.vertex_properties["RGBA"][t0]
  315 + G_cluster.vertex_properties["RGBA"][G_cluster.vertex(G.vertex_properties["clusters"][t1])] \
  316 + = G.vertex_properties["RGBA"][t1]
  317 + else:
  318 + G_cluster.edge_properties["weight"][G_cluster.edge(G_cluster.vertex(G.vertex_properties["clusters"][t1]), \
  319 + G_cluster.vertex(G.vertex_properties["clusters"][t0]))] += 1
  320 + G_cluster.edge_properties["volume"][G_cluster.edge(G_cluster.vertex(G.vertex_properties["clusters"][t1]), \
  321 + G_cluster.vertex(G.vertex_properties["clusters"][t0]))] += G.edge_properties["volume"][e]
  322 + G_cluster.vertex_properties["RGBA"][G_cluster.vertex(G.vertex_properties["clusters"][t1])] \
  323 + = G.vertex_properties["RGBA"][t1]
  324 + G_cluster.vertex_properties["RGBA"][G_cluster.vertex(G.vertex_properties["clusters"][t0])] \
  325 + = G.vertex_properties["RGBA"][t0]
  326 + else:
  327 + G_cluster.edge_properties["weight"][G_cluster.edge(G_cluster.vertex(G.vertex_properties["clusters"][t0]), \
  328 + G_cluster.vertex(G.vertex_properties["clusters"][t1]))] += 1
  329 + G_cluster.edge_properties["volume"][G_cluster.edge(G_cluster.vertex(G.vertex_properties["clusters"][t0]), \
  330 + G_cluster.vertex(G.vertex_properties["clusters"][t1]))] += G.edge_properties["volume"][e]
  331 + G_cluster.vertex_properties["RGBA"][G_cluster.vertex(G.vertex_properties["clusters"][t0])] \
  332 + = G.vertex_properties["RGBA"][t0]
  333 + G_cluster.vertex_properties["RGBA"][G_cluster.vertex(G.vertex_properties["clusters"][t1])] \
  334 + = G.vertex_properties["RGBA"][t1]
  335 + G_cluster.vertex_properties["degree"] = G_cluster.degree_property_map("total")
  336 + vbetweeness_centrality = G_cluster.new_vertex_property("double")
  337 + ebetweeness_centrality = G_cluster.new_edge_property("double")
  338 + nwt.gt.graph_tool.centrality.betweenness(G_cluster, vprop=vbetweeness_centrality, eprop=ebetweeness_centrality)
  339 + ebc = ebetweeness_centrality.get_array()/0.01
  340 + G_cluster.vertex_properties["bc"] = vbetweeness_centrality
  341 + G_cluster.edge_properties["bc"] = ebetweeness_centrality
  342 + G_cluster.edge_properties["bc_scaled"] = G_cluster.new_edge_property("double", vals=ebc)
  343 +
  344 + dg = G_cluster.vertex_properties["degree"].get_array()
  345 + dg = 2*max(dg) - dg
  346 + d = G_cluster.new_vertex_property("int", vals=dg)
  347 + G_cluster.vertex_properties["10-degree"] = d
  348 +
  349 + return G_cluster
  350 +
  351 +
  352 +
  353 +
  354 +def gen_clusters(G, bbl, bbu, n_c = 20, edge_metric = 'volume', vertex_metric = 'degree'):
  355 +
  356 + #Generate the clusters
  357 + labels = nwt.Network.spectral_clustering(G,'length', n_clusters = n_c)
  358 + #self.labels = nwt.Network.spectral_clustering(G,'length')
  359 +
  360 + #Add clusters as a vertex property
  361 + G.vertex_properties["clusters"] = G.new_vertex_property("int", vals=labels)
  362 + G.vertex_properties["idx"] = G.vertex_index
  363 +
  364 + #gen bc metric
  365 + vbetweeness_centrality = G.new_vertex_property("double")
  366 + ebetweeness_centrality = G.new_edge_property("double")
  367 + nwt.gt.graph_tool.centrality.betweenness(G, vprop=vbetweeness_centrality, eprop=ebetweeness_centrality, norm=True)
  368 + G.vertex_properties["bc"] = vbetweeness_centrality
  369 + G.edge_properties["bc"] = ebetweeness_centrality
  370 +
  371 + num_clusters = len(np.unique(labels))
  372 +
  373 + #add colormap
  374 + G.vertex_properties["RGBA"] = nwt.Network.map_property_to_color(G, G.vertex_properties["clusters"])
  375 + temp_pos = []
  376 + for i in range(num_clusters):
  377 + num_v_in_cluster = len(np.argwhere(labels == i))
  378 + vfilt = np.zeros([G.num_vertices(), 1], dtype="bool")
  379 + vfilt[np.argwhere(labels == i)] = 1
  380 + vfilt_prop = G.new_vertex_property("bool", vals = vfilt)
  381 + G.set_vertex_filter(vfilt_prop)
  382 +
  383 + #get the filtered properties
  384 + g = nwt.gt.Graph(G, prune=True, directed=False)
  385 + positions = g.vertex_properties["pos"].get_2d_array(range(3)).T
  386 + position = np.sum(positions, 0)/num_v_in_cluster
  387 + temp_pos.append(position)
  388 + G.clear_filters()
  389 +
  390 + return gen_cluster_graph(G, num_clusters, temp_pos), G
  391 +
  392 +
  393 +def gen_subclusters(G, G_cluster, i, reposition = False):
  394 + vfilt = np.zeros([G.num_vertices(), 1], dtype='bool')
  395 + labels = G.vertex_properties["clusters"].get_array()
  396 + num_v_in_cluster = len(np.argwhere(labels == i))
  397 + vfilt[np.argwhere(labels == i)] = 1
  398 + vfilt_prop = G.new_vertex_property("bool", vals = vfilt)
  399 + G.set_vertex_filter(vfilt_prop)
  400 +
  401 + g = nwt.gt.Graph(G, prune=True, directed=False)
  402 +
  403 +
  404 + if reposition == True:
  405 + vbetweeness_centrality = g.new_vertex_property("double")
  406 + ebetweeness_centrality = g.new_edge_property("double")
  407 + nwt.gt.graph_tool.centrality.betweenness(g, vprop=vbetweeness_centrality, eprop=ebetweeness_centrality, norm=True)
  408 + g.vertex_properties["bc"] = vbetweeness_centrality
  409 + g.edge_properties["bc"] = ebetweeness_centrality
  410 + g.vertex_properties["pos"] = nwt.gt.sfdp_layout(g, eweight = ebetweeness_centrality)
  411 +
  412 + positions = g.vertex_properties["pos"].get_2d_array(range(2)).T
  413 + center = np.sum(positions, 0)/num_v_in_cluster
  414 + G.clear_filters()
  415 + return g, center
  416 +
  417 +#def gen_hierarchical_layout(G, G_cluster):
  418 +
  419 +def gen_polygons(G_c, bb):
  420 + G_c.vertex_properties["region_idx"] = G_c.new_vertex_property("int")
  421 + pts = G_c.vertex_properties["pos"].get_2d_array(range(2)).T
  422 + bl = np.asarray([bb[0], bb[2]])
  423 + lx = bb[1]-bb[0]
  424 + ly = bb[3]-bb[2]
  425 + r = copy.deepcopy(bl)
  426 + t = copy.deepcopy(bl)
  427 + tr = copy.deepcopy(bl)
  428 + r[0] = r[0] + lx
  429 + t[1] = t[1] + ly
  430 + tr[0] = tr[0] + lx
  431 + tr[1] = tr[1] + ly
  432 +
  433 + boundary = np.asarray([bl, t, tr, r, bl])
  434 + diameter = np.linalg.norm(boundary.ptp(axis=0))
  435 + boundary_polygon = Polygon(boundary)
  436 + vor = Voronoi(pts)
  437 + polygons = []
  438 + idx = 0
  439 + for poly in voronoi_polygons(vor, diameter):
  440 + coords = np.array(poly.intersection(boundary_polygon).exterior.coords)
  441 + polygons.append([coords])
  442 + for v in G_c.vertices():
  443 + point = Point(G_c.vertex_properties["pos"][v])
  444 + if poly.contains(point):
  445 + G_c.vertex_properties["region_idx"][v] = idx
  446 + idx+=1
  447 + break
  448 + return G_c, polygons, vor
  449 +
  450 +def centroid_region(vertices):
  451 + # Polygon's signed area
  452 + A = 0
  453 + # Centroid's x
  454 + C_x = 0
  455 + # Centroid's y
  456 + C_y = 0
  457 + for i in range(0, len(vertices) - 1):
  458 + s = (vertices[i, 0] * vertices[i + 1, 1] - vertices[i + 1, 0] * vertices[i, 1])
  459 + A = A + s
  460 + C_x = C_x + (vertices[i, 0] + vertices[i + 1, 0]) * s
  461 + C_y = C_y + (vertices[i, 1] + vertices[i + 1, 1]) * s
  462 + A = 0.5 * A
  463 + C_x = (1.0 / (6.0 * A)) * C_x
  464 + C_y = (1.0 / (6.0 * A)) * C_y
  465 +
  466 + return np.array([C_x, C_y])
  467 +
  468 +
  469 +def gen_image(G, G_c, itr, bb_flag = False, bb = None, reposition = False):
  470 +#def gen_image(G, G_c, vor, vor_filtered):
  471 + #Draw the layout using graph-tool (for comparison)
  472 + title = "clusters.pdf"
  473 + nwt.gt.graph_draw(G_c, pos=G_c.vertex_properties["pos"], vertex_fill_color=G_c.vertex_index, output=title, bg_color=[0.0, 0.0, 0.0, 1.0], output_size=(1000,1000))
  474 +
  475 + #get points of the centers of every cluster
  476 + #generate voronoi region and plot it.
  477 + fig, ax = plt.subplots(3, 1, sharex='col', sharey='row')
  478 + fig.tight_layout()
  479 + grid = plt.GridSpec(3,1)
  480 + grid.update(wspace=0.025, hspace=0.2)
  481 + ax[0].axis('off')
  482 + ax[1].axis('off')
  483 + ax[2].axis('off')
  484 +
  485 + all_plots = fig.add_subplot(grid[0])
  486 + ax[0].set_title(itr)
  487 + no_links = fig.add_subplot(grid[1], sharey=all_plots, sharex=all_plots)
  488 + voronoi = fig.add_subplot(grid[2], sharey=all_plots, sharex=all_plots)
  489 + pts = G_c.vertex_properties["pos"].get_2d_array(range(2)).T
  490 + if bb_flag == False:
  491 + vor = Voronoi(pts)
  492 + voronoi_plot_2d(vor, all_plots)
  493 + voronoi_plot_2d(vor, no_links)
  494 + voronoi_plot_2d(vor, voronoi)
  495 + a = voronoi.get_ylim()
  496 + b = voronoi.get_xlim()
  497 + bb = np.array([b[0], b[1], a[0], a[1]])
  498 + G_c, regions, vor = gen_polygons(G_c, bb)
  499 + if bb_flag == True:
  500 + voronoi_plot_2d(vor, all_plots)
  501 + voronoi_plot_2d(vor, no_links)
  502 + voronoi_plot_2d(vor, voronoi)
  503 + #plot the top-level graph
  504 + pts = G_c.vertex_properties["pos"].get_2d_array(range(2)).T
  505 + all_plots.scatter(pts[:,0], pts[:, 1], s=20*G_c.vertex_properties["degree"].get_array(), marker="*")
  506 + no_links.scatter(pts[:,0], pts[:, 1], s=20*G_c.vertex_properties["degree"].get_array(), marker="*")
  507 + voronoi.scatter(pts[:,0], pts[:, 1], s=20*G_c.vertex_properties["degree"].get_array(), marker="*")
  508 + #plot the connections of the top level graph
  509 + for e in G_c.edges():
  510 + coord = G_c.vertex_properties["pos"][e.source()]
  511 + coord2 = G_c.vertex_properties["pos"][e.target()]
  512 + x = [coord[0], coord2[0]]
  513 + y = [coord[1], coord2[1]]
  514 + #all_plots.plot(x, y, 'go--', linewidth=1, markersize=1)
  515 + no_links.plot(x, y, 'go--', linewidth=1, markersize=1)
  516 + voronoi.plot(x, y, 'go--', linewidth=1, markersize=1)
  517 +
  518 + #for every subgraph generate a layout and plot the result
  519 + for i in range(num_clusters):
  520 + g, center = gen_subclusters(G, G_c, i, reposition)
  521 + d = G_c.vertex_properties["pos"][i] - center
  522 + t = Polygon_mass(g)
  523 + #t.distancefield()
  524 + for v in g.vertices():
  525 + G.vertex_properties["pos"][g.vertex_properties["idx"][v]] = g.vertex_properties["pos"][v] + d
  526 + #g.vertex_properties["pos"][g.vertex_properties["idx"][v]] = g.vertex_properties["pos"][v] + d
  527 + #sub_pts = g.vertex_properties["pos"].get_2d_array(range(2)).T
  528 + #all_plots.scatter(pts[:,0], pts[:, 1], marker="*")
  529 + # for e in g.edges():
  530 + # coord = g.vertex_properties["pos"][e.source()]
  531 + # coord2 = g.vertex_properties["pos"][e.target()]
  532 + # x = [coord[0], coord2[0]]
  533 + # y = [coord[1], coord2[1]]
  534 + # plt.plot(x, y, 'ro--', linewidth=1, markersize=1)
  535 +
  536 + for e in G.edges():
  537 + coord = G.vertex_properties["pos"][e.source()]
  538 + coord2 = G.vertex_properties["pos"][e.target()]
  539 + x = [coord[0], coord2[0]]
  540 + y = [coord[1], coord2[1]]
  541 + if (G.vertex_properties["clusters"][e.source()] == G.vertex_properties["clusters"][e.target()]):
  542 + all_plots.plot(x, y, 'ro--', linewidth=1, markersize=1)
  543 + no_links.plot(x, y, 'ro--', linewidth=1, markersize=1)
  544 + else:
  545 + all_plots.plot(x, y, 'bo--', linewidth=1, markersize=1)
  546 +
  547 + no_links.xaxis.set_visible(False)
  548 + all_plots.xaxis.set_visible(False)
  549 + for v in G_c.vertices():
  550 + region = regions[G_c.vertex_properties["region_idx"][v]]
  551 + centroid = centroid_region(region[0])
  552 + G_c.vertex_properties["pos"][v] = centroid
  553 +
  554 + pts_temp = G_c.vertex_properties["pos"].get_2d_array(range(2)).T
  555 + all_plots.scatter(pts_temp[:,0], pts_temp[:, 1], marker='.', color='r')
  556 + no_links.scatter(pts_temp[:,0], pts_temp[:, 1], marker='.', color='r')
  557 + voronoi.scatter(pts_temp[:,0], pts_temp[:, 1], marker='.', color='r')
  558 +
  559 + all_plots.set_xlim([bb[0], bb[1]])
  560 + all_plots.set_ylim([bb[2], bb[3]])
  561 +
  562 + no_links.set_xlim([bb[0], bb[1]])
  563 + no_links.set_ylim([bb[2], bb[3]])
  564 +
  565 + voronoi.set_xlim([bb[0], bb[1]])
  566 + voronoi.set_ylim([bb[2], bb[3]])
  567 +
  568 + plt.show()
  569 +
  570 +
  571 + return G, G_c, bb
  572 +
  573 +
  574 +
  575 +
  576 +
  577 +
  578 +
  579 +
  580 +#G_c.vertex_properties["pos"] = nwt.gt.fruchterman_reingold_layout(G_c, weight=G_c.edge_properties["weight"], r=G_c.num_vertices()*0.1, a = G_c.num_vertices()*500)
  581 +
  582 +G, bbl, bbu = load_nwt("/home/pavel/Documents/Python/GraphGuiQt/network_4.nwt")
  583 +G_c, G = gen_clusters(G, bbl, bbu)
  584 +num_clusters = 20
  585 +
  586 +#G_c.vertex_properties["pos"] = nwt.gt.radial_tree_layout(G_c, root=np.argwhere(G_c.vertex_properties["degree"].get_array() == max(G_c.vertex_properties["degree"].get_array())), node_weight = G_c.vertex_properties["10-degree"], r= 2.0)
  587 +G_c.vertex_properties["pos"] = nwt.gt.sfdp_layout(G_c, eweight=G_c.edge_properties["volume"], vweight=G_c.vertex_properties["degree"], C = 1.0, K = 10)
  588 +
  589 +G, G_c, bb = gen_image(G, G_c, "base", reposition = True)
  590 +itr = 0
  591 +
  592 +#for itr in range(5):
  593 +# G, G_c, bb = gen_image(G, G_c, itr, True, bb)
  594 +# itr+=1
  595 +
  596 +g, center = gen_subclusters(G, G_c)
  597 +#d = G_c.vertex_properties["pos"][0] - center
  598 +#for v in g.vertices():
  599 +# g.vertex_properties["pos"][v] = g.vertex_properties["pos"][v] + d
  600 +
  601 +#G_c = nwt.Network.gen_new_fd_layout(G_c)
  602 +#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)
  603 +
  604 +
  605 +
  606 +