Search code examples
c++boostcomputational-geometrycgalvoronoi

How to print the faces of a Voronoi diagram?


Code bellow assumes input is points, not line segments (which is wrong).


Following this 2D Voronoi Diagram Adaptor example, I am trying to write a program which takes as input line segments and prints the vertices of the faces of the Voronoi diagram.

Here is my attempt (keeping includes/typedefs of the example):

// standard includes
#include <iostream>
#include <fstream>
#include <cassert>
// includes for defining the Voronoi diagram adaptor
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Delaunay_triangulation_2.h>
#include <CGAL/Voronoi_diagram_2.h>
#include <CGAL/Delaunay_triangulation_adaptation_traits_2.h>
#include <CGAL/Delaunay_triangulation_adaptation_policies_2.h>
// typedefs for defining the adaptor
typedef CGAL::Exact_predicates_inexact_constructions_kernel                  K;
typedef CGAL::Delaunay_triangulation_2<K>                                    DT;
typedef CGAL::Delaunay_triangulation_adaptation_traits_2<DT>                 AT;
typedef CGAL::Delaunay_triangulation_caching_degeneracy_removal_policy_2<DT> AP;
typedef CGAL::Voronoi_diagram_2<DT,AT,AP>                                    VD;
// typedef for the result type of the point location
typedef AT::Site_2                    Site_2;
typedef AT::Point_2                   Point_2;
typedef VD::Locate_result             Locate_result;
typedef VD::Vertex_handle             Vertex_handle;
typedef VD::Face_handle               Face_handle;
typedef VD::Halfedge_handle           Halfedge_handle;
typedef VD::Ccb_halfedge_circulator   Ccb_halfedge_circulator;

int main()
{
    std::ifstream ifs("data.cin");
    assert( ifs );
    VD vd;
    Site_2 t;
    while ( ifs >> t ) { vd.insert(t); }
    ifs.close();
    assert( vd.is_valid() );
    Face_handle* f = boost::get<Face_handle>(vd);
    std::cout << "Exiting...\n";
    return 0;
}

This receives a compilation error:

/home/gsamaras/CGAL-4.7/code/voronoi_adaptor/voronoi_adaptor.cpp:46:48: error: no matching function for call to ‘get(VD&)’
     Face_handle* f = boost::get<Face_handle>(vd);
                                                ^
/home/gsamaras/CGAL-4.7/code/voronoi_adaptor/voronoi_adaptor.cpp:46:48: note: candidates are:
In file included from /usr/include/boost/variant.hpp:22:0,
                 from /home/gsamaras/CGAL-4.7/code/voronoi_adaptor/../../include/CGAL/Object.h:36,
                 from /home/gsamaras/CGAL-4.7/code/voronoi_adaptor/../../include/CGAL/kernel_basic.h:33,
                 from /home/gsamaras/CGAL-4.7/code/voronoi_adaptor/../../include/CGAL/basic.h:46,
                 from /home/gsamaras/CGAL-4.7/code/voronoi_adaptor/../../include/CGAL/Cartesian/Cartesian_base.h:28,
                 from /home/gsamaras/CGAL-4.7/code/voronoi_adaptor/../../include/CGAL/Simple_cartesian.h:28,
                 from /home/gsamaras/CGAL-4.7/code/voronoi_adaptor/../../include/CGAL/Exact_predicates_inexact_constructions_kernel.h:28,
                 from /home/gsamaras/CGAL-4.7/code/voronoi_adaptor/voronoi_adaptor.cpp:6:
