Boost::Geometry - 在 3d 空间中查找 2d 多边形的区域？

Question

我正在尝试在 3d 空间中获取 2d 多边形的面积。 Boost::Geometry 有没有办法做到这一点？这是我的实现，但它一直返回 0：

#include <iostream>

#include <boost/geometry.hpp>
#include <boost/geometry/geometries/point_xy.hpp>
#include <boost/geometry/geometries/polygon.hpp>
#include <boost/geometry/io/wkt/wkt.hpp>

namespace bg = boost::geometry;

typedef bg::model::point<double, 3, bg::cs::cartesian> point3d;

int main()
{
    bg::model::multi_point<point3d> square;
    bg::read_wkt("MULTIPOINT((0 0 0), (0 2 0), (0 2 2), (0 0 2), (0 0 0))", square);
    double area = bg::area(square);
    std::cout << "Area: " << area << std::endl;

    return 0;
}

UPD：实际上，我对简单的 2d 多点正方形也有同样的问题：

#include <iostream>

#include <boost/geometry.hpp>
#include <boost/geometry/geometries/point_xy.hpp>
#include <boost/geometry/geometries/polygon.hpp>
#include <boost/geometry/io/wkt/wkt.hpp>

namespace bg = boost::geometry;

typedef bg::model::point<double, 2, bg::cs::cartesian> point2d;

int main()
{
    bg::model::multi_point<point2d> square;
    bg::read_wkt("MULTIPOINT((0 0), (2 0), (2 2), (0 2))", square);
    double area = bg::area(square);
    std::cout << "Area: " << area << std::endl;

    return 0;
}

结果如下：

$ ./test_area
Area: 0

UPD：看起来像boost::geometry 中的面积计算仅适用于二维多边形。

Answer 1

我不熟悉 boost 的几何部分，但根据我的几何知识，我可以说 3D 和 2D 并没有太大区别。尽管可能已经有一些东西在 boost 中，但您可以编写自己的方法来相当容易地做到这一点。

编辑：

da code monkey 指出shoelace formula 这样会更高效，因为它更简单，速度更快。

原意如下：

---

为了计算这个，我首先将多边形细分为三角形，因为任何多边形都可以拆分为多个三角形。我会取每个三角形，并计算每个三角形的面积。要在 3d 空间中执行此操作，应用相同的概念。求底，取△ABC，任意指定-AB为底，-BC为高，-CA为斜边。只需执行 (-AB*-BC)/2 。只需将每个三角形的面积相加即可。

我不知道 boost 是否有内置的 tessellate 方法，这在 C++ 中实现起来相当困难，但您可能想要创建某种三角形风扇。 （注意：这只适用于凸多边形）。如果你确实有一个凹多边形，你应该看看这个：http://www.cs.unc.edu/~dm/CODE/GEM/chapter.html 我会把它作为一个练习放在 c++ 中，但过程相当简单。

Answer 2

我不希望点的集合有一个区域。您将需要 model::polygon<poind3d> 的等效项，但目前似乎不受支持。

如果保证点共面且线段不相交，则可以将多边形分解为一系列三角形，并根据以下公式使用一点线性代数计算面积对于三角形的面积：

在非凸多边形的情况下，需要调整面积的总和以减去多边形外的面积。实现这一点的最简单方法是使用三角形的符号区域，包括右侧三角形的正贡献和左侧三角形的负贡献：

请注意，似乎有一些计划在 Boost 中包含 cross_product 实现，但它似乎从 1.56 版开始不包含在内。以下替换应该可以为您的用例解决问题：

point3d cross_product(const point3d& p1, const point3d& p2)
{
  double x = bg::get<0>(p1);
  double y = bg::get<1>(p1);
  double z = bg::get<2>(p1);
  double u = bg::get<0>(p2);
  double v = bg::get<1>(p2);
  double w = bg::get<2>(p2);
  return point3d(y*w-z*v, z*u-x*w, x*v-y*u);
}
point3d cross_product(const bg::model::segment<point3d>& p1
                    , const bg::model::segment<point3d>& p2)
{
  point3d v1(p1.second);
  point3d v2(p2.second);
  bg::subtract_point(v1, p1.first);
  bg::subtract_point(v2, p2.first);

  return cross_product(v1, v2);
}

然后可以使用以下方法计算面积：

// compute the are of a collection of 3D points interpreted as a 3D polygon
// Note that there are no checks as to whether or not the points are
// indeed co-planar.
double area(bg::model::multi_point<point3d>& polygon)
{
  if (polygon.size()<3) return 0;

  bg::model::segment<point3d> v1(polygon[1], polygon[0]);
  bg::model::segment<point3d> v2(polygon[2], polygon[0]);
  // Compute the cross product for the first pair of points, to handle
  // shapes that are not convex.
  point3d n1 = cross_product(v1, v2);
  double normSquared = bg::dot_product(n1, n1);
  if (normSquared > 0)
  {
    bg::multiply_value(n1, 1.0/sqrt(normSquared));
  }
  // sum signed areas of triangles
  double result = 0.0;
  for (size_t i=1; i<polygon.size(); ++i)
  {
    bg::model::segment<point3d> v1(polygon[0], polygon[i-1]);
    bg::model::segment<point3d> v2(polygon[0], polygon[i]);

    result += bg::dot_product(cross_product(v1, v2), n1);
  }
  result *= 0.5;
  return abs(result);
}