How to fastest check if point (3D) is inside convex hull given by set of point

Question

I have set P of point (3D) which are vertices of convex hull (every one). I looking for method for checking if given point p0 is NOT outside of this convex hull.

I'll have to repeat the checking it several times (for different p0). So if is possible to reuse part of computation it would be great.

On stackoverflow pages i found this: Find if a point is inside a convex hull for a set of points without computing the hull itself There are 2 aproche: First based on convex hull property - set of linear equations. Secound based on observation: "The point lies outside of the convex hull of the other points if and only if the direction of all the vectors from it to those other points are on less than one half of a circle/sphere/hypersphere around it."

Unfortunately I don't know how exactly can do it. First give me insoluble system of equations - 3 equation with n unknown (n>3). How can I sovle it? Did I do some mistake? In second approach I don't know how check this assumption.

Answer 1

CGAL can easily do this test. Not only you could build the convex hull with CGAL, but you can quickly and easily determine whether a specific point is inside, on the surface or outside of the polyhedron.

A code snippet is shown below:

#include <CGAL/Simple_cartesian.h>
#include <CGAL/AABB_tree.h>
#include <CGAL/AABB_traits.h>
#include <CGAL/Polyhedron_3.h>
#include <CGAL/boost/graph/graph_traits_Polyhedron_3.h>
#include <CGAL/AABB_face_graph_triangle_primitive.h>
#include <CGAL/algorithm.h>
#include <CGAL/Side_of_triangle_mesh.h>

typedef CGAL::Simple_cartesian<double> K;
typedef K::Point_3 Point;
typedef CGAL::Polyhedron_3<K> Polyhedron;
typedef CGAL::AABB_face_graph_triangle_primitive<Polyhedron> Primitive;
typedef CGAL::AABB_traits<K, Primitive> Traits;
typedef CGAL::AABB_tree<Traits> Tree;
typedef CGAL::Side_of_triangle_mesh<Polyhedron, K> Point_inside;

bool pointInside(Polyhedron &polyhedron, Point &query) {
    // Construct AABB tree with a KdTree
    Tree tree(faces(polyhedron).first, faces(polyhedron).second, polyhedron);
    tree.accelerate_distance_queries();
    // Initialize the point-in-polyhedron tester
    Point_inside inside_tester(tree);

    // Determine the side and return true if inside!
    return inside_tester(query) == CGAL::ON_BOUNDED_SIDE;
}

Answer 2

Assuming P is large and there are many p0, you should compute a 3D triangulation in which to do point location. This is a CGAL demo.

Answer 3

You can write down the points in the hull as columns in a matrix and then have a vector which tells you what combination of points to take:

(X1 X2)  (a)      (X1a + X2b)
(Y1 Y2)  (b)   =  (Y1a + Y2b)
(Z1 Z2)           (Z1a + Z2b)

To generate your target point you want to find a vector that solves this, subject to the constraints that the elements of the vector are between 0 and 1, and that the elements of the vector add up to 1. You can solve this sort of problem - if there is a solution - via linear programming, which might involve using the http://en.wikipedia.org/wiki/Simplex_algorithm.

There are a variety of tricks to get this into the strict form. Adding another row to the matrix will allow you to say a + b = 1. To force b <= 1 you could have b + q = 1 and q >= 0, although as pointed out by Ted Hopp below, in this case b <= 1 because a + b = 1, and both a and b are >= 0.