Intersection of nD line with convex hull in Python

Question

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.

Answer 1

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 α:

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

Answer 2

As mentioned by Ante in the comments, you need to find the closest intersection of all the lines/planes/hyper-planes in the hull.

To find the intersection of the ray with the hyperplane, do a dot product of the normalized ray with the hyperplane normal, which will tell you how far in the direction of the hyperplane normal you move for each unit distance along the ray.

If the dot product is negative it means that the hyperplane is in the opposite direction of the ray, if zero it means the ray is parallel to it and won't intersect.

Once you have a positive dot product, you can work out how far away the hyperplane is in the direction of the ray, by dividing the distance of the plane in the direction of the plane normal by the dot product. For example if the plane is 3 units away, and the dot product is 0.5, then you only get 0.5 units closer for every unit you move along the ray, so the hyperplane is 3 / 0.5 = 6 units away in the direction of the ray.

Once you have calculated this distance for all the hyperplanes and found the closest one, the intersection point is just the ray multiplied by the closest distance.

Here is a solution in Python (normalize function is from 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)

Test code in 2D (the solution generalizes to higher dimensions):

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

output:

[ 0.66666667 -0.66666667]