CGAL users discussion list

[Shriramana Sharma <>]

Hello.

My previous mails have referred to my work on using CGAL for boolops
on bezier polygons and the difficulty in extracting the result of
those boolops. While intersection of beziers cannot always be exactly
computed, a good approximation to machine precision (actually it is
hardcoded to be around 2^-53) can indeed be computed, but I was able
to arrive at the actual code required for this only by myself diving
into the source code (the examples/API documentation being
insufficient).

To spare anyone else the effort, I'd hence like to contribute a
program for the CGAL examples. It is based on the existing
bezier_traits_adapter2 example but goes further and actually displays
the control points of the bezier polygons produced by the boolop
rather than just stopping at saying how many objects were produced.

Note that the file uses the de Casteljau algorithm from a separate
header file which I adapted from the internal CGAL private header to a
format which should hopefully make it usable in general with CGAL.

My contribution is under the GPLv3.

-- 
Shriramana Sharma
/*
Copyright (c) 2006,2007,2009,2010,2011 Tel-Aviv University (Israel).
All rights reserved.

This file is part of CGAL (www.cgal.org).

You can redistribute it and/or modify it under the terms of the GNU
General Public License as published by the Free Software Foundation,
either version 3 of the License, or (at your option) any later version.

This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.

Author(s)     : Ron Wein     <>
                Iddo Hanniel <>

This file has been adapted from the CGAL file
<CGAL_src_tree_root>/include/CGAL/Arr_geometry_traits/de_Casteljau_2.h
under the above licence.

Author: Shriramana Sharma <samjnaa at gmail dot com>
Copyright (c) 2013 all authors mentioned above
*/

# ifndef DE_CASTELJAU_H
# define DE_CASTELJAU_H

# include <vector>
# include <deque>
# include <iterator>
# include <stdexcept>

// input: a curve's control points given as start and end iterators, a parameter value t
// output: the control points of the first and second pieces when the curve is split at t
// return value: the point corresponding to the parameter t of the curve

template < class InputIterator, class BackInserter, class FrontInserter >
typename InputIterator :: value_type
splitBezier (
	InputIterator ctrlPtsBegin, InputIterator ctrlPtsEnd,
	const typename InputIterator :: value_type :: R :: FT &  t,
	// value_type will be a rational point, whose R member type is the rational kernel and whose FT member type is the number type
	BackInserter firstPieceBack, FrontInserter secondPieceFront )
{
	typedef typename InputIterator :: value_type  point ;
	typedef typename point :: R :: FT             numType ;

	const unsigned int n_pts = std :: distance ( ctrlPtsBegin, ctrlPtsEnd ) ;
	if ( n_pts == 0 ) throw std :: invalid_argument ( "The input curve is empty." ) ;
	
	std :: vector < point >  v ;
	std :: copy ( ctrlPtsBegin, ctrlPtsEnd, std :: back_inserter ( v ) ) ;

	const numType one_t = numType ( 1 ) - t ;
	unsigned int lastIndex = n_pts - 1 ;

	# define PUSH_FIRST_LAST_POINTS \
	* firstPieceBack   ++ = v [ 0         ] ; \
	* secondPieceFront ++ = v [ lastIndex ] ;

	PUSH_FIRST_LAST_POINTS
	while ( lastIndex > 0 )
	{
		// de Casteljau: The next level's i-th point is given by: (1 - t)*v[i] + t*v[i + 1]
		for ( unsigned int i = 0; i < lastIndex ; ++ i )
			v[i] = point ( one_t * v [ i ] . x () + t * v [ i + 1 ] . x (),
			               one_t * v [ i ] . y () + t * v [ i + 1 ] . y () ) ;

		-- lastIndex ;
		PUSH_FIRST_LAST_POINTS
	}
	return v [ 0 ] ; // this is the point on the curve at parameter t
}

// input: a curve's control points given as start and end iterators, parameter values t0 and t1
// output: the control points of the segment of the curve between t0 and t1

