Как определить, находится ли точка внутри многогранника

Срок действия ссылки в вашем вопросе истек, и я не смог понять алгоритм из вашего кода. Предполагая, что у вас есть выпуклый многогранник с гранями, ориентированными против часовой стрелки (если смотреть снаружи), этого должно быть достаточно, чтобы убедиться, что ваша точка находится за всеми гранями. Для этого вы можете взять вектор от точки к каждой грани и сверить знак скалярного произведения с нормалью грани. Если он положительный, точка находится за лицом; если он равен нулю, точка находится на грани; если он отрицательный, точка находится перед лицом.

Вот некоторый полный код C++11, который работает с гранями с 3 точками или простыми гранями с большим количеством точек (учитываются только первые 3 точки). Вы можете легко изменить bound, чтобы исключить границы.

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

Скомпилируйте его вашим любимым компилятором

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

и протестируй как

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

Если ваша сетка вогнутая и не обязательно водонепроницаемая, этого довольно сложно добиться.

В качестве первого шага найдите ближайшую к точке точку на поверхности сетки. Вам нужно отслеживать местоположение и конкретную особенность: находится ли ближайшая точка в середине грани, на краю сетки или в одной из вершин сетки.

Если особенность лица, вам повезло, вы можете использовать обмотки, чтобы определить, находится ли она внутри или снаружи. Вычислите нормаль к лицу (даже не нужно нормализовать его, подойдет не единичная длина), затем вычислите dot( normal, pt - tri[0] ), где pt — ваша точка, tri[0] — любая вершина лица. Если грани имеют последовательную обмотку, знак этого скалярного произведения скажет вам, находится ли она внутри или снаружи.

Если объект является краем, вычислите нормали к обеим граням (путем нормализации векторного произведения), сложите их вместе, используйте это как нормаль к сетке и вычислите то же скалярное произведение.

Самый сложный случай, когда вершина является ближайшим объектом. Чтобы вычислить нормаль сетки в этой вершине, вам нужно вычислить сумму нормалей граней, разделяющих эту вершину, взвешенных по двумерным углам этой грани в этой вершине. Например, для вершины куба с 3 соседними треугольниками веса будут Pi/2. Для вершины куба с 6 соседними треугольниками веса будут Pi/4. А для реальных сеток веса будут разными для каждой грани в диапазоне [ 0 .. +Pi ]. Это означает, что вам понадобится код обратной тригонометрии для этого случая, чтобы вычислить угол, возможно, acos().

Если вы хотите знать, почему это работает, см., например. «Создание полей расстояний со знаком из треугольных сеток», автор J. Andreas Bærentzen и Хенрик Аанес.

Оказывается, проблема заключалась в том, что я прочитал алгоритм, указанный в ссылке выше. Я читал:

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

as

N = - dot product of S and ni;

Изменив это, приведенный выше код теперь работает правильно. (Я также обновляю код в вопросе, чтобы отразить правильное решение).