Determining if a point is inside a polyhedron

Question

I'm attempting to determine if a specific point lies inside a polyhedron. In my current implementation, the method I'm working on take the point we're looking for an array of the faces of the polyhedron (triangles in this case, but it could be other polygons later). I've been trying to work from the info found here: http://softsurfer.com/Archive/algorithm_0111/algorithm_0111.htm

Below, you'll see my "inside" method. I know that the nrml/normal thing is kind of weird .. it's the result of old code. When I was running this it seemed to always return true no matter what input I give it. (This is solved, please see my answer below -- this code is working now).

bool Container::inside(Point* point, float* polyhedron[3], int faces) {
  Vector* dS = Vector::fromPoints(point->X, point->Y, point->Z,
                 100, 100, 100);
  int T_e = 0;
  int T_l = 1;

  for (int i = 0; i < faces; i++) {
    float* polygon = polyhedron[i];

    float* nrml = normal(&polygon[0], &polygon[1], &polygon[2]);
    Vector* normal = new Vector(nrml[0], nrml[1], nrml[2]);
    delete nrml;

    float N = -((point->X-polygon[0][0])*normal->X + 
                (point->Y-polygon[0][1])*normal->Y +
                (point->Z-polygon[0][2])*normal->Z);
    float D = dS->dot(*normal);

    if (D == 0) {
      if (N < 0) {
        return false;
      }

      continue;
    }

    float t = N/D;

    if (D < 0) {
      T_e = (t > T_e) ? t : T_e;
      if (T_e > T_l) {
        return false;
      }
    } else {
      T_l = (t < T_l) ? t : T_l;
      if (T_l < T_e) {
        return false;
      }
    }
  }

  return true;
}

This is in C++ but as mentioned in the comments, it's really very language agnostic.

Nicosia

You should update this question so that it is more language agnostic. What you are asking is not specific to openGL or C++. Once you have a general theory you can adapt it to what every language and 3D API you want – Gastroenterology

create a simple case where you can verify it's not inside the object, then start debugging it..the code looks about right after quickly skimming through it.. – Meridional

This doesn't look like a very robust method. First, it will only work with convex polyhedra. Second, it will fail for all kinds of boundary cases (chosen ray lies in the plane of one of the faces, etc). – Cacogenics

I'm not worried about concave polyhedra, so that doesn't matter. However, I am open to suggestions on how to catch more boundary cases. – Nicosia

@duedl0r, Thanks for what should have been an obvious approach. Heeding that simple advice is what lead me to my solution. – Nicosia

@greggory.hz Assuming variable names are related to geometric objects: Simplest convex polyhedron has 4 faces but your code has polyhedron structure with only 3 faces. – Villager

Answer 1

It turns out that the problem was my reading of the algorithm referenced in the link above. I was reading:

N = - dot product of (P0-Vi) and ni;

as

N = - dot product of S and ni;

