使用Boost将Voronoi图转换为Delaunay图时，缺少具有非整数点坐标的三角形。

Question

以下是两个资源：

• Boost基础教程
• SO问题

我使用boost编写了Delaunay三角剖分。如果点的坐标是整数，它可以正常工作（我生成了几个随机测试，并没有发现错误）。但是，如果点是非整数，我发现许多错误的三角剖分，缺少边缘或错误的边缘。

例如，这张图片使用四舍五入的值构建，是正确的（请参见下面的代码）

但是，这张图片使用原始值构建，是不正确的（请参见下面的代码）

此代码重现了这两个示例（不包括显示）。

#include <boost/polygon/voronoi.hpp>
using boost::polygon::voronoi_builder;
using boost::polygon::voronoi_diagram;

struct Point
{
  double a;
  double b;
  Point(double x, double y) : a(x), b(y) {}
};

namespace boost
{
  namespace polygon
  {
    template <>
    struct geometry_concept<Point>
    {
      typedef point_concept type;
    };

    template <>
    struct point_traits<Point>
    {
      typedef double coordinate_type;

      static inline coordinate_type get(const Point& point, orientation_2d orient)
      {
        return (orient == HORIZONTAL) ? point.a : point.b;
      }
    };
  }  // polygon
}  // boost

int main()
{ 

 std::vector<double> xx = {8.45619987, 9.96573889, 0.26574428, 7.63918524, 8.15187618, 0.09100718};
 std::vector<double> yy = {9.2452883, 7.0843455, 0.4811701, 3.1193826, 5.1336435, 5.5368374};

 // std::vector<double> xx = {8, 10, 0, 8, 8, 0};
 // std::vector<double> yy = {9, 7, 0, 3, 5, 6};

  std::vector<Point> points;

  for (int i = 0 ; i < xx.size() ; i++)
  {
    points.push_back(Point(xx[i], yy[i]));
  }

  // Construction of the Voronoi Diagram.
  voronoi_diagram<double> vd;
  construct_voronoi(points.begin(), points.end(), &vd);

  for (const auto& vertex: vd.vertices())
  {
    std::vector<Point> triangle;
    auto edge = vertex.incident_edge();
    do
    {
      auto cell = edge->cell();
      assert(cell->contains_point());

      triangle.push_back(points[cell->source_index()]);
      if (triangle.size() == 3)
      {   
        // process output triangles
        std::cout << "Got triangle: (" << triangle[0].a << " " << triangle[0].b << ") (" << triangle[1].a << " " << triangle[1].b << ") (" << triangle[2].a << " " << triangle[2].b << ")" << std::endl;
        triangle.erase(triangle.begin() + 1);
      }

      edge = edge->rot_next();
    } while (edge != vertex.incident_edge());
  }

  return 0;
}

- JRR

---

可能是boost库中的一个bug。他们的bug跟踪器中有一些票据看起来与你的问题相似：Voronoi图中缺少边和polygon/voronoi缺少边。 - Vincent Saulue-Laborde

“Decimal” 是什么意思？ - Marc Glisse

你能展示一下输入数据应该是什么样子吗？因为无论我怎么看这个输入（即使是整数输入），它看起来都像是无效的多边形。 - sehe

@MarcGlisse 对，这实际上是我的英语问题。我想说的是积分。谢谢。 - JRR

@JRR 不是数字长什么样子的问题。对应的形状会是什么样子呢？我已经看过这些数字了。但是我和我使用的软件都无法从中找出一个有效的多边形。 - sehe

Answer 1

如果点的坐标不是小数，它可以正常工作。

在尝试了您的示例后，我突然意识到您并不是指“当坐标不是小数时”。您是指“整数”。这是很大的区别。

文档：点概念

坐标类型应为整数。

和

并非所有算法都支持浮点坐标类型，并且通常不适合目前与库一起使用。

- sehe

---

谢谢，这就是我在文档中了解到的（特别是voronoi_builder）。实际上，我只需要2位数字精度，所以我将输入乘以100并转换为整数。 - JRR

Answer 2

我已经仔细查看了您的输入数据。从中，我只能“想象”出两个有效的多边形：

1. 环 { {0,0}, {8, 3}, {10, 7}, {8, 9}, {0, 6}, }
   
2. 环 { {0,0}, {8, 3}, {8, 5}, {10, 7}, {8, 9}, {0, 6}, }
   

