search

CGAL / cgal

The public CGAL repository, see the README below
https://github.com/CGAL/cgal#readme
Other
4.87k stars 1.38k forks source link

Voronoi_diagram_2 return a finite voronoi edge with a vertex having an invalid coordinate(super large). #4195

Closed diablohsh closed 5 years ago

diablohsh commented 5 years ago

Issue Details

Voronoi_diagram_2 return a finite voronoi edge with a vertex having a invalid coordinate(super large).

Source Code

#include <sstream>
#include <vector>

#define CGAL_HEADER_ONLY 1

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Segment_Delaunay_graph_filtered_traits_2.h>
#include <CGAL/Segment_Delaunay_graph_traits_2.h>
#include <CGAL/Segment_Delaunay_graph_2.h>
#include <CGAL/Segment_Delaunay_graph_adaptation_traits_2.h>
#include <CGAL/Segment_Delaunay_graph_adaptation_policies_2.h>
#include <CGAL/Voronoi_diagram_2.h>

using namespace std;

typedef CGAL::Exact_predicates_inexact_constructions_kernel             K;
typedef CGAL::Segment_Delaunay_graph_filtered_traits_2<K>               SDGT;
typedef CGAL::Segment_Delaunay_graph_2<SDGT>                            SDG;
typedef CGAL::Segment_Delaunay_graph_adaptation_traits_2<SDG>           SDGAT;
typedef CGAL::Segment_Delaunay_graph_degeneracy_removal_policy_2<SDG>   SDGAP;
typedef CGAL::Voronoi_diagram_2<SDG, SDGAT, SDGAP>                      VD;

typedef K::FT           Num_t;
typedef K::Point_2      Point_2;
typedef SDGAT::Site_2   Site_2;

int main()
{
    istringstream strm("2ftkqnron5d/2gosa7pa2gw 14o7lsvsnhx/18ce53un18g 24du0d5436p/2gosa7pa2gw oiv9ts2otv/m672jxbim8 1pj1zmqn3g3/2gosa7pa2gw 2g0t7s4y9vh/2gosa7pa2gw 1qvppf3jptx/2gosa7pa2gw p35x4z8lkr/m672jxbim8 1pmtfy222kr/4xdkkfek4xs 1fxltrt3fcj/18ce53un18g 1mxi0dc8tg3/4xdkkfek4xs bfx5ch9ua7/b33j9ynrb4 1iy0p7mruhd/9ur54ut49vk 1gfiq7inzsz/18ce53un18g 0 1 7v13lg71tl/jpia9pm8jr4 8whunyntav/b33j9ynrb4 194brmj62m5/9ur54ut49vk ral4ou2do3/18ce53un18g 3fj4l6xlwr/b33j9ynrb4 g557q50kj3/m672jxbim8 tqmz8g6vtn/18ce53un18g fga1rb5hd5/m672jxbim8 4lkriq10u1/5jjrmzbvnk 1chq45rxz7j/2gosa7pa2gw 1 40vl5jmjv1/5jjrmzbvnk");

    vector<Point_2> vecPoint;

    string line;
    while (true)
    {
        mpq_t x, y;
        mpq_init(x);
        mpq_init(y);

        if (!getline(strm, line, ' '))
            break;
        mpq_set_str(x, line.c_str(), 36);

        getline(strm, line, ' ');
        mpq_set_str(y, line.c_str(), 36);

        vecPoint.emplace_back(Num_t(mpq_get_d(x)), Num_t(mpq_get_d(y)));
    }

    vector<Site_2> vecSite;

    for (size_t i = 1; i < vecPoint.size(); i++)
        vecSite.push_back(Site_2::construct_site_2(vecPoint[i - 1], vecPoint[i]));

    vecSite.push_back(Site_2::construct_site_2(vecPoint.back(), vecPoint.front()));

    VD vd;
    vd.insert(vecSite.begin(), vecSite.end());

    for_each(vd.edges_begin(), vd.edges_end(), [&](const VD::Halfedge& e)
    {
        if (e.is_unbounded())
            return;
        auto s = e.source()->point();
        auto t = e.target()->point();
        cout << s << "," << t << endl;
    });

    return 0;
}