template < class InputIterator, class OutputIterator >
void
subcurveBezier (
	InputIterator ctrlPtsBegin, InputIterator ctrlPtsEnd,
	const typename InputIterator :: value_type :: R :: FT &  t0,
	const typename InputIterator :: value_type :: R :: FT &  t1,
	// value_type will be a rational point, whose R member type is the rational kernel and whose FT member type is the number type
	OutputIterator subcurve )
{
	typedef typename InputIterator :: value_type  point ;
	typedef typename point :: R :: FT             numType ;

	const unsigned int n_pts = std :: distance ( ctrlPtsBegin, ctrlPtsEnd ) ;
	if ( n_pts == 0 ) throw std :: invalid_argument ( "The input curve is empty." ) ;

	std :: deque < point >  p, q, r ; // temporaries

	bool swapped = false ;
	numType u0, u1 ;
	if ( t0 > t1 )
	{
		u0 = t1 ; u1 = t0 ;
		swapped = true ;
	}
	else
	{
		u0 = t0 ; u1 = t1 ;
	}

	splitBezier ( ctrlPtsBegin, ctrlPtsEnd,        u0               , std :: back_inserter ( p ), std :: front_inserter ( q ) ) ;
	splitBezier ( q . begin (), q . end (), ( u1 - u0 ) / ( 1 - u0 ), std :: back_inserter ( r ), std :: front_inserter ( p ) ) ;
	// note: contents of p now are not really meaningful but it is only a temporary anyway

	if ( swapped ) std :: reverse_copy ( r . begin (), r . end (), subcurve ) ;
	else           std ::         copy ( r . begin (), r . end (), subcurve ) ;
}

# endif // DE_CASTELJAU_H

/*! \file bezier_polygon_wh_output.cpp
 * Shows how to extract and print the boundaries and holes of a bezier polygon
 */

#include <CGAL/basic.h>

#ifndef CGAL_USE_CORE
#include <iostream>
int main ()
{
	std::cout << "Sorry, this example needs CORE ..." << std::endl;
	return 1 ;
}
#else

#include <CGAL/Cartesian.h>
#include <CGAL/CORE_algebraic_number_traits.h>
#include <CGAL/Arr_Bezier_curve_traits_2.h>
#include <CGAL/Gps_traits_2.h>
#include <CGAL/Boolean_set_operations_2.h>
#include <CGAL/General_polygon_set_2.h>
#include <CGAL/Timer.h>

#include <fstream>
#define NEWLINE std::cout << std::endl
#define APPROX_EPSILON 2e-16 // to be exact it is 2^-52

typedef CGAL::CORE_algebraic_number_traits              Nt_traits;
typedef Nt_traits::Rational                             Rational;
typedef Nt_traits::Algebraic                            Algebraic;
typedef CGAL::Cartesian<Rational>                       Rat_kernel;
typedef CGAL::Cartesian<Algebraic>                      Alg_kernel;
typedef Rat_kernel::Point_2                             Rat_point;
typedef CGAL::Arr_Bezier_curve_traits_2<Rat_kernel, Alg_kernel, Nt_traits>  Bezier_traits;
typedef Bezier_traits::Point_2                          Bezier_alg_point;
typedef Bezier_traits::Curve_2                          Bezier_curve;
typedef Bezier_traits::X_monotone_curve_2               Bezier_X_monotone_curve;
typedef CGAL::Gps_traits_2<Bezier_traits>               Bezier_polygon_traits;
typedef Bezier_polygon_traits::General_polygon_2                Bezier_polygon;
typedef Bezier_polygon_traits::General_polygon_with_holes_2     Bezier_polygon_with_holes;
typedef CGAL::General_polygon_set_2<Bezier_polygon_traits>      Bezier_polygon_set;
typedef std::vector<Bezier_polygon>                     Bezier_polygon_vector;

CGAL::Timer g_refineTimer ;
unsigned int g_totalRefinings = 0 ;

