Python中nD线与凸包的交点

Question

我使用 scipy.spatial.ConvexHull 创建了一个凸包。我需要计算凸包和射线之间的交点，从 0 开始，并在某个其他定义点的方向上。已知凸包包含 0，因此应保证相交。问题的维度可以在 2 到 5 之间变化。我尝试了一些谷歌搜索，但没有找到答案。我希望这是计算几何中已知解决方案的常见问题。谢谢。

Answer 1

根据qhull.org，凸包小平面的点x 验证V.x+b=0，其中V 和b 由hull.equations 给出。 （. 在这里代表点积。V 是长度为 1 的法向量。）

如果V是法线，b是偏移量，x是凸面内的点 船体，然后是 Vx+b

如果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

Answer 2

正如 Ante 在 cmets 中提到的，您需要找到船体中所有线/平面/超平面的最近交点。

要找到射线与超平面的交点，请对归一化射线与超平面法线进行点积，这将告诉您沿射线每单位距离在超平面法线方向上移动多远。

如果点积为负，则表示超平面与射线的方向相反，如果为零，则表示射线与其平行且不会相交。

一旦得到正点积，您就可以计算出超平面在射线方向上的距离，方法是将平面在平面法线方向上的距离除以点积。例如，如果平面距离为 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)

2D 测试代码（解决方案推广到更高维度）：

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]