Déterminer si un point est à l'intérieur d'un polyèdre

Question

J'essaie de déterminer si un point spécifique se trouve à l'intérieur d'un polyèdre.Dans mon implémentation actuelle, la méthode sur laquelle je travaille prend le point où nous recherchons un tableau des faces du polyèdre (triangles dans ce cas, mais cela pourrait être d'autres polygones plus tard).J'ai essayé de travailler à partir des informations trouvées ici: http://softsurfer.com/Archive/algorithm_0111 / algorithm_0111.htm

Ci-dessous, vous verrez ma méthode "inside".Je sais que le truc nrml / normal est un peu bizarre ... c'est le résultat d'un vieux code.Quand j'exécutais cela, il semblait toujours retourner vrai, quelle que soit l'entrée que je lui donne.(Ceci est résolu, s'il vous plaît voir ma réponse ci-dessous - ce code fonctionne maintenant).

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

C'est en C ++ mais comme mentionné dans les commentaires, c'est vraiment très indépendant du langage.

Answer 1

Il s'avère que le problème était ma lecture de l'algorithme référencé dans le lien ci-dessus.Je lisais:

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

comme

N = - dot product of S and ni;

Après avoir changé cela, le code ci-dessus semble maintenant fonctionner correctement.(Je mets également à jour le code de la question pour refléter la bonne solution).

Answer 2

Le lien dans votre question a expiré et je n'ai pas pu comprendre l'algorithme de votre code.En supposant que vous ayez un polyèdre convexe avec des faces orientées dans le sens antihoraire (vu de l'extérieur), il devrait être suffisant de vérifier que votre point est derrière toutes les faces.Pour ce faire, vous pouvez prendre le vecteur du point à chaque face et vérifier le signe du produit scalaire avec la normale du visage.S'il est positif, le point est derrière le visage;s'il est nul, le point est sur la face;s'il est négatif, le point est devant le visage.

Voici un code C ++ 11 complet, qui fonctionne avec des faces à 3 points ou des faces à plus de points simples (seuls les 3 premiers points sont considérés).Vous pouvez facilement modifier le bound pour exclure les limites.

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

Compilez-le avec votre compilateur préféré

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

et testez-le comme

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