Bezier_curve read_bezier_curve ( std::istream & is, bool format_is_double )
{
	// Read the number of control points.
	unsigned int n;
	is >> n;

	// Read the control points.
	std::vector<Rat_point> ctrl_pts;
	for ( unsigned int k = 0; k < n; k++ )
	{
		Rat_point p ;
		if ( format_is_double )
		{
			double x, y ;
			is >> x >> y ;
			p = Rat_point ( x, y );
		}
		else is >> p ;
		if ( k == 0 || ctrl_pts[k - 1] != p ) ctrl_pts.push_back ( p ) ;
	}
	
	return Bezier_curve ( ctrl_pts.begin(), ctrl_pts.end() );
}

bool read_bezier_polygon ( const char * filename, Bezier_polygon_set & S )
{
	bool success = false ;
	std::ifstream in_file ( filename );

	if ( not in_file ) return success ;

	try
	{
		std::cout << "Reading " << filename << std::endl ;
		std::string format ;
		std::getline ( in_file, format );
		bool format_is_double = ( format.length() >= 6 && format.substr ( 0, 6 ) == "DOUBLE" ) ;

		// Read the number of bezier polygon with holes (PWH)
		unsigned int n_pwhs ;
		in_file >> n_pwhs;

		for ( unsigned int r = 0 ; r < n_pwhs ; ++ r )
		{
			// Read the number of polygons (boundary with any holes) in this PWH
			unsigned int n_polygons;
			in_file >> n_polygons;

			// First we process each polygon and add it to a set of polygons
			Bezier_polygon_vector polygonSet ;
			for ( unsigned int b = 0 ; b < n_polygons ; ++ b )
			{
				// Read the number of bezier curves in this polygon
				unsigned int n_curves;
				in_file >> n_curves;

				// Process each curve in this polygon

				std::list<Bezier_X_monotone_curve> monotoneCurveSet;

				for ( unsigned int k = 0; k < n_curves; ++ k )
				{
					// Read the current curve and subdivide it into x-monotone subcurves.

					std::list<CGAL::Object>                  cgalObjects;
					std::list<CGAL::Object>::const_iterator  cgalObjIter;
					Bezier_X_monotone_curve                  monotoneCurve;
					Bezier_polygon_traits                    traits;
					Bezier_polygon_traits::Make_x_monotone_2 makeMonotone = traits.make_x_monotone_2_object();

					Bezier_curve curve = read_bezier_curve ( in_file, format_is_double );
					if ( curve.number_of_control_points() >= 2 )
					{
						makeMonotone ( curve, std::back_inserter ( cgalObjects ) );
						for ( cgalObjIter = cgalObjects.begin(); cgalObjIter != cgalObjects.end(); ++cgalObjIter )
						{
							if ( CGAL::assign ( monotoneCurve, *cgalObjIter ) ) monotoneCurveSet.push_back ( monotoneCurve );
						}
					}
				} // next curve in this polygon

				Bezier_polygon polygon ( monotoneCurveSet.begin(), monotoneCurveSet.end() );
				CGAL::Orientation  orient = polygon.orientation();
				std::cout << "Polygon with holes #" << r + 1 << " Polygon #" << b + 1 << std::endl ;
				std::cout << "  Orientation: " << orient ;

				if ( ( b == 0 && orient == CGAL::CLOCKWISE ) || ( b > 0 && orient == CGAL::COUNTERCLOCKWISE ) )
				{
					std::cout << " which is incorrect. Now reversing it." << std::endl  ;
					polygon.reverse_orientation();
				} else NEWLINE ;
				polygonSet.push_back ( polygon );
			} // next polygon in this PWH

			// Now we process the polygon set into a proper polygon with holes
			if ( polygonSet.size() > 0 )
			{
				Bezier_polygon_with_holes pwh ( polygonSet.front() );
				if ( polygonSet.size() > 1 )
				{
					Bezier_polygon_vector::const_iterator it;
					for ( it = CGAL::cpp0x::next ( polygonSet.begin() ); it != polygonSet.end(); ++ it ) pwh.add_hole ( *it );
				}

				if ( is_valid_polygon_with_holes ( pwh, S.traits() ) )
				{
					S.join ( pwh ) ;
					std::cout << "Successfully read bezier polygon with holes, totally comprising "
								<< polygonSet.size() << " polygons."
								<< std::endl ;
				}
				else std::cout << "**** Bezier polygon with holes is NOT VALID ****" << std::endl ;
			}
			success = true ; // TODO: needs some refining, since even if this PWH is not valid, this is set to true
		} // next PWH in this file
	}
	catch ( const std::exception & x )
	{
		std::cout << "An exception occurred during reading of this bezier polygon set:" << x.what() << std::endl ;
	}
	catch ( ... )
	{
		std::cout << "An exception occurred during reading of this bezier polygon set." << std::endl ;
	}
	return success ;
}

