segmentation / graph_gui

  #!/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<double>", val=color)
      color = [1.0, 0.0, 0.0, 0.9]
      G.vertex_properties["RGBA"] = G.new_vertex_property("vector<double>", 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<double>", val=np.zeros((3,1), dtype=np.float32))
      G_cluster.vertex_properties["RGBA"] = G_cluster.new_vertex_property("vector<double>", 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)