diff --git a/.gitattributes b/.gitattributes index 375421a17f8..63c04e580f5 100644 --- a/.gitattributes +++ b/.gitattributes @@ -1734,13 +1734,20 @@ Curved_kernel_via_analysis_2/include/CGAL/Arr_surfaces_intersecting_dupin_cyclid Curved_kernel_via_analysis_2/include/CGAL/Curved_kernel_via_analysis_2.h -text Curved_kernel_via_analysis_2/include/CGAL/Curved_kernel_via_analysis_2/Arc_2.h -text Curved_kernel_via_analysis_2/include/CGAL/Curved_kernel_via_analysis_2/Curve_interval_arcno_cache.h -text +Curved_kernel_via_analysis_2/include/CGAL/Curved_kernel_via_analysis_2/Curve_renderer_facade.h -text Curved_kernel_via_analysis_2/include/CGAL/Curved_kernel_via_analysis_2/Curved_kernel_via_analysis_2_functors.h -text Curved_kernel_via_analysis_2/include/CGAL/Curved_kernel_via_analysis_2/Generic_arc_2.h -text Curved_kernel_via_analysis_2/include/CGAL/Curved_kernel_via_analysis_2/Generic_point_2.h -text Curved_kernel_via_analysis_2/include/CGAL/Curved_kernel_via_analysis_2/Make_x_monotone_2.h -text Curved_kernel_via_analysis_2/include/CGAL/Curved_kernel_via_analysis_2/Non_x_monotone_arc_2.h -text Curved_kernel_via_analysis_2/include/CGAL/Curved_kernel_via_analysis_2/Point_2.h -text +Curved_kernel_via_analysis_2/include/CGAL/Curved_kernel_via_analysis_2/Qt_widget_Curve_renderer_2.h -text Curved_kernel_via_analysis_2/include/CGAL/Curved_kernel_via_analysis_2/Sweep_curves_adapter_2.h -text +Curved_kernel_via_analysis_2/include/CGAL/Curved_kernel_via_analysis_2/gfx/Curve_renderer_2.h -text +Curved_kernel_via_analysis_2/include/CGAL/Curved_kernel_via_analysis_2/gfx/Curve_renderer_internals.h -text +Curved_kernel_via_analysis_2/include/CGAL/Curved_kernel_via_analysis_2/gfx/Curve_renderer_traits.h -text +Curved_kernel_via_analysis_2/include/CGAL/Curved_kernel_via_analysis_2/gfx/Subdivision_1.h -text +Curved_kernel_via_analysis_2/include/CGAL/Curved_kernel_via_analysis_2/gfx/Subdivision_2.h -text Curved_kernel_via_analysis_2/include/CGAL/Curved_kernel_via_analysis_2/test/simple_models.h -text Curved_kernel_via_analysis_2/package_info/Curved_kernel_via_analysis_2/maintainer -text Curved_kernel_via_analysis_2/test/Curved_kernel_via_analysis_2/makefile -text diff --git a/Curved_kernel_via_analysis_2/include/CGAL/Curved_kernel_via_analysis_2/Curve_renderer_facade.h b/Curved_kernel_via_analysis_2/include/CGAL/Curved_kernel_via_analysis_2/Curve_renderer_facade.h new file mode 100644 index 00000000000..dd46d2bb615 --- /dev/null +++ b/Curved_kernel_via_analysis_2/include/CGAL/Curved_kernel_via_analysis_2/Curve_renderer_facade.h @@ -0,0 +1,286 @@ +// TODO: Add licence +// +// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE +// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. +// +// $URL:$ +// $Id: $ +// +// +// Author(s) : Pavel Emeliyanenko +// +// ============================================================================ + +/*!\file CGAL/Curved_kernel_via_analysis_2/Curve_renderer_facade.h + * \brief + * definition of \c Curve_renderer_facade<> + * + * high-level interface to the curve renderer + */ + +#ifndef CGAL_CKVA_CURVE_RENDERER_FACADE_H +#define CGAL_CKVA_CURVE_RENDERER_FACADE_H + +// do not compile curve renderer code (for fast debugging) +//#define CGAL_CKVA_DUMMY_RENDERER + +// whether to use multi-precision arithmetic +#define CGAL_CKVA_USE_MULTIPREC_ARITHMETIC +// whether to use exact rational arithmetic +#define CGAL_CKVA_USE_RATIONAL_ARITHMETIC + +#if AcX_SQRT_EXTENSION // not use rational arithmetic when sqrt-extended +#undef CGAL_CKVA_USE_RATIONAL_ARITHMETIC +#endif + +// this turns on a signleton curve renderer +// (not recommended for multi-threaded applications) +//#define CGAL_CKVA_USE_STATIC_RENDERER + +// prints out detailed info of the rendering process +// (only for internal debugging) +#ifndef Gfx_DETAILED_OUT +//#define Gfx_USE_DETAILED_OUT +#ifdef Gfx_USE_DETAILED_OUT +#define Gfx_DETAILED_OUT(x) std::cerr << x +#else +#define Gfx_DETAILED_OUT(x) static_cast(0) +#endif +#endif + +// prints out only high-level debug messages +#ifndef Gfx_OUT +//#define Gfx_USE_OUT +#ifdef Gfx_USE_OUT +#define Gfx_OUT(x) std::cerr << x +#else +#define Gfx_OUT(x) static_cast(0) +#endif +#endif + +#include +#include +#include + +#include + +CGAL_BEGIN_NAMESPACE + +/*! \brief + * CKvA interface to the curve renderer. Provides three levels of increasing + * arithmetic precision + * + * the curve renderer is instantiated with the base \c Float_ which can be + * either integral or user-defined floating-point number type, and optionally + * with multi-precision \c Bigfloat and exact rational number types as + * provided by the currently used \c Arithmetic_kernel + */ +template +class Curve_renderer_facade +{ +public: + //! \name Public typedefs + //!@{ + + //! this instance's first argument + typedef CurvedKernelViaAnalysis_2 Curved_kernel_via_analysis_2; + + //! base floating-point number type + typedef Float_ Float; + + //! type of a curve kernel + typedef typename Curved_kernel_via_analysis_2::Curve_kernel_2 + Curve_kernel_2; + + //! exact rational number type + typedef typename ::CGAL::Get_arithmetic_kernel< + typename Curve_kernel_2::Coefficient>::Arithmetic_kernel + Arithmetic_kernel; + + //! multi-precision float NT + //typedef typename Arithmetic_kernel::Bigfloat Bigfloat; + typedef typename Arithmetic_kernel::Bigfloat_interval BFI; + + typedef typename CGAL::Bigfloat_interval_traits::Boundary + Bigfloat; + + //! rational NT + typedef typename Arithmetic_kernel::Rational Rational; + + //! an arc of generic curve + typedef typename Curved_kernel_via_analysis_2::Arc_2 Arc_2; + + //! a point on curve + typedef typename Curved_kernel_via_analysis_2::Point_2 Point_2; + +#ifndef CGAL_CKVA_DUMMY_RENDERER + //! instance of a curve renderer + typedef Curve_renderer_2 + Default_renderer_2; +#ifdef CGAL_CKVA_USE_MULTIPREC_ARITHMETIC + //! the curve renderer instantiated with BigFloat + typedef Curve_renderer_2 + Bigfloat_renderer_2; +#endif +#ifdef CGAL_CKVA_USE_RATIONAL_ARITHMETIC + //! the curve renderer instantiated with Rationals + typedef Curve_renderer_2 + Exact_renderer_2; +#endif + +#ifdef CGAL_CKVA_USE_STATIC_RENDERER + static Default_renderer_2& renderer() + { + static Default_renderer_2 rend; + return rend; + } +#ifdef CGAL_CKVA_USE_MULTIPREC_ARITHMETIC + static Bigfloat_renderer_2& bigfloat_renderer() + { + static Bigfloat_renderer_2 rend; + return rend; + } +#endif +#ifdef CGAL_CKVA_USE_RATIONAL_ARITHMETIC + static Exact_renderer_2& exact_renderer() + { + static Exact_renderer_2 rend; + return rend; + } +#endif +#else // !CGAL_CKVA_USE_STATIC_RENDERER + Default_renderer_2 rend; + Default_renderer_2& renderer() + { + return rend; + } +#ifdef CGAL_CKVA_USE_MULTIPREC_ARITHMETIC + Bigfloat_renderer_2 bigfloat_rend; + Bigfloat_renderer_2& bigfloat_renderer() + { + return bigfloat_rend; + } +#endif +#ifdef CGAL_CKVA_USE_RATIONAL_ARITHMETIC + Exact_renderer_2 exact_rend; + Exact_renderer_2& exact_renderer() + { + return exact_rend; + } +#endif +#endif // !CGAL_CKVA_USE_STATIC_RENDERER + +#endif // !CGAL_CKVA_DUMMY_RENDERER + + //!@} +public: + //! \name Public methods and properties + //!@{ + + /*!\brief + * initializes renderer with drawing window dimensions + * and pixel resolution + * + * [x_min; y_min]x[x_max; y_max] - drawing window + * \c res_w and \c res_h - h/v pixel resolution + */ + inline void setup(const Bbox_2& bbox, int res_w, int res_h) + { +#ifndef CGAL_CKVA_DUMMY_RENDERER + renderer().setup(bbox, res_w, res_h); +#endif // !CGAL_CKVA_DUMMY_RENDERER + } + + /*!\brief + * rasterizes an x-monotone curve \c arc + * + * returns a list of sequences of pixel coordinates in \c points and + * end-point coordinats in \c end_points + */ + inline void draw(const Arc_2& arc, std::list& points, + std::pair& end_points) + { +#ifndef CGAL_CKVA_DUMMY_RENDERER + try { + renderer().draw(arc, points, end_points); + } + catch(CGALi::Insufficient_rasterize_precision_exception) { + std::cerr << "Switching to multi-precision arithmetic" << + std::endl; +#ifdef CGAL_CKVA_USE_MULTIPREC_ARITHMETIC + Bbox_2 bbox; + int res_w, res_h; + if(::boost::is_same::Is_exact, + CGAL::Tag_true>::value) + goto Lexit; + renderer().get_window(bbox); + renderer().get_resolution(res_w, res_h); + bigfloat_renderer().setup(bbox, res_w, res_h); + try { + points.clear(); + bigfloat_renderer().draw(arc, points, end_points); + return; + } + catch(CGALi::Insufficient_rasterize_precision_exception) { + + std::cerr << "Switching to exact arithmetic" << std::endl; +#ifdef CGAL_CKVA_USE_RATIONAL_ARITHMETIC + Bbox_2 bbox; + int res_w, res_h; + if(::boost::is_same::Is_exact, + CGAL::Tag_true>::value) + goto Lexit; + + renderer().get_window(bbox); + renderer().get_resolution(res_w, res_h); + exact_renderer().setup(bbox, res_w, res_h); + + try { + points.clear(); + exact_renderer().draw(arc, points, end_points); + return; + } + catch(CGALi::Insufficient_rasterize_precision_exception) { + goto Lexit; + } + return; +#endif // CGAL_CKVA_USE_RATIONAL_ARITHMETIC + } +Lexit: std::cerr << "Sorry, this does not work even with exact " + "arithmetic, bailing out..." << std::endl; + std::cerr << "polynomial: " << renderer().curve().polynomial_2() << + std::endl; + + renderer().get_window(bbox); + std::cerr << "window: " << bbox << std::endl; + +#endif // CGAL_CKVA_USE_MULTIPREC_ARITHMETIC + } +#endif // !CGAL_CKVA_DUMMY_RENDERER + } + + /*!\brief + * rasterizes a point on curve + */ + bool draw(const Point_2& point, CGALi::Coord_2& coord) { +#ifndef CGAL_CKVA_DUMMY_RENDERER + try { + return renderer().draw(point, coord); + } + catch(CGALi::Insufficient_rasterize_precision_exception) { + std::cerr << "Unable to rasterize point..\n"; + return false; + } +#else + return true; +#endif // !CGAL_CKVA_DUMMY_RENDERER + } + + //!@} +}; // Curve_renderer_facade + + +CGAL_END_NAMESPACE + +#endif // CGAL_CKVA_CURVE_RENDERER_FACADE_H + diff --git a/Curved_kernel_via_analysis_2/include/CGAL/Curved_kernel_via_analysis_2/Qt_widget_Curve_renderer_2.h b/Curved_kernel_via_analysis_2/include/CGAL/Curved_kernel_via_analysis_2/Qt_widget_Curve_renderer_2.h new file mode 100644 index 00000000000..a3484f81447 --- /dev/null +++ b/Curved_kernel_via_analysis_2/include/CGAL/Curved_kernel_via_analysis_2/Qt_widget_Curve_renderer_2.h @@ -0,0 +1,171 @@ +// TODO: Add licence +// +// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE +// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. +// +// $URL:$ +// $Id: $ +// +// +// Author(s) : Pavel Emeliyanenko +// +// ============================================================================ + +/*!\file CGAL/IO/Qt_widget_Curve_renderer_2.h + * \brief + * provides \c CGAL::Qt_widget interface for the curve renderer + */ + +#ifndef CGAL_QT_WIDGET_CURVE_RENDERER_2_H +#define CGAL_QT_WIDGET_CURVE_RENDERER_2_H + +#include + +#include +#include + +#include +#include +#include +#include + +#include + +CGAL_BEGIN_NAMESPACE + +namespace CGALi { + +/*!\brief + * singleton providing drawing routines for \c CurvedKernelViaAnalysis_2 + * + * @warning not recommended to use for multi-threaded applications + */ +template +class Renderer_instance +{ + Renderer_instance() { // private constructor + } + +public: + typedef typename CurvedKernelViaAnalysis_2::Curve_kernel_2::Coefficient + Coefficient; + + typedef CGAL::Interval_nt Interval_double; + + typedef Curve_renderer_facade + Curve_renderer_inst; + + static void resize_window(const Qt_widget& ws) + { + CGAL::Bbox_2 new_box(ws.x_min(), ws.y_min(), ws.x_max(), ws.y_max()); + if(bbox != new_box) { + bbox = new_box; + gfx_instance().setup(bbox, ws.width(), ws.height()); + } + } + + static void draw(Qt_widget& ws, + const typename Curve_renderer_inst::Arc_2& arc) + { + resize_window(ws); + QPainter *ppnt = &ws.get_painter(); + int height = ws.height(); + std::list points; + std::pair end_points; + + gfx_instance().draw(arc, points, end_points); + if(points.empty()) + return; + + std::list::const_iterator lit = points.begin(); + while(lit != points.end()) { + + const Coord_vec_2& vec = *lit; + Coord_vec_2::const_iterator vit = vec.begin(); + + if(vec.size() == 2) { + ppnt->moveTo(vit->x, height - vit->y); + vit++; + ppnt->lineTo(vit->x, height - vit->y); + + } else { + ppnt->moveTo(vit->x, height - vit->y); + //std::cout << "(" << (*it).x << "; " << height - (*it).y << + //")\n"; + while(vit != vec.end()) { + ppnt->lineTo(vit->x, height - vit->y); + //std::cout << "(" << (*it).x << "; " << height - (*it).y << + //")\n"; + vit++; + } + } + lit++; + } + + QPen old_pen = ppnt->pen(); + ppnt->setPen(QPen(Qt::NoPen)); // avoid drawing outlines + // draw with the current brush attributes + ppnt->drawEllipse(end_points.first.x-3,height-end_points.first.y-3, + 6, 6); + ppnt->drawEllipse(end_points.second.x-3,height-end_points.second.y-3, + 6, 6); + ppnt->setPen(old_pen); + } + + static void draw(Qt_widget& ws, + const typename Curve_renderer_inst::Point_2& pt) + { + resize_window(ws); + QPainter *ppnt = &ws.get_painter(); + int height = ws.height(); + + Coord_2 coord; + if(!gfx_instance().draw(pt, coord)) + return; + + QPen old_pen = ppnt->pen(); + ppnt->setPen(QPen(Qt::NoPen)); + ppnt->drawEllipse(coord.x - 3, height - coord.y - 3, 6, 6); + ppnt->setPen(old_pen); + } + +private: + static CGAL::Bbox_2 bbox; + + inline static Curve_renderer_inst& gfx_instance() { + static Curve_renderer_inst _instance; + return _instance; + } +}; + +template +CGAL::Bbox_2 Renderer_instance::bbox(0.0, 0.0, 0.0, 0.0); + +} // namespace CGALi + +/*! \brief + * outputs a curve arc to \c Qt_widget + */ +template +Qt_widget& operator << (Qt_widget& ws, + const CGALi::Arc_2& arc) { + + CGALi::Renderer_instance::draw(ws, arc); + return ws; +} + +/*! \brief + * outputs a curve point to \c Qt_widget + */ +template +Qt_widget& operator << (Qt_widget& ws, + const CGALi::Point_2& pt) { + + CGALi::Renderer_instance::draw(ws, pt); + return ws; +} + +CGAL_END_NAMESPACE + +#endif // CGAL_QT_WIDGET_CURVE_RENDERER_2_H + diff --git a/Curved_kernel_via_analysis_2/include/CGAL/Curved_kernel_via_analysis_2/gfx/Curve_renderer_2.h b/Curved_kernel_via_analysis_2/include/CGAL/Curved_kernel_via_analysis_2/gfx/Curve_renderer_2.h new file mode 100644 index 00000000000..3002981e709 --- /dev/null +++ b/Curved_kernel_via_analysis_2/include/CGAL/Curved_kernel_via_analysis_2/gfx/Curve_renderer_2.h @@ -0,0 +1,2653 @@ +// TODO: Add licence +// +// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE +// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. +// +// $URL:$ +// $Id: $ +// +// +// Author(s) : Pavel Emeliyanenko +// +// ============================================================================ + +/*!\file CGAL/Curved_kernel_via_analysis_2/gfx/Curve_renderer_2.h + * \brief definition of Curve_renderer_2<> + * rasterization of algebraic curves + */ + +#ifndef CGAL_CKVA_CURVE_RENDERER_2_H +#define CGAL_CKVA_CURVE_RENDERER_2_H + +#include +#include +#include +#include +#include +#include + +#include +#include + +#include +#include + +#include +#include + +using boost::multi_index::multi_index_container; +using boost::multi_index::get; +using boost::multi_index::project; + +CGAL_BEGIN_NAMESPACE + +#ifndef CGAL_CURVE_RENDERER_DEFS +#define CGAL_CURVE_RENDERER_DEFS + +// subdivision level beyond which visibly coincide branches are not further +// discriminated +#define CGAL_COINCIDE_LEVEL 4 + +// maximal recursion depth for derivative check +#define CGAL_DER_CHECK_DEPTH 5 + +// # of pixels to enlarge the drawing window in y-direction: required to skip +// closely located clip-points +#define CGAL_WINDOW_ENLARGE 15 + +// refine factor for intervals in x-direction (in pixel size) +#define CGAL_REFINE_X 100 + +// refine factor for intervals in y-direction (in pixel size) +#define CGAL_REFINE_Y 100000 + +// refine factor for clip-points +#define CGAL_REFINE_CLIP_POINTS 1000 + +#endif // CGAL_CURVE_RENDERER_DEFS + +/*! + * \brief The class template \c Curve_renderer_2 and its associate functions. + * + * The class implements rendering of distinct curve arcs and points as + * defined by \c CurvedKernelViaAnalysis_2::Arc_2 and + * \c CurvedKernelViaAnalysis_2::Point_2. \c Coeff_ template parameter + * defines an underlying number type to be used in polynomial and range + * evaluations. Valid instantiations are \c double, multi-precision float or + * exact rational number type. The main float-point number type is defined by + * \c Curve_renderer_traits::Float. + */ +//!@{ +template +class Curve_renderer_2 +{ +public: + //! \name public typedefs + //!@{ + + //! this instance's first template argument + typedef CurvedKernelViaAnalysis_2 Curved_kernel_via_analysis_2; + + //! this instance's second template argument + typedef Coeff_ Coeff; + + //! type of embedded curve kernel + typedef typename Curved_kernel_via_analysis_2::Curve_kernel_2 + Curve_kernel_2; + + //! type of x-monotone arc + typedef typename Curved_kernel_via_analysis_2::X_monotone_curve_2 Arc_2; + + //! type of a point on curve + typedef typename Curved_kernel_via_analysis_2::Point_2 Point_2; + + //! type of 1-curve analysis + typedef typename Curved_kernel_via_analysis_2::Curve_2 Curve_analysis_2; + + //!@} +private: + //! \name private typedefs + //!@{ + + //! type of X_coordinate + typedef typename Curve_kernel_2::X_coordinate_1 X_coordinate_1; + + //! type of xy-coordinate ? ;) + typedef typename Curve_kernel_2::Xy_coordinate_2 Xy_coordinate_2; + + //! rendering low-level routine + typedef CGALi::Curve_renderer_internals + Renderer_internals; + + //! exact rational number type + typedef typename Renderer_internals::Rational Rational; + + /// polynomial traits should be used whenever possible + typedef typename Renderer_internals::Polynomial_traits_2 + Polynomial_traits_2; + + //! curve renderer type conversion traits + typedef typename Renderer_internals::Renderer_traits Renderer_traits; + + //! coercion between rational and polynomial coefficients number type + typedef typename Renderer_internals::Rat_coercion_type Rat_coercion_type; + + //! instance of X_real_traits + typedef typename Curve_kernel_2::X_real_traits_1 X_real_traits_1; + + //! instance of Y_real_traits + typedef typename Curve_kernel_2::Y_real_traits_1 Y_real_traits_1; + + //! lower/upper boundary and one step refinement for \c X_coordinate_1 + typename X_real_traits_1::Lower_boundary lbound_x; + typename X_real_traits_1::Upper_boundary ubound_x; + typename X_real_traits_1::Refine refine_x; + + //! lower/upper boundary and one step refinement for \c Xy_coordinate_2 + typename Y_real_traits_1::Lower_boundary lbound_y; + typename Y_real_traits_1::Upper_boundary ubound_y; + typename Y_real_traits_1::Refine refine_y; + + //! specialized integer number type + typedef typename Renderer_traits::Integer Integer; + + //! event line instance type + typedef typename Curve_analysis_2::Status_line_1 Status_line_1; + + //! underlying bivariate polynomial type + typedef typename Renderer_internals::Polynomial_2 Polynomial_2; + + //! underlying univariate polynomial type + typedef typename Renderer_internals::Poly_dst_1 Poly_dst_1; + + //! basic number type used in all computations + typedef typename Renderer_internals::NT NT; + + //! instance of a univariate polynomial + typedef typename Renderer_internals::Poly_1 Poly_1; + //! instance of a bivariate polynomial + typedef typename Renderer_internals::Poly_2 Poly_2; + + //! conversion from the basic number type to doubles + typename CGAL::Real_embeddable_traits::To_double to_double; + + //! conversion from the basic number type to integers + typename Renderer_traits::Rat_to_integer rat2integer; + //! conversion from \c Integer type to built-in integer + typename Renderer_traits::Float_to_int float2int; + + //! conversion from \c Rational type to used number type + typename Renderer_traits::Rat_to_float rat2float; + //! makes the result exact after inexact operation (applicable only for + //! exact number types + typename Renderer_traits::Make_exact make_exact; + + //! returns \c true when the precision limit for a specified number type is + //! reached + typename Renderer_traits::Precision_limit limit; + //! maximum level of subdivision dependending on speficied number type + static const unsigned MAX_SUBDIVISION_LEVEL = + Renderer_traits::MAX_SUBDIVISION_LEVEL; + + //! pixel instance type + typedef CGALi::Pixel_2_templ Pixel_2; + //! seed point instance type + typedef CGALi::Seed_point_templ Seed_point; + //! support for multiple seed points + typedef std::stack Seed_stack; + + //! map container element's type for maintaining a list of cache instances + typedef std::pair LRU_entry; + //! a range of x-coordinates to define bottom/top clip points + typedef CGALi::Clip_point_templ + Clip_point_entry; + //! a container of bottom/top clip points + typedef std::vector Clip_points; + //! an integer index container + typedef std::vector index_vector; + + //! LRU list used for effective cache switching + typedef boost::multi_index::multi_index_container< + LRU_entry, boost::multi_index::indexed_by< + boost::multi_index::sequenced<>, + boost::multi_index::ordered_unique< + BOOST_MULTI_INDEX_MEMBER(LRU_entry,int,first) > > > + LRU_list; + + //! internal use only: a stripe defined by bounding polynomials + //! and coordinates on sides, index 0 - lower boundary, 1 - upper boundary + struct Stripe { + Poly_1 poly[2]; + NT key[2]; + unsigned level[2]; + } ; + + //!@} +public: + //! \name constructor + //!@{ + + //! default constructor: a curve segment is undefined + Curve_renderer_2() : cache_id(-1), initialized(false), one(1) { + } + + //!@} +public: + //! \name public methods + //!@{ + + //! sets up drawing window and pixel resolution + void setup(const Bbox_2& bbox_, int res_w_, int res_h_) { + + initialized = engine.setup(bbox_, res_w_, res_h_); + clip_pts_computed = false; // need to recompute clip points each time + // the window dimensions are changed + } + + //! \brief returns currently used supporting algebraic curve + Curve_analysis_2 curve() const + { + return *support; + } + + //! \brief returns the drawing window + void get_window(Bbox_2& bbox_) const { + bbox_ = Bbox_2(to_double(engine.x_min), to_double(engine.y_min), + to_double(engine.x_max), to_double(engine.y_max)); + } + + //! \brief returns the pixel resolution + void get_resolution(int& res_w_, int& res_h_) const { + res_w_ = engine.res_w; + res_h_ = engine.res_h; + } + + //! \brief the main rendering procedure for curve arcs + //! returns a sequence of integer pixel coordinates \c points + //! + //! An exception \c Insufficient_rasterize_precision_exception is + //! thrown whenever the precision of currently used NT is not enough + //! to correctly render a curve arc. The exception has to be caught + //! outside the renderer to switch to a higher-precision arithmetic + void draw(const Arc_2& arc, std::list& points, + std::pair& end_points); + + /*!\brief + * overloaded version to rasterize points on curves + */ + bool draw(const Point_2& pt, CGALi::Coord_2& coord); + + //!@} +private: + //! \name Private methods + //!@{ + + //! \brief switches to another cache instance depending on the + //! supporting curve of a segment + void select_cache_entry(const Arc_2& arc); + + //! \brief draws a piece of a semgment from pix_1 to pix_2 + void draw_lump(CGALi::Coord_vec_2& points, int& last_x, int arcno, + const Pixel_2& pix_1, const Pixel_2& pix_2); + + //! \brief copmutes an isolating interval for y-coordinate of an end-point, + //! uses caching if possible + Rational get_endpoint_y(const Arc_2& arc, const Rational& x, + CGAL::Arr_curve_end end, bool is_clipped); + + //! \brief returns whether a polynomial has zero over an interval, + //! we are not interested in concrete values + //! + //! a bit mask \c check indicates which boundaries are to be computed + //! 0th bit - of a polynomial itself (0th derivative, default) + //! 1st bit - of the first derivative (sets \c first_der flag) + //! 2nd bit - of the second derivative + bool get_range_1(int var, const NT& lower, const NT& upper, const NT& key, + const Poly_1& poly, int check = 1); + + //! \brief advances pixel's coordinates by given increments + void advance_pixel(Pixel_2& pix, int new_dir) + { + int x_inc = CGALi::directions[new_dir].x, + y_inc = CGALi::directions[new_dir].y; + if(pix.level == 0) { + pix.x += x_inc; + pix.y += y_inc; + } else { + Integer x = pix.sub_x + x_inc, + y = pix.sub_y + y_inc, pow = (one << pix.level) - 1; + (x < 0 ? pix.x-- : x > pow ? pix.x++ : x); + (y < 0 ? pix.y-- : y > pow ? pix.y++ : y); + pix.sub_x = x&pow; + pix.sub_y = y&pow; + } +// int taken = CGALi::DIR_TAKEN_MAP[new_dir]; +// if(taken != -1&&taken!=direction_taken) { +// Gfx_OUT("ERROR: direction reversed at pixel: " << +// pix << " new_dir = " << new_dir << std::endl); +// throw -1; +// } + } + + //! computes pixel coordinates from rational point + void get_pixel_coords(const Rational& x, const Rational& y, + Pixel_2& pix, Rational *ppix_x=NULL, Rational *ppix_y=NULL) + { + Rational p_x = (x - engine.x_min_r) / engine.pixel_w_r, + p_y = (y - engine.y_min_r) / engine.pixel_h_r; + pix.x = static_cast(std::floor(CGAL::to_double(p_x))); + pix.y = static_cast(std::floor(CGAL::to_double(p_y))); + + if(ppix_x != NULL && ppix_y != NULL) { +// NiX::simplify(p_x); +// NiX::simplify(p_y); + *ppix_x = p_x; + *ppix_y = p_y; + } + } + + //! refines y-coordinate of \c Xy_coordinate_2 to a certain bound + void refine_xy(const Xy_coordinate_2& xy, const Rational& bound) { + while(ubound_y(xy) - lbound_y(xy) >= bound) + refine_y(xy); + } + + //! returns h/v boundaries for 8-pixel neighbourhood + void get_boundaries(int var, const Pixel_2& pix, Stripe& stripe); + + //! returns univariate polynomials (depending on subdivision level cache + //! might be used) + void get_polynomials(int var, Stripe& stripe); + + //! checks 8-pixel neighbourhood of a pixel, returns \c true if + //! only one curve branch intersects pixel's neighbourhood, \c dir + //! defines backward direction, \c new_dir is a new tracking direction + bool test_neighbourhood(const Pixel_2& pix, int dir, int& new_dir, + int step_x = 2, int step_y = 2); + +//#ifdef Gfx_DEBUG_PRINT + // debug only! + void dump_neighbourhood(const Pixel_2& pix); +//#else +// #define dump_neighbourhood ((void)0) +//#endif + + //! returns a subpixel which is crossed by the curve branch, subpixels + //! are tested w.r.t. preferred direction (only for h/v directions) + int get_subpixel_hv(const Pixel_2& pix, int pref_dir); + + //! the same for diagonal directions + int get_subpixel_diag(const Pixel_2& pix, int pref_dir); + + //! recursively subdivides pixel into 4 subpixels, returns a new tracking + //! direction + bool subdivide(Pixel_2& pix, int dir, int& new_dir); + + //! returns the starting witness pixel and two directions from it + bool get_seed_point(const Rational& seed, int arcno, Pixel_2& start, + int *dir, int *b_taken, bool& b_coincide); + + //! tests whether a polynomial has only one root over an interval + bool recursive_check(int var, const NT& beg, const NT& end, + const NT& key, const Poly_1& poly, int depth = 0); + + //! if success returns two directions among 8, where a curve branch + //! crosses the pixel's neighbourhood + bool test_pixel(const Pixel_2& pix, int *dir, int *b_taken, + bool& b_coincide); + + //! computes an isolating box of an end-point + bool get_isolating_box(const Rational& x_s, const Rational& y_s, + Pixel_2& res); + + //! returns true if \c pix is encompassed into one of isolating rectangles + //! (stopping criteria) + bool is_isolated_pixel(const Pixel_2& pix); + + //! computes clip points on top and bottom window boundaries + void horizontal_clip(); + + //! refines an algebraic point w.r.t. certain criteria + void refine_alg_point(Rational& l, Rational& r, const Poly_dst_1& poly, + const Rational& criteria, int mode=0); + + //! computes bottom and top horizontal clip-points of a segment + void segment_clip_points(const Rational& x_lower, const Rational& x_upper, + const Rational& y_lower, const Rational& y_upper, + const Rational& y_clip, const Poly_dst_1& poly, int arcno, + Clip_points& clip_points, index_vector& clip_indices); + +public: + //! destructor + ~Curve_renderer_2() + { + cache_list.clear(); + } + + //!@} +private: + //! \name Private properties + //!@{ + + unsigned max_level; //! maximum subdivision level + unsigned current_level; + + bool clip_pts_computed; //! indicates whether clip points are computed + + Pixel_2 isolated_l, isolated_h; //! isolating rectangles for lower and + //! upper end-points + + //! an instance of rendering engine + CGALi::Curve_renderer_internals engine; + + Curve_analysis_2 *support; //! supporting 1-curve analysis + Curve_analysis_2 support_[CGAL_N_CACHES]; + + Clip_points btm_clip, top_clip; //! a set of bottom and top clip points + Poly_dst_1 btm_poly, top_poly; //! bottom and top polynomials + + int cache_id; //! index of currently used cache instance + LRU_list cache_list; //! list of indices of cache instances + + Seed_stack s_stack; //! a stack of seed points + Seed_point current_seed; //! current seed point + + bool initialized; //! indicates whether the renderer has been initialized + //! with correct parameters + const Integer one; //! just "one" + bool branches_coincide; //! indicates that there are several branches + //! passing through one neighbourhood pixel + int direction_taken; //! stores a direction taken from the seed point + //! during tracking, if it's possible to determine + //! 0 - towards lower point, 1 - towards upper + //!@} +}; // class Curve_renderer_2<> + +//! \brief main rasterization procedure, takes coordinates of event points and +//! the arc number and draws a curve segment +template +void Curve_renderer_2::draw( + const Arc_2& arc, std::list& points, + std::pair& end_points) +{ + if(!initialized) + return; + select_cache_entry(arc); // lookup for appropriate cache instance + + Gfx_OUT("\n////////////////////\n resolution: " << engine.res_w << " x " + << engine.res_h << std::endl); + Gfx_OUT("box: [" << engine.x_min << "; " << engine.y_min << "]x[" << + engine.x_max << "; " << engine.y_max << "]" << std::endl); + + Rational lower, upper, y_lower = 0, y_upper = 0; + bool infty_src = true, infty_tgt = true, clip_src=false, clip_tgt=false; + + CGAL::Arr_parameter_space loc_p1 = arc.location(CGAL::ARR_MIN_END), + loc_p2 = arc.location(CGAL::ARR_MAX_END); + + if(loc_p1 != CGAL::ARR_LEFT_BOUNDARY) { + const X_coordinate_1& x0 = arc.curve_end_x(CGAL::ARR_MIN_END); + + //while(x0.high() - x0.low() > engine.pixel_w_r/CGAL_REFINE_X) + // x0.refine(); + //lower = p1.finite().high(); + //if(!p1.finite().is_rational()&&!seg.is_vertical()) + // arcno_p1 = seg.arcno(); + + while(ubound_x(x0) - lbound_x(x0) > engine.pixel_w_r/CGAL_REFINE_X) + refine_x(x0); + + lower = ubound_x(x0); +// NiX::simplify(lower); + } else + lower = engine.x_min_r; + + // since end-points are sorted lexicographically, no need to check for + // lower boundary + if(loc_p2 == CGAL::ARR_TOP_BOUNDARY && arc.is_vertical()) + y_upper = engine.y_max_r; + + if(loc_p2 != CGAL::ARR_RIGHT_BOUNDARY) { + const X_coordinate_1& x0 = arc.curve_end_x(CGAL::ARR_MAX_END); + while(ubound_x(x0) - lbound_x(x0) > engine.pixel_w_r/CGAL_REFINE_X) + refine_x(x0); + + upper = lbound_x(x0); +// NiX::simplify(upper); + } else + upper = engine.x_max_r; + + if(upper <= engine.x_min_r||lower >= engine.x_max_r) + return; + + if(lower < engine.x_min_r) { + lower = engine.x_min_r; + clip_src = true; + } + if(upper > engine.x_max_r) { + upper = engine.x_max_r; + clip_tgt = true; + } + + Rational height_r = (engine.y_max_r-engine.y_min_r)*2; + switch(loc_p1) { + case CGAL::ARR_BOTTOM_BOUNDARY: + y_lower = (arc.is_vertical() ? engine.y_min_r : + engine.y_min_r - height_r); + break; + + case CGAL::ARR_TOP_BOUNDARY: + // endpoints must be sorted lexicographical + CGAL_precondition(!arc.is_vertical()); + y_lower = engine.y_max_r + height_r; + break; + + default: + infty_src = false; + y_lower = get_endpoint_y(arc, lower, CGAL::ARR_MIN_END, clip_src); + } + + switch(loc_p2) { + case CGAL::ARR_BOTTOM_BOUNDARY: + // endpoints must be sorted lexicographical + CGAL_precondition(!arc.is_vertical()); + y_upper = engine.y_min_r - height_r; + break; + + case CGAL::ARR_TOP_BOUNDARY: + y_upper = (arc.is_vertical() ? engine.y_max_r : + engine.y_max_r + height_r); + break; + default: + infty_tgt = false; + y_upper = get_endpoint_y(arc, upper, CGAL::ARR_MAX_END, clip_tgt); + } + +// NiX::simplify(y_lower); +// NiX::simplify(y_upper); + Pixel_2 pix_1, pix_2; + + Gfx_OUT("lower: " << CGAL::to_double(lower) << "; upper: " << + CGAL::to_double(upper) << "; y_lower: " << NiX::to_double(y_lower) << + "; y_upper: " << CGAL::to_double(y_upper) << "\n"); + + get_pixel_coords(lower, y_lower, pix_1); + get_pixel_coords(upper, y_upper, pix_2); + + end_points.first.x = (clip_src ? engine.res_w+4 : pix_1.x); + end_points.first.y = pix_1.y; + end_points.second.x = (clip_tgt ? engine.res_w+4 : pix_2.x); + end_points.second.y = pix_2.y; + + Gfx_OUT("lower pix: (" << pix_1.x << "; " << pix_1.y << + ") upper pix: (" << pix_2.x << "; " << pix_2.y << + ") pixel_w: " << engine.pixel_w << " pixel_h: " << engine.pixel_h << + std::endl); + + CGALi::Coord_vec_2 rev_points; + // reserve at least enough space for arc's x-length + rev_points.reserve(CGAL_ABS(pix_2.x - pix_1.x)); + + if(arc.is_vertical()) { + rev_points.push_back(CGALi::Coord_2(pix_1.x, pix_1.y)); + rev_points.push_back(CGALi::Coord_2(pix_2.x, pix_2.y)); + points.push_back(rev_points); + return; + } + + if(!clip_pts_computed) { + Gfx_OUT("computing clip points\n"); + horizontal_clip(); + clip_pts_computed = true; + } + + typedef std::pair Int_pair; + typedef std::vector index_pair_vector; + index_vector btm_idx, top_idx; + index_pair_vector seg_pts; + + Gfx_OUT("checking bottom clip: y = " << engine.y_min << "\n"); + segment_clip_points(lower, upper, y_lower, y_upper, engine.y_min_r, + btm_poly, arc.arcno(), btm_clip, btm_idx); + + Gfx_OUT("checking top clip: y = " << engine.y_max << "\n"); + segment_clip_points(lower, upper, y_lower, y_upper, engine.y_max_r, + top_poly, arc.arcno(), top_clip, top_idx); + + index_vector::iterator it1 = btm_idx.begin(), it2 = top_idx.begin(); + while(1) { + bool push_1 = true; + if(it1 != btm_idx.end()) { + if(it2 != top_idx.end()) + push_1 = ((btm_clip)[*it1].left < (top_clip)[*it2].left); + } else if(it2 != top_idx.end()) + push_1 = false; + else + break; + if(push_1) { // 0 - bottom point, 1 - top point + seg_pts.push_back(Int_pair(*it1,0)); + it1++; + } else { + seg_pts.push_back(Int_pair(*it2,1)); + it2++; + } + } + + max_level = 0; + Pixel_2 start; + int dir[2], b_taken[2], last_x = -engine.res_w; + bool b_coincide; + Rational pt, mid; + + branches_coincide = false; + // make a small tolerance + bool pt1_inside = (pix_1.y >= 0 && pix_1.y < engine.res_h), + pt2_inside = (pix_2.y >= 0 && pix_2.y < engine.res_h); + +//TODO: use Event1_info::multiplicity(i) - to gather information about roots ? + if(seg_pts.size() == 1 && pt1_inside && pt2_inside) + seg_pts.clear(); + + // easy case: no clip-points found + if(seg_pts.size() == 0) { + if(!pt1_inside&&!pt2_inside) { + Gfx_OUT("segment is outside\n"); + return; + } + /// WARNING: if x-interval is small while y coordinates are far away from + /// the window, we can get into the troubles.. + if(0){//pix_2.x - pix_1.x <= 1) {// it goes away right here + rev_points.push_back(CGALi::Coord_2(pix_1.x,pix_1.y)); + rev_points.push_back(CGALi::Coord_2(pix_2.x,pix_2.y)); + points.push_back(rev_points); + return; + } + Gfx_OUT("NO clip points\n"); + mid = (lower + upper)/2; +// NiX::simplify(mid); + pt = engine.x_min_r + (pix_1.x*2+5)*engine.pixel_w_r/2; +// NiX::simplify(pt); + // segment is almost vertical - start from a middle point + //if(pt > mid) + pt = mid; + if(!get_seed_point(pt, arc.arcno(), start, dir, b_taken, + b_coincide)) { + Gfx_OUT("generic error occurred..\n"); + return; + } + Gfx_OUT("starting pixel found: " << start << " directions: " << dir[0] + << " and " << dir[1] << std::endl); + // only one point list + draw_lump(rev_points, last_x, arc.arcno(), pix_1, pix_2); + points.push_back(rev_points); + + Gfx_OUT("exit normal, max_level: " << max_level << std::endl); + return; + } + + // clip-points presented + Rational l, r, y_clip; + X_coordinate_1 alpha; + Poly_dst_1 *ppoly; + Clip_point_entry *pclip, *ptmp; + + Gfx_OUT("segment is not completely inside the window..\n"); + index_pair_vector::iterator it = seg_pts.begin(), + eend = seg_pts.end(); + + Gfx_OUT("\n\nSEGMENT clip points: \n"); + Pixel_2 pix_beg, pix_end; + bool straight_segment; + + while(it != eend) { + if(it != seg_pts.begin() && it == eend-1 && !pt2_inside) + break; + straight_segment = false; + + int idx = it->first; + if(it->second == 0) {// bottom points + pclip = &btm_clip[idx]; + ppoly = &btm_poly; + y_clip = engine.y_min_r; + } else { // top points + pclip = &top_clip[idx]; + ppoly = &top_poly; + y_clip = engine.y_max_r; + } + l = pclip->left; + r = pclip->right; + + if(last_x != -engine.res_w) { // compute screen coordinates of a pixel + Rational last_pt = engine.x_min_r + static_cast(last_x)* + engine.pixel_w_r; + if(r <= last_pt) { // skip a clip-point if already drawn + Gfx_OUT("clip point skipped..\n"); + it++; + continue; + } + } + + rev_points.clear(); + if(it == seg_pts.begin() && pt1_inside) { + + Gfx_OUT("starting point is near the left end-point\n"); + //pt = engine.x_min_r + (pix_1.x*2+5)*engine.pixel_w_r/2; + pt = (lower + l)/2; +// NiX::simplify(pt); + if(pt <= lower) { + refine_alg_point(l, r, *ppoly, lower, 1); + pt = l; + if(l - lower <= engine.pixel_w_r*2) + straight_segment = true; + } + pix_beg = pix_1; + get_pixel_coords(l, y_clip, pix_end); + + } else if(it == eend - 1 && pt2_inside) { + + Gfx_OUT("starting point is near the right end-point\n"); + //pt = engine.x_min_r + (pix_2.x*2-3)*engine.pixel_w_r/2; + pt = (r + upper)/2; + Gfx_OUT("pt is: " << CGAL::to_double(pt) << "; r is: " << + CGAL::to_double(r) << ": pixel_w_r: " << + CGAL::to_double(engine.pixel_w_r) << "\n"); + +// NiX::simplify(pt); + if(pt >= upper) { + // ensure that clip-point interval is sufficiently small + refine_alg_point(l, r, *ppoly, upper, 1); + pt = r; + if(upper - r <= engine.pixel_w_r*2) + straight_segment = true; + } + pix_end = pix_2; + get_pixel_coords(r, y_clip, pix_beg); + + } else { + Gfx_OUT("starting point is between clip-points\n"); + it++; + + if(it == seg_pts.end()) { + Gfx_OUT("ERROR: clip point missed ?\n"); + break; + } + + idx = it->first; + ptmp = (it->second ? &top_clip[idx] : &btm_clip[idx]); + pt = pclip->alpha.rational_between(ptmp->alpha); + + get_pixel_coords(l, y_clip, pix_beg); + get_pixel_coords(ptmp->left, it->second ? engine.y_max_r : + engine.y_min_r,pix_end); + + if(CGAL_ABS(ptmp->left - l) <= engine.pixel_w_r*2) { + + Xy_coordinate_2 xy(X_coordinate_1(pt), *support, arc.arcno()); + refine_xy(xy, engine.pixel_h_r/CGAL_REFINE_Y); + mid = ubound_y(xy); + + rev_points.push_back(CGALi::Coord_2(pix_beg.x, pix_beg.y)); + get_pixel_coords(pt, mid, pix_beg); + straight_segment = true; + } + } + + if(straight_segment) { + Gfx_OUT("straight subsegment found\n"); + rev_points.push_back(CGALi::Coord_2(pix_beg.x, pix_beg.y)); + rev_points.push_back(CGALi::Coord_2(pix_end.x, pix_end.y)); + points.push_back(rev_points); + it++; + continue; + } + + if(!get_seed_point(pt, arc.arcno(), start, dir, b_taken, + b_coincide)) { + Gfx_OUT("get_seed_point: a problem occurred..\n"); + it++; + continue; + } + + Gfx_OUT("starting pixel found: " << start << " directions: " << + dir[0] << " and " << dir[1] << std::endl); + draw_lump(rev_points, last_x, arc.arcno(), pix_beg, pix_end); + points.push_back(rev_points); + if(it == eend) + break; + it++; + } + + /*if(!get_isolating_box(lower, y_lower, isolated_l)) + isolated_l.level = -1u; + if(!get_isolating_box(upper, y_upper, isolated_h)) + isolated_h.level = -1u;*/ +} + +//! draws a segment's piece from pix_1 to pix_2, the starting pixel is taken +//! from the seed point stack +template +void Curve_renderer_2::draw_lump( + CGALi::Coord_vec_2& rev_points, int& last_x, int arcno, + const Pixel_2& pix_1, const Pixel_2& pix_2) +{ + Pixel_2 start, pix, witness, prev_pix, stored_pix, stored_prev; + int back_dir, new_dir, stored_dir, dir[2], b_taken[2], ux, uy; + bool b_coincide, failed; + start = s_stack.top().start; + + bool pix_1_close = (CGAL_ABS(start.x - pix_1.x) <= 1&& + CGAL_ABS(start.y - pix_1.y) <= 1), + pix_2_close = (CGAL_ABS(start.x - pix_2.x) <= 1&& + CGAL_ABS(start.y - pix_2.y) <= 1); + + // stores an x-coordinate of the last traced pixel in --> direction: + // required to "jump" over clip-points + last_x = -engine.res_w; + + if(pix_1_close && pix_2_close) { // suppress drawing of trivial segments + rev_points.push_back(CGALi::Coord_2(pix_1.x,pix_1.y)); + rev_points.push_back(CGALi::Coord_2(pix_2.x,pix_2.y)); + return; + } + + CGALi::Coord_vec_2 points; + // reserve a list of point coordinates + points.reserve(CGAL_ABS(pix_2.x - pix_1.x)); + + direction_taken = -1; + bool overstep = (pix_1_close||pix_2_close); + bool ready = false; + while(!s_stack.empty()) + { + current_seed = s_stack.top(); + s_stack.pop(); + branches_coincide = current_seed.branches_coincide; + if(direction_taken == -1) + direction_taken = current_seed.direction_taken; + witness = current_seed.start; + pix = witness; + current_level = witness.level; + pix.sub_x = 0; + pix.sub_y = 0; + pix.level = 0; + back_dir = current_seed.back_dir; + stored_dir = back_dir; + stored_pix = pix; + stored_prev = pix; + + // store result in a list or reversed list depending on the orientation + CGALi::Coord_vec_2 *ppoints = (current_seed.orient == 0 ? + &rev_points : &points); + + Gfx_OUT("continuing from a seed point: " << witness << "; back_dir: " << + back_dir << " direction_taken: " << direction_taken << std::endl); + + bool is_exit = false; + while(1) { +// #ifdef CGAL_CKVA_SYNCHRONOUS_CANCEL +// pthread_testcancel(); // this is a cancellation point +// #endif + + if(pix.x < 0 || pix.x > engine.res_w || pix.y < -CGAL_WINDOW_ENLARGE|| + pix.y > engine.res_h + CGAL_WINDOW_ENLARGE) { + branches_coincide = false; + break; + } + ppoints->push_back(CGALi::Coord_2(pix.x, pix.y)); + + if((direction_taken == 0 && pix.x <= pix_1.x) || + (direction_taken == 1 && pix.x >= pix_2.x)) { + //Gfx_OUT("STOP: reached end-point x-coordinate\n"); + branches_coincide = false; + is_exit = true; + break; + } + + bool set_ready = false; + if(CGAL_ABS(pix.x - pix_1.x) <= 1 && CGAL_ABS(pix.y - pix_1.y) <= 1) { + if(!overstep||ready) { + pix.x = pix_1.x; + pix.y = pix_1.y; + branches_coincide = false; + break; + } + set_ready = true; + } + if(CGAL_ABS(pix.x - pix_2.x) <= 1 && CGAL_ABS(pix.y - pix_2.y) <= 1) { + if(!overstep||ready) { + pix.x = pix_2.x; + pix.y = pix_2.y; + branches_coincide = false; + break; + } + ready = true; + } + if(set_ready) + ready = true; + if(!test_neighbourhood(pix, back_dir, new_dir)) { + ux = pix.x; + uy = pix.y; + if(witness == pix) { // witness subpixel is a pixel itself + if(!subdivide(pix,back_dir,new_dir)) { + is_exit = true; + break; + } + prev_pix = pix; + advance_pixel(pix,new_dir); + back_dir = (new_dir+4)&7; + } else { + back_dir = stored_dir; + pix = witness; + prev_pix = witness; + } + stored_dir = -1; + failed = false; + while(ux == pix.x && uy == pix.y) { + if(!failed && ((prev_pix.sub_x & 2) != (pix.sub_x & 2)|| + (prev_pix.sub_y & 2) != (pix.sub_y & 2))) { + stored_pix = pix; + pix.sub_x >>= 1; + pix.sub_y >>= 1; + pix.level--; + stored_dir = back_dir; + stored_prev = prev_pix; + } + if(!test_neighbourhood(pix, back_dir, new_dir)) { + if(stored_dir != -1) { + pix = stored_pix; + prev_pix = stored_prev; + back_dir = stored_dir; + stored_dir = -1; + failed = true; + continue; + } + if(!subdivide(pix,back_dir,new_dir)) { + is_exit = true; + break; + } + } + //Gfx_OUT(pix << " dir = " << new_dir << std::endl); + prev_pix = pix; + advance_pixel(pix,new_dir); + if(is_isolated_pixel(pix)) { + branches_coincide = false; + is_exit = true; + break; + } + back_dir = (new_dir+4)&7; + } + if(is_exit) + break; + stored_dir = back_dir; + witness = pix; + pix.level = 0; + pix.sub_x = 0; + pix.sub_y = 0; + Gfx_OUT(witness << " " << prev_pix << std::endl); + } else { + ux = pix.x; + uy = pix.y; + advance_pixel(pix,new_dir); + witness = pix; + } + back_dir = (new_dir+4)&7; + } + + if(branches_coincide) { // oops, need another seed point + if(direction_taken == -1) { + Gfx_OUT("\n\nFATAL: unknown direction in coincide mode!\n\n"); + return; + } + + (direction_taken == 0) ? pix.x -= 1 : pix.x += 1; + if(pix.x >= pix_2.x-1||pix.x <= pix_1.x+1) + goto Lexit; + + Gfx_OUT("New seed point required at: " << pix << std::endl); + if(!get_seed_point(engine.x_min_r+pix.x*engine.pixel_w_r, arcno, + start, dir, + b_taken, b_coincide)) { + std::cerr << " wrong seed point found " << std::endl; + throw CGALi::Insufficient_rasterize_precision_exception(); + } + + Gfx_OUT("new seed point found: " << start << " directions: " << dir[0] + << " and " << dir[1] << "; taken = " << direction_taken << + std::endl); + + b_taken[0] = CGALi::DIR_TAKEN_MAP[dir[0]]; + b_taken[1] = CGALi::DIR_TAKEN_MAP[dir[1]]; + if(b_taken[0] != -1) + new_dir = dir[b_taken[0] != direction_taken ? 0: 1]; + else if(b_taken[1] != -1) + new_dir = dir[b_taken[1] != direction_taken ? 1: 0]; + else { + std::cerr << "ERROR: wrong backward dir after seed point: " << + dir[0] << " " << dir[1] << std::endl; + throw CGALi::Insufficient_rasterize_precision_exception(); + } + s_stack.push(Seed_point(start,new_dir, current_seed.orient, + direction_taken, b_coincide)); + continue; + } else if(is_exit) { +Lexit: if(direction_taken==0) { + pix.x = pix_1.x; + pix.y = pix_1.y; + } else if(direction_taken==1) { + pix.x = pix_2.x; + pix.y = pix_2.y; + } + } + ppoints->push_back(CGALi::Coord_2(pix.x, pix.y)); + + // we were tracing in --> direction: store the pixel where we stopped + if(direction_taken == 1) + last_x = pix.x; + if(direction_taken != -1) + direction_taken = 1 - direction_taken; + ready = false; + } + + std::reverse(rev_points.begin(), rev_points.end()); + // resize rev_points to accomodate the size of points vector + unsigned rsize = rev_points.size(); + rev_points.resize(rsize + points.size()); + std::copy(points.begin(), points.end(), rev_points.begin() + rsize); +} + +/*!\brief + * overloaded version to rasterize distinct points on curves + * + * \return \c false indicating that the point lies outside the drawing window + * \c true in case the point coordinates were successfully computed + */ +template +bool Curve_renderer_2::draw( + const Point_2& pt, CGALi::Coord_2& coord) { + + const X_coordinate_1& x0 = pt.x(); + while(ubound_x(x0) - lbound_x(x0) > engine.pixel_w_r/2) + refine_x(x0); + + Rational x_s = lbound_x(x0), y_s; +// NiX::simplify(x_s); + if(x_s < engine.x_min_r || x_s > engine.x_max_r) + return false; + + typename Curve_analysis_2::Internal_curve_2::Event1_info + event = pt.curve()._internal_curve().event_info_at_x(pt.x()); + event.refine_to(pt.arcno(), engine.pixel_h_r/CGAL_REFINE_X); + + y_s = event.lower_boundary(pt.arcno()); + + Pixel_2 pix; + get_pixel_coords(x_s, y_s, pix); + if(pix.x < 0 || pix.x >= engine.res_w || + pix.y < 0 || pix.y >= engine.res_h) + return false; + + coord.x = pix.x; + coord.y = pix.y; + return true; +} + +namespace CGALi +{ +// map from box sides to subpixel numbers (for h/v directions) +// these are old versions +static const int HV_SUBPIX_MAP[][4] = { + {1,3,0,2}, // right(0) + {3,2,1,0}, // top(1) + {2,0,3,1}, // left(2) + {0,1,2,3}}; // bottom(3) +static const int D_SUBPIX_MAP[][4] = { + {3,1,2,0}, // dir = 0 + {2,3,0,1}, // dir = 1 + {0,2,1,3}, // dir = 2 + {1,0,3,2}}; // dir = 3 +} // namespace CGALi + +//! recursively subdivides pixel into 4 subpixels, returns a new tracking +//! direction +template +bool Curve_renderer_2::subdivide( + Pixel_2& pix, int back_dir, int& new_dir) +{ + //Gfx_OUT("\n\nSubdivision of a pixel" << pix << " with back_dir = " << + //back_dir << std::endl); + NT inv = NT(1) / NT(one << pix.level); + make_exact(inv); + int idx, pref_dir = back_dir>>1; + + if(pix.level >= MAX_SUBDIVISION_LEVEL) { + std::cerr << "reached maximum subdivision level: " << pix << + std::endl; + throw CGALi::Insufficient_rasterize_precision_exception(); + } + + if(limit(engine.pixel_w * inv) || limit(engine.pixel_h * inv)) { + std::cerr << "too small subpixel size: " << pix << std::endl; + throw CGALi::Insufficient_rasterize_precision_exception(); + } + + // if several branches coincide withing this pixel we cannot perform + // a subdivision + if(branches_coincide) + return false; + + if(back_dir&1) { // diagonal direction + idx = get_subpixel_diag(pix,pref_dir); + if(idx == -1) { + std::cerr << "wrong diag subpixel: " << pix << " direction: " << + pref_dir << std::endl; + throw CGALi::Insufficient_rasterize_precision_exception(); + } + idx = CGALi::D_SUBPIX_MAP[pref_dir][idx]; + + } else { + idx = get_subpixel_hv(pix,pref_dir); + if(idx == -1) { + std::cerr << "wrong h/v subpixel" << pix << " direction: " << + pref_dir << std::endl; + throw CGALi::Insufficient_rasterize_precision_exception(); + } + idx = CGALi::HV_SUBPIX_MAP[pref_dir][idx]; + } + + pix.level++; + if(max_level < pix.level) + max_level = pix.level; + if(current_level < pix.level) + current_level = pix.level; + + pix.sub_x = (pix.sub_x<<1) + (idx&1); + pix.sub_y = (pix.sub_y<<1) + (idx>>1); + //Gfx_DETAILED_OUT("subpixel index: " << idx << " (" << pix.sub_x << "; " + // << pix.sub_y << ")" << std::endl); + if(!test_neighbourhood(pix, back_dir, new_dir)) + return subdivide(pix,back_dir,new_dir); + //Gfx_DETAILED_OUT("new direction found: " << new_dir << " at a pixel:" << + //pix << std::endl); + return true; +} + +//! returns a "seed" point of a segment and two possible directions from it +template +bool Curve_renderer_2::get_seed_point( + const Rational& seed, int arcno, Pixel_2& start, int *dir, + int *b_taken, bool& b_coincide) +{ + Rational x_s = seed, y_s; + + Xy_coordinate_2 xy(X_coordinate_1(seed), *support, arcno); + refine_xy(xy, engine.pixel_h_r/CGAL_REFINE_Y); + y_s = ubound_y(xy); + + Integer lvl = one; + Rational x_seed, y_seed; + get_pixel_coords(x_s, y_s, start, &x_seed, &y_seed); + + Gfx_OUT("y_seed = " << rat2float(y_seed) << std::endl); + + // we allow a small tolerance for seed point coordinates due to + // round off errors + if(start.y < -CGAL_WINDOW_ENLARGE || + start.y > engine.res_h + CGAL_WINDOW_ENLARGE) { + Gfx_OUT("get_seed_point: starting pixel does not fit the window " + "boundaries " << std::endl); + return false; + } + + start.level = 0; + start.sub_x = 0; + start.sub_y = 0; + current_level = 0; + + //Gfx_OUT("refining starting pixel: " << start << std::endl); +// Gfx_OUT("x/y seed: (" << NiX::to_double(x_seed) << "; " << +// NiX::to_double(y_seed) << ")" << std::endl); + b_taken[0] = b_taken[1] = -1; + b_coincide = false; + + int coincide_level = 0; + while(!test_pixel(start, (int *)dir, (int *)b_taken, b_coincide)) { + + //Gfx_OUT("refining starting pixel: " << start << std::endl); + if(start.level >= MAX_SUBDIVISION_LEVEL) { + std::cerr << "get_seed_point: reached maximum subdivision level " + << start.level << std::endl; + throw CGALi::Insufficient_rasterize_precision_exception(); + } + //dump_neighbourhood(start); + + if(limit(engine.pixel_w/NT(lvl))||limit(engine.pixel_h/NT(lvl))) { + std::cerr << "get_seed_point: too small subpixel size: " << + start.level << std::endl; + throw CGALi::Insufficient_rasterize_precision_exception(); + } + + start.level++; + if(start.level > max_level) + max_level = start.level; + if(start.level > current_level) + current_level = start.level; + lvl <<= 1; + + if(lvl > CGAL_REFINE_Y) { + refine_xy(xy, engine.pixel_h_r/(lvl*2)); + y_s = ubound_y(xy); +// event.refine_to(arcno, engine.pixel_h_r/(lvl*2)); + // y_s = event.upper_boundary(arcno); + + y_seed = (y_s - engine.y_min_r)/engine.pixel_h_r; +// NiX::simplify(y_seed); + } + + start.sub_x = rat2integer((x_seed - start.x)*lvl); + start.sub_y = rat2integer((y_seed - start.y)*lvl); + /*if(start.sub_x >= lvl || start.sub_x < 0 || + start.sub_y >= lvl || start.sub_y < 0) + Gfx_OUT << "bad subpixel found" << std::endl;*/ + + if(start.sub_x < 0) + start.sub_x = 0; + start.sub_x &= (lvl-1); + if(start.sub_y < 0) + start.sub_y = 0; + start.sub_y &= (lvl-1); + + if(current_level >= CGAL_COINCIDE_LEVEL) { + + Pixel_2 test = start; + test.sub_x = start.sub_x >> (start.level - coincide_level); + test.sub_y = start.sub_y >> (start.level - coincide_level); + test.level = coincide_level; + + if(test_pixel(test, (int *)dir, (int *)b_taken, b_coincide)) { + start = test; + break; + } + coincide_level++; + } + //Gfx_DETAILED_OUT("\nTesting pixel: " << start << std::endl); + } + + //Gfx_OUT("directions found: " << dir[0] << "; " << dir[1] << "\n"); + int t0 = CGALi::DIR_TAKEN_MAP[dir[0]], t1 = CGALi::DIR_TAKEN_MAP[dir[1]]; + if(t0 != -1&&t0 == t1) { + Gfx_OUT("get_seed_point: one-side directions found: " << + dir[0] << " and " << dir[1] << std::endl); + } + + if(!branches_coincide) { + if(t0 == -1) { + if(t1 == -1) { + // vx = -df/dy; vy = df/dx + typename Renderer_internals::Coercion::Cast rccast; + Rat_coercion_type xcs = rccast(x_s), ycs = rccast(y_s); + + Rat_coercion_type + vx = -engine.substitute_xy(*(engine.rational_fy), x_s, y_s), + vy = engine.substitute_xy(*(engine.rational_fx), x_s, y_s); + + /* typename Renderer_traits::Convert_poly convert_poly; + Poly_2 ffx, ffy; + ffx = convert_poly(*(engine.rational_fx)); + ffy = convert_poly(*(engine.rational_fy)); + + typename CGAL::Coercion_traits::Cast + castr; + + Coeff xx1 = castr(x_s), yy1 = castr(y_s); + Coeff vvx = NiX::substitute_xy(ffy, xx1, yy1); + Coeff vvy = NiX::substitute_xy(ffx, xx1, yy1); */ + + typename CGAL::Coercion_traits::Cast + cast; + +// Gfx_OUT("poly: " << *(engine.rational_fy) << +// "\n\nat: " << xcs << "; and " << ycs << "\n"); +// +// Gfx_OUT("vx before: " << cast(vx) << " " << (vy > 0) << +// "\nvvx before: " << vvx << " " << (vvy > 0) << "\n"); + +// vvx = -vvx; +// vx = -vx; + +// Gfx_OUT("vx after: " << cast(vx) << " " << (vy > 0) << "\nvy = " << +// (vy) << "\n"); +// Gfx_OUT("vvx after: " << vvx << " " << (vvy > 0) << "; vvy = " << vvy << "\n"); + + NT vvx = rat2float(vx), vvy = rat2float(vy); + //Coeff vvx = cast(vx), vvy = cast(vy); + if(vvy < 0) { + vvx = -vvx; + } + + Gfx_OUT("final vvx = " << vvx << " " << (vvx > 0) << "\n"); + + // vx > 0: 2 - right(1), 6 - left(0) + // vx < 0: 2 - left(0), 6 - right(1) + bool flag = (vvx > 0); + if(dir[0] == 2||dir[1] == 6) { + // taken 2 = 0 if vx>0, 1 if vx<0 + // taken 6 = 1 if vx>0, 0 if vx<0 + b_taken[0] = !flag; // !(vx > 0) + b_taken[1] = flag; // vx > 0 + } else if(dir[0] == 6||dir[1] == 2) { + b_taken[0] = flag; // vx > 0 + b_taken[1] = !flag; // !(vx > 0) + } + } else // t1 != -1 + b_taken[0] = 1 - b_taken[1]; + } else if(t1 == -1) + b_taken[1] = 1 - b_taken[0]; + + s_stack.push(Seed_point(start, dir[0], 1, b_taken[0], b_coincide)); + s_stack.push(Seed_point(start, dir[1], 0, b_taken[1], b_coincide)); + } + if(b_coincide) + Gfx_DETAILED_OUT("seed point with coincide branches found" << + std::endl); + return true; +} + +//! checks whether only one curve branch crosses the pixel, required to +//! compute a starting witness pixel +template +bool Curve_renderer_2::test_pixel( + const Pixel_2& pix, int *dir, int *b_taken, bool& b_coincide) +{ + Stripe box[2]; // 0 - left-right stripe, 1 - bottom-top stripe + NT lvl = NT(one << pix.level), inv = NT(1) / lvl; + make_exact(inv); + get_boundaries(CGAL_X_RANGE, pix, box[1]); + get_boundaries(CGAL_Y_RANGE, pix, box[0]); + NT bottom = box[1].key[0], //top = box[1].key[1], + left = box[0].key[0], lower;//, right = box[0].key[1], lower; + // bottom(2), top(3), left(0), right(1) + int n_sign = 0, i, j, n_dir = 0, shift, n_local, new_dir; + get_polynomials(CGAL_Y_RANGE, box[0]); + Gfx_DETAILED_OUT("\ncomputing left/right sides " << std::endl); + + b_coincide = false; + int n_corner_dir = 0, corner_dir[] = {-1, -1}; + // process subsegments: left/right + for(i = 0; i < 2; i++) { + for(j = 0, n_local = 0, lower = bottom; j < 3; j++, lower += inv) + + if(get_range_1(CGAL_Y_RANGE, lower, lower+inv, box[0].key[i], + box[0].poly[i], 3) || engine.zero_bounds) { + + Gfx_DETAILED_OUT("intersection detected at subsegment: " << j + << " side: " << i << "; left(0), right(1), bottom(2), top(3)" + << std::endl); + + if(engine.zero_bounds) { + + bool is_corner = false; + int diff = float2int((lower - pix.y)*lvl - pix.sub_y); + int low_sign = engine.evaluate_generic(CGAL_Y_RANGE, lower, + box[0].key[i], box[0].poly[i]); + + Gfx_DETAILED_OUT("zero bounds detected, low_sign: " + << low_sign << "\n"); + + if(diff != 0) { + if(diff == 1) { + is_corner = (engine.evaluate_generic(CGAL_Y_RANGE, + lower+inv, box[0].key[i], box[0].poly[i]) == 0); + // in case there is intersection at lower point => + // it is counted as normal intersection + if(!is_corner && low_sign != 0) + continue; + } else { + is_corner = (low_sign == 0); + if(!is_corner) + continue; + } + + if(is_corner) { + // compute corner direction + diff = CGALi::DIR_MAP[i][diff+1]; + if(n_corner_dir >= 2) { + Gfx_DETAILED_OUT("too many corners\n"); + return false; + } + // corner dirs on vertical segs cannot be identified + corner_dir[n_corner_dir++] = diff; + Gfx_DETAILED_OUT("corner direction saved: " + << diff << "\n"); + } + } else if(low_sign != 0) + continue; + } + + n_local++; + if(n_local > 1) { + Gfx_DETAILED_OUT("more than 1 intersection along vertical side " + << i << std::endl); + return false; + } + // no need to check the derivative if already coincide mode + if(!b_coincide) + if(engine.first_der && !recursive_check(CGAL_Y_RANGE, + lower, lower+inv, box[0].key[i],box[0].poly[i])) { + Gfx_DETAILED_OUT("\nrecursive_check failed" << std::endl); + if(current_level < CGAL_COINCIDE_LEVEL) + return false; + b_coincide = true; + } + shift = float2int((lower - pix.y)*lvl - pix.sub_y); + new_dir = CGALi::DIR_MAP[i][shift+1]; + // left side (i = 0) - direction taken = 1 + // right side (i = 1) - direction taken = 0 + b_taken[n_dir] = 1 - i; + Gfx_DETAILED_OUT("new dir found: " << new_dir); + dir[n_dir++] = new_dir; + } + n_sign += n_local; + } + + get_polynomials(CGAL_X_RANGE,box[1]); + Gfx_DETAILED_OUT("computing bottom/top sides" << std::endl); + for(i = 0; i < 2; i++) { + for(j = 0, n_local = 0, lower = left; j < 3; j++, lower += inv) { + + if(get_range_1(CGAL_X_RANGE,lower,lower+inv,box[1].key[i], + box[1].poly[i], 3) || engine.zero_bounds) { + + Gfx_DETAILED_OUT("intersection detected at subsegment: " << j << + " side: " << i + 2 << "; left(0), right(1), bottom(2), top(3) " + << std::endl); + + if(engine.zero_bounds) { + + bool is_corner = false; + int diff = float2int((lower - pix.x)*lvl - pix.sub_x); + int low_sign = engine.evaluate_generic(CGAL_X_RANGE, lower, + box[1].key[i], box[1].poly[i]); + + Gfx_DETAILED_OUT("zero bounds detected, low_sign: " + << low_sign << "\n"); + + if(diff != 0) { + if(diff == 1) { + is_corner = (engine.evaluate_generic(CGAL_X_RANGE, + lower+inv, box[1].key[i], box[1].poly[i]) == 0); + // in case there is intersection at lower point => + // it is counted as normal intersection + if(!is_corner && low_sign != 0) + continue; + } else { + is_corner = (low_sign == 0); + if(!is_corner) + continue; + } + + if(is_corner) { + Gfx_DETAILED_OUT("corner direction detected\n"); + // compute corner direction + diff = CGALi::DIR_MAP[i+2][diff+1]; + if(n_corner_dir > 0 && (corner_dir[0] == diff || + corner_dir[1] == diff)) { + Gfx_DETAILED_OUT("corner_dir identified: " + << diff << "\n"); + continue; + } + } + } else if(low_sign != 0) + continue; + } + + n_local++; + n_sign++; + if(n_local > 1||n_sign > 2) { + Gfx_DETAILED_OUT("more than 1 intersection along horizontal " + "side " << i+2 << std::endl); + return false; + } + + if(!b_coincide) + if(engine.first_der && !recursive_check(CGAL_X_RANGE, + lower, lower+inv, box[1].key[i],box[1].poly[i])) { + Gfx_DETAILED_OUT("\nrecursive_check failed" << std::endl); + if(current_level < CGAL_COINCIDE_LEVEL) + return false; + b_coincide = true; + } + shift = float2int((lower - pix.x)*lvl - pix.sub_x); + new_dir = CGALi::DIR_MAP[i+2][shift+1]; + b_taken[n_dir] = 1 - CGALi::DIR_TAKEN_MAP[new_dir]; + + Gfx_DETAILED_OUT(" new dir found: " << new_dir << "\n"); + dir[n_dir++] = new_dir; + } + } + } + if(n_sign < 2) { + Gfx_DETAILED_OUT("ERROR: not enough intersections found" << + std::endl); + return false; + } + return true; + + /*Stripe box[2]; // 0 - left-right stripe, 1 - bottom-top stripe + NT lvl = NT(one << pix.level), inv = NT(1) / lvl, llow[2]; + make_exact(inv); + get_boundaries(CGAL_X_RANGE,pix,box[1]); + get_boundaries(CGAL_Y_RANGE,pix,box[0]); + //NT bottom = box[1].key[0], //top = box[1].key[1], + //left = box[0].key[0], right = box[0].key[1]; + NT lower; + // bottom(2), top(3), left(0), right(1) + int f_der[2]; + int idx[2], ibox, ikey, var, n_sign = 0, i, j, n_dir = 0, shift, + n_local, new_dir; + get_polynomials(CGAL_Y_RANGE,box[0]); + b_coincide = false; + for(i = 0; i < 4; i++) { + if(i == 2) + get_polynomials(CGAL_X_RANGE,box[1]); + var = (i < 2 ? CGAL_Y_RANGE : CGAL_X_RANGE); + ibox = i >> 1; + ikey = i & 1; + lower = box[1-ibox].key[0]; + for(j = 0, n_local = 0; j < 3; j++, lower += inv) + if(get_range_1(var,lower,lower+inv,box[ibox].key[ikey], + box[ibox].poly[ikey])||zero_bounds) { + Gfx_DETAILED_OUT("intersection detected at subsegment: " << j + << " side: " << i << "; left(0), right(1), bottom(2), top(3)" + << std::endl); + if(zero_bounds&&evaluate_modular(var, lower, + box[ibox].key[ikey], true) != 0) + continue; + n_local++; + n_sign++; + if(n_local > 1||n_sign > 2) { + Gfx_DETAILED_OUT("more than 1 intersection along side " << i + << std::endl); + return false; + } + f_der[n_dir] = engine.first_der; + llow[n_dir] = lower; + idx[n_dir++] = i; + } + } + if(n_sign < 2) { + Gfx_DETAILED_OUT("not enough intersections found" << std::endl); + return false; + } + for(i = 0; i < 2; i++) { + var = (idx[i] < 2 ? CGAL_Y_RANGE : CGAL_X_RANGE); + ibox = idx[i] >> 1; + ikey = idx[i] & 1; + lower = llow[i]; + if(f_der[i] == -1) { // the 1st derivative needs to be computed + get_range_1(var,lower,lower+inv,box[ibox].key[ikey], + box[ibox].poly[ikey],2); + f_der[i] = engine.first_der; + } + // no need to check the derivative if already coincide mode + if(!b_coincide) + if(f_der[i]&&!recursive_check(var,lower,lower+inv,box[ibox].key[ikey], + box[ibox].poly[ikey])) { + Gfx_DETAILED_OUT("\nrecursive_check failed" << std::endl); + if(current_level < CGAL_COINCIDE_LEVEL) + return false; + b_coincide = true; + } + if(var == CGAL_Y_RANGE) + shift = float2int((lower - pix.y)*lvl - pix.sub_y); + else + shift = float2int((lower - pix.x)*lvl - pix.sub_x); + new_dir = CGALi::DIR_MAP[idx[i]][shift+1]; + b_taken[i] = 1 - CGALi::DIR_TAKEN_MAP[new_dir]; + dir[i] = new_dir; + Gfx_DETAILED_OUT(" direction " << i << " found: " << new_dir << + " taken: " << b_taken[i] << std::endl); + } + return true;*/ +} + +template +void Curve_renderer_2::horizontal_clip() +{ + typedef typename CGAL::Fraction_traits F_traits; + typename F_traits::Numerator_type num; + typename F_traits::Denominator_type denom; + typename F_traits::Decompose decompose; + + top_clip.clear(); + btm_clip.clear(); + Rational l, r; + + typename CGAL::Polynomial_traits_d< Poly_dst_1 >::Make_square_free msf; + + decompose(engine.y_min_r, num, denom); + btm_poly = msf(support->polynomial_2().evaluate_homogeneous(num, denom)); + + decompose(engine.y_max_r, num, denom); + top_poly = msf(support->polynomial_2().evaluate_homogeneous(num, denom)); + + CGAL::CGALi::Bitstream_descartes + isolator_btm(btm_poly), isolator_top(top_poly); + + int n_roots = isolator_btm.number_of_real_roots(), i; + Rational criteria = engine.pixel_w_r/CGAL_REFINE_CLIP_POINTS; +// NiX::simplify(criteria); + + for(i = 0; i < n_roots; i++) { + l = isolator_btm.left_boundary(i); + r = isolator_btm.right_boundary(i); + refine_alg_point(l, r, btm_poly, criteria); + if(l > engine.x_max_r || r < engine.x_min_r) + continue; + btm_clip.push_back(Clip_point_entry(l,r)); + } + + n_roots = isolator_top.number_of_real_roots(); + for(i = 0; i < n_roots; i++) { + l = isolator_top.left_boundary(i), + r = isolator_top.right_boundary(i); + refine_alg_point(l, r, top_poly, criteria); + if(l > engine.x_max_r || r < engine.x_min_r) + continue; + top_clip.push_back(Clip_point_entry(l,r)); + } +} + +// low_pix, up_pix and y_clip define segment end-points in pixel space +template +void Curve_renderer_2::segment_clip_points( + const Rational& x_lower, const Rational& x_upper, + const Rational& y_lower, const Rational& y_upper, + const Rational& y_clip, const Poly_dst_1& poly, int arcno, + Clip_points& clip_points, index_vector& clip_indices) +{ + //typename Curve_analysis_2::Internal_curve_2::Event1_info event; + + int n_arcs, i, j = 0; + typename Clip_points::iterator it = clip_points.begin(); + Rational low, high, d, l, r; + Gfx_DETAILED_OUT("segment: (" << rat2float(x_lower) << "; " << + rat2float(x_upper) << "); arcno: " << arcno << "\n"); + + while(it != clip_points.end()) { + l = it->left; + r = it->right; + if(r > x_lower && l < x_upper) { + X_coordinate_1 alpha = (l == r ? X_coordinate_1(r) : + X_coordinate_1(poly, l, r)); +// if(l == r) +// Gfx_DETAILED_OUT("checking rational clip-point\n "); +// else { +// Gfx_DETAILED_OUT("checking [" << l << "; " << r << +// "] clip-point; poly = " << poly << +// "; support: " << support->polynomial_2() << " \n\n"); +// +// if(!CGAL::CGALi::is_square_free(poly)) +// std::cerr << "polynomial is not square-free!\n"; +// +// if(poly.sign_at(l) == poly.sign_at(r)) +// std::cerr << "signs are the same!!\n\n"; +// } + + // need precise comparisons in case of tight boundaries + if(l < x_lower && alpha.compare(x_lower) != ::CGAL::LARGER) { + it++; j++; + continue; + } + if(r > x_upper && alpha.compare(x_upper) != ::CGAL::SMALLER) { + it++; j++; + continue; + } + if(it->arcno == -1) { +// event = support->_internal_curve().event_info_at_x( +// alpha); + Status_line_1 sline = support->status_line_for_x(alpha); + n_arcs = sline.number_of_events(); + for(i = 0; i < n_arcs; i++) { + + Xy_coordinate_2 xy(alpha, *support, i); + low = lbound_y(xy); + high = ubound_y(xy); + + //TODO: no need to refine y-intervals: you need only to check whether + // polynomial vanishes at y_clip, since algebraic real you specify is exact + if((low < y_clip && high > y_clip)|| + low == y_clip||high == y_clip) { + //event.refine_to(i, engine.pixel_h_r/CGAL_REFINE_Y); + refine_xy(xy, engine.pixel_h_r/CGAL_REFINE_Y); + if(lbound_y(xy) <= y_clip && ubound_y(xy) >= y_clip) { + it->arcno = i; + it->alpha = alpha; + Gfx_DETAILED_OUT("clip point assigned to arcno: " + << i << "\n"); + break; + } + } + } + } + if(it->arcno == arcno) { + bool set = true; + Xy_coordinate_2 xy(it->alpha, *support, arcno); + + if(r - x_lower <= engine.pixel_w_r) { + d = ubound_y(xy) - y_lower; + set = (CGAL_ABS(d) > engine.pixel_h_r); + } else if(x_upper - r <= engine.pixel_w_r) { + d = ubound_y(xy) - y_upper; + set = (CGAL_ABS(d) > engine.pixel_h_r); + } + if(set) { + Gfx_DETAILED_OUT("clip-point found\n"); + clip_indices.push_back(j); // store point index + } + } + } + it++; j++; + } +} + +// mode = 0: criteria specifies the length of the refined interval +// mode = 1: criteria specifies a point that must be lie outside the interval +template +void Curve_renderer_2::refine_alg_point( + Rational& l, Rational& r, const Poly_dst_1& poly, + const Rational& criteria, int mode) +{ + CGAL::Sign eval_l, eval_m; + Rational mid; + eval_l = poly.sign_at(l); + while(1) { + if(mode == 0) { + if(r - l < criteria) + break; + } else if(l > criteria || r < criteria) + break; + mid = (l+r)/2; +// NiX::simplify(mid); + eval_m = poly.sign_at(mid); + if(eval_m == CGAL::EQUAL) + l = r = mid; + else if(eval_m == eval_l) + l = mid; + else + r = mid; + } +} + +//! recursively checks whether only one curve branch crosses a line segment +template +bool Curve_renderer_2::recursive_check( + int var, const NT& beg_, const NT& end_, const NT& key, const Poly_1& poly, + int depth) +{ + // if polynomial is linear/quadratic and there is sign-change over interval + // => only one curve branch is guaranteed + if(poly.degree() < 3) + return true; + + int val_1, val_2, val_3; + NT mid = (beg_+end_)/2, key_1, key_2, beg = beg_, end = end_; + make_exact(mid); + + Gfx_DETAILED_OUT("executing recursive check" << std::endl); + val_1 = engine.evaluate_generic(var, beg, key, poly); + val_3 = engine.evaluate_generic(var, end, key, poly); + + Gfx_DETAILED_OUT("beg: " << val_1 << " end: " << val_3 << std::endl); + // no sing change: at least two curve branches inside + if((val_1^val_3)==0) + return false; + + val_2 = engine.evaluate_generic(var, mid, key, poly); + // select the half with even number of intersections + if((val_1^val_2)==0) + key_1 = beg; + else + key_1 = end; + if(get_range_1(var, key_1, mid, key, poly)) + return false; // at least three intersections over the interval: done + key_2 = beg + end - key_1; + + get_range_1(var, key_2, mid, key, poly, 3); + + // the first derivative does not straddle zero then only one intersection + if(!engine.first_der) + return true; + // if second derivative does not stranddle zero - only 1 first derivative + //if(engine.second_der == 0) + //return true; + if(depth > CGAL_DER_CHECK_DEPTH) + return true; + return recursive_check(var, key_2, mid, key, poly, depth+1); +} + +//! computes lower/upper boundaries for pixel's neighbourhood +template +void Curve_renderer_2::get_boundaries( + int var, const Pixel_2& pix, Stripe& stripe) +{ + int level = pix.level, val = pix.y; + Integer sub = pix.sub_y, cur_sub; + if(var == CGAL_Y_RANGE) { + sub = pix.sub_x; + val = pix.x; + } + cur_sub = sub; + if(level > 0) { + cur_sub += 2; // obtain local coordinates for upper boundary + // for even boundaries raise the level up + while(level > 0 && ((cur_sub & 1) == 0)) { + level--; + cur_sub >>= 1; + } + } + stripe.level[1] = level; + if(level == 0) + cur_sub = 2 - cur_sub; + stripe.key[1] = val + NT(cur_sub) / NT(one << level); + make_exact(stripe.key[1]); + level = pix.level; + cur_sub = sub - 1; + while(level > 0 && ((cur_sub & 1) == 0)) { + level--; + cur_sub >>= 1; + } + stripe.level[0] = level; + stripe.key[0] = val + NT(cur_sub) / NT(one << level); + make_exact(stripe.key[0]); +} + +//! shortcut to get precached polynomials +template + void Curve_renderer_2::get_polynomials( + int var, Stripe& stripe) +{ + engine.get_precached_poly(var, stripe.key[1], stripe.level[1], + stripe.poly[1]); + engine.get_precached_poly(var, stripe.key[0], stripe.level[0], + stripe.poly[0]); +} + +//! \brief checks 8-pixel neighbourhood of a pixel, returns \c true if +//! only one curve branch intersects pixel's neighbourhood, \c dir +//! defines backward direction, \c new_dir is a new tracking direction +template +bool Curve_renderer_2::test_neighbourhood( + const Pixel_2& pix, int dir, int& new_dir, int step_x, int step_y) +{ + NT lvl = NT(one << pix.level); + NT inv = NT(1.0) / lvl; + make_exact(inv); + int back_dir = dir, idx[5], n_pts = 5, n_sign; + int shift, i, j, res, var, s_shift = 0, tmp; + // box[0].key[0] - x-coordinate of left pixel boundary + // box[0].key[1] - x-coordinate of right pixel boundary + // box[1].key[0] - y-coordinate of bottom pixel boundary + // box[1].key[1] - y-coordinate of top pixel boundary + Stripe box[2]; // 0 - left-right stripe, 1 - bottom-top stripe + NT p1, p2, lower; + get_boundaries(CGAL_X_RANGE,pix,box[1]); + get_boundaries(CGAL_Y_RANGE,pix,box[0]); + + struct { // local point descriptor + NT key; // respective "key" + int flag; // 0th bit: 0: lower polynomial(left/bottom); + // 1: upper polynomial(right/top) + // 1st bit: 0: x-range poly; 1: y-range poly + // 2nd bit: 0: c is in lower poly; 1: c is in upper poly + // (for regular points only) + } pts[6] = {{box[0].key[0], 0}, {box[0].key[0], 1}, {box[0].key[1], 1|4}, + {box[0].key[1], 0|4}}; + + int ix = CGALi::directions[back_dir].x; + int e_dir = back_dir, _3and5 = 0, _3and7 = 0; + if(back_dir&1) {// 5 and 3: 4; 7 and 1: 0 + _3and5 = (e_dir == 5||e_dir == 3); + _3and7 = e_dir&2; + e_dir = _3and5 ? 4: 0; + } + int _4and6 = e_dir>>2, _2and6 = (e_dir&2)>>1; + int _2and4 = (e_dir==2)||(e_dir==4); + + pts[4].key = box[_2and6].key[_4and6] + (1-(_4and6<<1))*inv; + pts[4].flag = (_2and6<<1) + _2and4; + pts[5].key = pts[4].key; + pts[5].flag = pts[4].flag^1; //(_2and6<<1) + (!_2and4); // pts[4].flag^1 ? + idx[0] = 4; + idx[3] = 5; + if(back_dir&1) + (_3and7) ? pts[5].key += ix*inv : pts[4].key += ix*inv; + shift = (6 - e_dir) >> 1; + n_pts = 4; + idx[1] = (1 + shift) & 3; + idx[2] = (2 + shift) & 3; +/* NT xx, yy; + if(var == CGAL_X_RANGE) { + xx = pts[idx[i]].key; + yy = box[shift].key[fl&1]; + } else { + yy = pts[idx[i]].key; + xx = box[shift].key[fl&1]; + } + xx = xx*20+ofsxx; + yy = yy*20+100; + if(i == 0) + ppnt->moveTo(int(to_double(xx)),700-int(to_double(yy))); + else + ppnt->lineTo(int(to_double(xx)),700-int(to_double(yy))); + prev = curr; + }*/ + + Gfx_DETAILED_OUT("test_neigh for " << pix << "; dir = " << dir << "\n"); + + n_sign = 0; + get_polynomials(CGAL_Y_RANGE,box[0]); + get_polynomials(CGAL_X_RANGE,box[1]); + int f_der = -1, corner_dir = -1; + bool set_coincide = false; + new_dir = -1; + for(i = 0; i < n_pts - 1; i++) { + int f1 = pts[idx[i]].flag, f2 = pts[idx[i+1]].flag; + var = CGAL_Y_RANGE; + p1 = pts[idx[i]].key; + p2 = pts[idx[i+1]].key; + if((f1&2) + (f2&2) == 0) {// both can use x-range polynomials + // different polys in x-range - use y-range instead + if((f1&1)!=(f2&1)) { + shift = (f1&4)>>2; + p1 = box[1].key[0]; + p2 = box[1].key[1]; + } else { // same polys in x-range + shift = f1&1; + var = CGAL_X_RANGE; + shift += 2; + } + } else { + if((f1&2) == 0) { // first point is in x-range (convert to y-range) + shift = f2&1; + p1 = box[1].key[f1&1]; + } else { // second point is in x-range + shift = f1&1; + p2 = box[1].key[f2&1]; + } + } + if(p1 > p2) { + NT tmp1 = p1; + p1 = p2; + p2 = tmp1; + } + res = float2int((p2 - p1)*lvl); + int local_sign = 0; + for(j = 0; j < res; j++, p1 += inv) { + + NT lkey = box[shift>>1].key[shift&1]; + const Poly_1& lpoly = box[shift>>1].poly[shift&1]; + + if(get_range_1(var, p1, p1+inv, lkey, lpoly, 3) || + engine.zero_bounds) { + + Gfx_DETAILED_OUT("sign-change found at: " << i << "; subseg: " + << j << std::endl); + + if(engine.zero_bounds) { + Gfx_DETAILED_OUT("also with zero_bounds"); + + bool is_corner = false; + int diff = float2int(var == CGAL_X_RANGE ? + ((p1 - pix.x)*lvl - pix.sub_x) : + ((p1 - pix.y)*lvl - pix.sub_y)); + int low_sign = engine.evaluate_generic(var, p1, lkey, + lpoly); + + if(diff != 0) { + // diff == -1: test lower corner: p1 + // diff == 1: test upper corver: p1+inv + if(diff == 1) { + is_corner = (engine.evaluate_generic(var, p1+inv, + lkey, lpoly) == 0); + // in case there is intersection at lower point => + // it is counted as normal intersection + if(!is_corner && low_sign != 0) + continue; + } else { + is_corner = (low_sign == 0); + if(!is_corner) + continue; + } + + if(is_corner) { + // compute corner direction + tmp = CGALi::DIR_MAP[shift][diff+1]; + Gfx_DETAILED_OUT("corner direction detected: " + << tmp << "\n"); + + if(corner_dir == -1) + corner_dir = tmp; + else if(corner_dir != tmp) { + Gfx_DETAILED_OUT("different corner \ + directions found\n"); + return false; + } else { + // decrement since sign change occurs in the + // corner + local_sign--; + Gfx_DETAILED_OUT("corner found: sign: " + << n_sign << "; local: " << local_sign + << "\n"); + } + } + } else if(low_sign != 0) + continue; + } + local_sign++; + // Gfx_DETAILED_OUT("a curve branch detected at interval: (" << i + // << "; " << i+1 << ") segment: " << j << std::endl); + if(local_sign > 1) { + Gfx_DETAILED_OUT("more than 1 curve branch detected" << + std::endl); + return false; + } + lower = p1; + s_shift = shift; + f_der = engine.first_der; + } + } + n_sign += local_sign; + if(local_sign == 1) { + tmp = float2int(var == CGAL_X_RANGE ? // bottom/top + ((lower - pix.x)*lvl - pix.sub_x) : + ((lower - pix.y)*lvl - pix.sub_y)); + tmp = CGALi::DIR_MAP[s_shift][tmp+1]; + /*if(new_dir != -1&&new_dir != tmp) // more than 1 direction found + return false; + else { + if(!branches_coincide&& + (current_level < CGAL_COINCIDE_LEVEL)) + return false; // level is too small + set_coincide = true; // diagonal coincide mode + }*/ + new_dir = tmp; + if(n_sign > 1) { + Gfx_DETAILED_OUT("test_neigh: too many intersections found"); + return false; + } + } + } + + //Gfx_DETAILED_OUT("processing segment: [" << lower << "; " << lower+inv << + // "]" <>1, ikey = s_shift&1; + if(f_der == -1) { // need to compute the first derivative + get_range_1(var, lower, lower+inv, box[ibox].key[ikey], + box[ibox].poly[ikey], 3); + f_der = engine.first_der; + } + // if coincide already set - no need to check the derivative + if(!set_coincide && f_der && !recursive_check(var,lower,lower+inv, + box[ibox].key[ikey], box[ibox].poly[ikey])) { + + //Gfx_DETAILED_OUT("first derivative presented" << std::endl); + + if(!branches_coincide && + (current_level < CGAL_COINCIDE_LEVEL))//||direction_taken == -1)) + return false; + if(!branches_coincide) + Gfx_DETAILED_OUT("Branches coincide at pixel: " << pix << + std::endl); + set_coincide = true; + } + int taken = CGALi::DIR_TAKEN_MAP[new_dir]; + if(taken != -1&&taken != direction_taken) { + + Gfx_DETAILED_OUT("\nERROR: wrong direction " << new_dir << " at pixel: " + << pix << "; back_dir = " << back_dir << "; taken = " << + direction_taken << std::endl); + if(back_dir == 2) + new_dir = 6; + else if(back_dir == 6) + new_dir = 2; + else + return false; + } + if(set_coincide) + branches_coincide = true; + return true; +} + + +//! \brief returns whether a polynomial has zero at a given interval, +//! we are not interested in concrete values +//! +//! if \t der_check is set a range for the first derivative is also computed +template +inline bool Curve_renderer_2:: + get_range_1(int var, const NT& lower, const NT& upper, const NT& key, + const Poly_1& poly, int check) +{ + bool res = engine.get_range_QF_1(var, lower, upper, key, poly, check); + //engine.get_range_MAA_1(var, lower, upper, key, poly, check); + return res; +} + +//! \brief copmutes an isolating interval for y-coordinate of an end-point, +//! uses caching if possible +template +typename Curve_renderer_2::Rational +Curve_renderer_2::get_endpoint_y( + const Arc_2& arc, const Rational& x, CGAL::Arr_curve_end end, + bool is_clipped) { + +// typename Curve_analysis_2::Internal_curve_2::Event1_info event; + int arcno; + Xy_coordinate_2 xy; + + if(arc.location(end) != CGAL::ARR_LEFT_BOUNDARY && + arc.location(end) != CGAL::ARR_RIGHT_BOUNDARY && !is_clipped) { + //event = arc.curve_end(end).curve()._internal_curve(). + // event_info_at_x(arc.curve_end_x(end)); + arcno = arc.curve_end(end).arcno(); + xy = Xy_coordinate_2(arc.curve_end_x(end), + arc.curve_end(end).curve(), arcno); + } else { + //sline = support->status_line_for_x(X_coordinate_1(x)); + //event = support->_internal_curve().event_info_at_x( + // X_coordinate_1(x)); + arcno = arc.arcno(); + xy = Xy_coordinate_2(X_coordinate_1(x), *support, arcno); + } + + // refine the y-interval until its size is smaller than the pixel size + ////////////////////////// CHANGED REFINE_Y by REFINE_X + //event.refine_to(arcno, engine.pixel_h_r/CGAL_REFINE_X); + refine_xy(xy, engine.pixel_h_r/CGAL_REFINE_X); + ////////////////////////// CHANGED REFINE_Y by REFINE_X + return ubound_y(xy);//event.upper_boundary(arcno); +} + +//! \brief switches to a certain cache instance depending on currently used +//! algebraic curve +template +void Curve_renderer_2::select_cache_entry( + const Arc_2& arc) +{ + typename boost::multi_index::nth_index::type& + idx = cache_list.get<1>(); + typename boost::multi_index::nth_index_iterator::type + it = idx.find(arc.curve().id()); + + int new_id; + if(it == idx.end()) { + new_id = cache_list.size(); + if(new_id >= CGAL_N_CACHES) { + new_id = cache_list.back().second; + cache_list.pop_back(); + } + cache_list.push_front(LRU_entry(arc.curve().id(),new_id)); + } else { // mark this element as most recently used + new_id = (*it).second; + cache_list.relocate(cache_list.begin(), cache_list.project<0>(it)); + } + + if(cache_id != new_id) + clip_pts_computed = false; + cache_id = new_id; + + support = support_ + cache_id; + engine.select_cache_entry(cache_id); + + while(!s_stack.empty()) + s_stack.pop(); + + if(it == idx.end()) { + *support = arc.curve(); + engine.precompute(support->polynomial_2()); + } +} + +//! \brief returns one of 4 subpixels of a pixel which is crossed by the curve +//! branch +//! +//! subpixels are chosen with the priority \c dir which defines a preferred +//! direction, i.e. right(0), top(1), left(2) or bottom(3) this is only for h/v +//! directions +template +int Curve_renderer_2::get_subpixel_hv( + const Pixel_2& pix, int dir) +{ + Poly_1 box[4]; + NT inv = NT(1) / NT(one << pix.level), inv_2; + make_exact(inv); + inv_2 = inv/2; + make_exact(inv_2); + NT pix_x = pix.x + pix.sub_x * inv; + NT pix_y = pix.y + pix.sub_y * inv; + NT s_c, s_key, inc_c, inc_key; + int var, inv_var; + if(dir&1) { // 1 or 3 (vertical direction) + var = CGAL_X_RANGE; + s_c = pix_x; + s_key = pix_y; + } else { // 0 or 2 (horizontal direction) + var = CGAL_Y_RANGE; + s_c = pix_y; + s_key = pix_x; + } + inc_key = ((dir&2)-1)*inv_2; + inc_c = (1-2*((dir>>1)^(dir&1)))*inv_2; + inv_var = 1 - var; + s_c = s_c + inv_2 - inc_c; + s_key = s_key + inv_2 - inc_key; + engine.get_precached_poly(var, s_key, pix.level, box[0]); + + //Gfx_DETAILED_OUT("hv subdpixel" << std::endl); + int p0 = engine.evaluate_generic(var, s_c, s_key, box[0]); // 0 + int p1 = engine.evaluate_generic(var, s_c + inc_c, s_key, box[0]); // 1 + if(p0^p1) + return 0; + + int p2 = engine.evaluate_generic(var, s_c + inc_c*2, s_key, box[0]); // 2 + if(p2^p1) + return 1; + + engine.get_precached_poly(inv_var, s_c, pix.level, box[1]); + int p3 = engine.evaluate_generic(inv_var, s_key+inc_key, s_c, box[1]); // 3 + if(p3^p0) + return 0; + + engine.get_precached_poly(inv_var, s_c + inc_c*2, pix.level, box[2]); + int p4 = engine.evaluate_generic(inv_var, s_key+inc_key, s_c + inc_c*2, + box[2]); + + if(p4^p2) + return 1; + + int p5 = engine.evaluate_generic(inv_var, s_key+inc_key*2, s_c, box[1]); + if(p3^p5) + return 2; + + int p7 = engine.evaluate_generic(inv_var, s_key+inc_key*2, s_c + inc_c*2, + box[2]); // 7 + if(p4^p7) + return 3; + + engine.get_precached_poly(var, s_key+inc_key*2, pix.level, box[3]); + int p6 = engine.evaluate_generic(var, s_c + inc_c, s_key+inc_key*2, + box[3]); // 6 + if(p5^p6) + return 2; + if(p6^p7) + return 3; + return -1; +} + +//! \brief the same for diagonal directions +//! +//! preferred direction: NE(0), NW(1), SW(2), SE(3) +template +int Curve_renderer_2::get_subpixel_diag( + const Pixel_2& pix, int dir) +{ + Poly_1 box[4]; + NT inv = NT(1) / NT(one << pix.level); + make_exact(inv); + NT inv_2 = inv/2; + make_exact(inv_2); + NT pix_x = pix.x + pix.sub_x * inv; + NT pix_y = pix.y + pix.sub_y * inv; + NT key_1, key_2, inc_1, inc_2; + int var, inv_var; + if(dir&1) { // 1 or 3 (vertical direction) + var = CGAL_X_RANGE; + key_1 = pix_x; + key_2 = pix_y; + } else { // 0 or 2 (horizontal direction) + var = CGAL_Y_RANGE; + key_1 = pix_y; + key_2 = pix_x; + } + inc_1 = (2*((dir>>1)^(dir&1))-1)*inv_2; + inc_2 = ((dir&2)-1)*inv_2; + inv_var = 1 - var; + key_1 = key_1 + inv_2 - inc_1; + key_2 = key_2 + inv_2 - inc_2; + + engine.get_precached_poly(var, key_2, pix.level, box[0]); + int p0 = engine.evaluate_generic(var, key_1, key_2, box[0]); // 0 + int p1 = engine.evaluate_generic(var, key_1 + inc_1, key_2, box[0]); // 1 + //Gfx_DETAILED_OUT("p0: " << p0 << "; p1: " << p1 << std::endl); + + if(p0^p1) + return 0; + + engine.get_precached_poly(inv_var, key_1, pix.level, box[1]); + int p2 = engine.evaluate_generic(inv_var, key_2+inc_2, key_1, box[1]); // 2 + //Gfx_DETAILED_OUT << "p2: " << p2 << std::endl; + if(p2^p0) + return 0; + + int p3 = engine.evaluate_generic(var, key_1 + inc_1*2, key_2, box[0]); // 3 + //Gfx_DETAILED_OUT << "p3: " << p3 << std::endl; + if(p3^p1) + return 1; + + int p4 = engine.evaluate_generic(inv_var, key_2+inc_2*2, key_1, box[1]); + //Gfx_DETAILED_OUT << "p4: " << p4 << std::endl; + if(p4^p2) + return 2; + + engine.get_precached_poly(inv_var, key_1 + inc_1*2, pix.level, box[2]); + int p5 = engine.evaluate_generic(inv_var, key_2+inc_2, key_1 + inc_1*2, + box[2]); + //Gfx_DETAILED_OUT << "p5: " << p5 << std::endl; + if(p3^p5) + return 1; + + engine.get_precached_poly(var, key_2+inc_2*2, pix.level, box[3]); + int p6 = engine.evaluate_generic(var, key_1 + inc_1, key_2+inc_2*2, + box[3]); // 6 + //Gfx_DETAILED_OUT << "p6: " << p6 << std::endl; + if(p4^p6) + return 2; + + int p7 = engine.evaluate_generic(var, key_1+inc_1*2, key_2 + inc_2*2, + box[3]); // 7 + //Gfx_DETAILED_OUT << "p7: " << p7 << std::endl; + if(p5^p7) + return 3; + + if(p6^p7) + return 3; + return -1; +} + +//! attempts to find an isolating box for an end-point whose upper/lower +//! boundaries do not intersect with a curve +template +bool Curve_renderer_2::get_isolating_box( + const Rational& x_s, const Rational& y_s, Pixel_2& res) +{ + Integer lvl = one; + Rational x_seed, y_seed; + get_pixel_coords(x_s, y_s, res, &x_seed, &y_seed); + res.level = 0; + res.sub_x = 0; + res.sub_y = 0; + + NT inv = NT(1), bottom = NT(res.y), top = bottom + inv, left = NT(res.x); + + Poly_1 poly_btm, poly_top; + engine.get_precached_poly(CGAL_X_RANGE, bottom, 0, poly_btm); + engine.get_precached_poly(CGAL_X_RANGE, top, 0, poly_top); + + while(1) { + inv = NT(1) / NT(lvl); + make_exact(inv); + left = NT(res.x) + NT(res.sub_x) * inv; + + Gfx_DETAILED_OUT("processing pixel: " << res << std::endl); + if(get_range_1(CGAL_X_RANGE,left,left+inv,bottom,poly_btm)&& + !engine.zero_bounds) { + + Gfx_DETAILED_OUT("bottom boundary intersection found at pixel: " << + res << std::endl); + } else { + if(!get_range_1(CGAL_X_RANGE,left,left+inv,top,poly_top)) { + + Gfx_DETAILED_OUT("no bottom/top intersections found at pixel:" + << res << std::endl); + break; + } else { + Gfx_DETAILED_OUT("top boundary intersection found at pixel: " + << res << std::endl); + } + } + lvl <<= 1; + /*if(lvl > CGAL_REFINE_Y) { + event.refine_to(arcno, pixel_h_r/(lvl*2)); + y_seed = (event.lower_boundary(arcno) + + event.upper_boundary(arcno))/2; + }*/ + if(res.level >= MAX_SUBDIVISION_LEVEL) { + std::cerr("get_isolating_box: reached maximum subdivision level " + << res.level << std::endl); + return false; + } + res.level++; + res.sub_x = rat2integer((x_seed - res.x)*lvl); + res.sub_y = rat2integer((y_seed - res.y)*lvl); + if(res.sub_x < 0) + res.sub_x = 0; + res.sub_x &= (lvl-1); + if(res.sub_y < 0) + res.sub_y = 0; + res.sub_y &= (lvl-1); + } + return true; +} + +//! returns true if \c pix is encompassed into one of isolating rectangles +//! (stopping criteria) +template +inline bool +Curve_renderer_2::is_isolated_pixel( + const Pixel_2& pix) +{ + /*Integer sub_x; + if(isolated_l.level != -1u&&isolated_l.y == pix.y&&isolated_l.x == pix.x&& + pix.level >= isolated_l.level) { + sub_x = pix.sub_x >> (pix.level - isolated_l.level); + if(CGAL_ABS(sub_x - isolated_l.sub_x) <= 1) + return true; + } + if(isolated_h.level != -1u&&isolated_h.y == pix.y&&isolated_h.x == pix.x&& + pix.level >= isolated_h.level) { + sub_x = pix.sub_x >> (pix.level - isolated_h.level); + if(CGAL_ABS(sub_x - isolated_h.sub_x) <= 1) + return true; + }*/ + return false; +} + +#ifdef Gfx_DEBUG_PRINT +// DEBUG ONLY +template +void Curve_renderer_2::dump_neighbourhood( + const Pixel_2& pix) +{ + CGAL::set_mode(std::cerr, CGAL::IO::PRETTY); + CGAL::set_mode(std::cout, CGAL::IO::PRETTY); + + Stripe box[2]; // 0 - left-right stripe, 1 - bottom-top stripe + //NT inv = NT(1) / NT(one << pix.level); + get_boundaries(CGAL_X_RANGE,pix,box[1]); + get_boundaries(CGAL_Y_RANGE,pix,box[0]); + NT bottom = box[1].key[0], top = box[1].key[1], + left = box[0].key[0], right = box[0].key[1]; + + Gfx_OUT("\n\nComputing nighbourhood for pixel: " << pix << std::endl); + get_polynomials(CGAL_X_RANGE, box[1]); + get_polynomials(CGAL_Y_RANGE, box[0]); + NT a,b,range,inc,val; + /////////////////////////////////////////////////////////////////////// + Gfx_OUT("\nevaluate RIGHT side:" << std::endl); + range = top - bottom; + inc = range / 3; + val = bottom; + + if(get_range_1(CGAL_Y_RANGE, bottom, top, right, box[0].poly[1])) + Gfx_OUT("segment RIGHT registered" << std::endl); + if(get_range_1(CGAL_Y_RANGE, val, val + inc, right, box[0].poly[1])) + Gfx_OUT("segment 0 registered" << std::endl); + + a = engine.evaluate_generic(CGAL_Y_RANGE, val, right, box[0].poly[1]); + Gfx_OUT("val = " << a << std::endl); + val += inc; + + b = engine.evaluate_generic(CGAL_Y_RANGE, val, right, box[0].poly[1]); + Gfx_OUT("val = " << b << std::endl); + if(a*b < 0) + Gfx_OUT("sign change at segment 0" << std::endl); + a = b; + if(get_range_1(CGAL_Y_RANGE, val, val + inc, right, box[0].poly[1])) + Gfx_OUT("segment 1 registered" << std::endl); + val += inc; + b = engine.evaluate_generic(CGAL_Y_RANGE, val, right, box[0].poly[1]); + + Gfx_OUT("val = " << b << std::endl); + if(a*b < 0) + Gfx_OUT("sign change at segment 1" << std::endl); + if(get_range_1(CGAL_Y_RANGE, val, val + inc, right, box[0].poly[1])) + Gfx_OUT("segment 2 registered" << std::endl); + a = b; + val = top; + + Gfx_OUT("poly right: " << box[0].poly[1] << " at (" << val << "; " << + right << ")\n"); + + b = engine.evaluate_generic(CGAL_Y_RANGE, val, right, box[0].poly[1]); + + Gfx_OUT("val = " << b << std::endl); + if(a*b < 0) + Gfx_OUT("sign change at segment 2" << std::endl); + /////////////////////////////////////////////////////////////////////// + Gfx_OUT("\nevaluate BOTTOM side:" << std::endl); + range = right - left; + inc = range / 3; + val = left; + if(get_range_1(CGAL_X_RANGE, left, right, bottom, box[1].poly[0])) + Gfx_OUT("segment BOTTOM registered" << std::endl); + if(get_range_1(CGAL_X_RANGE, val, val + inc, bottom, box[1].poly[0])) + Gfx_OUT("segment 0 registered" << std::endl); + + Gfx_OUT("poly bottom: " << box[1].poly[0] << " at (" << val << "; " << + bottom << ")\n"); + + engine.show_dump = true; + a = engine.evaluate_generic(CGAL_X_RANGE, val, bottom, box[1].poly[0]); + engine.show_dump = false; + + Gfx_OUT("val = " << a << std::endl); + val += inc; + b = engine.evaluate_generic(CGAL_X_RANGE, val, bottom, box[1].poly[0]); + + Gfx_OUT("val = " << b << std::endl); + if(a*b < 0) + Gfx_OUT("sign change at segment 0" << std::endl); + a = b; + if(get_range_1(CGAL_X_RANGE, val, val + inc, bottom, box[1].poly[0])) + Gfx_OUT("segment 1 registered" << std::endl); + val += inc; + b = engine.evaluate_generic(CGAL_X_RANGE, val, bottom, box[1].poly[0]); + + Gfx_OUT("val = " << b << std::endl); + if(a*b < 0) + Gfx_OUT("sign change at segment 1" << std::endl); + a = b; + if(get_range_1(CGAL_X_RANGE, val, val + inc, bottom, box[1].poly[0])) + Gfx_OUT("segment 2 registered" << std::endl); + val = right; + b = engine.evaluate_generic(CGAL_X_RANGE, val, bottom, box[1].poly[0]); + + Gfx_OUT("val = " << b << std::endl); + if(a*b < 0) + Gfx_OUT("sign change at segment 2" << std::endl); + /////////////////////////////////////////////////////////////////////// + Gfx_OUT("\nevaluate LEFT side:" << std::endl); + range = top - bottom; + inc = range / 3; + val = bottom; + if(get_range_1(CGAL_Y_RANGE, bottom, top, left, box[0].poly[0])) + Gfx_OUT("segment LEFT registered" << std::endl); + if(get_range_1(CGAL_Y_RANGE, val, val + inc, left, box[0].poly[0])) + Gfx_OUT("segment 0 registered" << std::endl); + + Gfx_OUT("poly left: " << box[0].poly[0] << " at (" << val << "; " << + left << ")\n"); + engine.show_dump = true; + a = engine.evaluate_generic(CGAL_Y_RANGE, val, left, box[0].poly[0]); + engine.show_dump = false; + + + Gfx_OUT("val = " << a << std::endl); + val += inc; + b = engine.evaluate_generic(CGAL_Y_RANGE, val, left, box[0].poly[0]); + + Gfx_OUT("val = " << b << std::endl); + if(a*b < 0) + Gfx_OUT("sign change at segment 0" << std::endl); + a = b; + if(get_range_1(CGAL_Y_RANGE, val, val + inc, left, box[0].poly[0])) + Gfx_OUT("segment 1 registered" << std::endl); + val += inc; + b = engine.evaluate_generic(CGAL_Y_RANGE, val, left, box[0].poly[0]); + + Gfx_OUT("val = " << b << std::endl); + if(a*b < 0) + Gfx_OUT("sign change at segment 1" << std::endl); + a = b; + if(get_range_1(CGAL_Y_RANGE, val, val + inc, left, box[0].poly[0])) + Gfx_OUT("segment 2 registered" << std::endl); + val = top; + b = engine.evaluate_generic(CGAL_Y_RANGE, val, left, box[0].poly[0]); + + Gfx_OUT("val = " << b << std::endl); + if(a*b < 0) + Gfx_OUT("sign change at segment 2" << std::endl); + /////////////////////////////////////////////////////////////////////// + Gfx_OUT("\nevaluate TOP side:" << std::endl); + range = right - left; + inc = range / 3; + val = left; + if(get_range_1(CGAL_X_RANGE, left, right, top, box[1].poly[1])) + Gfx_OUT("segment TOP registered" << std::endl); + + if(get_range_1(CGAL_X_RANGE, val, val + inc, top, box[1].poly[1])) + Gfx_OUT("segment 0 registered" << std::endl); + a = engine.evaluate_generic(CGAL_X_RANGE, val, top, box[1].poly[1]); + + Gfx_OUT("val = " << a << std::endl); + val += inc; + b = engine.evaluate_generic(CGAL_X_RANGE, val, top, box[1].poly[1]); + + Gfx_OUT("val = " << b << std::endl); + if(a*b < 0) + Gfx_OUT("sign change at segment 0" << std::endl); + a = b; + if(get_range_1(CGAL_X_RANGE, val, val + inc, top, box[1].poly[1])) + Gfx_OUT("segment 1 registered" << std::endl); + val += inc; + b = engine.evaluate_generic(CGAL_X_RANGE, val, top, box[1].poly[1]); + + Gfx_OUT("val = " << b << std::endl); + if(a*b < 0) + Gfx_OUT("sign change at segment 1" << std::endl); + a = b; + + if(get_range_1(CGAL_X_RANGE, val, val + inc, top, box[1].poly[1])) + Gfx_OUT("segment 2 registered" << std::endl); + + val = right; + + Gfx_OUT("poly top: " << box[1].poly[1] << " at (" << val << "; " << + top << ")\n"); + b = engine.evaluate_generic(CGAL_X_RANGE, val, top, box[1].poly[1]); + + Gfx_OUT("val = " << b << std::endl); + if(a*b < 0) + Gfx_OUT("sign change at segment 2" << std::endl); +} +#endif + +//!@} +CGAL_END_NAMESPACE + +#endif // CGAL_CKVA_CURVE_RENDERER_2_H diff --git a/Curved_kernel_via_analysis_2/include/CGAL/Curved_kernel_via_analysis_2/gfx/Curve_renderer_internals.h b/Curved_kernel_via_analysis_2/include/CGAL/Curved_kernel_via_analysis_2/gfx/Curve_renderer_internals.h new file mode 100644 index 00000000000..e3c271bafcc --- /dev/null +++ b/Curved_kernel_via_analysis_2/include/CGAL/Curved_kernel_via_analysis_2/gfx/Curve_renderer_internals.h @@ -0,0 +1,1072 @@ +// TODO: Add licence +// +// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE +// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. +// +// $URL:$ +// $Id: $ +// +// +// Author(s) : Pavel Emeliyanenko +// +// ============================================================================ + +/*!\file CGAL/Curved_kernel_via_analysis_2/gfx/Curve_renderer_internals.h + * \brief + * contains various low-level routines used by \c Curve_renderer_2 and + * \c Subdivision_1 + * + * provide caching for polynomials and polynomial evaluations; 1D range + * analysis using First Quadratic Affine Form (QF) and Modified Affine + * Arithmetic (MAA), both equipped with recursive derivative information + */ + +#ifndef CGAL_CKVA_CURVE_RENDERER_INTERNALS_H +#define CGAL_CKVA_CURVE_RENDERER_INTERNALS_H 1 + +#include +#include +#include +#include +#include +#include + +#include +#include + +#include + +using boost::multi_index::multi_index_container; +using boost::multi_index::get; +using boost::multi_index::project; + +CGAL_BEGIN_NAMESPACE + +namespace CGALi { + +#define CGAL_N_CACHES 2 // maximal number of cache instances + +#define CGAL_X_RANGE 0 // polynomial in x-range + +#define CGAL_Y_RANGE 1 // polynomial in y-range + +// defines subdivision level beyond which univariate polynomials are not +// cached (if 0 cache is off) +#define CGAL_MAX_POLY_CACHE_LEVEL 12 + +// maximal number of entries in univariate polynomial cache containers +#define CGAL_POLY_CACHE_SIZE 2*1024*1024 + +// defines subdivision level beyond which polynomial evaluations are not +// cached (if 0 cache is off) +#define CGAL_MAX_EVAL_CACHE_LEVEL 12 + +// maximal number of entries in polynomial evaluation cache container +#define CGAL_EVAL_CACHE_SIZE 3*1024*1024 + +// maximal degree of the derivative taken into account during recursive +// derivative range analysis +#define CGAL_RECURSIVE_DER_MAX_DEGREE 7 + +// 8-pixel neighbouthood directions +static const struct { int x; int y; } directions[] = { + { 1, 0}, { 1, 1}, { 0, 1}, {-1, 1}, + {-1, 0}, {-1,-1}, { 0,-1}, { 1,-1}}; +// map from rectangle sides to pixel directions +// 0 - left side, 1 - right side, 2 - bottom side, 3 - top side +static const int DIR_MAP[][3] = + {{5, 4, 3}, {7, 0, 1}, {5, 6, 7}, {3, 2, 1}}; +// map from 8 pixel directions to "left" (0) or "right" (1) direction of motion +// vertical directions (2 and 6) are mapped to -1 +static const int DIR_TAKEN_MAP[] = + {1, 1, -1, 0, 0, 0, -1, 1}; + +template +struct Pixel_2_templ +{ + int x, y; // pixel coordinates relative to the drawing area + unsigned level; // subdivision level: 0 - for pixels, + // 1 - for 1/2 pixels, and so on) + Integer sub_x, sub_y; // subpixel coordinates relative to pixel's boundary + // (always 0 for pixels) + bool operator ==(const Pixel_2_templ& pix) const + { + if(memcmp(this, &pix, sizeof(int)*3)) + return false; + return (sub_x==pix.sub_x&&sub_y==pix.sub_y); + } +}; + +// structure describing the seed point and backward direction, support for +// multiple seed points +template +struct Seed_point_templ +{ + Seed_point_templ() + { } + Seed_point_templ(const Pixel_2_templ& start_, int dir_, + int orient_, int taken_, bool coincide_) : start(start_), + back_dir(dir_), orient(orient_), direction_taken(taken_), + branches_coincide(coincide_) + { } + Pixel_2_templ start; // starting pixel + int back_dir; // backward direction + int orient; // 0 - push_back, 1 - push_front + int direction_taken; + bool branches_coincide; +}; + +template +struct Clip_point_templ // describes a bottom/top clip point +{ + Clip_point_templ() { } + Clip_point_templ(Rational left_, Rational right_, int arcno_=-1) : + left(left_), right(right_), arcno(arcno_) { } + Rational left, right; // isolating interval boundaries + int arcno; // arcno of a segment this point belongs to + AlgebraicReal alpha; // this clip point event line +}; + +template +std::ostream& operator <<(std::ostream& os, const Pixel_2_templ& pix) +{ + os << " (" << pix.x << "/" << pix.sub_x << "; " << pix.y << "/" << + pix.sub_y << ") level = " << pix.level; + return os; +} + +//! this exception is thrown whenever the precision of used number +//! type is not sufficient +class Insufficient_rasterize_precision_exception +{ }; + +struct Coord_2 { + + Coord_2() { + } + Coord_2(int x_, int y_) : x(x_), y(y_) { + } + int x, y; +}; + +typedef std::vector Coord_vec_2; + +/*! \brief defines class \c Curve_renderer_internals + * + * provides an interface to low-level range analysis and polynomials evaluation + * methods along with caching mechanism. \c Coeff_ defines internal polynomial + * coefficients used by the renderer. + * \c CurveKernel_2 specifies underlying 2D curve kernel + */ +template +class Curve_renderer_internals +{ +public: + //! \name typedefs + //!@{ + + //! this instance's first template argument + typedef CurveKernel_2 Curve_kernel_2; + + //! this instance's second template argument + typedef Coeff_ Coeff; + + //! type of 1-curve analysis + typedef typename Curve_kernel_2::Curve_analysis_2 Curve_analysis_2; + + //! type of supporting polynomial + typedef typename Curve_analysis_2::Polynomial_2 Polynomial_2; + + //! underlying univariate polynomial type + typedef typename Polynomial_2::NT Poly_dst_1; + + //! rational number type + typedef typename ::CGAL::Get_arithmetic_kernel< + typename Curve_kernel_2::Coefficient>::Arithmetic_kernel::Rational + Rational; + + typedef Curve_renderer_traits Renderer_traits; + + /// polynomial traits should be used whenever possible + typedef CGAL::Polynomial_traits_d Polynomial_traits_2; + + typedef CGAL::Coercion_traits Coercion; + + //! coercion between rational and polynomial coefficients number type + typedef typename Coercion::Type Rat_coercion_type; + + //! base number type for all internal computations + typedef typename Renderer_traits::Float NT; + + //! instance of a univariate polynomial + typedef CGAL::Polynomial Poly_1; + //! instance of a bivariate polynomial + typedef CGAL::Polynomial Poly_2; + //! container's const iterator (random access) + typedef typename Poly_1::const_iterator const_iterator_1; + //! container's const iterator (random access) + typedef typename Poly_2::const_iterator const_iterator_2; + + //! a univariate rational polynomial + typedef CGAL::Polynomial Rational_poly_1; + //! a bivariate rational polynomial + typedef CGAL::Polynomial Rational_poly_2; + + //! conversion from \c Rational type to used number type + typename Renderer_traits::Rat_to_float rat2float; + //! conversion from the basic NT to \c Rational + typename Renderer_traits::Float_to_rat float2rat; + //! makes the result exact after inexact operation (applicable only for + //! exact number types + typename Renderer_traits::Make_exact make_exact; + + //!@} +private: + //!\name private typedefs + //!@{ + + //! modular traits for bivariate polynomial + typedef CGAL::Modular_traits MT_poly_2; + //! a modular image of a bivariate polynomial + typedef typename MT_poly_2::Modular_NT Modular_poly_2; + //! a modular image of a univariate polynomial + typedef typename Modular_poly_2::NT Modular_poly_1; + //! modular traits for rationals + typedef CGAL::Modular_traits MT_rational; + //! a modular image for rationals + typedef typename MT_rational::Modular_NT Modular_rat; + //! a modular converter for rationals + typename MT_rational::Modular_image image_rat; + + //! container to store derivative coefficients + typedef std::vector Derivative_1; + //! instance of polynomial derivatives type + typedef std::vector > Derivative_2; + //! const_iterator_2 for polynomial derivatives + typedef typename Derivative_2::const_iterator der_iterator_2; + //! const_iterator_1 for polynomial derivatives + typedef typename std::vector::const_iterator der_iterator_1; + + //! hashed container element's type for univariate polynomials + typedef std::pair Poly_entry; + //! hashed key type for function evaluations + typedef std::pair Eval_hash_key; + //! hashed container element's type for function evaluations + typedef std::pair Eval_entry; + + //! hash function to cache univariate polynomials + struct hash_func_poly { + typename Renderer_traits::Hash_function hash; + std::size_t operator()(const NT& key) const { + return hash(key); + } + }; + + //! hash function to cache function evaluations + struct hash_func_eval { + typename Renderer_traits::Hash_function hash; + std::size_t operator()(const Eval_hash_key& key) const { + std::size_t h1 = hash(key.first), h2 = hash(key.second); + return (h1 - h2); + } + }; + + //! hashed map container with LRU capabilities used to store precomputed + //! univariate polynomials at fixed x or y coordinates + typedef boost::multi_index::multi_index_container< + Poly_entry, + boost::multi_index::indexed_by< + boost::multi_index::sequenced<>, + boost::multi_index::hashed_unique< + BOOST_MULTI_INDEX_MEMBER(Poly_entry, NT, first), + hash_func_poly > > > + Poly_cache; + + //! hashed map used to store precomputed polynomial evaluations + typedef boost::multi_index::multi_index_container< + Eval_entry, + boost::multi_index::indexed_by< + boost::multi_index::sequenced<>, + boost::multi_index::hashed_unique< + BOOST_MULTI_INDEX_MEMBER(Eval_entry, Eval_hash_key, first), + hash_func_eval > > > + Eval_cache; + + //!@} +public: + //!\name public interface + //!@{ + + //! default constructor + Curve_renderer_internals() { + } + + //! sets up drawing window and pixel resolution + //! + //! returns \c false if parameters are incorrect + bool setup(const ::CGAL::Bbox_2& box_, int res_w_, int res_h_); + + //! activates a certain cache entry + void select_cache_entry(int cache_id); + + //! precomputes polynomials and derivative coefficients + void precompute(const Polynomial_2& in); + + //! \brief fixes one coordinate of a bivariate polynomial, uses caching + //! if appropriate + void get_precached_poly(int var, const NT& key, int level, Poly_1& poly); + + //! \brief computes the sign of polynomial at (x; y) + //! + //! if computation with the current precision is not reliable, the sign is + //! recomputed with modular or exact rational arithmetic + int evaluate_generic(int var, const NT& c, const NT& key, + const Poly_1& poly); + + //! \brief evaluates a polynomial at (x; y) using modular arithmetic + //! + //! this is a fast way to test for exact zero if the test fails we obtain + //! a correct sign by evaluating rational polynomial + int evaluate_modular(int var, const NT& c, const NT& key); + + //! \brief evaluates a polynomial at (x, y) using exact rational arithmetic + int evaluate_rational(int var, const NT& c, const NT& key); + + //! \brief checks whether interval contains zero + inline bool is_zero(const NT& low, const NT& up) + { + return (low*up < 0); //(low < 0&&up > 0); + } + + //! \brief evalutates a certain polynomial derivative at x + //! + //! \c der_coeffs is a set of derivative coefficients, + //! \c poly - polynomial coefficients + NT evaluate_der(const CGAL::Polynomial& der_coeffs, const Poly_1& poly, + const NT& x) + { + typename Renderer_traits::Extract_eval extract; + const_iterator_1 poly_it = poly.end() - 1; + der_iterator_1 der_it = der_coeffs.end() - 1; + NT y(extract(*poly_it--) * (*der_it)); + while((der_it--) != der_coeffs.begin()) { + y = y * x + extract(*poly_it--) * (*der_it); + } + return y; + } + + //! \brief the same as \c evaluate but arguments are passed by value + //! (needed to substitute variables in bivariate polynomial) + inline static NT binded_eval(Poly_1 poly, NT x) + { + return evaluate(poly, x); + } + + //! \brief evalutates a polynomial at certain x-coordinate + static NT evaluate(const Poly_1& poly, const NT& x, + bool *error_bounds_ = NULL) + { + typename Renderer_traits::Extract_eval extract; + int n = poly.degree()+1, m = (n-1)>>1, odd = n&1; + if(error_bounds_ != NULL) + *error_bounds_ = false; + if(n == 1) + return extract(poly.lcoeff(), error_bounds_); + Coeff cc = static_cast(x); + const_iterator_1 it1 = poly.end()-1, it2 = it1 - (n>>1), + beg = poly.begin()+odd; + Coeff y1 = *it1, y2 = *it2, mul = cc, y; + // unrolled loop for better instruction pairing + while(m-- > odd) { + it1--; it2--; + y1 = y1*cc + (*it1); + y2 = y2*cc + (*it2); + mul = mul*cc; + } + y = y2 + y1*mul; + if(odd) + y = poly[0] + y*cc; + //Gfx_OUT("eval results x = " << x << " res: [" << y.lower() << + // "; " << y.upper() << "\n"); + + return extract(y, error_bounds_); + } + + //! \brief evaluates polynomial at a point (x; y) + //! + //! y is the outermost variable, x is the innermost + static Rat_coercion_type substitute_xy( + const Rational_poly_2& poly, const Rational& x, const Rational& y) { + + typename Rational_poly_2::const_iterator rit = poly.end()-1; + Rat_coercion_type r = rit->evaluate(x), + yc = typename Coercion::Cast()(y); + + while((rit--) != poly.begin()) + r = r * yc + rit->evaluate(x); + return r; + } + + //! First Quadratic Form for univariate case (with recursive + //! derivative technique) + bool get_range_QF_1(int var, const NT& lower, const NT& upper, + const NT& key, const Poly_1& poly, int check = 1); + + //! Modified Affine Arithmetic for univariate case (with recursive + //! derivative technique) + bool get_range_MAA_1(int var, const NT& lower, const NT& upper, + const NT& key, const Poly_1& poly, int check = 1); + + //! empties all cache instances + void clear_caches(); + + //! destructor + ~Curve_renderer_internals() + { + clear_caches(); + } + + //!@} +public: + //!\name public data (can be accessed from hosting curve renderer) + //!@{ + + NT x_min, x_max, y_min, y_max; //! drawing window boundaries + NT pixel_w, pixel_h; //! pixel dimensions w.r.t. resolution + + Rational x_min_r, y_min_r, x_max_r, y_max_r; //! rational versions + Rational pixel_w_r, pixel_h_r; + + int res_w, res_h; //! pixel resolution + + // TODO: make them pointers to const ? + + Poly_2 *coeffs_x, *coeffs_y; //! f(x(y)) / f(y(x)) + Derivative_2 *der_x, *der_y; //! derivative coefficients df/dx (df/dy) + + //! caches of precomputed univariate polynomials in x and y directions + Poly_cache *cached_x, *cached_y; + + Eval_cache *eval_cached; //! cache of polynomial evaluations + + Modular_poly_2 *modular_x, *modular_y; //! modular images of a polynomial + + Rational_poly_2 *rational_x, *rational_y; //! poly with rational coeffs + + Rational_poly_2 *rational_fx, *rational_fy; //! partial derivatives + + //! 0 - the 1st (2nd) derivative does not have zero over an + //! interval; 1 - does have 0; -1 - not computed + int first_der, second_der; + + bool zero_bounds; //! indicates that the result of the last range + //! evaluation has at least one zero boundary + + static bool show_dump; //! for debugging + + //!@} +private: + //!\name private members for caching support + //!@{ + + //! data structures grouped into arrays to facilitate caching + Poly_2 coeffs_x_[CGAL_N_CACHES], coeffs_y_[CGAL_N_CACHES]; + Derivative_2 der_x_[CGAL_N_CACHES], der_y_[CGAL_N_CACHES]; + Poly_cache cached_x_[CGAL_N_CACHES], cached_y_[CGAL_N_CACHES]; + Eval_cache eval_cached_[CGAL_N_CACHES]; + Modular_poly_2 modular_x_[CGAL_N_CACHES], modular_y_[CGAL_N_CACHES]; + Rational_poly_2 rational_x_[CGAL_N_CACHES], rational_y_[CGAL_N_CACHES]; + Rational_poly_2 rational_fx_[CGAL_N_CACHES], rational_fy_[CGAL_N_CACHES]; + + //!@} +}; + +template +bool Curve_renderer_internals::show_dump = false; + +//! sets up drawing window and pixel resolution +//! +//! returns \c false if parameters are incorrect +template +bool Curve_renderer_internals::setup( + const ::CGAL::Bbox_2& box_, int res_w_, int res_h_) +{ + x_min = static_cast(box_.xmin()); + y_min = static_cast(box_.ymin()); + x_max = static_cast(box_.xmax()); + y_max = static_cast(box_.ymax()); + res_w = res_w_; + res_h = res_h_; + + if(x_min >= x_max||y_min >= y_max||res_w < 5||res_h < 5||res_w > 1024|| + res_h > 1024) { + Gfx_OUT("Incorrect setup parameters" << std::endl); + return false; + } + + x_min_r = Rational(box_.xmin()); + y_min_r = Rational(box_.ymin()); + x_max_r = Rational(box_.xmax()); + y_max_r = Rational(box_.ymax()); + + pixel_w_r = (x_max_r - x_min_r) / res_w; + pixel_h_r = (y_max_r - y_min_r) / res_h; +// NiX::simplify(pixel_w_r); +// NiX::simplify(pixel_h_r); + + pixel_w = rat2float(pixel_w_r); + pixel_h = rat2float(pixel_h_r); + make_exact(pixel_w); + make_exact(pixel_h); + + show_dump = false; + + return true; +} + +//! \brief evaluates a univariate polynomial over an interval using +//! First Quadratic Affine Form +template +bool Curve_renderer_internals::get_range_QF_1( + int var, const NT& l_, const NT& r_, const NT& key, const Poly_1& poly, + int check) +{ + Derivative_2 *der = (var == CGAL_X_RANGE ? der_x : der_y); + NT l(l_), r(r_), l1, r1, low, up; + NT v1, v2; + int eval1, eval2; + der_iterator_2 der_it_2 = der->end()-1; + der_iterator_1 der_it, der_begin; + const_iterator_1 cache_it, begin; + + first_der = false; + if(poly.degree()==0) { + zero_bounds = false; + return (poly.lcoeff()==NT(0.0)); + } + if(l_ == r_) { + zero_bounds = false; + return false; + } + if(l > r) { + l = r_; + r = l_; + } + + eval1 = evaluate_generic(var, l, key, poly); + eval2 = evaluate_generic(var, r, key, poly); + bool sign_change = (eval1*eval2 < 0); + + zero_bounds = ((eval1&eval2) == 0); + if((sign_change||zero_bounds)&&check==1) + return true; + + if(var == CGAL_X_RANGE) { + l1 = x_min + l*pixel_w; + r1 = x_min + r*pixel_w; + } else { + l1 = y_min + l*pixel_h; + r1 = y_min + r*pixel_h; + } + + typename Renderer_traits::Extract_eval extract; + unsigned index = CGAL_RECURSIVE_DER_MAX_DEGREE; + if(index >= der->size()) { + low = up = extract(poly.lcoeff()) * (*der_it_2).lcoeff(); + } else { + der_it_2 = der->begin() + index; + low = 1; + up = -1; + } + + NT x0 = (l1 + r1)/2, x1 = (r1 - l1)/2; + make_exact(x0); + make_exact(x1); + NT x0_abs = CGAL_ABS(x0), x1_abs = CGAL_ABS(x1); + + while((der_it_2--)!=der->begin()) { + + // iterate through derivative coefficients + der_it = (*der_it_2).end()-1; + der_begin = (*der_it_2).begin(); + cache_it = poly.end()-1; // iterate through precomputed y-values + + // if a derivative does not straddle zero we can obtain exact bounds + // by evaluating a polynomial at end-points + if(low * up >= 0) { + v1 = v2 = extract(*cache_it--)* (*der_it); + // calculate the ith derivative at xa and xb + while((der_it--) != der_begin) { + NT cc1 = extract(*cache_it--)* (*der_it); + v1 = v1 * l1 + cc1; + v2 = v2 * r1 + cc1; + } + low = v1; + up = v2; + } else { // use Quadratic Form to compute the bounds + + NT y0 = extract(*cache_it) * (*der_it), y1(0), z1(0), e1(0); + cache_it--; + while((der_it--)!=der_begin) { + e1 = x0_abs*e1 + x1_abs*e1 + CGAL_ABS(x1*z1); + z1 = x0*z1 + x1*y1; + y1 = y1*x0 + x1*y0; + y0 = x0*y0 + extract(*cache_it)*(*der_it); + cache_it--; + } + NT spread = CGAL_ABS(y1) + e1; + low = spread; + up = spread; + if(z1 > 0) + up = up + z1; + else + low = low - z1; + low = y0 - low; + up = y0 + up; + } + + if(der_it_2 == der->begin() && check==3) { + first_der = //(eval1*eval2 < 0);// + is_zero(low, up); + if(sign_change||zero_bounds) + return true; + } + } + + if(low * up < 0) { + cache_it = poly.end()-1; + begin = poly.begin(); + NT y0 = extract(*cache_it), y1(0), z1(0), e1(0); + while((cache_it--) != begin) { + e1 = x0_abs*e1 + x1_abs*e1 + CGAL_ABS(x1*z1); + z1 = x0*z1 + x1*y1; + y1 = y1*x0 + x1*y0; + y0 = x0*y0 + extract(*cache_it); + } + NT spread = CGAL_ABS(y1) + e1; + low = spread; + up = spread; + if(z1 > 0) + up = up + z1; + else + low = low - z1; + low = y0 - low; + up = y0 + up; + eval1 = CGAL_SGN(low); + eval2 = CGAL_SGN(up); + } + + zero_bounds = ((eval1 & eval2) == 0); + return (eval1*eval2 < 0); +} + +//! \brief evaluates a univariate polynomial over an interval using +//! Modified Affine Arithmetic +template +bool Curve_renderer_internals::get_range_MAA_1( + int var, const NT& l_, const NT& r_, const NT& key, const Poly_1& poly, + int check) +{ + Derivative_2 *der = (var == CGAL_X_RANGE) ? der_x : der_y; + // stores precomputed polynomial derivatives and binominal coeffs + Derivative_1 der_cache + //(der->size()+1, NT(0)) + , binom;//(der->size()+1, NT(0)); + NT l(l_), r(r_), l1, r1, low = NT(1), up = NT(-1); + NT v1, v2, v; + int eval1, eval2; + + first_der = false; + if(poly.degree()==0) { + zero_bounds = false; + return (poly.lcoeff()==NT(0.0)); + } + if(l_ == r_) { + first_der = false; + return false; + } + if(l > r) { + l = r_; + r = l_; + } + eval1 = evaluate_generic(var, l, key, poly); + eval2 = evaluate_generic(var, r, key, poly); + + bool sign_change = (eval1*eval2 < 0); + zero_bounds = ((eval1&eval2) == 0); + + if((sign_change || zero_bounds) && check == 1) + return true; + + if(var == CGAL_X_RANGE) { + l1 = x_min + l*pixel_w; + r1 = x_min + r*pixel_w; + } else { + l1 = y_min + l*pixel_h; + r1 = y_min + r*pixel_h; + } + NT x0 = (l1 + r1)/2, x1 = (r1 - l1)/2; + make_exact(x0); + make_exact(x1); + + int d = 0; + v1 = evaluate(poly, x0); + //der_cache[d] = v1; + der_cache.push_back(v1); + v = x1, d++; + der_iterator_2 der_it_2; + for(der_it_2 = der->begin(); der_it_2 != der->end(); der_it_2++) { + v1 = evaluate_der((*der_it_2), poly, x0); + der_cache.push_back(v1); + //der_cache[d] = v1; // replace by push_backs ? + //binom[d] = v; + binom.push_back(v); + d++; + //v *= x1/d; + //make_exact(v); + v *= x1; // ONLY when derivative coefficients are normalized + } + + typename Renderer_traits::Extract_eval extract; + unsigned index = CGAL_RECURSIVE_DER_MAX_DEGREE; + if(index >= der->size()) { + low = up = extract(poly.lcoeff()) * (*der_it_2).lcoeff(); + } else { + der_it_2 = der->begin()+index; + low = 1; + up = -1; + } + // assume we have an array of derivatives: + // der_cache: {f^(0); f^(1); f^(2); ...} + // and binominal coefficients: [h; h^2/2; h^3/6; ... h^d/d!] + der_iterator_1 eval_it = der_cache.end()-1, local_it, binom_it, + eval_end = der_cache.end(); + d = poly.degree(); + der_it_2 = der->end()-1; + while(1) {//der_it!=end) { + if(low * up < 0) { + v2 = *eval_it; + local_it = eval_it + 1; + binom_it = binom.begin(); + low = v2, up = v2; + while((local_it) != eval_end) {// calculate derivative bounds + v1 = (*local_it++) * (*binom_it++); + if(v1 >= 0) { // derivative index is odd + up += v1; + low -= v1; + } else { + up -= v1; + low += v1; + } + if(local_it == eval_end) + break; + // derivative index is even + v1 = (*local_it++) * (*binom_it++); + (v1 > 0) ? up += v1 : low += v1; + } + eval1 = CGAL_SGN(low); + eval2 = CGAL_SGN(up); + + } else if(d > 0) { + low = evaluate_der((*der_it_2), poly, l1); + up = evaluate_der((*der_it_2), poly, r1); + eval1 = CGAL_SGN(low); + eval2 = CGAL_SGN(up); + + } else if(d == 0) + ;//Gfx_DETAILED_OUT("MAA bounds: sign change\n"); + + if(d == 1&&check == 3) { + first_der = (eval1*eval2 < 0); + if(sign_change||zero_bounds) + return true; + } else if(d == 0) { + + zero_bounds = ((eval1&eval2) == 0); + return (eval1*eval2 < 0); + } + d--; der_it_2--; + if((eval_it--) == der_cache.begin()) + break; + } + return true; +} + +//! \brief fixes one coordinate of a bivariate polynomial, uses caching +//! if appropriate +template +void Curve_renderer_internals::get_precached_poly( + int var, const NT& key, int level, Poly_1& poly) +{ + NT key1; + Poly_cache *cached = cached_x; + Poly_2 *coeffs = coeffs_x; + typename boost::multi_index::nth_index_iterator::type it; + bool not_cached = (level >= CGAL_MAX_POLY_CACHE_LEVEL), not_found = false; + + if(var == CGAL_Y_RANGE) { + cached = cached_y; + coeffs = coeffs_y; + key1 = x_min + key*pixel_w; + } else + key1 = y_min + key*pixel_h; + + if(!not_cached) { + typename boost::multi_index::nth_index::type& + idx = cached->get<1>(); + it = idx.find(key1); //*4 + not_found = (it == idx.end()); + } + + if(not_cached||not_found) { + poly = Poly_1(::boost::make_transform_iterator(coeffs->begin(), + std::bind2nd(std::ptr_fun(binded_eval), key1)), + ::boost::make_transform_iterator(coeffs->end(), + std::bind2nd(std::ptr_fun(binded_eval), key1))); + if(not_cached) + return; + // all available space consumed: drop the least recently used entry + if(cached->size() >= CGAL_POLY_CACHE_SIZE) + cached->pop_back(); + cached->push_front(Poly_entry(key1,poly)); //*4 + return; + } + cached->relocate(cached->begin(), cached->project<0>(it)); + poly = (*it).second; +} + +//! \brief computes the sign of univariate polynomial at (x; y). +//! +//! if computation with the current precision is not reliable, the sign is +//! recomputed with modular or exact rational arithmetic +template +int Curve_renderer_internals::evaluate_generic( + int var, const NT& c, const NT& key, const Poly_1& poly) +{ + NT x = key, y = c, c1; + Eval_hash_key hash_key; + typename boost::multi_index::nth_index_iterator::type it; + bool not_cached = true, + //(max_level > CGAL_MAX_EVAL_CACHE_LEVEL), + not_found = false; + + // TODO define max_level somehow + + if(var == CGAL_X_RANGE) { + x = c; + y = key; + c1 = x_min + c*pixel_w; + } else + c1 = y_min + c*pixel_h; + + if(show_dump) { + //Gfx_OUT("evaluate " << poly << " at: " << c1 << "\n"); + } + + if(!not_cached) { + hash_key.first = x_min + x*pixel_w; + hash_key.second = y_min + y*pixel_h; + typename boost::multi_index::nth_index::type& + idx = eval_cached->get<1>(); + it = idx.find(hash_key); + not_found = (it == idx.end()); + } + + if(not_cached||not_found) { + + bool error_bounds_; + int sign; + NT res = evaluate(poly, c1, &error_bounds_); + + if(error_bounds_) { + //Gfx_OUT("error_bounds_\n"); + sign = evaluate_modular(var, c, key); + if(sign != 0) + sign = evaluate_rational(var, c, key); + } else + sign = CGAL_SGN(res); + if(not_cached) + return sign; + // drop the least recently used entry if cache is full + if(eval_cached->size() >= CGAL_EVAL_CACHE_SIZE) + eval_cached->pop_back(); + eval_cached->push_front(Eval_entry(hash_key, sign)); + return sign; + } + + eval_cached->relocate(eval_cached->begin(), + eval_cached->project<0>(it)); + return (*it).second; +} + +//! \brief evaluates a polynomial at (x, y) using modular arithmetic +//! +//! this is a fast way to test for exact zero if the test fails we obtain +//! a correct sign by evaluating with rational polynomial, returns -1, 0 or 1 +//! depending on the sign of the evaluation +template +int Curve_renderer_internals::evaluate_modular( + int var, const NT& c, const NT& key) +{ +#if !AcX_SQRT_EXTENSION + Modular_poly_1 poly; + Rational c_r = float2rat(c), key_r = float2rat(key); + c_r = (var == CGAL_X_RANGE) ? (x_min_r + c_r*pixel_w_r) : + (y_min_r + c_r*pixel_h_r); + + Modular_poly_2 *mod = modular_x; + if(var == CGAL_Y_RANGE) { + mod = modular_y; + key_r = x_min_r + key_r*pixel_w_r; + } else + key_r = y_min_r + key_r*pixel_h_r; + + //!@todo: get rid of NiX::substitute_x !!! + poly = NiX::substitute_x(*mod, image_rat(key_r)); + Modular_rat res = poly.evaluate(image_rat(c_r)); + return CGAL_SGN(res.x()); +#else //!@todo: decide by Algebraic_structure_traits::Field_with_sqrt + return -1; // modular arithmetic is disabled with sqrt extension +#endif +} + +//! \brief evaluates a polynomial at (x, y) using exact rational arithmetic +template +int Curve_renderer_internals::evaluate_rational( + int var, const NT& c, const NT& key) +{ + //Rational_poly_1 poly; + Rational c_r = float2rat(c), key_r = float2rat(key); + c_r = (var == CGAL_X_RANGE) ? (x_min_r + c_r*pixel_w_r) : + (y_min_r + c_r*pixel_h_r); + + Rational_poly_2 *rat = rational_x; + if(var == CGAL_Y_RANGE) { + rat = rational_y; + key_r = x_min_r + key_r*pixel_w_r; + } else + key_r = y_min_r + key_r*pixel_h_r; + + Rat_coercion_type res = substitute_xy(*rat, key_r, c_r); + return CGAL_SGN(res); +} + +//! precomputes polynomials and derivative coefficients +template +void Curve_renderer_internals::precompute( + const Polynomial_2& in) +{ + typedef typename Polynomial_traits_2::Innermost_coefficient Coeff_src; + + Max_coeff max_coeff; + Coeff_src max_c = max_coeff(in); + /////// magic symbol //////////// + // somehow relates to double precision fix + std::cerr << ' '; + + typedef Reduce_by Reduce_op; + Transform transform; + Reduce_op op(max_c); + + typedef CGAL::Polynomial_traits_d RP_traits; + *rational_y = transform(in, op); + *rational_x = typename RP_traits::Swap()(*rational_y, 0, 1); + + // rational fx and fy must have y outermost var and x innermost + *rational_fx = typename RP_traits::Derivative()(*rational_y, 0); + *rational_fy = typename RP_traits::Derivative()(*rational_y, 1); + + // modular polynomials are not used with Field_with_sqrt +//#if !AcX_SQRT_EXTENSION + typename MT_poly_2::Modular_image image_poly_2; + *modular_y = image_poly_2(*rational_y); + *modular_x = typename CGAL::Polynomial_traits_d:: + Swap()(*modular_y, 0, 1); +//#endif + + typename Renderer_traits::Convert_poly convert_poly; + //////////////////////////////////////////////////////// + /////// ATTENTION: need to call makeExact for bigfloats after conversion + //////////////////////////////////////////////////////// + *coeffs_y = convert_poly(*rational_y); + *coeffs_x = typename CGAL::Polynomial_traits_d:: + Swap()(*coeffs_y, 0, 1); + + int degree_x = coeffs_x->degree(), + degree_y = coeffs_y->degree(); + cached_x->clear(); + cached_y->clear(); + eval_cached->clear(); + der_x->clear(); + der_y->clear(); + int i, j, maxdeg = (degree_x > degree_y ? degree_x : degree_y); + std::vector X(maxdeg, NT(0)); + der_x->reserve(degree_x); + der_y->reserve(degree_y); + + NT det(1.0); + for(i = 0; i < degree_x; i++) { + if(i != 0) + det = X[0]; + for(j = 1; j <= degree_x - i; j++) { + // divide by the lowest coefficient ? + X[j - 1] = (i == 0 ? NT(j) : X[j] * NT(j) / det); + make_exact(X[j-1]); + } + der_x->push_back(CGAL::Polynomial(X.begin(), + (X.begin() + degree_x - i))); + } + + det = NT(1.0); + for(i = 0; i < degree_y; i++) { + if(i != 0) + det = X[0]; + for(j = 1; j <= degree_y - i; j++) { + // divide by the lowest coefficient ? + X[j - 1] = (i == 0 ? NT(j) : X[j] * NT(j) / det); + make_exact(X[j-1]); + } + der_y->push_back(CGAL::Polynomial(X.begin(), + (X.begin() + degree_y - i))); + } +} + +//! \brief activates the cache entry \c cache_id +template +void Curve_renderer_internals::select_cache_entry( + int cache_id) +{ + coeffs_x = coeffs_x_ + cache_id; + coeffs_y = coeffs_y_ + cache_id; + modular_x = modular_x_ + cache_id; + modular_y = modular_y_ + cache_id; + rational_x = rational_x_ + cache_id; + rational_y = rational_y_ + cache_id; + rational_fx = rational_fx_ + cache_id; + rational_fy = rational_fy_ + cache_id; + der_x = der_x_ + cache_id; + der_y = der_y_ + cache_id; + cached_x = cached_x_ + cache_id; + cached_y = cached_y_ + cache_id; + eval_cached = eval_cached_ + cache_id; +} + +//! \brief empties all cache instances +template +void Curve_renderer_internals::clear_caches() +{ + for(unsigned i = 0; i < CGAL_N_CACHES; i++) { + cached_x_[i].clear(); + cached_y_[i].clear(); + eval_cached_[i].clear(); + der_x_[i].clear(); + der_y_[i].clear(); + } +} + +} // namespace CGALi + +CGAL_END_NAMESPACE + +#endif // CGAL_CKVA_CURVE_RENDERER_INTERNALS_H diff --git a/Curved_kernel_via_analysis_2/include/CGAL/Curved_kernel_via_analysis_2/gfx/Curve_renderer_traits.h b/Curved_kernel_via_analysis_2/include/CGAL/Curved_kernel_via_analysis_2/gfx/Curve_renderer_traits.h new file mode 100644 index 00000000000..7d7fab7d244 --- /dev/null +++ b/Curved_kernel_via_analysis_2/include/CGAL/Curved_kernel_via_analysis_2/gfx/Curve_renderer_traits.h @@ -0,0 +1,575 @@ +// TODO: Add licence +// +// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE +// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. +// +// $URL:$ +// $Id: $ +// +// +// Author(s) : Pavel Emeliyanenko +// +// ============================================================================ + +#ifndef CGAL_CKVA_CURVE_RENDERER_TRAITS_H +#define CGAL_CKVA_CURVE_RENDERER_TRAITS_H + +#include +#include + +/*! \file CGAL/Curved_kernel_via_analysis_2/gfx/Curve_renderer_traits.h + * \brief + * defines class Curve_renderer_traits. + * + * provides specialisations of Curve_renderer_traits for different number + * types compatible with the curve renderer +*/ + +#ifndef CGAL_MAX_COEFF_TRANSFORM +#define CGAL_MAX_COEFF_TRANSFORM + +#ifndef CGAL_ABS +#define CGAL_ABS(x) ((x) < 0 ? -(x): (x)) +#endif + +#ifndef CGAL_SGN +#define CGAL_SGN(x) ((x) > 0 ? 1 : ((x) < 0 ? -1 : 0)) +#endif + +CGAL_BEGIN_NAMESPACE + +// transformation routine +namespace CGALi { + +template +struct Max_coeff +{ + template + NT operator()(const CGAL::Polynomial& p) const + { + typename CGAL::Polynomial::const_iterator it = p.begin(); + Max_coeff max_coeff; + NT max(max_coeff(*it)); + while(++it != p.end()) { + NT tmp(max_coeff(*it)); + if(max < CGAL_ABS(tmp)) + max = CGAL_ABS(tmp); + } + return max; + } + NT operator()(const NT& x) const + { return CGAL_ABS(x); } +}; + +/*!\brief + * divides an input value by a contant + * + * provided that there is a coercion between \c Input and \c Result types + */ +template +struct Reduce_by { + + typedef Input argument_type; + typedef Result result_type; + + Reduce_by(const Input& denom_) : + denom(cast(denom_)) { + } + + Result operator()(const Input& x) { + return (cast(x)/denom); + } + + typename CGAL::Coercion_traits::Cast cast; + Result denom; +}; + +/*!\brief + * transforms bivaritate polynomial of type \c InputPoly_2 to + * \c OutputPoly_2 by recursively applying operation \c Op to all of its + * coefficients + * + * Op: InputPoly_2::Innermost_coefficient -> + * OutputPoly_2::Innermost_coefficient + */ +template +struct Transform { + + typedef InputPoly_2 first_argument_type; + typedef Op second_argument_type; + typedef OutputPoly_2 result_type; + + template + OutputPoly_2 operator()(const CGAL::Polynomial& p, Op op = Op()) const { + + Transform tr; + return OutputPoly_2( + ::boost::make_transform_iterator(p.begin(), std::bind2nd(tr, op)), + ::boost::make_transform_iterator(p.end(), std::bind2nd(tr, op))); + } + + OutputPoly_2 operator()( + const typename Innermost_coefficient::Type& x, Op op) + const { + + return static_cast(op(x)); + } +}; + +#endif // CGAL_MAX_COEFF_TRANSFORM + +/*!\brief + * class template \c Curve_renderer_traits + * + * this traits class prodives various number type conversions for the + * curve renderer + */ + +template +struct Curve_renderer_traits_base +{ + //! type of innermost polynomial coefficients + typedef Coeff_ Coeff; + + //! an integer number type + typedef Integer_ Integer; + + typedef Rational_ Rational; + + //! conversion from rational to floating-point + template + struct Rat_to_float { + typedef Float result_type; + + template + Float operator()(const Sqrt_extension& x) const { + typename CGAL::Coercion_traits, Float>::Cast + cast; + return cast(x); + } + + Float operator()(const Rational& x) const { + return static_cast(x); + } + }; + + //! provided for convenience when there exists an implicit coercion + //! between number types + template + struct Implicit_coercion { + typedef To result_type; + + template + To operator()(const From& x) const { + return static_cast(x); + } + }; + + struct Float_to_int { + typedef int result_type; + + template + int operator()(const Float& x) const + { return static_cast(std::floor(CGAL::to_double(x))); } + //return static_cast(std::floor(x)); } + }; + + /*!\brief + * conversion from bivariate polynomial over extended number type to + * polynomial with coefficients of type \c Coeff + * + * valid instantiations of \c Extended number type are \c Rational , + * \c FieldWithSqrt , etc. Provided that there is a type coercion between + * \c Extended and \c Coeff + */ + struct Convert_poly { + typedef CGAL::Polynomial > result_type; + + template + inline result_type operator()(const + CGAL::Polynomial >& poly) const { + + //!@todo: use Rat_to_float functor instead of coercion traits ? + //! need some sort of polymorphism.. + typedef typename CGAL::Coercion_traits::Cast + Cast; + Transform >, Cast> + transform; + return transform(poly); + } + }; + + //! converts polynomial coefficients to floating-point representation + //! \c error_bounds is set if the result is not reliable + struct Extract_eval { + typedef Coeff argument_type; + typedef Coeff result_type; + + Coeff operator()(const Coeff& x, + bool *error_bounds = NULL) const { + if(error_bounds != NULL) + *error_bounds = false; + return x; + } + }; + + //! computes a 32-bit hash value from floating-point argument + struct Hash_function { + typedef std::size_t result_type; + struct long_long { + long low, high; + }; + + template + std::size_t operator()(const Float& key) const { + const long_long *hk = reinterpret_cast(&key); + return (hk->low ^ hk->high); + } + }; + + //! makes result exact after inexact operation such as div, sqrt, etc. + //! (required for multi-precision arithmetic) + struct Make_exact { + typedef void result_type; + + template + void operator()(const Float& x) const + { } + }; + + //! compares a given quantity with the precision limit a floating-point + //! number type, returns \c true if this limit is exceeded + struct Precision_limit { + typedef bool result_type; + + template + bool operator()(const Float& x) const + { return false;/*(CGAL_ABS(x) <= 1e-16)*/; } + }; + + //! maximum subdivision level for the curve renderer by exceeding which + //! an exception of type \c Insufficient_rasterize_precision_exception + //! will be thrown, this is also limited by \c Integer number type, since + //! the integer must be able to store numbers up to 2^MAX_SUBDIVISION_LEVEL + static const unsigned MAX_SUBDIVISION_LEVEL = 12; +}; + +template +struct Curve_renderer_traits +{ }; + +#ifdef CGAL_USE_CORE +//! Specialization for \c CGAL::Interval_nt +template <> +struct Curve_renderer_traits, CORE::BigRat > : + public Curve_renderer_traits_base, int, + CORE::BigRat> { + + typedef double Float; + + struct Rat_to_float { + typedef Float result_type; + + template + Float operator()(const Extended& x) const { + return CGAL::to_double(x); + } + }; + + typedef Implicit_coercion Float_to_rat; + + struct Rat_to_integer { + typedef Rational argument_type; + typedef Integer result_type; + + Integer operator()(const Rational& x) const { + return static_cast(std::floor(CGAL::to_double(x))); + } + }; + + struct Extract_eval { + typedef Coeff argument_type; + typedef Float result_type; + + Float operator()(const Coeff& x, + bool *error_bounds = NULL) const { + bool err_bnd; + Float l = x.inf(), u = x.sup(), mid = (l+u)/2; + //err_bnd = ((l <= 0 && u >= 0)); + err_bnd = (CGAL_ABS(l) < 1E-15 || CGAL_ABS(u) < 1E-15) || + ((l <= 0 && u >= 0)); + if(error_bounds != NULL) + *error_bounds = err_bnd; + if(err_bnd)// && CGAL_ABS(mid) < 1E-15) + return 0; + return mid; + } + }; + + //! compares a given quantity with the precision limit a floating-point + //! number type, returns \c true if this limit is exceeded + struct Precision_limit { + typedef bool result_type; + + template + bool operator()(const Float& x) const + { return (CGAL_ABS(x) <= 1e-16); } + }; +}; + +//! Specialization for \c CORE::BigFloat +template <> +struct Curve_renderer_traits + : public Curve_renderer_traits_base { + + typedef CORE::BigFloat Float; + + struct Rat_to_integer { + typedef Rational argument_type; + typedef Integer result_type; + + Integer operator()(const Rational& x) const { + return x.BigIntValue(); + } + }; + + typedef Rat_to_float Rat_to_float; + + struct Float_to_rat { + typedef Float argument_type; + typedef Rational result_type; + + Rational operator()(const Float& x) const + { return x.BigRatValue(); } + }; + + struct Hash_function { + typedef Float argument_type; + typedef std::size_t result_type; + + inline result_type operator()(const Float& key) const { + const BigFloatRep& rep = key.getRep(); + std::size_t ret = reinterpret_cast(&rep); + return ret; + } + }; + + struct Make_exact { + typedef Float argument_type; + typedef void result_type; + + inline void operator()(Float& x) const + { x.makeExact(); } + }; +}; + +//! Specialization for \c CORE::BigRat +template <> +struct Curve_renderer_traits : + public Curve_renderer_traits_base { + + typedef CORE::BigRat Float; + + typedef Rat_to_float Rat_to_float; + + typedef Implicit_coercion Float_to_rat; + + struct Rat_to_integer { + typedef Rational argument_type; + typedef Integer result_type; + + Integer operator()(const Rational& x) const { + return x.BigIntValue(); + } + }; + + struct Hash_function { + typedef Float argument_type; + typedef std::size_t result_type; + + inline result_type operator()(const Float& key) const { + const BigRatRep& rep = key.getRep(); + std::size_t ret = reinterpret_cast(&rep); + return ret; + } + }; +}; +#endif // CGAL_USE_CORE + +#ifdef CGAL_USE_LEDA +template <> +struct Curve_renderer_traits, leda::rational > : + public Curve_renderer_traits_base, int, + leda::rational> { + + typedef double Float; + + struct Rat_to_float { + typedef Float result_type; + + template + Float operator()(const Extended& x) const { + return CGAL::to_double(x); + } + }; + + typedef Implicit_coercion Float_to_rat; + + struct Rat_to_integer { + typedef Rational argument_type; + typedef Integer result_type; + + Integer operator()(const Rational& x) const { + return static_cast(std::floor(CGAL::to_double(x))); + //return static_cast(x.to_double()); + } + }; + + struct Extract_eval { + typedef Coeff argument_type; + typedef Float result_type; + + Float operator()(const Coeff& x, + bool *error_bounds = NULL) const { + bool err_bnd; + Float l = x.inf(), u = x.sup(), mid = (l+u)/2; + //err_bnd = ((l <= 0 && u >= 0)); + err_bnd = (CGAL_ABS(l) < 1E-15 || CGAL_ABS(u) < 1E-15) || + ((l <= 0 && u >= 0)); + if(error_bounds != NULL) + *error_bounds = err_bnd; + if(err_bnd)// && CGAL_ABS(mid) < 1E-15) + return 0; + return mid; + } + }; + + //! compares a given quantity with the precision limit a floating-point + //! number type, returns \c true if this limit is exceeded + struct Precision_limit { + typedef bool result_type; + + template + bool operator()(const Float& x) const + { return (CGAL_ABS(x) <= 1e-16); } + }; +}; + +//! Specialization for \c leda::bigfloat +template <> +struct Curve_renderer_traits + : public Curve_renderer_traits_base { + + typedef leda::bigfloat Float; + + struct Rat_to_integer { + typedef Rational argument_type; + typedef Integer result_type; + + Integer operator()(const Rational& x) const { + return leda::floor(x); + } + }; + + struct Rat_to_float { + typedef Float result_type; + + template + Float operator()(const Sqrt_extension& x) const { + typename CGAL::Coercion_traits, Float>::Cast + cast; + return cast(x); + } + // no implicit coercion between leda rational and floats, therefore + // decompose rational to compute the result + Float operator()(const Rational& x) const { + return (static_cast(x.numerator()) / + static_cast(x.denominator())); + } + }; + + struct Float_to_rat { + typedef Float argument_type; + typedef Rational result_type; + + Rational operator()(const Float& x) const + { return x.to_rational(); } + }; + + struct Convert_poly { + typedef CGAL::Polynomial > result_type; + + template + inline result_type operator()(const + CGAL::Polynomial >& poly) const { + + Transform >, + Rat_to_float> transform; + return transform(poly); + } + }; + + struct Hash_function { + typedef Float argument_type; + typedef std::size_t result_type; + + inline result_type operator()(const Float& key) const { + return static_cast( + key.get_significant_length()); + } + }; +}; + +//! Specialization for \c leda::rational +template <> +struct Curve_renderer_traits : + public Curve_renderer_traits_base { + + typedef leda::rational Float; + + typedef Rat_to_float Rat_to_float; + + typedef Implicit_coercion Float_to_rat; + + struct Rat_to_integer { + typedef Rational argument_type; + typedef Integer result_type; + + Integer operator()(const Rational& x) const { + return leda::floor(x); + } + }; + + struct Hash_function { + typedef Float argument_type; + typedef std::size_t result_type; + + inline result_type operator()(const Float& key) const { + std::size_t ret = reinterpret_cast( + key.numerator().word_vector()); + return ret; + } + }; + + struct Make_exact { + typedef Float argument_type; + typedef void result_type; + + inline void operator()(Float& x) const + { x.normalize(); } + }; +}; + +#endif // CGAL_USE_LEDA + +} // namespace CGALi + +CGAL_END_NAMESPACE + +#endif // CGAL_CKVA_CURVE_RENDERER_TRAITS_H diff --git a/Curved_kernel_via_analysis_2/include/CGAL/Curved_kernel_via_analysis_2/gfx/Subdivision_1.h b/Curved_kernel_via_analysis_2/include/CGAL/Curved_kernel_via_analysis_2/gfx/Subdivision_1.h new file mode 100644 index 00000000000..c778bfa1174 --- /dev/null +++ b/Curved_kernel_via_analysis_2/include/CGAL/Curved_kernel_via_analysis_2/gfx/Subdivision_1.h @@ -0,0 +1,432 @@ +// TODO: Add licence +// +// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE +// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. +// +// $URL:$ +// $Id: $ +// +// +// Author(s) : Pavel Emeliyanenko +// +// ============================================================================ + +/*!\file CGAL/Curved_kernel_via_analysis_2/gfx//Subdivision_1.h + * \brief definition of Subdivision_1<> + * 1D space subdivision for rasterization of planar curves + */ + +#ifndef CGAL_CKVA_SUBDIVISION_1_H +#define CGAL_CKVA_SUBDIVISION_1_H 1 + +#include +#include +#include +//#include +#include +#include + +#include +#include + +#include +#include + +#include + +using boost::multi_index::multi_index_container; +using boost::multi_index::get; +using boost::multi_index::project; + +CGAL_BEGIN_NAMESPACE + +namespace CGALi { + +#ifndef SoX_CURVE_RENDERER_DEFS +#define SoX_CURVE_RENDERER_DEFS + +#define SoX_WINDOW_ENLARGE 15 // # of pixels by which the drawing window + // is enlarged in y-direction: this is used + // to skip closely located clip-points + +#define SoX_REFINE_X 1000 // refine factor for X-intervals + // (in pixel size) +#define SoX_REFINE_Y 100000 // refine factor for Y-intervals + // (in pixel size) + +#endif // SoX_CURVE_RENDERER_DEFS + +#define SoX_REFINE_ISOLATED_POINTS 2 // refine factor for clip points + +template +struct Range_templ +{ + Range_templ() { } + Range_templ(const NT& l_, const NT& r_) : left(l_), right(r_) + { + sign_change = true; + } + NT left, right; + bool sign_change; +}; + +} // namespace CGALi + +/*!\brief + * The class template \c Subdivision_1 and its associate functions. + * + * The class implements a space method to plot algebraic curves, we use Affine + * Arithmetic with recursive derivative information + */ +template +class Subdivision_1 +{ +public: + //! \name public typedefs + //!@{ + + //! this instance's first template argument + typedef Coeff_ Coeff; + //! this instance's second template argument + typedef Algebraic_curve_2_ Algebraic_curve_2; + + //!@} +private: + //! \name private typedefs + //!@{ + //! rational number type + typedef typename Algebraic_curve_2::Rational Rational; + + //! special number type traits dependent on polynomial coefficient + typedef SoX::Curve_renderer_traits Renderer_traits; + + //! specialized integer number type + typedef typename Renderer_traits::Integer Integer; + //! supporting bivariate polynomial type + typedef typename Algebraic_curve_2::Poly_d Poly_dst_2; + //! supporting univariate polynomial type + typedef typename Algebraic_curve_2::Poly_d::NT Poly_dst_1; + + //! a univariate rational polynomial + typedef CGAL::Polynomial Rational_poly_1; + //! a bivariate rational polynomial + typedef CGAL::Polynomial Rational_poly_2; + + //! basic number type used in all computations + typedef typename Renderer_traits::Float NT; + //! instance of a univariate polynomial + typedef CGAL::Polynomial Poly_1; + //! instance of a bivariate polynomial + typedef CGAL::Polynomial Poly_2; + + //! conversion from the basic number type to doubles + typename NiX::NT_traits::To_double to_double; + //! conversion from the basic number type to integers + typename Renderer_traits::To_integer to_integer; + //! conversion from \c Integer type to built-in integer + typename Renderer_traits::To_machine_int to_int; + //! conversion from \c Rational type to used number type + typename Renderer_traits::From_rational from_rat; + //! conversion from the basic NT to \c Rational + typename Renderer_traits::To_rational to_rat; + //! makes the result exact after inexact operation (applicable only for + //! exact number types + typename Renderer_traits::Make_exact make_exact; + + //! returns \c true whenever a precision limit of used number type is + //! reached + typename Renderer_traits::Precision_limit limit; + //! maximum level of subdivision, dependent on used data type + static const unsigned MAX_SUBDIVISION_LEVEL = + Renderer_traits::MAX_SUBDIVISION_LEVEL; + //! isolated point + typedef Intern::Range_templ Isolated_point; + //! a list of isolated points + typedef std::list Isolated_points; + + //! map container element's type for maintaining a list of cache instances + typedef std::pair LRU_entry; + + //! LRU list used for effective cache switching + typedef boost::multi_index::multi_index_container< + LRU_entry, boost::multi_index::indexed_by< + boost::multi_index::sequenced<>, + boost::multi_index::ordered_unique< + BOOST_MULTI_INDEX_MEMBER(LRU_entry,int,first) > > > + LRU_list; + + //!@} +public: + //! \name Constructors + //!@{ + + //! default constructor + Subdivision_1() : initialized(false), one(1), polynomial_set(false) {} + //!@} + +public: + //! \name public methods + //!@{ + + //! sets up drawing window and pixel resolution + void setup(const double& x_min_,const double& y_min_, + const double& x_max_,const double& y_max_, + int res_w_, int res_h_) { + initialized = engine.setup(x_min_, y_min_, x_max_, y_max_, res_w_, + res_h_); + } + + //! sets up the underlying polynomial + void set_polynomial(const Poly_dst_2& poly) + { + input_poly = poly; + polynomial_set = true; + select_cache_entry(); + } + //! returns the underlying polynomial + Poly_dst_2 get_polynomial() const + { + return input_poly; + } + + //! \brief returns drawing window boundaries + void get_window(double& x_min_, double& y_min_, double& x_max_, + double& y_max_) const + { + x_min_ = to_double(engine.x_min); + x_max_ = to_double(engine.x_max); + y_min_ = to_double(engine.y_min); + y_max_ = to_double(engine.y_max); + } + + //! \brief returns pixel resolution + void get_resolution(int& res_w_, int& res_h_) const + { + res_w_ = engine.res_w; + res_h_ = engine.res_h; + } + + //! \brief the main rendering procedure, copies a set of pixel coordinates + //! to the output iterator \c oi + template + OutputIterator draw(OutputIterator oi); + + ~Subdivision_1() + { + cache_list.clear(); + } + + //!@} +private: + //! \name Private methods + //!@{ + + //! makes certain cache entry active + void select_cache_entry(); + + //! refines isolated point intervals until pixel size + void refine_points(int var, const NT& key, const Poly_1& poly, + Isolated_points& points); + + //! isolates curve points along a sweep-line + bool isolate_recursive(int var, const NT& beg, const NT& end, + const NT& clip, const Poly_1& poly, Isolated_points& points); + + //! \brief returns whether a polynomial has zero over an interval, + //! we are not interested in concrete values + //! + //! a bit mask \c check indicates which boundaries are to be computed + //! 0th bit - of a polynomial itself (0th derivative, default) + //! 1st bit - of the first derivative (sets \c first_der flag) + //! 2nd bit - of the second derivative + bool get_range_1(int var, const NT& lower, const NT& upper, const NT& key, + const Poly_1& poly, int check = 1); + + //!@} +private: + //! \name Private properties + //!@{ + + Poly_dst_2 input_poly; //! supporting polynomial + + //! an instance of rendering engine + Intern::Curve_renderer_internals engine; + + int cache_id; //! index of currently used cache instance + LRU_list cache_list; //! list of indices of cache instances + + bool initialized; //! indicates whether the renderer has been initialized + //! with correct parameters + const Integer one; //! just "one" + + bool polynomial_set; //! true if input polynomial was set + + //!@} +}; // class Subdivision_1<> + +//! \brief main rasterization procedure, copies in the the output iterator +//! \c oi a set of pixel coordinates + +template +template +OutputIterator Subdivision_1::draw( + OutputIterator oi) +{ + if(!initialized||!polynomial_set) + return oi; + //std::cout << "resolution: " << res_w << " x " << res_h << std::endl; + //std::cout << "box: [" << x_min << "; " << y_min << "]x[" << x_max << "; " + // << y_max << "]" << std::endl; + + Isolated_points points; + Poly_1 poly; + int x, y; + + for(x = 0; x < engine.res_w; x++) { + NT key = NT(x); + points.clear(); + + engine.get_precached_poly(SoX_Y_RANGE, key, 0, poly); + isolate_recursive(SoX_Y_RANGE, NT(0), NT(engine.res_h), key, poly, + points); + refine_points(SoX_Y_RANGE, key, poly, points); + + typename Isolated_points::iterator it = points.begin(); + while(it != points.end()) { + // std::cout << "(" << x << "; " << y << ") "; + y = engine.res_h - (int)floor((*it).left); + //if((*it).sign_change) + //painter->drawPoint(x, y); + //else { + *oi++ = std::make_pair(x, y); + //int yend = engine.res_h - (int)floor((*it).right); + //painter->moveTo(x,y); + //painter->lineTo(x,yend); + it++; + } + } + + for(y = 0; y < engine.res_h; y++) { + NT key = NT(y); + points.clear(); + + engine.get_precached_poly(SoX_X_RANGE, key, 0, poly); + isolate_recursive(SoX_X_RANGE, NT(0), NT(engine.res_h), key, poly, + points); + refine_points(SoX_X_RANGE, key, poly, points); + + typename Isolated_points::iterator it = points.begin(); + while(it != points.end()) { + // std::cout << "(" << x << "; " << y << ") "; + x = (int)floor((*it).left); + *oi++ = std::make_pair(x, engine.res_h - y); + + //int yend = res_h - (int)floor((*it).right); + //painter->moveTo(x,y); + //painter->lineTo(x,yend); + + it++; + } + } + //std::cout << "exit normal" << std::endl; + return oi; +} + +template +void Subdivision_1::refine_points(int var, + const NT& key, const Poly_1& poly, Isolated_points& points) +{ + NT dist = NT(1)/NT(SoX_REFINE_ISOLATED_POINTS); + make_exact(dist); + typename Isolated_points::iterator it = points.begin(); + while(it != points.end()) { + NT l = (*it).left, r = (*it).right; + //Gfx_OUT << l << " " << r << "\n"; + int eval_l = engine.evaluate_generic(var, l, key, poly), + eval_r = engine.evaluate_generic(var, r, key, poly); + if((eval_l^eval_r)==0) { // no sign change - interval is too small + /*if(eval_l==0) + (*it).left = (*it).right = l; + else if(eval_r==0) + (*it).left = (*it).right = r;*/ + (*it).sign_change = false; + //std::cout << "no sign change\n"; + it++; + continue; + } + while(r - l > dist) { + NT mid = (l+r)/2; + make_exact(mid); + int eval_m = engine.evaluate_generic(var, mid, key, poly); + if(eval_m == 0) + l = r = mid; + else if(eval_m == eval_l) + l = mid; + else + r = mid; + } + (*it).left = l; + (*it).right = r; + it++; + } +} + +template +bool Subdivision_1::isolate_recursive + (int var, const NT& beg, const NT& end, const NT& key, const Poly_1& poly, + Isolated_points& points) +{ + // y_clip: 0 for top boundary, res_h-1 for bottom boundary + // beg, end: 0 and res_w-1 respectively + + if(!get_range_1(var, beg, end, key, poly, 3)) + return false; + + if(!engine.first_der) { + //std::cout << "interval found: " << (x_end - x_beg) << std::endl; + points.push_back(Isolated_point(beg, end)); + return true; + } + + NT dist = NT(1)/NT(SoX_REFINE_ISOLATED_POINTS); + make_exact(dist); + if(end - beg < dist) {// pixel size is reached + //std::cout << "WARNING: roots are too close..\n"; + points.push_back(Isolated_point(beg, end)); + return true; + } + NT mid = (beg + end)/2; + make_exact(mid); + + isolate_recursive(var, beg, mid, key, poly, points); + isolate_recursive(var, mid, end, key, poly, points); + return true; +} + +//! \brief returns whether a polynomial has zero at a given interval, +//! we are not interested in concrete values +//! +//! if \t der_check is set a range for the first derivative is also computed +template +inline bool Subdivision_1::get_range_1(int var, + const NT& lower, const NT& upper, const NT& key, const Poly_1& poly, + int check) +{ + bool res = engine.get_range_QF_1(var, lower, upper, key, poly, check); + return res; +} + +//! \brief switches to a certain cache instance depending on currently used +//! algebraic curve +template +void Subdivision_1::select_cache_entry() +{ + cache_id = 0; // no cache is currently used + + engine.select_cache_entry(cache_id); + engine.precompute(input_poly); +} + +CGAL_END_NAMESPACE + +#endif // CGAL_CKVA_SUBDIVISION_1_H diff --git a/Curved_kernel_via_analysis_2/include/CGAL/Curved_kernel_via_analysis_2/gfx/Subdivision_2.h b/Curved_kernel_via_analysis_2/include/CGAL/Curved_kernel_via_analysis_2/gfx/Subdivision_2.h new file mode 100644 index 00000000000..675f9020b18 --- /dev/null +++ b/Curved_kernel_via_analysis_2/include/CGAL/Curved_kernel_via_analysis_2/gfx/Subdivision_2.h @@ -0,0 +1,549 @@ +// TODO: Add licence +// +// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE +// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. +// +// $URL:$ +// $Id: $ +// +// +// Author(s) : Pavel Emeliyanenko +// +// ============================================================================ + +/*!\file CGAL/Curved_kernel_via_analysis_2/gfx//Subdivision_2.h + * \brief definition of Subdivision_2<> + * 2D space subdivision for rasterization of planar curves + */ + +#ifndef CGAL_CKVA_SUBDIVISION_2_H +#define CGAL_CKVA_SUBDIVISION_2_H 1 +#warning "this file is considered obsolete" + +#include +#include +#include + +#include + +CGAL_BEGIN_NAMESPACE + +namespace CGALi { + +template +struct Range_templ +{ + Range_templ() { } + Range_templ(const NT& l_, const NT& u_) : lower(l_), upper(u_) { } + NT lower, upper; +}; + +} // namespace CGALi + +/*!\brief + * The class template \c Subdivision_2 and its associate functions. + * + * The class implements a space method to plot algebraic curves, we use Affine + * Arithmetic with recursive derivative information + */ +template +class Subdivision_2 +{ +public: + //! \name public typedefs + //!@{ + //! this instance's first template argument + typedef NT_ NT; + //! this instance's second template argument + typedef Algebraic_curve_2_ Algebraic_curve_2; + //! specialized integer number type + typedef typename SoX::Curve_renderer_traits::Integer Integer; + //! instance of univariate polynomial + typedef CGAL::Polynomial Poly_1; + //! instance of bivariate polynomial + typedef CGAL::Polynomial Poly_2; + //! rational number type instance + typedef typename Algebraic_curve_2::Rational Rational; + //! instance of input bivariate polynomial + typedef typename Algebraic_curve_2::Poly_d Poly_input_2; + //! instance of input univariate polynomuial + typedef typename Poly_input_2::NT Poly_input_1; + //! instance of coefficient number type + typedef typename Poly_input_1::NT Coeff; + //! container used to store coefficient sequence + typedef typename Poly_1::Vector Vector_1; + //! container's const iterator (random access) + typedef typename Poly_1::const_iterator const_iterator_1; + //! container used to store univariate polynomials sequence + typedef typename Poly_2::Vector Vector_2; + //! container's const iterator (random access) + typedef typename Poly_2::const_iterator const_iterator_2; + //@} +private: + //! \name private typedefs + //@{ + //! conversion from the basic number type to doubles + typename NiX::NT_traits::To_double to_double; + //! conversion from the basic number type to integers + typename SoX::Curve_renderer_traits::To_integer to_integer; + //! conversion from \c Integer type to built-in integer + typename SoX::Curve_renderer_traits::To_machine_int to_int; + //! conversion from \c Rational type to used number type + typename SoX::Curve_renderer_traits::From_rational from_rat; + //! makes the result exact after inexact operation (applicable only for + //! exact number types + typename SoX::Curve_renderer_traits::Make_exact make_exact; + //! an interval + typedef Intern::Range_templ Range; + //! affine form instance + typedef SoX::Affine_form Affine_form; + //@} +public: + //! \name Constructors + //@{ + //! default constructor + Subdivision_2() : initialized(false), polynomial_set(false) {} + //@} +public: + //! \name public methods + //@{ + //! specifies drawing window and pixel resolution + void setup(const double& x_min_,const double& y_min_, + const double& x_max_,const double& y_max_, + int res_w_, int res_h_) + { + x_min = static_cast(x_min_); + x_max = static_cast(x_max_); + y_min = static_cast(y_min_); + y_max = static_cast(y_max_); + res_w = res_w_; + res_h = res_h_; + if(x_min >= x_max||y_min >= y_max||res_w < 5||res_h < 5||res_w > 1024|| + res_h > 1024) { + std::cout << "Incorrect setup parameters" << std::endl; + initialized = false; + return; + } + pixel_w = (x_max - x_min) / res_w; + pixel_h = (y_max - y_min) / res_h; + make_exact(pixel_w); + make_exact(pixel_h); + initialized = true; + } + //! sets up the underlying polynomial + void set_polynomial(const Poly_input_2& poly) + { + input_poly = poly; + precompute(); + } + //! returns the underlying polynomial + Poly_input_2 get_polynomial() const + { + return input_poly; + } + //! \brief returns drawing window boundaries + void get_window(double& x_min_, double& y_min_, double& x_max_, + double& y_max_) const + { + x_min_ = to_double(x_min); + x_max_ = to_double(x_max); + y_min_ = to_double(y_min); + y_max_ = to_double(y_max); + } + //! \brief returns pixel resolution + void get_resolution(int& res_w_, int& res_h_) const + { + res_w_ = res_w; + res_h_ = res_h; + } + //! \brief the main rendering procedure + void draw(QPainter *painter_); + //@} +private: + //! \name Private methods + //@{ + //! precomputes polynomials and derivative coefficients + void precompute(); + //! \brief switches to another cache instance depending on the + //! supporting curve of a segment + //! \brief evalutates the ith derivative at certain x + //! + //! \c cache_it - an intetator pointing to the end of an array of + //! polynomial coefficients, \c der_it - an iterator for derivative + //! coefficients + NT evaluate_der(const_iterator_1 der_it, const_iterator_1 begin, + const_iterator_1 cache_it, const NT& x) + { + NT val((*cache_it--) * (*der_it)); + while((der_it--)!=begin) + val = val * x + (*cache_it--) * (*der_it); + return val; + } + //! evalutates a function at a certain x + NT evaluate(const Poly_1& poly, const NT& x) + { + const_iterator_1 it = poly.end() - 1, begin = poly.begin(); + NT val(*it); + while((it--)!=begin) + val = val * x + (*it); + return val; + } + //! \brief computes the ranges of univariate polynomial values over an interval + void get_range_1(int var, const NT& lower, const NT& upper, + const Poly_1& poly, NT& l, NT& h); + //! Affine Arithmetic with Recursive Derivative information + //! for univariate case + void get_range_AARD_1(int var, const NT& lower, const NT& upper, + const Poly_1& poly, NT& l, NT& h); + //! \brief Recursive Taylor, bivariate case + //! + //! returns a range of polynomial values as Affine_form + void get_range_RT_2(const NT& x_low, const NT& x_high, const NT& y_low, + const NT& y_high, int depth, int index, Affine_form& res); + //! checks a rectangular area with 2D range analysis: either discrads it, + //! subdivides further or draws a pixel + void quad_tree(const NT& x_low, const NT& x_high, const NT& y_low, + const NT& y_high); + //! \brief recursive quad tree subdivision + void subdivide(const NT& x_low, const NT& x_high, const NT& y_low, + const NT& y_high); +public: + //! destructor + ~Subdivision_2() + { + + } + //@} +private: + //! \name Private properties + //@{ + NT x_min, x_max, y_min, y_max; //! drawing window boundaries + int res_w, res_h; //! pixel resolution + NT pixel_w, pixel_h; //! pixel dimensions w.r.t. resolution + Poly_2 coeffs_x, coeffs_y; //! f(x(y)) / f(y(x)) + Poly_input_2 input_poly; + Vector_2 der_x, der_y; //! derivative coefficients df/dx (df/dy) + std::vector mixed_derivatives; + int max_deg; //! the maximal degree w.r.t. x and y + bool initialized, polynomial_set; + QPainter *painter; + //@} +}; // class Subdivision_2<> + +//! \brief main rasterization procedure +template +void Subdivision_2::draw(QPainter *painter_) +{ + if(!initialized||!polynomial_set||painter_==NULL) + return; + painter = painter_; + //std::cout << " P(x(y)): " << coeffs_x << std::endl; + //std::cout << " P(y(x)): " << coeffs_y << std::endl; + std::cout << "resolution: " << res_w << " x " << res_h << std::endl; + //std::cout << "box: [" << x_min << "; " << y_min << "]x[" << x_max << "; " + // << y_max << "]" << std::endl; + + //quad_tree(x_min/15, x_max/15, y_min/15, y_max/15); + quad_tree(x_min, x_max, y_min, y_max); + //std::cout << "exit normal" << std::endl; +} + +//! \brief checks a rectangular area with 2D range analysis: either discrads it, +//! subdivides further or draws a pixel +template +void Subdivision_2::quad_tree(const NT& x_low, + const NT& x_high, const NT& y_low, const NT& y_high) +{ + Affine_form res; + NT lower, upper; + get_range_RT_2(x_low, x_high, y_low, y_high, 0, 0, res); + res.convert(lower, upper); + if(lower*upper > 0) + return; + if(x_high - x_low <= pixel_w&&y_high - y_low <= pixel_h) { + NT x = (x_low), y = (y_low); + int pix_x = static_cast(NiX::to_double((x - x_min) / pixel_w)), + pix_y = static_cast(NiX::to_double((y - y_min) / pixel_h)); + painter->drawPoint(pix_x, res_h - pix_y); + //painter->drawEllipse(pix_x-2,res_h-pix_y-2,4,4); + } + else + subdivide(x_low, x_high, y_low, y_high); +} + +//! \brief recursive quad tree subdivision +template +void Subdivision_2::subdivide(const NT& x_low, + const NT& x_high, const NT& y_low, const NT& y_high) +{ + NT x_mid = (x_low + x_high)/2, y_mid = (y_low + y_high)/2; + quad_tree(x_low,x_mid,y_low,y_mid); + quad_tree(x_low,x_mid,y_mid,y_high); + quad_tree(x_mid,x_high,y_mid,y_high); + quad_tree(x_mid,x_high,y_low,y_mid); +} + +//! \brief Recursive Taylor, bivariate case +template +void Subdivision_2::get_range_RT_2( + const NT& x_low, const NT& x_high, const NT& y_low, const NT& y_high, + int depth, int index, Affine_form& res) +{ + //std::cout << "range for [" << x_low << "; " << y_low << "]x[" << + //x_high << + //"; " << y_high << "]: (" << low << "; " << high << ")" << std::endl; + typename std::vector::const_iterator der_it = + mixed_derivatives.begin() + index; + if((*der_it).degree()==0) { + NT c = (*der_it).lcoeff().lcoeff(); + res = Affine_form(c,c); + return; + } else if((*der_it).degree()==1) { + Poly_1 p_x_low = NiX::substitute_x((*der_it), x_low), + p_x_high = NiX::substitute_x((*der_it), x_high); + NT eval_1 = p_x_low.evaluate(y_low), + eval_2 = p_x_low.evaluate(y_high), + eval_3 = p_x_high.evaluate(y_low), + eval_4 = p_x_high.evaluate(y_high); + NT min1 = eval_1, max1 = eval_2; + if(min1 > eval_2) { + min1 = eval_2; + max1 = eval_1; + } + NT min2 = eval_3, max2 = eval_4; + if(min2 > eval_4) { + min2 = eval_4; + max2 = eval_3; + } + if(min1 > min2) + min1 = min2; + if(max1 < max2) + max1 = max2; + res = Affine_form(min1,max1); + return; + } + + if(depth >= 4) + { + NT l, h; + std::vector x_forms; + const_iterator_2 it_2 = (*der_it).begin(); + while(it_2 != (*der_it).end()) + { + // it_2 - a poly in x_range + if((*it_2).is_zero()) + l = h = 0; + else + get_range_1(X_RANGE, x_low, x_high, *it_2, l, h); + x_forms.push_back(Affine_form(l, h)); + it_2++; + } + typename std::vector::const_iterator it = + x_forms.end()-1; + Affine_form y_form(y_low, y_high); + res = Affine_form(*it); + while((it--) != x_forms.begin()) + res = res * y_form + (*it); + //res.convert(lower, upper); + return; + } + + NT x0 = (x_low+x_high)/2, y0 = (y_low+y_high)/2, x1 = (x_high-x_low)/2, + y1 = (y_high-y_low)/2; + Affine_form one1(-1,1), one2(-1,1), one3(-1,1), + zero1(0,1), zero2(0,1), fxx, fxy, fyy; + NT eval_f = NiX::substitute_xy(*der_it, x0, y0), + eval_fx = NiX::substitute_xy(*(der_it+depth+1), x0, y0)*x1, + eval_fy = NiX::substitute_xy(*(der_it+depth+2), x0, y0)*y1; + res = eval_f + eval_fx*one1 + eval_fy*one2; + + int idx = index+2*depth+3; + if(idx < (int)mixed_derivatives.size()) { + get_range_RT_2(x_low, x_high, y_low, y_high, depth+2, idx, fxx); + get_range_RT_2(x_low, x_high, y_low, y_high, depth+2, idx+1, fxy); + get_range_RT_2(x_low, x_high, y_low, y_high, depth+2, idx+2, fyy); + res = res + ((x1*x1/2*zero1)*fxx) + ((y1*y1/2*zero2)*fyy) + + ((x1*y1*one3)*fxy); + } + + //res.convert(lower, upper); + //std::cout << "range for depth = " << depth << " index = " << index << + //" [" << lower << "; " << upper << "]" << std::endl; +} + +//! \brief computes the range of polynomial values \c f([lower,upper]) using +//! Affine Arithmetic with Recursive Derivative information +//! +//! \c var = \c X_RANGE: \c y is fixed - polynomial is given in \c x +//! coordinates, otherwise: \c x is fixed - polynomial is given in \c y +//! coordinates +template +void Subdivision_2::get_range_AARD_1(int var, + const NT& l_, const NT& r_, const Poly_1& poly, NT& l1, NT& h1) +{ + Vector_2 *der = &der_y; + if(var == X_RANGE) + der = &der_x; + NT low, up, l(l_), r(r_); + const_iterator_2 der_it_2 = der->end()-1; + const_iterator_1 der_it, cache_it, begin; + if(poly.degree()==0) { + l1 = h1 = poly.lcoeff(); + return; + } + low = up = poly.lcoeff() * (*der_it_2).lcoeff(); + if(l > r) { + l = r_; + r = l_; + } + Affine_form f(l,r); + /*while((der_it_2--)!=der->begin()) { + // iterate through derivative coefficients + der_it = (*der_it_2).end()-1; + begin = (*der_it_2).begin(); + cache_it = poly.end()-1; // iterate through precomputed y-values + // if a derivative does not straddle zero we can + // calculate the exact boundaries for f(x) + if(low * up > 0) { + v1 = v2 = (*cache_it--) * (*der_it); + // calculate the ith derivative at xa and xb + while((der_it--)!=begin) { + v1 = v1 * l + (*cache_it) * (*der_it); + v2 = v2 * r + (*cache_it--) * (*der_it); + } + if(low < 0 && up < 0) { + v = v1; + v1 = v2; + v2 = v; + } + low = v1; + up = v2; + } else { // use affine arithmetic to compute bounds + Affine_form range(((*cache_it--) * (*der_it))); + // calculate the ith derivative using affine arithmetic + while((der_it--)!=begin) + range = (*cache_it--) * (*der_it) + (range * f); + range.convert(low,up); + } + } + if(low * up >= 0) { + Gfx_OUT << "+ " << std::endl; + v1 = evaluate(poly, l); + v2 = evaluate(poly, r); + //Gfx_OUT << "v1 = " << v1 << " v2 = " << v2 << std::endl; + if(low < 0 && up < 0) { + v = v1; + v1 = v2; + v2 = v; + } + low = v1; + up = v2; + } else */{ // use affine arithmetic to compute bounds + cache_it = poly.end()-1; + begin = poly.begin(); + Affine_form res(*cache_it); + while((cache_it--)!=begin) + res = (*cache_it) + (res * f); + res.convert(low,up); + } + l1 = low; + h1 = up; +//std::cout << "AARD bounds: [" << low << "; " << up << "]" << std::endl; +} + +//! \brief returns whether a polynomial has zero at a given interval, +//! since we are not interested in concrete values +//! +//! flag \t der_check forces to calculate boundaries for the first +//! derivative instead of function itself. Caching is used if +//! appropriate +template +void Subdivision_2::get_range_1(int var, + const NT& lower, const NT& upper, const Poly_1& poly, NT& l, NT& h) +{ + //std::cout << "range for: [" << lower << "; " << upper << "] poly: " << + //poly << std::endl; + get_range_AARD_1(var, lower, upper, poly, l, h); +} + +//! precomputes some initial data, necessary for subsequent computations +template +void Subdivision_2::precompute() +{ + Intern::Max_coeff max_coeff; + Coeff mx1 = max_coeff(input_poly); + Intern::Transform::From_rational, + typename SoX::Curve_renderer_traits::Make_exact> tr; + tr.max = Rational(mx1); + coeffs_y = tr(input_poly); + // x - outer variable, y - inner + coeffs_x = transpose_bivariate_polynomial(coeffs_y); + int degree_x = coeffs_x.degree(), + degree_y = coeffs_y.degree(); + der_x.clear(); + der_y.clear(); + int i, j; + max_deg = degree_x; + if(degree_y > max_deg) + max_deg = degree_y; + NT *X = new NT[max_deg]; + NT det(1.0); + std::cout << "start" << std::endl; + for(i = 0; i < degree_x; i++) { + if(i != 0) + det = X[0]; + for(j = 1; j <= degree_x - i; j++) + if(i == 0) + X[j-1] = j; + else { + X[j-1] = X[j] * j / det; // divide by the lowest coefficient ? + make_exact(X[j-1]); + } + der_x.push_back(Poly_1(X,(X + degree_x - i))); + } + for(i = 0; i < degree_y; i++) { + if(i != 0) + det = X[0]; + for(j = 1; j <= degree_y - i; j++) + if(i == 0) + X[j-1] = j; // divide by the lowest coefficient ? + else { + X[j-1] = X[j] * j / det; + make_exact(X[j-1]); + } + der_y.push_back(Poly_1(X,(X + degree_y - i))); + } + delete []X; + mixed_derivatives.clear(); + // f, fx, fy, fxx, fxy, fyy, fxxx, fxxy, fxyy, fyyy, ... + mixed_derivatives.push_back(coeffs_y); + int idx = 0; + for(i = 1; i <= max_deg; i++) + { + mixed_derivatives.push_back(NiX::diff_x(mixed_derivatives[idx])); + for(j = 0; j < i; j++) +// compute fx^(i-j)y^(j), i.e. (i-j) times derivate by x; and j times +// derivate by y + { + Poly_2 p = mixed_derivatives[idx+j]; + p.diff(); + mixed_derivatives.push_back(p); + } + idx += i; + } + std::cout << "finished" << std::endl; + polynomial_set = true; + /*typename std::vector::const_iterator der_it = + mixed_derivatives.end()-1; + for(i = max_deg; i >= 0; i--) + { + std::cout << i << "th mixed derivatives: " << std::endl; + for(j = 0; j < i+1; j++, der_it--) + std::cout << *der_it << std::endl; + }*/ +} + +CGAL_END_NAMESPACE + +#endif // CGAL_CKVA_SUBDIVISION_2_H