output with an invalid finite edge vertex having extremely large coordinate(2.78205e+10 2.22564e+11):

0.773244 1.12446,0.798476 1.19867
0.861265 1.10608,0.798476 1.19867
0.273577 1.17925,0.347474 1.17114
0.471248 1.01689,0.347474 1.17114
0.439592 0.181246,0.309526 0.728072
0.322629 0.87997,0.309526 0.728072
0.210585 0.888338,0.309526 0.728072
0.184313 0.857639,0.0110809 0.80305
0.332233 1.03213,0.471248 1.01689
2.78205e+10 2.22564e+11,0.471248 1.01689
0.814059 0.857826,0.861265 1.10608
0.709009 1.13151,0.773244 1.12446
0.693769 0.992491,0.773244 1.12446
0.127187 0.615504,0.184313 0.857639
0.200711 0.889078,0.184313 0.857639
0.554752 1.00773,0.693769 0.992491
0.800249 0.835478,0.693769 0.992491
0.502608 0.917517,2.78205e+10 2.22564e+11
0.554752 1.00773,2.78205e+10 2.22564e+11
0.174677 0.941349,0.154882 1.18237
0.258641 1.30485,0.154882 1.18237
0.709009 1.13151,0.554752 1.00773
0.174677 0.941349,0 1
0.200711 0.889078,0.174677 0.941349
0.273577 1.17925,0.258641 1.30485
0.502608 0.917517,0.322629 0.87997
0.210585 0.888338,0.322629 0.87997
0.273577 1.17925,0.332233 1.03213
0.210435 0.888349,0.332233 1.03213
0.210435 0.888349,0.200711 0.889078
0.210585 0.888338,0.210435 0.888349
0.683401 0.844613,0.502608 0.917517
0.670661 0.69692,0.683401 0.844613
0.800249 0.835478,0.683401 0.844613
0.434114 0.12756,0.433321 0.119626
0.127187 0.615504,0.433321 0.119626
0.814059 0.857826,0.990225 0.9172
0.829755 0.546772,0.434114 0.12756
0.439592 0.181246,0.434114 0.12756
0.801346 0.835391,0.800249 0.835478
0.823062 0.791469,1 0.726024
0.670661 0.69692,0.801346 0.835391
0.802831 0.835271,0.801346 0.835391
0.439592 0.181246,0.670661 0.69692
0.823062 0.791469,0.829755 0.546772
0.802831 0.835271,0.814059 0.857826
0.823062 0.791469,0.802831 0.835271

image

Environment

diablohsh commented 5 years ago

using the same vertex data from the about code, Segment_Delaunay_graph_2 never return from draw_dual function.

Source Code


#include <sstream>
#include <vector>

#define CGAL_HEADER_ONLY 1

#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
#include <CGAL/Segment_Delaunay_graph_filtered_traits_2.h>
#include <CGAL/Segment_Delaunay_graph_traits_2.h>
#include <CGAL/Segment_Delaunay_graph_2.h>

using namespace std;

typedef CGAL::Exact_predicates_exact_constructions_kernel               K;
typedef CGAL::Segment_Delaunay_graph_filtered_traits_2<K>               SDGT;
typedef CGAL::Segment_Delaunay_graph_2<SDGT>                            SDG;

typedef K::FT           Num_t;
typedef K::Point_2      Point_2;
typedef SDGT::Site_2    Site_2;

struct Stream_t
{
    Stream_t& operator<<(const K::Segment_2& e)
    {
        auto s = e.source();
        auto t = e.target();
        cout << s << "," << t << endl;
        return *this;
    }
    template<typename T>Stream_t& operator<<(const T&r) { return *this; }
};

