代码之家  ›  专栏  ›  技术社区  ›  JRR

带boost的voronoi中的delaunay:具有非积分点坐标的缺失三角形

  •  3
  • JRR  · 技术社区  · 7 年前

    遵循这两种资源:

    我写了一个Delaunay三角测量 boost . 如果点坐标是积分的(我生成了几个随机测试,但没有观察到误差),它就可以正常工作。但是,如果这些点是非积分的,我发现许多不正确的三角剖分缺少边或错误边。

    例如,此图像是用舍入值构建的,并且是正确的(请参见下面的代码)

    enter image description here

    但是这个图像是用原始值构建的,并且是不正确的(参见下面的代码)

    enter image description here

    这段代码复制了这两个示例(没有显示)。

    #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;
    }
    
    2 回复  |  直到 6 年前
        1
  •  5
  •   sehe    7 年前

    如果点坐标不是十进制的,它就可以工作了

    玩了你的样本之后,我突然意识到你不是说“坐标不是十进制的时候”。你的意思是“积分”。差别很大。

    Documentation: Point Concept

    坐标类型应为整数

    所有算法都不支持浮点坐标类型,目前一般不适合与库一起使用。

        2
  •  1
  •   sehe    7 年前

    我已经看了很多你的输入数据。我只能从中“想象”出两个有效的多边形:

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

    注释掉的结束点是在您有一个多边形模型需要关闭该多边形的情况下。

    在本例中,我选择了boost geometries多边形模型,并将其参数化为不闭合:

    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. 要创建不使用整数的测试用例,让我们通过将多边形从(0,0)移动到(1,1)来转换多边形,并将每个维度的比例缩放为。

    2. 我们还要检查输入的有效性(并尝试更正错误):

      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. 最后,让我们保存一些输入和三角剖分的SVG可视化

    演示程序

    Live On 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文件:

    enter image description here