Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

Determining if a point is inside a polyhedron

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.

like image 423
gregghz Avatar asked Jan 16 '12 09:01

gregghz


People also ask

How do you know if a point is inside a shape?

Draw a horizontal line to the right of each point and extend it to infinity. Count the number of times the line intersects with polygon edges. A point is inside the polygon if either count of intersections is odd or point lies on an edge of polygon.

How do you determine if a point is inside or outside a polygon?

One simple way of finding whether the point is inside or outside a simple polygon is to test how many times a ray, starting from the point and going in any fixed direction, intersects the edges of the polygon. If the point is on the outside of the polygon the ray will intersect its edge an even number of times.

How do you know if a point is inside a triangle?

The simplest way to determine if a point lies inside a triangle is to check the number of points in the convex hull of the vertices of the triangle adjoined with the point in question. If the hull has three points, the point lies in the triangle's interior; if it is four, it lies outside the triangle.

Can a polyhedron have a hole in it?

The faces of a polyhedron can have different sizes and shapes, just as long as each one is a polygon. The polyhedron itself can have a hole (or two or more).


1 Answers

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
like image 181
John Avatar answered Oct 25 '22 07:10

John