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.spatial import Delaunay, delaunay_plot_2d

from triangle import triangulate
from triangle import plot as triangle_plot
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
        print(nwt.gt.graph_tool.topology.is_planar(G))
        self.get_aabb()
        self.gen_polygon()
        self.torque = []
        self.forces_r = np.zeros(2)
        self.forces_a = np.zeros(2)
        self.vel = 0.0
        self.aa = 0.0
        self.degree = 0
      
        
    def clear_torques(self):
        self.torque = []        
        
    def clear_forces(self):
        self.forces_r = np.zeros(2)
        self.forces_a = np.zeros(2)
        
    def set_degree(self, degree):
        self.degree = degree
        
    def add_torque(self, p, f):
        #direction of the torque = cross of (r, f)
        #magnitude = ||r||*||f||*sin(theta)
        #r = level arm vector
        d = self.CoM - p
        #r = np.linalg.norm(self.CoM - p)
        value = np.dot(d, f)/np.dot(d, d)/np.dot(f, f)
        if value < 1.0 and value > -1.0:
            theta = math.acos(value)
            torque = math.sin(theta) * np.sqrt(f[0]*f[0]+f[1]*f[1]) * np.sqrt(d[0]*d[0]+d[1]*d[1])
        else:
            #print("value = ", value)
            torque = 0.0
        
        #if < 0 then clockwise, else counter
        direction = np.cross(d, f)
        if direction < 0:
            self.torque.append([torque, f, p, "counterclock"])
        else:
            self.torque.append([torque, f, p, "clockwise"])
        
            
    def calculate_moment(self, use_graph=False):
        
        #returns the area of a triangle defined by two points
        def area(t):
            output = 1.0/2.0*abs(t[0,0]*(t[1,1]-t[2,1]) + t[1,0]*(t[2,1]-t[0,1]) + t[2,0]*(t[0,1]-t[1,1]))
            return output
        
        def center(t):
            output = np.asarray([(t[0,0]+t[1,0]+t[2,0])/3.0, (t[0,1]+t[1,1]+t[2,1])/3.0])
            return output
        
        segs = []
        pts = np.asarray(self.polygon.exterior.xy).T
        pts = pts[:-1, :]
        for k in range(pts.shape[0]-1):
            segs.append([k, k+1])
        segs.append([pts.shape[0]-1, 0])
        if use_graph:
            points = self.G.vertex_properties["pos"].get_2d_array(range(2)).T
            segs2 = []
            n_pts = pts.shape[0]
            for e in self.G.edges():
                segs2.append([int(e.source())+n_pts, int(e.target())+n_pts])
            pts = np.concatenate((pts, points))
            segs = segs + segs2
        mesh = dict(vertices=pts, segments=segs)
        #print(self.polygon.area())
        tri = triangulate(mesh, 'pq20Ds')
        self.mesh = mesh
        self.tri = tri
        moment = 0.0
        #NEED TO ADD MASS maybe?
        for i in range(tri['triangles'].shape[0]):
            t = tri['vertices'][tri['triangles'][i]]
            A = area(t)
            C = center(t)
            Moi = A+A*np.linalg.norm(C-self.CoM)**2
            moment += Moi
            
        self.MoI = abs(math.log(moment))
        #self.MoI = 10
        print(self.MoI)
