Determinare se un punto è all'interno di un poliedro
-
28-10-2019 - |
Domanda
Sto tentando di determinare se un punto specifico si trova all'interno di un poliedro. Nella mia attuale implementazione, il metodo a cui sto lavorando prendendo il punto in cui stiamo cercando una serie di facce del poliedro (triangoli in questo caso, ma potrebbero essere altri poligoni in seguito). Ho provato a lavorare dalle informazioni trovate qui: http://softsurfer.com/archive/algorithm_0111/algorithm_0111.htm
Di seguito, vedrai il mio metodo "dentro". So che la cosa NRML/normale è un po 'strana .. è il risultato del vecchio codice. Quando stavo eseguendo questo sembrava restituire sempre vero, indipendentemente dall'input che gli do. (Questo è risolto, consulta la mia risposta di seguito: questo codice ora funziona).
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;
}
Questo è in C ++ ma come menzionato nei commenti, è davvero molto agnostico.
Soluzione 2
Si scopre che il problema era la mia lettura dell'algoritmo a cui si fa riferimento nel link sopra. Stavo leggendo:
N = - dot product of (P0-Vi) and ni;
come
N = - dot product of S and ni;
Avendo cambiato questo, il codice sopra ora sembra funzionare correttamente. (Sto anche aggiornando il codice nella domanda per riflettere la soluzione corretta).
Altri suggerimenti
Il link nella tua domanda è scaduto e non sono riuscito a capire l'algoritmo dal tuo codice. Supponendo che tu abbia un convesso poliedro con Antiorario Facce orientate (visto dall'esterno), dovrebbe essere sufficiente verificare che il tuo punto sia dietro tutti i volti. Per fare ciò, puoi prendere il vettore dal punto a ciascuna faccia e controllare il segno del prodotto scalare con il normale. Se è positivo, il punto è dietro il viso; Se è zero, il punto è sul viso; Se è negativo, il punto è di fronte al viso.
Ecco un codice C ++ 11 completo, che funziona con facce a 3 punti o facce più semplici (vengono considerati solo i primi 3 punti). Puoi facilmente cambiare bound
per escludere i confini.
#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;
}
Compilalo con il tuo compilatore preferito
clang++ -Wall -std=c++11 code.cpp -o inpoly
e testalo come
$ ./inpoly 0.5 0.5 0.5
inside
$ ./inpoly 1 1 1
inside
$ ./inpoly 2 2 2
outside