Union of polygons with segments and circular arcs

five.txt

-2 2 -2 14 -4 14 -4 2
-2 11 8 11 8 13 -2 13
0 2 12 2 12 8 0 8
2 9 5 12 3 14 0 11
9 0 9 12 7 12 7 0

union.cpp

#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
#include <CGAL/Gps_circle_segment_traits_2.h>
#include <CGAL/Boolean_set_operations_2.h>
#include <CGAL/Timer.h>

#include <list>
#include <cstdlib>
#include <string>

typedef CGAL::Exact_predicates_exact_constructions_kernel Kernel;
typedef Kernel::FT                                        FT;
typedef Kernel::Point_2                                   Point_2;
typedef Kernel::Circle_2                                  Circle_2;

typedef CGAL::Gps_circle_segment_traits_2<Kernel>         Traits_2;
typedef typename Traits_2::Point_2                        Arc_point_2;
typedef typename Arc_point_2::CoordNT                     CoordNT;
typedef typename Traits_2::Curve_2                        Curve_2;
typedef typename Traits_2::X_monotone_curve_2             X_monotone_curve_2;
typedef Traits_2::Polygon_2                               Polygon_2;
typedef Polygon_2::Curve_const_iterator                   Curve_const_iterator;
typedef Traits_2::Polygon_with_holes_2                    Polygon_with_holes_2;

typedef CGAL::General_polygon_set_2<Traits_2>             Polygon_set_2;




void construct_arc(const Point_2& ps, const Point_2& pt, Polygon_2& pgn)
{
  Traits_2 traits;
  Traits_2::Make_x_monotone_2 make_x_monotone = traits.make_x_monotone_2_object();
  CGAL::Object obj_vec[3];
  CGAL::Object *obj_begin = (obj_vec + 0);
  CGAL::Object *obj_end;
  X_monotone_curve_2 cv1;
  Circle_2 supp_circ;
  
  const FT x_coord = ((ps.x() + pt.x())/2);
  const FT y_coord = ((ps.y() + pt.y())/2);
  const FT sqr_rad = CGAL::squared_distance(ps, pt)  / 4; 
  
  supp_circ = Circle_2(Point_2(x_coord, y_coord), sqr_rad,
                       CGAL::COUNTERCLOCKWISE);
  
  Curve_2 circ_arc(supp_circ, 
                   Arc_point_2(ps.x(), ps.y()),
                   Arc_point_2(pt.x(), pt.y()));
  
  // Break the arc into x-monotone subarcs (there can be at most
  // three subarcs) and add them to the polygon.
  obj_end = make_x_monotone(circ_arc, obj_begin);
  std::size_t n_subarcs = (obj_end - obj_begin);
  
  for (std::size_t i = 0; i < n_subarcs; i++)
    {
      if (CGAL::assign(cv1, obj_vec[i]))
        pgn.push_back(cv1);
    }
}


// Constructs a polygon with round caps from p1 to p2, and from p3 to p0
Polygon_2 construct_polygon (const Point_2& p0, const Point_2& p1,
                             const Point_2& p2, const Point_2& p3)
{
  Polygon_2 pgn;
  X_monotone_curve_2 s1(p0, p1);
  pgn.push_back(s1);
  construct_arc(p1, p2, pgn);
  X_monotone_curve_2 s3(p2, p3);
  pgn.push_back(s3);
  construct_arc(p3, p0, pgn);
  std::cout << pgn.size() << std::endl;
  return pgn;
}


void print(const Polygon_2& gp)
{
  Arc_point_2 first = gp.curves_begin()->source();
  for(Curve_const_iterator cci = gp.curves_begin(), end = gp.curves_end();
      cci != end;
      ++cci){
    X_monotone_curve_2 curve = *cci;
    std::string ori = (curve.orientation() == CGAL::COLLINEAR)?"seg"
      : (curve.orientation() == CGAL::NEGATIVE)? "cw":"ccw";
    
    std::cout << curve.source();
    CoordNT x = curve.source().x();
    if(x.is_extended()){
      std::cout << " (x: " << x.a0() << " "<< x.a1() << " "<< x.root() << ")";
    }
    CoordNT y = curve.source().y();
    if(y.is_extended()){
      std::cout << " (y: " << y.a0() << " "<< y.a1() << " "<< y.root() << ")";
      
    }
    std::cout << " -" << ori << "-> " << std::flush; 
  }
  std::cout << first;
}


void print(const Polygon_with_holes_2& p)
{
  const Polygon_2& gp =  p.outer_boundary();
  std::cout << "Outer boundary:\n";
  print(gp);
  std::cout << std::endl;
  for(auto it = p.holes_begin(); it != p.holes_end(); ++it){
    std::cout << "Hole:\n";
    print(*it);
  }
}


// The main program:
int main (int argc, char* argv[])
{
  std::list<Polygon_2> polygons;
#if 1
  std::ifstream in(argv[1]);
  Point_2 p,q,r,s;
  while(in >> p >> q >> r >> s){
    polygons.push_back (construct_polygon (p,q,r,s));
  }
#else  

  unsigned int k;
  for (k = 0; k < 2; k++) {
    Point_2 ll(0, k-0.1);
    Point_2 lr(1000, k-0.2);
    Point_2 tr(1000, k+0.2);
    Point_2 tl(0, k+0.1);
    polygons.push_back (construct_polygon (ll,lr,tr,tl));
  }
  
  for (k = 0; k < 2; k++) {
    Point_2 ll(k-0.1, 0);
    Point_2 lr(k+0.1, 0);
    Point_2 tr(k+0.1, 1000);
    Point_2 tl(k-0.1, 1000);
    polygons.push_back (construct_polygon (lr,tr,tl, ll));

  }
#endif 
  CGAL::Timer t;
  t.start();
  // Compute the union aggregately.
  std::list<Polygon_with_holes_2> res;

  std::cout << "Join " << polygons.size() << " polygons"<< std::endl;
  CGAL::join (polygons.begin(), polygons.end(), std::back_inserter (res));

  std::cout << t.time() << " sec." << std::endl;
  t.stop();
  // Print the result.
  std::cout << res.size() << " polygons with holes" << std::endl;
  
  for(const Polygon_with_holes_2& p : res){
    print(p);
  }
  return 0;
}