int main()
{
    istringstream strm("2ftkqnron5d/2gosa7pa2gw 14o7lsvsnhx/18ce53un18g 24du0d5436p/2gosa7pa2gw oiv9ts2otv/m672jxbim8 1pj1zmqn3g3/2gosa7pa2gw 2g0t7s4y9vh/2gosa7pa2gw 1qvppf3jptx/2gosa7pa2gw p35x4z8lkr/m672jxbim8 1pmtfy222kr/4xdkkfek4xs 1fxltrt3fcj/18ce53un18g 1mxi0dc8tg3/4xdkkfek4xs bfx5ch9ua7/b33j9ynrb4 1iy0p7mruhd/9ur54ut49vk 1gfiq7inzsz/18ce53un18g 0 1 7v13lg71tl/jpia9pm8jr4 8whunyntav/b33j9ynrb4 194brmj62m5/9ur54ut49vk ral4ou2do3/18ce53un18g 3fj4l6xlwr/b33j9ynrb4 g557q50kj3/m672jxbim8 tqmz8g6vtn/18ce53un18g fga1rb5hd5/m672jxbim8 4lkriq10u1/5jjrmzbvnk 1chq45rxz7j/2gosa7pa2gw 1 40vl5jmjv1/5jjrmzbvnk");

    vector<Point_2> vecPoint;

    string line;
    while (true)
    {
        mpq_t x, y;
        mpq_init(x);
        mpq_init(y);

        if (!getline(strm, line, ' '))
            break;
        mpq_set_str(x, line.c_str(), 36);

        getline(strm, line, ' ');
        mpq_set_str(y, line.c_str(), 36);

        vecPoint.emplace_back(Num_t(x), Num_t(y));
    }

    vector<Site_2> vecSite;

    for (size_t i = 1; i < vecPoint.size(); i++)
        vecSite.push_back(Site_2::construct_site_2(vecPoint[i - 1], vecPoint[i]));

    vecSite.push_back(Site_2::construct_site_2(vecPoint.back(), vecPoint.front()));

    SDG sdg;

    sdg.insert(vecSite.begin(), vecSite.end());
    sdg.draw_dual(Stream_t());

    return 0;
}
diablohsh commented 5 years ago

After some digestion into the code, i narrow down the cause of the problem.there is probably a failure on parallel detection.

  if ( parallel || c_ == cq_ ) {

of CGAL\Segment_Delaunay_graph_2\Voronoi_vertex_sqrt_field_new_C2.h:228

where c and cq have extremly close value c_: _inf=0.050863068750012563 sup=0.050863068750012841 cq :_inf=0.050863068750011564 _sup=0.050863068750011842 but still failed on the equivalence predication. after i change the code into:

if ( parallel || CGAL::abs(c_ - cq_)<0.000001 ) {

The problem is gone,but i don't think of it a good modification. any suggestion for a fast fix? i'd rather not to wait for the next release.

purple line of the image is the Segment Delaunay Graph edge.

vdbug

diablohsh commented 5 years ago

futher information, i have tried all kinds of traits, including:

Segment_Delaunay_graph_traits_without_intersections_2<K, Integral_domain_without_division_tag> Segment_Delaunay_graph_traits_without_intersections_2<K, Field_with_sqrt_tag> Segment_Delaunay_graph_traits_without_intersections_2<K, Field_with_sqrt_tag> with CGAL_SDG_USE_OLD_INCIRCLE compile definition

and K defined to different kinds of kernels including:

CGAL::Simple_cartesian CGAL::Simple_cartesian CGAL::Simple_cartesian CGAL::Simple_cartesian<CGAL::Lazy_exact_nt> CGAL::Simple_cartesian CGAL::Simple_cartesian< CGAL::Lazy_exact_nt<CGAL::Quotient>> CGAL::Cartesian<CGAL::Quotient > CGAL::Exact_predicates_inexact_constructions_kernel

computation will go down to three class:

Voronoi_vertex_ring_C2 Voronoi_vertex_sqrt_field_C2 Voronoi_vertex_sqrt_field_new_C2

neither of them can get correct result.

diablohsh commented 5 years ago

Hi all, with the following kernel and traits, the problem is gone.

typedef CGAL::Exact_predicates_exact_constructions_kernel_with_sqrt K; typedef CGAL::Segment_Delaunay_graph_traits_2<K,CGAL::Field_with_sqrt_tag> SDGT; typedef CGAL::Segment_Delaunay_graph_2 SDG; typedef CGAL::Segment_Delaunay_graph_adaptation_traits_2 SDGAT; typedef CGAL::Segment_Delaunay_graph_degeneracy_removal_policy_2 SDGAP; typedef CGAL::Voronoi_diagram_2<SDG, SDGAT, SDGAP> VD;