#        triangle_plot(plt.gca(), **tri)
#        plt.gca().set_title(str(self.polygon.area))
     
    def translate(self, step):
        d = self.forces_a + self.forces_r
        #print(self.forces_a, self.forces_r)
        d0 = step*d
        self.CoM = self.CoM + d0
        pos = self.G.vertex_properties["pos"].get_2d_array(range(2)).T
        pos = pos + d0
        self.G.vertex_properties["pos"] = self.G.new_vertex_property("vector<double>", vals = pos)
    
    def rotate(self, phi, direction = "counterclock"):
        if("counterclock"):
            for v in self.G.vertices():
                p_prime = copy.deepcopy(self.G.vertex_properties["pos"][v])
                p = copy.deepcopy(self.G.vertex_properties["pos"][v])
                p_prime[0] = self.CoM[0] + math.cos(phi) * (p[0] - self.CoM[0]) - math.sin(phi) * (p[1] - self.CoM[1])
                p_prime[1] = self.CoM[1] + math.sin(phi) * (p[0] - self.CoM[0]) + math.cos(phi) * (p[1] - self.CoM[1])
                self.G.vertex_properties["pos"][v] = p_prime
            #rotate points in mesh
            points = copy.deepcopy(self.mesh['vertices'])
            for v in range(points.shape[0]):
                points[v][0] = self.CoM[0] + math.cos(phi) * (self.mesh['vertices'][v][0] - self.CoM[0]) - math.sin(phi) * (self.mesh['vertices'][v][1] - self.CoM[1])
                points[v][1] = self.CoM[1] + math.sin(phi) * (self.mesh['vertices'][v][0] - self.CoM[0]) + math.cos(phi) * (self.mesh['vertices'][v][1] - self.CoM[1])
            self.mesh['vertices'] = points
        else:
            for v in self.G.vertices():
                p_prime = copy.deepcopy(self.G.vertex_properties["pos"][v])
                p = copy.deepcopy(self.G.vertex_properties["pos"][v])
                p_prime[0] = self.CoM[0] + math.cos(phi) * (p[0] - self.CoM[0]) + math.sin(phi) * (p[1] - self.CoM[1])
                p_prime[1] = self.CoM[1] - math.sin(phi) * (p[0] - self.CoM[0]) + math.cos(phi) * (p[1] - self.CoM[1])
                self.G.vertex_properties["pos"][v] = p_prime
            #rotate points in mesh
            points = copy.deepcopy(self.mesh['vertices'])
            for v in range(points.shape[0]):
                points[v][0] = self.CoM[0] + math.cos(phi) * (self.mesh['vertices'][v][0] - self.CoM[0]) + math.sin(phi) * (self.mesh['vertices'][v][1] - self.CoM[1])
                points[v][1] = self.CoM[1] - math.sin(phi) * (self.mesh['vertices'][v][0] - self.CoM[0]) + math.cos(phi) * (self.mesh['vertices'][v][1] - self.CoM[1])
            self.mesh['vertices'] = points
                
        
    
     
    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='*')
        
        #mesh = dict(vertices=pts, segments=segs)
        #print(self.polygon.area())
        #tri = triangulate(mesh, 'pq20Ds')
        triangle_plot(plt.gca(), **self.mesh)
        
        #plot polygon
        #plt.plot(*self.polygon.exterior.xy, color = 'r')
        
        
#        for i in range(len(segs)):
#            plt.plot((pts[segs[i][0]][0], pts[segs[i][1]][0]), (pts[segs[i][0]][1], pts[segs[i][1]][1]), color='b')
        plt.gca().set_title(str(self.polygon.area))
        
#        for e in self.torque:
#            plt.quiver(e[2][0], e[2][1], e[1][0], e[1][1], color='r')
        #tri = Delaunay(np.asarray(self.polygon.exterior.coords.xy).T)
        #tri = triangulate(mesh, 'pq20a' + str(self.polygon.area/100.0)+'D')
        #delaunay_plot_2d(tri)
        
#        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.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)
        self.calculate_moment(True)
        #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)
        self.D = D
        self.x = x
        self.y = y
        
        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])  
    
#Graph G
#List of Polygonmass objects with the same cluster index.
def get_torques(G, masses):
    for i in masses:
        i.clear_torques()
    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()]):
            #index of the cluster
            t0 = G.vertex_properties["clusters"][e.target()]
            t1 = G.vertex_properties["clusters"][e.source()]
            #index of the vertex outside of the subgraph
            v0_index = G.vertex_properties["idx"][e.target()]
            v1_index = G.vertex_properties["idx"][e.source()]
            #location of torque arm in the subgraph
            p0 = masses[t0].G.vertex_properties["pos"][np.argwhere(masses[t0].G.vertex_properties["idx"].get_array() == v0_index)]
            p1 = masses[t1].G.vertex_properties["pos"][np.argwhere(masses[t1].G.vertex_properties["idx"].get_array() == v1_index)]
            
            f0 = np.subtract(p0, p1)
            f1 = np.subtract(p1, p0)
            masses[t0].add_torque(p0, f1)
            masses[t1].add_torque(p1, f0)
'''
    c1 scales the attractive force before log.
    c2 scales the attractive force inside log.
    c3 scales the repulsive force.
'''
def get_forces(G, masses, c1 = 50.0, c2 = 1.0, c3 = 1.0):
    for i in range(len(masses)):
        masses[i].clear_forces()
        f_total = np.zeros(2)
        for j in range(len(masses)):
            if i != j:
                f0 = np.subtract(masses[i].CoM, masses[j].CoM)
                f1 = np.power(f0, 3.0)
                f1[0] = c3/f1[0]*np.sign(f0[0])
                f1[1] = c3/f1[1]*np.sign(f0[1])
                f_total = np.add(f_total, f1)
        masses[i].forces_r = f_total
            
    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()]):
            #index of the cluster
            t0 = G.vertex_properties["clusters"][e.target()]
            t1 = G.vertex_properties["clusters"][e.source()]
            #index of the vertex outside of the subgraph
            v0_index = G.vertex_properties["idx"][e.target()]
            v1_index = G.vertex_properties["idx"][e.source()]
            #location of torque arm in the subgraph
            p0 = masses[t0].G.vertex_properties["pos"][np.argwhere(masses[t0].G.vertex_properties["idx"].get_array() == v0_index)]
            p1 = masses[t1].G.vertex_properties["pos"][np.argwhere(masses[t1].G.vertex_properties["idx"].get_array() == v1_index)]
            
            f0 = np.subtract(p1, p0)
            f0_1 = abs(f0)
            f0_1 = c1*np.log(f0_1/c2)/masses[t0].degree
            f0_1 = f0_1*np.sign(f0)
            f1 = np.subtract(p0, p1)
            f1_1 = abs(f1)
            f1_1 = c1*np.log(f1_1/c2)/masses[t1].degree
            f1_1 = f1_1*np.sign(f1)
            masses[t0].forces_a = np.add(masses[t0].forces_a, f1_1)
            masses[t1].forces_a = np.add(masses[t1].forces_a, f0_1)
            
    
    
        

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)
    G_cluster.edge_properties["log"] = G_cluster.new_edge_property("double", vals=abs(np.log(G_cluster.edge_properties["volume"].get_array())))    
    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 find_equlibrium(masses, t = 0.01):
    for m in masses:
        sum_torque = 0
        for torque in m.torque:
            if torque[3] == "clockwise":
                sum_torque -= torque[0]
            else:
                sum_torque += torque[0]
        m.vel = m.aa * t + m.vel
        m.aa = sum_torque/m.MoI
        #print(m.G.vertex_properties["clusters"][0], m.vel)
        if m.vel != 0.0:
            if m.vel < 0.0:
                m.rotate(abs(m.vel * t),"counterclock")
            else:
                m.rotate(abs(m.vel * t), "clockwise")


