From 193cb4c6c21380f415e8a5ef6c16c4b42980e090 Mon Sep 17 00:00:00 2001 From: Pavel Govyadinov Date: Thu, 12 Sep 2019 18:48:34 -0500 Subject: [PATCH] need this to test --- GraphCanvas.py | 4 ++-- GuiVisPy_tube.py | 2 +- TubeCanvas.py | 49 +++++++++++++++++++++++++++++++++++++++++++++++-- TubeWidget.py | 28 +++++++++++++++++++++++++++- graph_shaders.py | 4 ++-- voronoi_test.py | 606 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 6 files changed, 685 insertions(+), 8 deletions(-) create mode 100644 voronoi_test.py diff --git a/GraphCanvas.py b/GraphCanvas.py index 15ca17b..c1f96f6 100644 --- a/GraphCanvas.py +++ b/GraphCanvas.py @@ -828,8 +828,8 @@ class GraphCanvas(scene.SceneCanvas): #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 = []; + 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) diff --git a/GuiVisPy_tube.py b/GuiVisPy_tube.py index ddfafec..48ca01b 100644 --- a/GuiVisPy_tube.py +++ b/GuiVisPy_tube.py @@ -93,7 +93,7 @@ def load_nwt(filepath): 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") +#G, bbl, bbu = load_nwt("/home/pavel/Documents/Python/GraphGuiQt/net`````````work_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) diff --git a/TubeCanvas.py b/TubeCanvas.py index 60ffac1..b430c1b 100644 --- a/TubeCanvas.py +++ b/TubeCanvas.py @@ -53,6 +53,10 @@ class TubeDraw(scene.SceneCanvas): high=(4-1)).astype(np.uint32) self.vbo = gloo.VertexBuffer(self.cylinder_data) self.triangles = gloo.IndexBuffer(self.triangle_data) + + #generate an index buffer for the caps. + self.cap_data = np.random.randint(size=(5, 3), low=0, high=(4-1)).astype(np.uint32) + self.caps = gloo.IndexBuffer(self.cap_data) self.program.bind(self.vbo) self.scale = [1,1,1] self.r1 = np.eye(4, dtype=np.float32) @@ -61,7 +65,7 @@ class TubeDraw(scene.SceneCanvas): #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'), cull_face='back') + 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.] @@ -93,7 +97,9 @@ class TubeDraw(scene.SceneCanvas): #create program self.gen_cylinder_vbo(self.G, self.num_sides) self.vbo = gloo.VertexBuffer(self.cylinder_data) + #self.triangle_data = np.append(self.triangle_data, self.cap_data, axis=0) self.triangles = gloo.IndexBuffer(self.triangle_data) + self.caps = gloo.IndexBuffer(self.cap_data) #self.view = np.eye(4, dtype=np.float32) self.model = np.eye(4, dtype=np.float32) @@ -153,6 +159,14 @@ class TubeDraw(scene.SceneCanvas): set_state(clear_color='white', depth_test=True, blend=True, blend_func=('src_alpha', 'one_minus_src_alpha'), depth_func = ('lequal'), cull_face='back') self.program.draw('triangles', self.triangles) + if self.select: + set_state(clear_color='white', depth_test=True, blend=True, + blend_func=('src_alpha', 'one_minus_src_alpha'), depth_func = ('always'), cull_face=False) + else: + set_state(clear_color='white', depth_test=True, blend=True, + blend_func=('src_alpha', 'one_minus_src_alpha'), depth_func = ('lequal'), cull_face=False) + self.program.draw('triangles', self.caps) + self.projection = perspective(90.0, self.physical_size[0]/self.physical_size[1], 1.0, 1000.0) self.program['u_projection'] = self.projection @@ -161,9 +175,12 @@ class TubeDraw(scene.SceneCanvas): i = 0 num_pts = 0 num_tri = 0 + num_tri_caps = 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 + num_tri_caps += (2*(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), @@ -171,8 +188,12 @@ class TubeDraw(scene.SceneCanvas): ]) self.triangle_data = np.random.randint(size=(num_tri, 3), low=0, high=(G.num_edges()-1)).astype(np.uint32) + + self.cap_data = np.random.randint(size = (num_tri_caps, 3), + low = 0, high=(G.num_edges()-2)).astype(np.uint32) index = 0 t_index = 0 + c_index = 0 #for each edge generate a cylinder. for e in G.edges(): #print("done") @@ -266,9 +287,26 @@ class TubeDraw(scene.SceneCanvas): #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) - + cap_triangles = np.random.randint(size=((num_sides-2)*2 , 3), low=0, high=(num_sides-2)).astype(np.uint32) t_index_temp = 0 + c_index_temp = 0 + #for every ring in the fiber + for ring in range(0, pts.shape[0], pts.shape[0]-1): + #if we have the first ring or the last ring add a cap. + idx_ori = index+ring*num_sides + for side in range(0, num_sides-2): + idx_current_point = index+ring*num_sides + side + if(side < num_sides-2): + cap_tri = [idx_ori, idx_current_point+1, idx_current_point+2] + else: + cap_tri = [idx_ori, idx_current_point+1, idx_ori] + cap_triangles[c_index_temp] = cap_tri + self.cap_data[c_index] = cap_tri + c_index_temp += 1 + c_index += 1 + for ring in range(pts.shape[0]-1): + #otherwise generate the sides. for side in range(num_sides): if(side < num_sides-1): idx_current_point = index+ring*num_sides + side @@ -322,6 +360,13 @@ class TubeDraw(scene.SceneCanvas): 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(cap_triangles.shape[0]): + tri = np.zeros((3,4)) + tri[:,0] = self.cylinder_data['a_position'][cap_triangles[j][0]] + tri[:,1] = self.cylinder_data['a_position'][cap_triangles[j][1]] + tri[:,2] = self.cylinder_data['a_position'][cap_triangles[j][2]] + tri[:,3] = self.cylinder_data['a_position'][cap_triangles[j][0]] + ax.plot(tri[0,:], tri[1,:], tri[2,:], c='b') for j in range(triangles.shape[0]): tri = np.zeros((3,4)) tri[:,0] = self.cylinder_data['a_position'][triangles[j][0]] diff --git a/TubeWidget.py b/TubeWidget.py index b1df505..b451563 100644 --- a/TubeWidget.py +++ b/TubeWidget.py @@ -13,6 +13,8 @@ Created on Mon Aug 5 15:53:16 2019 from pyqtgraph.Qt import QtCore, QtGui, QtWidgets from TubeCanvas import TubeDraw +import numpy as np +import trimesh as tm DEBUG = False @@ -50,7 +52,13 @@ class TubeWidget(QtGui.QWidget): #Handles the mouse press events def on_mouse_press(self, event): self.down = True - n = 3 + if event.button == 2: + menu = QtWidgets.QMenu(self) + NS = menu.addAction('Export Model') + action = menu.exec_(event.native.globalPos()) + if action == NS: + self.save_stl_mesh() + #Handles the mouse wheel event def on_mouse_wheel(self, event): @@ -86,3 +94,21 @@ class TubeWidget(QtGui.QWidget): """ def connect(self, signal_object): signal_object.select.connect(self.select) + + def save_stl_mesh(self): + triangles = self.canvas.triangle_data + caps = self.canvas.cap_data + points = self.canvas.cylinder_data['a_position'] + colors = self.canvas.cylinder_data['a_fg_color'] + normals = self.canvas.cylinder_data['a_normal'] + + mesh = tm.Trimesh(vertices=points, faces=np.append(triangles, caps, axis=0), + vertex_colors=colors, vertex_normals=normals) + mesh.export('with_caps.obj') + mesh.export('with_caps.ply') + + + + + + diff --git a/graph_shaders.py b/graph_shaders.py index cc799b2..d2b6a5f 100644 --- a/graph_shaders.py +++ b/graph_shaders.py @@ -94,7 +94,7 @@ void main() // 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)) + if( r > (v_linewidth/3.0+v_antialias)) { discard; } @@ -113,7 +113,7 @@ void main() 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; + gl_FragColor = vec4(v_fg_color.rgb, v_bg_color.a); } } else diff --git a/voronoi_test.py b/voronoi_test.py new file mode 100644 index 0000000..98857e1 --- /dev/null +++ b/voronoi_test.py @@ -0,0 +1,606 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Created on Tue Sep 3 13:15:54 2019 + +@author: pavel +""" + +from scipy.spatial import Voronoi, voronoi_plot_2d +from scipy.interpolate import interp1d +#import matplotlib._cntr as cntr +from shapely.geometry import Point +from shapely.geometry import Polygon +import numpy as np +import scipy as sp +import math +import matplotlib.pyplot as plt +import sys +import copy + +from skimage import measure + +from collections import defaultdict + +import network_dep as nwt + +class Polygon_mass: + def __init__(self, G): + self.G = G + self.get_aabb() + self.gen_polygon() + self.torque = [] + + def add_torque(self, p, f): + d = self.CoM - p + r = np.linalg.norm(self.CoM - p) + theta = math.acos(np.dot(d, f)/np.dot(d, d)/np.dot(f, f)) + torque = r * math.sin(theta) * f + self.torque.append(torque) + + def rotate(self, phi, direction = "counterclock"): + if("counterclock"): + for v in self.G.vertices: + p = self.G.vertex_properties["pos"][v] + p[0] = self.CoM[0] + math.cos(phi) * (self.CoM[0] - p[0]) - math.sin(phi) * (p[1] - self.CoM[1]) + p[1] = self.CoM[1] + math.sin(phi) * (self.CoM[0] - p[0]) + math.cos(phi) * (p[1] - self.CoM[1]) + self.G.vertex_properties["pos"][v] = p + else: + for v in self.G.vertices: + p = self.G.vertex_properties["pos"][v] + p[0] = self.CoM[0] + math.cos(phi) * (self.CoM[0] + p[0]) + math.sin(phi) * (p[1] - self.CoM[1]) + p[1] = self.CoM[1] + math.sin(phi) * (self.CoM[0] + p[0]) - math.cos(phi) * (p[1] - self.CoM[1]) + self.G.vertex_properties["pos"][v] = p + + def plot_graph(self, D, x, y): + plt.figure() + ext = [self.a[0], self.b[0], self.a[1], self.b[1]] + plt.imshow(D, origin = 'lower', extent=ext) + p = self.G.vertex_properties["pos"].get_2d_array(range(2)).T + plt.scatter(p[:,0], p[:,1], color='r') + plt.scatter(self.CoM[0], self.CoM[1], marker='*') + +# for n, contour in enumerate(self.cn): +# X = interp1d(np.arange(0, x.shape[0]), x) +# Y = interp1d(np.arange(0, y.shape[0]), y) +# contour[:, 1] = X(contour[:, 1]) +# contour[: ,0] = Y(contour[:, 0]) +# plt.plot(contour[:, 1], contour[:, 0]) +# mx = np.amax(D) +# mn = np.amin(D) +# level = (mx-mn)/5.5 +# cn = plt.contour(x, y, D, levels = [level]) + + for e in self.G.edges(): + coord = self.G.vertex_properties["pos"][e.source()] + coord2 = self.G.vertex_properties["pos"][e.target()] + X = [coord[0], coord2[0]] + Y = [coord[1], coord2[1]] + #all_plots.plot(x, y, 'go--', linewidth=1, markersize=1) + plt.plot(X, Y, 'go--', linewidth=1, markersize=1) + + plt.plot(*self.polygon.exterior.xy, color = 'r') + plt.show() + + def get_aabb(self): + pts = self.G.vertex_properties["pos"].get_2d_array(range(2)).T + a = np.asarray([100000.0, 100000.0]) + b = np.asarray([-100000.0, -100000.0]) + + #Find the bounding box based on the vertices. + for i in pts: + if(i[0] < a[0]): + a[0] = i[0] + if(i[1] < a[1]): + a[1] = i[1] + if(i[0] > b[0]): + b[0] = i[0] + if(i[1] > b[1]): + b[1] = i[1] + + #add 50% of the bounding box as padding on each side + d = 0.5*abs(a-b) + self.a = a - d + self.b = b + d + + def line(self, p1, p2, step1, step2): + 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))) + + def gen_polygon(self): + D, x, y = self.distancefield() + mx = np.amax(D) + mn = np.amin(D) + level = (mx-mn)/5.5 + cn = measure.find_contours(D, level) + contour = copy.deepcopy(cn[0]) + X = interp1d(np.arange(0, x.shape[0]), x) + Y = interp1d(np.arange(0, y.shape[0]), y) + contour[:, 0] = X(cn[0][:, 1]) + contour[: ,1] = Y(cn[0][:, 0]) + self.polygon = Polygon(contour) + self.CoM = self.centroid_com(contour) + + #cn = plt.contour(x, y, D, levels = [level]) + #cn = plt.contour(x, y, D, levels = [level]) +# plt.close() + #p = cn.collections[0].get_paths()[0] +# for i in range(len(cn.allsegs[0])): +# +# self.p = p +# v = p.vertices +# x = v[:, 0] +# y = v[:, 1] +# pts = np.array(zip(x, y)) +# #nlist = c.trace(level, level, 0) +# #segs = nlist[:len(nlist)//2] + #self.polygon = Polygon(pts) + self.plot_graph(D, x, y) + + + + def distancefield(self): + + #generate a meshgrid of the appropriate size and resolution to surround the network + #get the space occupied by the network + lower = self.a + upper = self.b + R = np.asarray(np.floor(abs(lower-upper)), dtype=np.int) + if(R[0] < 10): + R[0] = 10 + if(R[1] < 10): + R[1] = 10 + 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]) + X, Y = np.meshgrid(x, y) + #Z = 150 * numpy.ones(X.shape) + + Q = np.stack((X, Y), 2) + d_x = abs(x[1]-x[0]); + d_y = abs(y[1]-y[0]); + dis1 = math.sqrt(pow(d_x,2)+pow(d_y,2)) + #dx = abs(x[1]-x[0]) + + #dy = abs(y[1]-y[0]) + #dz = abs(z[1]-z[0]) + #get a list of all node positions in the network + P = [] + + for e in self.G.edges(): #12-17 + start = self.G.vertex_properties["pos"][e.source()] + end = self.G.vertex_properties["pos"][e.target()] + l = self.line(start, end, 10, 10) + P = P + l + + for j in range(len(l)-1): + d_t = l[j+1]-l[j] + dis2 = math.sqrt(pow(d_t[0],2)+pow(d_t[1],2)) + ins = max(int(d_t[0]/d_x), int(d_t[1]/d_y)) + if(ins > 0): + ins = ins+1 + for k in range(ins): + p_ins =l[j]+(k+1)*(l[j+1]-l[j])/ins + P.append(p_ins) + #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) + + #specify the resolution of the ouput grid + # R = (200, 200, 200) + + D, I = tree.query(Q) + + return D, x, y + + def centroid_com(self, vertices): + # Polygon's signed area + A = 0 + # Centroid's x + C_x = 0 + # Centroid's y + C_y = 0 + for i in range(0, len(vertices) - 1): + s = (vertices[i, 0] * vertices[i + 1, 1] - vertices[i + 1, 0] * vertices[i, 1]) + A = A + s + C_x = C_x + (vertices[i, 0] + vertices[i + 1, 0]) * s + C_y = C_y + (vertices[i, 1] + vertices[i + 1, 1]) * s + A = 0.5 * A + C_x = (1.0 / (6.0 * A)) * C_x + C_y = (1.0 / (6.0 * A)) * C_y + + return np.array([C_x, C_y]) + + + +def voronoi_polygons(voronoi, diameter): + """Generate shapely.geometry.Polygon objects corresponding to the + regions of a scipy.spatial.Voronoi object, in the order of the + input points. The polygons for the infinite regions are large + enough that all points within a distance 'diameter' of a Voronoi + vertex are contained in one of the infinite polygons. + + """ + centroid = voronoi.points.mean(axis=0) + + # Mapping from (input point index, Voronoi point index) to list of + # unit vectors in the directions of the infinite ridges starting + # at the Voronoi point and neighbouring the input point. + ridge_direction = defaultdict(list) + for (p, q), rv in zip(voronoi.ridge_points, voronoi.ridge_vertices): + u, v = sorted(rv) + if u == -1: + # Infinite ridge starting at ridge point with index v, + # equidistant from input points with indexes p and q. + t = voronoi.points[q] - voronoi.points[p] # tangent + n = np.array([-t[1], t[0]]) / np.linalg.norm(t) # normal + midpoint = voronoi.points[[p, q]].mean(axis=0) + direction = np.sign(np.dot(midpoint - centroid, n)) * n + ridge_direction[p, v].append(direction) + ridge_direction[q, v].append(direction) + + for i, r in enumerate(voronoi.point_region): + region = voronoi.regions[r] + if -1 not in region: + # Finite region. + yield Polygon(voronoi.vertices[region]) + continue + # Infinite region. + inf = region.index(-1) # Index of vertex at infinity. + j = region[(inf - 1) % len(region)] # Index of previous vertex. + k = region[(inf + 1) % len(region)] # Index of next vertex. + if j == k: + # Region has one Voronoi vertex with two ridges. + dir_j, dir_k = ridge_direction[i, j] + else: + # Region has two Voronoi vertices, each with one ridge. + dir_j, = ridge_direction[i, j] + dir_k, = ridge_direction[i, k] + + # Length of ridges needed for the extra edge to lie at least + # 'diameter' away from all Voronoi vertices. + length = 2 * diameter / np.linalg.norm(dir_j + dir_k) + + # Polygon consists of finite part plus an extra edge. + finite_part = voronoi.vertices[region[inf + 1:] + region[:inf]] + extra_edge = [voronoi.vertices[j] + dir_j * length, + voronoi.vertices[k] + dir_k * length] + yield Polygon(np.concatenate((finite_part, extra_edge))) + + +def load_nwt(filepath): + net = nwt.Network(filepath) + G = net.createFullGraph_gt() + G = net.filterDisconnected(G) + 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 + +def gen_cluster_graph(G, num_clusters, cluster_pos): + #create a graph that stores the edges of between the clusters + G_cluster = nwt.gt.Graph(directed=False) + G_cluster.vertex_properties["pos"] = G_cluster.new_vertex_property("vector", val=np.zeros((3,1), dtype=np.float32)) + G_cluster.vertex_properties["RGBA"] = G_cluster.new_vertex_property("vector", val=np.zeros((4,1), dtype=np.float32)) + for v in range(num_clusters): + G_cluster.add_vertex() + G_cluster.vertex_properties["pos"][G_cluster.vertex(v)] = np.asarray(cluster_pos[v], dtype=np.float32) + G_cluster.edge_properties["weight"] = G_cluster.new_edge_property("int", val = 0) + G_cluster.edge_properties["volume"] = G_cluster.new_edge_property("float", val = 0.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()]): + t0 = e.source() + t1 = e.target() + ct0 = G_cluster.vertex(G.vertex_properties["clusters"][t0]) + ct1 = G_cluster.vertex(G.vertex_properties["clusters"][t1]) + if(G_cluster.edge(ct0, ct1) == None): + if(G_cluster.edge(ct1, ct0) == None): + #temp_e.append([G.vertex_properties["clusters"][e.source()], G.vertex_properties["clusters"][e.target()]]) + G_cluster.add_edge(G_cluster.vertex(G.vertex_properties["clusters"][t0]), \ + G_cluster.vertex(G.vertex_properties["clusters"][t1])) + G_cluster.edge_properties["weight"][G_cluster.edge(G_cluster.vertex(G.vertex_properties["clusters"][t0]), \ + G_cluster.vertex(G.vertex_properties["clusters"][t1]))] += 1 + G_cluster.edge_properties["volume"][G_cluster.edge(G_cluster.vertex(G.vertex_properties["clusters"][t0]), \ + G_cluster.vertex(G.vertex_properties["clusters"][t1]))] += G.edge_properties["volume"][e] + G_cluster.vertex_properties["RGBA"][G_cluster.vertex(G.vertex_properties["clusters"][t0])] \ + = G.vertex_properties["RGBA"][t0] + G_cluster.vertex_properties["RGBA"][G_cluster.vertex(G.vertex_properties["clusters"][t1])] \ + = G.vertex_properties["RGBA"][t1] + else: + G_cluster.edge_properties["weight"][G_cluster.edge(G_cluster.vertex(G.vertex_properties["clusters"][t1]), \ + G_cluster.vertex(G.vertex_properties["clusters"][t0]))] += 1 + G_cluster.edge_properties["volume"][G_cluster.edge(G_cluster.vertex(G.vertex_properties["clusters"][t1]), \ + G_cluster.vertex(G.vertex_properties["clusters"][t0]))] += G.edge_properties["volume"][e] + G_cluster.vertex_properties["RGBA"][G_cluster.vertex(G.vertex_properties["clusters"][t1])] \ + = G.vertex_properties["RGBA"][t1] + G_cluster.vertex_properties["RGBA"][G_cluster.vertex(G.vertex_properties["clusters"][t0])] \ + = G.vertex_properties["RGBA"][t0] + else: + G_cluster.edge_properties["weight"][G_cluster.edge(G_cluster.vertex(G.vertex_properties["clusters"][t0]), \ + G_cluster.vertex(G.vertex_properties["clusters"][t1]))] += 1 + G_cluster.edge_properties["volume"][G_cluster.edge(G_cluster.vertex(G.vertex_properties["clusters"][t0]), \ + G_cluster.vertex(G.vertex_properties["clusters"][t1]))] += G.edge_properties["volume"][e] + G_cluster.vertex_properties["RGBA"][G_cluster.vertex(G.vertex_properties["clusters"][t0])] \ + = G.vertex_properties["RGBA"][t0] + G_cluster.vertex_properties["RGBA"][G_cluster.vertex(G.vertex_properties["clusters"][t1])] \ + = G.vertex_properties["RGBA"][t1] + G_cluster.vertex_properties["degree"] = G_cluster.degree_property_map("total") + vbetweeness_centrality = G_cluster.new_vertex_property("double") + ebetweeness_centrality = G_cluster.new_edge_property("double") + nwt.gt.graph_tool.centrality.betweenness(G_cluster, vprop=vbetweeness_centrality, eprop=ebetweeness_centrality) + ebc = ebetweeness_centrality.get_array()/0.01 + G_cluster.vertex_properties["bc"] = vbetweeness_centrality + G_cluster.edge_properties["bc"] = ebetweeness_centrality + G_cluster.edge_properties["bc_scaled"] = G_cluster.new_edge_property("double", vals=ebc) + + dg = G_cluster.vertex_properties["degree"].get_array() + dg = 2*max(dg) - dg + d = G_cluster.new_vertex_property("int", vals=dg) + G_cluster.vertex_properties["10-degree"] = d + + return G_cluster + + + + +def gen_clusters(G, bbl, bbu, n_c = 20, edge_metric = 'volume', vertex_metric = 'degree'): + + #Generate the clusters + 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=labels) + G.vertex_properties["idx"] = G.vertex_index + + #gen bc metric + vbetweeness_centrality = G.new_vertex_property("double") + ebetweeness_centrality = G.new_edge_property("double") + nwt.gt.graph_tool.centrality.betweenness(G, vprop=vbetweeness_centrality, eprop=ebetweeness_centrality, norm=True) + G.vertex_properties["bc"] = vbetweeness_centrality + G.edge_properties["bc"] = ebetweeness_centrality + + num_clusters = len(np.unique(labels)) + + #add colormap + G.vertex_properties["RGBA"] = nwt.Network.map_property_to_color(G, G.vertex_properties["clusters"]) + temp_pos = [] + for i in range(num_clusters): + num_v_in_cluster = len(np.argwhere(labels == i)) + vfilt = np.zeros([G.num_vertices(), 1], dtype="bool") + vfilt[np.argwhere(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 + temp_pos.append(position) + G.clear_filters() + + return gen_cluster_graph(G, num_clusters, temp_pos), G + + +def gen_subclusters(G, G_cluster, i, reposition = False): + vfilt = np.zeros([G.num_vertices(), 1], dtype='bool') + labels = G.vertex_properties["clusters"].get_array() + num_v_in_cluster = len(np.argwhere(labels == i)) + vfilt[np.argwhere(labels == i)] = 1 + vfilt_prop = G.new_vertex_property("bool", vals = vfilt) + G.set_vertex_filter(vfilt_prop) + + g = nwt.gt.Graph(G, prune=True, directed=False) + + + if reposition == True: + vbetweeness_centrality = g.new_vertex_property("double") + ebetweeness_centrality = g.new_edge_property("double") + nwt.gt.graph_tool.centrality.betweenness(g, vprop=vbetweeness_centrality, eprop=ebetweeness_centrality, norm=True) + g.vertex_properties["bc"] = vbetweeness_centrality + g.edge_properties["bc"] = ebetweeness_centrality + g.vertex_properties["pos"] = nwt.gt.sfdp_layout(g, eweight = ebetweeness_centrality) + + positions = g.vertex_properties["pos"].get_2d_array(range(2)).T + center = np.sum(positions, 0)/num_v_in_cluster + G.clear_filters() + return g, center + +#def gen_hierarchical_layout(G, G_cluster): + +def gen_polygons(G_c, bb): + G_c.vertex_properties["region_idx"] = G_c.new_vertex_property("int") + pts = G_c.vertex_properties["pos"].get_2d_array(range(2)).T + bl = np.asarray([bb[0], bb[2]]) + lx = bb[1]-bb[0] + ly = bb[3]-bb[2] + r = copy.deepcopy(bl) + t = copy.deepcopy(bl) + tr = copy.deepcopy(bl) + r[0] = r[0] + lx + t[1] = t[1] + ly + tr[0] = tr[0] + lx + tr[1] = tr[1] + ly + + boundary = np.asarray([bl, t, tr, r, bl]) + diameter = np.linalg.norm(boundary.ptp(axis=0)) + boundary_polygon = Polygon(boundary) + vor = Voronoi(pts) + polygons = [] + idx = 0 + for poly in voronoi_polygons(vor, diameter): + coords = np.array(poly.intersection(boundary_polygon).exterior.coords) + polygons.append([coords]) + for v in G_c.vertices(): + point = Point(G_c.vertex_properties["pos"][v]) + if poly.contains(point): + G_c.vertex_properties["region_idx"][v] = idx + idx+=1 + break + return G_c, polygons, vor + +def centroid_region(vertices): + # Polygon's signed area + A = 0 + # Centroid's x + C_x = 0 + # Centroid's y + C_y = 0 + for i in range(0, len(vertices) - 1): + s = (vertices[i, 0] * vertices[i + 1, 1] - vertices[i + 1, 0] * vertices[i, 1]) + A = A + s + C_x = C_x + (vertices[i, 0] + vertices[i + 1, 0]) * s + C_y = C_y + (vertices[i, 1] + vertices[i + 1, 1]) * s + A = 0.5 * A + C_x = (1.0 / (6.0 * A)) * C_x + C_y = (1.0 / (6.0 * A)) * C_y + + return np.array([C_x, C_y]) + + +def gen_image(G, G_c, itr, bb_flag = False, bb = None, reposition = False): +#def gen_image(G, G_c, vor, vor_filtered): + #Draw the layout using graph-tool (for comparison) + title = "clusters.pdf" + 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)) + + #get points of the centers of every cluster + #generate voronoi region and plot it. + fig, ax = plt.subplots(3, 1, sharex='col', sharey='row') + fig.tight_layout() + grid = plt.GridSpec(3,1) + grid.update(wspace=0.025, hspace=0.2) + ax[0].axis('off') + ax[1].axis('off') + ax[2].axis('off') + + all_plots = fig.add_subplot(grid[0]) + ax[0].set_title(itr) + no_links = fig.add_subplot(grid[1], sharey=all_plots, sharex=all_plots) + voronoi = fig.add_subplot(grid[2], sharey=all_plots, sharex=all_plots) + pts = G_c.vertex_properties["pos"].get_2d_array(range(2)).T + if bb_flag == False: + vor = Voronoi(pts) + voronoi_plot_2d(vor, all_plots) + voronoi_plot_2d(vor, no_links) + voronoi_plot_2d(vor, voronoi) + a = voronoi.get_ylim() + b = voronoi.get_xlim() + bb = np.array([b[0], b[1], a[0], a[1]]) + G_c, regions, vor = gen_polygons(G_c, bb) + if bb_flag == True: + voronoi_plot_2d(vor, all_plots) + voronoi_plot_2d(vor, no_links) + voronoi_plot_2d(vor, voronoi) + #plot the top-level graph + pts = G_c.vertex_properties["pos"].get_2d_array(range(2)).T + all_plots.scatter(pts[:,0], pts[:, 1], s=20*G_c.vertex_properties["degree"].get_array(), marker="*") + no_links.scatter(pts[:,0], pts[:, 1], s=20*G_c.vertex_properties["degree"].get_array(), marker="*") + voronoi.scatter(pts[:,0], pts[:, 1], s=20*G_c.vertex_properties["degree"].get_array(), marker="*") + #plot the connections of the top level graph + for e in G_c.edges(): + coord = G_c.vertex_properties["pos"][e.source()] + coord2 = G_c.vertex_properties["pos"][e.target()] + x = [coord[0], coord2[0]] + y = [coord[1], coord2[1]] + #all_plots.plot(x, y, 'go--', linewidth=1, markersize=1) + no_links.plot(x, y, 'go--', linewidth=1, markersize=1) + voronoi.plot(x, y, 'go--', linewidth=1, markersize=1) + + #for every subgraph generate a layout and plot the result + for i in range(num_clusters): + g, center = gen_subclusters(G, G_c, i, reposition) + d = G_c.vertex_properties["pos"][i] - center + t = Polygon_mass(g) + #t.distancefield() + for v in g.vertices(): + G.vertex_properties["pos"][g.vertex_properties["idx"][v]] = g.vertex_properties["pos"][v] + d + #g.vertex_properties["pos"][g.vertex_properties["idx"][v]] = g.vertex_properties["pos"][v] + d + #sub_pts = g.vertex_properties["pos"].get_2d_array(range(2)).T + #all_plots.scatter(pts[:,0], pts[:, 1], marker="*") + # for e in g.edges(): + # coord = g.vertex_properties["pos"][e.source()] + # coord2 = g.vertex_properties["pos"][e.target()] + # x = [coord[0], coord2[0]] + # y = [coord[1], coord2[1]] + # plt.plot(x, y, 'ro--', linewidth=1, markersize=1) + + for e in G.edges(): + coord = G.vertex_properties["pos"][e.source()] + coord2 = G.vertex_properties["pos"][e.target()] + x = [coord[0], coord2[0]] + y = [coord[1], coord2[1]] + if (G.vertex_properties["clusters"][e.source()] == G.vertex_properties["clusters"][e.target()]): + all_plots.plot(x, y, 'ro--', linewidth=1, markersize=1) + no_links.plot(x, y, 'ro--', linewidth=1, markersize=1) + else: + all_plots.plot(x, y, 'bo--', linewidth=1, markersize=1) + + no_links.xaxis.set_visible(False) + all_plots.xaxis.set_visible(False) + for v in G_c.vertices(): + region = regions[G_c.vertex_properties["region_idx"][v]] + centroid = centroid_region(region[0]) + G_c.vertex_properties["pos"][v] = centroid + + pts_temp = G_c.vertex_properties["pos"].get_2d_array(range(2)).T + all_plots.scatter(pts_temp[:,0], pts_temp[:, 1], marker='.', color='r') + no_links.scatter(pts_temp[:,0], pts_temp[:, 1], marker='.', color='r') + voronoi.scatter(pts_temp[:,0], pts_temp[:, 1], marker='.', color='r') + + all_plots.set_xlim([bb[0], bb[1]]) + all_plots.set_ylim([bb[2], bb[3]]) + + no_links.set_xlim([bb[0], bb[1]]) + no_links.set_ylim([bb[2], bb[3]]) + + voronoi.set_xlim([bb[0], bb[1]]) + voronoi.set_ylim([bb[2], bb[3]]) + + plt.show() + + + return G, G_c, bb + + + + + + + + +#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) + +G, bbl, bbu = load_nwt("/home/pavel/Documents/Python/GraphGuiQt/network_4.nwt") +G_c, G = gen_clusters(G, bbl, bbu) +num_clusters = 20 + +#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) +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) + +G, G_c, bb = gen_image(G, G_c, "base", reposition = True) +itr = 0 + +#for itr in range(5): +# G, G_c, bb = gen_image(G, G_c, itr, True, bb) +# itr+=1 + +g, center = gen_subclusters(G, G_c) +#d = G_c.vertex_properties["pos"][0] - center +#for v in g.vertices(): +# g.vertex_properties["pos"][v] = g.vertex_properties["pos"][v] + d + +#G_c = nwt.Network.gen_new_fd_layout(G_c) +#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) + + + + -- libgit2 0.21.4