template <typename Iter>
void print_double_points ( Iter begin, Iter end )
{
	for ( Iter i = begin ; i != end ; ++ i )
		std::cout << "(" << i->x().doubleValue() << "," << i->y().doubleValue() << ")," ;
	NEWLINE ;
}

void print_time_bounds ( const Bezier_alg_point * pt, const Bezier_curve * curve, unsigned int xid )
{
	// Adapted from include/CGAL/Arr_geometry_traits/Bezier_x_monotone_2.h parameter_range() function
	
	Bezier_alg_point::Originator_iterator  i = pt->get_originator (*curve, xid);
	CGAL_assertion (i != pt->originators_end());
	
	double min_t = CGAL::to_double(i->point_bound().t_min),
	       max_t = CGAL::to_double(i->point_bound().t_max) ;

	if ( max_t - min_t < APPROX_EPSILON ) std::cout << "a half-epsilon span near " << min_t << " " ;
	else std::cout << min_t << " and " << max_t << " spanning " << max_t - min_t << std::endl ;
}

bool refine ( const Bezier_alg_point * pt ) // returns whether refining happened
{
	unsigned int refineIterations = 0 ;
	g_refineTimer . start () ;
	while ( pt->refine() ) ++ refineIterations ;
	g_refineTimer . stop () ;
	g_totalRefinings += refineIterations ;
	if ( refineIterations == 0 ) std::cout << "which is already refined." << std::endl ;
	else std::cout << "Refined these in " << refineIterations << " iterations to: " ;
	return refineIterations != 0 ;
}

# include "de_casteljau.h"
void print_bezier_polygon ( const Bezier_polygon & b )
{
	std::vector<Rat_point> subcurveCtrlPts ;
	Bezier_polygon::Curve_const_iterator i ;
	unsigned int k ;
	std::pair<double, double> params ;
	double t0, t1 ;
	for ( i = b.curves_begin(), k = 1 ;
	      i != b.curves_end() ;
	      ++i, ++k ) // i points to each monotone bezier forming the polygon
	{
		const Bezier_curve * supportCurvePtr = & ( i->supporting_curve() ) ;
		const Bezier_alg_point * sourcePtr = & ( i->source() ) ;
		const Bezier_alg_point * targetPtr = & ( i->target() ) ;
		unsigned int xid = i->xid() ;

		# define GETPARAMS  params = i->parameter_range() ; t0 = params.first ; t1 = params.second
		GETPARAMS ;

		std::cout << "Segment " << k << std::endl ;

		std::cout << "Supporting curve control points " ;
		print_double_points ( supportCurvePtr->control_points_begin(), supportCurvePtr->control_points_end() );

		std::cout << "Segment from " << *sourcePtr << " to " << *targetPtr << std::endl ;

		# define PRINTSOURCE  print_time_bounds ( sourcePtr, supportCurvePtr, xid )
		# define PRINTTARGET  print_time_bounds ( targetPtr, supportCurvePtr, xid )
		
		if ( sourcePtr->is_exact() )
			std::cout << "Exact start time " << t0 << std::endl ;
		else
		{
			std::cout << "Start time bounds " ; PRINTSOURCE ; if ( refine ( sourcePtr ) ) { PRINTSOURCE ; NEWLINE ; }
		}

		if ( targetPtr->is_exact() )
			std::cout << "Exact end time "   << t1 << std::endl ;
		else
		{
			std::cout << "End time bounds "   ; PRINTTARGET ; if ( refine ( targetPtr ) ) { PRINTTARGET ; NEWLINE ; }
		}
		
		GETPARAMS ; // refreshing
		
		subcurveCtrlPts . clear () ;
		subcurveBezier ( supportCurvePtr->control_points_begin(), supportCurvePtr->control_points_end (),
		                                     t0, t1, std :: back_inserter ( subcurveCtrlPts ) ) ;
		std::cout << "Subcurve " << ( sourcePtr->is_exact() and targetPtr->is_exact() ? "exact" : "refined" ) << " control points " ;
		print_double_points ( subcurveCtrlPts . begin (), subcurveCtrlPts . end () ) ;
		NEWLINE ;
	}
}

