Revisión | c30a7f5c33c935d863797a46d6bef292f1531eb7 (tree) |
---|---|
Tiempo | 2016-04-14 18:58:14 |
Autor | Lorenzo Isella <lorenzo.isella@gmai...> |
Commiter | Lorenzo Isella |
A self-contained example about plotting some point in 3D, calculating their convex hull and visualizing them.
@@ -0,0 +1,61 @@ | ||
1 | +#! /usr/bin/env python | |
2 | + | |
3 | + | |
4 | +import numpy as np | |
5 | +import pylab as pl | |
6 | +import scipy as sp | |
7 | +from scipy.spatial import ConvexHull | |
8 | +from scipy.spatial.distance import euclidean | |
9 | +import matplotlib.pyplot as plt | |
10 | +import mpl_toolkits.mplot3d as a3 | |
11 | + | |
12 | + | |
13 | +# see http://bit.ly/1qJlWkP | |
14 | + | |
15 | +aspect = 0 | |
16 | +while aspect == 0: | |
17 | + | |
18 | + # Generate random points & convex hull | |
19 | + points = np.random.rand(20,3) | |
20 | + hull = ConvexHull(points) | |
21 | + | |
22 | + # Check aspect ratios of surface facets | |
23 | + aspectRatio = [] | |
24 | + for simplex in hull.simplices: | |
25 | + a = euclidean(points[simplex[0],:], points[simplex[1],:]) | |
26 | + b = euclidean(points[simplex[1],:], points[simplex[2],:]) | |
27 | + c = euclidean(points[simplex[2],:], points[simplex[0],:]) | |
28 | + circRad = (a*b*c)/(np.sqrt((a+b+c)*(b+c-a)*(c+a-b)*(a+b-c))) | |
29 | + inRad = 0.5*np.sqrt(((b+c-a)*(c+a-b)*(a+b-c))/(a+b+c)) | |
30 | + aspectRatio.append(inRad/circRad) | |
31 | + | |
32 | + # Threshold for minium allowable aspect raio of surface facets | |
33 | + if np.amin(aspectRatio) > 0.3: | |
34 | + aspect = 1 | |
35 | + | |
36 | + | |
37 | + | |
38 | + | |
39 | +ax = a3.Axes3D(pl.figure()) | |
40 | +facetCol = sp.rand(3) #[0.0, 1.0, 0.0] | |
41 | + | |
42 | +# Plot hull's vertices | |
43 | +#for vert in hull.vertices: | |
44 | +# ax.scatter(points[vert,0], points[vert,1], zs=points[vert,2]) | |
45 | + | |
46 | +# Plot surface traingulation | |
47 | +for simplex in hull.simplices: | |
48 | + vtx = [points[simplex[0],:], points[simplex[1],:], points[simplex[2],:]] | |
49 | + tri = a3.art3d.Poly3DCollection([vtx], linewidths = 2, alpha = 0.8) | |
50 | + tri.set_color(facetCol) | |
51 | + tri.set_edgecolor('k') | |
52 | + ax.add_collection3d(tri) | |
53 | + | |
54 | +# mp = Axes3D.scatter(points[:,0], points[:,1], point[:,2]) | |
55 | + | |
56 | +ax.scatter(points[:,0], points[:,1], points[:,2], s=100) | |
57 | + | |
58 | +plt.axis('off') | |
59 | +plt.show() | |
60 | + | |
61 | +print "So far so good" |