Search code examples
ralgorithmintersectionnumerical-methodsconvex-hull

Algorithm intersecting polytope and half line


I am looking for an algorithm in R to intersect a convex polytope with a line segment. I found several post here on stack exchange for in the plane but I am wondering if this algorithms exists in higher dimensions. My google searches didn't really produce a lot of answers.

The line segment is composed of a point inside and a point outside the convex polytope. Are there algorithms in R available that can do this in dimension N<=10 ? Or does anyone know a reference so I can implement the algorithm myself? Is there information on the complexity of finding the polytope and the intersection ?


Solution

  • For problems in computational geometry, dimension d > 3 usually might as well be d arbitrary. If you have the polytope as a collection of intersected halfspaces, then likely the only sensible thing to do is to intersect the line segment with each of the separating hyperplanes (by solving a system of d linear equations) and take the intersection closest to the point inside.

    If you have only the vertices of the polytope or even just a set of vertices whose convex closure is the polytope, then the easiest approach given R's libraries probably is linear programming. (Conceivably you could compute the facets using an algorithm to find high-dimensional convex hulls, but there could be Theta(n^floor(d/2)) of them, where n is the number of vertices.) I'm not familiar with LP solvers in R, so I'll write down the program mathematically. It shouldn't be too hard to translate. Let p_0 be the point outside and p_1 be the point inside and v_i be the ith point defining the polytope.

    maximize alpha_0
    subject to
    
    for 1 <= j <= d,
        p_0[j] alpha_0 + p_1[j] alpha_1 - sum_{1 <= i <= n} v_i[j] beta_i = 0
    alpha_0 + alpha_1 = 1
    sum_{1 <= i <= n} beta_i = 1
    
    alpha_0 >= 0
    alpha_1 >= 0
    for 1 <= i <= n,
        beta_i >= 0
    

    The intersection is defined by the point p_0 alpha_0 + p_1 alpha_1 (unless the program is infeasible, in which case there is no intersection).