/usr/include/boost/variant/get.hpp:141:1: note: template<class U, class T0, class T1, class T2, class T3, class T4, class T5, class T6, class T7, class T8, class T9, class T10, class T11, class T12, class T13, class T14, class T15, class T16, class T17, class T18, class T19> typename boost::add_pointer<T>::type boost::get(boost::variant<T0, T1, T2, T3, T4, T5, T6, T7, T8, T9, T10, T11, T12, T13, T14, T15, T16, T17, T18, T19>*)
 get(
 ^
/usr/include/boost/variant/get.hpp:141:1: note:   template argument deduction/substitution failed:
/home/gsamaras/CGAL-4.7/code/voronoi_adaptor/voronoi_adaptor.cpp:46:48: note:   mismatched types ‘boost::variant<T0, T1, T2, T3, T4, T5, T6, T7, T8, T9, T10, T11, T12, T13, T14, T15, T16, T17, T18, T19>*’ and ‘CGAL::Voronoi_diagram_2<CGAL::Delaunay_triangulation_2<CGAL::Epick>, CGAL::Delaunay_triangulation_adaptation_traits_2<CGAL::Delaunay_triangulation_2<CGAL::Epick> >, CGAL::Delaunay_triangulation_caching_degeneracy_removal_policy_2<CGAL::Delaunay_triangulation_2<CGAL::Epick> > >’
     Face_handle* f = boost::get<Face_handle>(vd);
                                                ^
...

Solution

  • I needed to print the faces of the Voronoi diagram as well-known text polygons. This entailed iterating over the edges of the faces and giving finite representations to the infinite points. I did so as follows.

    //Generate WKT polygons from Voronoi cells using CGAL
    //Compile with: g++ main.cpp -Wall -lCGAL -lgmp
    //Author: Richard Barnes (rbarnes.org)
    #include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
    #include <CGAL/Regular_triangulation_filtered_traits_2.h>
    #include <CGAL/Regular_triangulation_adaptation_traits_2.h>
    #include <CGAL/Regular_triangulation_adaptation_policies_2.h>
    #include <CGAL/Regular_triangulation_2.h>
    #include <CGAL/Voronoi_diagram_2.h>
    #include <CGAL/intersections.h>
    #include <CGAL/bounding_box.h>
    #include <CGAL/Polygon_2.h>
    #include <iostream>
    #include <cstdint>
    
    typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
    typedef CGAL::Regular_triangulation_filtered_traits_2<K>  Traits;
    
    typedef CGAL::Regular_triangulation_2<Traits> RT2;
    typedef CGAL::Regular_triangulation_adaptation_traits_2<RT2>         AT;
    typedef CGAL::Regular_triangulation_degeneracy_removal_policy_2<RT2> DRP;
    typedef CGAL::Voronoi_diagram_2<RT2, AT, DRP> VD;
    
    int main(int argc, char **argv){
      std::vector<RT2::Weighted_point> wpoints;
    
      std::cout.precision(4);
      std::cout.setf(std::ios::fixed);
    
      //Generated random points
      for(int i=0;i<100;i++)
        //Weight of 0 gives a Voronoi diagram. Non-zero weight gives a power diagram
        wpoints.push_back(RT2::Weighted_point(K::Point_2(rand()%100,rand()%100), 0)); 
    
      //Create a Regular Triangulation from the points
      RT2 rt(wpoints.begin(), wpoints.end());
      rt.is_valid();
    
      //Wrap the triangulation with a Voronoi diagram adaptor. This is necessary to
      //get the Voronoi faces.
      VD vd(rt);
    
      //CGAL often returns objects that are either segments or rays. This converts
      //these objects into segments. If the object would have resolved into a ray,
      //that ray is intersected with the bounding box defined above and returned as
      //a segment.
      const auto ConvertToSeg = [&](const CGAL::Object seg_obj, bool outgoing) -> K::Segment_2 {
        //One of these will succeed and one will have a NULL pointer
        const K::Segment_2 *dseg = CGAL::object_cast<K::Segment_2>(&seg_obj);
        const K::Ray_2     *dray = CGAL::object_cast<K::Ray_2>(&seg_obj);
        if (dseg) { //Okay, we have a segment
          return *dseg;
        } else {    //Must be a ray
          const auto &source = dray->source();
          const auto dsx     = source.x();
          const auto dsy     = source.y();
          const auto &dir    = dray->direction();
          const auto tpoint  = K::Point_2(dsx+1000*dir.dx(),dsy+1000*dir.dy());
          if(outgoing)
            return K::Segment_2(
              dray->source(),
              tpoint
            );
          else
            return K::Segment_2(
              tpoint,
              dray->source()
            );
        }
      };
    
      std::cout<<"\"id\",\"geom\"\n";
    
      int fnum = 0;
      //Loop over the faces of the Voronoi diagram in some arbitrary order
      for(VD::Face_iterator fit = vd.faces_begin(); fit!=vd.faces_end();++fit,fnum++){
        CGAL::Polygon_2<K> pgon;
    
        //Edge circulators traverse endlessly around a face. Make a note of the
        //starting point so we know when to quit.
        VD::Face::Ccb_halfedge_circulator ec_start = fit->ccb();
    
        //Find a bounded edge to start on
        for(;ec_start->is_unbounded();ec_start++){}
    
        //Current location of the edge circulator
        VD::Face::Ccb_halfedge_circulator ec = ec_start;
    
        do {
          //A half edge circulator representing a ray doesn't carry direction
          //information. To get it, we take the dual of the dual of the half-edge.
          //The dual of a half-edge circulator is the edge of a Delaunay triangle.
          //The dual of the edge of Delaunay triangle is either a segment or a ray.
          // const CGAL::Object seg_dual = rt.dual(ec->dual());
          const CGAL::Object seg_dual = vd.dual().dual(ec->dual());
    
          //Convert the segment/ray into a segment
          const auto this_seg = ConvertToSeg(seg_dual, ec->has_target());
    
          pgon.push_back(this_seg.source());
    
          //If the segment has no target, it's a ray. This means that the next
          //segment will also be a ray. We need to connect those two rays with a
          //segment. The following accomplishes this.
          if(!ec->has_target()){
            const CGAL::Object nseg_dual = vd.dual().dual(ec->next()->dual());
            const auto next_seg = ConvertToSeg(nseg_dual, ec->next()->has_target());
            pgon.push_back(next_seg.target());
          }
        } while ( ++ec != ec_start ); //Loop until we get back to the beginning
    
        std::cout<<fnum<<", "
        "\"POLYGON ((";
        for(auto v=pgon.vertices_begin();v!=pgon.vertices_end();v++)
          std::cout<<v->x()<<" "<<v->y()<<", ";
        std::cout<<pgon.vertices_begin()->x()<<" "<<pgon.vertices_begin()->y()<<"))\"\n";
      }
    
      return 0;
    }