确定一个点是否在多面体内部答案

【问题标题】：Determining if a point is inside a polyhedron确定一个点是否在多面体内部
【发布时间】：2012-02-11 06:25:55
【问题描述】：

我正在尝试确定一个特定点是否位于多面体内。在我当前的实现中，我正在研究的方法是我们正在寻找多面体的面数组（在这种情况下是三角形，但以后可能是其他多边形）。我一直在尝试使用此处找到的信息：http://softsurfer.com/Archive/algorithm_0111/algorithm_0111.htm

下面，您将看到我的“内部”方法。我知道 nrml/normal 有点奇怪.. 这是旧代码的结果。当我运行它时，无论我给它什么输入，它似乎总是返回 true。（这已经解决了，请看下面我的回答——这段代码现在可以工作了）。

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++ 编写的，但正如 cmets 中所提到的，它确实与语言无关。

【问题讨论】：

您应该更新此问题，使其与语言无关。您要问的内容并非特定于 openGL 或 C++。一旦你有了一个通用的理论，你就可以让它适应你想要的每种语言和 3D API
创建一个简单的案例，您可以在其中验证它不在对象内，然后开始调试它。快速浏览后代码看起来差不多......
这看起来不是一个非常健壮的方法。首先，它只适用于凸多面体。其次，对于各种边界情况（选择的射线位于其中一个面的平面等），它都会失败。
我不担心凹多面体，所以没关系。但是，我愿意接受有关如何捕获更多边界情况的建议。
@duedl0r，感谢您提供的本应是显而易见的方法。听取这个简单的建议是我找到解决方案的原因。

标签： c++ 3d computational-geometry polyhedra

【解决方案1】：

您问题中的链接已过期，我无法从您的代码中理解算法。假设您有一个凸多面体，其面为 逆时针（从外部看），那么检查您的点是否在所有面的后面就足够了。为此，您可以将向量从点带到每个面，并检查标量积的符号与面的法线。如果为正，则该点在脸的后面；如果为零，则该点在脸上；如果是负数，则点在人脸前面。

这是一些完整的 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

【讨论】：

感谢您的简单解释，但我认为有问题。逆时针顶点排序基本上意味着法向量朝向多面体的外侧。因此，对于在后面（内部）的点，叉积应该是负数（连接点和三角形的向量相对于法向量大于 90 度）。还是我错过了什么？
因此，确实一个面的法线向量指向外部，但对于从多面体内部任意点p 到面f 的线段p2f 也是如此。标量（或点）积则为正，因为它小于 90°。也许你错过了p2f。看看isInConvexPoly。有p2f（而不是p）用于正常的点积。
哦！现在我明白了，当我进行数学计算时，我考虑了从面到点的向量（f2p），但是，是的，如果你假设这个约定，你是对的 :) 谢谢！

【解决方案2】：

如果你的网格是凹的，而且不一定是防水的，那就很难做到。

第一步，在网格表面上找到最接近该点的点。您需要跟踪位置和特定特征：最近的点是在面的中间，在网格的边缘，还是在网格的顶点之一。

如果特征是脸，你很幸运，可以使用绕组来确定是在里面还是在外面。计算法线到面（甚至不需要对其进行归一化，非单位长度即可），然后计算 dot( normal, pt - tri[0] ) 其中 pt 是您的点，tri[0] 是面的任何顶点。如果面具有一致的缠绕，则该点积的符号会告诉您它是在里面还是在外面。

如果特征是边缘，则计算两个面的法线（通过对叉积进行归一化），将它们相加，将其用作网格的法线，然后计算相同的点积。

最困难的情况是顶点是最近的特征。要计算该顶点处的网格法线，您需要计算共享该顶点的面的法线之和，由该顶点处该面的 2D 角度加权。例如，对于具有 3 个相邻三角形的立方体的顶点，权重将为 Pi/2。对于具有 6 个相邻三角形的立方体的顶点，权重将为 Pi/4。对于现实生活中的网格，每个面的权重会有所不同，在 [0 .. +Pi ] 范围内。这意味着您需要一些反三角代码来计算角度，可能是acos()。

如果您想知道为什么会这样，请参阅例如“Generating Signed Distance Fields From Triangle Meshes”，作者：J. Andreas Bærentzen 和 Henrik Aanæs。

【讨论】：

【解决方案3】：

事实证明，问题在于我阅读了上面链接中引用的算法。我正在阅读：

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

作为

N = - dot product of S and ni;

更改后，上面的代码现在似乎可以正常工作了。（我也在更新问题中的代码以反映正确的解决方案）。

【讨论】：