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

Question

Я пытаюсь определить, лежит ли конкретная точка внутри многогранника. В моей текущей реализации метод, над которым я работаю, принимает точку, в которой мы ищем массив граней многогранника (в данном случае треугольники, но позже это могут быть другие многоугольники). Я пытался работать с информацией, найденной здесь: http://softsurfer.com/Archive/algorithm_0111/algorithm_0111.htm

Ниже вы увидите мой метод "изнутри". Я знаю, что nrml / normal довольно странная вещь.. это результат старого кода. Когда я запускал эту программу, казалось, что она всегда возвращает истину, независимо от того, какой вклад я ей даю. (Это решено, пожалуйста, смотрите мой ответ ниже - этот код работает сейчас).

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++, но, как уже упоминалось в комментариях, это действительно очень не зависит от языка.

Answer 1

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

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

как

N = - dot product of S and ni;

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

Answer 2

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

Вот некоторый полный код 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

Answer 3

Я уже отвечал на этот вопрос пару лет назад. Но с тех пор я обнаружил намного лучший алгоритм. Его изобрели в 2018 году, вот ссылка.

Идея довольно проста. Учитывая эту конкретную точку, вычислите сумму телесных углов со знаком всех граней многогранника, если смотреть из этой точки. Если точка находится снаружи, эта сумма должна быть равна нулю. Если точка находится внутри, эта сумма должна быть ±4·π стерадиан, + или - зависит от порядка намотки граней многогранника.

Этот конкретный алгоритм упаковывает многогранник в дерево, что значительно повышает производительность, когда вам нужно несколько внутренних/внешних запросов для одного и того же многогранника. Алгоритм вычисляет телесные углы только для отдельных граней, когда грань находится очень близко к точке запроса. Для больших наборов граней, удаленных от точки запроса, алгоритм вместо этого использует аппроксимацию этих наборов, используя некоторые числа, которые они сохраняют в узлах того дерева BVH, которое они строят из исходной сетки.

При ограниченной точности математики FP и при использовании аппроксимированных потерь дерева BVH из аппроксимации этот угол никогда не будет точно равен 0 или ±4·π. Но тем не менее порог 2·π на практике работает довольно хорошо, по крайней мере, по моему опыту. Если абсолютное значение этой суммы телесных углов меньше 2·π, считают, что точка находится снаружи.

Answer 4

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

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

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

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

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

Если вы хотите знать, почему это работает, см., Например, " Генерация полей расстояний со знаком из треугольных сеток " Дж. Андреаса Бёрентцена и Хенрика Аанеса.