diff --git a/voronoi_test.py b/voronoi_test.py index c32242f..34f5f65 100644 --- a/voronoi_test.py +++ b/voronoi_test.py @@ -7,6 +7,10 @@ Created on Tue Sep 3 13:15:54 2019 """ 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 @@ -30,13 +34,59 @@ class Polygon_mass: self.get_aabb() self.gen_polygon() self.torque = [] + self.vel = 0.0 + self.aa = 0.0 + + + def clear_torques(self): + self.torque = [] 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) 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) + torque = math.sin(theta) * np.sqrt(f[0]*f[0]+f[1]*f[1]) * np.sqrt(r[0]*r[0]+r[1]*r[1]) + #if < 0 then clockwise, else counter + direction = np.cross(r, f) + if direction < 0: + self.torque.append([torque, "clockwise"]) + else: + self.torque.append([torque, "counterclock"]) + + def calculate_moment(self): + + #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]) + + 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]) + mesh = dict(vertices=pts, segments=segs) + #print(self.polygon.area()) + tri = triangulate(mesh, 'pq20D') + 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 = moment +# triangle_plot(plt.gca(), **tri) +# plt.gca().set_title(str(self.polygon.area)) def rotate(self, phi, direction = "counterclock"): if("counterclock"): @@ -62,6 +112,35 @@ class Polygon_mass: plt.scatter(p[:,0], p[:,1], color='r') plt.scatter(self.CoM[0], self.CoM[1], marker='*') + + 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]) + 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, 'pq20D') + triangle_plot(plt.gca(), **mesh) + for i in range(len(segs)): + plt.plot((pts[segs[i][0]][0], pts[segs[i][0]][1]), (pts[segs[i][1]][0], pts[segs[i][1]][1]), color='b') + plt.gca().set_title(str(self.polygon.area)) + + #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) @@ -73,14 +152,14 @@ class Polygon_mass: # 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) - +# 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() @@ -211,8 +290,33 @@ class Polygon_mass: 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]) - + 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()]): + pf0 = G.vertex_properties["pos"][e.target()].get_2D_array(range(2)) + pf1 = G.vertex_properties["pos"][e.source()].get_2D_array(range(2)) + + t0 = G.vertex_properties["clusters"][e.source()] + t1 = G.vertex_properties["clusters"][e.target()] + + v0_index = G.vertex_properties["idx"][e.source()] + v1_index = G.vertex_properties["idx"][e.target()] + + p0 = masses[t0].G.vertex_properties["pos"][masses[t0].G.vertex(v0_index)].get_2D_array(range(2)) + p1 = masses[t1].G.vertex_properties["pos"][masses[t1].G.vertex(v1_index)].get_2D_array(range(2)) + + f0 = pf0-p1 + f1 = pf1-p0 + masses[t0].add_torque(p0, f0) + masses[t1].add_torque(p1, f1) def voronoi_polygons(voronoi, diameter): @@ -517,15 +621,19 @@ def gen_image(G, G_c, itr, bb_flag = False, bb = None, reposition = False): 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 every subgraph generate a layout and plot the resulting Polygon_mass + masses = [] 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) + masses.append(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 + + + #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(): -- libgit2 0.21.4