// The main program.
int main ( int argc, char * argv[] )
{
	const char * filename1 = ( argc > 1 ) ? argv[1] : "char_g.bps";
	const char * filename2 = ( argc > 2 ) ? argv[2] : "char_m.bps";
	const char * bop       = ( argc > 3 ) ? argv[3] : "i";

	// Read the general polygons from the input files.
	CGAL::Timer           timer;
	Bezier_polygon_set S, S2 ;

	timer.start();

	if ( not read_bezier_polygon ( filename1, S ) )
	{
		std::cerr << "Failed to read " << filename1 << " ..." << std::endl;
		return 1;
	}

	if ( not read_bezier_polygon ( filename2, S2 ) )
	{
		std::cerr << "Failed to read " << filename2 << " ..." << std::endl;
		return 1;
	}

	timer.stop();
	std::cout << "Constructed the input polygons in " << timer.time() << " seconds." << std::endl ;
	timer.reset();

	std::cout << "Starting boolean operation..." << std::endl ;
	timer.start();
	try
	{
		switch ( bop[0] )
		{
		case 'i':
			S.intersection ( S2 );
			break ;
		case 'u':
			S.join        ( S2 );
			break ;
		}
	}
	catch ( std::exception const & x )
	{
		std::cout << "An exception occurred during the boolean operation:" << x.what() << std::endl ;
	}
	catch ( ... )
	{
		std::cout << "An exception occurred during the boolean operation." << std::endl;
	}
	timer.stop();

	std::cout << "The " << ( bop[0] == 'i' ? "intersection" : "union" ) << " computation took "
	          << timer.time() << " seconds." << std::endl;

	std::list<Bezier_polygon_with_holes> pwh ;
	std::list<Bezier_polygon_with_holes>::const_iterator pwhit;
	Bezier_polygon_with_holes::Hole_const_iterator hit;

	std::cout << "The result contains " << S.number_of_polygons_with_holes()
	          << " components:" << std::endl << std::endl ;

	std::cout << std::setprecision(15);
	unsigned int polynum, holenum;
	S.polygons_with_holes ( std::back_inserter ( pwh ) );
	for ( pwhit = pwh.begin(), polynum = 1; pwhit != pwh.end(); ++pwhit, ++polynum )
	{
		std::cout << "Polygon #" << polynum << " {" << std::endl;
		if ( not pwhit -> is_unbounded() )
		{
			std::cout << "Outer boundary = " << std::endl;
			print_bezier_polygon ( pwhit->outer_boundary() );
		}
		else std::cout << "No outer boundary." << std::endl;
		
		std::cout << "  " << pwhit->number_of_holes() << " holes" << std::endl;
		for ( hit = pwhit->holes_begin(), holenum = 1; hit != pwhit->holes_end(); ++hit, ++holenum )
		{
			std::cout << "    Hole #" << holenum << " = ";
			print_bezier_polygon ( *hit );
		}
		std::cout << "}" << std::endl << std::endl;
	}
	std::cout << "Refining the bezier vertices in " << g_totalRefinings
	          << " refine() calls took " << g_refineTimer.time() << " seconds" << std::endl ;
	return 0;
}

#endif