Search code examples
pythoncomputational-geometryintersectionconvex-hull

Intersection of nD line with convex hull in Python


I have created a convex hull using scipy.spatial.ConvexHull. I need to compute the intersection point between the convex hull and a ray, starting at 0 and in the direction of some other defined point. The convex hull is known to contain 0 so the intersection should be guaranteed. The dimension of the problem can vary between 2 and 5. I have tried some google searching but haven't found an answer. I am hoping this is a common problem with known solutions in computational geometry. Thank you.


Solution

  • According to qhull.org, the points x of a facet of the convex hull verify V.x+b=0, where V and b are given by hull.equations. (. stands for the dot product here. V is a normal vector of length one.)

    If V is a normal, b is an offset, and x is a point inside the convex hull, then Vx+b <0.

    If U is a vector of the ray starting in O, the equation of the ray is x=αU, α>0. so the intersection of ray an facet is x = αU = -b/(V.U) U. The unique intersection point with the hull corresponds to the min of the positive values of α: enter image description here

    The next code give it :

    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
    

    It is a pure numpy solution so it is fast. An example for 1 million points in the [-1,1]^3 cube :

    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