Как преобразовать диаграмму Вороного, обрезанную прямоугольником (граничной рамкой), в набор многоугольников?

Здесь есть экспериментальный код: http://code.google.com/p/cgal-voronoi-cropping, которые обрезают диаграмму Вороного до прямоугольника, в результате чего получается HDS. См. main.cpp в тестовом каталоге.

Следующее будет генерировать случайное облако точек, находить его диаграмму Вороного, обрезать эту диаграмму до ограничивающей рамки облака и генерировать хорошо известные текстовые полигоны.

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