Bestimmen, ob ein Punkt innerhalb eines Polyeders liegt

Question

Ich versuche zu bestimmen, ob ein bestimmter Punkt innerhalb eines Polyeders liegt. In meiner aktuellen Implementierung nimmt die Methode, an der ich arbeite, den Punkt, den wir suchen, ein Array der Flächen des Polyeders (Dreiecke in diesem Fall, aber es könnte später andere Polygone sein). Ich habe versucht, von den Informationen hier gefunden zu arbeiten: http://softsurfer.com/Archive/algorithm_0111/algorithm_0111.htm

Unten sehen Sie meine "innere" Methode. Ich weiß, dass die nrml/normal Sache ist irgendwie seltsam .. es ist das Ergebnis der alten Code. Als ich diesen Code ausführte, schien er immer true zurückzugeben, egal, welche Eingabe ich ihm gebe. (Dies ist gelöst, siehe meine Antwort unten - dieser Code funktioniert jetzt).

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

Dies ist in C++, aber wie in den Kommentaren erwähnt, ist es wirklich sehr Sprache agnostisch.

Answer 1

Der Link in Ihrer Frage ist abgelaufen und ich konnte den Algorithmus in Ihrem Code nicht verstehen. Angenommen, Sie haben eine konvex Polyeder mit gegen den Uhrzeigersinn orientierten Flächen (von außen gesehen), sollte es ausreichen, zu prüfen, ob Ihr Punkt hinter allen Flächen liegt. Dazu können Sie den Vektor vom Punkt zu jeder Fläche nehmen und das Vorzeichen des Skalarprodukts mit der Flächennormalen überprüfen. Ist es positiv, liegt der Punkt hinter der Fläche; ist es Null, liegt der Punkt auf der Fläche; ist es negativ, liegt der Punkt vor der Fläche.

Hier ist ein vollständiger C++11-Code, der mit 3-Punkt-Flächen oder einfachen Mehr-Punkt-Flächen (nur die ersten 3 Punkte werden berücksichtigt) funktioniert. Sie können leicht ändern bound um die Grenzen auszuschließen.

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

Kompilieren Sie es mit Ihrem bevorzugten Compiler

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

und testen Sie es wie

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

Answer 2

Wenn Ihr Netz konkav und nicht unbedingt wasserdicht ist, ist das ziemlich schwer zu bewerkstelligen.

Finden Sie zunächst den Punkt auf der Oberfläche des Netzes, der dem Punkt am nächsten liegt. Sie müssen sich die Position und die Besonderheit merken: ob der nächstgelegene Punkt in der Mitte der Fläche, am Rand des Netzes oder an einem der Scheitelpunkte des Netzes liegt.

Wenn das Merkmal eine Fläche ist, können Sie mit etwas Glück anhand der Windungen feststellen, ob es sich um eine Innen- oder Außenfläche handelt. Berechnen Sie die Normale zur Fläche (Sie brauchen sie nicht einmal zu normalisieren, die Nicht-Einheitslänge reicht aus), dann berechnen Sie dot( normal, pt - tri[0] ) wobei pt Ihr Punkt und tri[0] ein beliebiger Scheitelpunkt der Fläche ist. Wenn die Flächen konsistent gewickelt sind, sagt das Vorzeichen des Punktprodukts aus, ob sie innen oder außen liegen.

Handelt es sich bei dem Feature um eine Kante, berechnen Sie die Normalen beider Flächen (durch Normalisierung eines Kreuzprodukts), addieren Sie diese, verwenden Sie dies als Normalen für das Netz und berechnen Sie dasselbe Punktprodukt.

Der schwierigste Fall ist, wenn ein Scheitelpunkt das nächstgelegene Merkmal ist. Um die Netznormale an diesem Scheitelpunkt zu berechnen, müssen Sie die Summe der Normalen der Flächen berechnen, die diesen Scheitelpunkt teilen, gewichtet nach 2D-Winkeln dieser Fläche an diesem Scheitelpunkt. Für eine Würfelspitze mit 3 Nachbardreiecken beträgt die Gewichtung beispielsweise Pi/2. Für die Spitze eines Würfels mit 6 Nachbardreiecken sind die Gewichte Pi/4. Bei realen Netzen werden die Gewichte für jede Fläche unterschiedlich sein, im Bereich [ 0 .. +Pi ]. Das bedeutet, dass Sie für diesen Fall einen Code für inverse Trigonometrie benötigen, um den Winkel zu berechnen, wahrscheinlich acos() .

Wenn Sie wissen wollen, warum das funktioniert, lesen Sie z. B. " Erzeugung von Abstandsfeldern mit Vorzeichen aus Dreiecksnetzen " von J. Andreas Bærentzen und Henrik Aanæs.

Answer 3

Ich habe diese Frage bereits vor einigen Jahren beantwortet. Aber seither habe ich einen viel besseren Algorithmus entdeckt. Er wurde im Jahr 2018 erfunden, hier ist der Link .

Die Idee ist recht einfach. Berechnen Sie für diesen bestimmten Punkt eine Summe von vorzeichenbehafteten Raumwinkel aller Flächen des Polyeders, von diesem Punkt aus gesehen. Liegt der Punkt außerhalb, muss die Summe gleich Null sein. Befindet sich der Punkt im Inneren, muss die Summe ±4- Steradian betragen, wobei + oder - von der Reihenfolge der Windungen der Flächen des Polyeders abhängt.

Dieser spezielle Algorithmus packt das Polyeder in einen Baum, was die Leistung erheblich verbessert, wenn Sie mehrere Innen/Außen-Abfragen für dasselbe Polyeder benötigen. Der Algorithmus berechnet Raumwinkel für einzelne Flächen nur dann, wenn die Fläche sehr nahe am Abfragepunkt liegt. Bei großen Mengen von Flächen, die weit vom Abfragepunkt entfernt sind, verwendet der Algorithmus stattdessen eine Annäherung an diese Mengen, indem er einige Zahlen verwendet, die er in den Knoten des BVH-Baum die sie aus dem Quellennetz erstellen.

Bei der begrenzten Genauigkeit der FP-Mathematik und bei Verwendung der approximierten BVH-Baumverluste aus der Annäherung wird dieser Winkel niemals genau 0 oder ±4- sein. Dennoch funktioniert der Schwellenwert von 2 in der Praxis recht gut, zumindest nach meiner Erfahrung. Wenn der absolute Wert der Summe der Raumwinkel kleiner als 2- ist, ist der Punkt als außerhalb zu betrachten.

Answer 4

Wie sich herausstellte, lag das Problem bei meiner Lektüre des Algorithmus, auf den im obigen Link verwiesen wird. Ich habe gelesen:

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

als

N = - dot product of S and ni;

Nachdem ich dies geändert habe, scheint der obige Code nun korrekt zu funktionieren. (Ich aktualisiere auch den Code in der Frage, um die richtige Lösung wiederzugeben).