songsenand
/
cgal
mirror of https://github.com/CGAL/cgal
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

Curve renderer migrated

This commit is contained in:
Pavel Emeliyanenko 2008-07-02 17:01:02 +00:00
parent 3a6da39537
commit be8878919f
8 changed files with 5745 additions and 0 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

7
.gitattributes vendored
View File

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

286
Curved_kernel_via_analysis_2/include/CGAL/Curved_kernel_via_analysis_2/Curve_renderer_facade.h Normal file
View File

@ -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 <asm@mpi-sb.mpg.de>
//
// ============================================================================
/*!\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<void>(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<void>(0)
#endif
#endif
#include <CGAL/basic.h>
#include <CGAL/Bbox_2.h>
#include <CGAL/Arithmetic_kernel.h>
#include <CGAL/Curved_kernel_via_analysis_2/gfx/Curve_renderer_2.h>
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 CurvedKernelViaAnalysis_2, class Float_>
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<BFI>::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<Curved_kernel_via_analysis_2, Float>
Default_renderer_2;
#ifdef CGAL_CKVA_USE_MULTIPREC_ARITHMETIC
//! the curve renderer instantiated with BigFloat
typedef Curve_renderer_2<Curved_kernel_via_analysis_2, Bigfloat>
Bigfloat_renderer_2;
#endif
#ifdef CGAL_CKVA_USE_RATIONAL_ARITHMETIC
//! the curve renderer instantiated with Rationals
typedef Curve_renderer_2<Curved_kernel_via_analysis_2, Rational>
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
*
* <tt>[x_min; y_min]x[x_max; y_max]</tt> - 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<CGALi::Coord_vec_2>& points,
std::pair<CGALi::Coord_2, CGALi::Coord_2>& 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<typename NiX::NT_traits<Float>::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<typename NiX::NT_traits<Float>::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

171
Curved_kernel_via_analysis_2/include/CGAL/Curved_kernel_via_analysis_2/Qt_widget_Curve_renderer_2.h Normal file
View File

@ -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 <asm@mpi-sb.mpg.de>
//
// ============================================================================
/*!\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 <vector>
#include <qpainter.h>
#include <CGAL/IO/Qt_widget.h>
#include <CGAL/Interval_traits.h>
#include <CGAL/Bigfloat_interval_traits.h>
#include <CGAL/convert_to_bfi.h>
#include <CGAL/Arithmetic_kernel.h>
#include <CGAL/Curved_kernel_via_analysis_2/Curve_renderer_facade.h>
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 CurvedKernelViaAnalysis_2>
class Renderer_instance
{
Renderer_instance() { // private constructor
}
public:
typedef typename CurvedKernelViaAnalysis_2::Curve_kernel_2::Coefficient
Coefficient;
typedef CGAL::Interval_nt<true> Interval_double;
typedef Curve_renderer_facade<CurvedKernelViaAnalysis_2, Interval_double>
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<Coord_vec_2> points;
std::pair<Coord_2, Coord_2> end_points;
gfx_instance().draw(arc, points, end_points);
if(points.empty())
return;
std::list<Coord_vec_2>::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 <class CKvA>
CGAL::Bbox_2 Renderer_instance<CKvA>::bbox(0.0, 0.0, 0.0, 0.0);
} // namespace CGALi
/*! \brief
* outputs a curve arc to \c Qt_widget
*/
template <class CurvedKernelViaAnalysis_2>
Qt_widget& operator << (Qt_widget& ws,
const CGALi::Arc_2<CurvedKernelViaAnalysis_2>& arc) {
CGALi::Renderer_instance<CurvedKernelViaAnalysis_2>::draw(ws, arc);
return ws;
}
/*! \brief
* outputs a curve point to \c Qt_widget
*/
template <class CurvedKernelViaAnalysis_2>
Qt_widget& operator << (Qt_widget& ws,
const CGALi::Point_2<CurvedKernelViaAnalysis_2>& pt) {
CGALi::Renderer_instance<CurvedKernelViaAnalysis_2>::draw(ws, pt);
return ws;
}
CGAL_END_NAMESPACE
#endif // CGAL_QT_WIDGET_CURVE_RENDERER_2_H

