nD 线与 Python 中凸包的交点
Intersection of nD line with convex hull in Python
我使用 scipy.spatial.ConvexHull 创建了一个凸包。我需要计算凸包和射线之间的交点,从 0 开始并沿其他定义点的方向。已知凸包包含 0,因此应保证相交。问题的维度可以在 2 到 5 之间变化。我尝试了一些 google 搜索但没有找到答案。我希望这是计算几何中已知解决方案的常见问题。谢谢。
正如Ante在评论中提到的,您需要找到船体中所有lines/planes/hyper-planes最近的交集。
要找到射线与超平面的交点,请对归一化射线与超平面法线进行点积,这将告诉您沿着超平面法线的每个单位距离移动多远雷.
如果点积为负则表示超平面与光线方向相反,如果为零则表示光线与其平行且不相交。
一旦得到正点积,就可以计算出超平面在光线方向上的距离,方法是将平面在平面法线方向上的距离除以点积。例如,如果平面距离为 3 个单位,并且点积为 0.5,那么沿着射线每移动一个单位,您只能靠近 0.5 个单位,因此超平面在射线方向上距离为 3 / 0.5 = 6 个单位.
一旦你为所有超平面计算了这个距离并找到了最近的一个,交点就是射线乘以最近的距离。
下面是Python中的一个解决方案(归一化函数来自here):
def normalize(v):
norm = np.linalg.norm(v)
if norm == 0:
return v
return v / norm
def find_hull_intersection(hull, ray_point):
# normalise ray_point
unit_ray = normalize(ray_point)
# find the closest line/plane/hyperplane in the hull:
closest_plane = None
closest_plane_distance = 0
for plane in hull.equations:
normal = plane[:-1]
distance = plane[-1]
# if plane passes through the origin then return the origin
if distance == 0:
return np.multiply(ray_point, 0) # return n-dimensional zero vector
# if distance is negative then flip the sign of both the
# normal and the distance:
if distance < 0:
np.multiply(normal, -1);
distance = distance * -1
# find out how much we move along the plane normal for
# every unit distance along the ray normal:
dot_product = np.dot(normal, unit_ray)
# check the dot product is positive, if not then the
# plane is in the opposite direction to the rayL
if dot_product > 0:
# calculate the distance of the plane
# along the ray normal:
ray_distance = distance / dot_product
# is this the closest so far:
if closest_plane is None or ray_distance < closest_plane_distance:
closest_plane = plane
closest_plane_distance = ray_distance
# was there no valid plane? (should never happen):
if closest_plane is None:
return None
# return the point along the unit_ray of the closest plane,
# which will be the intersection point
return np.multiply(unit_ray, closest_plane_distance)
二维测试代码(解决方案推广到更高维度):
from scipy.spatial import ConvexHull
import numpy as np
points = np.array([[-2, -2], [2, 0], [-1, 2]])
h = ConvexHull(points)
closest_point = find_hull_intersection(h, [1, -1])
print closest_point
输出:
[ 0.66666667 -0.66666667]
根据qhull.org,凸包的一个面的点x
验证V.x+b=0
,其中V
和b
由hull.equations
。 (这里.
代表点积,V
是长度为1的法向量。)
If V is a normal, b is an offset, and x is a point inside the convex
hull, then Vx+b <0.
如果U
是从O
开始的射线向量,则射线方程为x=αU, α>0
。所以射线的交点是 x = αU = -b/(V.U) U
。与船体的唯一交点对应于 α
正值的最小值:
下一个代码给它:
import numpy as np
from scipy.spatial import ConvexHull
def hit(U,hull):
eq=hull.equations.T
V,b=eq[:-1],eq[-1]
alpha=-b/np.dot(V,U)
return np.min(alpha[alpha>0])*U
这是一个纯粹的 numpy 解决方案,所以速度很快。 [-1,1]^3
立方体中 100 万个点的示例:
In [13]: points=2*np.random.rand(1e6,3)-1;hull=ConvexHull(points)
In [14]: %timeit x=hit(np.ones(3),hull)
#array([ 0.98388702, 0.98388702, 0.98388702])
10000 loops, best of 3: 30 µs per loop
我使用 scipy.spatial.ConvexHull 创建了一个凸包。我需要计算凸包和射线之间的交点,从 0 开始并沿其他定义点的方向。已知凸包包含 0,因此应保证相交。问题的维度可以在 2 到 5 之间变化。我尝试了一些 google 搜索但没有找到答案。我希望这是计算几何中已知解决方案的常见问题。谢谢。
正如Ante在评论中提到的,您需要找到船体中所有lines/planes/hyper-planes最近的交集。
要找到射线与超平面的交点,请对归一化射线与超平面法线进行点积,这将告诉您沿着超平面法线的每个单位距离移动多远雷.
如果点积为负则表示超平面与光线方向相反,如果为零则表示光线与其平行且不相交。
一旦得到正点积,就可以计算出超平面在光线方向上的距离,方法是将平面在平面法线方向上的距离除以点积。例如,如果平面距离为 3 个单位,并且点积为 0.5,那么沿着射线每移动一个单位,您只能靠近 0.5 个单位,因此超平面在射线方向上距离为 3 / 0.5 = 6 个单位.
一旦你为所有超平面计算了这个距离并找到了最近的一个,交点就是射线乘以最近的距离。
下面是Python中的一个解决方案(归一化函数来自here):
def normalize(v):
norm = np.linalg.norm(v)
if norm == 0:
return v
return v / norm
def find_hull_intersection(hull, ray_point):
# normalise ray_point
unit_ray = normalize(ray_point)
# find the closest line/plane/hyperplane in the hull:
closest_plane = None
closest_plane_distance = 0
for plane in hull.equations:
normal = plane[:-1]
distance = plane[-1]
# if plane passes through the origin then return the origin
if distance == 0:
return np.multiply(ray_point, 0) # return n-dimensional zero vector
# if distance is negative then flip the sign of both the
# normal and the distance:
if distance < 0:
np.multiply(normal, -1);
distance = distance * -1
# find out how much we move along the plane normal for
# every unit distance along the ray normal:
dot_product = np.dot(normal, unit_ray)
# check the dot product is positive, if not then the
# plane is in the opposite direction to the rayL
if dot_product > 0:
# calculate the distance of the plane
# along the ray normal:
ray_distance = distance / dot_product
# is this the closest so far:
if closest_plane is None or ray_distance < closest_plane_distance:
closest_plane = plane
closest_plane_distance = ray_distance
# was there no valid plane? (should never happen):
if closest_plane is None:
return None
# return the point along the unit_ray of the closest plane,
# which will be the intersection point
return np.multiply(unit_ray, closest_plane_distance)
二维测试代码(解决方案推广到更高维度):
from scipy.spatial import ConvexHull
import numpy as np
points = np.array([[-2, -2], [2, 0], [-1, 2]])
h = ConvexHull(points)
closest_point = find_hull_intersection(h, [1, -1])
print closest_point
输出:
[ 0.66666667 -0.66666667]
根据qhull.org,凸包的一个面的点x
验证V.x+b=0
,其中V
和b
由hull.equations
。 (这里.
代表点积,V
是长度为1的法向量。)
If V is a normal, b is an offset, and x is a point inside the convex hull, then Vx+b <0.
如果U
是从O
开始的射线向量,则射线方程为x=αU, α>0
。所以射线的交点是 x = αU = -b/(V.U) U
。与船体的唯一交点对应于 α
正值的最小值:
下一个代码给它:
import numpy as np
from scipy.spatial import ConvexHull
def hit(U,hull):
eq=hull.equations.T
V,b=eq[:-1],eq[-1]
alpha=-b/np.dot(V,U)
return np.min(alpha[alpha>0])*U
这是一个纯粹的 numpy 解决方案,所以速度很快。 [-1,1]^3
立方体中 100 万个点的示例:
In [13]: points=2*np.random.rand(1e6,3)-1;hull=ConvexHull(points)
In [14]: %timeit x=hit(np.ones(3),hull)
#array([ 0.98388702, 0.98388702, 0.98388702])
10000 loops, best of 3: 30 µs per loop