How to convert a voronoi diagram cropped by a rectangle (boundaing box) to a set of polygons?

Question

I am very new to CGAL and recently I have been trying to use CGAL to compute Voronoi diagram cropped by a bounding box (rectangle) in a C++ code and I was successful with that. I took advantage of an available example code in CGAL documentation. Here is the code:

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Delaunay_triangulation_2.h>
#include <iterator>
typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
typedef K::Point_2 Point_2;
typedef K::Iso_rectangle_2 Iso_rectangle_2;
typedef K::Segment_2 Segment_2;
typedef K::Ray_2 Ray_2;
typedef K::Line_2 Line_2;
typedef CGAL::Delaunay_triangulation_2<K> Delaunay_triangulation_2;
//A class to recover Voronoi diagram from stream.
//Rays, lines and segments are cropped to a rectangle
//so that only segments are stored
struct Cropped_voronoi_from_delaunay{
  std::list<Segment_2> m_cropped_vd;
  Iso_rectangle_2 m_bbox;
  Cropped_voronoi_from_delaunay(const Iso_rectangle_2& bbox):m_bbox(bbox){}
  template <class RSL>
  void crop_and_extract_segment(const RSL& rsl){
    CGAL::Object obj = CGAL::intersection(rsl,m_bbox);
    const Segment_2* s=CGAL::object_cast<Segment_2>(&obj);
    if (s) m_cropped_vd.push_back(*s);
  }
  void operator<<(const Ray_2& ray) { crop_and_extract_segment(ray); }
  void operator<<(const Line_2& line) { crop_and_extract_segment(line); }
  void operator<<(const Segment_2& seg){ crop_and_extract_segment(seg); }
};
int main(){
  //consider some points
  std::vector<Point_2> points;
  points.push_back(Point_2(0,0));
  points.push_back(Point_2(1,1));
  points.push_back(Point_2(0,1));
  Delaunay_triangulation_2 dt2;
  //insert points into the triangulation
  dt2.insert(points.begin(),points.end());
  //construct a rectangle
  Iso_rectangle_2 bbox(-1,-1,2,2);
  Cropped_voronoi_from_delaunay vor(bbox);
  //extract the cropped Voronoi diagram
  dt2.draw_dual(vor);
  //print the cropped Voronoi diagram as segments
  std::copy(vor.m_cropped_vd.begin(),vor.m_cropped_vd.end(),
            std::ostream_iterator<Segment_2>(std::cout,"\n"));
}

Now I intend to generate the voronoi faces and convert them to polygons in order to use CGAL::intersection boolean operations on polygons. A similar question have been previously asked but no CGAL solution was provided. Two sets of polygons need to be considered; first a set of complete voronoi cells within the bounding box that have no intersections with the cropping rectangle. Second set would consist of voronoi cells that are actually clipped by the bounding box. Any comments or hints would be really appreciated.

Answer 1

The following will generate a random point cloud, find its Voronoi diagram, crop that diagram to the cloud's bounding box, and generate well-known text polygons.

//Finds the cropped Voronoi diagram of a set of points and saves it as WKT
//Compile with: g++ main.cpp -Wall -lCGAL -lgmp
//Author: Richard Barnes (rbarnes.org)
#include <CGAL/Exact_predicates_exact_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/Boolean_set_operations_2.h>
#include <CGAL/bounding_box.h>
#include <CGAL/Polygon_2.h>
#include <iostream>
#include <cstdint>

//Used to convert otherwise infinite rays into looooong line segments
const int RAY_LENGTH = 1000;

typedef CGAL::Exact_predicates_exact_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)); 

  //Find the bounding box of the points. This will be used to crop the Voronoi
  //diagram later.
  const K::Iso_rectangle_2 bbox = CGAL::bounding_box(wpoints.begin(), wpoints.end());

  //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+RAY_LENGTH*dir.dx(),dsy+RAY_LENGTH*dir.dy());
      if(outgoing)
        return K::Segment_2(
          dray->source(),
          tpoint
        );
      else
        return K::Segment_2(
          tpoint,
          dray->source()
        );
    }
  };

  //First line of WKT CSV output
  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();

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

    //In order to crop the Voronoi diagram, we need to convert the bounding box
    //into a polygon. You'd think there'd be an easy way to do this. But there
    //isn't (or I haven't found it).
    CGAL::Polygon_2<K> bpoly;
    bpoly.push_back(K::Point_2(bbox.xmin(),bbox.ymin()));
    bpoly.push_back(K::Point_2(bbox.xmax(),bbox.ymin()));
    bpoly.push_back(K::Point_2(bbox.xmax(),bbox.ymax()));
    bpoly.push_back(K::Point_2(bbox.xmin(),bbox.ymax()));

    //Perform the intersection. Since CGAL is very general, it believes the
    //result might be multiple polygons with holes.
    std::list<CGAL::Polygon_with_holes_2<K>> isect;
    CGAL::intersection(pgon, bpoly, std::back_inserter(isect));

    //But we know better. The intersection of a convex polygon and a box is
    //always a single polygon without holes. Let's assert this.
    assert(isect.size()==1);

    //And recover the polygon of interest
    auto &poly_w_holes = isect.front();
    auto &poly_outer   = poly_w_holes.outer_boundary();

    //Print the polygon as a WKT polygon
    std::cout<<fnum<<", "
    "\"POLYGON ((";
    for(auto v=poly_outer.vertices_begin();v!=poly_outer.vertices_end();v++)
      std::cout<<v->x()<<" "<<v->y()<<", ";
    std::cout<<poly_outer.vertices_begin()->x()<<" "<<poly_outer.vertices_begin()->y()<<"))\"\n";
  }

  return 0;
}

Answer 2

There is some experimental code here: http://code.google.com/p/cgal-voronoi-cropping that crop a voronoi diagram to a rectangle, the result being a HDS. See main.cpp in the test directory.