2653
Curved_kernel_via_analysis_2/include/CGAL/Curved_kernel_via_analysis_2/gfx/Curve_renderer_2.h Normal file
View File

File diff suppressed because it is too large Load Diff

1072
Curved_kernel_via_analysis_2/include/CGAL/Curved_kernel_via_analysis_2/gfx/Curve_renderer_internals.h Normal file
View File

File diff suppressed because it is too large Load Diff

575
Curved_kernel_via_analysis_2/include/CGAL/Curved_kernel_via_analysis_2/gfx/Curve_renderer_traits.h Normal file
View File

@ -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 <asm@mpi-sb.mpg.de>
//
// ============================================================================
#ifndef CGAL_CKVA_CURVE_RENDERER_TRAITS_H
#define CGAL_CKVA_CURVE_RENDERER_TRAITS_H
#include <CGAL/basic.h>
#include <CGAL/function_objects.h>
/*! \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 <class NT>
struct Max_coeff
{
template <class X>
NT operator()(const CGAL::Polynomial<X>& p) const
{
typename CGAL::Polynomial<X>::const_iterator it = p.begin();
Max_coeff<NT> 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 <class Result, class Input>
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<Input, Result>::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
*
* <tt>Op: InputPoly_2::Innermost_coefficient ->
* OutputPoly_2::Innermost_coefficient</tt>
*/
template <class OutputPoly_2, class InputPoly_2, class Op>
struct Transform {
typedef InputPoly_2 first_argument_type;
typedef Op second_argument_type;
typedef OutputPoly_2 result_type;
template <class X>
OutputPoly_2 operator()(const CGAL::Polynomial<X>& p, Op op = Op()) const {
Transform<typename OutputPoly_2::NT, typename InputPoly_2::NT, Op> 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<InputPoly_2>::Type& x, Op op)
const {
return static_cast<OutputPoly_2>(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 <class Coeff_, class Integer_, class Rational_>
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 <class Float>
struct Rat_to_float {
typedef Float result_type;
template <class X, class Y>
Float operator()(const Sqrt_extension<X, Y>& x) const {
typename CGAL::Coercion_traits<Sqrt_extension<X, Y>, Float>::Cast
cast;
return cast(x);
}
Float operator()(const Rational& x) const {
return static_cast<Float>(x);
}
};
//! provided for convenience when there exists an implicit coercion
//! between number types
template <class To>
struct Implicit_coercion {
typedef To result_type;
template <class From>
To operator()(const From& x) const {
return static_cast<To>(x);
}
};
struct Float_to_int {
typedef int result_type;
template <class Float>
int operator()(const Float& x) const
{ return static_cast<int>(std::floor(CGAL::to_double(x))); }
//return static_cast<int>(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<CGAL::Polynomial<Coeff> > result_type;
template <class Extended>
inline result_type operator()(const
CGAL::Polynomial<CGAL::Polynomial<Extended> >& poly) const {
//!@todo: use Rat_to_float functor instead of coercion traits ?
//! need some sort of polymorphism..
typedef typename CGAL::Coercion_traits<Extended, Coeff>::Cast
Cast;
Transform<result_type,
CGAL::Polynomial<CGAL::Polynomial<Extended> >, 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 <class Float>
std::size_t operator()(const Float& key) const {
const long_long *hk = reinterpret_cast<const long_long *>(&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 <class Float>
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 <class Float>
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 <class Float, class Rational>
struct Curve_renderer_traits
{ };
#ifdef CGAL_USE_CORE
//! Specialization for \c CGAL::Interval_nt<true>
template <>
struct Curve_renderer_traits<CGAL::Interval_nt<true>, CORE::BigRat > :
public Curve_renderer_traits_base<CGAL::Interval_nt<true>, int,
CORE::BigRat> {
typedef double Float;
struct Rat_to_float {
typedef Float result_type;
template <class Extended>
Float operator()(const Extended& x) const {
return CGAL::to_double(x);
}
};
typedef Implicit_coercion<double> Float_to_rat;
struct Rat_to_integer {
typedef Rational argument_type;
typedef Integer result_type;
Integer operator()(const Rational& x) const {
return static_cast<int>(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 <class Float>
bool operator()(const Float& x) const
{ return (CGAL_ABS(x) <= 1e-16); }
};
};
//! Specialization for \c CORE::BigFloat
template <>
struct Curve_renderer_traits<CORE::BigFloat, class CORE::BigRat>
: public Curve_renderer_traits_base<CORE::BigFloat, CORE::BigInt,
CORE::BigRat> {
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<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<std::size_t>(&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<CORE::BigRat, CORE::BigRat> :
public Curve_renderer_traits_base<CORE::BigRat, CORE::BigInt,
CORE::BigRat> {
typedef CORE::BigRat Float;
typedef Rat_to_float<Float> Rat_to_float;
typedef Implicit_coercion<Float> 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<std::size_t>(&rep);
return ret;
}
};
};
#endif // CGAL_USE_CORE
#ifdef CGAL_USE_LEDA
template <>
struct Curve_renderer_traits<CGAL::Interval_nt<true>, leda::rational > :
public Curve_renderer_traits_base<CGAL::Interval_nt<true>, int,
leda::rational> {
typedef double Float;
struct Rat_to_float {
typedef Float result_type;
template <class Extended>
Float operator()(const Extended& x) const {
return CGAL::to_double(x);
}
};
typedef Implicit_coercion<double> Float_to_rat;
struct Rat_to_integer {
typedef Rational argument_type;
typedef Integer result_type;
Integer operator()(const Rational& x) const {
return static_cast<int>(std::floor(CGAL::to_double(x)));
//return static_cast<int>(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 <class Float>
bool operator()(const Float& x) const
{ return (CGAL_ABS(x) <= 1e-16); }
};
};
//! Specialization for \c leda::bigfloat
template <>
struct Curve_renderer_traits<leda::bigfloat, class leda::rational>
: public Curve_renderer_traits_base<leda::bigfloat, leda::integer,
leda::rational> {
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 <class X, class Y>
Float operator()(const Sqrt_extension<X, Y>& x) const {
typename CGAL::Coercion_traits<Sqrt_extension<X, Y>, 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<Float>(x.numerator()) /
static_cast<Float>(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<CGAL::Polynomial<Coeff> > result_type;
template <class Extended>
inline result_type operator()(const
CGAL::Polynomial<CGAL::Polynomial<Extended> >& poly) const {
Transform<result_type,
CGAL::Polynomial<CGAL::Polynomial<Extended> >,
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<std::size_t>(
key.get_significant_length());
}
};
};
//! Specialization for \c leda::rational
template <>
struct Curve_renderer_traits<leda::rational, leda::rational> :
public Curve_renderer_traits_base<leda::rational, leda::integer,
leda::rational> {
typedef leda::rational Float;
typedef Rat_to_float<Float> Rat_to_float;
typedef Implicit_coercion<Float> 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<std::size_t>(
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

432
Curved_kernel_via_analysis_2/include/CGAL/Curved_kernel_via_analysis_2/gfx/Subdivision_1.h Normal file
View File

@ -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 <asm@mpi-sb.mpg.de>
//
// ============================================================================
/*!\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 <vector>
#include <boost/multi_index_container.hpp>
#include <boost/multi_index/member.hpp>
//#include <boost/multi_index/hashed_index.hpp>
#include <boost/multi_index/ordered_index.hpp>
#include <boost/multi_index/sequenced_index.hpp>
#include <CGAL/Polynomial.h>
#include <CGAL/Interval_nt.h>
#include <CGAL/Arithmetic_kernel.h>
#include <CGAL/Algebraic_kernel_d/Bitstream_descartes.h>
#include <CGAL/Curved_kernel_via_analysis_2/gfx/Curve_renderer_internals.h>
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 <class NT>
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 Coeff_, class Algebraic_curve_2_>
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<Coeff> 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> Rational_poly_1;
//! a bivariate rational polynomial
typedef CGAL::Polynomial<Rational_poly_1> Rational_poly_2;
//! basic number type used in all computations
typedef typename Renderer_traits::Float NT;
//! instance of a univariate polynomial
typedef CGAL::Polynomial<Coeff> Poly_1;
//! instance of a bivariate polynomial
typedef CGAL::Polynomial<Poly_1> Poly_2;
//! conversion from the basic number type to doubles
typename NiX::NT_traits<NT>::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<NT> Isolated_point;
//! a list of isolated points
typedef std::list<Isolated_point> Isolated_points;
//! map container element's type for maintaining a list of cache instances
typedef std::pair<int,int> 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 <class OutputIterator>
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<Coeff, Algebraic_curve_2> 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 <class Coeff_, class Algebraic_curve_2_>
template <class OutputIterator>
OutputIterator Subdivision_1<Coeff_, Algebraic_curve_2_>::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 <class Coeff_, class Algebraic_curve_2_>
void Subdivision_1<Coeff_, Algebraic_curve_2_>::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 <class Coeff_, class Algebraic_curve_2_>
bool Subdivision_1<Coeff_, Algebraic_curve_2_>::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 <class Coeff_, class Algebraic_curve_2_>
inline bool Subdivision_1<Coeff_, Algebraic_curve_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);
return res;
}
//! \brief switches to a certain cache instance depending on currently used
//! algebraic curve
template <class Coeff_, class Algebraic_curve_2_>
void Subdivision_1<Coeff_, Algebraic_curve_2_>::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

549
Curved_kernel_via_analysis_2/include/CGAL/Curved_kernel_via_analysis_2/gfx/Subdivision_2.h Normal file
View File

@ -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 <asm@mpi-sb.mpg.de>
//
// ============================================================================
/*!\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 <vector>
#include <CGAL/Polynomial.h>
#include <CGAL/Interval_nt.h>
#include <CGAL/Curved_kernel_via_analysis_2/gfx/Curve_renderer_internals.h>
CGAL_BEGIN_NAMESPACE
namespace CGALi {
template <class NT>
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 NT_, class Algebraic_curve_2_>
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<NT>::Integer Integer;
//! instance of univariate polynomial
typedef CGAL::Polynomial<NT> Poly_1;
//! instance of bivariate polynomial
typedef CGAL::Polynomial<Poly_1> 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<NT>::To_double to_double;
//! conversion from the basic number type to integers
typename SoX::Curve_renderer_traits<NT>::To_integer to_integer;
//! conversion from \c Integer type to built-in integer
typename SoX::Curve_renderer_traits<NT>::To_machine_int to_int;
//! conversion from \c Rational type to used number type
typename SoX::Curve_renderer_traits<NT>::From_rational from_rat;
//! makes the result exact after inexact operation (applicable only for
//! exact number types
typename SoX::Curve_renderer_traits<NT>::Make_exact make_exact;
//! an interval
typedef Intern::Range_templ<NT> Range;
//! affine form instance
typedef SoX::Affine_form<NT> 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<NT>(x_min_);
x_max = static_cast<NT>(x_max_);
y_min = static_cast<NT>(y_min_);
y_max = static_cast<NT>(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<Poly_2> 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 <class NT, class Algebraic_curve_2_>
void Subdivision_2<NT, Algebraic_curve_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 <class NT_, class Algebraic_curve_2_>
void Subdivision_2<NT_, Algebraic_curve_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<int>(NiX::to_double((x - x_min) / pixel_w)),
pix_y = static_cast<int>(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 <class NT_, class Algebraic_curve_2_>
void Subdivision_2<NT_, Algebraic_curve_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 <class NT_, class Algebraic_curve_2_>
void Subdivision_2<NT_, Algebraic_curve_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<Poly_2>::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<Affine_form> 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<Affine_form>::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 <class NT_, class Algebraic_curve_2_>
void Subdivision_2<NT_, Algebraic_curve_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<NT> 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 <class NT_, class Algebraic_curve_2_>
void Subdivision_2<NT_, Algebraic_curve_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 <class NT_, class Algebraic_curve_2_>
void Subdivision_2<NT_, Algebraic_curve_2_>::precompute()
{
Intern::Max_coeff<Coeff> max_coeff;
Coeff mx1 = max_coeff(input_poly);
Intern::Transform<Poly_2, Coeff, Rational,
typename SoX::Curve_renderer_traits<NT>::From_rational,
typename SoX::Curve_renderer_traits<NT>::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<Poly_2>::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