def convex_hull(self):
"""Return a 3D mesh that represents the convex hull of the mesh.
"""
hull = ss.ConvexHull(self.vertices_)
hull_tris = hull.simplices
if self.normals_ is None:
cvh_mesh = Mesh3D(self.vertices_.copy(), hull_tris.copy(), center_of_mass=self.center_of_mass_)
else:
cvh_mesh = Mesh3D(self.vertices_.copy(), hull_tris.copy(), normals=self.normals_.copy(), center_of_mass=self.center_of_mass_)
cvh_mesh.remove_unreferenced_vertices()
return cvh_mesh
python类ConvexHull()的实例源码
def determine_convex_hull(formation_en_grid):
"""
Wraps the scipy.spatial ConvexHull algo for our purposes.
For now only for 2D phase diagrams
Adds the points [1.0, 0.0] and [0.0, 1.0], because in material science these
are always there.
:params: formation_en_grid: list of points in phase space [[x, formation_energy]]
:returns: a hul datatype
"""
import numpy as np
from scipy.spatial import ConvexHull
# TODO multi d
# check if endpoints are in
if [1.0, 0.0] not in formation_en_grid:
formation_en_grid.append([1.0, 0.0])
if [0.0, 0.0] not in formation_en_grid:
formation_en_grid.append([0.0, 0.0])
points = np.array(formation_en_grid)
hull = ConvexHull(points)
return hull
def compute_vc(self, points):
"""compute the Voronoi cell"""
#compute S and V
S = ConvexHull(points).area
V = ConvexHull(points).volume
#voronoi cell
eta = S ** 3 / (36 * np.pi * V ** 2)
return eta
def convex_hull_area(pts):
"""
Calculates the surface area from a given point cloud using simplices of
its convex hull. For the estimation of the synapse contact area, divide by
a factor of two, in order to get the area of only one face (we assume that
the contact site is sufficiently thin represented by the points).
:param pts: np.array of coordinates in nm (scaled)
:return: Area of the point cloud (nm^2)
"""
if len(pts) < 4:
return 0
area = 0
try:
ch = ConvexHull(pts)
triangles = ch.points[ch.simplices]
for triangle in triangles:
area += poly_area(triangle)
except QhullError as e:
# warnings.warn("%s encountered during calculation of convex hull "
# "area with %d points. Returning 0 nm^2." %
# (e, len(pts)), RuntimeWarning)
pass
return area
def fit_convex_hull(points):
""" Creates a feasible set by taking a convex hull of the points given. Returns P = { x : Ax >= b }
Args:
points (list): Set of numpy points.
Returns:
A (numpy): constraint matrix
b (numpy): constraint vector
"""
hull = ConvexHull(points)
m,n = hull.equations.shape
A = -1 * hull.equations[:,0:n-1]
b = hull.equations[:,n-1]
return np.mat(A), np.mat(b).T
def tri_normals(self, align_to_hull=False):
"""Returns a list of the triangle normals.
Parameters
----------
align_to_hull : bool
If true, we re-orient the normals to point outward from
the mesh by using the convex hull.
Returns
-------
:obj:`numpy.ndarray` of float
A #triangles by 3 array of floats, where each 3-ndarray
represents the 3D normal vector of the corresponding triangle.
"""
# compute normals
v0 = self.vertices_[self.triangles_[:,0],:]
v1 = self.vertices_[self.triangles_[:,1],:]
v2 = self.vertices_[self.triangles_[:,2],:]
n = np.cross(v1 - v0, v2 - v0)
normals = n / np.tile(np.linalg.norm(n, axis=1)[:,np.newaxis], [1,3])
# reverse normal based on alignment with convex hull
if align_to_hull:
tri_centers = self.tri_centers()
hull = ss.ConvexHull(tri_centers)
hull_tris = hull.simplices
hull_vertex_ind = hull_tris[0][0]
hull_vertex = tri_centers[hull_vertex_ind]
hull_vertex_normal = normals[hull_vertex_ind]
v = hull_vertex.reshape([1,3])
n = hull_vertex_normal
ip = (tri_centers - np.tile(hull_vertex,
[tri_centers.shape[0], 1])).dot(n)
if ip[0] > 0:
normals = -normals
return normals
def _rubberband(bands, intensities, num_ranges):
'''Basic rubberband method,
from p.77 of "IR and Raman Spectroscopy" (OPUS manual)'''
# create n ranges of equal size in the spectrum
range_size = len(intensities) // num_ranges
y = intensities[:range_size * num_ranges].reshape((num_ranges, range_size))
# find the smallest intensity point in each range
idx = np.arange(num_ranges) * range_size + np.argmin(y, axis=1)
# add in the start and end points as well, to avoid weird edge effects
if idx[0] != 0:
idx = np.append(0, idx)
if idx[-1] != len(intensities) - 1:
idx = np.append(idx, len(intensities) - 1)
baseline_pts = np.column_stack((bands[idx], intensities[idx]))
# wrap a rubber band around the baseline points
hull = ConvexHull(baseline_pts)
hidx = idx[hull.vertices]
# take only the bottom side of the hull
left = np.argmin(bands[hidx])
right = np.argmax(bands[hidx])
mask = np.ones(len(hidx), dtype=bool)
for i in range(len(hidx)):
if i > right and (i < left or right > left):
mask[i] = False
elif i < left and i < right:
mask[i] = False
hidx = hidx[mask]
hidx = hidx[np.argsort(bands[hidx])]
# interpolate a baseline
return np.interp(bands, bands[hidx], intensities[hidx])
def compute_polyhedron_shadow(n, light_poly, ground_poly, inf_dist=1000):
light_hull = ConvexHull(light_poly)
ground_hull = ConvexHull(ground_poly)
light_vertices = [light_poly[i] for i in light_hull.vertices]
ground_vertices = [ground_poly[i] for i in ground_hull.vertices]
rays = [gv - lv for gv in ground_vertices for lv in light_vertices]
return ground_vertices + [v + inf_dist * r
for v in ground_vertices for r in rays]
def polygons(latitudes, longitudes, clusters, maptype=MAPTYPE):
"""Plot clusters of points on map, including them in a polygon defining their convex hull.
:param pandas.Series latitudes: series of sample latitudes
:param pandas.Series longitudes: series of sample longitudes
:param pandas.Series clusters: marker clusters, as integers
:param string maptype: type of maps, see GoogleStaticMapsAPI docs for more info
:return: None
"""
width = SCALE * MAX_SIZE
img, pixels = background_and_pixels(latitudes, longitudes, MAX_SIZE, maptype)
polygons = []
for c in clusters.unique():
in_polygon = clusters == c
if in_polygon.sum() < 3:
print '[WARN] Cannot draw polygon for cluster {} - only {} samples.'.format(c, in_polygon.sum())
continue
cluster_pixels = pixels.loc[clusters == c]
polygons.append(Polygon(cluster_pixels.iloc[ConvexHull(cluster_pixels).vertices], closed=True))
plt.figure(figsize=(10, 10))
ax = plt.subplot(111)
plt.imshow(np.array(img)) # Background map
p = PatchCollection(polygons, cmap='jet', alpha=0.15) # Collection of polygons
p.set_array(clusters.unique())
ax.add_collection(p)
plt.scatter( # Scatter plot
pixels['x_pixel'],
pixels['y_pixel'],
c=clusters,
s=width / 40,
linewidth=0,
alpha=0.25,
)
plt.gca().invert_yaxis() # Origin of map is upper left
plt.axis([0, width, width, 0]) # Remove margin
plt.axis('off')
plt.tight_layout()
plt.show()
def get_edges(a, convex=False):
a = checkma(a)
#Need to deal with RGB images here
#Need to be careful, probably want to take minimum value from masks
if a.ndim == 3:
#Assume that the same mask applies to each band
#Only create edges for one band
b = a[:,:,0]
#Could convert to HSL and do edges for L channel
#b = a[:,:,0].mask + a[:,:,1].mask + a[:,:,2].mask
else:
b = a
#Compute edges along both axes, need both to handle undercuts
#These are inclusive, indices indicate position of unmasked data
edges0 = np.ma.notmasked_edges(b, axis=0)
edges1 = np.ma.notmasked_edges(b, axis=1)
edges = np.array([np.concatenate([edges0[0][0], edges0[1][0], edges1[1][0], edges1[0][0]]), np.concatenate([edges0[0][1], edges0[1][1], edges1[1][1], edges1[0][1]])])
#This is a rough outline - needs testing
if convex:
from scipy.spatial import ConvexHull
#hull = ConvexHull(edges.T)
#edges = edges.T[hull.simplices]
#This is in scipy v0.14
#edges0 = edges1 = hull.vertices
return edges0, edges1, edges
def convex_hull(self, model):
train = np.array(model.data)
global t1,t2
t1,t2 = [],[]
for i,example in enumerate(train):
if example[2]==1:
t1.append(example[0:2])
else:
t2.append(example[0:2])
global ch1,ch2,has_margin
if (has_margin):
self.ax.lines.pop(len(t1)+len(t2)-1)
for i in range(ch1+ch2):
self.ax.lines.pop(len(t1)+len(t2)-1)
if len(t1) > 2:
X1 = np.matrix(np.array(t1))
hull1 = ConvexHull(X1);
for simplex in hull1.simplices:
self.ax.plot(X1[simplex,0],X1[simplex,1], color='grey',linewidth=1)
ch1 = len(hull1.simplices)
if len(t2) > 2:
X2 = np.matrix(np.array(t2))
hull2 = ConvexHull(X2);
for simplex in hull2.simplices:
self.ax.plot(X2[simplex,0],X2[simplex,1], color='0.25')
#self.ax.plot(cx2,cy2,'+',ms=10, color="blue")
ch2 = len(hull2.simplices)
global isFitted
#print len(train),"\t"
#if isFitted:
# self.draw_distance(model,model.clf.gamma)
self.set_dim()
def build_from_hull(self):
"""
Test oriented bounding box algorithm using convex hull points.
Raises:
None
Returns:
EigenVectors(OpenMaya.MVector)
CenterPoint(OpenMaya.MVector)
BoundingExtents(OpenMaya.MVector)
"""
npPointList = [[self.points[i].x, self.points[i].y, self.points[i].z]
for i in xrange(self.points.length())]
try:
hull = ConvexHull(npPointList)
except NameError:
raise RuntimeError(
"From hull method unavailable because"
" scipy cannot be imported."
"Please install it if you need it.")
indices = hull.simplices
vertices = npPointList[indices]
hullPoints = list(vertices.flatten())
hullTriPoints = list(indices.flatten())
hullArray = OpenMaya.MVectorArray()
for ind in xrange(0, len(hullPoints), 3):
hullArray.append(
OpenMaya.MVector(
hullPoints[ind], hullPoints[ind+1], hullPoints[ind+2]))
triPoints = OpenMaya.MIntArray()
for tri in xrange(len(hullTriPoints)):
triPoints.append(tri)
return self.build_from_triangles(points=hullArray, triangles=triPoints)
def generate_one_example(n_nodes):
points = numpy.random.rand(n_nodes, 2)
hull = ConvexHull(points) # scipy.spatial.ConvexHull will generate points in CCW order
v = hull.vertices
# v = numpy.roll(v, -list(v).index(numpy.min(v))) # start from the smallest indice
return points, v + 1
def oriented_bounds_2D(points):
'''
Find an oriented bounding box for a set of 2D points.
Arguments
----------
points: (n,2) float, 2D points
Returns
----------
transform: (3,3) float, homogenous 2D transformation matrix to move the input set of
points to the FIRST QUADRANT, so no value is negative.
rectangle: (2,) float, size of extents once input points are transformed by transform
'''
c = ConvexHull(np.asanyarray(points))
# (n,2,3) line segments
hull = c.points[c.simplices]
# (3,n) points on the hull to check against
dot_test = c.points[c.vertices].reshape((-1,2)).T
edge_vectors = unitize(np.diff(hull, axis=1).reshape((-1,2)))
perp_vectors = np.fliplr(edge_vectors) * [-1.0,1.0]
bounds = np.zeros((len(edge_vectors), 4))
for i, edge, perp in zip(range(len(edge_vectors)),
edge_vectors,
perp_vectors):
x = np.dot(edge, dot_test)
y = np.dot(perp, dot_test)
bounds[i] = [x.min(), y.min(), x.max(), y.max()]
extents = np.diff(bounds.reshape((-1,2,2)), axis=1).reshape((-1,2))
area = np.product(extents, axis=1)
area_min = area.argmin()
offset = -bounds[area_min][0:2]
theta = np.arctan2(*edge_vectors[area_min][::-1])
transform = transformation_2D(offset, theta)
rectangle = extents[area_min]
return transform, rectangle
def convex_hull(mesh, clean=True):
'''
Get a new Trimesh object representing the convex hull of the
current mesh. Requires scipy >.12.
Argments
--------
clean: boolean, if True will fix normals and winding
to be coherent (as qhull/scipy outputs are not)
Returns
--------
convex: Trimesh object of convex hull of current mesh
'''
type_trimesh = type_named(mesh, 'Trimesh')
c = ConvexHull(mesh.vertices.view(np.ndarray).reshape((-1,3)))
vid = np.sort(c.vertices)
mask = np.zeros(len(c.points), dtype=np.int64)
mask[vid] = np.arange(len(vid))
faces = mask[c.simplices]
vertices = c.points[vid].copy()
convex = type_trimesh(vertices = vertices,
faces = faces,
process = True)
if clean:
# the normals and triangle winding returned by scipy/qhull's
# ConvexHull are apparently random, so we need to completely fix them
convex.fix_normals()
return convex
def planar_hull(points, normal, origin=None, input_convex=False):
'''
Find the convex outline of a set of points projected to a plane.
Arguments
-----------
points: (n,3) float, input points
normal: (3) float vector, normal vector of plane
origin: (3) float, location of plane origin
input_convex: bool, if True we assume the input points are already from
a convex hull which provides a speedup.
Returns
-----------
hull_lines: (n,2,2) set of unordered line segments
T: (4,4) float, transformation matrix
'''
if origin is None:
origin = np.zeros(3)
if not input_convex:
pass
planar, T = project_to_plane(points,
plane_normal = normal,
plane_origin = origin,
return_planar = False,
return_transform = True)
hull_edges = ConvexHull(planar[:,0:2]).simplices
hull_lines = planar[hull_edges]
planar_z = planar[:,2]
height = np.array([planar_z.min(),
planar_z.max()])
return hull_lines, T, height
def calculate_boundary(vertices):
# find the points that actually comprise a convex hull
hull = ConvexHull(list(vertices))
# keep only vertices that actually comprise a convex hull and arrange in CCW order
vertices = vertices[hull.vertices]
# get the real number of vertices
nVertices = vertices.shape[0]
# initialize normals array
unit_normals = np.zeros([nVertices, 2])
# determine if point is inside or outside of each face, and distance from each face
for j in range(0, nVertices):
# calculate the unit normal vector of the current face (taking points CCW)
if j < nVertices - 1: # all but the set of point that close the shape
normal = np.array([vertices[j+1, 1]-vertices[j, 1],
-(vertices[j+1, 0]-vertices[j, 0])])
unit_normals[j] = normal/np.linalg.norm(normal)
else: # the set of points that close the shape
normal = np.array([vertices[0, 1]-vertices[j, 1],
-(vertices[0, 0]-vertices[j, 0])])
unit_normals[j] = normal/np.linalg.norm(normal)
return vertices, unit_normals
def get_normals(hull, number_fitnodes=12):
"""
Calculate normals from given hull points using local convex hull fitting.
Orientation of normals is found using local center of mass.
:param hull: 3D coordinates of points representing cell hull
:type hull: np.array
:return: normals for each hull point
"""
normals = np.zeros_like(hull, dtype=np.float)
hull_tree = spatial.cKDTree(hull)
dists, nearest_nodes_ixs = hull_tree.query(hull, k=number_fitnodes,
distance_upper_bound=1000)
for ii, nearest_ixs in enumerate(nearest_nodes_ixs):
nearest_nodes = hull[nearest_ixs[dists[ii] != np.inf]]
ch = ConvexHull(nearest_nodes, qhull_options='QJ Pp')
triangles = ch.points[ch.simplices]
normal = np.zeros((3), dtype=np.float)
# average normal
for triangle in triangles:
cnt = 0
n_help = unit_normal(triangle[0], triangle[1], triangle[2])
if not np.any(np.isnan(n_help)):
normal += np.abs(n_help)
normal /= np.linalg.norm(normal)
normal_sign = (hull[ii] - np.mean(nearest_nodes, axis=0))/\
np.abs(hull[ii] - np.mean(nearest_nodes, axis=0))
normals[ii] = normal * normal_sign
return normals
def convexhull_example(self, length, scales):
points = np.random.uniform(0, 1, [length, self.input_size])
target = -1 * np.ones([length])
ch = ConvexHull(points).vertices
argmin = np.argsort(ch)[0]
ch = list(ch[argmin:]) + list(ch[:argmin])
target[:len(ch)] = np.array(ch)
target += 1
return points, target
def _in_volume_convex(points, volume, remote_instance=None, approximate=False, ignore_axis=[]):
""" Uses scipy to test if points are within a given CATMAID volume.
The idea is to test if adding the point to the cloud would change the
convex hull.
"""
if remote_instance is None:
if 'remote_instance' in sys.modules:
remote_instance = sys.modules['remote_instance']
elif 'remote_instance' in globals():
remote_instance = globals()['remote_instance']
if type(volume) == type(str()):
volume = fetch.get_volume(volume, remote_instance)
verts = volume['vertices']
if not approximate:
intact_hull = ConvexHull(verts)
intact_verts = list(intact_hull.vertices)
if isinstance(points, list):
points = np.array(points)
elif isinstance(points, pd.DataFrame):
points = points.to_matrix()
return [list(ConvexHull(np.append(verts, list([p]), axis=0)).vertices) == intact_verts for p in points]
else:
bbox = [(min([v[0] for v in verts]), max([v[0] for v in verts])),
(min([v[1] for v in verts]), max([v[1] for v in verts])),
(min([v[2] for v in verts]), max([v[2] for v in verts]))
]
for a in ignore_axis:
bbox[a] = (float('-inf'), float('inf'))
return [False not in [bbox[0][0] < p.x < bbox[0][1], bbox[1][0] < p.y < bbox[1][1], bbox[2][0] < p.z < bbox[2][1], ] for p in points]
def get_hull(self):
points = self.get_points_hull()
return ConvexHull(points)
def run_hpr(grid, viewpoint, gamma, boundary):
"""Runs the HPR operator for a single viewpoint and its neighborhood
Runs the HPR operator for a single viewpoint and its neighborhood.
Args:
grid: The support grid containing the data
viewpoint: The viewpoint to be used by the HPR operator.
gamma: The exponent used to flip the points.
boundary: An array of flags used as output.
"""
candidates = grid.get_candidates(viewpoint)
if (candidates.shape[0] == 0):
return
if (candidates.shape[0] <= grid.dimension):
boundary[candidates] = True
return
flipped = exponential_flip(grid.points[candidates], viewpoint, gamma)
# add the viewpoint to the end of the list
flipped = np.vstack([flipped, np.zeros([1, grid.dimension])])
hull = ConvexHull(flipped)
visible_idx = hull.vertices
# remove the index corresponding to the viewpoint
visible_idx.sort()
visible_idx = np.delete(visible_idx, -1)
visible_idx = candidates[visible_idx]
boundary[visible_idx] = True
def compute_polygon_shadow(n, light_poly, ground_poly):
"""Polygons are described by their convex hulls.
n -- normal vector used in ZMP computations
light_poly -- polygon acting as light source
ground_poly -- polygon whose shadow is projected under the light one
"""
t = [n[2] - n[1], n[0] - n[2], n[1] - n[0]]
n = array(n) / pymanoid.toolbox.norm(n)
t = array(t) / pymanoid.toolbox.norm(t)
b = cross(n, t)
light_proj = [array([dot(t, p), dot(b, p)]) for p in light_poly]
light_hull = ConvexHull(light_proj)
ground_proj = [array([dot(t, p), dot(b, p)]) for p in ground_poly]
ground_hull = ConvexHull(ground_proj)
vertex2poly = {i: j for (i, j) in enumerate(ground_hull.vertices)}
light_vertices = [light_proj[i] for i in light_hull.vertices]
ground_vertices = [ground_proj[i] for i in ground_hull.vertices]
mink_diff = [gv - lv for gv in ground_vertices for lv in light_vertices]
try:
u_low, u_high = pymanoid.draw.pick_2d_extreme_rays(mink_diff)
except pymanoid.exceptions.UnboundedPolyhedron:
big_dist = 1000 # [m]
vertices = [
array([-big_dist, -big_dist, light_poly[0][2]]),
array([-big_dist, +big_dist, light_poly[0][2]]),
array([+big_dist, -big_dist, light_poly[0][2]]),
array([+big_dist, +big_dist, light_poly[0][2]])]
return vertices, []
nb_vertices = len(ground_vertices)
vertex_indices = range(len(ground_vertices))
def f_low(i):
return cross(u_low, ground_vertices[i])
def f_high(i):
return cross(u_high, ground_vertices[i])
v_low = min(vertex_indices, key=f_low)
v_high = max(vertex_indices, key=f_high)
vertices = [ground_poly[vertex2poly[vertex_index % nb_vertices]]
for vertex_index in xrange(v_high, (v_low + nb_vertices + 1))]
rays = [u_low[0] * t + u_low[1] * b,
u_high[0] * t + u_high[1] * b]
return vertices, rays
def sortCorners(corners):
'''
sort the corners of a given quadrilateral of the type
corners : [ [xi,yi],... ]
to an anti-clockwise order starting with the bottom left corner
or (if plotted as image where y increases to the bottom):
clockwise, starting top left
'''
corners = np.asarray(corners)
# bring edges in order:
corners2 = corners[ConvexHull(corners).vertices]
if len(corners2) == 3:
# sometimes ConvexHull one point is missing because it is
# within the hull triangle
# find the right position of set corner as the minimum perimeter
# built with that point as different indices
for c in corners:
if c not in corners2:
break
perimeter = []
for n in range(0, 4):
corners3 = np.insert(corners2, n, c, axis=0)
perimeter.append(
np.linalg.norm(
np.diff(
corners3,
axis=0),
axis=1).sum())
n = np.argmin(perimeter)
corners2 = np.insert(corners2, n, c, axis=0)
# find the edge with the right angle to the quad middle:
mn = corners2.mean(axis=0)
d = (corners2 - mn)
ascent = np.arctan2(d[:, 1], d[:, 0])
bl = np.abs(BL_ANGLE + ascent).argmin()
# build a index list starting with bl:
i = list(range(bl, 4))
i.extend(list(range(0, bl)))
return corners2[i]
def lambda_cone(accuracy, runtime, ax, c, conesize, lines=1, aspect_equal=1):
#ax.scatter(runtime, accuracy, c=tradeoff, lw=0)
P = zip(runtime, accuracy)
#P.extend([(0,0), (runtime.max()+0.5e7, accuracy.max())])
P.extend([(runtime.min()-.1*runtime.ptp(), 0),
(runtime.max()+.1*runtime.ptp(), accuracy.max())])
# tradeoff = np.array(list(tradeoff) + [np.inf, 0.0])
hull = ConvexHull(P)
v = hull.points[hull.vertices]
ddd = []
V = len(hull.vertices)
for i in range(V):
x1,y1 = v[(i-1) % V]
x,y = v[i]
x3,y3 = v[(i+1) % V]
# slopes give a range for lambda values to try.
m12 = (y-y1)/(x-x1)
m23 = (y3-y)/(x3-x)
# Note: This filter skips cases where the convex hull contains a
# segments (edge) that is beneath the fronter (we know this because of
# the orientation of the edge)..
if m12 >= m23:
continue
ddd.append([m12, m23, x, y])
ax.add_artist(Polygon([[x,y],
point(x, y, m12, -conesize),
point(x, y, m23, -conesize),
[x,y]],
linewidth=0,
alpha=0.5,
color=c,
))
if lines:
ax.plot([x1,x,x3], [y1,y,y3], c=c, alpha=0.5)
if aspect_equal:
ax.set_aspect('equal')
ax.grid(True)
return ddd