Having changed this, the code above now seems to work correctly. (I'm also updating the code in the question to reflect the correct solution).

Nicosia

Answer 2

The link in your question has expired and I could not understand the algorithm from your code. Assuming you have a convex polyhedron with counterclockwise oriented faces (seen from outside), it should be sufficient to check that your point is behind all faces. To do that, you can take the vector from the point to each face and check the sign of the scalar product with the face's normal. If it is positive, the point is behind the face; if it is zero, the point is on the face; if it is negative, the point is in front of the face.

Here is some complete C++11 code, that works with 3-point faces or plain more-point faces (only the first 3 points are considered). You can easily change bound to exclude the boundaries.

#include <vector>
#include <cassert>
#include <iostream>
#include <cmath>

struct Vector {
  double x, y, z;

  Vector operator-(Vector p) const {
    return Vector{x - p.x, y - p.y, z - p.z};
  }

  Vector cross(Vector p) const {
    return Vector{
      y * p.z - p.y * z,
      z * p.x - p.z * x,
      x * p.y - p.x * y
    };
  }

  double dot(Vector p) const {
    return x * p.x + y * p.y + z * p.z;
  }

  double norm() const {
    return std::sqrt(x*x + y*y + z*z);
  }
};

using Point = Vector;

struct Face {
  std::vector<Point> v;

  Vector normal() const {
    assert(v.size() > 2);
    Vector dir1 = v[1] - v[0];
    Vector dir2 = v[2] - v[0];
    Vector n  = dir1.cross(dir2);
    double d = n.norm();
    return Vector{n.x / d, n.y / d, n.z / d};
  }
};

bool isInConvexPoly(Point const& p, std::vector<Face> const& fs) {
  for (Face const& f : fs) {
    Vector p2f = f.v[0] - p;         // f.v[0] is an arbitrary point on f
    double d = p2f.dot(f.normal());
    d /= p2f.norm();                 // for numeric stability

    constexpr double bound = -1e-15; // use 1e15 to exclude boundaries
    if (d < bound)
      return false;
  }

  return true;
}

int main(int argc, char* argv[]) {
  assert(argc == 3+1);
  char* end;
  Point p;
  p.x = std::strtod(argv[1], &end);
  p.y = std::strtod(argv[2], &end);
  p.z = std::strtod(argv[3], &end);

  std::vector<Face> cube{ // faces with 4 points, last point is ignored
    Face{{Point{0,0,0}, Point{1,0,0}, Point{1,0,1}, Point{0,0,1}}}, // front
    Face{{Point{0,1,0}, Point{0,1,1}, Point{1,1,1}, Point{1,1,0}}}, // back
    Face{{Point{0,0,0}, Point{0,0,1}, Point{0,1,1}, Point{0,1,0}}}, // left
    Face{{Point{1,0,0}, Point{1,1,0}, Point{1,1,1}, Point{1,0,1}}}, // right
    Face{{Point{0,0,1}, Point{1,0,1}, Point{1,1,1}, Point{0,1,1}}}, // top
    Face{{Point{0,0,0}, Point{0,1,0}, Point{1,1,0}, Point{1,0,0}}}, // bottom
  };

  std::cout << (isInConvexPoly(p, cube) ? "inside" : "outside") << std::endl;

  return 0;
}

Compile it with your favorite compiler

clang++ -Wall -std=c++11 code.cpp -o inpoly

and test it like

$ ./inpoly 0.5 0.5 0.5
inside
$ ./inpoly 1 1 1
inside
$ ./inpoly 2 2 2
outside

Enjoy

Thanks for the easy explanation, but I think there is something wrong. Counterclockwise vertices ordering basically means that the normal vector is towards the outside of the polyhedron. Therefore, for the point to be behind (inside) the cross product should be negative (the vector connecting the point and the triangle is more than 90 degrees with respect the normal vector). Or am I missing something? – Graecize

So, indeed the normal vector of a face points to the outside, but the same holds for the segment p2f from any point p inside the polyhedron to the face f. The scalar (or dot) product is then positive, since it is less then 90°. Maybe you missed p2f. Have a look at isInConvexPoly. There p2f is used (and not p) for the dot product with the normal. – Enjoy

Oh! Now I see, when I did the math I considered the vector from the face to the point (f2p), but yes, you are right if you assume that convention :) Thanks! – Graecize

Answer 3

If your mesh is concave, and not necessarily watertight, that’s rather hard to accomplish.

As a first step, find the point on the surface of the mesh closest to the point. You need to keep track the location, and specific feature: whether the closest point is in the middle of face, on the edge of the mesh, or one of the vertices of the mesh.

If the feature is face, you’re lucky, can use windings to find whether it’s inside or outside. Compute normal to face (don't even need to normalize it, non-unit-length will do), then compute dot( normal, pt - tri[0] ) where pt is your point, tri[0] is any vertex of the face. If the faces have consistent winding, the sign of that dot product will tell you if it’s inside or outside.

If the feature is edge, compute normals to both faces (by normalizing a cross-product), add them together, use that as a normal to the mesh, and compute the same dot product.

The hardest case is when a vertex is the closest feature. To compute mesh normal at that vertex, you need to compute sum of the normals of the faces sharing that vertex, weighted by 2D angles of that face at that vertex. For example, for vertex of cube with 3 neighbor triangles, the weights will be Pi/2. For vertex of a cube with 6 neighbor triangles the weights will be Pi/4. And for real-life meshes the weights will be different for each face, in the range [ 0 .. +Pi ]. This means you gonna need some inverse trigonometry code for this case to compute the angle, probably acos().

If you want to know why that works, see e.g. “Generating Signed Distance Fields From Triangle Meshes” by J. Andreas Bærentzen and Henrik Aanæs.

Northington

Answer 4

I have already answered this question couple years ago. But since that time I’ve discovered much better algorithm. It was invented in 2018, here’s the link.

The idea is rather simple. Given that specific point, compute a sum of signed solid angles of all faces of the polyhedron as viewed from that point. If the point is outside, that sum gotta be zero. If the point is inside, that sum gotta be ±4·π steradians, + or - depends on the winding order of the faces of the polyhedron.

That particular algorithm is packing the polyhedron into a tree, which dramatically improves performance when you need multiple inside/outside queries for the same polyhedron. The algorithm only computes solid angles for individual faces when the face is very close to the query point. For large sets of faces far away from the query point, the algorithm is instead using an approximation of these sets, using some numbers they keep in the nodes of that BVH tree they build from the source mesh.

With limited precision of FP math, and if using that approximated BVH tree losses from the approximation, that angle will never be exactly 0 nor ±4·π. But still, the 2·π threshold works rather well in practice, at least in my experience. If the absolute value of that sum of solid angles is less than 2·π, consider the point to be outside.

Northington

This is the most helpful answer, by far – Tired