Determinar si un punto está dentro de un poliedro

Question

Estoy intentando determinar si un punto específico se encuentra dentro de un poliedro. En mi implementación actual, el método en el que estoy trabajando Toma el punto que estamos buscando una variedad de caras del poliedro (triángulos en este caso, pero podría ser otros polígonos más adelante). He estado tratando de trabajar a partir de la información que se encuentra aquí: http://softsurfer.com/archive/algorithm_0111/algorithm_0111.htm

A continuación, verá mi método "interior". Sé que la cosa NRML/Normal es un poco raro ... es el resultado del código antiguo. Cuando estaba ejecutando esto, parecía volver siempre verdadero sin importar qué aporte le diera. (Esto se resuelve, consulte mi respuesta a continuación: este código funciona ahora).

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;
}

Esto está en C ++, pero como se menciona en los comentarios, es realmente un lenguaje agnóstico.

Answer 1

Resulta que el problema era mi lectura del algoritmo a los que se hace referencia en el enlace anterior. Estaba leyendo:

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

como

N = - dot product of S and ni;

Habiendo cambiado esto, el código anterior ahora parece funcionar correctamente. (También estoy actualizando el código en la pregunta para reflejar la solución correcta).

Answer 2

El enlace en su pregunta ha expirado y no pude entender el algoritmo de su código. Asumiendo que tienes un convexo poliedro con en sentido anti-horario caras orientadas (vistas desde el exterior), debe ser suficiente verificar que su punto esté detrás de todas las caras. Para hacer eso, puede tomar el vector desde el punto hasta cada cara y verificar el signo del producto escalar con la cara normal. Si es positivo, el punto está detrás de la cara; Si es cero, el punto está en la cara; Si es negativo, el punto está frente a la cara.

Aquí hay un código C ++ 11 completo, que funciona con caras de 3 puntos o caras más pequeñas (solo se consideran los primeros 3 puntos). Puedes cambiar fácilmente bound para excluir los límites.

#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;
}

Compilarlo con su compilador favorito

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

y probarlo como

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