让我们在代码中定义它们：

Ring const inputs[] = {
            Ring { {0,0}, {8, 3}, {10, 7}, {8, 9}, {0, 6}, }, // {0, 0},
            Ring { {0,0}, {8, 3}, {8, 5}, {10, 7}, {8, 9}, {0, 6}, } // {0, 0},
        };

The closing points commented out are in case you have a polygon model that requires the polygon be closed.

In this case, I've opted fro Boost Geometries polygon model, and parameterized it to be not-closed:

static constexpr bool closed_polygons = false;
using bgPoly  = bgm::polygon<Point, false, closed_polygons>;
using bgMulti = bgm::multi_polygon<bgPoly>;
using Ring    = bgPoly::ring_type;

让我们进行测试

1. To create the test-cases that are not using integral numbers, let's transform the polygon by shifting it from (0,0) to (1,1) and also scaling every dimension by a factor of π.

2. let's also check for input validity (and optionally attempt to correct for errors):

   template <typename G> void validate(std::string name, G& geom) {
       std::cout << name << ": " << bg::wkt(geom) << "\n";
   
       std::string reason;
       if (!bg::is_valid(geom, reason)) {
           std::cout << name << ": " << reason << "\n";
   
           bg::correct(geom);
   
           std::cout << bg::wkt(geom) << "\n";
           if (!bg::is_valid(geom, reason)) {
               std::cout << name << " corrected: " << reason << "\n";
           }
       }
   }
   

3. Finally, let's save some SVG visualizations of the input and triangulations

演示程序

在 Coliru 上运行

#include <boost/polygon/voronoi.hpp>
#include <cassert>
#include <iostream>
using boost::polygon::voronoi_builder;
using boost::polygon::voronoi_diagram;

struct Point
{
    double a;
    double b;
    Point(double x = 0, double y = 0) : a(x), b(y) {}
};

namespace boost { namespace polygon {
    template <> struct geometry_concept<Point> { typedef point_concept type; };

    template <> struct point_traits<Point> {
        typedef double coordinate_type;

        static inline coordinate_type get(const Point& point, orientation_2d orient) {
            return (orient == HORIZONTAL) ? point.a : point.b;
        }
    };
} }

#include <boost/geometry.hpp>
#include <boost/geometry/geometries/point_xy.hpp>
#include <boost/geometry/algorithms/convex_hull.hpp>
#include <boost/geometry/algorithms/transform.hpp>
#include <boost/geometry/strategies/transform.hpp>
#include <boost/geometry/geometries/polygon.hpp>
#include <boost/geometry/geometries/multi_polygon.hpp>
#include <boost/geometry/geometries/register/point.hpp>
#include <boost/geometry/io/io.hpp>
#include <fstream>

namespace bg  = boost::geometry;
namespace bgm = bg::model;
namespace bgs = bg::strategy;

BOOST_GEOMETRY_REGISTER_POINT_2D(Point, double, bg::cs::cartesian, a, b)

static constexpr bool closed_polygons = false;
using bgPoly  = bgm::polygon<Point, false, closed_polygons>;
using bgMulti = bgm::multi_polygon<bgPoly>;
using Ring    = bgPoly::ring_type;

template <typename G> void validate(std::string name, G& geom) {
    std::cout << name << ": " << bg::wkt(geom) << "\n";

    std::string reason;
    if (!bg::is_valid(geom, reason)) {
        std::cout << name << ": " << reason << "\n";

        bg::correct(geom);

        std::cout << bg::wkt(geom) << "\n";
        if (!bg::is_valid(geom, reason)) {
            std::cout << name << " corrected: " << reason << "\n";
        }
    }
}