def gen_Eades(G, masses, M = 10):
    for i in range(M):
        get_forces(G, masses)
        for j in masses:
            j.translate(0.001)
        
#def onion_springs(G, masses, min_length):
#    for v in G.vertices():
        

    
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(4, 1, sharex='col', sharey='row')
    fig.tight_layout()
    grid = plt.GridSpec(4,1)
    grid.update(wspace=0.025, hspace=0.2)
    ax[0].axis('off')
    ax[1].axis('off')
    ax[2].axis('off')
    ax[3].axis('off')
    
    
    #Add plots to the axes and get their handles
    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)
    rotated = fig.add_subplot(grid[3], sharey=all_plots, sharex=all_plots)
    
    #Get the points and generate the voronoi region
    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]])
        
    #generate the polygons based on the voronoi regions
    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 resulting Polygon_mass
    #These polygons are not the same as the voronoi polygons and instead surround
    #graph.
    masses = []
    for i in range(num_clusters):
        g, center = gen_subclusters(G, G_c, i, reposition)
        d = G_c.vertex_properties["pos"][i] - center
        for v in g.vertices():
            G.vertex_properties["pos"][g.vertex_properties["idx"][v]] = g.vertex_properties["pos"][v] + d
            g.vertex_properties["pos"][v] = g.vertex_properties["pos"][v] + d
        t = Polygon_mass(g)
        t.set_degree(G_c.vertex_properties["degree"][i])
        masses.append(t)
        #t.distancefield()
            
    #get the torques generated by the positioning of the graphs.
    get_torques(G, masses)
#    for i in masses:
#        i.plot_graph(i.D, i.x, i.y)
        #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)
    
    
    #Plot the cluster level connections and the vertex level connections.
    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)
    
    #Update the centroids based on the voronoi polygons
    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
    
    
    print(G_c.num_vertices(), G_c.num_edges())
    
    #Plots the vertices of the subgraphs
    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')
    
    #set the limits of the plots.
    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]])
    
    #show the plots.
    plt.show()

#    for j in [7, 10, 15]:
#        masses[j].plot_graph(masses[j].D, masses[j].x, masses[j].y)

    for j in range(10):
        for i in range(100):
            get_torques(G, masses)
            find_equlibrium(masses)
        gen_Eades(G, masses)
#        for j in [7, 10]:
#            masses[j].plot_graph(masses[j].D, masses[j].x, masses[j].y)
#    for j in [7, 10, 15]:
#        masses[j].plot_graph(masses[j].D, masses[j].x, masses[j].y)
        
    for i in masses:
        for v in i.G.vertices():
            G.vertex_properties["pos"][i.G.vertex_properties["idx"][v]] = i.G.vertex_properties["pos"][v]

    
    #Plot the cluster level connections and the vertex level connections.
    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()]):
            rotated.plot(x, y, 'ro--', linewidth=1, markersize=1)
        else:
            rotated.plot(x, y, 'bo--', linewidth=1, markersize=1)
    
    pts_temp = G.vertex_properties["pos"].get_2d_array(range(2)).T
    rotated.scatter(pts_temp[:,0], pts_temp[:, 1], marker='.', color='r')
    
    

    return G, G_c, bb, masses

            


        
        
        
            
#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, masses = gen_image(G, G_c, "base", reposition = True)

print("Planarity test: G_c, G = " , nwt.gt.graph_tool.topology.is_planar(G_c), nwt.gt.graph_tool.topology.is_planar(G))

itr = 0

#for itr in range(5):
#    G, G_c, bb, masses = 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)