Delaunay from Voronoi with boost: missing triangle with non-integral point coordinates

Question

Following this two resources:

• Boost basic tutorial
• SO Question

I wrote a Delaunay triangulation with boost. It works fine if the points coordinates are integral (I generated several random tests and I did not observed error). However if the points are non integral I found many incorrect triangulations with missing edges or wrong edges.

For example this image has been build with rounded value and is correct (see code below)

But this image as been build with raw values and is incorrect (see code below)

This code reproduce this two examples (without the display).

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

Susannsusanna

Might be a bug in the boost library. Some tickets in their bugtracker look similar to your problem: Missing edges in Voronoi Diagram and polygon/voronoi missing edges – Kaolinite

Probably. Bad new because these bug reports are respectively one year and a half and two years old... – Susannsusanna

What do you mean by "decimal"? – Rosary

Computational geometry with floating point numbers is very hard. Don't expect boost to do it for you. Use integral or rational coordinates. – Hierarchy

Can you show us what the input data should look like? Because any way I look at the input (even the integer input) it looks like invalid polygons anyway. – Koheleth

@MarcGlisse right, it is a problem with my English actually. I meant integral. Thanks – Susannsusanna

@Koheleth the input data really look like the one in the code. – Susannsusanna

@Susannsusanna No no. Not what the NUMBERS look like. What would a corresponding shape look like? I've seen the numbers. And neither me, nor the software I used can figure out a valid polygon from it. – Koheleth

Answer 1

It works fine if the points coordinates are not decimal

After playing around with your sample I suddenly realized you didn't mean "when the coordinates are not decimal". You meant "integral". Big difference.

Documentation: Point Concept

The coordinate type is expected to be integral

and

Floating point coordinate types are not supported by all the algorithms and generally not suitable for use with the library at present.

Koheleth

Thank you, this is what I understood digging in the doc (especially the voronoi_builder. Actually I only need a 2 digits accuracy so I multiplied the input by 100 and I transtyped to integers. – Susannsusanna

Answer 2

I have looked a lot more at your input data. There's only two valid polygons I can "imagine" from that:

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

Let's define them in code:

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;

Let's do the tests

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

Demo Program

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
}

The output:

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)))

And the corresponding SVG files:

Koheleth

imgur.com/a/eGtSwlf has the svg renders in 150dpi detail, SVGs in the live Coliru output – Koheleth

Scaling down by 100x does run into the problems, leaving the output empty imgur.com/a/akTzxmM – Koheleth

Answer 3

I have the same issue when the points are too small, like (0.411, 0.633), (0.142, 0.453), etc.

So I enlarge the points with point.x * 100, point.y * 100, and then the boost::polygon::voronoi works well, maybe this is due to the floating precision.

Chirp