int main()
{
    int count = 0;

    Ring const inputs[] = {
                Ring { {0,0}, {8, 3}, {10, 7}, {8, 9}, {0, 6}, }, // {0, 0},
                Ring { {0,0}, {8, 3}, {8, 5}, {10, 7}, {8, 9}, {0, 6}, } // {0, 0},
            };

    bgs::transform::matrix_transformer<double, 2, 2> const transformations[] = {
            { 1,    0,    0, // identity transformation
              0,    1,    0,
              0,    0,    1 },
            { M_PI, 0,    1, // just to get nice non-integral numbers everywhere
              0,    M_PI, 1, // shift to (1,1) and scale everything by π
              0,    0,    1 },
            };

    for (auto transformation : transformations) {
        for (auto input : inputs) {
            validate("Input", input);

            Ring transformed_input;
            bg::transform(input, transformed_input, transformation);

            validate("transformed_input", transformed_input);

            // Construction of the Voronoi Diagram.
            voronoi_diagram<double> vd;
            construct_voronoi(transformed_input.begin(), transformed_input.end(), &vd);

            bgMulti out;
            Ring triangle;

            for (const auto& vertex: vd.vertices()) {
                triangle.clear();
                for(auto edge = vertex.incident_edge(); triangle.empty() || edge != vertex.incident_edge(); edge = edge->rot_next()) {
                    triangle.push_back(transformed_input[edge->cell()->source_index()]);

                    if (triangle.size() == 3) {

#if 0
                        std::cout << " -- found \n";
                        bgPoly t{triangle};
                        validate("Triangle", t);
                        out.push_back(t);
#else
                        out.push_back({ triangle });
#endif

                        triangle.erase(triangle.begin() + 1);
                    }
                }
            }

            std::cout << "Out " << bg::wkt(out) << "\n";
            {
                std::ofstream svg("/tmp/svg" + std::to_string(++count) + ".svg");
                boost::geometry::svg_mapper<Point> mapper(svg, 600, 600);

                mapper.add(out);
                mapper.map(out, R"(fill-opacity:0.5;fill:rgb(153,204,0);stroke:rgb(153,204,0);stroke-dasharray=5,5;stroke-width:2)");

                mapper.add(transformed_input);
                mapper.map(transformed_input, R"(fill-opacity:0.1;fill:rgb(204,153,0);stroke:red;stroke-width:3)");
            }
        } // inputs
    } // transformations
}

输出结果：

Input: POLYGON((0 0,8 3,10 7,8 9,0 6))
transformed_input: POLYGON((0 0,8 3,10 7,8 9,0 6))
Out MULTIPOLYGON(((0 6,0 0,8 3,0 6)),((8 9,0 6,8 3,8 9)),((10 7,8 9,8 3,10 7)))
Input: POLYGON((0 0,8 3,8 5,10 7,8 9,0 6))
transformed_input: POLYGON((0 0,8 3,8 5,10 7,8 9,0 6))
Out MULTIPOLYGON(((0 6,0 0,8 3,0 6)),((8 5,0 6,8 3,8 5)),((8 9,0 6,8 5,8 9)),((8 9,8 5,10 7,8 9)),((10 7,8 5,8 3,10 7)))
Input: POLYGON((0 0,8 3,10 7,8 9,0 6))
transformed_input: POLYGON((1 1,26.1327 10.4248,32.4159 22.9911,26.1327 29.2743,1 19.8496))
Out MULTIPOLYGON(((1 19.8496,1 1,26.1327 10.4248,1 19.8496)),((26.1327 29.2743,1 19.8496,26.1327 10.4248,26.1327 29.2743)),((32.4159 22.9911,26.1327 29.2743,26.1327 10.4248,32.4159 22.9911)))
Input: POLYGON((0 0,8 3,8 5,10 7,8 9,0 6))
transformed_input: POLYGON((1 1,26.1327 10.4248,26.1327 16.708,32.4159 22.9911,26.1327 29.2743,1 19.8496))
Out MULTIPOLYGON(((1 19.8496,1 1,26.1327 10.4248,1 19.8496)),((26.1327 16.708,1 19.8496,26.1327 10.4248,26.1327 16.708)),((26.1327 29.2743,1 19.8496,26.1327 16.708,26.1327 29.2743)),((26.1327 29.2743,26.1327 16.708,32.4159 22.9911,26.1327 29.2743)),((32.4159 22.9911,26.1327 16.708,26.1327 10.4248,32.4159 22.9911)))

对应的SVG文件：

- sehe

---

https://imgur.com/a/eGtSwlf 包含150dpi细节的SVG渲染，SVG在实时Coliru输出中。 - sehe

缩小100倍规模（http://coliru.stacked-crooked.com/a/e6e257ec975bad7d）会遇到问题，导致输出为空 https://imgur.com/a/akTzxmM - sehe

Answer 3

当点的尺寸太小，比如(0.411, 0.633)，(0.142, 0.453)等，我也遇到了同样的问题。

所以我使用point.x * 100, point.y * 100来放大这些点，然后boost::polygon::voronoi就可以正常工作，可能是由于浮点精度造成的。

- miyan