diff --git a/Curved_kernel_via_analysis_2/include/CGAL/Arrangement_2l/Adjacencies_3.h b/Curved_kernel_via_analysis_2/include/CGAL/Arrangement_2l/Adjacencies_3.h new file mode 100644 index 00000000000..b0c5d6b9c4c --- /dev/null +++ b/Curved_kernel_via_analysis_2/include/CGAL/Arrangement_2l/Adjacencies_3.h @@ -0,0 +1,214 @@ +// ============================================================================ +// +// Copyright (c) 2001-2008 Max-Planck-Institut Saarbruecken (Germany). +// All rights reserved. +// +// This file is part of EXACUS (http://www.mpi-inf.mpg.de/projects/EXACUS/); +// you may redistribute it under the terms of the Q Public License version 1.0. +// See the file LICENSE.QPL distributed with EXACUS. +// +// Licensees holding a valid commercial license may use this file in +// accordance with the commercial license agreement provided with the software. +// +// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE +// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. +// +// ---------------------------------------------------------------------------- +// +// Library : SoX +// File : include/SoX/GAPS/Adjacencies_3.h +// SoX_release : $Name: $ +// Revision : $Revision$ +// Revision_date : $Date$ +// +// Author(s) : Eric Berberich +// Michael Kerber +// +// ============================================================================ + +#ifndef SoX_GAPS_ADJACENCIES_3_H +#define SoX_GAPS_ADJACENCIES_3_H 1 + +/*! \file SoX/GAPS/Adjacencies_3 + \brief Contains Adjacencies_3 class +*/ + +#include +#include + +#include +#include +#include + +namespace SoX { + + +namespace Intern { + +class Adjacencies_3_rep { + +public: + + //! the class itself + typedef Adjacencies_3_rep Self; + + //! a pair of adjacencies + typedef std::pair< int, int > Adjacency_pair; + + //! an interval + typedef std::pair< int, int > Adjacency_interval; + + //! an adjacency vector + typedef std::vector < Adjacency_pair > Adjacency_vector; + + //! Default constructor + Adjacencies_3_rep() { + } + + //! Constructor with data + Adjacencies_3_rep(Adjacency_vector data) : + _m_data(data) { + } + + Adjacency_vector _m_data; + +}; + +} // namespace Intern + +//! The Adjacency output object +class Adjacencies_3 : + public CGAL::Handle_with_policy< Intern::Adjacencies_3_rep > { + +public: + //! the class itself + typedef Adjacencies_3 Self; + + //! the rep type + typedef Intern::Adjacencies_3_rep Rep; + + //! the base type + typedef CGAL::Handle_with_policy< Rep > Base; + + //! a pair of adjacencies + typedef std::pair< int, int > Adjacency_pair; + + //! an interval + typedef std::pair< int, int > Adjacency_interval; + + //! a vector of adjacency + typedef std::vector < Adjacency_pair > Adjacency_vector; + + //! const version + typedef Adjacency_vector::const_iterator Const_adjacency_iterator; + +public: + + //! Default constructor + Adjacencies_3() : + Base(Rep()) { + } + + //! Constructor with data + Adjacencies_3(Adjacency_vector data) : + Base(Rep(data)) { + } + + //! Returns an interval [i,j] such that all sheets between + //! i and j are adjacenct to index + Adjacency_interval interval(int index) const { + + std::vector adjacent_features; + + for (int i = 0; + i < static_cast(this->ptr()->_m_data.size()); + i++) { + if (this->ptr()->_m_data[i].first == index) { + adjacent_features.push_back(this->ptr()->_m_data[i].second); + } + } + + std::sort(adjacent_features.begin(), adjacent_features.end()); + + int n = static_cast(adjacent_features.size()); + + if (n == 0) { + return std::make_pair(1,0); + } + // Make sure that the adjacent features form an interval + for (int i = 0; i < n-1; i++) { + + CGAL_assertion( adjacent_features[i] == + adjacent_features[i+1]-1 ); + + } + + return std::make_pair(adjacent_features[0], + adjacent_features[n-1]); + } + + //! beginning of adjacencies + Const_adjacency_iterator begin() const { + return this->ptr()->_m_data.begin(); + } + + //! end of adjacencies + Const_adjacency_iterator end() const { + return this->ptr()->_m_data.end(); + } + + //! The number of adjacency pairs + unsigned int size() const { + return this->ptr()->_m_data.size(); + } + + //! print (for testing only) + void print() const { + for(int i = 0; + i < static_cast(this->ptr()->_m_data.size()); + i++) { + std::cout << "(" << this->ptr()->_m_data[i].first << ", " + << this->ptr()->_m_data[i].second << ")" << std::endl; + } + } + + //! swap two instances + Self swap() const { + Self copy; + for(Const_adjacency_iterator it = begin(); it!=end(); it++) { + copy.ptr()->_m_data.push_back + (std::make_pair(it->second,it->first)); + } + return copy; + } + + //! output operator + void pretty_print(std::ostream& os) const { + for (int i = 0; + i < static_cast(this->ptr()->_m_data.size()); + i++) { + os << "<" << this->ptr()->_m_data[i].first + << "," << this->ptr()->_m_data[i].second << ">"; + if (i < static_cast(this->ptr()->_m_data.size()) - 1) { + os << ","; + } + } + } + +}; // end of Adjacencies data structure + +/*!\relates Adjacencies_3 + * \brief + * outputs Adjacencies_3 object to stream \c os + */ +inline +std::ostream& operator<<( + std::ostream& os, + const Adjacencies_3& adj) { + adj.pretty_print(os); + return os; +} + +} // namespace SoX + +#endif // SoX_GAPS_ADJACENCIES_3 diff --git a/Curved_kernel_via_analysis_2/include/CGAL/Arrangement_2l/Arr_p_dcel_info_overlay_traits.h b/Curved_kernel_via_analysis_2/include/CGAL/Arrangement_2l/Arr_p_dcel_info_overlay_traits.h new file mode 100644 index 00000000000..a5d5d03a88b --- /dev/null +++ b/Curved_kernel_via_analysis_2/include/CGAL/Arrangement_2l/Arr_p_dcel_info_overlay_traits.h @@ -0,0 +1,471 @@ +// ============================================================================ +// +// Copyright (c) 2001-2008 Max-Planck-Institut Saarbruecken (Germany). +// All rights reserved. +// +// This file is part of EXACUS (http://www.mpi-inf.mpg.de/projects/EXACUS/); +// you may redistribute it under the terms of the Q Public License version 1.0. +// See the file LICENSE.QPL distributed with EXACUS. +// +// Licensees holding a valid commercial license may use this file in +// accordance with the commercial license agreement provided with the software. +// +// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE +// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. +// +// ---------------------------------------------------------------------------- +// +// Library : SoX +// File : include/SoX/GAPS/Arr_p_dcel_info_overlay_traits.h +// SoX_release : $Name: $ +// Revision : $Revision$ +// Revision_date : $Date$ +// +// Author(s) : Eric Berberich +// +// ============================================================================ + +/*! \file SoX/GAPS/Arr_p_dcel_info_overlay_traits.h + * \brief definition of Arr_p_dcel_info_overlay_traits class template + */ + +#ifndef SoX_GAPS_ARR_P_DCEL_INFO_OVERLAY_TRAITS_H +#define SoX_GAPS_ARR_P_DCEL_INFO_OVERLAY_TRAITS_H 1 + +#include + +#include + +#include + +#include +#include + +namespace SoX { + +/*!\brief + * Model of CGAL::ArrangementOverlayTraits_2 for Arrangements + * attached with SoX::P_dcel_info<> data. + */ +template < class RestrictedCad_3 > +class Arr_p_dcel_info_overlay_traits : public +CGAL::Arr_default_overlay_traits< typename RestrictedCad_3::Rep > { + +public: + //! this instance's first template parameter + typedef RestrictedCad_3 Restricted_cad_3; + + //! type of arrangement + typedef typename Restricted_cad_3::Rep Arrangement_2; + + //! type of surface + typedef typename Restricted_cad_3::Surface_3 Surface_3; + + //! base type + //typedef CGAL::Arr_default_overlay_traits< Arrangement_2 > Base; + + typedef Restricted_cad_3 RS_3; + + //! type of accesor + typedef Restricted_cad_3_accessor< Restricted_cad_3 > Accessor; + + //! the class itself + typedef Arr_p_dcel_info_overlay_traits< Arrangement_2 > Self; + + typedef typename RS_3::Vertex_const_handle Vertex_const_handle_A; + typedef typename RS_3::Vertex_const_handle Vertex_const_handle_B; + typedef typename RS_3::Vertex_const_handle Vertex_const_handle_R; + + typedef typename RS_3::Halfedge_const_handle Halfedge_const_handle_A; + typedef typename RS_3::Halfedge_const_handle Halfedge_const_handle_B; + typedef typename RS_3::Halfedge_const_handle Halfedge_const_handle_R; + + typedef typename RS_3::Face_const_handle Face_const_handle_A; + typedef typename RS_3::Face_const_handle Face_const_handle_B; + typedef typename RS_3::Face_const_handle Face_const_handle_R; + + + /*!\brief + * Constructs a new overlay traits + */ + Arr_p_dcel_info_overlay_traits(const Surface_3& surface, + const Restricted_cad_3& cad, + SoX::Nk::Value_type type) : + _m_surface_mode(true), + _m_surface(surface), + _m_cad(cad), + _m_acc(cad), + _m_type(type), + _m_factors_of_same_curve(false) + { + + } + + + /*!\brief + * Constructs a new overlay traits where \c factors_of_same_curve defines + * whether the two parts originate from two factors of the same curve + */ + Arr_p_dcel_info_overlay_traits(const Restricted_cad_3& cad, + bool factors_of_same_curve) : + _m_surface_mode(false), + _m_cad(cad), + _m_acc(cad), + _m_factors_of_same_curve(factors_of_same_curve) { + } + +private: + + // for vertex/face + inline void merge(Nk nk1, Nk nk2, const Nk& nk, + const Dcel_feature& feature) const { + + if (this->_m_type == SoX::Nk::MULT) { + // mults! + if (feature == SoX::VERTEX) { + nk.set_mult(0); + } else if (feature == SoX::EDGE) { + nk.set_mult( + (nk1.mult() > 0 ? nk1.mult() : 0) + + (nk2.mult() > 0 ? nk2.mult() : 0) + ); + CGAL_assertion(nk.mult() > 0); + } + } else { + // copy already set values! + if (nk1.mult() >= 0) { + if (feature == SoX::VERTEX) { + nk.set_mult(0); + } else { + CGAL_assertion(feature == SoX::EDGE); + nk.set_mult(nk1.mult()); + } + } + if (nk1.n() >= -1) { + nk.set_n(nk1.n()); + } + if (nk1.k() >= -1) { + nk.set_k(nk1.k()); + } + } + } + +public: + + + /*! + * Create a vertex v that corresponds to the coinciding vertices v1 and v2. + */ + virtual void create_vertex (Vertex_const_handle_A v1, + Vertex_const_handle_B v2, + Vertex_const_handle_R v) + { + CGAL_precondition(v1->point() == v->point()); + CGAL_precondition(v2->point() == v->point()); + CGAL_assertion(v1->data()); + CGAL_assertion(v2->data()); + + if (_m_surface_mode) { + CGAL_precondition(!v->data()); + _m_cad._init_feature(_m_surface, v, SoX::VERTEX); + + v1->data()->_nk(_m_surface); + v2->data()->_nk(_m_surface); + v->data()->_nk(_m_surface); + + merge(v1->data()->_nk(_m_surface), + v2->data()->_nk(_m_surface), + v->data()->_nk(_m_surface), SoX::VERTEX); + CGAL_postcondition(v->data()); + + v->data()->_nk(_m_surface).set_feature1(SoX::VERTEX); + v->data()->_nk(_m_surface).set_feature2(SoX::VERTEX); + } else { + _m_acc.non_const_handle(v)->data() = + v1->data()->_overlay_with(*v2->data(), + _m_factors_of_same_curve); + } + } + + /*! + * Create a vertex v that mathces v1, which lies of the edge e2. + */ + virtual void create_vertex (Vertex_const_handle_A v1, + Halfedge_const_handle_B e2, + Vertex_const_handle_R v) + { + CGAL_precondition(v1->point() == v->point()); + CGAL_assertion(v1->data()); + CGAL_assertion(e2->data()); + + if (_m_surface_mode) { + CGAL_precondition(!v->data()); + _m_cad._init_feature(_m_surface, v, SoX::VERTEX); + merge(v1->data()->_nk(_m_surface), + e2->data()->_nk(_m_surface), + v->data()->_nk(_m_surface), SoX::VERTEX); + CGAL_postcondition(v->data()); + + v->data()->_nk(_m_surface).set_feature1(SoX::VERTEX); + v->data()->_nk(_m_surface).set_feature2(SoX::EDGE); + } else { + _m_acc.non_const_handle(v)->data() = + v1->data()->_overlay_with(*e2->data(), + _m_factors_of_same_curve); + } + } + + /*! + * Create a vertex v that mathces v1, contained in the face f2. + */ + virtual void create_vertex (Vertex_const_handle_A v1, + Face_const_handle_B f2, + Vertex_const_handle_R v) + { + CGAL_precondition(v1->point() == v->point()); + CGAL_assertion(v1->data()); + CGAL_assertion(f2->data()); + + if (_m_surface_mode) { + CGAL_precondition(!v->data()); + _m_cad._init_feature(_m_surface, v, SoX::VERTEX); + merge(v1->data()->_nk(_m_surface), + f2->data()->_nk(_m_surface), + v->data()->_nk(_m_surface), SoX::VERTEX); + CGAL_postcondition(v->data()); + + v->data()->_nk(_m_surface).set_feature1(SoX::VERTEX); + v->data()->_nk(_m_surface).set_feature2(SoX::FACE); + + } else { + _m_acc.non_const_handle(v)->data() = + v1->data()->_overlay_with(*f2->data(), + _m_factors_of_same_curve); + } + } + + /*! + * Create a vertex v that mathces v2, which lies of the edge e1. + */ + virtual void create_vertex (Halfedge_const_handle_A e1, + Vertex_const_handle_B v2, + Vertex_const_handle_R v) + { + CGAL_precondition(v2->point() == v->point()); + CGAL_assertion(e1->data()); + CGAL_assertion(v2->data()); + + if (_m_surface_mode) { + CGAL_precondition(!v->data()); + _m_cad._init_feature(_m_surface, v, SoX::VERTEX); + merge(e1->data()->_nk(_m_surface), + v2->data()->_nk(_m_surface), + v->data()->_nk(_m_surface), SoX::VERTEX); + CGAL_postcondition(v->data()); + + v->data()->_nk(_m_surface).set_feature1(SoX::EDGE); + v->data()->_nk(_m_surface).set_feature2(SoX::VERTEX); + } else { + _m_acc.non_const_handle(v)->data() = + v2->data()->_overlay_with(*e1->data(), + _m_factors_of_same_curve); + } + } + + /*! + * Create a vertex v that mathces v2, contained in the face f1. + */ + virtual void create_vertex (Face_const_handle_A f1, + Vertex_const_handle_B v2, + Vertex_const_handle_R v) + { + CGAL_precondition(v2->point() == v->point()); + CGAL_assertion(f1->data()); + CGAL_assertion(v2->data()); + + if (_m_surface_mode) { + CGAL_precondition(!v->data()); + _m_cad._init_feature(_m_surface, v, SoX::VERTEX); + merge(f1->data()->_nk(_m_surface), + v2->data()->_nk(_m_surface), + v->data()->_nk(_m_surface), SoX::VERTEX); + CGAL_postcondition(v->data()); + + v->data()->_nk(_m_surface).set_feature1(SoX::FACE); + v->data()->_nk(_m_surface).set_feature2(SoX::VERTEX); + } else { + _m_acc.non_const_handle(v)->data() = + v2->data()->_overlay_with(*f1->data(), + _m_factors_of_same_curve); + } + } + + /*! + * Create a vertex v that mathces the intersection of the edges e1 and e2. + */ + virtual void create_vertex (Halfedge_const_handle_A e1, + Halfedge_const_handle_B e2, + Vertex_const_handle_R v) + { + CGAL_assertion(e1->data()); + CGAL_assertion(e2->data()); + + if (_m_surface_mode) { + CGAL_precondition(!v->data()); + _m_cad._init_feature(_m_surface, v, SoX::VERTEX); + merge(e1->data()->_nk(_m_surface), + e2->data()->_nk(_m_surface), + v->data()->_nk(_m_surface), SoX::VERTEX); + CGAL_postcondition(v->data()); + + v->data()->_nk(_m_surface).set_feature1(SoX::EDGE); + v->data()->_nk(_m_surface).set_feature2(SoX::EDGE); + } else { + _m_acc.non_const_handle(v)->data() = + e1->data()->_overlay_with(*e2->data(), + _m_factors_of_same_curve); + } + } + + /*! + * Create an edge e that matches the overlap between e1 and e2. + */ + virtual void create_edge (Halfedge_const_handle_A e1, + Halfedge_const_handle_B e2, + Halfedge_const_handle_R e) + { + CGAL_assertion(e1->data()); + CGAL_assertion(e2->data()); + + if (_m_surface_mode) { + CGAL_precondition(!e->data()); + _m_cad._init_feature(_m_surface, e, SoX::EDGE); + merge(e1->data()->_nk(_m_surface), + e2->data()->_nk(_m_surface), + e->data()->_nk(_m_surface), SoX::EDGE); + CGAL_postcondition(e->data()); + + e->data()->_nk(_m_surface).set_feature1(SoX::EDGE); + e->data()->_nk(_m_surface).set_feature2(SoX::EDGE); + } else { + _m_acc.non_const_handle(e)->data() = + e1->data()->_overlay_with(*e2->data(), + _m_factors_of_same_curve); + } + CGAL_assertion(!e->twin()->data()); + _m_acc.non_const_handle(e->twin())->data() = *e->data(); + } + + /*! + * Create an edge e that matches the edge e1, contained in the face f2. + */ + virtual void create_edge (Halfedge_const_handle_A e1, + Face_const_handle_B f2, + Halfedge_const_handle_R e) + { + CGAL_assertion(e1->data()); + CGAL_assertion(f2->data()); + + if (_m_surface_mode) { + CGAL_precondition(!e->data()); + _m_cad._init_feature(_m_surface, e, SoX::EDGE); + merge(e1->data()->_nk(_m_surface), + f2->data()->_nk(_m_surface), + e->data()->_nk(_m_surface), SoX::EDGE); + CGAL_postcondition(e->data()); + + e->data()->_nk(_m_surface).set_feature1(SoX::EDGE); + e->data()->_nk(_m_surface).set_feature2(SoX::FACE); + } else { + _m_acc.non_const_handle(e)->data() = + e1->data()->_overlay_with(*f2->data(), + _m_factors_of_same_curve); + } + CGAL_assertion(!e->twin()->data()); + _m_acc.non_const_handle(e->twin())->data() = *e->data(); + } + + /*! + * Create an edge e that matches the edge e2, contained in the face f1. + */ + virtual void create_edge (Face_const_handle_A f1, + Halfedge_const_handle_B e2, + Halfedge_const_handle_R e) + { + CGAL_assertion(f1->data()); + CGAL_assertion(e2->data()); + + if (_m_surface_mode) { + CGAL_precondition(!e->data()); + _m_cad._init_feature(_m_surface, e, SoX::EDGE); + merge(f1->data()->_nk(_m_surface), + e2->data()->_nk(_m_surface), + e->data()->_nk(_m_surface), SoX::EDGE); + CGAL_postcondition(e->data()); + + e->data()->_nk(_m_surface).set_feature1(SoX::FACE); + e->data()->_nk(_m_surface).set_feature2(SoX::EDGE); + } else { + _m_acc.non_const_handle(e)->data() = + e2->data()->_overlay_with(*f1->data(), + _m_factors_of_same_curve); + } + CGAL_assertion(!e->twin()->data()); + _m_acc.non_const_handle(e->twin())->data() = *e->data(); + } + + /*! + * Create a face f that matches the overlapping region between f1 and f2. + */ + virtual void create_face (Face_const_handle_A f1, + Face_const_handle_B f2, + Face_const_handle_R f) + { + CGAL_assertion(f1->data()); + CGAL_assertion(f2->data()); + + if (_m_surface_mode) { + CGAL_precondition(!f->data()); + _m_cad._init_feature(_m_surface, f, SoX::FACE); + CGAL_assertion(f1->data()->_nk(_m_surface).mult() == -1); + CGAL_assertion(f2->data()->_nk(_m_surface).mult() == -1); + merge(f1->data()->_nk(_m_surface), + f2->data()->_nk(_m_surface), + f->data()->_nk(_m_surface), SoX::FACE); + CGAL_postcondition(f->data()); + + f->data()->_nk(_m_surface).set_feature1(SoX::FACE); + f->data()->_nk(_m_surface).set_feature2(SoX::FACE); + } else { + _m_acc.non_const_handle(f)->data() = + f1->data()->_overlay_with(*f2->data(), + _m_factors_of_same_curve); + } + } + +private: + // members + //! surface mode? + bool _m_surface_mode; + + //! surface + Surface_3 _m_surface; + + //! cad + Restricted_cad_3 _m_cad; + + //! accessor + Accessor _m_acc; + + //! overlay type + SoX::Nk::Value_type _m_type; + + //! factors of same curve + bool _m_factors_of_same_curve; +}; + + +} // namespace SoX + +#endif // SoX_GAPS_ARR_P_DCEL_INFO_OVERLAY_TRAITS_H +// EO diff --git a/Curved_kernel_via_analysis_2/include/CGAL/Arrangement_2l/Curve_3.h b/Curved_kernel_via_analysis_2/include/CGAL/Arrangement_2l/Curve_3.h new file mode 100644 index 00000000000..317cf64e702 --- /dev/null +++ b/Curved_kernel_via_analysis_2/include/CGAL/Arrangement_2l/Curve_3.h @@ -0,0 +1,725 @@ +// ============================================================================ +// +// Copyright (c) 2001-2008 Max-Planck-Institut Saarbruecken (Germany). +// All rights reserved. +// +// This file is part of EXACUS (http://www.mpi-inf.mpg.de/projects/EXACUS/); +// you may redistribute it under the terms of the Q Public License version 1.0. +// See the file LICENSE.QPL distributed with EXACUS. +// +// Licensees holding a valid commercial license may use this file in +// accordance with the commercial license agreement provided with the software. +// +// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE +// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. +// +// ---------------------------------------------------------------------------- +// +// Library : SoX +// File : include/SoX/GAPS/Curve_3.h +// SoX_release : $Name: $ +// Revision : $Revision$ +// Revision_date : $Date$ +// +// Author(s) : Eric Berberich +// +// ============================================================================ + +#ifndef SoX_GAPS_CURVE_3_H +#define SoX_GAPS_CURVE_3_H 1 + +/*!\file SoX/GAPS/Curve_3.h + * \brief + * definition of \c Curve_3<> + */ + +#include + +#include + +#include +#include +#include + +#include + +namespace SoX { + +namespace Intern { + +template < class SurfaceZAtXyIsolatorTraits > +class Curve_3_rep : + public SoX::Intern::Surface_pair_3_rep< SurfaceZAtXyIsolatorTraits > { + +public: + //! this instance's template parameter + typedef SurfaceZAtXyIsolatorTraits Surface_z_at_xy_isolator_traits; + + SoX_SURFACE_Z_AT_XY_ISOLATOR_TRAITS_SNAP_TYPEDEFS( + Surface_z_at_xy_isolator_traits + ); + + //! this instance itself + typedef Curve_3_rep< Surface_z_at_xy_isolator_traits > Self; + + //! the base classe + typedef SoX::Intern::Surface_pair_3_rep< Surface_z_at_xy_isolator_traits > + Base; + + //! type of restricted cad + typedef typename Base::Restricted_cad_3 Restricted_cad_3; + + //! type of creator + typedef typename Base::Creator Creator; + + //! type of creator + typedef typename Base::Overlayer Overlayer; + +public: + //!\name Constructors + //!@{ + + Curve_3_rep(const Surface_3& surface1, const Surface_3& surface2) : + Base(surface1, surface2) { + } + + //!@} + + // add special members?? + +}; + +} // namespace Intern + +template < +class SurfaceZAtXyIsolatorTraits, +class Rep_ = SoX::Intern::Curve_3_rep< SurfaceZAtXyIsolatorTraits > +> +class Curve_3 : + public Surface_pair_3 { + +public: + //! this instance's first template parameter + typedef SurfaceZAtXyIsolatorTraits Surface_z_at_xy_isolator_traits; + + //! this instance's second template parameter + typedef Rep_ Rep; + + SoX_SURFACE_Z_AT_XY_ISOLATOR_TRAITS_SNAP_TYPEDEFS( + Surface_z_at_xy_isolator_traits + ); + + //! type of Base + typedef SoX::Surface_pair_3< Surface_z_at_xy_isolator_traits, Rep > Base; + + //! type of restricted cad + typedef typename Base::Restricted_cad_3 Restricted_cad_3; + + //!\name Constructors + //!@{ + + Curve_3(const Surface_3& surface1, const Surface_3& surface2) : + Base(surface1, surface2) { + } + + //!@} + +private: // make public members inaccessible + //!\name Cads + //!@{ + + //! rscad of first silhouette + Restricted_cad_3 silhouette1() const { + return Base::silhouette1(); + } + + //! rscad of second silhouette + Restricted_cad_3 silhouette2() const { + return Base::silhouette1(); + } + + //! rscad of cut + Restricted_cad_3 cut() const { + return Base::cut(); + } + + //! rscad of first silhouette with cut + Restricted_cad_3 silhouettes_cut() const { + return Base::silhouettes_cut(); + } + + //!@} + +private: // make public members inaccessible + //!\name Curves + //!@{ + + /*!\brief + * returns silhouettes curves and points for a given \c surface in pair + */ + template < class CurveOutputIterator, class PointOutputIterator > + void silhouette_objects( + const Surface_3& surface, + CurveOutputIterator coi, PointOutputIterator poi) const { + return Base::silhouette_objects(surface, coi, poi); + + } + + /*!\brief + * returns projected intersection curves and points of the pair + */ + template < class CurveOutputIterator, class PointOutputIterator > + void cut_objects( + CurveOutputIterator coi, PointOutputIterator poi) const { + return Base::cut_objects(coi, poi); + } + + //!@} + +public: + + //!\name Space Curve Objects + //!@{ + /*!\brief + * returns spatial curves/points as combination of curve/sheet or + * point/sheet, where sheet is an int wrt to smaller-degree surface + */ + template < + class SurfaceArc_3, + class CurveOutputIterator, + class PointOutputIterator + > + void objects( + SurfaceArc_3 dummy, + Surface_3& support, + CurveOutputIterator coi, PointOutputIterator poi + ) const { + + typedef SurfaceArc_3 Surface_arc_3; + typedef typename Surface_arc_3::Surface_point_2l Surface_point_3; + + typedef typename Restricted_cad_3::Vertex_const_iterator + Vertex_const_iterator; + typedef typename Restricted_cad_3::Edge_const_iterator + Edge_const_iterator; + typedef typename Restricted_cad_3::Z_stack Z_stack; + + // Remark (1): + // The simple case is an edge with mult == 1, ask member for it + // -> Refine until only one z-cell of f remains (f is simpler than g!) + // Else: Multiple intersection can occur. + // Either use local gcd (point_on_curve for SR_i, i > 0!!!, as point + // already lies on SR_0, + // or use silhouette of g (costly?!?!) to determine isolator :-( + + Surface_3 other; + + // TODO select smaller surface out of the two by traits?!?! + if (this->surface2().f().degree() < this->surface1().f().degree()) { + support = this->surface2(); + other = this->surface1(); + } else { + support = this->surface1(); + other = this->surface2(); + } + + // TODO support is not allowed to be vertical! + CGAL_assertion(support.f().degree() > 0); + if (support.f().degree() != CGAL::total_degree(support.f())) { + // TODO what if support is not z-regular + // and what if other is not z-regular??! + std::cerr << "Warning: Support is not z-regular, " + << "we might encounter problems with vertical " + << "objects, or asympotes!" + << std::endl; + } + + Restricted_cad_3 cad = + (support.id() == this->surface1().id() ? + this->_silhouette1_cut() : + this->_silhouette2_cut()); + // we won't access z-stacks of cad directly, only for support! + + SoX::Intern::Refineable_interval_helper< Z_at_xy_isolator > + iv_helper; + + // iterate over all features of cut + for (Edge_const_iterator eit = cad.edges_begin(); + eit != cad.edges_end(); eit++) { + + if (cad.has_cut(eit)) { + + boost::optional< int > mult = + eit->data()->multiplicity_of_cut(support, other); + CGAL_assertion(mult); + + Point_2 point = cad.sample_point(eit); + + typename + Surface_z_at_xy_isolator_traits::Construct_isolator + construct_isolator; // TODO single instance + + Z_at_xy_isolator support_isolator = + construct_isolator(support, point, + cad.nk(eit, support), + (cad.has_silhouette(eit) ? + SoX::EDGE : SoX::FACE) + ); + + int n = support_isolator.number_of_real_roots(); + if (n == 0) { + continue; + } + + std::list< int > sheets; + + if (*mult == 1) { + + Z_at_xy_isolator other_isolator = + construct_isolator( + other, point, + cad.nk(cad.faces_begin(), support), + SoX::FACE + ); + std::cout << "RRs" + << support_isolator.number_of_real_roots() + << std::endl; + + + std::cout << "RRo" + << other_isolator.number_of_real_roots() + << std::endl; + + CGAL_assertion(other_isolator.number_of_real_roots() > 0); + + CGAL_assertion(point == other_isolator.traits().point()); + + std::list< std::pair< int, int > > overlaps; + + iv_helper.refined_overlaps( + support_isolator, other_isolator, + std::back_inserter(overlaps) + ); + + CGAL_assertion(!overlaps.empty()); + + std::pair< int, int > overlap = + iv_helper.unique_overlap( + support_isolator, other_isolator, + overlaps.begin(), overlaps.end() + ); + + sheets.push_back(overlap.first); + + } else { + + typename + Surface_z_at_xy_isolator_traits::Equal_z + equal_z; // TODO single instance + + typename + Surface_z_at_xy_isolator_traits::Point_on_curve_2 + point_on_curve; // TODO single instance + + // can start with 1 as it lies on k = 0 + int k = 1; + + // TODO what happens if surfaces are not z-regular??! + Polynomial_3 local_gcd = + equal_z.local_gcd(support, other, + support_isolator.traits(), k); + + // TODO make use of k + + // TODO be careful if a surface has a vertical line! + + Surface_3 gcd_surface = + Surface_3::surface_cache()(local_gcd); + + Polynomial_2 sil_local_gcd = gcd_surface.resultant_f_fz(); + + bool on_gcd_sil = point_on_curve(point, sil_local_gcd); + + Z_at_xy_isolator gcd_isolator = + construct_isolator( + local_gcd, point, + cad.nk(eit, support), + (on_gcd_sil ? SoX::EDGE : SoX::FACE) + ); + + // compute overlaps + std::vector< std::pair< int, int > > overlaps; + + iv_helper.refined_overlaps( + support_isolator, gcd_isolator, + std::back_inserter(overlaps) + ); + + // refine until a gcd-interval is fully include in + // support interval or has moved away + for (int k = 0; k < static_cast< int >(overlaps.size()); + k++) { + + while (true) { + + if (iv_helper.overlap_or_order( + gcd_isolator, overlaps[k].second, + support_isolator, overlaps[k].first + ) != CGAL::EQUAL) { + break; + } + // else + + if (iv_helper.is_included( + gcd_isolator, overlaps[k].second, + support_isolator, overlaps[k].first, + true + ) + ) { + sheets.push_back(overlaps[k].first); + break; + } + + // else + gcd_isolator.refine_interval(overlaps[k].second); + } + } + } + + for (std::list< int >::const_iterator sit = sheets.begin(); + sit != sheets.end(); sit++) { + + //std::cout << "cv: " << eit->curve() << std::endl; + + bool min_finite = + eit->curve().is_finite(CGAL::ARR_MIN_END); + bool max_finite = + eit->curve().is_finite(CGAL::ARR_MAX_END); + + bool src_at_min = + (eit->direction() == CGAL::ARR_LEFT_TO_RIGHT); + + Surface_point_3 min, max; + + int sheet = *sit; + int sheet_at_min = -2; + int sheet_at_max = -2; + bool z_inf_at_min = false; + bool z_inf_at_max = false; + + if (min_finite) { + Point_2 pmin = + eit->curve().curve_end(CGAL::ARR_MIN_END); + + std::pair< Z_stack, std::pair< int, int > > adjmin = + cad.surface_adjacency(eit, support, sheet, + (src_at_min ? + eit->source() : + eit->target())); + + sheet_at_min = adjmin.second.first; + z_inf_at_min = + (sheet_at_min == -1 || sheet_at_min == n); + if (!z_inf_at_min) { + min = Surface_point_3(pmin, support, sheet_at_min); + } + } + + if (max_finite) { + Point_2 pmax = + eit->curve().curve_end(CGAL::ARR_MAX_END); + + std::pair< Z_stack, std::pair< int, int > > adjmax = + cad.surface_adjacency(eit, support, sheet, + (src_at_min ? + eit->target() : + eit->source())); + + sheet_at_max = adjmax.second.first; + z_inf_at_max = + (sheet_at_max == -1 || sheet_at_max == n); + if (!z_inf_at_max) { + max = Surface_point_3(pmax, support, sheet_at_max); + } + } + + if (min_finite && max_finite) { + // arc is bounded + + if (z_inf_at_min && z_inf_at_max) { + // no z-asymptote + *coi++ = Surface_arc_3( + eit->curve(), min, max, support, + sheet, sheet_at_min, sheet_at_max + ); + } else if (!z_inf_at_min && !z_inf_at_max) { + // z-asymptotes at both sides + *coi++ = Surface_arc_3( + eit->curve(), + (sheet_at_min == -1 ? + CGAL::ARR_MIN_END : CGAL::ARR_MAX_END), + (sheet_at_max == -1 ? + CGAL::ARR_MIN_END : CGAL::ARR_MAX_END), + support, + sheet + ); + } else { + // z-asymptote at one-side + int sheet_at_point = sheet_at_min; + if (!z_inf_at_min) { + sheet_at_point = sheet_at_max; + } + *coi++ = Surface_arc_3( + eit->curve(), + (z_inf_at_min ? max : min), + (sheet_at_point == -1 ? + CGAL::ARR_MIN_END : CGAL::ARR_MAX_END), + support, + sheet, + (z_inf_at_min ? sheet_at_max : sheet_at_min) + ); + } + } else if (!min_finite && !max_finite) { + // planar arc is a branch + *coi++ = Surface_arc_3( + eit->curve(), + support, + sheet + ); + } else { + // arc is ray + if (!z_inf_at_min && !z_inf_at_max) { + // usual one + *coi++ = Surface_arc_3( + eit->curve(), + (min_finite ? min : max), + support, + sheet, + (min_finite ? sheet_at_min : sheet_at_max) + ); + + } else { + int sheet_at_point = sheet_at_min; + if (!z_inf_at_min) { + sheet_at_point = sheet_at_max; + } + // z-asymptote at finite-side + *coi++ = Surface_arc_3 ( + eit->curve(), + (sheet_at_point == -1 ? + CGAL::ARR_MIN_END : CGAL::ARR_MAX_END), + support, + sheet + ); + } + } + } + } + } + + // search only for isolated points in support and + // check only them with the help of local gcd! + + // iterate over all features of cut + for (Vertex_const_iterator vit = cad.vertices_begin(); + vit != cad.vertices_end(); vit++) { + if (vit->is_at_infinity()) { + continue; + } + if (cad.has_cut(vit)) { + + // check whether vertex is isolated or whether support + // has an isolated event here + + // actually we test whether each sheet k at vit is connected + // to something + + Point_2 point = cad.sample_point(vit); + + typename + Surface_z_at_xy_isolator_traits::Construct_isolator + construct_isolator; // TODO single instance + + Z_at_xy_isolator support_isolator = + // edge is ok, as we don't want to activate artifical + // x-interval, but other isolators (m-k..) + construct_isolator( + support, point, + cad.nk(vit, support), + (cad.has_silhouette(vit) ? + SoX::EDGE : SoX::FACE) + ); + + int n = support_isolator.number_of_real_roots(); + + if (cad.nk(vit, support).n() == -1) { // vertical line + // TODO sufficient condition? actually both surfaces + // need to be vertical, right? + // better criterion: both polynomials vanish? + + if (n == 0) { + // whole line + *coi++ = Surface_arc_3(point, support); + } else { + CGAL_assertion(n > 0); + Surface_point_3 min, max; + min = Surface_point_3(point, support, 0); + for (int sheet = 0; sheet < n; sheet++) { + if (sheet == 0) { + // ray from z=-oo + *coi++ = Surface_arc_3(min, + CGAL::ARR_MIN_END, + support); + } + if (sheet + 1 == n) { + // ray to z==oo + *coi++ = Surface_arc_3(min, + CGAL::ARR_MAX_END, + support); + break; + } + CGAL_assertion(sheet + 1 < n); + max = Surface_point_3(point, support, sheet + 1); + // bounded vertical arc + *coi++ = Surface_arc_3(min, max, support); + min = max; + } + } + return; // as no further isolated vertices will occur + } + + // else + + if (n == 0) { + continue; + } + + std::set< int > sheets; + + if (!vit->is_isolated()) { + for (int k = 0; k < n; k++) { + typename Restricted_cad_3:: + Halfedge_around_vertex_const_circulator circ, + start = + vit->incident_halfedges(); + circ = start; + bool isolated = true; + do { + ++circ; + + std::pair< Z_stack, std::pair< int, int > > + adjmin = + cad.surface_adjacency(vit, support, k, + circ); + + if (adjmin.second.first <= adjmin.second.second) { + isolated = false; + break; + } + + } while (circ != start); + + if (isolated) { + sheets.insert(k); + } + } + + if (sheets.empty()) { + continue; + } + } else { + for (int k = 0; k < n; k++) { + sheets.insert(k); + } + } + + typename + Surface_z_at_xy_isolator_traits::Equal_z + equal_z; // TODO single instance + + typename + Surface_z_at_xy_isolator_traits::Point_on_curve_2 + point_on_curve; // TODO single instance + + int k = 1; + + Polynomial_3 local_gcd = + equal_z.local_gcd(support, other, + support_isolator.traits(), k); + + // TODO make use of k + + // TODO be carefull if one of the surfaces has a vertical line + + Surface_3 gcd_surface = + Surface_3::surface_cache()(local_gcd); + + Polynomial_2 sil_local_gcd = gcd_surface.resultant_f_fz(); + + bool on_gcd_sil = point_on_curve(point, sil_local_gcd); + + Z_at_xy_isolator gcd_isolator = + construct_isolator( + local_gcd, point, + cad.nk(vit, support), + (on_gcd_sil ? SoX::EDGE : SoX::FACE) + ); + + // compute overlaps + std::vector< std::pair< int, int > > overlaps; + + iv_helper.refined_overlaps( + support_isolator, gcd_isolator, + std::back_inserter(overlaps) + ); + + // refine until a gcd-interval is fully include in + // support interval or has moved away + for (int k = 0; k < static_cast< int >(overlaps.size()); + k++) { + + if (sheets.find(overlaps[k].first) == sheets.end()) { + // the overlaps[k].first sheets is not isolated! + continue; + } + + // here we are only faces with sheet levels + // that form an isolated point of support :-) + + while (true) { + + if (iv_helper.overlap_or_order( + gcd_isolator, overlaps[k].second, + support_isolator, overlaps[k].first + ) != CGAL::EQUAL) { + break; + } + // else + + if (iv_helper.is_included( + gcd_isolator, overlaps[k].second, + support_isolator, overlaps[k].first, + true + ) + ) { + *poi++ = Surface_point_3(vit->point(),support, + overlaps[k].first); + break; + } + + // else + gcd_isolator.refine_interval(overlaps[k].second); + } + } + } + } + } + + //!@} + +}; // Curve_3 + +} // namespace SoX + +#endif // SoX_GAPS_CURVE_3_H +// EOF diff --git a/Curved_kernel_via_analysis_2/include/CGAL/Arrangement_2l/P_dcel_info.h b/Curved_kernel_via_analysis_2/include/CGAL/Arrangement_2l/P_dcel_info.h new file mode 100644 index 00000000000..b869ad476fc --- /dev/null +++ b/Curved_kernel_via_analysis_2/include/CGAL/Arrangement_2l/P_dcel_info.h @@ -0,0 +1,1725 @@ +// ============================================================================ +// +// Copyright (c) 2001-2008 Max-Planck-Institut Saarbruecken (Germany). +// All rights reserved. +// +// This file is part of EXACUS (http://www.mpi-inf.mpg.de/projects/EXACUS/); +// you may redistribute it under the terms of the Q Public License version 1.0. +// See the file LICENSE.QPL distributed with EXACUS. +// +// Licensees holding a valid commercial license may use this file in +// accordance with the commercial license agreement provided with the software. +// +// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE +// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. +// +// ---------------------------------------------------------------------------- +// +// Library : SoX +// File : include/SoX/GAPS/P_dcel_info.h +// SoX_release : $Name: $ +// Revision : $Revision$ +// Revision_date : $Date$ +// +// Author(s) : Eric Berberich +// +// ============================================================================ + +/*! \file SoX/GAPS/P_dcel_info.h + * \brief definition of P_dcel_info class template + */ + +#ifndef SoX_GAPS_P_DCEL_INFO_H +#define SoX_GAPS_P_DCEL_INFO_H 1 + +#include + +#include +#include +#include + +#include +#include + +#include +#include +#include + +#include + +#include +#include + +#include +#include +#include +#include +#include + +namespace SoX { + +namespace Intern { + +/*!\brief + * Representation class for P_dcel_info + */ +template < class SurfaceZAtXyIsolatorTraits > +class P_dcel_info_rep { + +public: + // types + //! this instance's first template parameter + typedef SurfaceZAtXyIsolatorTraits Surface_z_at_xy_isolator_traits; + + SoX_SURFACE_Z_AT_XY_ISOLATOR_TRAITS_SNAP_TYPEDEFS( + Surface_z_at_xy_isolator_traits + ); + + //! the class itself + typedef P_dcel_info_rep< Surface_z_at_xy_isolator_traits > Self; + + //! type of Less for surfaces + typedef typename Surface_3::Surface_less_than Surface_less_than; + + //! type of surface container + typedef std::set< Surface_3, Surface_less_than > Surface_container; + + //! type of nk map for silhouettes + typedef std::map< Surface_3, SoX::Nk, Surface_less_than > + Silhouettes_nk_map; + + //! type of Silhouette map + typedef std::map< Surface_3, int, Surface_less_than > Silhouettes_map; + + //! type of Surface_pair + typedef std::pair< Surface_3, Surface_3 > Surface_pair; + + //! type of Less of surface pair + typedef CGAL::Pair_lexicographical_less_than< Surface_3, Surface_3, + Surface_less_than, Surface_less_than > + Surface_pair_less; + + //! type of Cuts map + typedef std::map< Surface_pair, int, Surface_pair_less > Cuts_map; + + //! type of Dcel_data + typedef SoX::P_dcel_info< Surface_z_at_xy_isolator_traits > Dcel_data; + + //! type of Z_stack + typedef SoX::Z_stack< Surface_z_at_xy_isolator_traits, Dcel_data > + Z_stack; + + //! type of id + typedef std::ptrdiff_t Rs_id; + +private: + //! originating dcel-handles for single surfaces + typedef std::map< Surface_3, std::pair< SoX::Dcel_feature, CGAL::Object >, + Surface_less_than > Surface_dcel_handle_map; + + //! originating dcel-handles for single surfaces + typedef std::map< Surface_pair, + std::pair< SoX::Dcel_feature,CGAL::Object >, Surface_pair_less > + Surface_pair_dcel_handle_map; + + + //! Standard constructor + P_dcel_info_rep (const SoX::Dcel_feature feature, + const CGAL::Object dcel_handle) : + _m_rs_id(0), + _m_dcel_feature(feature), + _m_dcel_handle(dcel_handle) { + } + +private: + // data members + //! stored id of handle + mutable Rs_id _m_rs_id; + + //! stored feature + mutable SoX::Dcel_feature _m_dcel_feature; + + //! stored for this dcel-handle + mutable CGAL::Object _m_dcel_handle; + + //! stored point + mutable boost::optional< Point_2 > _m_point; + + //! stored z-stack + mutable boost::optional < Z_stack > _m_z_stack; + + //! list of surfaces + mutable Surface_container _m_surfaces; + + //! map for Nk objects for silhouettes + mutable Silhouettes_nk_map _m_silhouettes_nk; + + //! indication whether silhouette is stored + mutable bool _m_has_silhouette_curve; + + //! indication whether vertical surface exists + mutable bool _m_has_vertical_surface; + + //! indication whether vertical line exists + mutable bool _m_has_vertical_line; + + //! list of involved cut curves + mutable Cuts_map _m_cuts; + + //! stored boundary class in x + mutable boost::optional < CGAL::Arr_boundary_type > _m_boundary_type_in_x; + + //! stored boundary class in y + mutable boost::optional < CGAL::Arr_boundary_type > _m_boundary_type_in_y; + + //! stored data + mutable std::list< Z_stack > _m_z_stacks; + + // TODO _m_intermediate_stack_cells + + //! stores originating dcel handles for surface + mutable Surface_dcel_handle_map _m_surface_dcel_handle_map; + + //! stores originating dcel handles for surface pairs + mutable Surface_pair_dcel_handle_map _m_surface_pair_dcel_handle_map; + + //! friends + friend class P_dcel_info< Surface_z_at_xy_isolator_traits >; + friend class Restricted_cad_3< Surface_z_at_xy_isolator_traits>; +}; + +} // namespace Intern + + +/*!\brief + * Attach information about surface to planar dcel. + */ +template < class SurfaceZAtXyIsolatorTraits > +class P_dcel_info : public +::CGAL::Handle_with_policy< Intern::P_dcel_info_rep< SurfaceZAtXyIsolatorTraits > > { +public: + + //! this instance's first template parameter + typedef SurfaceZAtXyIsolatorTraits Surface_z_at_xy_isolator_traits; + + SoX_SURFACE_Z_AT_XY_ISOLATOR_TRAITS_SNAP_TYPEDEFS( + Surface_z_at_xy_isolator_traits + ); + + //! type of rep; + typedef Intern::P_dcel_info_rep< Surface_z_at_xy_isolator_traits > Rep; + + //! type of Base + typedef ::CGAL::Handle_with_policy< Rep > Base; + + //! the class itself + typedef P_dcel_info< Surface_z_at_xy_isolator_traits > Self; + + //! type of Z_stack + typedef typename Rep::Z_stack Z_stack; + + //! type of Rs_id + typedef typename Rep::Rs_id Rs_id; + +protected: + + //! type of Silhouettes map + typedef typename Rep::Silhouettes_map Silhouettes_map; + + //! type of Silhouettes nk map + typedef typename Rep::Silhouettes_nk_map Silhouettes_nk_map; + + //! type of Cuts map + typedef typename Rep::Cuts_map Cuts_map; + + //! type of Surface_pair + typedef typename Rep::Surface_pair Surface_pair; + + //! type of Surface_dcel_handle_map + typedef typename Rep::Surface_dcel_handle_map Surface_dcel_handle_map; + + //! type of Surface_pair_dcel_handle_map + typedef typename Rep::Surface_pair_dcel_handle_map + Surface_pair_dcel_handle_map; + + //! pair canonicalizer + struct Canonicalize_pair { + + //! first in pair in always less than second + Surface_pair operator()(Surface_pair pair) { + typename Rep::Surface_less_than less; + if (!less(pair.first, pair.second)) { + std::swap(pair.first, pair.second); + } + return pair; + } + }; + +private: + //! type of X_coordinate + typedef typename Curve_kernel_2::X_coordinate_1 X_coordinate_1; + + //! type of Boundary + typedef typename Curve_kernel_2::Boundary Boundary; + + //! type of Event info + typedef typename Curve_analysis_2::Status_line_1 Status_line_1; + + //! type of Polynomial 1 + typedef typename Polynomial_2::NT Polynomial_1; + +private: + //! Standard constructor + P_dcel_info(const SoX::Dcel_feature feature, + const CGAL::Object& dcel_handle) : + Base(Rep(feature, dcel_handle)) { + } + + //!\name Id of cads + //!@{ + + //! returns id of underlying cad + Rs_id _rs_id() const { + return this->ptr()->_m_rs_id; + } + + //! sets id of underlying cad + void _set_rs_id(Rs_id rs_id) const { + this->ptr()->_m_rs_id = rs_id; + } + + //!@} + + //!\name Dcel Handles + //!@{ + + //! the dcel-handle + CGAL::Object _dcel_handle() const { + return this->ptr()->_m_dcel_handle; + } + + //! sets the dcel-handle + void _set_dcel(const SoX::Dcel_feature feature, + const CGAL::Object& dcel_handle) const { + this->ptr()->_m_dcel_feature = feature; + this->ptr()->_m_dcel_handle = dcel_handle; + } + + //! return (if exiting) a pair of dcel-feature + handle for \c surface + boost::optional< std::pair< SoX::Dcel_feature, CGAL::Object > > + _dcel(const Surface_3& surface) const { + typedef + boost::optional< std::pair< SoX::Dcel_feature, CGAL::Object > > + Return_type; + typename Surface_dcel_handle_map::iterator it = + this->ptr()->_m_surface_dcel_handle_map.find(surface); + if (it != this->ptr()->_m_surface_dcel_handle_map.end()) { + Return_type tmp = it->second; + return tmp; + } + // else + return boost::none; + } + + + /*!\brief + * stores for given \c surface the dcel in originating rscad + */ + void _set_dcel(const Surface_3& surface, + const SoX::Dcel_feature feature, + CGAL::Object obj) const { + + CGAL_assertion_code(( + { + typedef Restricted_cad_3< Surface_z_at_xy_isolator_traits> + Restricted_cad_3; + switch (feature) { + case SoX::FACE: { + typename Restricted_cad_3::Face_const_iterator fit; + CGAL_assertion_code(bool check = ) + CGAL::assign (fit, obj); + CGAL_assertion(check); + break; + } + case SoX::EDGE: { + typename Restricted_cad_3::Halfedge_const_iterator hit; + CGAL_assertion_code(bool check = ) + CGAL::assign (hit, obj); + CGAL_assertion(check); + break; + } + case SoX::VERTEX: { + typename Restricted_cad_3::Vertex_const_iterator vit; + CGAL_assertion_code(bool check = ) + CGAL::assign (vit, obj); + CGAL_assertion(check); + break; + } + } + }) + ); + + typename Surface_dcel_handle_map::iterator it = + this->ptr()->_m_surface_dcel_handle_map.find(surface); + CGAL_precondition(it == + this->ptr()->_m_surface_dcel_handle_map.end()); + this->ptr()->_m_surface_dcel_handle_map.insert( + std::make_pair(surface,std::make_pair(feature,obj)) + ); + } + + //! return (if exiting) a pair of dcel-feature + handle for \c surfaces + boost::optional< std::pair< SoX::Dcel_feature, CGAL::Object > > + _dcel(const Surface_3& surface1, const Surface_3& surface2) const { + typedef + boost::optional< std::pair< SoX::Dcel_feature, CGAL::Object > > + Return_type; + static Canonicalize_pair canonicalize; + Surface_pair pair = canonicalize(std::make_pair(surface1, surface2)); + typename Surface_pair_dcel_handle_map::iterator it = + this->ptr()->_m_surface_pair_dcel_handle_map.find(pair); + if (it != this->ptr()->_m_surface_pair_dcel_handle_map.end()) { + Return_type tmp = it->second; + return tmp; + } + // else + return boost::none; + } + + /*!\brief + * stores for given \c surface the dcel in originating rscad + */ + void _set_dcel(const Surface_3& surface1, + const Surface_3& surface2, + const SoX::Dcel_feature feature, + CGAL::Object obj) const { + static Canonicalize_pair canonicalize; + Surface_pair pair = canonicalize(std::make_pair(surface1, surface2)); + typename Surface_pair_dcel_handle_map::iterator it = + this->ptr()->_m_surface_pair_dcel_handle_map.find(pair); + CGAL_precondition( + it == this->ptr()->_m_surface_pair_dcel_handle_map.end() + ); + this->ptr()->_m_surface_pair_dcel_handle_map.insert( + std::make_pair(pair,std::make_pair(feature,obj)) + ); + } + + //!@} + + //!\name Adding + //!@{ + +private: + /*! adds surface + */ + // TODO remove with adding init_nk(s1,s2,nk) + void _add_surface(const Surface_3& surface) const { + this->ptr()->_m_surfaces.insert(surface); + } + + + //! initialized nk-object for surface + void _init_nk(const Surface_3& surface) const { + this->ptr()->_m_surfaces.insert(surface); + + SoX::Nk nk; + + typename Silhouettes_nk_map::iterator it = + this->ptr()->_m_silhouettes_nk.find(surface); + if (it == this->ptr()->_m_silhouettes_nk.end()) { + this->ptr()->_m_silhouettes_nk.insert( + it, std::make_pair(surface, nk) + ); + } + } + + /*! finalizes internal structures + */ + void _finalize_nk(const Surface_3& surface) const { + + SoX::Nk nk = _nk(surface); + +#if !NDEBUG + std::cout << "Finalizing Nk for " << this->ptr()->_m_dcel_feature + << " with " << nk << " ... " << std::flush; +#endif + + if (nk.mult() != -1) { + if (!this->ptr()->_m_has_silhouette_curve) { + this->ptr()->_m_has_silhouette_curve = true; + } + if (nk.n() == -1) { + if (this->ptr()->_m_dcel_feature == SoX::VERTEX) { + // TODO be carefull if point is at infinity + if (!this->ptr()->_m_has_vertical_line) { + this->ptr()->_m_has_vertical_line = true; + } + } else if (this->ptr()->_m_dcel_feature == SoX::EDGE) { + if (!this->ptr()->_m_has_vertical_surface) { + this->ptr()->_m_has_vertical_surface = true; + } + } + } + } +#if !NDEBUG + std::cout << "done." << std::endl; +#endif + } + + /*!\brief + * Add cut of \c surface1 and \c surface2 with \c multiplicity + */ + void _add_cut(const Surface_3& surface1, + const Surface_3& surface2, + int multiplicity, + bool factors_of_same_curve) const { + static Canonicalize_pair canonicalize; + Surface_pair pair = canonicalize(std::make_pair(surface1, surface2)); + typename Cuts_map::iterator it = + this->ptr()->_m_cuts.find(pair); + if (it == this->ptr()->_m_cuts.end()) { + this->ptr()->_m_cuts.insert( + it, std::make_pair(pair, multiplicity) + ); + } else { + if (factors_of_same_curve) { + it->second += multiplicity; + } + } + } + + //!@} + +private: + //!\name Query members + //!@{ + + /*!\brief + * returns \c true if carries positive number of stored silhouette curves + */ + bool _has_silhouette() const { + return (this->ptr()->_m_has_silhouette_curve); + } + + /*!\brief + * returns \c true if carries vertical silhouettes + */ + bool _has_vertical_surface() const { + return (this->ptr()->_m_has_vertical_surface); + } + + /*!\brief + * returns \c true if carries positive number of stored cut curves + */ + bool _has_cut() const { + return (!this->ptr()->_m_cuts.empty()); + } + + /*!\brief + * returns \c true if handle is silhouette of \c surface + */ + bool _has_silhouette(const Surface_3& surface) const { + if (!this->ptr()->_m_has_silhouette_curve) { + return false; + } + // else + typename Silhouettes_nk_map::iterator it = + this->ptr()->_m_silhouettes_nk.find(surface); + if (it == this->ptr()->_m_silhouettes_nk.end()) { + return false; + } + // else + return it->second.mult() != -1; + } + + /*! return nk-instance of given surface + */ + const SoX::Nk& _nk(const Surface_3& surface) const { + + typename Silhouettes_nk_map::iterator it = + this->ptr()->_m_silhouettes_nk.find(surface); + CGAL_assertion(it != this->ptr()->_m_silhouettes_nk.end()); + return (it->second); + } + + /*!\brief + * returns \c true of projected curve belongs to vertical part of + * \c surface + */ + // TODO new name with n + bool _is_surface_vertical( + const Surface_3& surface + ) const { + typename Silhouettes_nk_map::iterator it = + this->ptr()->_m_silhouettes_nk.find(surface); + CGAL_assertion(it != this->ptr()->_m_silhouettes_nk.end()); + return (it->second.n() == -1); + } + + /*!\brief + * returns whether supports vertical line + */ + // TODO new name with n + bool _supports_vertical_line(const Surface_3& surface) const { + typename Silhouettes_nk_map::iterator it = + this->ptr()->_m_silhouettes_nk.find(surface); + CGAL_assertion(it != this->ptr()->_m_silhouettes_nk.end()); + return (it->second.n() == -1); + } + +private: + /*!\brief + * returns \c true if carries positive number of stored cut curves + */ + bool _has_cut(const Surface_3& surface1, + const Surface_3& surface2) const { + static Canonicalize_pair canonicalize; + Surface_pair pair = canonicalize(std::make_pair(surface1, surface2)); + typename Cuts_map::iterator it = + this->ptr()->_m_cuts.find(pair); + return (it != this->ptr()->_m_cuts.end()); + } + +public: + + /*!\brief + * returns multiplicity of silhouette curve of \c surface + * + * \pre Is info for an edge. + */ + boost::optional< int > multiplicity_of_silhouette( + const Surface_3& surface + ) const { + CGAL_precondition(this->ptr()->_m_dcel_feature == SoX::EDGE); + typename Silhouettes_nk_map::iterator it = + this->ptr()->_m_silhouettes_nk.find(surface); + if (it == this->ptr()->_m_silhouettes_nk.end()) { + return boost::none; + } + // else + return boost::optional< int >(it->second.mult()); + } + + /*!\brief + * returns multiplicity of cut curve of \c surface1 and + * \c surface2 + * + * \pre Is info for an edg.e + */ + boost::optional< int > multiplicity_of_cut( + const Surface_3& surface1, + const Surface_3& surface2 + ) const { + CGAL_precondition(this->ptr()->_m_dcel_feature == SoX::EDGE); + static Canonicalize_pair canonicalize; + Surface_pair pair = canonicalize(std::make_pair(surface1, surface2)); + typename Cuts_map::iterator it = + this->ptr()->_m_cuts.find(pair); + if (it == this->ptr()->_m_cuts.end()) { + return boost::none; + } + // else + return boost::optional< int >(it->second); + } + + //!@} + + // TOOD make private and accessible only through RS_3 + //!\name Boundary classes + //!@{ + + /*!\brief + * returns known boundary class for object + */ + CGAL::Arr_boundary_type boundary_type_in_x() const { + if (this->ptr()->_m_boundary_type_in_x) { + return *this->ptr()->_m_boundary_type_in_x; + } + // else + return CGAL::ARR_NUMBER_OF_BOUNDARY_TYPES; + } + + /*!\brief + * returns known boundary class for object + */ + CGAL::Arr_boundary_type boundary_type_in_y() const { + if (this->ptr()->_m_boundary_type_in_y) { + return *this->ptr()->_m_boundary_type_in_y; + } + // else + return CGAL::ARR_NUMBER_OF_BOUNDARY_TYPES; + } + + /*!\brief + * Sets boundary class for data to given value \c cl + */ + void set_boundary_type_in_x(CGAL::Arr_boundary_type cl) const { + this->ptr()->_m_boundary_type_in_x = cl; + } + + /*!\brief + * Sets boundary class for data to given value \c cl + */ + void set_boundary_type_in_y(CGAL::Arr_boundary_type cl) const { + this->ptr()->_m_boundary_type_in_y = cl; + } + + //!@} + + +private: + // TASK move the following members to RS_3?; make static? make "global" + + //!\name Stored point + //!@{ + + //! returns the sample point for a vertex handle + template < class VertexHandle > + Point_2 _sample_point_for_vertex_handle(VertexHandle vh) const { + CGAL_precondition(!vh->is_at_infinity()); + + typedef Restricted_cad_3< Surface_z_at_xy_isolator_traits> + Restricted_cad_3; + +#if 0 && SoX_Z_STACK_Z_TIMERS + + bool running = z_stack_time[0].is_running(); + if (!running) { + z_stack_time[0].start(); + } +#endif + + Restricted_cad_3 cad = + Restricted_cad_3::cad_cache()( + vh->data()->ptr()->_m_rs_id + ); + + // check existing stacks + boost::optional< Z_stack > helper; + + for (typename std::list< Z_stack >::iterator sit = + this->ptr()->_m_z_stacks.begin(); + sit != this->ptr()->_m_z_stacks.end(); + sit++) { + + if (cad._point_on_dcel_handle(sit->point(), vh)) { + // if helper is not set or + // number of z_cell in current z_stack is greater + // than number of z_stack cells in current helper + if (!helper || + (sit->number_of_z_cells() > + helper->number_of_z_cells())) { + // store new values. + helper = *sit; + } + } // TASK implement assert + } + + // to obtain a final candidate + this->ptr()->_m_z_stacks.clear(); + if (helper) { + this->ptr()->_m_z_stacks.push_front(*helper); + } + + if (!this->ptr()->_m_point) { + CGAL_assertion(!this->ptr()->_m_z_stack); + + Point_2 pt; + // either helper gives point + if (helper) { + pt = helper->point(); + // _m_z_stacks store also z-stack candidate + } else { + // or we compute it + pt = vh->point(); + } + + _set_point(pt); + } else { + // second case: point is already set + if (this->ptr()->_m_z_stack) { + // and also z-stack + if (helper) { + // -> delete candidate + this->ptr()->_m_z_stacks.clear(); + } + } else { + if (helper) { + // reset point as helper is set and z_stack is unset + Point_2 pt = helper->point(); + _set_point(pt); + } + } + } + +#if 0 && SoX_Z_STACK_Z_TIMERS + if (!running) { + z_stack_time[0].stop(); + } +#endif + CGAL_postcondition(this->ptr()->_m_point); + CGAL_postcondition( + cad._point_on_dcel_handle(*this->ptr()->_m_point, vh) + ); + + // return stored z_stack + return *(this->ptr()->_m_point); + } + + //! returns the sample point for a halfedge handle + template < class HalfedgeHandle > + Point_2 _sample_point_for_halfedge_handle(HalfedgeHandle heh) const { + + typedef Restricted_cad_3< Surface_z_at_xy_isolator_traits> + Restricted_cad_3; + +#if 0 && SoX_Z_STACK_Z_TIMERS + bool running = z_stack_time[1].is_running(); + if (!running) { + z_stack_time[1].start(); + } +#endif + Restricted_cad_3 cad = + Restricted_cad_3::cad_cache()( + heh->data()->ptr()->_m_rs_id + ); + + // check existing stacks + boost::optional< Z_stack > helper; + + for (typename std::list< Z_stack >::iterator sit = + this->ptr()->_m_z_stacks.begin(); + sit != this->ptr()->_m_z_stacks.end(); + sit++) { + if (cad._point_on_dcel_handle(sit->point(), heh)) { + // if helper is not set or + // number of z_cell in current z_stack is greater + // than number of z_stack cells in current helper + if (!helper || + (sit->number_of_z_cells() > + helper->number_of_z_cells())) { + // store new values. + helper = *sit; + } + } // TASK implement assert + } + + // to obtain a final candidate + this->ptr()->_m_z_stacks.clear(); + if (helper) { + this->ptr()->_m_z_stacks.push_front(*helper); + } + + + if (!this->ptr()->_m_point) { + CGAL_assertion(!this->ptr()->_m_z_stack); + + Point_2 pt; + // either helper gives point + if (helper) { + pt = helper->point(); + // _m_z_stacks store also z-stack candidate + } else { + // or we compute it + pt = Restricted_cad_3::_point_in_interior(heh); + } + + _set_point(pt); + } else { + // second case: point is already set + if (this->ptr()->_m_z_stack) { + // and also z-stack + if (helper) { + // -> delete candidate + this->ptr()->_m_z_stacks.clear(); + } + } else { + if (helper) { + // reset point as helper is set and z_stack is unset + Point_2 pt = helper->point(); + _set_point(pt); + } + } + } + +#if 0 && SoX_Z_STACK_Z_TIMERS + if (!running) { + z_stack_time[1].stop(); + } +#endif + CGAL_postcondition(this->ptr()->_m_point); + CGAL_postcondition( + cad._point_on_dcel_handle(*this->ptr()->_m_point, heh) + ); + // return stored z_stack + return *(this->ptr()->_m_point); + } + + //! returns sample point for face handle + template < class FaceHandle > + Point_2 _sample_point_for_face_handle(FaceHandle fh) const { + + typedef Restricted_cad_3< Surface_z_at_xy_isolator_traits> + Restricted_cad_3; + +#if 0 && SoX_Z_STACK_Z_TIMERS + bool running = z_stack_time[2].is_running(); + if (!running) { + z_stack_time[2].start(); + } +#endif + + Restricted_cad_3 cad = + Restricted_cad_3::cad_cache()( + fh->data()->ptr()->_m_rs_id + ); + + // check existing stacks + boost::optional< Z_stack > helper; + + for (typename std::list< Z_stack >::iterator sit = + this->ptr()->_m_z_stacks.begin(); + sit != this->ptr()->_m_z_stacks.end(); + sit++) { + + if (cad._point_on_dcel_handle(sit->point(), fh)) { + // if helper is not set or + // number of z_cell in current z_stack is greater + // than number of z_stack cells in current helper + if (!helper || + (sit->number_of_z_cells() > + helper->number_of_z_cells())) { + // store new values. + helper = *sit; + } + } // TASK implement assert + } + + // to obtain a final candidate + this->ptr()->_m_z_stacks.clear(); + if (helper) { + this->ptr()->_m_z_stacks.push_front(*helper); + } + + + if (!this->ptr()->_m_point) { + CGAL_assertion(!this->ptr()->_m_z_stack); + + Point_2 pt; + // either helper gives point + if (helper) { + pt = helper->point(); + // _m_z_stacks store also z-stack candidate + } else { + // or we compute it + pt = _rational_point_in_face(fh, cad); + } + + _set_point(pt); + } else { + // second case: point is already set + if (this->ptr()->_m_z_stack) { + // and also z-stack + if (helper) { + // -> delete candidate + this->ptr()->_m_z_stacks.clear(); + } + } else { + if (helper) { + // reset point as helper is set and z_stack is unset + Point_2 pt = helper->point(); + _set_point(pt); + } + } + } +#if 0 && SoX_Z_STACK_Z_TIMERS + if (!running) { + z_stack_time[2].stop(); + } +#endif + CGAL_postcondition(this->ptr()->_m_point); + CGAL_postcondition( + cad._point_on_dcel_handle(*this->ptr()->_m_point, fh) + ); + + // return stored z_stack + return *(this->ptr()->_m_point); + } + + //!@} + +private: + //!\name Stored z-stack + //!@{ + + /*!\brief + * returns \c true if z-stack is known + */ + bool _knows_z_stack() const { + return this->ptr()->_m_z_stack; + } + + /*!\brief + * compute Z_stack for a given Vertex_handle + */ + template < class VertexHandle > + Z_stack _z_stack_for_vertex_handle(VertexHandle vh) const { + CGAL_precondition(!vh->is_at_infinity()); +#if 0 && SoX_Z_STACK_Z_TIMERS + bool running = z_stack_time[0].is_running(); + if (!running) { + z_stack_time[0].start(); + } +#endif + if (!this->ptr()->_m_z_stack) { + if (!this->ptr()->_m_point) { + _sample_point_for_vertex_handle(vh); + } + + // store z_stack + _set_z_stack(SoX::VERTEX); + } + + // delete all old z_stacks + this->ptr()->_m_z_stacks.clear(); + +#if 0 && SoX_Z_STACK_Z_TIMERS + if (!running) { + z_stack_time[0].stop(); + } +#endif + + // return stored z_stack + return *(this->ptr()->_m_z_stack); + } + + /*!\brief + * compute Z_stack for a given Halfedge_handle + */ + template < class HalfedgeHandle > + Z_stack _z_stack_for_halfedge_handle(HalfedgeHandle heh) const { + + typedef Restricted_cad_3< Surface_z_at_xy_isolator_traits> + Restricted_cad_3; + +#if 0 && SoX_Z_STACK_Z_TIMERS + bool running = z_stack_time[1].is_running(); + if (!running) { + z_stack_time[1].start(); + } +#endif + if (!this->ptr()->_m_z_stack) { + if (!this->ptr()->_m_point) { + _sample_point_for_halfedge_handle(heh); + } + + // store z_stack + _set_z_stack(SoX::EDGE); + } + + // delete all old z_stacks + this->ptr()->_m_z_stacks.clear(); + +#if 0 && SoX_Z_STACK_Z_TIMERS + if (!running) { + z_stack_time[1].stop(); + } +#endif + + // return stored z_stack + return *(this->ptr()->_m_z_stack); + } + + /*!\brief + * compute Z_stack for a given Face_handle + */ + template < class FaceHandle > + Z_stack _z_stack_for_face_handle(FaceHandle fh) const { + + typedef Restricted_cad_3< Surface_z_at_xy_isolator_traits> + Restricted_cad_3; + +#if 0 && SoX_Z_STACK_Z_TIMERS + bool running = z_stack_time[2].is_running(); + if (!running) { + z_stack_time[2].start(); + } +#endif + if (!this->ptr()->_m_z_stack) { + if (!this->ptr()->_m_point) { + _sample_point_for_face_handle(fh); + } + + // store z_stack + _set_z_stack(SoX::FACE); + } + + // delete all old z_stacks + this->ptr()->_m_z_stacks.clear(); + +#if 0 && SoX_Z_STACK_Z_TIMERS + if (!running) { + z_stack_time[2].stop(); + } +#endif + + // return stored z_stack + return *(this->ptr()->_m_z_stack); + } + + //! returns for given \c surface and \cad the store z_stack in the sil-cad + std::pair< Z_stack, SoX::Dcel_feature > _z_stack_of_surface( + const Surface_3& surface + ) const { + typedef Restricted_cad_3< Surface_z_at_xy_isolator_traits> + Restricted_cad_3; + + typename Surface_dcel_handle_map::iterator it = + this->ptr()->_m_surface_dcel_handle_map.find(surface); + CGAL_precondition(it != this->ptr()->_m_surface_dcel_handle_map.end()); + + if (it->second.first == SoX::FACE) { + typename Restricted_cad_3::Face_const_iterator fit; + CGAL_assertion_code(bool check = ) + CGAL::assign (fit, it->second.second); + CGAL_assertion(check); + return std::make_pair( + fit->data()->_z_stack_for_face_handle(fit), + SoX::FACE + ); + } else if (it->second.first == SoX::EDGE) { + typename Restricted_cad_3::Halfedge_const_iterator hit; + CGAL_assertion_code(bool check = ) + CGAL::assign (hit, it->second.second); + CGAL_assertion(check); + return std::make_pair( + hit->data()->_z_stack_for_halfedge_handle(hit), + SoX::EDGE + ); + } + // else + CGAL_assertion(it->second.first == SoX::VERTEX); + typename Restricted_cad_3::Vertex_const_iterator vit; + CGAL_assertion_code(bool check =) + CGAL::assign (vit, it->second.second); + CGAL_assertion(check); + return std::make_pair(vit->data()->_z_stack_for_vertex_handle(vit), + SoX::VERTEX); + } + + /*!\brief + * scans a connected component of the given ccb indicated by \c start + */ + template < class HalfedgeHandle > + static + std::pair< HalfedgeHandle, HalfedgeHandle > + _scan_ccb(HalfedgeHandle start) { + + HalfedgeHandle invalid_he; + HalfedgeHandle he; + HalfedgeHandle he_v; + + HalfedgeHandle curr = start; + do { + ++curr; + if (!curr->is_fictitious()) { + //std::cout << "curve: " << curr->curve() << std::endl; + //std::cout << "magickey" << std::endl; + if (curr->curve().is_vertical()) { + if (he_v == invalid_he) { + he_v = curr; + } + } else { + if (he == invalid_he) { + he = curr; + } + } + } + if (he != invalid_he) { + break; + } + } while (curr != start); + + // finally return + return std::make_pair(he, he_v); + } + + /*!\brief + * computes rational point in given face + */ + template < class FaceHandle, class RestrictedCad_3> + Point_2 _rational_point_in_face(FaceHandle fh, + RestrictedCad_3& cad) const { + CGAL_precondition(fh != FaceHandle()); + + typedef RestrictedCad_3 Restricted_cad_3; + + typedef typename Restricted_cad_3::Ccb_halfedge_const_circulator + Ccb_halfedge_const_circulator; + typedef typename Restricted_cad_3::Outer_ccb_const_iterator + Outer_ccb_const_iterator; + typedef typename Restricted_cad_3::Inner_ccb_const_iterator + Inner_ccb_const_iterator; + + Ccb_halfedge_const_circulator invalid_he; + + std::pair< Ccb_halfedge_const_circulator, + Ccb_halfedge_const_circulator > ocbpair, icbpair; + + Ccb_halfedge_const_circulator he, he_v; + + // select halfedge on boundary of face + // avoid vertical halfedges + for (Outer_ccb_const_iterator ocb = fh->outer_ccbs_begin(); + ocb != fh->outer_ccbs_end(); + ocb++) { + + ocbpair = _scan_ccb(*ocb); + + if (ocbpair.first != invalid_he && he == invalid_he) { + he = ocbpair.first; + break; + } + if (ocbpair.second != invalid_he && he_v == invalid_he) { + he_v = ocbpair.second; + } + } + + if (he == invalid_he) { + // search on inner ccbs + for (Inner_ccb_const_iterator icb = fh->inner_ccbs_begin(); + icb != fh->inner_ccbs_end(); + icb++) { + + icbpair = _scan_ccb(*icb); + + if (icbpair.first != invalid_he && he == invalid_he) { + he = icbpair.first; + break; + } + if (icbpair.second != invalid_he && he_v == invalid_he) { + he_v = icbpair.second; + } + } + + } + + if (he == invalid_he) { + // try whether vertical one has been found + if (he_v != invalid_he) { + he = he_v; + } + } + // he is either invalid, or we've found a normal edge or a vertical one + + // helping variables + int type = -1; + Boundary xr(0), y0(0); + X_coordinate_1 x0(xr); + int count = 0; + int arc = -1; + Status_line_1 sl; + + // type init phase + if (he == invalid_he) { + type = 0; + } else if (!he->curve().is_vertical()) { + type = 1; + Boundary xr1 = he->curve().boundary_in_x_range_interior(); + x0 = X_coordinate_1(xr1); + + arc = he->curve().arcno(); + sl = he->curve().curve().status_line_at_exact_x(x0); + } else { + CGAL_assertion(he->curve().is_vertical()); + type = 2; + } + + + typename Curve_kernel_2::Lower_boundary_y_2 lower_boundary_y = + Arrangement_traits_2::instance().kernel(). + lower_boundary_y_2_object(); + typename Curve_kernel_2::Upper_boundary_y_2 upper_boundary_y = + Arrangement_traits_2::instance().kernel(). + upper_boundary_y_2_object(); + typename Curve_kernel_2::Refine_y_2 refine_y = + Arrangement_traits_2::instance().kernel(). + refine_y_2_object(); + + while (true) { + // desired object + Point_2 pt; + + Point_2 ptv; + + // compute value phase + switch (type) { + case 0: { + x0 = X_coordinate_1(xr); + y0 = xr * xr; + break; + } + case 1: + // we found a normal halfedge + // (needs direction) + y0 = (he->direction() == CGAL::ARR_LEFT_TO_RIGHT ? + upper_boundary_y(sl.algebraic_real_2(arc)) : + lower_boundary_y(sl.algebraic_real_2(arc)) + ); + break; + case 2: { + ++count; + // in this case we could only find a + // vertical halfedge for the face + // that has to go from y=-oo to y =+00 + CGAL_assertion(he->curve().is_vertical()); + X_coordinate_1 x1 = he->curve().x(); + Boundary xt2 = (he->direction() == CGAL::ARR_LEFT_TO_RIGHT ? + x1.low() : x1.high()); + if (x1.is_rational()) { + xr = xt2 + + (he->direction() == CGAL::ARR_LEFT_TO_RIGHT ? + Boundary(-1) : Boundary(1))/ + CGAL::ipower(Boundary(2),count); + } else { + // TODO use ak + xr = x1.rational_between(X_coordinate_1(xt2)); + } + x0 = X_coordinate_1(xr); + + ptv = Restricted_cad_3::_point_in_interior(he); + + break; + } + default: + CGAL_error_msg("Only 0,1,2 are allowed as types"); + } + + if (type == 2) { + pt = Point_2(X_coordinate_1(x0),ptv.curve(),0); + } else { + // construct point phase + pt = Restricted_cad_3::_construct_point_with_rational_y( + x0, y0 + ); + } + + // check phase + // test whether point lies in the face with point location + if (cad._point_on_dcel_handle(pt, fh)) { + // found face matches given one, + // i.e., constructed point lies within correct face + return pt; + } + + // refine phase + switch (type) { + case 0: + xr += Boundary(1); + break; + case 1: + CGAL_assertion(arc != -1); + // refine y-coordinate + refine_y(sl.algebraic_real_2(arc)); + break; + case 2: + CGAL_assertion(arc == -1); + // refine x-coordinate + x0.strong_refine(xr); + break; + } + } + } + +private: + + /*!\brief + * sets point for z-stack construction + */ + void _set_point(const Point_2& pt) const { + CGAL_precondition(!this->ptr()->_m_point); + this->ptr()->_m_point = pt; + } + + /*!\brief + * sets z_stack for object + */ + void _set_z_stack(SoX::Dcel_feature feature) const { + CGAL_precondition(this->ptr()->_m_point); + CGAL_precondition(!this->ptr()->_m_z_stack); + + Point_2 pt = *this->ptr()->_m_point; + +#if !NDEBUG + std::cout << "feature: " << feature << std::endl; + std::cout << "Computing z-stack for " << pt << " ... " << std::flush; +#endif + + // tasks: + // - check whether all surface in stored z_stack exist in _m_surfaces + if (static_cast< int >(this->ptr()->_m_z_stacks.size()) > 1) { + CGAL_assertion( + static_cast< int >(this->ptr()->_m_z_stacks.size()) == 2 + ); + CGAL_assertion(pt == this->ptr()->_m_z_stacks.begin()->point()); + CGAL_assertion( + pt == (++this->ptr()->_m_z_stacks.begin())->point() + ); + Z_stack z_stack( + this->ptr()->_m_z_stacks.begin()->_merge( + *(this->ptr()->_m_z_stacks.begin()++), + this + ), + feature, this + ); + + this->ptr()->_m_z_stack = z_stack; + } else { + + CGAL_assertion_code( + if (static_cast< int >(this->ptr()->_m_z_stacks.size()) + == 1) { + CGAL_assertion( + pt == this->ptr()->_m_z_stacks.begin()->point() + ); + } + ); + const Z_stack& z_stack = + (static_cast< int >(this->ptr()->_m_z_stacks.size()) == 1 ? + Z_stack(*this->ptr()->_m_z_stacks.begin(), feature, this) : + Z_stack(feature,this,pt)); + + typename Rep::Surface_container erase; + + for (typename Rep::Surface_container::iterator it = + this->ptr()->_m_surfaces.begin(); + it != this->ptr()->_m_surfaces.end(); it++) { + + if (this->ptr()->_m_dcel_feature == SoX::EDGE && + _is_surface_vertical(*it)) { + continue; + } + + bool empty_isolator; + if (!z_stack._knows_isolator(*it, empty_isolator)) { + // construct traits + typename Surface_z_at_xy_isolator_traits:: + Construct_isolator + construct_isolator; + // TODO (construct_isolator_object()) + // construct isolator + boost::optional< + std::pair< SoX::Dcel_feature, CGAL::Object > > + dcelinfo = _dcel(*it); + CGAL_assertion(dcelinfo); + Z_at_xy_isolator isolator = + construct_isolator( + *it, pt, _nk(*it), dcelinfo->first + ); + // add roots in z_stack sequence + z_stack._add_surface(*it, isolator); + if (isolator.number_of_real_roots() == 0) { + // if isolator has no root remove + // surface from _m_surfaces + erase.insert(*it); + } + } else { + CGAL_assertion_code( + Z_at_xy_isolator isolator = z_stack._isolator(*it); + ); + CGAL_assertion(isolator.traits().point() == pt); + } + } + + for (typename Rep::Surface_container::iterator it = + erase.begin(); + it != erase.end(); it++) { + this->ptr()->_m_surfaces.erase(*it); + } + + this->ptr()->_m_z_stack = z_stack; + } + + this->ptr()->_m_z_stacks.clear(); + + CGAL_postcondition(this->ptr()->_m_z_stack); + + // check that each involved surface shows up in z_stack + CGAL_postcondition_code(( + { + std::list< int > levels; + for (typename Rep::Surface_container::iterator it = + this->ptr()->_m_surfaces.begin(); + it != this->ptr()->_m_surfaces.end(); it++){ + if (!_is_surface_vertical(*it)) { + levels.clear(); + typename + Surface_z_at_xy_isolator_traits::Construct_isolator + construct_isolator; + // TODO (construct_isolator_object()) + boost::optional< + std::pair< SoX::Dcel_feature, CGAL::Object > > + dcelinfo = _dcel(*it); + CGAL_assertion(dcelinfo); + Z_at_xy_isolator isolator = + construct_isolator( + *it, pt, _nk(*it), dcelinfo->first + ); + if (isolator.number_of_real_roots() > 0) { + this->ptr()->_m_z_stack->levels_of( + *it, std::back_inserter(levels) + ); + CGAL_postcondition(!levels.empty()); + } + } + } + } + )); + +#if !NDEBUG + std::cout << "done." << std::endl; +#endif + } + + //!@} + + //!\name Overlay + //!@{ + + Self _overlay_with(const Self& info, bool factors_of_same_curve) const { + Self tmp(info); + tmp.copy_on_write(); + + // set rs_id to zero + tmp.ptr()->_m_rs_id = 0; + + // empty object + CGAL::Object obj; + tmp.ptr()->_m_dcel_handle = obj; + + // involved surfaces + tmp.ptr()->_m_surfaces.insert(this->ptr()->_m_surfaces.begin(), + this->ptr()->_m_surfaces.end()); + + // silhouettes and ... + for (typename Silhouettes_nk_map::iterator it = + this->ptr()->_m_silhouettes_nk.begin(); + it != this->ptr()->_m_silhouettes_nk.end(); it++) { + tmp.ptr()->_m_silhouettes_nk.insert(*it); + } + + if (this->ptr()->_m_has_silhouette_curve) { + tmp.ptr()->_m_has_silhouette_curve = true; + } + + // cuts .. note that add_cut may also be additive + for (typename Cuts_map::iterator it = + this->ptr()->_m_cuts.begin(); + it != this->ptr()->_m_cuts.end(); it++) { + tmp._add_cut(it->first.first, it->first.second, + it->second, factors_of_same_curve); + } + + // overlay of boundary classes in x + if (tmp.ptr()->_m_boundary_type_in_x) { + if (this->ptr()->_m_boundary_type_in_x) { + CGAL_assertion(tmp.boundary_type_in_x() == + this->boundary_type_in_x()); + } + } else { + if (this->ptr()->_m_boundary_type_in_x) { + tmp.ptr()->_m_boundary_type_in_x = + this->boundary_type_in_x(); + } + } + // overlay of boundary classes in y + if (tmp.ptr()->_m_boundary_type_in_y) { + if (this->ptr()->_m_boundary_type_in_y) { + CGAL_assertion(tmp.boundary_type_in_y() == + this->boundary_type_in_y()); + } + } else { + if (this->ptr()->_m_boundary_type_in_y) { + tmp.ptr()->_m_boundary_type_in_y = + this->boundary_type_in_y(); + } + } + + // new construction may benefit from original z_stacks + // store to _m_z_stacks + if (tmp.ptr()->_m_z_stack) { + tmp.ptr()->_m_z_stacks.push_front(*(tmp.ptr()->_m_z_stack)); + } + if (this->ptr()->_m_z_stack) { + tmp.ptr()->_m_z_stacks.push_front(*(this->ptr()->_m_z_stack)); + } + // this list will be deleted when asking for z_stack of *this + + // but first, we invalidate current z_stack + tmp.ptr()->_m_point = boost::none; + tmp.ptr()->_m_z_stack = boost::none; + + if (!factors_of_same_curve) { + // surface dcel handles + for (typename Surface_dcel_handle_map::iterator it = + this->ptr()->_m_surface_dcel_handle_map.begin(); + it != this->ptr()->_m_surface_dcel_handle_map.end(); it++) { + tmp._set_dcel(it->first, it->second.first, it->second.second); + } + + // surface pair dcel handles + for (typename Surface_pair_dcel_handle_map::iterator it = + this->ptr()->_m_surface_pair_dcel_handle_map.begin(); + it != this->ptr()->_m_surface_pair_dcel_handle_map.end(); + it++) { + tmp._set_dcel( + it->first.first, it->first.second, + it->second.first, it->second.second + ); + } + } + + return tmp; + } + + //!@} + +public: + //!\name IO members + //!@{ + + /*!\brief + * prints pretty-formated version of P_dcel_info + */ + void pretty_print(std::ostream& os) const { + os << "PDcelInfo(id: " << this->id() + << ", rsid: " << this->_rs_id() + << ", feat=" << this->ptr()->_m_dcel_feature << "): " + << std::endl; + if (this->ptr()->_m_dcel_feature == SoX::VERTEX) { + typedef Restricted_cad_3< Surface_z_at_xy_isolator_traits> + Restricted_cad_3; + typename Restricted_cad_3::Vertex_const_iterator vit; + CGAL_assertion_code(bool check = ) + CGAL::assign(vit, this->ptr()->_m_dcel_handle); + CGAL_assertion(check); + os << "Point: " << vit->point() << std::endl; + } + if (this->ptr()->_m_dcel_feature == SoX::EDGE) { + typedef Restricted_cad_3< Surface_z_at_xy_isolator_traits> + Restricted_cad_3; + typename Restricted_cad_3::Halfedge_const_iterator hit; + CGAL_assertion_code(bool check = ) + CGAL::assign(hit, this->ptr()->_m_dcel_handle); + CGAL_assertion(check); + os << "Curve: " << hit->curve() << std::endl; + } + os << "Surfaces "; + for (typename Rep::Surface_container::iterator it = + this->ptr()->_m_surfaces.begin(); + it != this->ptr()->_m_surfaces.end(); it++){ + os << "<" << it->id() << ">" << std::flush; + } + if (this->ptr()->_m_dcel_feature != SoX::FACE) { + os << std::endl; + } + if (this->_has_silhouette()) { + os << "Silhouettes "; + for (typename Silhouettes_nk_map::iterator it = + this->ptr()->_m_silhouettes_nk.begin(); + it != this->ptr()->_m_silhouettes_nk.end(); it++) { + if (it->second.mult() != -1) { + os << "<" << it->first.id() << "," << it->second; + if (it->second.n() == 0) { + if (this->ptr()->_m_dcel_feature == SoX::EDGE) { + os << ",V"; + } else if (this->ptr()->_m_dcel_feature == + SoX::VERTEX) { + os << ",VL"; + } + } + os << "> "; + } + } + os << " "; + } + if (this->_has_cut()) { + if (this->_has_silhouette()) { + os << " "; + } + os << "Cuts "; + for (typename Cuts_map::iterator it = + this->ptr()->_m_cuts.begin(); + it != this->ptr()->_m_cuts.end(); it++) { + os << "<(" << it->first.first.id() << "," + << it->first.second.id() << ")," << it->second + //<< " (" << it->first.first.f() << "," + //<< it->first.second.f() << ")" + << "> "; + } + } + if (this->ptr()->_m_boundary_type_in_x || + this->ptr()->_m_boundary_type_in_y) { + if (this->ptr()->_m_boundary_type_in_x) { + os << "BdryX " << *this->ptr()->_m_boundary_type_in_x; + } + if (this->ptr()->_m_boundary_type_in_y) { + if (this->ptr()->_m_boundary_type_in_x) { + os << " "; + } + os << "BdryY " << *this->ptr()->_m_boundary_type_in_y; + } + if (this->ptr()->_m_dcel_feature != SoX::FACE) { + os << std::endl; + } + } + if (this->ptr()->_m_z_stack) { + os << *this->ptr()->_m_z_stack; + } + } + + //!@} + +private: + + // friends + // for z_stacks + friend class Restricted_cad_3< Surface_z_at_xy_isolator_traits>; + friend class Restricted_cad_3_accessor< Restricted_cad_3< + Surface_z_at_xy_isolator_traits> >; + + // for add members + friend class Create_restricted_cad_3< Surface_z_at_xy_isolator_traits >; + // for _set_rs_id + friend class Overlay_restricted_cad_3< Surface_z_at_xy_isolator_traits >; + + // for _overlay_with + friend class Arr_p_dcel_info_overlay_traits< + Restricted_cad_3< Surface_z_at_xy_isolator_traits> >; + + // for _z_stack_of_surface + friend class SoX::Z_stack< Surface_z_at_xy_isolator_traits, Self>; + + // data members +}; + + +/*!\relates P_dcel_info + * \brief + * outputs P_dcel_instance object to stream \c os + */ +template < class SurfaceZAtXyIsolatorTraits > +std::ostream& operator<<( + std::ostream& os, + const P_dcel_info< SurfaceZAtXyIsolatorTraits >& data +) { + data.pretty_print(os); + return os; +} + + +} // namespace SoX + +#endif // SoX_GAPS_P_DCEL_INFO_H +// EOF diff --git a/Curved_kernel_via_analysis_2/include/CGAL/Arrangement_2l/Restricted_cad_3.h b/Curved_kernel_via_analysis_2/include/CGAL/Arrangement_2l/Restricted_cad_3.h new file mode 100644 index 00000000000..59b61d1aa36 --- /dev/null +++ b/Curved_kernel_via_analysis_2/include/CGAL/Arrangement_2l/Restricted_cad_3.h @@ -0,0 +1,1629 @@ +// ============================================================================ +// +// Copyright (c) 2001-2008 Max-Planck-Institut Saarbruecken (Germany). +// All rights reserved. +// +// This file is part of EXACUS (http://www.mpi-inf.mpg.de/projects/EXACUS/); +// you may redistribute it under the terms of the Q Public License version 1.0. +// See the file LICENSE.QPL distributed with EXACUS. +// +// Licensees holding a valid commercial license may use this file in +// accordance with the commercial license agreement provided with the software. +// +// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE +// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. +// +// ---------------------------------------------------------------------------- +// +// Library : SoX +// File : include/SoX/GAPS/Restricted_cad_3.h +// SoX_release : $Name: $ +// Revision : $Revision$ +// Revision_date : $Date$ +// +// Author(s) : Eric Berberich +// +// ============================================================================ + +#ifndef SoX_GAPS_RESTRICTED_CAD_3_H +#define SoX_GAPS_RESTRICTED_CAD_3_H 1 + +/*! \file SoX/GAPS/Restricted_cad_3.h + \brief Contains class template Restricted_cad_3 +*/ + +#include + +#include +#include + +#include +#include + +#include +#include + +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include + +namespace SoX { + +namespace Intern { + +/*!\brief + * Functor class that provides creation and overlaying of + * restriced_cads. + */ +template < class SurfaceZAtXyIsolatorTraits > +class Restricted_cad_3_cache { +public: + + //! this instance's first template parameter + typedef SurfaceZAtXyIsolatorTraits Surface_z_at_xy_isolator_traits; + + //! type of instantiated class + typedef Restricted_cad_3_cache< Surface_z_at_xy_isolator_traits > Self; + + //! type of Restricted_cad_3 + typedef Restricted_cad_3< Surface_z_at_xy_isolator_traits > + Restricted_cad_3; + + //! type of id + typedef typename Restricted_cad_3::Id_type Id_type; + + //! type of Surface_3 + typedef typename Surface_z_at_xy_isolator_traits::Surface_3 Surface_3; + +private: + //! type of creator + typedef SoX::Create_restricted_cad_3< Surface_z_at_xy_isolator_traits > + Creator; + + //! less for cad + typedef CGAL::Handle_id_less_than< Restricted_cad_3 > Less; + + //! type of map for restricted cad + typedef std::map< Id_type, Restricted_cad_3, std::less< Id_type > > + Restricted_cad_map; + + //! type of Less for surfaces + typedef typename Surface_3::Surface_less_than Surface_less_than; + + //! type of map + typedef std::map< Surface_3, Id_type, Surface_less_than > + Surface_cache_3; + + //! type of Surface_pair + typedef std::pair< Surface_3, Surface_3 > Surface_pair_3; + + //! type of Less of surface pair + typedef CGAL::Pair_lexicographical_less_than< Surface_3, Surface_3, + Surface_less_than, Surface_less_than > + Surface_pair_less; + + //! type of pair map + typedef std::map< Surface_pair_3, Id_type, Surface_pair_less > + Surface_pair_cache_3; + + //! pair canonicalizer + struct Canonicalize_pair { + + //! first in pair in always less than second + Surface_pair_3 operator()(Surface_pair_3 pair) { + Surface_less_than less; + if (!less(pair.first, pair.second)) { + std::swap(pair.first, pair.second); + } + return pair; + } + }; + +public: + //!\name Caching for cads with ids + //!@{ + + //! returns cad of given \c id + Restricted_cad_3 operator()(Id_type id) { + typename Restricted_cad_map::iterator it = _m_cads.find(id); + CGAL_assertion(it != _m_cads.end()); + return it->second; + } + + //! cached \c cad + Restricted_cad_3 operator()(Restricted_cad_3 cad) { + typename Restricted_cad_map::iterator it = _m_cads.find(cad.id()); + CGAL_assertion(it == _m_cads.end()); + _m_cads.insert(it, std::make_pair(cad.id(), cad)); + return cad; + } + + //!@} + + //!\name Caching for cads representing the silhouettes of single surface + //!@{ + +public: + //! returns \c cad of silhouette-curve for \c surface + Restricted_cad_3 operator()(const Surface_3& surface) { + typename Surface_cache_3::iterator it = _m_map.find(surface); + if (it == _m_map.end()) { + Creator creator; // TODO single instance + typename Creator::Non_cached_construction_tag t; + Restricted_cad_3 tmp = creator(surface, t); + _m_cads.insert(std::make_pair(tmp.id(), tmp)); + it = _m_map.insert(it, std::make_pair(surface, tmp.id())); + } + return _m_cads.find(it->second)->second; + } + +private: + //! return \c true if cad for \c surface's silhouette-curve is known + bool is_cached(const Surface_3& surface) { + return (_m_map.find(surface) != _m_map.end()); + } + + //!@} + + //!\name Caching for cads representing the cuts of surface pairs + //!@{ + +public: + //! return stored \c cad of cut-curve for \c surface1 and \c surface2 + Restricted_cad_3 operator()(const Surface_3& surface1, + const Surface_3& surface2) { + static Canonicalize_pair canonicalize; + Surface_pair_3 pair = canonicalize(std::make_pair(surface1, surface2)); + typename Surface_pair_cache_3::iterator it = _m_pair_map.find(pair); + if (it == _m_pair_map.end()) { + Creator creator; // TODO single instance + typename Creator::Non_cached_construction_tag t; + Restricted_cad_3 tmp = creator(surface1, surface2, t); + _m_cads.insert(std::make_pair(tmp.id(), tmp)); + it = _m_pair_map.insert(it, std::make_pair(pair, tmp.id())); + } + return _m_cads.find(it->second)->second; + } + +private: + //! return \c true if cad for given surfaces cut-curve is known + bool is_cached(const Surface_3& surface1, const Surface_3& surface2) { + static Canonicalize_pair canonicalize; + Surface_pair_3 pair = canonicalize(std::make_pair(surface1, surface2)); + return (_m_pair_map.find(pair) != _m_pair_map.end()); + } + + //!@} + + // Remark: We omit to cache sil+sil+cut - as this can be done in + // Surface_pair_3, and other classes can derive from it + +private: + + // members + //! cache for for cads + Restricted_cad_map _m_cads; + + //! cache for silhouettes + Surface_cache_3 _m_map; + + //! cache for cuts + Surface_pair_cache_3 _m_pair_map; +}; + +} // namespace Intern + +// TODO document +template < class SurfaceZAtXyIsolatorTraits > +class Restricted_cad_3 : + public CGAL::Handle_with_policy < CGAL::Arrangement_2< +typename SurfaceZAtXyIsolatorTraits::Arrangement_traits_2, + CGAL::Arr_extended_dcel< + typename SurfaceZAtXyIsolatorTraits::Arrangement_traits_2, + boost::optional< SoX::P_dcel_info< SurfaceZAtXyIsolatorTraits > >, + boost::optional< SoX::P_dcel_info< SurfaceZAtXyIsolatorTraits > >, + boost::optional< SoX::P_dcel_info< SurfaceZAtXyIsolatorTraits > > +> +> +> { +public: + //! this instance' template parameter + typedef SurfaceZAtXyIsolatorTraits Surface_z_at_xy_isolator_traits; + + SoX_SURFACE_Z_AT_XY_ISOLATOR_TRAITS_SNAP_TYPEDEFS( + Surface_z_at_xy_isolator_traits + ); + +#if DOXYGEN_RUNNING + //! type of Surface_3 + typedef typename Surface_z_at_xy_isolator_traits::Surface_3 Surface_3; +#endif + + //! data type + typedef SoX::P_dcel_info< Surface_z_at_xy_isolator_traits > Data; + +private: + //! optionalized data + typedef boost::optional < Data > Opt_data; + +#if DOXYGEN_RUNNING + //! arrangment traits + typedef typename Surface_z_at_xy_isolator_traits::Arrangement_traits_2 + Arrangement_traits_2; +#endif + + //! dcel + typedef CGAL::Arr_extended_dcel< Arrangement_traits_2, + Opt_data, Opt_data, Opt_data > Dcel; + +public: + //! rep class + typedef CGAL::Arrangement_2< Arrangement_traits_2, Dcel > Rep; + +private: + //! base class + typedef CGAL::Handle_with_policy< Rep > Base; + + typedef Restricted_cad_3< Surface_z_at_xy_isolator_traits > Self; + + //! type of cache + typedef Intern::Restricted_cad_3_cache< Surface_z_at_xy_isolator_traits > + Restricted_cad_3_cache; + +public: + + //! Z_Stack type + typedef typename Data::Z_stack Z_stack; + + // from arrangement + typedef typename Dcel::Size Size; + +#if DOXYGEN_RUNNING + typedef typename Rep::Point_2 Point_2; + typedef typename Rep::X_monotone_curve_2 X_monotone_curve_2; +#endif + +private: + typedef typename Rep::Vertex_iterator Vertex_iterator; + typedef typename Rep::Edge_iterator Edge_iterator; + typedef typename Rep::Halfedge_iterator Halfedge_iterator; + typedef typename Rep::Face_iterator Face_iterator; + +public: + typedef typename Rep::Vertex_const_iterator Vertex_const_iterator; + + typedef typename Rep::Edge_const_iterator Edge_const_iterator; +private: + typedef typename Rep::Halfedge_const_iterator Halfedge_const_iterator; + +public: + typedef typename Rep::Face_const_iterator Face_const_iterator; + + typedef typename Rep::Ccb_halfedge_const_circulator + Ccb_halfedge_const_circulator; + + typedef typename Rep::Outer_ccb_const_iterator Outer_ccb_const_iterator; + + typedef typename Rep::Inner_ccb_const_iterator Inner_ccb_const_iterator; + + typedef typename Rep::Isolated_vertex_const_iterator + Isolated_vertex_const_iterator; + +public: + typedef typename Rep::Halfedge_around_vertex_const_circulator + Halfedge_around_vertex_const_circulator; + +private: + typedef Vertex_iterator Vertex_handle; + typedef Halfedge_iterator Halfedge_handle; + typedef Edge_iterator Edge_handle; + typedef Face_iterator Face_handle; + +public: + typedef Vertex_const_iterator Vertex_const_handle; +private: + typedef Halfedge_const_iterator Halfedge_const_handle; +public: + typedef Edge_const_iterator Edge_const_handle; + typedef Face_const_iterator Face_const_handle; + + // TODO choose better point location + typedef CGAL::Arr_naive_point_location< Rep > Point_location; + +public: + //! type of X_coordinate + typedef typename Curve_kernel_2::X_coordinate_1 X_coordinate_1; + + //! type of Boundary + typedef typename Curve_kernel_2::Boundary Boundary; + + //!\name Static members + //!@{ + + //! returns static instance of cad_cache for silhouette- and cut-curves + static Restricted_cad_3_cache& cad_cache() { + static Restricted_cad_3_cache cache; + return cache; + } + + //!@} + +private: + + //!\name Global + //!@{ + + //! make id unique and clear the rep + void _renew() { + this->copy_on_write(); + this->ptr()->clear(); + } + + //! construct cad for given curve + static + Restricted_cad_3 _construct_for_curve(const Curve_analysis_2& curve) { + + Restricted_cad_3 cad; + cad._renew(); + + // split curve into segments + + typename + Surface_z_at_xy_isolator_traits::Arrangement_traits_2 + geo_traits; + ; + + std::list< Curve_analysis_2 > input; + input.push_front(curve); + + std::list< X_monotone_curve_2 > xcurves; + std::list< Point_2 > points; + CGAL::make_x_monotone( + input.begin(), input.end(), + std::back_inserter(xcurves), + std::back_inserter(points), + &geo_traits + ); + + // TODO use non-intersecting + cad._insert_empty(xcurves.begin(), xcurves.end(), + points.begin(), points.end()); + + return cad; + } + + //! initializes a single feature + template < class DcelConstHandle > + void _init_feature(const Surface_3& surface, + DcelConstHandle handle, + SoX::Dcel_feature feature) { + + CGAL::Object obj = CGAL::make_object(handle); + this->non_const_handle(handle)->data() = Data(feature, obj); + handle->data()->_set_dcel(feature, obj); + handle->data()->_set_rs_id(this->id()); + + handle->data()->_set_dcel(surface, feature, obj); + + handle->data()->_init_nk(surface); + } + + //! initialize data objects + void _init(const Surface_3& surface) { + for (Vertex_const_iterator vit = + this->ptr()->vertices_begin(); + vit != this->ptr()->vertices_end(); vit++) { + + CGAL_precondition(!vit->data()); + + _init_feature(surface, vit, SoX::VERTEX); + + CGAL_postcondition(vit->data()); + } + + for (Halfedge_const_iterator hit = + this->ptr()->halfedges_begin(); + hit != this->ptr()->halfedges_end(); hit++) { + + if (!hit->data()) { + + CGAL_precondition(!hit->twin()->data()); + + + _init_feature(surface, hit, SoX::EDGE); + this->non_const_handle(hit->twin())->data() = *hit->data(); + } + + CGAL_postcondition(hit->data()); + CGAL_postcondition(hit->twin()->data()); + } + + for (Face_const_iterator fit = + this->ptr()->faces_begin(); + fit != this->ptr()->faces_end(); fit++) { + + CGAL_precondition(!fit->data()); + + _init_feature(surface, fit, SoX::FACE); + + CGAL_postcondition(fit->data()); + } + } + + //! set final data + void _finalize(const Surface_3& surface) { + + for (Vertex_const_iterator vit = + this->ptr()->vertices_begin(); + vit != this->ptr()->vertices_end(); vit++) { + CGAL_precondition(nk(vit,surface).mult() != -1); + vit->data()->_finalize_nk(surface); + } + + for (Edge_const_iterator eit = + this->ptr()->edges_begin(); + eit != this->ptr()->edges_end(); eit++) { + CGAL_precondition(nk(eit,surface).mult() != -1); + eit->data()->_finalize_nk(surface); + } + + for (Face_const_iterator fit = + this->ptr()->faces_begin(); + fit != this->ptr()->faces_end(); fit++) { + CGAL_precondition(nk(fit,surface).mult() == -1); + fit->data()->_finalize_nk(surface); + } + } + + //!@} + +private: + //!\name Arrangement methods + //!@{ + + /*! insert a range of curves */ + template < class InputIterator > + void _insert_non_intersecting_curves( + InputIterator begin, InputIterator end + ) { + CGAL::insert_non_intersecting_curves(*(this->ptr()), begin, end); + } + + /*! insert a range of curves */ + template < class CurveInputIterator, class PointInputIterator > + void _insert_empty( + CurveInputIterator cbegin, CurveInputIterator cend, + PointInputIterator pbegin, PointInputIterator pend + + ) { + CGAL::insert_empty(*(this->ptr()), + cbegin, cend, + pbegin, pend); + } + + /*! inserts an isolated point */ + void _insert_point(Point_2 pt) { + CGAL::insert_point(*(this->ptr()),pt); + } + + /*! remove edge from arrangement */ + void _remove_edge(Halfedge_handle he) { + CGAL::remove_edge(*(this->ptr()), he); + } + + /*! remove edge from arrangement */ + void _remove_vertex(Vertex_handle vh) { + CGAL::remove_vertex(*(this->ptr()), vh); + } + + /*! overlays two instances */ + void _overlay(Self rscad1, Self rscad2, bool factors_of_same_curve) { + // use correct overlay instance + SoX::Arr_p_dcel_info_overlay_traits< Self > + ovltraits(*this, factors_of_same_curve); + this->ptr()->clear(); + CGAL::overlay(*(rscad1.ptr()), *(rscad2.ptr()), + *(this->ptr()), ovltraits); + } + + /*! overlays two instances */ + void _overlay(Self rscad1, Self rscad2, + const Surface_3& surface, + SoX::Nk::Value_type type) { + // use correct overlay instance + SoX::Arr_p_dcel_info_overlay_traits< Self > + ovltraits(surface, *this, type); + this->ptr()->clear(); + CGAL::overlay(*(rscad1.ptr()), *(rscad2.ptr()), + *(this->ptr()), ovltraits); + } + + + //!@} + + //!\name Geometry and the Dcel + //!@{ + +public: + + /*! perform point location */ + CGAL::Object locate(Point_2 pt) const { + Point_location pl(*(this->ptr())); + return pl.locate(pt); + } + +private: + + /*!\brief + * use point location to check whether \c pt is part of given + * \c handle + */ + template < class DcelConstHandle > + bool _point_on_dcel_handle(const Point_2& pt, + DcelConstHandle handle) const { + + typedef DcelConstHandle Dcel_const_handle; + // Perform the point-location query. + CGAL::Object obj = this->locate(pt); + + // check if correct type + Dcel_const_handle handle2; + if (CGAL::assign (handle2, obj)) { + if (handle->data()->id() == handle2->data()->id()) { + return true; + } else { + return false; + } + } + // else + return false; + } + + /*!\brief + * constructs a rational point at coordinate \c x0, \c y0 + */ + static + Point_2 _construct_point_with_rational_y(X_coordinate_1 x0, Boundary y0) { + + //! type of Polynomial 1 + typedef typename Polynomial_2::NT Polynomial_1; + + typedef CGAL::Fraction_traits< Boundary > FT; + typedef typename FT::Numerator_type Integer; + Integer yn; + Integer yd; + typename FT::Decompose decompose; + decompose(y0, yn, yd); + + Polynomial_2 poly(Polynomial_1(-yn), Polynomial_1(yd)); + + // construct curve analysis + typename + Arrangement_traits_2::Curve_kernel_2::Construct_curve_2 + construct_curve = + Arrangement_traits_2::instance().kernel(). + construct_curve_2_object(); + typename + Arrangement_traits_2::Curve_kernel_2::Curve_analysis_2 hcurve = + construct_curve(poly); + Point_2 pt(x0, hcurve, 0); + return pt; + } + + //! construct point, with rational x (or y) in interior of \c heh's curve + static + Point_2 _point_in_interior(const Halfedge_const_handle& heh) { + Point_2 pt = _point_in_interior(heh->curve()); + CGAL_postcondition( + Restricted_cad_3::cad_cache()( + heh->data()->ptr()->_m_rs_id + )._point_on_dcel_handle(pt, heh) + ); + return pt; + } + + //! construct point, with rational x (or y) in interior of \c heh's curve + static + Point_2 _point_in_interior(const X_monotone_curve_2& cv) { + + Point_2 pt; +#if 1 + typename + Surface_z_at_xy_isolator_traits::Arrangement_traits_2 + geo_traits; + + typename Arrangement_traits_2::Construct_interior_vertex_2 + construct_interior_vertex = + geo_traits.construct_interior_vertex_2_object(); + + pt = construct_interior_vertex(cv); +#else + // TODO remove old code + //! type of Event info + typedef typename + Arrangement_traits_2::Curve_kernel_2::Curve_analysis_2:: + Status_line_1 + Status_line_1; + + // find point in halfedge interior + if (cv.is_vertical()) { + // handle also vertical curves + // find y0 between source.y and target.y + + Boundary y0(0); + + CGAL::Arr_parameter_space ps_y_min = + cv.location(CGAL::ARR_MIN_END); + CGAL::Arr_parameter_space ps_y_max = + cv.location(CGAL::ARR_MAX_END); + + Point_2 min_p; + if (ps_y_min == CGAL::ARR_INTERIOR) { + min_p = cv.curve_end(CGAL::ARR_MIN_END); + } + Point_2 max_p; + if (ps_y_max == CGAL::ARR_INTERIOR) { + max_p = cv.curve_end(CGAL::ARR_MAX_END); + } + + X_coordinate_1 x0 = cv.x(); + + typename Arrangement_traits_2::Curve_kernel_2::Lower_boundary_y_2 + lower_boundary_y = + Arrangement_traits_2::instance().kernel(). + lower_boundary_y_2_object(); + typename Arrangement_traits_2::Curve_kernel_2::Upper_boundary_y_2 + upper_boundary_y = + Arrangement_traits_2::instance().kernel(). + upper_boundary_y_2_object(); + + if (ps_y_min == CGAL::ARR_BOTTOM_BOUNDARY) { // -oo + if (ps_y_max == CGAL::ARR_TOP_BOUNDARY) { // +oo + // nothing to do since y0 == 0 already + } else { + Status_line_1 max_sl = + max_p.curve().status_line_at_exact_x(x0); + int max_i = max_p.arcno(); + y0 = lower_boundary_y(max_sl.algebraic_real_2(max_i)) - + Boundary(1); + } + } else { + if (ps_y_max == CGAL::ARR_TOP_BOUNDARY) { // +oo + Status_line_1 min_sl = + min_p.curve().status_line_at_exact_x(x0); + int min_i = min_p.arcno(); + y0 = upper_boundary_y(min_sl.algebraic_real_2(min_i)) + + Boundary(1); + } else { + // both endpoints are finite + Status_line_1 min_sl = + min_p.curve().status_line_at_exact_x(x0); + int min_i = min_p.arcno(); + + Status_line_1 max_sl = + max_p.curve().status_line_at_exact_x(x0); + int max_i = max_p.arcno(); + + typename + Arrangement_traits_2::Curve_kernel_2:: + Boundary_between_y_2 + boundary_between_y = + Arrangement_traits_2::instance().kernel(). + boundary_between_y_2_object(); + + y0 = boundary_between_y(min_sl.algebraic_real_2(min_i), + max_sl.algebraic_real_2(max_i)); + } + } + + pt = _construct_point_with_rational_y(x0, y0); + } else { + // construct point on curve + X_coordinate_1 x0(cv.boundary_in_x_range_interior()); + pt = Point_2(x0, cv.curve(), cv.arcno()); + } +#endif + return pt; + } + + //!@} + +public: + //!\name Counting + //!@{ + /*! Check whether the arrangement is empty. */ + bool is_empty () const + { + return (this->ptr()->is_empty()); + } + + /*! Get the number of arrangement vertices. */ + Size number_of_vertices () const + { + return (this->ptr()->number_of_vertices()); + } + + /*! Get the number of isolated arrangement vertices. */ + Size number_of_isolated_vertices () const + { + return (this->ptr()->number_of_isolated_vertices()); + } + +private: + /*! Get the number of arrangement halfedges (the result is always even). */ + Size number_of_halfedges () const + { + return (this->ptr()->number_of_halfedges()); + } + +public: + /*! Get the number of arrangement edges. */ + Size number_of_edges () const + { + return (this->ptr()->number_of_halfedges() / 2); + } + + /*! Get the number of arrangement faces. */ + Size number_of_faces () const + { + return (this->ptr()->number_of_faces()); + } + + /*! Get the number of arrangement faces. */ + Size number_of_unbounded_faces () const + { + return (this->ptr()->number_of_unbounded_faces()); + } + + //!@} + + //!\name Traversal functions for the arrangement vertices. + //!@{ + + /*! Get a const iterator for the first vertex in the arrangement. */ + Vertex_const_iterator vertices_begin() const + { + return (this->ptr()->vertices_begin()); + } + + /*! Get a past-the-end const iterator for the arrangement vertices. */ + Vertex_const_iterator vertices_end() const + { + return (this->ptr()->vertices_end()); + } + //!@} + +private: + //\name Traversal functions for the arrangement halfedges. + //!@{ + + /*! Get a const iterator for the first halfedge in the arrangement. */ + Halfedge_const_iterator halfedges_begin() const + { + return (this->ptr()->halfedges_begin()); + } + + /*! Get a past-the-end const iterator for the arrangement halfedges. */ + Halfedge_const_iterator halfedges_end() const + { + return (this->ptr()->halfedges_end()); + } + //!@} + +public: + //\name Traversal functions for the arrangement edges. + //!@{ + + /*! Get a const iterator for the first edge in the arrangement. */ + Edge_const_iterator edges_begin() const + { + return (this->ptr()->edges_begin()); + } + + /*! Get a past-the-end const iterator for the arrangement halfedges. */ + Edge_const_iterator edges_end() const + { + return (this->ptr()->edges_end()); + } + //!@} + + //!\name Traversal functions for the arrangement faces. + //!@{ + + /*! Get a const iterator for the first face in the arrangement. */ + Face_const_iterator faces_begin() const + { + return (this->ptr()->faces_begin()); + } + + /*! Get a past-the-end const iterator for the arrangement faces. */ + Face_const_iterator faces_end() const + { + return (this->ptr()->faces_end()); + } + //!@} + + //!\name Conversions of special halfedges to edges + //!{ + + //! convert to Edge_const_handle + Edge_const_handle convert_to_edge_const_handle( + const Halfedge_around_vertex_const_circulator& circ) const { + Halfedge_const_handle heh(circ); + Edge_const_handle eh = edges_begin(); + std::advance(eh, std::distance(halfedges_begin(),heh)/2); + return eh; + } + + //! convert to Edge_const_handle + Edge_const_handle convert_to_edge_const_handle( + const Ccb_halfedge_const_circulator& circ) const { + Halfedge_const_handle heh(circ); + Edge_const_handle eh = edges_begin(); + std::advance(eh, std::distance(halfedges_begin(),heh)/2); + return eh; + } + + //!@} + +private: + //!\name Casting away constness for handle types. + //!@{ + //! converts const vertex handle to non-const vertex handle + typename Rep::Vertex_handle non_const_handle (Vertex_const_handle vh) + { + return (this->ptr()->non_const_handle(vh)); + } + + //! converts const halfedge handle to non-const halfedge handle + typename Rep::Halfedge_handle non_const_handle (Halfedge_const_handle hh) + { + return (this->ptr()->non_const_handle(hh)); + } + + //! converts const face handle to non-const face handle + typename Rep::Face_handle non_const_handle (Face_const_handle fh) + { + return (this->ptr()->non_const_handle(fh)); + } + //!@} + + //!\name Nk + //!@{ + +private: + //! set value for vertex handle + void _set_nk_value(const Vertex_const_handle& vh, + const Surface_3& surface, + SoX::Nk::Value_type type, + int value) const { + switch (type) { + case SoX::Nk::MULT: + vh->data()->_nk(surface).set_mult(value); + break; + case SoX::Nk::N: + vh->data()->_nk(surface).set_n(value); + break; + case SoX::Nk::K: + vh->data()->_nk(surface).set_k(value); + break; + } + } + + //! set value for halfedge handle + void _set_nk_value(const Halfedge_const_handle& heh, + const Surface_3& surface, + SoX::Nk::Value_type type, + int value) const { + switch (type) { + case SoX::Nk::MULT: + heh->data()->_nk(surface).set_mult(value); + break; + case SoX::Nk::N: + heh->data()->_nk(surface).set_n(value); + break; + case SoX::Nk::K: + heh->data()->_nk(surface).set_k(value); + break; + } + } + + //! set value for edge handle + void _set_nk_value(const Edge_const_handle& eh, + const Surface_3& surface, + SoX::Nk::Value_type type, + int value) const { + switch (type) { + case SoX::Nk::MULT: + eh->data()->_nk(surface).set_mult(value); + break; + case SoX::Nk::N: + eh->data()->_nk(surface).set_n(value); + break; + case SoX::Nk::K: + eh->data()->_nk(surface).set_k(value); + break; + } + } + + //! set value for face handle + void _set_nk_value(const Face_const_handle& fh, + const Surface_3& surface, + SoX::Nk::Value_type type, + int value) const { + switch (type) { + case SoX::Nk::MULT: + fh->data()->_nk(surface).set_mult(value); + break; + case SoX::Nk::N: + fh->data()->_nk(surface).set_n(value); + break; + case SoX::Nk::K: + fh->data()->_nk(surface).set_k(value); + break; + } + } + +public: + //! nk values for vertex handle + const SoX::Nk& nk(const Vertex_const_handle& vh, + const Surface_3& surface) const { + return vh->data()->_nk(surface); + } + +private: + //! nk values for halfedge handle + const SoX::Nk& nk(const Halfedge_const_handle& heh, + const Surface_3& surface) const { + return heh->data()->_nk(surface); + } + +public: + //! nk values for edge handle + const SoX::Nk& nk(const Edge_const_handle& eh, + const Surface_3& surface) const { + return eh->data()->_nk(surface); + } + + //! nk values for face handle + const SoX::Nk& nk(const Face_const_handle& fh, + const Surface_3& surface) const { + return fh->data()->_nk(surface); + } + + //!@} + +public: + //!\name Silhouette/Cut + //!@{ + + //! returns whether silhouette at given vertex handle exists + bool has_silhouette(const Vertex_const_handle& vh) const { + return vh->data()->_has_silhouette(); + } +private: + //! returns whether silhouette at given halfedge handle exists + bool has_silhouette(const Halfedge_const_handle& heh) const { + return heh->data()->_has_silhouette(); + } +public: + //! returns whether silhouette at given edge handle exists + bool has_silhouette(const Edge_const_handle& eh) const { + return eh->data()->_has_silhouette(); + } + + //! returns whether silhouette at given face handle exists + bool has_silhouette(const Face_const_handle& fh) const { + return fh->data()->_has_silhouette(); + } + + + //! returns whether cut at given vertex handle exists + bool has_cut(const Vertex_const_handle& vh) const { + return vh->data()->_has_cut(); + } +private: + //! returns whether cut at given halfedge handle exists + bool has_cut(const Halfedge_const_handle& heh) const { + return heh->data()->_has_cut(); + } +public: + //! returns whether cut at given edge handle exists + bool has_cut(const Edge_const_handle& eh) const { + return eh->data()->_has_cut(); + } + + //! returns whether cut at given face handle exists + bool has_cut(const Face_const_handle& fh) const { + return fh->data()->_has_cut(); + } + + //! returns whether silhouette of \c surface at \c vh exisst + bool has_silhouette(const Vertex_const_handle& vh, + const Surface_3& surface) const { + return vh->data()->_has_silhouette(surface); + } +private: + //! returns whether silhouette of \c surface at \c heh exists + bool has_silhouette(const Halfedge_const_handle& heh, + const Surface_3& surface) const { + return heh->data()->_has_silhouette(surface); + } +public: + //! returns whether silhouette of \c surface at \c eh exisst + bool has_silhouette(const Edge_const_handle& eh, + const Surface_3& surface) const { + return eh->data()->_has_silhouette(surface); + } + + //! returns whether silhouette of \c surface at \c fh exisst + bool has_silhouette(const Face_const_handle& fh, + const Surface_3& surface) const { + return fh->data()->_has_silhouette(surface); + } + + //! returns whether \c vh supports vertical line of \c surface + bool supports_vertical_line(const Vertex_const_handle& vh, + const Surface_3& surface) const { + return vh->data()->_supports_vertical_line(surface); + } + + //! returns whether \c surface supports a vertical line at \c point + bool supports_vertical_line_at(const Point_2& point, + const Surface_3& surface) { + CGAL::Object obj = this->locate(point); + Face_const_handle fh; + if (CGAL::assign(fh,obj)) { + return false; + } + Halfedge_const_handle heh; + if (CGAL::assign(heh, obj)) { + return (this->nk(heh, surface).n() == -1); + } + Vertex_const_handle vh; + CGAL_assertion_code(bool check =) + CGAL::assign (vh, obj); + CGAL_assertion(check); + return (this->nk(vh, surface).n() == -1); + }; + + //! returns whether cut the two surface at \c vh exisst + bool has_cut(const Vertex_const_handle& vh, + const Surface_3& surface1, + const Surface_3& surface2) const { + return vh->data()->_has_cut(surface1, surface2); + } +private: + //! returns whether cut the two surface at \c heh exisst + bool has_cut(const Halfedge_const_handle& heh, + const Surface_3& surface1, + const Surface_3& surface2) const { + return heh->data()->_has_cut(surface1, surface2); + } +public: + //! returns whether cut the two surface at \c eh exisst + bool has_cut(const Edge_const_handle& eh, + const Surface_3& surface1, + const Surface_3& surface2) const { + return eh->data()->_has_cut(surface1, surface2); + } + + //! returns whether cut the two surface at \c fh exisst + bool has_cut(const Face_const_handle& fh, + const Surface_3& surface1, + const Surface_3& surface2) const { + return fh->data()->_has_cut(surface1, surface2); + } + + //!@} + + //!\name Sample Points + //!@{ + + //! returns sample point for vertex handle + Point_2 sample_point(const Vertex_const_handle& vh) const { + CGAL_precondition(vh->data()->_rs_id() == this->id()); + return vh->data()->_sample_point_for_vertex_handle(vh); + } +private: + //! returns sample point for given halfedge handle + Point_2 sample_point(const Halfedge_const_handle& heh) const { + CGAL_precondition(heh->data()->_rs_id() == this->id()); + return heh->data()->_sample_point_for_halfedge_handle(heh); + } +public: + //! returns sample point for given edge handle + Point_2 sample_point(const Edge_const_handle& eh) const { + CGAL_precondition(eh->data()->_rs_id() == this->id()); + Halfedge_const_handle heh = eh; + return heh->data()->_sample_point_for_halfedge_handle(heh); + } + + //! returns sample point for given edge handle + Point_2 sample_point(const Face_const_handle& fh) const { + CGAL_precondition(fh->data()->_rs_id() == this->id()); + return fh->data()->_sample_point_for_face_handle(fh); + } + + //!@} + + //!\name Isolators + //!@{ + + //! returns isolator (if existing) for given surface at given handle + boost::optional< Z_at_xy_isolator> + isolator(const Vertex_const_handle& vh, + const Surface_3& surface) const { + Z_stack z_stackp = z_stack(vh); + bool empty; + if (z_stackp._knows_isolator(surface, empty)) { + return z_stackp._isolator(surface); + } + // else + return boost::none; + } + +private: + //! returns isolator (if existing) for given surface at given handle + boost::optional< Z_at_xy_isolator> + isolator(const Halfedge_const_handle& heh, + const Surface_3& surface) const { + Z_stack z_stackp = z_stack(heh); + bool empty; + if (z_stackp._knows_isolator(surface, empty)) { + return z_stackp._isolator(surface); + } + // else + return boost::none; + } + +public: + //! returns isolator (if existing) for given surface at given handle + boost::optional< Z_at_xy_isolator> + isolator(const Edge_const_handle& eh, + const Surface_3& surface) const { + Z_stack z_stackp = z_stack(eh); + bool empty; + if (z_stackp._knows_isolator(surface, empty)) { + return z_stackp._isolator(surface); + } + // else + return boost::none; + } + + //! returns isolator (if existing) for given surface at given handle + boost::optional< Z_at_xy_isolator> + isolator(const Face_const_handle& fh, + const Surface_3& surface) const { + Z_stack z_stackp = z_stack(fh); + bool empty; + if (z_stackp._knows_isolator(surface, empty)) { + return z_stackp._isolator(surface); + } + // else + return boost::none; + } + + //! returns isolator (if existing) for \c surface for given point \c pt + boost::optional< Z_at_xy_isolator > + isolator_for(const Point_2& pt, + const Surface_3& surface) const { + Z_stack z_stackp = z_stack_for(pt); + bool empty; + if (z_stackp._knows_isolator(surface, empty)) { + return z_stackp._isolator(surface); + } + // else + return boost::none; + } + +#if 1 + // TODO reimplement avoiding cache! + //! returns isolator (if existing) for \c surface for given point \c pt + boost::optional< Z_at_xy_isolator > + isolator_at(const Point_2& pt, + const Surface_3& surface) const { + + + typename Surface_z_at_xy_isolator_traits:: + Construct_isolator + construct_isolator; + // TODO (construct_isolator_object()) + // construct isolator + + CGAL::Object obj = (cad_cache()(surface)).locate(pt); + + SoX::Nk nk; + SoX::Dcel_feature feature; + + Face_const_handle fh; + Halfedge_const_handle heh; + Vertex_const_handle vh; + if (CGAL::assign(fh,obj)) { + nk = fh->data()->_nk(surface); + feature = SoX::FACE; + } else if (CGAL::assign(heh, obj)) { + nk = heh->data()->_nk(surface); + feature = SoX::EDGE; + } else { + CGAL_assertion_code(bool check =) + CGAL::assign (vh, obj); + CGAL_assertion(check); + nk = vh->data()->_nk(surface); + feature = SoX::VERTEX; + } + + Z_at_xy_isolator isolator = + construct_isolator(surface, pt, nk, feature); + + return isolator; + } +#endif + + //!@} + + //!\name Z_Stacks + //!@{ + + //! returns z_stack for given vertex handle + Z_stack z_stack(const Vertex_const_handle& vh) const { + CGAL_precondition(vh->data()->_rs_id() == this->id()); + return vh->data()->_z_stack_for_vertex_handle(vh); + } +private: + //! returns z_stack for given halfedge handle + Z_stack z_stack(const Halfedge_const_handle& heh) const { + CGAL_precondition(heh->data()->_rs_id() == this->id()); + return heh->data()->_z_stack_for_halfedge_handle(heh); + } +public: + //! returns z_stack for given edge handle + Z_stack z_stack(const Edge_const_handle& eh) const { + CGAL_precondition(eh->data()->_rs_id() == this->id()); + Halfedge_const_handle heh = eh; + return heh->data()->_z_stack_for_halfedge_handle(heh); + } + + //! returns z_stack for given edge handle + Z_stack z_stack(const Face_const_handle& fh) const { + CGAL_precondition(fh->data()->_rs_id() == this->id()); + return fh->data()->_z_stack_for_face_handle(fh); + } + + //! returns z_stack for given point + Z_stack z_stack_for(const Point_2& point) const { + CGAL::Object obj = this->locate(point); + Face_const_handle fh; + if (CGAL::assign(fh,obj)) { + return this->z_stack(fh); + } + Halfedge_const_handle heh; + if (CGAL::assign(heh, obj)) { + return this->z_stack(heh); + } + Vertex_const_handle vh; + CGAL_assertion_code(bool check =) + CGAL::assign (vh, obj); + CGAL_assertion(check); + return this->z_stack(vh); + }; + +#if 0 + //! returns z_stack for given point + std::pair< Z_stack, SoX::Dcel_feature > + z_stack_at(const Point_2& point) const { + CGAL::Object obj = this->locate(point); + Face_const_handle fh; + if (CGAL::assign(fh,obj)) { + return std::make_pair(fh->data()->_z_stack(point, SoX::FACE), + SoX::FACE); + } + Halfedge_const_handle heh; + if (CGAL::assign(heh, obj)) { + return std::make_pair(heh->data()->_z_stack(point, SoX::EDGE), + SoX::EDGE); + } + Vertex_const_handle vh; + CGAL_assertion_code(bool check =) + CGAL::assign (vh, obj); + CGAL_assertion(check); + return std::make_pair(this->z_stack(vh), SoX::VERTEX); + }; +#endif + + //! returns z_stack of silhouette-curve of \c surface at given vertex handle + std::pair< Z_stack, SoX::Dcel_feature > + z_stack(const Vertex_const_handle& vh, + const Surface_3& surface) const { + return vh->data()->_z_stack_of_surface(surface); + } +private: + //! returns z_stack of silhouette-curve of \c surface at given edge handle + std::pair< Z_stack, SoX::Dcel_feature > + z_stack(const Halfedge_const_handle& heh, + const Surface_3& surface) const { + return heh->data()->_z_stack_of_surface(surface); + } +public: + //! returns z_stack of silhouette-curve of \c surface at given edge handle + std::pair < Z_stack, SoX::Dcel_feature > + z_stack(const Edge_const_handle& eh, + const Surface_3& surface) const { + return eh->data()->_z_stack_of_surface(surface); + } + + //! returns z_stack of silhouette-curve of \c surface at given face handle + std::pair < Z_stack, SoX::Dcel_feature > + z_stack(const Face_const_handle& fh, + const Surface_3& surface) const { + return fh->data()->_z_stack_of_surface(surface); + } + + //! returns z_stack for given point + std::pair< Z_stack, SoX::Dcel_feature > + z_stack_for(const Point_2& point, const Surface_3& surface) const { + CGAL::Object obj = this->locate(point); + Face_const_handle fh; + if (CGAL::assign(fh,obj)) { + return this->z_stack(fh, surface); + } + Halfedge_const_handle heh; + if (CGAL::assign(heh, obj)) { + return this->z_stack(heh, surface); + } + Vertex_const_handle vh; + CGAL_assertion_code(bool check =) + CGAL::assign (vh, obj); + CGAL_assertion(check); + return this->z_stack(vh, surface); + }; + +#if 0 + //! returns z_stack for given point + std::pair< Z_stack, SoX::Dcel_feature > + z_stack_at(const Point_2& point, const Surface_3& surface) const { + CGAL::Object obj = this->locate(point); + Face_const_handle fh; + if (CGAL::assign(fh,obj)) { + return fh->data()->_z_stack_of_surface_at(point, surface); + } + Halfedge_const_handle heh; + if (CGAL::assign(heh, obj)) { + return heh->data()->_z_stack_of_surface_at(point, surface); + } + Vertex_const_handle vh; + CGAL_assertion_code(bool check =) + CGAL::assign (vh, obj); + CGAL_assertion(check); + return this->z_stack(vh, surface); + }; +#endif + //!@} + + //!\name Adjacencies + //!{@ + + /*\brief Computes adjacency when going from \c surface sheet \c sheet in + * \c from to \c to, and return the corresponding z-cell as \c .first + * (std::pair< Z_stack, std::pair >) and the sheet numbers + * of \c surface at \c to as \c .second . + */ + template < class DcelConstHandle1, class DcelConstHandle2 > + std::pair< std::pair< Z_stack, std::pair< int, int > >, + std::pair< int, int > > adjacency( + DcelConstHandle1 from, + const Surface_3& surface, int sheet, + DcelConstHandle2 to) const { + + std::pair< Z_stack, SoX::Dcel_feature > z_stack_from = + z_stack(from, surface); + std::pair< Z_stack, SoX::Dcel_feature > z_stack_to = + z_stack(to, surface); + + CGAL_assertion(from->data()->_rs_id() == to->data()->_rs_id()); + + boost::optional< + std::pair< SoX::Dcel_feature,CGAL::Object > + > sffrom = from->data()->_dcel(surface); + CGAL_assertion(sffrom); + + boost::optional< + std::pair< SoX::Dcel_feature,CGAL::Object > + > sfto = to->data()->_dcel(surface); + CGAL_assertion(sfto); + + CGAL::Object fromobj = sffrom->second; + CGAL::Object toobj = sfto->second; + + std::pair< int, int > sheet_to = + z_stack_from.first.adjacency( + surface, sheet, + z_stack_from.second, + fromobj, + z_stack_to.first, + z_stack_to.second, + toobj + ); + + Z_stack goal = z_stack(to); + + if (sheet_to.first > sheet_to.second) { + // not involved + return std::make_pair(std::make_pair(goal, sheet_to), sheet_to); + } + if (sheet_to.first == -1) { + if (sheet_to.second == -1) { + // TODO handling at minus inf + return std::make_pair(std::make_pair(goal, sheet_to), + sheet_to); + } else if (sheet_to.second == + z_stack_to.first.number_of_z_cells()) { + int z = goal.number_of_z_cells(); + return std::make_pair( + std::make_pair(goal, + std::make_pair(-1, z)), + sheet_to + ); + } else { + int z_to_high = + goal.z_level_of_sheet(surface, sheet_to.second); + return std::make_pair( + std::make_pair(goal, + std::make_pair(-1, z_to_high)), + sheet_to + ); + } + } else if (sheet_to.second == z_stack_to.first.number_of_z_cells()) { + int z = goal.number_of_z_cells(); + if (sheet_to.first == z_stack_to.first.number_of_z_cells()) { + // TODO handling at plus inf + return std::make_pair( + std::make_pair(goal, std::make_pair(z,z)), + sheet_to + ); + } else { + int z_to_low = goal.z_level_of_sheet(surface, sheet_to.first); + + return std::make_pair( + std::make_pair(goal, std::make_pair(z_to_low, z)), + sheet_to + ); + } + } + + // else + int z_to_low = goal.z_level_of_sheet(surface, sheet_to.first); + int z_to_high = goal.z_level_of_sheet(surface, sheet_to.second); + + return std::make_pair( + std::make_pair(goal, std::make_pair(z_to_low, z_to_high)), + sheet_to + ); + } + + //!@} + + //!\name Surface Adjacency + //!@{ + + /*!\brief Computes adjacency when going from \c surface sheet \c sheet in + * \c from to \c to, and return the corresponding in sil-arr as .first + * and the sheet number of \c surface + * at \c to as \c .second . + */ + template < class DcelConstHandle1, class DcelConstHandle2 > + std::pair< Z_stack, std::pair< int, int > > surface_adjacency( + DcelConstHandle1 from, + const Surface_3& surface, int sheet, + DcelConstHandle2 to) const { + + std::pair< Z_stack, SoX::Dcel_feature > z_stack_from = + z_stack(from, surface); + std::pair< Z_stack, SoX::Dcel_feature > z_stack_to = + z_stack(to, surface); + + CGAL_assertion(from->data()->_rs_id() == to->data()->_rs_id()); + + boost::optional< + std::pair< SoX::Dcel_feature,CGAL::Object > + > sffrom = from->data()->_dcel(surface); + CGAL_assertion(sffrom); + + boost::optional< + std::pair< SoX::Dcel_feature,CGAL::Object > + > sfto = to->data()->_dcel(surface); + CGAL_assertion(sfto); + + CGAL::Object fromobj = sffrom->second; + CGAL::Object toobj = sfto->second; + + std::pair< int, int > sheet_to = + z_stack_from.first.adjacency( + surface, sheet, + z_stack_from.second, + fromobj, + z_stack_to.first, + z_stack_to.second, + toobj + ); + + return std::make_pair(z_stack_to.first, sheet_to); + } + + template < class DcelConstHandle1, class DcelConstHandle2 > + SoX::Adjacencies_3 adjacency(DcelConstHandle1 from, + const Surface_3& surface, + DcelConstHandle2 to) const { + + std::pair< Z_stack, SoX::Dcel_feature > z_stack_from = + z_stack(from, surface); + std::pair< Z_stack, SoX::Dcel_feature > z_stack_to = + z_stack(to, surface); + + CGAL_assertion(from->data()->_rs_id() == to->data()->_rs_id()); + + boost::optional< + std::pair< SoX::Dcel_feature,CGAL::Object > > sffrom = + from->data()->_dcel(surface); + CGAL_assertion(sffrom); + + boost::optional< + std::pair< SoX::Dcel_feature,CGAL::Object > > sfto = + to->data()->_dcel(surface); + CGAL_assertion(sfto); + + CGAL::Object fromobj = sffrom->second; + CGAL::Object toobj = sfto->second; + + typename Surface_z_at_xy_isolator_traits::Adjacency + compute_adj; // TODO single instance + + SoX::Adjacencies_3 adj = + compute_adj(surface, + // is valid as we just constructed z-stack + *isolator(from, surface), + z_stack_from.second, fromobj, + from->data()->_supports_vertical_line(surface), + // is valid as we just constructed z-stack + *isolator(to, surface), + z_stack_to.second, toobj, + to->data()->_supports_vertical_line(surface) + ); + return adj; + } + + //!@} + + // friends + // for non_const_handle, _insert*, _overlay, halfedge-stuff + friend class + SoX::Create_restricted_cad_3< Surface_z_at_xy_isolator_traits >; + friend class + SoX::Overlay_restricted_cad_3< Surface_z_at_xy_isolator_traits >; + friend class SoX::P_dcel_info< Surface_z_at_xy_isolator_traits >; + friend class SoX::Arr_p_dcel_info_overlay_traits< Self >; + // for types + friend class SoX::Z_stack< Surface_z_at_xy_isolator_traits, Data >; + // for all private members + friend class SoX::Restricted_cad_3_accessor< Self >; + friend class SoX::Intern::Construct_nk_decomposition_2< + Surface_z_at_xy_isolator_traits >; + +}; + +} // namespace SoX + +#endif // SoX_GAPS_RESTRICTED_CAD_3_H +// EOF diff --git a/Curved_kernel_via_analysis_2/include/CGAL/Arrangement_2l/Restricted_cad_3_accessor.h b/Curved_kernel_via_analysis_2/include/CGAL/Arrangement_2l/Restricted_cad_3_accessor.h new file mode 100644 index 00000000000..bf7dd5a2fea --- /dev/null +++ b/Curved_kernel_via_analysis_2/include/CGAL/Arrangement_2l/Restricted_cad_3_accessor.h @@ -0,0 +1,441 @@ +// ============================================================================ +// +// Copyright (c) 2001-2008 Max-Planck-Institut Saarbruecken (Germany). +// All rights reserved. +// +// This file is part of EXACUS (http://www.mpi-inf.mpg.de/projects/EXACUS/); +// you may redistribute it under the terms of the Q Public License version 1.0. +// See the file LICENSE.QPL distributed with EXACUS. +// +// Licensees holding a valid commercial license may use this file in +// accordance with the commercial license agreement provided with the software. +// +// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE +// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. +// +// ---------------------------------------------------------------------------- +// +// Library : SoX +// File : include/SoX/GAPS/Restricted_cad_3_accessor.h +// SoX_release : $Name: $ +// Revision : $Revision$ +// Revision_date : $Date$ +// +// Author(s) : Eric Berberich +// +// ============================================================================ + +#ifndef SoX_GAPS_RESTRICTED_CAD_3_ACCESSOR_H +#define SoX_GAPS_RESTRICTED_CAD_3_ACCESSOR_H 1 + +/*! \file SoX/GAPS/Restricted_cad_3_accessor.h + \brief Contains class template Restricted_cad_3_accessor +*/ + +#include + +#include + +namespace SoX { + +/*!\brief + * Accessor class for private types and members of RestrictedCad_3 + */ +template < class RestrictedCad_3 > +class Restricted_cad_3_accessor { +public: + + //! this instance's first template parameter + typedef RestrictedCad_3 Restricted_cad_3; + + //! type of instantiated class + typedef Restricted_cad_3_accessor< Restricted_cad_3 > Self; + + //! standard constructor + Restricted_cad_3_accessor(const Restricted_cad_3& rscad) : + _m_rscad(rscad) { + } + +private: + typedef Restricted_cad_3 RSC_3; + +public: + //! rep class + typedef typename RSC_3::Rep Rep; + + //! base class + typedef typename RSC_3::Base Base; + + //! type of Surface + typedef typename RSC_3::Surface_3 Surface_3; + + //! Z_Stack type + typedef typename RSC_3::Z_stack Z_stack; + + //! Z_at_xy_isolator type + typedef typename RSC_3::Z_at_xy_isolator Z_at_xy_isolator; + + // from arrangement + typedef typename RSC_3::Size Size; + + typedef typename RSC_3::Curve_analysis_2 Curve_analysis_2; + + typedef typename RSC_3::Point_2 Point_2; + typedef typename RSC_3::X_monotone_curve_2 X_monotone_curve_2; + + typedef typename RSC_3::Vertex_const_iterator Vertex_const_iterator; + + typedef typename RSC_3::Edge_const_iterator Edge_const_iterator; + + typedef typename RSC_3::Halfedge_const_iterator Halfedge_const_iterator; + + typedef typename RSC_3::Face_const_iterator Face_const_iterator; + + typedef typename RSC_3::Ccb_halfedge_const_circulator + Ccb_halfedge_const_circulator; + + typedef typename RSC_3::Outer_ccb_const_iterator Outer_ccb_const_iterator; + + typedef typename RSC_3::Inner_ccb_const_iterator Inner_ccb_const_iterator; + + typedef typename RSC_3::Vertex_const_handle Vertex_const_handle; + typedef typename RSC_3::Halfedge_const_handle Halfedge_const_handle; + typedef typename RSC_3::Edge_const_handle Edge_const_handle; + typedef typename RSC_3::Face_const_handle Face_const_handle; + + typedef typename RSC_3::Vertex_handle Vertex_handle; + typedef typename RSC_3::Halfedge_handle Halfedge_handle; + typedef typename RSC_3::Edge_handle Edge_handle; + typedef typename RSC_3::Face_handle Face_handle; + + typedef typename RSC_3::Point_location Point_location; + + //! type of X_coordinate + typedef typename RSC_3::X_coordinate_1 X_coordinate_1; + + //! type of Boundary + typedef typename RSC_3::Boundary Boundary; + + //!\name Global + //!@{ + + //! Returns a pointer to the representation of the cad + inline const Rep* rep() const { + return _m_rscad.ptr(); + } + + //! make id unique and clear the rep + inline void renew() { + _m_rscad._renew(); + } + + //! construct cad for curve + static + Restricted_cad_3 construct_for_curve(const Curve_analysis_2& curve) { + return Restricted_cad_3::_construct_for_curve(curve); + } + + //! initialize data objects + inline void init(const Surface_3& surface) { + _m_rscad._init(surface); + } + + //! set final data + inline void finalize(const Surface_3& surface) { + _m_rscad._finalize(surface); + } + + //!@} + + //!\name Arrangement methods + //!@{ + + /*! insert a range of curves */ + template < class InputIterator > + inline + void insert_non_intersecting_curves( + InputIterator begin, InputIterator end + ) { + _m_rscad._insert_non_intersecting_curves(begin, end); + } + + /*! insert a range of curves */ + template < class CurveInputIterator, class PointInputIterator > + inline + void insert_empty( + CurveInputIterator cbegin, CurveInputIterator cend, + PointInputIterator pbegin, PointInputIterator pend + + ) { + _m_rscad._insert_empty(cbegin, cend, pbegin, pend); + } + + + /*! inserts an isolated point */ + inline + void insert_point(Point_2 pt) { + _m_rscad._insert_point(pt); + } + + /*! remove edge from arrangement */ + void remove_edge(Halfedge_handle he) { + _m_rscad._remove_edge(he); + } + + /*! remove edge from arrangement */ + void remove_vertex(Vertex_handle vh) { + _m_rscad._remove_vertex(vh); + } + + /*! overlays two instances */ + inline + void overlay(Restricted_cad_3 rscad1, + Restricted_cad_3 rscad2, bool factors_of_same_curve) { + _m_rscad._overlay(rscad1, rscad2, factors_of_same_curve); + } + + /*! overlays two instances */ + inline + void overlay(Restricted_cad_3 rscad1, + Restricted_cad_3 rscad2, + Surface_3 surface, SoX::Nk::Value_type type) { + _m_rscad._overlay(rscad1, rscad2, surface, type); + } + + //!@} + + //!\name Geometry and the Dcel + //!@{ + + /*!\brief + * use point location to check whether \c pt is part of given + * \c handle + */ + template < class DcelConstHandle > + inline + bool point_on_dcel_handle(const Point_2& pt, + DcelConstHandle handle) const { + return _m_rscad._point_on_dcel_handle(pt, handle); + } + + /*!\brief + * constructs a rational point at coordinate \c x0, \c y0 + */ + static + Point_2 construct_point_with_rational_y(X_coordinate_1 x0, Boundary y0) { + return RSC_3::_construct_point_with_rational_y(x0, y0); + } + + //! construct point, with rational x (or y) in interior of \c heh's curve + static + Point_2 point_in_interior(const Halfedge_const_handle& heh) { + return RSC_3::_point_in_interior(heh); + } + + //! construct point, with rational x (or y) in interior of \c cv + static + Point_2 point_in_interior(const X_monotone_curve_2& cv) { + return RSC_3::_point_in_interior(cv); + } + + //!@} + + //!\name Counting + //!@{ + /*! Check whether the arrangement is empty. */ + + /*! Get the number of arrangement halfedges (the result is always even). */ + inline + Size number_of_halfedges () const + { + return (_m_rscad.number_of_halfedges()); + } + + //!@} + + //\name Traversal functions for the arrangement halfedges. + //!@{ + + /*! Get a const iterator for the first halfedge in the arrangement. */ + inline + Halfedge_const_iterator halfedges_begin() const + { + return (_m_rscad.halfedges_begin()); + } + + /*! Get a past-the-end const iterator for the arrangement halfedges. */ + inline + Halfedge_const_iterator halfedges_end() const + { + return (_m_rscad.halfedges_end()); + } + //!@} + + //!\name Casting away constness for handle types. + //!@{ + + //! converts const vertex handle to non-const vertex handle + inline + typename Rep::Vertex_handle non_const_handle (Vertex_const_handle vh) + { + return (_m_rscad.non_const_handle(vh)); + } + + //! converts const halfedge handle to non-const halfedge handle + inline + typename Rep::Halfedge_handle non_const_handle (Halfedge_const_handle hh) + { + return (_m_rscad.non_const_handle(hh)); + } + + //! converts const face handle to non-const face handle + inline + typename Rep::Face_handle non_const_handle (Face_const_handle fh) + { + return (_m_rscad.non_const_handle(fh)); + } + //!@} + + //!\name Nk + //!@{ + + + //! set value for vertex handle + inline void set_nk_value(const Vertex_const_handle& vh, + const Surface_3& surface, + SoX::Nk::Value_type type, + int value) const { + _m_rscad._set_nk_value(vh, surface, type, value); + } + + //! set value for halfedge handle + inline void set_nk_value(const Halfedge_const_handle& heh, + const Surface_3& surface, + SoX::Nk::Value_type type, + int value) const { + _m_rscad._set_nk_value(heh, surface, type, value); + } + + //! set value for edge handle + inline void set_nk_value(const Edge_const_handle& eh, + const Surface_3& surface, + SoX::Nk::Value_type type, + int value) const { + _m_rscad._set_nk_value(eh, surface, type, value); + } + + //! set value for face handle + inline void set_nk_value(const Face_const_handle& fh, + const Surface_3& surface, + SoX::Nk::Value_type type, + int value) const { + _m_rscad._set_nk_value(fh, surface, type, value); + } + + //! nk for halfedge handle + const SoX::Nk& nk(const Halfedge_const_handle& heh, + const Surface_3& surface) const { + return _m_rscad.nk(heh, surface); + } + + //!@} + + + //!\name Silhouette/Cut + //!@{ + + //! returns whether silhouette at given halfedge handle exists + bool has_silhouette(const Halfedge_const_handle& heh) const { + return (_m_rscad.has_silhouette(heh)); + } + + //! returns whether cut at given halfedge handle exists + bool has_cut(const Halfedge_const_handle& heh) const { + return (_m_rscad.has_cut(heh)); + } + + //! returns whether silhouette of \c surface at \c heh exists + bool has_silhouette(const Halfedge_const_handle& heh, + const Surface_3& surface) const { + return (_m_rscad.has_silhouette(heh,surface)); + } + + //! returns whether cut the two surface at \c heh exisst + bool has_cut(const Halfedge_const_handle& heh, + const Surface_3& surface1, + const Surface_3& surface2) const { + return (_m_rscad.has_cut(heh,surface1, surface2)); + } + + //!@} + + //!@{ + + //! returns sample point for halfedge handle + inline + Point_2 sample_point(const Halfedge_const_handle& heh) const { + return _m_rscad.sample_point(heh); + } + + //!@} + + //!\name Isolators + //!@{ + + //! returns isolator (if existing) for given surface at given handle + boost::optional< Z_at_xy_isolator> + isolator(const Halfedge_const_handle& heh, + const Surface_3& surface) const { + return _m_rscad_isolator(heh, surface); + } + + //!@} + + //!\name Z_Stacks + //!@{ + + //! returns whether z_stack is known + bool knows_z_stack(const Vertex_const_handle& vh) const { + return vh->data()->_knows_z_stack(); + } + + //! returns whether z_stack is known + bool knows_z_stack(const Halfedge_const_handle& heh) const { + return heh->data()->_knows_z_stack(); + } + + //! returns whether z_stack is known + bool knows_z_stack(const Edge_const_handle& eh) const { + return eh->data()->_knows_z_stack(); + } + + //! returns whether z_stack is known + bool knows_z_stack(const Face_const_handle& fh) const { + return fh->data()->_knows_z_stack(); + } + + //! returns z_stack for given edge handle + inline + Z_stack z_stack(const Halfedge_const_handle& heh) const { + return _m_rscad.z_stack(heh); + } + + //! returns z_stack of silhouette-curve of \c surface at halfedge handle + inline + std::pair< Z_stack, SoX::Dcel_feature > + z_stack(const Halfedge_const_handle& heh, + const Surface_3& surface) const { + return _m_rscad.z_stack(heh, surface); + } + + //!@} + +private: + Restricted_cad_3 _m_rscad; + +}; + +} // namespace SoX + +#endif // SoX_GAPS_RESTRICTED_CAD_3_ACCESSOR_H +// EOF diff --git a/Curved_kernel_via_analysis_2/include/CGAL/Arrangement_2l/Restricted_cad_3_enums.h b/Curved_kernel_via_analysis_2/include/CGAL/Arrangement_2l/Restricted_cad_3_enums.h new file mode 100644 index 00000000000..448c43e35b8 --- /dev/null +++ b/Curved_kernel_via_analysis_2/include/CGAL/Arrangement_2l/Restricted_cad_3_enums.h @@ -0,0 +1,171 @@ +// ============================================================================ +// +// Copyright (c) 2001-2008 Max-Planck-Institut Saarbruecken (Germany). +// All rights reserved. +// +// This file is part of EXACUS (http://www.mpi-inf.mpg.de/projects/EXACUS/); +// you may redistribute it under the terms of the Q Public License version 1.0. +// See the file LICENSE.QPL distributed with EXACUS. +// +// Licensees holding a valid commercial license may use this file in +// accordance with the commercial license agreement provided with the software. +// +// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE +// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. +// +// ---------------------------------------------------------------------------- +// +// Library : SoX +// File : include/SoX/GAPS/Restricted_cad_3_enums.h +// SoX_release : $Name: $ +// Revision : $Revision$ +// Revision_date : $Date$ +// +// Author(s) : Eric Berberich +// +// ============================================================================ + +#ifndef SoX_GAPS_RESTRICTED_CAD_3_ENUMS_H +#define SoX_GAPS_RESTRICTED_CAD_3_ENUMS_H 1 + +/*! \file SoX/GAPS/Restricted_cad_3_enums.h + \brief Contains enumeration related to + \link Restricted_cad_3 \endlink +*/ + +#include + +#include + +namespace SoX { + +//! distinguishes between the three possible dcel-features +enum Dcel_feature { + VERTEX = 0, //!< identifies a vertex + EDGE = 1, //!< identifies an edge + FACE = 2 //!< identifies a face +}; + +/*!\relates Dcel_feature + \brief output operator + + The output is intended for debugging purposes only and hence always + is in a human-readable pretty-print format. + + +*/ +// \pre ::CGAL::is_pretty(os) +inline std::ostream& operator<< (std::ostream& os, Dcel_feature feat) { + //CGAL_precondition(::CGAL::is_pretty(os)); + static const char* names[] = { "VERTEX", "EDGE", "FACE" }; + + CGAL_assertion(feat >= 0 && + feat < static_cast(sizeof names / sizeof *names)); + return os << names[feat]; +} + + +//! container for integers +class Nk { +public: + + enum Value_type { + MULT = 1, + N = 2, + K = 3 + }; + + Nk() : + _m_feat1(SoX::VERTEX), + _m_feat2(SoX::VERTEX), + _m_mult(-1), + // TODO other initial value? + _m_n(-2), + _m_k(-2), + _m_fixed(false) { + } + +public: + + SoX::Dcel_feature feature1() const { + return _m_feat1; + } + + SoX::Dcel_feature feature2() const { + return _m_feat2; + } + + void set_feature1(SoX::Dcel_feature feature) const { + CGAL_precondition(!_m_fixed); + _m_feat1 = feature; + } + + void set_feature2(SoX::Dcel_feature feature) const { + CGAL_precondition(!_m_fixed); + _m_feat2 = feature; + } + + int mult() const { + return _m_mult; + } + + int n() const { + return _m_n; + } + + int k() const { + return _m_k; + } + + void set_mult(int mult) const { + CGAL_precondition(mult >= 0); + CGAL_precondition(!_m_fixed); + _m_mult = mult; + } + + void set_n(int n) const { + CGAL_precondition(n >= -1); + CGAL_precondition(!_m_fixed); + _m_n = n; + } + + void set_k(int k) const { + CGAL_precondition(k >= -1); + CGAL_precondition(!_m_fixed); + _m_k = k; + } + + void fix() const { + _m_fixed = true; + } + + bool same_nk(const Nk& nk) const { + return + (this->_m_n == nk._m_n) && + (this->_m_k == nk._m_k); + } + +private: + // members + mutable SoX::Dcel_feature _m_feat1; + mutable SoX::Dcel_feature _m_feat2; + + + mutable int _m_mult; + mutable int _m_n; + mutable int _m_k; + + mutable bool _m_fixed; +}; + + +inline +std::ostream& operator<<(std::ostream& out, const SoX::Nk& nk) { + out << "NK(mult=" << nk.mult() << ",n=" << nk.n() + << ",k=" << nk.k() << ")"; + return out; +} + +} // namespace SoX + +#endif // SoX_GAPS_RESTRICTED_CAD_3_ENUMS_H diff --git a/Curved_kernel_via_analysis_2/include/CGAL/Arrangement_2l/Restricted_cad_3_functors.h b/Curved_kernel_via_analysis_2/include/CGAL/Arrangement_2l/Restricted_cad_3_functors.h new file mode 100644 index 00000000000..454f8786912 --- /dev/null +++ b/Curved_kernel_via_analysis_2/include/CGAL/Arrangement_2l/Restricted_cad_3_functors.h @@ -0,0 +1,1519 @@ +// ============================================================================ +// +// Copyright (c) 2001-2008 Max-Planck-Institut Saarbruecken (Germany). +// All rights reserved. +// +// This file is part of EXACUS (http://www.mpi-inf.mpg.de/projects/EXACUS/); +// you may redistribute it under the terms of the Q Public License version 1.0. +// See the file LICENSE.QPL distributed with EXACUS. +// +// Licensees holding a valid commercial license may use this file in +// accordance with the commercial license agreement provided with the software. +// +// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE +// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. +// +// ---------------------------------------------------------------------------- +// +// Library : SoX +// File : include/SoX/GAPS/Restricted_cad_3_functors.h +// SoX_release : $Name: $ +// Revision : $Revision$ +// Revision_date : $Date$ +// +// Author(s) : Eric Berberich +// +// ============================================================================ + +#ifndef SoX_GAPS_RESTRICTED_CAD_3_FUNCTORS_H +#define SoX_GAPS_RESTRICTED_CAD_3_FUNCTORS_H 1 + +/*! \file SoX/GAPS/Restricted_cad_3_functor.h + \brief Contains functor classes related to + \link Restricted_cad_3 \endlink +*/ + +#include + +#include +#include + +// TODO remove Poly_traits_d +#include + +#include +#include +#include +#include + +#include +#include +#include +#include + +namespace SoX { + +namespace Intern { + +/*!\brief + * Computes projected cut curves for two two given surfaces + */ +template < class SurfaceZAtXyIsolatorTraits > +class Construct_nk_decomposition_2 { + public: + + //! this instance's first template parameter + typedef SurfaceZAtXyIsolatorTraits Surface_z_at_xy_isolator_traits; + + SoX_SURFACE_Z_AT_XY_ISOLATOR_TRAITS_SNAP_TYPEDEFS( + Surface_z_at_xy_isolator_traits + ); + + typedef SoX::Restricted_cad_3< Surface_z_at_xy_isolator_traits > + Restricted_cad_3; + typedef SoX::Restricted_cad_3_accessor< Restricted_cad_3 > + Accessor; + + /*!\brief + * returns through \c oi projected + * intersection curve(s) of \c s1 and \c s2 + */ + Restricted_cad_3 operator()(const Surface_3& surface) const { + Restricted_cad_3 nkd; + + typename Surface_z_at_xy_isolator_traits:: + Construct_projected_surface_curves_2 + construct_projected_surface_curves + // TODO (traits.construct_projected_surface_curves_2_object()) + ; + + // STEP 1: Construct Silhouette-Arr with mults + +#if !NDEBUG + std::cout << "Construct A(S) for surface: " + << surface.id() << " ... " << std::flush; +#endif + + std::list< std::pair< Curve_analysis_2, int > > curves; + + construct_projected_surface_curves( + surface, + std::back_inserter(curves) + ); +#if !NDEBUG + std::cout << "#sil-curves: " << curves.size() << std::endl; + + for (typename std::list< std::pair< Curve_analysis_2, int > >:: + const_iterator it = curves.begin(); it != curves.end(); + it++) { + std::cout << "sil-curve: " + << it->first.polynomial_2() << std::endl; + } +#endif + Accessor acc(nkd); + + if (curves.empty()) { + // empty arrangement + acc.init(surface); + for (typename Restricted_cad_3::Face_const_iterator fit = + nkd.faces_begin(); + fit != nkd.faces_end(); fit++) { + nkd._set_nk_value(fit, surface, + SoX::Nk::N, surface.f().degree()); + nkd._set_nk_value(fit, surface, + SoX::Nk::K, 0); + } + return nkd; + } + + std::list< Restricted_cad_3 > nkds; + + // otherwise create for each part of silhouette curve + for (typename std::list< std::pair< Curve_analysis_2, int > >:: + const_iterator + cit = curves.begin(); cit != curves.end(); cit++) { + // split curve and insert points/segments into the arrangement + nkds.push_back( + _arrangement_of_curve( + surface, cit->first, + SoX::Nk::MULT, cit->second + ) + ); + } + + curves.clear(); + + Overlay_restricted_cad_3< Surface_z_at_xy_isolator_traits > + overlay; // TODO single instance + + nkd = overlay(nkds.begin(), nkds.end(), surface, SoX::Nk::MULT); + +#if !NDEBUG + std::cout << "done." << std::endl; +#endif + + CGAL_assertion_code( + unsigned int num_faces = nkd.number_of_faces(); + unsigned int num_unb_faces = nkd.number_of_unbounded_faces(); + ); + + CGAL_assertion_code(( + { + for (typename Restricted_cad_3::Vertex_const_iterator vit = + nkd.vertices_begin(); + vit != nkd.vertices_end(); vit++) { + CGAL_assertion(nkd.nk(vit, surface).mult() != -1); + } + for (typename Restricted_cad_3::Edge_const_iterator eit = + nkd.edges_begin(); + eit != nkd.edges_end(); eit++) { + CGAL_assertion(nkd.nk(eit, surface).mult() != -1); + } + for (typename Restricted_cad_3::Face_const_iterator fit = + nkd.faces_begin(); + fit != nkd.faces_end(); fit++) { + CGAL_assertion(nkd.nk(fit, surface).mult() == -1); + } + }) + ); + + // STEP 1-2: Set n and k for faces + for (typename Restricted_cad_3::Face_const_iterator fit = + nkd.faces_begin(); + fit != nkd.faces_end(); fit++) { + nkd._set_nk_value(fit, surface, + SoX::Nk::N, surface.f().degree()); + nkd._set_nk_value(fit, surface, + SoX::Nk::K, 0); + } + + CGAL_assertion_code(( + { + for (typename Restricted_cad_3::Vertex_const_iterator vit = + nkd.vertices_begin(); + vit != nkd.vertices_end(); vit++) { + CGAL_assertion(nkd.nk(vit, surface).mult() != -1); + } + for (typename Restricted_cad_3::Edge_const_iterator eit = + nkd.edges_begin(); + eit != nkd.edges_end(); eit++) { + CGAL_assertion(nkd.nk(eit,surface).mult() != -1); + } + for (typename Restricted_cad_3::Face_const_iterator fit = + nkd.faces_begin(); + fit != nkd.faces_end(); fit++) { + CGAL_assertion(nkd.nk(fit,surface).mult() == -1); + } + }) + ); + + // STEP 2: Refine edges of nkd with respect to a_i + + std::list< int > used_i; + + int max_i = surface.f().degree(); // == degree_z + + bool stop = false; + + // take care of a_i! + for (int i = max_i; i >= 0; i--) { + +#if !NDEBUG + std::cout << "Refine A(S) wrt a_" + << i + << " ... " << std::flush; + bool skip = false; +#endif + + if (surface.f(i).degree() != i) { +#if !NDEBUG + std::cout << "skipped." << std::endl; + skip = true; +#else + continue; +#endif + } + + curves.clear(); + + construct_projected_surface_curves( + surface, i, + std::back_inserter(curves) + ); + + bool whole_plane = + construct_projected_surface_curves(surface, i); + +#if !NDEBUG + std::cout << "#alpha-curves: " << curves.size() << std::endl; + for (typename std::list< std::pair< Curve_analysis_2, int > >:: + const_iterator it = curves.begin(); it != curves.end(); + it++) { + + std::cout << "alpha-curve: " + << it->first.polynomial_2() << std::endl; + } + if (skip) { + continue; + } +#endif + bool set_i = false; + stop = _refine_with_respect_to( + surface, nkd, curves.begin(), curves.end(), + whole_plane, + SoX::Nk::N, set_i, i + ); + + if (set_i) { + // store set i + used_i.push_back(i); + } +#if !NDEBUG + std::cout << "done." << std::endl; +#endif + if (stop) { +#if !NDEBUG + std::cout << "All n set." << std::endl; +#endif + break; + } + } + + if (!stop) { + for (typename Restricted_cad_3::Vertex_const_iterator vit = + nkd.vertices_begin(); + vit != nkd.vertices_end(); vit++) { + const SoX::Nk& nk = nkd.nk(vit, surface); + if (nk.n() == -2) { + nkd._set_nk_value(vit, surface, SoX::Nk::N, -1); + } + } + for (typename Restricted_cad_3::Edge_const_iterator eit = + nkd.edges_begin(); + eit != nkd.edges_end(); eit++) { + const SoX::Nk& nk = nkd.nk(eit, surface); + if (nk.n() == -2) { + nkd._set_nk_value(eit, surface, SoX::Nk::N, -1); + } + } + } + + CGAL_assertion(nkd.number_of_faces() == num_faces); + CGAL_assertion(nkd.number_of_unbounded_faces() == num_unb_faces); + + CGAL_assertion_code(( + { + for (typename Restricted_cad_3::Vertex_const_iterator vit = + nkd.vertices_begin(); + vit != nkd.vertices_end(); vit++) { + const SoX::Nk& nk = nkd.nk(vit, surface); + CGAL_assertion(nk.mult() != -1); + CGAL_assertion(nk.n() != -2); + } + for (typename Restricted_cad_3::Edge_const_iterator eit = + nkd.edges_begin(); + eit != nkd.edges_end(); eit++) { + const SoX::Nk& nk = nkd.nk(eit, surface); + CGAL_assertion(nk.mult() != -1); + CGAL_assertion(nk.n() != -2); + } + for (typename Restricted_cad_3::Face_const_iterator fit = + nkd.faces_begin(); + fit != nkd.faces_end(); fit++) { + const SoX::Nk& nk = nkd.nk(fit, surface); + CGAL_assertion(nk.mult() == -1); + CGAL_assertion(nk.n() != -2); + } + }) + ); + + // STEP 3: Refine Sil-Arr with respect to stha_i^n + // iterate over set a_i + stop = false; + + + for (typename std::list< int >::const_iterator nit = + used_i.begin(); nit != used_i.end(); nit++) { + + +#if !NDEBUG + std::cout << "*nit: " << *nit << std::endl; + std::cout << "f(n): " << surface.f(*nit).degree() << std::endl; +#endif + + // start with 0, as this can happen for asymptotic faces + for (int k = 0; k < surface.f(*nit).degree(); k++) { + +#if !NDEBUG + std::cout << "Refine A(S) wrt stha_" + << k << "^F(" << *nit << ")" + << " ... " << std::flush; + bool skip = false; +#endif + + // quick exit for silhouette itself + if (k == 0 && + *nit == typename CGAL::Polynomial_traits_d< Polynomial_3 >:: + Total_degree()(surface.f())) { +#if !NDEBUG + std::cout << "skipped." << std::endl; + skip = true; +#else + continue; +#endif + } + + curves.clear(); + + construct_projected_surface_curves( + surface, k, *nit, + std::back_inserter(curves) + ); + + bool whole_plane = + construct_projected_surface_curves(surface, k, *nit); + +#if !NDEBUG + std::cout << "#sigma-curves_" << k + << ": " << curves.size() << std::endl; + for (typename std::list< std::pair< Curve_analysis_2, int > >:: + const_iterator it = curves.begin(); it != curves.end(); + it++) { + + std::cout << "sigma-curve: " + << it->first.polynomial_2() << std::endl; + } + if (skip) { + continue; + } +#endif + + bool set_k = false; + stop = _refine_with_respect_to( + surface, nkd, curves.begin(), curves.end(), + whole_plane, + SoX::Nk::K, set_k, k, *nit + ); + +#if !NDEBUG + std::cout << "done." << std::endl; +#endif + + if (stop) { +#if !NDEBUG + std::cout << "All k set." << std::endl; +#endif + break; + } + } + } + + CGAL_assertion(nkd.number_of_faces() == num_faces); + CGAL_assertion(nkd.number_of_unbounded_faces() == num_unb_faces); + + if (!stop) { + for (typename Restricted_cad_3::Vertex_const_iterator vit = + nkd.vertices_begin(); + vit != nkd.vertices_end(); vit++) { + const SoX::Nk& nk = nkd.nk(vit, surface); + if (nk.k() == -2) { + if (nk.n() <= 0) { + CGAL_assertion(nk.n() > -2); + nkd._set_nk_value(vit, surface, SoX::Nk::K, nk.n()); + } + } + } + for (typename Restricted_cad_3::Edge_const_iterator eit = + nkd.edges_begin(); + eit != nkd.edges_end(); eit++) { + const SoX::Nk& nk = nkd.nk(eit, surface); + if (nk.k() == -2) { + if (nk.n() <= 0) { + CGAL_assertion(nk.n() > -2); + nkd._set_nk_value(eit, surface, SoX::Nk::K, nk.n()); + } + } + } + } + + // STEP 4: Delete unneccessary vertices + + for (typename Restricted_cad_3::Vertex_const_iterator vt, vit = + nkd.vertices_begin(); + vit != nkd.vertices_end();) { + + bool remove = false; + + if (vit->degree() == 2) { + const SoX::Nk& nk = nkd.nk(vit, surface); + //std::cout << "vnk: " << nk << std::endl; + + + typename + Restricted_cad_3::Halfedge_around_vertex_const_circulator + e1 = vit->incident_halfedges(); + + const SoX::Nk& nk1 = nkd.nk(e1, surface); + + if (nk1.same_nk(nk)) { + + typename + Restricted_cad_3:: + Halfedge_around_vertex_const_circulator e2 = + e1->next(); + CGAL_assertion(e1 != e2); + + const SoX::Nk& nk2 = nkd.nk(e2, surface); + + if (nk2.same_nk(nk)) { + if (e2->curve().are_mergeable(e1->curve())) { + remove = true; + } + } + } + } + + if (remove) { + vt = vit; + vit++; + acc.remove_vertex(acc.non_const_handle(vt)); + } else { + vit++; + } + } + + // FINISHED! + CGAL_assertion_code(( + { + for (typename Restricted_cad_3::Face_const_iterator fit = + nkd.faces_begin(); + fit != nkd.faces_end(); fit++) { + const SoX::Nk& nk = nkd.nk(fit, surface); + //std::cout << "fnk: " << nk << std::endl; + CGAL_assertion(nk.mult() == -1); + CGAL_assertion(nk.n() == surface.f().degree()); + CGAL_assertion(nk.k() == 0); + } + + for (typename Restricted_cad_3::Edge_const_iterator + eit = + nkd.edges_begin(); + eit != nkd.edges_end(); eit++) { + const SoX::Nk& nk = nkd.nk(eit, surface); + //std::cout << "hnk: " << nk << std::endl; + CGAL_assertion(nk.mult() > 0); + CGAL_assertion(nk.n() != -2); + CGAL_assertion(nk.k() != -2); + CGAL_assertion( + (nk.n() == -1 && nk.k() == -1) || + (nk.n() == 0 && nk.k() == 0) || + (nk.n() > 0 && nk.k() >= 0 && nk.k() < nk.n()) + ); + } + + for (typename Restricted_cad_3::Vertex_const_iterator vit = + nkd.vertices_begin(); + vit != nkd.vertices_end(); vit++) { + const SoX::Nk& nk = nkd.nk(vit, surface); + //std::cout << "vnk: " << nk << std::endl; + CGAL_assertion(nk.mult() == 0); + CGAL_assertion(nk.n() != -2); + CGAL_assertion(nk.k() != -2); + CGAL_assertion( + (nk.n() == 0 && nk.k() == 0) || + (nk.n() > 0 && nk.k() >= 0 && nk.k() < nk.n()) || + (nk.n() == -1 && nk.k() == -1) + ); + } + }) + ); + + return nkd; + } + +private: + + /*!\brief + * creates segments and points for a given curve and inserts them + * into an arrangement + */ + Restricted_cad_3 + _arrangement_of_curve(const Surface_3& surface, + const Curve_analysis_2& curve, + SoX::Nk::Value_type type, + int value) const { + + Restricted_cad_3 nkd = Accessor::construct_for_curve(curve); + + Accessor acc(nkd); + + // construct data objects + acc.init(surface); + + // set nk-data + if (type == SoX::Nk::MULT) { + // set mult + // TODO store multiplicity of singular point??!?! + for (typename Restricted_cad_3::Vertex_const_iterator vit = + nkd.vertices_begin(); + vit != nkd.vertices_end(); vit++) { + nkd._set_nk_value(vit, surface, SoX::Nk::MULT, 0); + // set it to 0 to indicate + // that it belongs to sil + } + for (typename Restricted_cad_3::Edge_const_iterator eit = + nkd.edges_begin(); + eit != nkd.edges_end(); eit++) { + nkd._set_nk_value(eit, surface, SoX::Nk::MULT, value); + } + for (typename Restricted_cad_3::Face_const_iterator fit = + nkd.faces_begin(); + fit != nkd.faces_end(); fit++) { + CGAL_assertion(nkd.nk(fit, surface).mult() == -1); + } + } + + CGAL_assertion_code(( + { + for (typename Restricted_cad_3::Vertex_const_iterator vit = + nkd.vertices_begin(); + vit != nkd.vertices_end(); vit++) { + CGAL_postcondition(vit->data()); + } + for (typename Restricted_cad_3::Edge_const_iterator eit = + nkd.edges_begin(); + eit != nkd.edges_end(); eit++) { + CGAL_postcondition(eit->data()); + CGAL_postcondition(eit->twin()->data()); + } + for (typename Restricted_cad_3::Face_const_iterator fit = + nkd.faces_begin(); + fit != nkd.faces_end(); fit++) { + CGAL_postcondition(fit->data()); + } + }) + ); + + return nkd; + } + + template < class InputIterator > + bool _refine_with_respect_to(const Surface_3& surface, + Restricted_cad_3& nkd, + InputIterator begin, InputIterator end, + bool whole_plane, + SoX::Nk::Value_type type, + bool& set_value, + int value1, + int value2 = -1) const { + + Overlay_restricted_cad_3< Surface_z_at_xy_isolator_traits > + overlay; // TODO single instance + + Restricted_cad_3 nkrefine; + + set_value = false; + + if (begin == end) { + + Accessor acc(nkrefine); + acc.init(surface); + + } else { + + CGAL_precondition(!whole_plane); + + std::list< Restricted_cad_3 > nkds; + + for (InputIterator it = begin; it != end; it++) { + // split curve and insert points/segments + // into the arrangement + nkds.push_back( + _arrangement_of_curve( + surface, it->first, type, -1 + ) + ); + } + + nkrefine = overlay(nkds.begin(), nkds.end(), surface, type); + } + + // overlay nkd with nkrefine + Restricted_cad_3 tmp = overlay(nkd, nkrefine, surface, type); + + Accessor acc(tmp); + + if (begin != end) { + + // clean wrong halfedges within a face of nkd + for (typename Restricted_cad_3::Edge_const_iterator + ti, eit = + tmp.edges_begin(); + eit != tmp.edges_end();) { + if (tmp.nk(eit, surface).feature1() == SoX::FACE) { + ti = eit; + eit++; + acc.remove_edge(acc.non_const_handle(ti)); + } else { + eit++; + } + } + + // clean wrong vertices within a face of nkd + for (typename Restricted_cad_3::Vertex_const_iterator + vt, vit = + tmp.vertices_begin(); + vit != tmp.vertices_end();) { + if (tmp.nk(vit, surface).feature1() == SoX::FACE) { + vt = vit; + vit++; + acc.remove_vertex(acc.non_const_handle(vt)); + } else { + vit++; + } + } + + // clean vertices with already stored data + for (typename Restricted_cad_3::Vertex_const_iterator + vt, vit = + tmp.vertices_begin(); + vit != tmp.vertices_end(); ) { + const SoX::Nk& nk = tmp.nk(vit, surface); + bool remove = false; + if (nk.feature1() == SoX::EDGE) { + if (type == SoX::Nk::N) { + if (nk.n() == -2) { + remove = true; + } + } else if (type == SoX::Nk::K) { + if (nk.k() != -2 || + nk.n() != value2) { + remove = true; + } + } + } + if (remove) { + vt = vit; + vit++; + acc.remove_vertex(acc.non_const_handle(vt)); + } else { + vit++; + } + } + } + + bool stop = true; + // set values for correct halfedges ... + for (typename Restricted_cad_3::Edge_const_iterator + eit = + tmp.edges_begin(); + eit != tmp.edges_end(); eit++) { + SoX::Nk nk = tmp.nk(eit, surface); + //std::cout << "H1: " << nk << std::endl; + if (nk.feature2() == SoX::FACE) { + CGAL_assertion(nk.mult() != -1); + if (type == SoX::Nk::N) { + if (!whole_plane && + nk.n() == -2) { + //std::cout << "HN: " << value1 << std::endl; + nkd._set_nk_value(eit, surface, SoX::Nk::N, value1); + set_value = true; + } + } else if (type == SoX::Nk::K) { + if (nk.n() == value2 && + !whole_plane && + nk.k() == -2) { + //std::cout << "HK: " << value1 << std::endl; + nkd._set_nk_value(eit, surface, SoX::Nk::K, value1); + set_value = true; + } + } + } + nk = tmp.nk(eit, surface); + if (stop) { + if (type == SoX::Nk::N) { + if (nk.n() == -2) { + stop = false; + } + } else if (type == SoX::Nk::K) { + if (nk.k() == -2) { + stop = false; + } + } + } + } + + // ... and vertices + for (typename Restricted_cad_3::Vertex_const_iterator vit = + tmp.vertices_begin(); + vit != tmp.vertices_end(); vit++) { + SoX::Nk nk = tmp.nk(vit, surface); + //std::cout << "V1: " << nk << std::endl; + if (nk.feature2() == SoX::FACE) { + CGAL_assertion(nk.mult() != -1); + if (type == SoX::Nk::N) { + if (!whole_plane && + nk.n() == -2) { + //std::cout << "VN: " << value1 << std::endl; + nkd._set_nk_value(vit, surface, SoX::Nk::N, value1); + set_value = true; + } + } else if (type == SoX::Nk::K) { + if (!whole_plane && + nk.n() == value2 && + nk.k() == -2) { + //std::cout << "VK: " << value1 << std::endl; + nkd._set_nk_value(vit, surface, SoX::Nk::K, value1); + set_value = true; + } + } + } + nk = tmp.nk(vit, surface); + if (stop) { + if (type == SoX::Nk::N) { + if (nk.n() == -2) { + stop = false; + } + } else if (type == SoX::Nk::K) { + if (nk.k() == -2) { + stop = false; + } + } + } + } + + nkd = tmp; + + return stop; + } +}; + + + +} // namespace Intern + +/*!\brief + * Functor class that provides creation of restriced_cads. + */ +template < class SurfaceZAtXyIsolatorTraits > +class Create_restricted_cad_3 { +public: + + //! this instance's first template parameter + typedef SurfaceZAtXyIsolatorTraits Surface_z_at_xy_isolator_traits; + + SoX_SURFACE_Z_AT_XY_ISOLATOR_TRAITS_SNAP_TYPEDEFS( + Surface_z_at_xy_isolator_traits + ); + + //! type of instantiated class + typedef Create_restricted_cad_3< Surface_z_at_xy_isolator_traits > Self; + +#if DOXYGEN_RUNNING + //! type of projected curve + typedef typename Surface_z_at_xy_isolator_traits::Curve_analysis_2 + Curve_analysis_2; + + //! type of projected point + typedef typename Surface_z_at_xy_isolator_traits::Point_2 Point_2; + + //! type of projected segment + typedef typename Surface_z_at_xy_isolator_traits::X_monotone_curve_2 + X_monotone_curve_2; +#endif + + //! type of Restricted_cad_3 + typedef Restricted_cad_3< Surface_z_at_xy_isolator_traits > + Restricted_cad_3; + +private: + //! type of data + typedef typename Restricted_cad_3::Data Data; + + +public: + /*!\brief + * creates cad for silhouette-curve of \c surface + */ + Restricted_cad_3 operator()(const Surface_3& surface) const { + return Restricted_cad_3::cad_cache()(surface); + } + +private: + //! used to indicate special constructors + class Non_cached_construction_tag {}; + + /*!\brief + * creates cad for silhouette-curve of \c surface (called in cache only) + */ + Restricted_cad_3 operator()(const Surface_3& surface, + Non_cached_construction_tag t) const { + + typedef Intern::Construct_nk_decomposition_2< + Surface_z_at_xy_isolator_traits > + Construct_nk_decomposition_2; + + Construct_nk_decomposition_2 construct_nk_decomposition + // TODO (traits.construct_nk_decomposition_2_object()) + ; + + Restricted_cad_3 nkd = construct_nk_decomposition(surface); + + nkd._finalize(surface); + + // and return it + return nkd; + } + +public: + /*!\brief + * creates cad for cut-curve of \c surface1 and \c surface2 + */ + Restricted_cad_3 operator()(const Surface_3& surface1, + const Surface_3& surface2) const { + return Restricted_cad_3::cad_cache()(surface1, surface2); + } + +private: + /*!\brief + * creates cad for cut-curve of \c surface1 and \c surface2 + * (called in cache only) + */ + Restricted_cad_3 operator()(const Surface_3& surface1, + const Surface_3& surface2, + Non_cached_construction_tag t) const { + + typename + Surface_z_at_xy_isolator_traits:: + Construct_projected_cuts_2 + construct_projected_cuts + // TODO (traits.construct_projected_cuts_2_object()) + ; + + // store components of curve + std::list< std::pair< Curve_analysis_2, int > > curves; + + construct_projected_cuts( + surface1, surface2, + std::back_inserter(curves) + ); + + if (curves.empty()) { + // empty arrangement + Restricted_cad_3 cad; + cad._renew(); + + // for all faces + for (typename Restricted_cad_3::Face_const_iterator fit = + cad.faces_begin(); + fit != cad.faces_end(); fit++) { + if (!fit->data()) { + CGAL::Object obj = CGAL::make_object(fit); + cad.non_const_handle(fit)->data() = Data(SoX::FACE, obj); + fit->data()->_set_rs_id(cad.id()); + fit->data()->_add_surface(surface1); + fit->data()->_add_surface(surface2); + fit->data()->_set_dcel(surface1, surface2, SoX::FACE, obj); + } + } + // and return it + return cad; + } + + // for each part create RS_CAD_3 + std::vector< Restricted_cad_3 > + cads(std::distance(curves.begin(), curves.end()), + Restricted_cad_3() + ); + int i = 0; + + typedef SoX::Restricted_cad_3_accessor< Restricted_cad_3 > + Accessor; + + for (typename std::list< std::pair< Curve_analysis_2, int > >:: + const_iterator + cit = curves.begin(); cit != curves.end(); cit++, i++) { + + // make rep unique + cads[i]._renew(); + + // split curve and insert points/segments into the arrangement + cads[i] = Accessor::construct_for_curve(cit->first); + + // assign values + _set_cut(surface1, surface2, cit->second, cads[i]); + } + + // merge the RS_CAD_3s + Overlay_restricted_cad_3< Surface_z_at_xy_isolator_traits > + overlay; // TODO single instance + + CGAL_assertion(cads.begin() != cads.end()); + Restricted_cad_3 tmp = overlay(cads.begin(), cads.end(), true); + + // it remains to set handles correctly + for (typename Restricted_cad_3::Vertex_const_iterator vit = + tmp.vertices_begin(); + vit != tmp.vertices_end(); vit++) { + vit->data()->_set_rs_id(tmp.id()); + CGAL::Object obj = CGAL::make_object(vit); + vit->data()->_set_dcel(SoX::VERTEX, obj); + vit->data()->_set_dcel(surface1, surface2, SoX::VERTEX, obj); + } + for (typename Restricted_cad_3::Halfedge_const_iterator hit = + tmp.halfedges_begin(); + hit != tmp.halfedges_end(); hit++) { + bool correct_dir = (hit->direction() == CGAL::ARR_LEFT_TO_RIGHT); + if (correct_dir) { + hit->data()->_set_rs_id(tmp.id()); + CGAL::Object obj = CGAL::make_object(hit); + hit->data()->_set_dcel(SoX::EDGE, obj); + hit->data()->_set_dcel(surface1, surface2, SoX::EDGE, obj); + } + } + for (typename Restricted_cad_3::Face_const_iterator fit = + tmp.faces_begin(); + fit != tmp.faces_end(); fit++) { + fit->data()->_set_rs_id(tmp.id()); + CGAL::Object obj = CGAL::make_object(fit); + fit->data()->_set_dcel(SoX::FACE, obj); + fit->data()->_set_dcel(surface1, surface2, SoX::FACE, obj); +#if 0 + // TASK benchmark: check whether this saves time + // it may be even problematic, as we have to isolate + // over non-face-features of the "virtual" silhouette-curve + cad.slice(fit); +#endif + } + + // and return it + return tmp; + } + +public: + /*!\brief + * Compute cad for a range of surfaces defined by \c [begin,end) + */ + // TASK move to Surface_cad?? + template < class InputIterator > + Restricted_cad_3 operator()(InputIterator begin, InputIterator end) const { + + CGAL_precondition((boost::is_same< typename InputIterator::value_type, + Surface_3 >::value + )); + + std::vector< Restricted_cad_3 > boundaries; + boundaries.reserve(std::distance(begin,end)); + for (InputIterator it = begin; it != end; it++) { + boundaries.push_back(this->operator()(*it)); + } + + std::vector< Restricted_cad_3 > cuts; + int d = std::distance(begin,end); + cuts.reserve(d*(d-1)); + for (InputIterator oit = begin; oit != end; oit++) { + for (InputIterator iit = oit; iit != end; iit++) { + if (oit != iit) { + cuts.push_back(this->operator()(*oit,*iit)); + } + } + } + + Overlay_restricted_cad_3< Surface_z_at_xy_isolator_traits > + overlay; // TODO single instance + + // TASK stream +#if NDEBUG + std::cout << "Boundaries ... " << std::flush; +#endif + Restricted_cad_3 bdry_overlay = overlay(boundaries.begin(), + boundaries.end()); +#if NDEBUG + std::cout << "done." << std::endl; +#endif + Restricted_cad_3 out; + out._renew(); + + if (static_cast< int >(boundaries.size()) > 1) { +#if NDEBUG + std::cout << "Cuts ... " << std::flush; +#endif + Restricted_cad_3 itxs_overlay = overlay(cuts.begin(), + cuts.end()); +#if NDEBUG + std::cout << "done." << std::endl; +#endif + +#if NDEBUG + std::cout << "Final ... " << std::flush; +#endif + out = overlay(bdry_overlay, itxs_overlay); +#if NDEBUG + std::cout << "done." << std::endl; +#endif + } else { + out = bdry_overlay; +#if NDEBUG + std::cout << "Final consists of a single silhouette ... done." + << std::endl; +#endif + } + + return out; + } + +private: + + /*!\brief + * sets features of dcel with silhouette of \c surface + */ + void _set_silhouette(const Surface_3& surface, + int multiplicity, bool vertical, + Restricted_cad_3& rscad) const { + + // for all vertices + for (typename Restricted_cad_3::Vertex_const_iterator vit = + rscad.vertices_begin(); + vit != rscad.vertices_end(); vit++) { + CGAL_assertion(!vit->is_at_infinity()); + if (!vit->data()) { + rscad.non_const_handle(vit)->data() = + Data(SoX::VERTEX, CGAL::make_object(vit)); + vit->data()->_add_surface(surface); + } + // store multiplicity of singular point +#if 0 + if (vit->point().curve().may_be_singular()) { + int vmult = vit->point().curve().event_info_at_x( + vit->point().x().number().finite() + ).multiplicity( + vit->point().arcno() + ); + vit->data()->_add_silhouette(surface, multiplicity * vmult, + false, vertical); + } else { + vit->data()->_add_silhouette(surface, multiplicity, + false, vertical); + } +#else + vit->data()->_add_silhouette(surface, -1, + false, vertical); +#endif + } + + // for all edges + for (typename Restricted_cad_3::Halfedge_const_iterator hit = + rscad.halfedges_begin(); + hit != rscad.halfedges_end(); hit++) { + if (!hit->data()) { + rscad.non_const_handle(hit)->data() = + Data(SoX::EDGE, CGAL::make_object(hit)); + hit->data()->_add_silhouette(surface, multiplicity, + false, vertical); + hit->data()->_add_surface(surface); + } + if (!hit->twin()->data()) { + rscad.non_const_handle(hit)->twin()->data() = *hit->data(); + } + } + + // for all faces + for (typename Restricted_cad_3::Face_const_iterator fit = + rscad.faces_begin(); + fit != rscad.faces_end(); fit++) { + if (!fit->data()) { + rscad.non_const_handle(fit)->data() = + Data(SoX::FACE, CGAL::make_object(fit)); + fit->data()->_add_surface(surface); + } + } + } + + /*!\brief + * sets features of dcel with cut of \c surface1 and \c surface2 + */ + void _set_cut(const Surface_3& surface1, + const Surface_3& surface2, + int multiplicity, + Restricted_cad_3& rscad) const { + + // for all vertices + for (typename Restricted_cad_3::Vertex_const_iterator vit = + rscad.vertices_begin(); + vit != rscad.vertices_end(); vit++) { + CGAL_assertion(!vit->is_at_infinity()); + if (!vit->data()) { + rscad.non_const_handle(vit)->data() = + Data(SoX::VERTEX, CGAL::make_object(vit)); + vit->data()->_add_surface(surface1); + vit->data()->_add_surface(surface2); + } + // store multiplicity of singular point +#if 0 + if (vit->point().curve().may_be_singular()) { + int vmult = vit->point().curve().event_info_at_x( + vit->point().x().number().finite() + ).multiplicity( + vit->point().arcno() + ); + vit->data()->_add_cut( + surface1, surface2, multiplicity * vmult, false + ); + } else { + vit->data()->_add_cut( + surface1, surface2, multiplicity, false + ); + } +#else + vit->data()->_add_cut( + surface1, surface2, -1, false + ); +#endif + } + + // for all edges + for (typename Restricted_cad_3::Halfedge_const_iterator hit = + rscad.halfedges_begin(); + hit != rscad.halfedges_end(); hit++) { + if (!hit->data()) { + rscad.non_const_handle(hit)->data() = + Data(SoX::EDGE, CGAL::make_object(hit)); + hit->data()->_add_cut(surface1, surface2, + multiplicity, false); + hit->data()->_add_surface(surface1); + hit->data()->_add_surface(surface2); + } + if (!hit->twin()->data()) { + rscad.non_const_handle(hit)->twin()->data() = *hit->data(); + } + } + + // for all faces + for (typename Restricted_cad_3::Face_const_iterator fit = + rscad.faces_begin(); + fit != rscad.faces_end(); fit++) { + if (!fit->data()) { + rscad.non_const_handle(fit)->data() = + Data(SoX::FACE, CGAL::make_object(fit)); + fit->data()->_add_surface(surface1); + fit->data()->_add_surface(surface2); + } + } + } + +private: + + // friend + // for non-cached construction + friend class + Intern::Restricted_cad_3_cache< Surface_z_at_xy_isolator_traits >; + +}; // Create + + +// TODO move Overlay_restricted_cad_3 to namespace Intern?!?! +/*!\brief + * Functor class that provides overlaying of restriced_cads. + */ +template < class SurfaceZAtXyIsolatorTraits > +class Overlay_restricted_cad_3 { +public: + + //! this instance's first template parameter + typedef SurfaceZAtXyIsolatorTraits Surface_z_at_xy_isolator_traits; + + SoX_SURFACE_Z_AT_XY_ISOLATOR_TRAITS_SNAP_TYPEDEFS( + Surface_z_at_xy_isolator_traits + ); + + //! type of instantiated class + typedef Overlay_restricted_cad_3< Surface_z_at_xy_isolator_traits > Self; + + //! type of Restricted_cad_3 + typedef Restricted_cad_3< Surface_z_at_xy_isolator_traits > + Restricted_cad_3; + + /*!\brief + * merges two restricted cads + */ + Restricted_cad_3 operator()(const Restricted_cad_3& rscad1, + const Restricted_cad_3& rscad2) const { + return this->operator()(rscad1, rscad2, false); + } + +private: + /*!\brief + * merges two restricted cads + */ + Restricted_cad_3 operator()(const Restricted_cad_3& rscad1, + const Restricted_cad_3& rscad2, + bool factors_of_same_curve) const { + + Restricted_cad_3 out; + out._renew(); + + CGAL_precondition_code( + for (typename Restricted_cad_3::Vertex_const_iterator vit = + rscad1.vertices_begin(); + vit != rscad1.vertices_end(); vit++) { + CGAL_precondition(vit->data()); + } + for (typename Restricted_cad_3::Halfedge_const_iterator hit = + rscad1.halfedges_begin(); + hit != rscad1.halfedges_end(); hit++) { + CGAL_precondition(hit->data()); + } + for (typename Restricted_cad_3::Face_const_iterator fit = + rscad1.faces_begin(); + fit != rscad1.faces_end(); fit++) { + CGAL_precondition(fit->data()); + } + for (typename Restricted_cad_3::Vertex_const_iterator vit = + rscad2.vertices_begin(); + vit != rscad2.vertices_end(); vit++) { + CGAL_precondition(vit->data()); + } + for (typename Restricted_cad_3::Halfedge_const_iterator hit = + rscad2.halfedges_begin(); + hit != rscad2.halfedges_end(); hit++) { + CGAL_precondition(hit->data()); + } + for (typename Restricted_cad_3::Face_const_iterator fit = + rscad2.faces_begin(); + fit != rscad2.faces_end(); fit++) { + CGAL_precondition(fit->data()); + } + ); + + // overlay computation (unique place to do it) + out._overlay(rscad1, rscad2, factors_of_same_curve); + + // TASK benchmark: + // compute slices already here vs. delay their computation + + for (typename Restricted_cad_3::Vertex_const_iterator vit = + out.vertices_begin(); + vit != out.vertices_end(); vit++) { + vit->data()->_set_rs_id(out.id()); + CGAL::Object obj = CGAL::make_object(vit); + vit->data()->_set_dcel(SoX::VERTEX, obj); + CGAL_postcondition(vit->data()); + } + for (typename Restricted_cad_3::Halfedge_const_iterator hit = + out.halfedges_begin(); + hit != out.halfedges_end(); hit++) { + bool correct_dir = (hit->direction() == CGAL::ARR_LEFT_TO_RIGHT); + if (correct_dir) { + hit->data()->_set_rs_id(out.id()); + CGAL::Object obj = CGAL::make_object(hit); + hit->data()->_set_dcel(SoX::EDGE, obj); + } + CGAL_postcondition(hit->data()); + } + for (typename Restricted_cad_3::Face_const_iterator fit = + out.faces_begin(); + fit != out.faces_end(); fit++) { + fit->data()->_set_rs_id(out.id()); + CGAL::Object obj = CGAL::make_object(fit); + fit->data()->_set_dcel(SoX::FACE, obj); + CGAL_postcondition(fit->data()); + } + + // cache id + return Restricted_cad_3::cad_cache()(out); + } + + + Restricted_cad_3 operator()(const Restricted_cad_3& rscad1, + const Restricted_cad_3& rscad2, + const Surface_3& surface, + SoX::Nk::Value_type type) const { + + Restricted_cad_3 out; + out._renew(); + + CGAL_precondition_code( + for (typename Restricted_cad_3::Vertex_const_iterator vit = + rscad1.vertices_begin(); + vit != rscad1.vertices_end(); vit++) { + CGAL_precondition(vit->data()); + } + for (typename Restricted_cad_3::Halfedge_const_iterator hit = + rscad1.halfedges_begin(); + hit != rscad1.halfedges_end(); hit++) { + CGAL_precondition(hit->data()); + } + for (typename Restricted_cad_3::Face_const_iterator fit = + rscad1.faces_begin(); + fit != rscad1.faces_end(); fit++) { + CGAL_precondition(fit->data()); + } + for (typename Restricted_cad_3::Vertex_const_iterator vit = + rscad2.vertices_begin(); + vit != rscad2.vertices_end(); vit++) { + CGAL_precondition(vit->data()); + } + for (typename Restricted_cad_3::Halfedge_const_iterator hit = + rscad2.halfedges_begin(); + hit != rscad2.halfedges_end(); hit++) { + CGAL_precondition(hit->data()); + } + for (typename Restricted_cad_3::Face_const_iterator fit = + rscad2.faces_begin(); + fit != rscad2.faces_end(); fit++) { + CGAL_precondition(fit->data()); + } + ); + + // overlay computation (unique place to do it) + out._overlay(rscad1, rscad2, surface, type); + + for (typename Restricted_cad_3::Vertex_const_iterator vit = + out.vertices_begin(); + vit != out.vertices_end(); vit++) { + CGAL_postcondition(vit->data()); + } + for (typename Restricted_cad_3::Halfedge_const_iterator hit = + out.halfedges_begin(); + hit != out.halfedges_end(); hit++) { + CGAL_postcondition(hit->data()); + } + for (typename Restricted_cad_3::Face_const_iterator fit = + out.faces_begin(); + fit != out.faces_end(); fit++) { + CGAL_postcondition(fit->data()); + } + + // cache id + return Restricted_cad_3::cad_cache()(out); + } + +public: + /*!\brief + * Computes overlay of a range of restricted cads + */ + template < class InputIterator > + Restricted_cad_3 operator()(InputIterator begin, InputIterator end) const { + return this->operator()(begin, end, false); + } + +private: + /*!\brief + * Computes overlay of a range of restricted cads + * + * \precond range contains at least one element + */ + template < class InputIterator > + Restricted_cad_3 operator()(InputIterator begin, InputIterator end, + bool factors_of_same_curve) const { + + CGAL_precondition(begin != end); + + CGAL_precondition((boost::is_same< typename InputIterator::value_type, + Restricted_cad_3 >::value + )); + + std::vector< Restricted_cad_3 > cads; + std::copy(begin, end, std::back_inserter(cads)); + + int n = static_cast< int >(cads.size()); + CGAL_assertion(n > 0); + while (n > 1) { + // TASK better random generator + std::random_shuffle(cads.begin(),cads.end()); + + // TASK stream +#if NDEBUG + std::cout << "Size " << n << std::flush; +#endif + std::vector< Restricted_cad_3 > tmp; + int i = 0; + for (; i + 1 < n; i += 2) { + tmp.push_back(this->operator()(cads[i],cads[i+1])); + } + // odd number of items + if (i + 1 == n) { + tmp.push_back(cads[i]); + } + cads.clear(); + cads = tmp; + n = static_cast< int >(cads.size()); + if (n > 1) { +#if NDEBUG + std::cout << "," << std::flush; +#endif + } +#if NDEBUG + std::cout << " " << std::flush; +#endif + } + CGAL_assertion(n == 1); + return cads[0]; + } + + + template < class InputIterator > + Restricted_cad_3 operator()(InputIterator begin, InputIterator end, + const Surface_3& surface, + SoX::Nk::Value_type type) const { + + CGAL_precondition(begin != end); + + CGAL_precondition((boost::is_same< typename InputIterator::value_type, + Restricted_cad_3 >::value + )); + + std::vector< Restricted_cad_3 > cads; + std::copy(begin, end, std::back_inserter(cads)); + + int n = static_cast< int >(cads.size()); + CGAL_assertion(n > 0); + while (n > 1) { + // TASK better random generator + std::random_shuffle(cads.begin(),cads.end()); + + // TASK stream +#if NDEBUG + std::cout << "Size " << n << std::flush; +#endif + std::vector< Restricted_cad_3 > tmp; + int i = 0; + for (; i + 1 < n; i += 2) { + tmp.push_back(this->operator()(cads[i], cads[i+1], + surface, type)); + } + // odd number of items + if (i + 1 == n) { + tmp.push_back(cads[i]); + } + cads.clear(); + cads = tmp; + n = static_cast< int >(cads.size()); + if (n > 1) { +#if NDEBUG + std::cout << "," << std::flush; +#endif + } +#if NDEBUG + std::cout << " " << std::flush; +#endif + } + CGAL_assertion(n == 1); + return cads[0]; + } + + +private: + + // friends + // for overlay with set factors_of_same_curve + friend class Create_restricted_cad_3< Surface_z_at_xy_isolator_traits >; + + friend class + Intern::Construct_nk_decomposition_2< Surface_z_at_xy_isolator_traits >; + +}; // Overlay + +} // namespace SoX + +#endif // SoX_GAPS_RESTRICTED_CAD_3_FUNCTORS_H diff --git a/Curved_kernel_via_analysis_2/include/CGAL/Arrangement_2l/Surface_3_envelope_traits.h b/Curved_kernel_via_analysis_2/include/CGAL/Arrangement_2l/Surface_3_envelope_traits.h new file mode 100644 index 00000000000..c8f34b1520c --- /dev/null +++ b/Curved_kernel_via_analysis_2/include/CGAL/Arrangement_2l/Surface_3_envelope_traits.h @@ -0,0 +1,841 @@ +// ============================================================================ +// +// Copyright (c) 2001-2008 Max-Planck-Institut Saarbruecken (Germany). +// All rights reserved. +// +// This file is part of EXACUS (http://www.mpi-inf.mpg.de/projects/EXACUS/); +// you may redistribute it under the terms of the Q Public License version 1.0. +// See the file LICENSE.QPL distributed with EXACUS. +// +// Licensees holding a valid commercial license may use this file in +// accordance with the commercial license agreement provided with the software. +// +// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE +// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. +// +// ---------------------------------------------------------------------------- +// +// Library : QdX +// File : include/SoX/GAPS/Surface_3_envelope_traits.h +// QdX_release : $Name: $ +// Revision : $Revision$ +// Revision_date : $Date$ +// +// Author(s) : Eric Berberich +// +// ============================================================================ + +#ifndef SoX_GAPS_SURFACE_3_ENVELOPE_TRAITS +#define SoX_GAPS_SURFACE_3_ENVELOPE_TRAITS 1 + +/*!\file SoX/GAPS/Surface_3_envelope_traits.h + * \brief Model for CGAL's EnvelopeTraits_3 concept. + */ + +#ifndef CGAL_ENVELOPE_3_USE_EDGE_HANDLES +#define CGAL_ENVELOPE_3_USE_EDGE_HANDLES 0 +#endif + +#include + +#include + +#include +#include +#include + +#include + +namespace SoX { + +namespace Intern { + +#if CGAL_ENVELOPE_3_USE_EDGE_HANDLES +template < class SurfacePair_3 > +struct Intersection_info : public CGAL::Handle_with_policy< +std::map< SurfacePair_3, + typename SurfacePair_3::Restricted_cad_3::Edge_const_handle, + CGAL::Handle_id_less_than< SurfacePair_3 > +> +> { + + typedef SurfacePair_3 Surface_pair_3; + + typedef typename Surface_pair_3::Restricted_cad_3::Edge_const_handle + Edge_const_handle; + + typedef CGAL::Handle_id_less_than< Surface_pair_3 > Less; + + typedef std::map< SurfacePair_3, Edge_const_handle, Less > Rep; + + typedef CGAL::Handle_with_policy< Rep > Base; + + typedef Intersection_info< Surface_pair_3 > Self; + + typedef typename Rep::const_iterator const_iterator; + + typedef typename Rep::value_type value_type; + +private: + typedef typename Rep::iterator iterator; + +public: + + Intersection_info() { + } + + Intersection_info(const Surface_pair_3& surface, + const Edge_const_handle& eh) { + CGAL_precondition(this->ptr()->empty()); + this->ptr()->insert(std::make_pair(surface, eh)); + } + + const_iterator begin() const { + return this->ptr()->begin(); + } + + const_iterator end() const { + return this->ptr()->end(); + } + + const_iterator find(const Surface_pair_3& pair) const { + return this->ptr()->find(pair); + } + + Self operator() (Self d1, Self d2) { + Self tmp = d2; + tmp.copy_on_write(); + + for (iterator it1 = d1.ptr()->begin(); + it1 != d1.ptr()->end(); it1++) { + iterator it2 = tmp.ptr()->find(it1->first); + if (it2 == tmp.ptr()->end()) { + tmp.ptr()->insert(it2, *it1); + } + } + return tmp; + } + + bool operator== (Self self) const { + return this->id() == self.id(); + } +}; +#endif + +} // namespace Intern + + +// TODO implement Apollonius-Mode +// TODO implicitly requires Surface_pair_3::surface_pair_cache() +template < class SurfacePair_3 > +class Surface_3_envelope_traits : public +#if CGAL_ENVELOPE_3_USE_EDGE_HANDLES +CGAL::Arr_curve_data_traits_2< +typename SurfacePair_3::Surface_z_at_xy_isolator_traits::Arrangement_traits_2, +Intern::Intersection_info< SurfacePair_3 >, +Intern::Intersection_info< SurfacePair_3 > +> +#else +SurfacePair_3::Surface_z_at_xy_isolator_traits::Arrangement_traits_2 +#endif +{ +public: + //! this instance template parameter + typedef SurfacePair_3 Surface_pair_3; + + //! the class itself + typedef Surface_3_envelope_traits< Surface_pair_3 > Self; + + // types for EnvelopeTraits + //! type of surfaces + typedef typename Surface_pair_3::Surface_3 Surface_3; + + //! type of Xy_monotone_surface_3 + typedef Surface_3 Xy_monotone_surface_3; + + //! type of multiplicity + typedef unsigned int Multiplicity; + +private: + //! type of Restricted_cad + typedef typename Surface_pair_3::Restricted_cad_3 Restricted_cad_3; + + //! type of surface traits + typedef typename Surface_pair_3::Surface_z_at_xy_isolator_traits + Surface_z_at_xy_isolator_traits; + + //! type of basic traits + typedef typename Surface_z_at_xy_isolator_traits::Arrangement_traits_2 + Arrangement_traits_2; + + //! type of Edge_const_handle + typedef typename Restricted_cad_3::Edge_const_handle Edge_const_handle; + +#if CGAL_ENVELOPE_3_USE_EDGE_HANDLES + typedef Intern::Intersection_info< Surface_pair_3 > Info; + + //! type of base class + typedef + CGAL::Arr_curve_data_traits_2< Arrangement_traits_2, Info, Info > Base; +#else + //! type of base class + typedef Arrangement_traits_2 Base; +#endif + + //! type of Z_Stack + typedef typename Restricted_cad_3::Z_stack Z_stack; + + +public: + + //! class of point + typedef typename Base::Point_2 Point_2; + + //! class of x-monotone curve + typedef typename Base::X_monotone_curve_2 X_monotone_curve_2; + +public: + + /*!\brief + * Subdivide the given surface into envelope relevant xy-monotone + * parts, and insert them into the output iterator. + * + * The iterator value-type is Xy_monotone_surface_3 + */ + class Make_xy_monotone_3 + { + protected: + const Self *parent; + public: + Make_xy_monotone_3(const Self* p) + : parent(p) + {} + + + template + OutputIterator operator()( + const Surface_3& s, + bool is_lower, + OutputIterator oi) + { + //parent->total_timer.start(); + + // TODO surfaces must be coprime? + // TASK Ask RW, MM for triangles, and planes + + // we just apply the sophisticated stuff in sqff_3 + // that also deals with the case of vertical components + typename + Surface_z_at_xy_isolator_traits::Square_free_factorization_3 + sqff_3; // TODO unique instance! + + std::list< int > mults; + sqff_3(s.f(), oi, std::back_inserter(mults)); + + //parent->total_timer.stop(); + + return oi; + } + }; + + /*! Get a Make_xy_monotone_3 functor object. */ + Make_xy_monotone_3 make_xy_monotone_3_object() + { + return Make_xy_monotone_3(this); + } + + /*!\brief + * Insert all 2D curves, which form of the boundary of the + * vertical projection of s onto the xy-plane, into the output iterator. + * The iterator value-type is Curve_2. + */ + class Construct_projected_boundary_2 { + protected: + const Self *parent; + public: + Construct_projected_boundary_2(const Self* p) + : parent(p) + {} + + template + OutputIterator operator()( + const Xy_monotone_surface_3& s, + OutputIterator oi) { + +#if !NDEBUG + std::cout << "Construct_projected_boundary ... " << std::flush; +#endif + //parent->total_timer.start(); + //parent->pboundary_timer.start(); + + // create cached instance + Restricted_cad_3 cad = Restricted_cad_3::cad_cache()(s); + + // we run over all segments/isolated points of cad and + // check whether they belong to boundary of s + + for (typename Restricted_cad_3::Vertex_const_iterator vit = + cad.vertices_begin(); + vit != cad.vertices_end(); vit++) { + if (vit->is_isolated() && cad.has_silhouette(vit)) { + Z_stack z_stack = cad.z_stack(vit); + // we are only interested in the lowest boundary! + if (!z_stack.is_empty()) { + if (z_stack.level_of_surface_in_z_cell(s,0) == 0) { + *oi++ = CGAL::make_object(vit->point()); + } + } + } + } + for (typename Restricted_cad_3::Edge_const_iterator eit = + cad.edges_begin(); + eit != cad.edges_end(); eit++) { + if (cad.has_silhouette(eit)) { + Z_stack z_stack = cad.z_stack(eit); + if (!z_stack.is_empty()) { + // we are only interested in the lowest boundary + // TODO check whether it is cheaper to check z_stacks + // of incident faces first, instead of + // checking z_stacks for eit + if (z_stack.level_of_surface_in_z_cell(s,0) == 0) { + CGAL::Arr_halfedge_direction dir = + eit->direction(); + // check both incident faces of eit: + // 1) both stacks are empty: + // return ON_ORIENTED_BOUNDARY + // 2) both stacks are non-empty: + // this part of the curve isn't + // a true boundary of s + // 3) one is empty, the other not + // we indicate the non-empty side! + int k = cad.z_stack( + eit->face() + ).number_of_z_cells(); + int l = cad.z_stack( + eit->twin()->face() + ).number_of_z_cells(); + CGAL_assertion(k >= 0); + CGAL_assertion(l >= 0); + CGAL::Oriented_side side = + CGAL::ON_ORIENTED_BOUNDARY; + if (k > 0) { + if (l > 0) { + continue; + } else { + CGAL_assertion(l == 0); + side = CGAL::ON_NEGATIVE_SIDE; + } + } else { + CGAL_assertion(k == 0); + if (l > 0) { + side = CGAL::ON_POSITIVE_SIDE; + } + } + if (dir == CGAL::ARR_LEFT_TO_RIGHT) { + side = -side; + } + Edge_const_handle invalid_eh; + *oi++ = CGAL::make_object( + std::make_pair( + X_monotone_curve_2(eit->curve()), + side + ) + ); + } + } + } + } + + //parent->total_timer.stop(); + //parent->pboundary_timer.stop(); + +#if !NDEBUG + std::cout << "done." << std::endl; +#endif + return oi; + } + }; + + /*! Get a Construct_projected_boundary_2 functor object. */ + Construct_projected_boundary_2 + construct_projected_boundary_2_object() const + { + return Construct_projected_boundary_2(this); + } + + + /*!\brief + * Insert all the 2D projections (onto the xy-plane) of the + * intersection objects between s1 and s2 into the output iterator. + * + * The iterator value-type is Object. An Object may be: + * 1. A pair, where the intersection + * type is an enumeration that can take the values + * {Transversal, Tangency, Unknown}. + * 2. A Point_2 instance (in degenerate cases). + */ + class Construct_projected_intersections_2 { + protected: + const Self *parent; + public: + Construct_projected_intersections_2(const Self* p) + : parent(p) + {} + + template + OutputIterator operator()( + const Xy_monotone_surface_3& s1, + const Xy_monotone_surface_3& s2, + OutputIterator oi) { +#if !NDEBUG + std::cout << "Construct_projected_intersections ... " + << std::flush; +#endif + //parent->total_timer.start(); + //parent->intersection_timer.start(); + + Surface_pair_3 pair = + Surface_pair_3::surface_pair_cache()(std::make_pair(s1, s2)); + + Restricted_cad_3 cad = pair.silhouettes_cut(); + + // run over all segments/isolated points of cad and + // check whether they belong to intersection of s1 and s2 + for (typename Restricted_cad_3::Vertex_const_iterator vit = + cad.vertices_begin(); + vit != cad.vertices_end(); vit++) { + if (vit->is_isolated() && cad.has_cut(vit)) { + Z_stack z_stack = cad.z_stack(vit); + if (!z_stack.is_empty()) { + if (z_stack.level_of_surface_in_z_cell(s1,0) == 0 && + z_stack.level_of_surface_in_z_cell(s2,0) == 0) { + *oi++ = CGAL::make_object(vit->point()); + } + } + } + } + + for (typename Restricted_cad_3::Edge_const_iterator eit = + cad.edges_begin(); + eit != cad.edges_end(); eit++) { + if (cad.has_cut(eit)) { + Z_stack z_stack = cad.z_stack(eit); + if (!z_stack.is_empty()) { + if (z_stack.level_of_surface_in_z_cell(s1,0) == 0 && + z_stack.level_of_surface_in_z_cell(s2,0) == 0) { + // TODO compute multiplicity of intersection + // if mult = 1 -> give 1 + // if mult = 2 + // -> check if unique intersection - give 2 + // if mult > 2 + // -> complex intersections destroy picture + Multiplicity multiplicity = 0; +#if CGAL_ENVELOPE_3_USE_EDGE_HANDLES + Info info(pair, eit); +#endif + *oi++ = CGAL::make_object( + std::make_pair( + // store eit in info +#if CGAL_ENVELOPE_3_USE_EDGE_HANDLES + X_monotone_curve_2(eit->curve(), + info), +#else + X_monotone_curve_2(eit->curve()), +#endif + multiplicity + ) + ); + } + } + } + } + + //parent->total_timer.stop(); + //parent->intersection_timer.stop(); + +#if !NDEBUG + std::cout << "done." << std::endl; +#endif + + return oi; + } + }; + + /*! Get a Construct_projected_intersections_2 functor object. */ + Construct_projected_intersections_2 + construct_projected_intersections_2_object() const + { + return Construct_projected_intersections_2(this); + } + +private: + + //! compares the vertical alignment of \c s1 and \c s2 along \c z_stack + static + CGAL::Comparison_result _intern_compare(const Xy_monotone_surface_3& s1, + const Xy_monotone_surface_3& s2, + const Z_stack& z_stack, + bool reverse) { + CGAL_assertion(!z_stack.is_empty()); + int i1 = z_stack.level_of_surface_in_z_cell(s1,0); + int i2 = z_stack.level_of_surface_in_z_cell(s2,0); + + CGAL_assertion(i1 >= -1 || i2 >= -1); + CGAL_assertion(i1 < 1); + CGAL_assertion(i2 < 1); + if (i1 == i2) { + CGAL_assertion(i1 == 0); + CGAL_assertion(i2 == 0); + return CGAL::EQUAL; + } else if (i1 == 0) { + if (reverse) { + return CGAL::LARGER; + } + return CGAL::SMALLER; + } + // else + CGAL_assertion(i2 == 0); + if (reverse) { + return CGAL::SMALLER; + } + return CGAL::LARGER; + } + +public: + + /*!\brief + * Check if the surface s1 is closer/equally distanced/farther + * from the envelope with respect to s2 at the xy-coordinates of p/c. + */ + class Compare_z_at_xy_3 { + protected: + const Self *parent; + public: + Compare_z_at_xy_3(const Self* p) + : parent(p) + {} + + //! compare over point + CGAL::Comparison_result operator()( + const Point_2& p, + const Xy_monotone_surface_3& s1, + const Xy_monotone_surface_3& s2) const { +#if !NDEBUG + std::cout << "Compare over point ... " << std::flush; +#endif + //parent->total_timer.start(); + //parent->compare_timer.start(); + + Surface_pair_3 pair = + Surface_pair_3::surface_pair_cache()(std::make_pair(s1, s2)); + + bool reverse = (s1 == pair.surface2()); + + Restricted_cad_3 cad = pair.silhouettes_cut(); + + // locate point + // TODO develop strategy to avoid point location! + Z_stack z_stack = cad.z_stack_for(p); + + // and compare lowest z-cells of s1 and s2 + CGAL::Comparison_result result = + _intern_compare(s1, s2, z_stack, reverse); + + //parent->total_timer.stop(); + //parent->compare_timer.stop(); +#if !NDEBUG + std::cout << "done." << std::endl; +#endif + + return result; + } + + //! compare over curve + CGAL::Comparison_result operator() ( + const X_monotone_curve_2& cv, + const Xy_monotone_surface_3& s1, + const Xy_monotone_surface_3& s2) const { + +#if !NDEBUG + std::cout << "Compare over curve ... " << std::flush; +#endif + //parent->total_timer.start(); + //parent->compare_on_cv_timer.start(); + + Surface_pair_3 pair = + Surface_pair_3::surface_pair_cache()(std::make_pair(s1, s2)); + + bool reverse = (s1 == pair.surface2()); + + Restricted_cad_3 cad = pair.silhouettes_cut(); + + typedef Restricted_cad_3_accessor< Restricted_cad_3 > Accessor; + Accessor acc(cad); + + // construct point on c + Point_2 p = acc.point_in_interior(cv); + + // locate point + // TODO develop strategy to avoid construction/point location! + Z_stack z_stack = cad.z_stack_for(p); + + // and compare lowest z-cells of s1 and s2 + CGAL::Comparison_result result = + _intern_compare(s1, s2, z_stack, reverse); + + // the two surfaces are not allowed to intersect ove c + CGAL_postcondition(result != CGAL::EQUAL); + + //parent->total_timer.stop(); + //parent->compare_on_cv_timer.stop(); + +#if !NDEBUG + std::cout << "done." << std::endl; +#endif + + return result; + } + + //! compare over full surfaces + CGAL::Comparison_result operator() ( + const Xy_monotone_surface_3& s1, + const Xy_monotone_surface_3& s2) const { +#if !NDEBUG + std::cout << "Compare over face ... " << std::flush; +#endif + //parent->total_timer.start(); + //parent->compare_on_face_timer.stare(); + + Surface_pair_3 pair = + Surface_pair_3::surface_pair_cache()(std::make_pair(s1, s2)); + + bool reverse = (s1 == pair.surface2()); + + Restricted_cad_3 cad = pair.silhouettes_cut(); + + // use z_stack of only existing face + Z_stack z_stack = cad.z_stack(cad.faces_begin()); + + // and compare lowest z-cells of s1 and s2 + CGAL::Comparison_result result = + _intern_compare(s1, s2, z_stack, reverse); + + // the two surfaces are not allowed to intersect at all + CGAL_postcondition(result != CGAL::EQUAL); + + //parent->total_timer.stop(); + //parent->compare_on_face_timer.stop(); + +#if !NDEBUG + std::cout << "done." << std::endl; +#endif + return result; + } + }; + + /*! Get a Compare_z_at_xy_3 functor object. */ + Compare_z_at_xy_3 + compare_z_at_xy_3_object() const + { + return Compare_z_at_xy_3(this); + } + + /*!\brief + * Check if the surface s1 is closer/equally distanced/farther + * from the envelope with + * respect to s2 immediately below the curve c. + */ + class Compare_z_at_xy_below_3 { + protected: + const Self *parent; + public: + Compare_z_at_xy_below_3(const Self* p) + : parent(p) + {} + CGAL::Comparison_result operator() ( + const X_monotone_curve_2& c, + const Xy_monotone_surface_3& s1, + const Xy_monotone_surface_3& s2) const { +#if !NDEBUG + std::cout << "Compare below curve ... " << std::flush; +#endif + //parent->total_timer.start(); + //parent->compare_side_timer.start(); + + Surface_pair_3 pair = + Surface_pair_3::surface_pair_cache()(std::make_pair(s1, s2)); + + bool reverse = (s1 == pair.surface2()); + + Restricted_cad_3 cad = pair.silhouettes_cut(); + + // not invalid, i.e., formed by an intersection and not a boundary +#if CGAL_ENVELOPE_3_USE_EDGE_HANDLES + typename Info::const_iterator dit = c.data().find(pair); + CGAL_assertion(dit != c.data().end()); + Edge_const_handle eh = dit->second; + + CGAL_assertion(cad.has_cut(eh,s1,s2)); + + // locate halfedge of c in cad using stored information + // as we want to compare "below" we chose the face incident the + // the halfegde of the pair that is directed from RIGHT_TO_LEFT + // and compare the lowest z-cells over this face + // Remark: This also works when c is vertical! + typename Restricted_cad_3::Face_const_handle fh = + (eh->direction() == CGAL::ARR_RIGHT_TO_LEFT ? + eh->face() : eh->twin()->face()); + Z_stack z_stack = cad.z_stack(fh); +#else + typedef Restricted_cad_3_accessor< Restricted_cad_3 > Accessor; + Accessor acc(cad); + + // construct point on c + Point_2 p = acc.point_in_interior(c); + + // TODO develop strategy to avoid point location! + CGAL::Object obj = cad.locate(p); + typename Accessor::Halfedge_const_handle heh; + CGAL_assertion_code(bool check = ) + CGAL::assign(heh, obj); + CGAL_assertion(check); + + // locate halfedge of c in cad using stored information + // as we want to compare "below" we chose the face incident the + // the halfegde of the pair that is directed from RIGHT_TO_LEFT + // and compare the lowest z-cells over this face + // Remark: This also works when c is vertical! + typename Restricted_cad_3::Face_const_handle fh = + (heh->direction() == CGAL::ARR_RIGHT_TO_LEFT ? + heh->face() : heh->twin()->face()); + Z_stack z_stack = cad.z_stack(fh); +#endif + + // and compare lowest z-cells of s1 and s2 + CGAL::Comparison_result result = + _intern_compare(s1, s2, z_stack, reverse); + + // the two surfaces are not allowed to intersect in a + // two-dimensional patch + CGAL_postcondition(result != CGAL::EQUAL); + + //parent->total_timer.stop(); + //parent->compare_side_timer.stop(); + +#if !NDEBUG + std::cout << "done." << std::endl; +#endif + return result; + } + + }; + + /*! Get a Compare_z_at_xy_below_3 functor object. */ + Compare_z_at_xy_below_3 + compare_z_at_xy_below_3_object() const + { + return Compare_z_at_xy_below_3(this); + } + + /*!\brief + * Check if the surface s1 is closer/equally distanced/farther + * from the envelope with + * respect to s2 immediately above the curve c. + */ + class Compare_z_at_xy_above_3 { + protected: + const Self *parent; + public: + Compare_z_at_xy_above_3(const Self* p) + : parent(p) + {} + + CGAL::Comparison_result operator()( + const X_monotone_curve_2& c, + const Xy_monotone_surface_3& s1, + const Xy_monotone_surface_3& s2) { +#if !NDEBUG + std::cout << "Compare above curve ... " << std::flush; +#endif + //parent->total_timer.start(); + //parent->compare_side_timer.start(); + + Surface_pair_3 pair = + Surface_pair_3::surface_pair_cache()(std::make_pair(s1, s2)); + + bool reverse = (s1 == pair.surface2()); + + Restricted_cad_3 cad = pair.silhouettes_cut(); + +#if CGAL_ENVELOPE_3_USE_EDGE_HANDLES + // not invalid, i.e., formed by an intersection and not a boundary + typename Info::const_iterator dit = c.data().find(pair); + CGAL_assertion(dit != c.data().end()); + Edge_const_handle eh = dit->second; + + CGAL_assertion(cad.has_cut(eh,s1,s2)); + + // locate halfedge of c in cad using stored information + // as we want to compare "above" we chose the face incident the + // the halfegde of the pair that is directed from LEFT_TO_RIGHT + // and compare the lowest z-cells over this face + // Remark: This also works when c is vertical! + typename Restricted_cad_3::Face_const_handle fh = + (eh->direction() == CGAL::ARR_LEFT_TO_RIGHT ? + eh->face() : eh->twin()->face()); + Z_stack z_stack = cad.z_stack(fh); +#else + typedef Restricted_cad_3_accessor< Restricted_cad_3 > Accessor; + Accessor acc(cad); + + // construct point on c + Point_2 p = acc.point_in_interior(c); + + // TODO develop strategy to avoid point location! + CGAL::Object obj = cad.locate(p); + typename Accessor::Halfedge_const_handle heh; + CGAL_assertion_code(bool check = ) + CGAL::assign(heh, obj); + CGAL_assertion(check); + + // locate halfedge of c in cad using stored information + // as we want to compare "below" we chose the face incident the + // the halfegde of the pair that is directed from RIGHT_TO_LEFT + // and compare the lowest z-cells over this face + // Remark: This also works when c is vertical! + typename Restricted_cad_3::Face_const_handle fh = + (heh->direction() == CGAL::ARR_LEFT_TO_RIGHT ? + heh->face() : heh->twin()->face()); + Z_stack z_stack = cad.z_stack(fh); +#endif + // and compare lowest z-cells of s1 and s2 + CGAL::Comparison_result result = + _intern_compare(s1, s2, z_stack, reverse); + + // the two surfaces are not allowed to intersect in a + // two-dimensional patch + CGAL_postcondition(result != CGAL::EQUAL); + + //parent->total_timer.stop(); + //parent->compare_side_timer.stop(); + +#if !NDEBUG + std::cout << "done." << std::endl; +#endif + return result; + } + }; + + /*! Get a Compare_z_at_xy_above_3 functor object. */ + Compare_z_at_xy_above_3 + compare_z_at_xy_above_3_object() const + { + return Compare_z_at_xy_above_3(this); + } + +}; // Surface_3_envelope_traits + +} // namespace SoX + +#endif // SoX_GAPS_SURFACE_3_ENVELOPE_TRAITS +// EOF diff --git a/Curved_kernel_via_analysis_2/include/CGAL/Arrangement_2l/Surface_pair_3.h b/Curved_kernel_via_analysis_2/include/CGAL/Arrangement_2l/Surface_pair_3.h new file mode 100644 index 00000000000..ac8bc10b829 --- /dev/null +++ b/Curved_kernel_via_analysis_2/include/CGAL/Arrangement_2l/Surface_pair_3.h @@ -0,0 +1,324 @@ +// ============================================================================ +// +// Copyright (c) 2001-2008 Max-Planck-Institut Saarbruecken (Germany). +// All rights reserved. +// +// This file is part of EXACUS (http://www.mpi-inf.mpg.de/projects/EXACUS/); +// you may redistribute it under the terms of the Q Public License version 1.0. +// See the file LICENSE.QPL distributed with EXACUS. +// +// Licensees holding a valid commercial license may use this file in +// accordance with the commercial license agreement provided with the software. +// +// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE +// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. +// +// ---------------------------------------------------------------------------- +// +// Library : SoX +// File : include/SoX/GAPS/Surface_pair_3.h +// SoX_release : $Name: $ +// Revision : $Revision$ +// Revision_date : $Date$ +// +// Author(s) : Eric Berberich +// +// ============================================================================ + +#ifndef SoX_GAPS_SURFACE_PAIR_3_H +#define SoX_GAPS_SURFACE_PAIR_3_H 1 + +/*!\file SoX/GAPS/Surface_pair_3.h + * \brief + * definition of \c Surface_pair_3<> + */ + +#include + +#include + +#include + +#include +#include +#include +#include + +namespace SoX { + +namespace Intern { + +template < class SurfaceZAtXyIsolatorTraits > +class Surface_pair_3_rep { + +public: + //! this instance's template parameter + typedef SurfaceZAtXyIsolatorTraits Surface_z_at_xy_isolator_traits; + + SoX_SURFACE_Z_AT_XY_ISOLATOR_TRAITS_SNAP_TYPEDEFS( + Surface_z_at_xy_isolator_traits + ); + + //! this instance itself + typedef Surface_pair_3_rep< Surface_z_at_xy_isolator_traits > Self; + + //! type of restricted cad + typedef Restricted_cad_3< Surface_z_at_xy_isolator_traits > + Restricted_cad_3; + + //! type of creator + typedef Create_restricted_cad_3< Surface_z_at_xy_isolator_traits > Creator; + + //! type of creator + typedef Overlay_restricted_cad_3< Surface_z_at_xy_isolator_traits > + Overlayer; + +public: + //!\name Constructors + //!@{ + + Surface_pair_3_rep(const Surface_3& surface1, const Surface_3& surface2) : + _m_surface1(surface1), + _m_surface2(surface2) { + } + + //!@} + +// TODO make again private: + + //! surface1 + mutable Surface_3 _m_surface1; + + //! surface2 + mutable Surface_3 _m_surface2; + + //! restricted cad of surface1 + mutable boost::optional< Restricted_cad_3 > _m_silhouette1; + + //! restricted cad of surface2 + mutable boost::optional< Restricted_cad_3 > _m_silhouette2; + + //! restricted cad of surface2 + mutable boost::optional< Restricted_cad_3 > _m_silhouettes; + + //! restricted cad of intersection + mutable boost::optional< Restricted_cad_3 > _m_cut; + + //! restricted cad of first silhouette with x + mutable boost::optional< Restricted_cad_3 > _m_silhouette1cut; + + //! restricted cad of second silhouette with x + mutable boost::optional< Restricted_cad_3 > _m_silhouette2cut; + + //! restricted cad + mutable boost::optional< Restricted_cad_3 > _m_silhouettescut; + + //! friends + friend class Surface_pair_3< Surface_z_at_xy_isolator_traits, Self > ; +}; + +} // namespace Intern + +template < +class SurfaceZAtXyIsolatorTraits, +class Rep_ = Intern::Surface_pair_3_rep < SurfaceZAtXyIsolatorTraits > +> +class Surface_pair_3 : + public ::CGAL::Handle_with_policy< Rep_ > { + +public: + //! this instance's first template parameter + typedef SurfaceZAtXyIsolatorTraits Surface_z_at_xy_isolator_traits; + + //! this instance's second template parameter + typedef Rep_ Rep; + + SoX_SURFACE_Z_AT_XY_ISOLATOR_TRAITS_SNAP_TYPEDEFS( + Surface_z_at_xy_isolator_traits + ); + + //! type of Base + typedef ::CGAL::Handle_with_policy< Rep > Base; + + //! type of restricted cad + typedef typename Rep::Restricted_cad_3 Restricted_cad_3; + + //! type of creator + typedef typename Rep::Creator Creator; + + //! type of overlayer + typedef typename Rep::Overlayer Overlayer; + + //! type of Z_stack + typedef typename Restricted_cad_3::Z_stack Z_stack; + + //!\name Constructors + //!@{ + + Surface_pair_3(const Surface_3& surface1, const Surface_3& surface2) : + Base(Rep(surface1, surface2)) { + } + + //!@} + + //!\name Access members + //!@{ + + //! the first surface + Surface_3 surface1() const { + return this->ptr()->_m_surface1; + } + + //! the second surface + Surface_3 surface2() const { + return this->ptr()->_m_surface2; + } + + //!@} + + //!\name Cads + //!@{ + + //! rscad of first silhouette + Restricted_cad_3 silhouette1() const { + if (!this->ptr()->_m_silhouette1) { + Creator creator; // TODO static?? + // create rs_cad_3 for surface1 + this->ptr()->_m_silhouette1 = creator(this->ptr()->_m_surface1); + } + return *this->ptr()->_m_silhouette1; + } + + //! rscad of second silhouette + Restricted_cad_3 silhouette2() const { + if (!this->ptr()->_m_silhouette2) { + Creator creator; // TODO static?? + // create rs_cad_3 for surface2 + this->ptr()->_m_silhouette2 = creator(this->ptr()->_m_surface2); + } + return *this->ptr()->_m_silhouette2; + } + +protected: + //! rscad of both silhouettes + Restricted_cad_3 _silhouettes() const { + if (!this->ptr()->_m_silhouettes) { + Overlayer overlay; // TODO static? + this->ptr()->_m_silhouettes = + overlay(silhouette1(), silhouette2()); + } + return *this->ptr()->_m_silhouettes; + } + +public: + //! rscad of cut + Restricted_cad_3 cut() const { + if (!this->ptr()->_m_cut) { + // create rs_cad_3 for intersection of surface1 and surface2 + Creator creator; // TODO static?? + this->ptr()->_m_cut = + creator(this->ptr()->_m_surface1, this->ptr()->_m_surface2); + } + return *this->ptr()->_m_cut; + } + +protected: + //! rscad of first silhouette with cut + Restricted_cad_3 _silhouette1_cut() const { + if (!this->ptr()->_m_silhouette1cut) { + Overlayer overlay; // TODO static? + this->ptr()->_m_silhouette1cut = overlay(silhouette1(), cut()); + } + return *this->ptr()->_m_silhouette1cut; + } + + //! rscad of first silhouette with cut + Restricted_cad_3 _silhouette2_cut() const { + if (!this->ptr()->_m_silhouette2cut) { + Overlayer overlay; // TODO static? + this->ptr()->_m_silhouette2cut = overlay(silhouette2(), cut()); + } + return *this->ptr()->_m_silhouette2cut; + } + +public: + //! rscad of first silhouette with cut + Restricted_cad_3 silhouettes_cut() const { + if (!this->ptr()->_m_silhouettescut) { + Overlayer overlay; // TODO static? + this->ptr()->_m_silhouettescut = overlay(_silhouettes(), cut()); + } + return *this->ptr()->_m_silhouettescut; + } + + //!@} + + //!\name Curves + //!@{ + + /*!\brief + * returns boundary curves and points for a given \c surface in pair + */ + template < class CurveOutputIterator, class PointOutputIterator > + void silhouette_objects( + const Surface_3& surface, + CurveOutputIterator coi, PointOutputIterator poi) const { + + const Restricted_cad_3& rsc = (surface == this->surface1() ? + this->silhouette1() : + this->silhouette2()); + + for (typename Restricted_cad_3::Vertex_const_iterator vit = + rsc.vertices_begin(); + vit != rsc.vertices_end(); vit++) { + if (vit->is_isolated()) { + *poi++ = vit->point(); + } + } + for (typename Restricted_cad_3::Edge_const_iterator eit = + rsc.edges_begin(); + eit != rsc.edges_end(); eit++) { + *coi++ = eit->curve(); + } + } + + /*!\brief + * returns projected intersection curves and points of the pair + */ + template < class CurveOutputIterator, class PointOutputIterator > + void cut_objects( + CurveOutputIterator coi, PointOutputIterator poi) const { + + const Restricted_cad_3& rsc = this->cut(); + + for (typename Restricted_cad_3::Vertex_const_iterator vit = + rsc.vertices_begin(); + vit != rsc.vertices_end(); vit++) { + if (vit->is_isolated()) { + *poi++ = vit->point(); + } + } + for (typename Restricted_cad_3::Edge_const_iterator eit = + rsc.edges_begin(); + eit != rsc.edges_end(); eit++) { + *coi++ = eit->curve(); + } + } + + //!@} + + //!\name Z_stack for compare_xyz + //!@{ + + //! returns z-stack of silhouttes_cut() for \c point + Z_stack z_stack_for(const Point_2& point) const { + return silhouettes_cut().z_stack_for(point); + } + + //!@} +}; + +} // namespace SoX + +#endif // SoX_GAPS_SURFACE_PAIR_3_H +// EOF diff --git a/Curved_kernel_via_analysis_2/include/CGAL/Arrangement_2l/Z_stack.h b/Curved_kernel_via_analysis_2/include/CGAL/Arrangement_2l/Z_stack.h new file mode 100644 index 00000000000..8e35561aa5c --- /dev/null +++ b/Curved_kernel_via_analysis_2/include/CGAL/Arrangement_2l/Z_stack.h @@ -0,0 +1,1629 @@ +// ============================================================================ +// +// Copyright (c) 2001-2008 Max-Planck-Institut Saarbruecken (Germany). +// All rights reserved. +// +// This file is part of EXACUS (http://www.mpi-inf.mpg.de/projects/EXACUS/); +// you may redistribute it under the terms of the Q Public License version 1.0. +// See the file LICENSE.QPL distributed with EXACUS. +// +// Licensees holding a valid commercial license may use this file in +// accordance with the commercial license agreement provided with the software. +// +// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE +// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. +// +// ---------------------------------------------------------------------------- +// +// Library : SoX +// File : include/SoX/GAPS/Z_stack.h +// SoX_release : $Name: $ +// Revision : $Revision$ +// Revision_date : $Date$ +// +// Author(s) : Eric Berberich +// +// ============================================================================ + +#ifndef SoX_GAPS_Z_STACK_H +#define SoX_GAPS_Z_STACK_H 1 + +/*! \file SoX/GAPS/Z_stack.h + * \brief definition of Z_stack class template + */ + +#include + +#include + +#include + +#include + +#include +#include +#include +#include +#include + +namespace SoX { + +namespace Intern { + +template < class SurfaceZAtXyIsolatorTraits, class DcelData > +class Z_cell_rep { +public: + + //! this instance's first template parameter + typedef SurfaceZAtXyIsolatorTraits Surface_z_at_xy_isolator_traits; + + //! this instance's second template parameter + typedef DcelData Dcel_data; + + //! the class itself + typedef Z_cell_rep< Surface_z_at_xy_isolator_traits, Dcel_data > Self; + + SoX_SURFACE_Z_AT_XY_ISOLATOR_TRAITS_SNAP_TYPEDEFS( + Surface_z_at_xy_isolator_traits + ); + + //! type of z_stack + typedef Z_stack< Surface_z_at_xy_isolator_traits, Dcel_data > Z_stack; + + //! type of Surface less + typedef typename Surface_3::Surface_less_than Surface_less_than; + + //! type of sheet map + typedef std::map< Surface_3, int, Surface_less_than > Sheet_map; + + //! type of Boundary + typedef typename Z_at_xy_isolator::Boundary Boundary; + + Z_cell_rep(const Z_stack& z_stack, int z) : + _m_z_stack(z_stack), + _m_z(z), + _m_min_length(0) { + } + +private: + //! the z_stack + Z_stack _m_z_stack; + + //! z-position in stacke + int _m_z; + + //! sheet of surface + Sheet_map _m_sheets; + + //! current min_length + Boundary _m_min_length; + + //! surface that induces min_length + Surface_3 _m_min_length_surface; + + //! surface sheet that induces min_length + int _m_min_length_surface_sheet; + +#if 0 + //! + std::map< CGAL::Object, int, Object_less_than > _m_adjacent_items; +#endif + + friend class Z_cell< Surface_z_at_xy_isolator_traits, Dcel_data >; + +}; + +template < class SurfaceZAtXyIsolatorTraits, class DcelData > +class Z_cell : public +CGAL::Handle_with_policy< Z_cell_rep > +{ +public: + + //! this instance's first template parameter + typedef SurfaceZAtXyIsolatorTraits Surface_z_at_xy_isolator_traits; + + //! this instance's second template parameter + typedef DcelData Dcel_data; + + //! type of rep + typedef Z_cell_rep< Surface_z_at_xy_isolator_traits, Dcel_data > Rep; + + //! type of base + typedef CGAL::Handle_with_policy< Rep > Base; + + //! the class itself + typedef Z_cell< Surface_z_at_xy_isolator_traits, Dcel_data > Self; + + SoX_SURFACE_Z_AT_XY_ISOLATOR_TRAITS_SNAP_TYPEDEFS( + Surface_z_at_xy_isolator_traits + ); + + //! type of sheet map + typedef typename Rep::Sheet_map Sheet_map; + + + //! type of sheet map iterator + typedef typename Sheet_map::const_iterator Sheet_const_iterator; + + //! type of Boundary + typedef typename Rep::Boundary Boundary; + + //! type of z_stack + typedef Z_stack< Surface_z_at_xy_isolator_traits, Dcel_data > Z_stack; + + //!\name Constructor + //!@{ + + //! standard constructor for given \c z-stack and \c int + Z_cell(const Z_stack& z_stack, int z) : + Base(Rep(z_stack, z)) { + } + + //!@} + + //!\name Surfaces, sheets, and levels + //!@{ + + //! returns z-entry of cell + int z() { + return this->ptr()->_m_z; + } + + //! returns number of involved sheets + int number_of_sheets() const { + return static_cast< int >(this->ptr()->_m_sheets.size()); + } + + //! returns sheet number of surface in given cell + int sheet_number(const Surface_3& surface) const { + typename Sheet_map::const_iterator it = + this->ptr()->_m_sheets.find(surface); + if (it != this->ptr()->_m_sheets.end()) { + CGAL_postcondition(it->second >= 0); + return it->second; + } + return -1; + } + + //! returns beginning of sheets + Sheet_const_iterator sheets_begin() const { + return this->ptr()->_m_sheets.begin(); + } + + //! returns past-the-end iterator for sheets + Sheet_const_iterator sheets_end() const { + return this->ptr()->_m_sheets.end(); + } + + // TODO random surface + + //! return surface whose z-interval is minimal in cell + const Surface_3& minimal_length_surface() const { + CGAL_assertion(!this->ptr()->_m_sheets.empty()); + return this->ptr()->_m_min_length_surface; + } + + //! return sheet of surface whose z-interval is minimal in cell + int minimal_length_surface_sheet() const { + CGAL_assertion(!this->ptr()->_m_sheets.empty()); + return this->ptr()->_m_min_length_surface_sheet; + } + + //!@} + + //! Approximation + //!@{ + + //! returns double approximation of z-cell + std::pair< double, double > to_double() const { + CGAL_assertion(!this->ptr()->_m_sheets.empty()); + Surface_3 surface = this->ptr()->_m_min_length_surface; + int sheet = this->ptr()->_m_min_length_surface_sheet; + return std::make_pair( + CGAL::to_double( + this->ptr()->_m_z_stack._isolator(surface). + left_boundary(sheet) + ), + CGAL::to_double( + this->ptr()->_m_z_stack._isolator(surface). + right_boundary(sheet) + ) + ); + } + + //!@} + +private: + + //!\name Modifying + //!@{ + + //! add sheet \c sheet of \c surface to cell + void _add_sheet(const Surface_3& surface, int sheet, + Boundary interval_length) { + CGAL_precondition(sheet >= 0); + this->ptr()->_m_sheets.insert(std::make_pair(surface, sheet)); + if (interval_length < this->ptr()->_m_min_length || + this->ptr()->_m_min_length == Boundary(0)) { + this->ptr()->_m_min_length = interval_length; + this->ptr()->_m_min_length_surface = surface; + this->ptr()->_m_min_length_surface_sheet = sheet; + } + + } + + + //! makes cell unique for \c z's cell of \c z_stack + void _make_unique_for(const Z_stack& z_stack, int z) { + if (z_stack.id() != this->ptr()->_m_z_stack.id()) { + this->copy_on_write(); + this->ptr()->_m_z_stack = z_stack; + } + this->ptr()->_m_z = z; + } + + //!@} + +public: + //!\name IO + //!@{ + + /*!\brief + * prints pretty-formated version of z_stack + */ + void pretty_print(std::ostream& os) const { + os << "(z=" << this->ptr()->_m_z << ","; + std::pair< double, double > dp = to_double(); + os << "[" << dp.first << "," << dp.second << "],"; + for (typename Sheet_map::const_iterator iit = + this->ptr()->_m_sheets.begin(); + iit != this->ptr()->_m_sheets.end(); iit++) { + os << "<" << iit->first.id() << "," << iit->second << ">"; + typename Sheet_map::const_iterator tmp = iit; + tmp++; + if (tmp != this->ptr()->_m_sheets.end()) { + os << ","; + } + } + os << ")"; + } + + // friends + friend class SoX::Z_stack< Surface_z_at_xy_isolator_traits, Dcel_data >; + +}; // Z_cell + + + +template < class SurfaceZAtXyIsolatorTraits, class DcelData > +class Z_stack_rep { + +public: + //! this instance's first template parameter + typedef SurfaceZAtXyIsolatorTraits Surface_z_at_xy_isolator_traits; + + //! this instance's second template parameter + typedef DcelData Dcel_data; + + SoX_SURFACE_Z_AT_XY_ISOLATOR_TRAITS_SNAP_TYPEDEFS( + Surface_z_at_xy_isolator_traits + ); + + //! the class itself + typedef Z_stack_rep< Surface_z_at_xy_isolator_traits, Dcel_data > + Self; + + //! type of Surface less + typedef typename Surface_3::Surface_less_than Surface_less_than; + + //! type of z-isolator + typedef std::map< Surface_3, Z_at_xy_isolator, Surface_less_than > + Z_at_xy_isolators; + + //! type of z-cell + typedef Z_cell< Surface_z_at_xy_isolator_traits, Dcel_data > Z_cell; + + //! type of z-cell container + typedef std::list< Z_cell > Z_cell_container; + + //! type of const_iterator + typedef typename Z_cell_container::const_iterator Z_cell_const_iterator; + + //!\name Constructors + //!@{ + + + //! default constructor + Z_stack_rep() {} + + /*!\brief + * standard constructor with the help of a point \c xy + */ + Z_stack_rep(SoX::Dcel_feature feature, + const Dcel_data* data, const Point_2& pt) : + _m_dcel_feature(feature), + _m_dcel_data(data), + _m_point(pt) { + _m_z_at_xy_isolators.clear(); + _m_z_at_xy_empty_isolators.clear(); + _m_z_cells.clear(); + } + + //!@} + +private: + // data members + //! type of dcel_feature + mutable SoX::Dcel_feature _m_dcel_feature; + + //! storage for isolators + mutable Z_at_xy_isolators _m_z_at_xy_isolators; + + //! storage for empty isolators + mutable Z_at_xy_isolators _m_z_at_xy_empty_isolators; + + //! z-cells + // taken a list since we need to update the cells + // actually we also require RAM with the help of operator[] + mutable Z_cell_container _m_z_cells; + + //! dcel data + const Dcel_data* _m_dcel_data; + + //! stored point + mutable Point_2 _m_point; + + // friends + friend class Z_stack< Surface_z_at_xy_isolator_traits, Dcel_data >; +}; + +} // namespace Intern + +/*!\brief + * Class to describe the order of surfaces of a planar (refinable) point. + */ +template < class SurfaceZAtXyIsolatorTraits, class DcelData > +class Z_stack : public +::CGAL::Handle_with_policy< +Intern::Z_stack_rep< SurfaceZAtXyIsolatorTraits, DcelData > > { +public: + //! this instance's first template parameter + typedef SurfaceZAtXyIsolatorTraits Surface_z_at_xy_isolator_traits; + + //! this instance's second template parameter + typedef DcelData Dcel_data; + + SoX_SURFACE_Z_AT_XY_ISOLATOR_TRAITS_SNAP_TYPEDEFS( + Surface_z_at_xy_isolator_traits + ); + + //! type of rep + typedef Intern::Z_stack_rep< Surface_z_at_xy_isolator_traits, Dcel_data > + Rep; + + //! base type + typedef ::CGAL::Handle_with_policy< Rep > Base; + + //! the class itself + typedef Z_stack< Surface_z_at_xy_isolator_traits, Dcel_data > Self; + + //! type of isolator map + typedef typename Rep::Z_at_xy_isolators Z_at_xy_isolators; + + //! type of Z_stack_cell + typedef typename Rep::Z_cell Z_cell; + + //! type of Z_cell_container + typedef typename Rep::Z_cell_container Z_cell_container; + + //! type of const_iterator + typedef typename Rep::Z_cell_const_iterator Z_cell_const_iterator; + + //!\name Constructors + //!@{ + + //! defualt constructor + Z_stack() {}; + + /*!\brief + * standard constructor from a refinable point + */ + Z_stack(SoX::Dcel_feature feature, + const Dcel_data* data, const Point_2& pt) : + Base(Rep(feature,data,pt)) { + }; + + /*!\brief + * copies z-stack \c self, but replaces feature + */ + Z_stack(const Self& self, + SoX::Dcel_feature feature, const Dcel_data* data) : + Base(*self.ptr()) { + this->copy_on_write(); + this->ptr()->_m_dcel_feature = feature; + this->ptr()->_m_dcel_data = data; + }; + + //!@} + +private: + //!\name Query isolators + //!@{ + + /*!\brief + * returns \c true if z_stack already knows an isolator + * at stored point for \c surface + */ + bool _knows_isolator(const Surface_3& surface, bool& empty) const { + empty = true; + typename Z_at_xy_isolators::const_iterator it = + this->ptr()->_m_z_at_xy_isolators.find(surface); + if (it != this->ptr()->_m_z_at_xy_isolators.end()) { + return true; + } + // else + empty = false; + return (this->ptr()->_m_z_at_xy_empty_isolators.find(surface) != + this->ptr()->_m_z_at_xy_empty_isolators.end()); + } + + /*!\brief returns isolator for surface + * + * \pre _knows_isolator(surface) + */ + Z_at_xy_isolator _isolator(const Surface_3& surface) const { + CGAL_precondition_code(bool empty;); + CGAL_precondition(_knows_isolator(surface, empty)); + typename Z_at_xy_isolators::const_iterator it = + this->ptr()->_m_z_at_xy_isolators.find(surface); + if (it != this->ptr()->_m_z_at_xy_isolators.end()) { + return it->second; + } + // else + return this->ptr()->_m_z_at_xy_empty_isolators.find(surface)->second; + } + +public: + + /*!\brief + * returns \c true if z_stack already knows an isolator + * at stored point for \c surface + */ + inline + bool knows_isolator(const Surface_3& surface, bool& empty) const { + return _knows_isolator(surface, empty); + } + + /*!\brief returns isolator for surface + * + * \pre knows_isolator(surface) + */ + inline + Z_at_xy_isolator isolator(const Surface_3& surface) const { + return _isolator(surface); + } + + //!@} + +private: + //!\name Adding surfaces + //!@{ + + /*!\brief + * add surface to z_stack and maintain sequence of cells + */ + void _add_surface(const Surface_3& surface, + const Z_at_xy_isolator& isolator) const { + + // check whether surface influences z_stack + if (isolator.number_of_real_roots() == 0) { + typename Z_at_xy_isolators::iterator it = + this->ptr()->_m_z_at_xy_empty_isolators.find(surface); + if (it == this->ptr()->_m_z_at_xy_empty_isolators.end()) { + this->ptr()->_m_z_at_xy_empty_isolators.insert( + it, std::make_pair(surface,isolator) + ); + } + // surface has no cut with line parallel to z-axis + // -> nothing else to do + return; + } + // else + + // first check whether z_stack aldready knows something about the curve + typename Z_at_xy_isolators::iterator it = + this->ptr()->_m_z_at_xy_isolators.find(surface); + if (it == this->ptr()->_m_z_at_xy_isolators.end()) { + this->ptr()->_m_z_at_xy_isolators.insert( + it, std::make_pair(surface,isolator) + ); + } else { + // surface already exists in z_stack -> no further changes on stack + return; + } + + // maintenance + _update_cells(surface, isolator); + } + +private: + + /*!\brief + * updates the cells, i.e., find intersections! + */ + void _update_cells(const Surface_3& surface2, + const Z_at_xy_isolator& isolator2) const { + + int roots2 = isolator2.number_of_real_roots(); + + if (this->ptr()->_m_z_cells.empty()) { + // simple case, fast exit + for (int i = roots2 - 1; i >= 0; i--) { + Z_cell cell(*this, i); + cell._add_sheet(surface2, i, isolator2.length(i)); + this->ptr()->_m_z_cells.push_front(cell); + } + // exit function + return; + } + + // else + // for each 0 <= i < roots2: + // find position of (isol2.left_boundary(i), isol2.right_boundary(i)) + // in sequence of cells + + // we try to apply filters + // FACE: No intersection can take place + // EDGE/VERTEX: No projected intersection -> no equal_z-test + // EDGE: If mult == 1 => only one intersection + // VERTEX: not isolated -> adjacency + + // there are situation where filters do not work + // EDGE: If mult > 1 => equal-z + // VERTEX: isolated + + SoX::Dcel_feature feature = this->ptr()->_m_dcel_feature; + + typedef typename Surface_3::Surface_less_than Surface_less_than; + + typedef Restricted_cad_3< Surface_z_at_xy_isolator_traits > + Restricted_cad_3; + + Restricted_cad_3 cad = + Restricted_cad_3::cad_cache()( + this->ptr()->_m_dcel_data->_rs_id() + ); + + + SoX::Intern::Refineable_interval_helper< Z_at_xy_isolator > iv_helper; + + std::set< Surface_3, Surface_less_than > projected_sil1s; + std::map< Surface_3, std::pair< int, + std::list< std::pair < int, int > > + >, Surface_less_than > projected_cut1s; + + bool projected_sil2 = false; + bool projected_cut2 = false; + + typedef std::pair< Surface_3, int > Sheet; + + //! type of Less of sheet + typedef CGAL::Pair_lexicographical_less_than< Surface_3, int, + Surface_less_than, std::less< int > > Sheet_less; + + typedef std::map< Sheet, int, Sheet_less > Adj_x_map; + + Adj_x_map adj_x_map; + + if (feature != FACE) { + + // check existing projected intersections and boundaries + projected_sil2 = + this->ptr()->_m_dcel_data->_has_silhouette(surface2); + + for (typename Z_at_xy_isolators::iterator it = + this->ptr()->_m_z_at_xy_isolators.begin(); + it != this->ptr()->_m_z_at_xy_isolators.end(); it++) { + if (this->ptr()->_m_dcel_data->_has_silhouette(it->first)) { + projected_sil1s.insert(it->first); + } + if (this->ptr()->_m_dcel_data->_has_cut(it->first, surface2)) { + // compute overlaps + std::list< std::pair< int, int > > overlaps; + iv_helper.refined_overlaps(it->second, isolator2, + std::back_inserter(overlaps)); + + // and cache them + projected_cut1s.insert( + std::make_pair( + it->first, + std::make_pair( + // make sense for EDGE + // in the test *mult == 1 + (feature == SoX::VERTEX ? + -1 : + *this->ptr()-> + _m_dcel_data->multiplicity_of_cut( + it->first, surface2 + ) + ), + overlaps + ) + ) + ); + projected_cut2 = true; + } + } + + // for VERTEX: + if (feature == SoX::VERTEX && projected_cut2) { + // if there is a projected cut over a vertex + + + // Create adjacency information at z_stack-location from + // neighboring cells and for one non-isolated + // surfaces of each cell + + // search for neighboring cells that form an intersection + // of surface2 with at least one other surface + typename Restricted_cad_3::Vertex_const_handle vh; + + CGAL_assertion_code(bool check = ) + CGAL::assign(vh,this->ptr()->_m_dcel_data->_dcel_handle()); + CGAL_assertion(check); + + if (!vh->is_isolated()) { + + std::pair< Self, SoX::Dcel_feature > z_stack2to = + vh->data()->_z_stack_of_surface(surface2); + + boost::optional< + std::pair< SoX::Dcel_feature,CGAL::Object > + > opt2to = vh->data()->_dcel(surface2); + CGAL_assertion(opt2to); + CGAL_assertion( + opt2to->first == z_stack2to.second + ); + + typename Restricted_cad_3:: + Halfedge_around_vertex_const_circulator circ, start = + vh->incident_halfedges(); + circ = start; + + do { + ++circ; + Self z_stack = + circ->data()->_z_stack_for_halfedge_handle(circ); + + for (typename Z_cell_container::iterator it = + z_stack.z_cells_begin(); + it != z_stack.z_cells_end(); it++) { + // going over all cells + + // check whether surface2 and at least some + // other curve is involved + int sheet2 = it->sheet_number(surface2); + + if (it->number_of_sheets() > 1 && + sheet2 >= 0) { + + // compute adjacencies of surface2 towards goal + std::pair< Self, SoX::Dcel_feature > + z_stack2from = + circ->data()->_z_stack_of_surface( + surface2 + ); + boost::optional< + std::pair< SoX::Dcel_feature,CGAL::Object > + > opt2from = circ->data()->_dcel(surface2); + CGAL_assertion(opt2from); + CGAL_assertion( + opt2from->first == z_stack2from.second + ); + + std::pair< int, int > adj2 = + z_stack2from.first.adjacency( + surface2, sheet2, + opt2from->first, + opt2from->second, + z_stack2to.first, + opt2to->first, + opt2to->second + ); + + // TODO check adj [1,0] + CGAL_assertion(adj2.first >= 0 && + adj2.first <= adj2.second); + + // and adjacancies of other surfaces + // towards goal + for (typename Z_cell::Sheet_const_iterator + sit = it->sheets_begin(); + sit != it->sheets_end(); sit++) { + + Surface_3 surface1 = sit->first; + int sheet1 = sit->second; + + if (surface1.id() == surface2.id()) { + continue; + } + + std::pair< Self, SoX::Dcel_feature > + z_stack1from = + circ->data()->_z_stack_of_surface( + surface1 + ); + boost::optional< + std::pair< SoX::Dcel_feature, + CGAL::Object > + > opt1from = + circ->data()->_dcel(surface1); + CGAL_assertion(opt1from); + CGAL_assertion( + opt1from->first == + z_stack1from.second + ); + + std::pair< Self, SoX::Dcel_feature > + z_stack1to = + vh->data()->_z_stack_of_surface( + surface1 + ); + boost::optional< + std::pair< SoX::Dcel_feature, + CGAL::Object > + > opt1to = + vh->data()->_dcel(surface1); + CGAL_assertion(opt1to); + CGAL_assertion( + opt1to->first == + z_stack1to.second + ); + + std::pair< int, int > adj1 = + z_stack1from.first.adjacency( + surface1, sheet1, + opt1from->first, + opt1from->second, + z_stack1to.first, + opt1to->first, + opt1to->second); + + // TODO check adj [1,0] + CGAL_assertion(adj1.first >= 0 && + adj1.first <= adj1.second); + + // and store other surface-sheets as key +#if 0 + std::cout << "s1: " << surface1.id() + << "t1: " << adj1.first + << std::endl; + std::cout << "s2: " << surface2.id() + << "t2: " << adj2.first + << std::endl; +#endif + adj_x_map.insert( + std::make_pair( + std::make_pair( + surface1, + adj1.first), + adj2.first) + ); + + // whenever we later see the stored + // surface-sheet-key, we've actually + // detected an intersection + } + } + } + } while (circ != start); + } + } + } + // else: if z_stack is for FACE, we will not detect any intersection, + // as we require surface to be "coprime" + + // setup for filters finished + // actual merge starts now + int z = 0; + int root2 = 0; + + typename Z_cell_container::iterator it = + this->ptr()->_m_z_cells.begin(); + + bool still_possible_equal = true; + + while (root2 < roots2 && it != this->ptr()->_m_z_cells.end()) { + // read isolator and correspondent root + // TASK make selection of surface1 random! + const Surface_3& surface1 = it->minimal_length_surface(); + const int root1 = it->minimal_length_surface_sheet(); + + const Z_at_xy_isolator& isolator1 = + this->ptr()->_m_z_at_xy_isolators.find(surface1)->second; + + switch (iv_helper.overlap_or_order(isolator1, root1, + isolator2, root2)) { + case CGAL::SMALLER: { + //std::cout << "SMALLER root1: " << root1 << " root2: " + // << root2 << std::endl; + it->_make_unique_for(*this, z); + + // go to next existing cell in z_stack + it++; + z++; + + // the next combination can be equal again + still_possible_equal = true; + break; + } + case CGAL::EQUAL: { + //std::cout << "EQUAL root1: " << root1 << " root2: " << root2 + // << std::endl; + + bool surely_equal = false; + bool surely_unequal = false; + + bool projected_sil1 = + (projected_sil1s.find(surface1) != + projected_sil1s.end()); + + // still_possible_equal can already be *false*, i.e., + // if in last iterator CGAL::EQAUL was detected and + // refinement was triggered. Then CGAL:EQUAL is still + // a wrong overlap ..so directly refine! + + // General //////////////////////////////////////////////////// + + // filter: intersections can only happen for non-FACES + if (!surely_equal && !surely_unequal && still_possible_equal) { + still_possible_equal = (feature != SoX::FACE); + //std::cout << "filterA" << std::endl; + if (feature == SoX::FACE) { + //std::cout << "filterA2" << std::endl; + surely_unequal = true; + } + } + + // filter: is surface2 involved in any projected cut? + if (!surely_equal && !surely_unequal && still_possible_equal) { + still_possible_equal = projected_cut2; + //std::cout << "filterB" << std::endl; + if (!projected_cut2) { + //std::cout << "filterB2" << std::endl; + surely_equal = false; + surely_unequal = true; + } + } + + typename + std::map< Surface_3, + std::pair< int, + std::list< std::pair< int, int > > + >, Surface_less_than >::iterator sit; + + int mult = -1; + + // filter: is surface1 involved in any projected cut? + if (!surely_equal && !surely_unequal && still_possible_equal) { + //std::cout << "filterC" << std::endl; + sit = projected_cut1s.find(surface1); + + if (sit == projected_cut1s.end()) { + //std::cout << "filterC2" << std::endl; + still_possible_equal = false; + surely_equal = false; + surely_unequal = true; + } else { + //std::cout << "filterC3" << std::endl; + if (sit->second.second.empty()) { + // no overlapping found + still_possible_equal = false; + surely_equal = false; + surely_unequal = true; + } else { + mult = sit->second.first; + } + } + } + + // Edge /////////////////////////////////////////////////////// + + // filter: EDGE: at most one intersection? + if (!surely_equal && !surely_unequal && still_possible_equal) { + //std::cout << "filterD" << std::endl; + if (feature == SoX::EDGE) { + //std::cout << "filterD2" << "mult=" << mult << std::endl; + CGAL_assertion(mult > 0); + if (mult == 1) { + //std::cout << "filterD3" << std::endl; + std::pair< int, int > intersection = + iv_helper.unique_overlap( + isolator1, isolator2, + sit->second.second.begin(), + sit->second.second.end()); + if (intersection.first == root1 && + intersection.second == root2) { + //std::cout << "filterD4" << std::endl; + surely_equal = true; + surely_unequal = false; + } else { + //std::cout << "filterD5" << std::endl; + surely_equal = false; + surely_unequal = true; + } + } // else complex intersections can create + // the projection -> equal_z + } + } + + // the other case (mult > 1) must be handled by equal_z + + // Vertex ///////////////////////////////////////////////////// + // we now turn to vertices only + + // filter: VERTEX using adjacencies + if (!surely_equal && !surely_unequal && still_possible_equal) { + //std::cout << "filterE" << std::endl; + if (feature == SoX::VERTEX) { + //std::cout << "filterE2" << std::endl; + if (!adj_x_map.empty()) { + //std::cout << "filterE3" << std::endl; + // it is possible that we can derive the + // intersection by adjacencies, + // so let's give it a try + typename Adj_x_map::const_iterator + iit = adj_x_map.find( + std::make_pair(surface1, root1) + ); + if (iit != adj_x_map.end() && + iit->second == root2) { + //std::cout << "filterE4" << std::endl; + // adj-intersection detected + surely_equal = true; + surely_unequal = false; + } + } + } + } + + // Remark: When reching this point, two z-intervals are + // overlapping and the adjacency does not found + // an intersection. If no silhouette of + // surface1 and surface2 exists, there will be + // no further intersection + // The other way around: There can only be an + // intersection if an isolated object meets another one + // of "lies within a sheet that locally looks like a + // plane". + + // Filter: VERTEX further intersections require silhouette + if (!surely_equal && !surely_unequal && still_possible_equal) { + if (feature == SoX::VERTEX) { + if (!projected_sil2) { + + if (! projected_sil1) { + surely_equal = false; + surely_unequal = true; + } + } + } + } + + // Remark: We cannot claim that the existence of a single + // overlap induces a real intersection, + // this is only true oif the point is non-singular, + // i.e., mult = 1 (but we do not have this information) + // otherwise complex intersection can create the + // projected intersection. Therefore it is not possible + // to use the following filter + /* + // filter: only one overlap? + if (!surely_equal && !surely_unequal && still_possible_equal) { + // find possible pairings first: unique_overlap! + } + */ + + // So, here can assume that the point is singular + // and at least one silhouette exists + // Filter: + if (!surely_equal && !surely_unequal && still_possible_equal) { + if (feature == SoX::VERTEX) { + CGAL_assertion(projected_sil1 || projected_sil2); + } + } + + // If *at least one of) the current events is isolated, + // we can also do not deduce whether the overlapping + // intervals will seperate or not. So the following + // filter is also not possible -> need answer from equal_z + /* + // Filter: Surface-"sheet" at vertex isolated and + // NOT involved in overlaps -> no intersection + if (!surely_equal && !surely_unequal && still_possible_equal) { + + } + */ + + // TODO move filter before adjacenceny? + // Filter: If vertex is not a genuine vertex in cut-curve arr + if (!surely_equal && !surely_unequal && still_possible_equal) { + if (feature == SoX::VERTEX) { + boost::optional< + std::pair< SoX::Dcel_feature, CGAL::Object > > + dcel_pair = + this->ptr()->_m_dcel_data->_dcel( + surface1, surface2 + ); + CGAL_assertion(dcel_pair); + if (dcel_pair->first == SoX::EDGE) { + typename Restricted_cad_3::Halfedge_const_handle + hh; + CGAL_assertion_code(bool check = ) + CGAL::assign(hh, dcel_pair->second); + CGAL_assertion(check); + boost::optional< int > m = + hh->data()->multiplicity_of_cut( + surface1, surface2 + ); + CGAL_assertion(m); + if (*m == 1) { + std::pair< int, int > intersection = + iv_helper.unique_overlap( + isolator1, isolator2, + sit->second.second.begin(), + sit->second.second.end()); + if (intersection.first == root1 && + intersection.second == root2) { + //std::cout << "filterD4" << std::endl; + surely_equal = true; + surely_unequal = false; + } else { + //std::cout << "filterD5" << std::endl; + surely_equal = false; + surely_unequal = true; + } + } + } + } + } + + +#if 0 + // empty filter + if (!surely_equal && !surely_unequal && still_possible_equal) { + + } +#endif + + // final check (costly one!) + if (!surely_equal && !surely_unequal && still_possible_equal) { + // reaching here, no filter applied + + // what are the possibilities: + // (1) vertex and at least one silhouette + CGAL_assertion_code( + bool poss1 = ( + feature == VERTEX && + (projected_sil1 || projected_sil2) + ); + + ); + // (2) EDGE with mult > 1 + CGAL_assertion_code( + bool poss2 = ( + feature == EDGE && + (mult > 1) + ); + ); + // TODO (2a) an isolated segment? Something special? + CGAL_assertion_code( + bool poss2a = false; + ); + CGAL_assertion(poss1 || poss2 || poss2a); + + // Remark: It is CRUCIAL that we can decide (somehow) + // whether two roots are equal + typename Surface_z_at_xy_isolator_traits::Equal_z equal_z + // TODO (object isolator1.equal_z_object()) + ; + + // we have to ask the traits class using equal_z + // remove mult as upper bound of possible intersections + surely_equal = equal_z( + surface1, isolator1, root1, projected_sil1, + surface2, isolator2, root2, projected_sil2); + surely_unequal = !surely_equal; + } + + // final result + if (surely_equal) { + CGAL_assertion(!surely_unequal); + // reach here, it's not only possible equal -> it is equal! + + // update z-id + it->_make_unique_for(*this, z); + + // insert it to current cell + it->_add_sheet(surface2, root2, isolator2.length(root2)); + + // proceed to next root and next cell! + it++; + root2++; + z++; + + // the next combination can be equal again + still_possible_equal = true; + + } else { + CGAL_assertion(!still_possible_equal || surely_unequal); + CGAL_assertion(!surely_equal); + + // otherwise, we refine interval + isolator1.refine_interval(root1); + isolator2.refine_interval(root2); + + // the next combination is equal to this one, + // which is cannot be equal! + still_possible_equal = false; + } + break; + } + case CGAL::LARGER: { + //std::cout << "LARGER root1: " << root1 << " root2: " << root2 + // << std::endl; + // create new Z_cell + Z_cell cell(*this, z); + cell._add_sheet(surface2, root2, isolator2.length(root2)); + + // and insert it before current position + this->ptr()->_m_z_cells.insert(it, cell); + + // go to next root + root2++; + z++; + + // the next combination can be equal again + still_possible_equal = true; + + break; + } + } + } + while (it != this->ptr()->_m_z_cells.end()) { + it->_make_unique_for(*this, z); + + // go to next existing cell in z_stack + it++; + z++; + } + while (root2 < roots2) { + // append cells of isolator2 at the end + Z_cell cell(*this,z); + cell._add_sheet(surface2, root2, isolator2.length(root2)); + + // and insert it at the end + this->ptr()->_m_z_cells.push_back(cell); + + // go to next root + root2++; + z++; + } + + } + + //!@} + +private: + //!\name Merging + //!@{ + + /*!\brief + * merges \c *this and \c z_stack to a new z_stack + */ + Self _merge(const Self& z_stack, const Dcel_data* data) { + Self tmp(z_stack); + tmp.copy_on_write(); + tmp.ptr()->_m_dcel_data = data; + + for (typename Z_at_xy_isolators::iterator it = + this->ptr()->_m_z_at_xy_isolators.begin(); + it != this->ptr()->_m_z_at_xy_isolators.end(); it++) { + tmp._add_surface(it->first, it->second); + } + + for (typename Z_at_xy_isolators::iterator it = + this->ptr()->_m_z_at_xy_empty_isolators.begin(); + it != this->ptr()->_m_z_at_xy_empty_isolators.end(); it++) { + tmp._add_surface(it->first, it->second); + } + + return tmp; + } + + //!@} + +public: + //!\name Information + //!@{ + + /*\brief + * returns one of stored projected point for z_stack + */ + Point_2 point() const { + return this->ptr()->_m_point; + } + + /*!\brief + * returns the number of z_stack cells + */ + bool is_empty() const { + return this->ptr()->_m_z_cells.empty(); + } + + /*!\brief + * returns the number of z_stack cells + */ + int number_of_z_cells() const { + return static_cast< int >(this->ptr()->_m_z_cells.size()); + } + + //!@} + + //!\name Accessors + //!@{ + + /*!\brief + * returns \c i-th z_stack cell + */ + Z_cell z_cell(int i) const { + CGAL_precondition(i >= 0 && i < this->number_of_z_cells()); + typename Z_cell_container::const_iterator it = + this->ptr()->_m_z_cells.begin(); + std::advance(it, i); + return *it; + } + + /*\brief + * return level of \c surface in \c i-th z-cell + */ + int level_of_surface_in_z_cell(const Surface_3& surface, int i) const { + CGAL_precondition(i >= 0); + CGAL_precondition(i < this->number_of_z_cells()); + const Z_cell& cell = this->z_cell(i); + return cell.sheet_number(surface); + } + + /*!\brief + * beginning of z-cells + */ + typename Z_cell_container::iterator z_cells_begin() { + return this->ptr()->_m_z_cells.begin(); + } + + /*!\brief + * past-the-end of z-cells + */ + typename Z_cell_container::iterator z_cells_end() { + return this->ptr()->_m_z_cells.end(); + } + + /*!\brief + * beginning of z-cells + */ + Z_cell_const_iterator z_cells_begin() const { + return this->ptr()->_m_z_cells.begin(); + } + + /*!\brief + * past-the-end of z-cells + */ + Z_cell_const_iterator z_cells_end() const { + return this->ptr()->_m_z_cells.end(); + } + + /*\brief + * returns all levels of \c surface + */ + template < class OutputIterator > + OutputIterator levels_of( + const Surface_3& surface, + OutputIterator oi + ) const { + int level = 0; + for (typename Z_cell_container::const_iterator it = + this->z_cells_begin(); it != this->z_cells_end(); + it++, level++) { + if (it->sheet_number(surface) != -1) { + *oi++ = level; + } + } + return oi; + } + + + /*!\brief + * returns z-level for a given surface sheet + */ + int z_level_of_sheet(const Surface_3& surface, int sheet) const { + CGAL_precondition(sheet < this->number_of_z_cells()); + // it suffices to start at "sheet" + for (int i = 0; i < this->number_of_z_cells(); i++) { + if (this->level_of_surface_in_z_cell(surface, i) == sheet) { + return i; + } + } + CGAL_error_msg("Not allowed to reach here"); + /* Should not be reached */ return -1; + } + +#if 0 + // TASK implement other intersections predicates + /*\brief + * returns all pairs of intersection \c surface1 and \c surface2 + * are involved + */ + template < class OutputIterator > + OutputIterator intersections( + const Surface_3& surface1, const Surface_3& surface2, + OutputIterator oi + ) const { + typedef Z_cell::value_type Surface_id; + typedef std::pair< Surface_id, Surface_id > Surface_id_pair; + + for (typename Z_cell_container::const_iterator it = + z_cells_begin(); it != z_cells.end(); it++) { + Surface_id si1 = it->find(surface1); + Surface_id si2 = it->find(surface2); + if (si1 != it->end() && si2 != it->end()) { + *oi++ = std::make_pair(si1, si2); + } + } + return oi; + } + + /*\brief + * returns all pairs of intersection \c surface1 and \c surface2 + * are involved + */ + template < class OutputIterator > + OutputIterator intersections( + const Surface_3& surface1, + OutputIterator oi + ) const { + typedef Z_cell::value_type Surface_id; + typedef std::pair< Surface_id, Surface_id > Surface_id_pair; + + for (typename Z_cell_container::const_iterator it = + z_cells_begin(); it != z_cells.end(); it++) { + Surface_id si1 = it->find(surface1); + Surface_id si2 = it->find(surface2); + if (si1 != it->end()) { + for (typename Z_cell::const_iterator sit = + it->begin(); sit != it->end(); sit++) { + if (sit != si1) { + *oi++ = std::make_pair(si1, *sit); + } + } + } + } + return oi; + } + +#endif + + /*\brief + * returns all levels of \c surface that is involved in some intersection + */ + template < class OutputIterator > + OutputIterator intersection_levels_of( + const Surface_3& surface, + OutputIterator oi + ) const { + if (this->ptr()->_m_dcel_data->_has_cut()) { + int level = 0; + for (typename Z_cell_container::const_iterator it = + this->z_cells_begin(); it != this->z_cells_end(); + it++, level++) { + if (it->number_of_sheets() > 1) { + int sheet = it->sheet_number(surface); + if (sheet != -1) { + *oi++ = sheet; + } + } + } + } + return oi; + } + + + //!@} + +public: + //!\name Adjacency + //!@{ + + /*!\brief + * Compute adjacency information when going from isolator \c *this and + * sheet defined by \c surface and \c sheet to \c to. Return the final + * sheet number at \c to. + */ + std::pair< int, int > adjacency(const Surface_3& surface, int sheet, + SoX::Dcel_feature feat_from, + CGAL::Object dcel_handle_from, + Self to, + SoX::Dcel_feature feat_to, + CGAL::Object dcel_handle_to + ) const { + typename Z_at_xy_isolators::const_iterator thisit = + this->ptr()->_m_z_at_xy_isolators.find(surface); + CGAL_assertion(thisit != this->ptr()->_m_z_at_xy_isolators.end()); + + typename Z_at_xy_isolators::const_iterator thatit = + to.ptr()->_m_z_at_xy_isolators.find(surface); + if (thatit == to.ptr()->_m_z_at_xy_isolators.end()) { + // if surface is not involved in "that" dcel feature + // -> either not connected or infinity + thatit = to.ptr()->_m_z_at_xy_empty_isolators.find(surface); + CGAL_assertion( + thatit != to.ptr()->_m_z_at_xy_empty_isolators.end() + ); + } + + typedef typename Surface_z_at_xy_isolator_traits::Adjacency Adjacency; + + typedef typename SoX::Adjacencies_3::Adjacency_interval + Adjacency_interval; + + if (this->id() == to.id()) { + return std::make_pair(sheet, sheet); + } + // else + + typedef Restricted_cad_3< Surface_z_at_xy_isolator_traits > + Restricted_cad_3; + + CGAL_assertion_code( + Restricted_cad_3 cad = Restricted_cad_3::cad_cache()(surface); + ); + + bool has_vertical_line_from = false; + if (feat_from == SoX::VERTEX) { + typename Restricted_cad_3::Vertex_const_handle vh; + CGAL_assertion_code(bool check = ) + CGAL::assign(vh,dcel_handle_from); + CGAL_assertion(check); + has_vertical_line_from = + vh->data()->_supports_vertical_line(surface); + + CGAL_assertion(vh->data()->_rs_id() == cad.id()); + } + + CGAL_assertion_code(( + { + if (feat_from == SoX::EDGE) { + typename Restricted_cad_3::Halfedge_const_handle hh1, hh2; + CGAL_assertion_code(CGAL::assign(hh1, dcel_handle_from)); + CGAL_assertion(hh1->data()->_rs_id() == cad.id()); + CGAL::Object obj = + cad.locate(thisit->second.traits().point()); + CGAL_assertion(CGAL::assign(hh2, obj)); + CGAL_assertion(hh1 == hh2); + } + }) + ); + CGAL_assertion_code(( + { + if (feat_from == SoX::FACE) { + typename Restricted_cad_3::Face_const_handle fh1, fh2; + CGAL_assertion(CGAL::assign(fh1, dcel_handle_from)); + CGAL_assertion(fh1->data()->_rs_id() == cad.id()); + CGAL::Object obj = + cad.locate(thisit->second.traits().point()); + CGAL_assertion(CGAL::assign(fh2, obj)); + CGAL_assertion(fh1 == fh2); + } + }) + ); + + bool has_vertical_line_to = false; + if (feat_to == SoX::VERTEX) { + typename Restricted_cad_3::Vertex_const_handle vh; + CGAL_assertion_code(bool check = ) + CGAL::assign(vh,dcel_handle_to); + CGAL_assertion(check); + has_vertical_line_to = + vh->data()->_supports_vertical_line(surface); + + CGAL_assertion(vh->data()->_rs_id() == cad.id()); + } + + CGAL_assertion_code(( + { + if (feat_to == SoX::EDGE) { + typename Restricted_cad_3::Halfedge_const_handle hh1, hh2; + CGAL_assertion_code(CGAL::assign(hh1, dcel_handle_to)); + CGAL_assertion(hh1->data()->_rs_id() == cad.id()); + CGAL::Object obj = + cad.locate(thatit->second.traits().point()); + CGAL_assertion(CGAL::assign(hh2, obj)); + CGAL_assertion(hh1 == hh2 || hh1->twin() == hh2); + } + }) + ); + CGAL_assertion_code(( + { + if (feat_to == SoX::FACE) { + typename Restricted_cad_3::Face_const_handle fh1, fh2; + CGAL_assertion(CGAL::assign(fh1, dcel_handle_to)); + CGAL_assertion(fh1->data()->_rs_id() == cad.id()); + CGAL::Object obj = + cad.locate(thatit->second.traits().point()); + CGAL_assertion(CGAL::assign(fh2, obj)); + CGAL_assertion(fh1 == fh2); + } + }) + ); + + + Adjacency compute_adj; // TODO single instance + SoX::Adjacencies_3 adj = + compute_adj(surface, + thisit->second, feat_from, dcel_handle_from, + has_vertical_line_from, + thatit->second, feat_to, dcel_handle_to, + has_vertical_line_to); + Adjacency_interval intv = adj.interval(sheet); + return intv; + } + + //!@} + + //!\name IO + //!@{ + + /*!\brief + * prints pretty-formated version of z_stack + */ + void pretty_print(std::ostream& os) const { + os << "#ZS @ "; + os << point() << ":" << std::endl; + os << "[" << number_of_z_cells() << ","; + for (typename Z_cell_container::const_iterator sit = + z_cells_begin(); sit != z_cells_end(); sit++) { + sit->pretty_print(os); + } + os << " id=" << this->id() << ">"; + os << "]"; + + } + + //!@} + + // friend + friend class SoX::P_dcel_info< Surface_z_at_xy_isolator_traits >; + // for _isolator + friend class + SoX::Intern::Z_cell< Surface_z_at_xy_isolator_traits, Dcel_data >; + friend class SoX::Restricted_cad_3< Surface_z_at_xy_isolator_traits >; +}; + + +/*!\relates Z_stack + * \brief + * outputs Z_stack object to stream \c os + */ +template < class SurfaceZAtXyIsolatorTraits, class DcelData > +std::ostream& operator<<( + std::ostream& os, + const Z_stack< SurfaceZAtXyIsolatorTraits, DcelData >& z_stack +) { + z_stack.pretty_print(os); + return os; +} + + +} // namespace SoX + +#endif // SoX_GAPS_Z_STACK_H +// EOF diff --git a/Curved_kernel_via_analysis_2/include/CGAL/Arrangement_2l/Z_stack_helpers.h b/Curved_kernel_via_analysis_2/include/CGAL/Arrangement_2l/Z_stack_helpers.h new file mode 100644 index 00000000000..5cc566477e0 --- /dev/null +++ b/Curved_kernel_via_analysis_2/include/CGAL/Arrangement_2l/Z_stack_helpers.h @@ -0,0 +1,309 @@ +// ============================================================================ +// +// Copyright (c) 2001-2008 Max-Planck-Institut Saarbruecken (Germany). +// All rights reserved. +// +// This file is part of EXACUS (http://www.mpi-inf.mpg.de/projects/EXACUS/); +// you may redistribute it under the terms of the Q Public License version 1.0. +// See the file LICENSE.QPL distributed with EXACUS. +// +// Licensees holding a valid commercial license may use this file in +// accordance with the commercial license agreement provided with the software. +// +// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE +// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. +// +// ---------------------------------------------------------------------------- +// +// Library : SoX +// File : include/SoX/GAPS/Z_stack_helpers.h +// SoX_release : $Name: $ +// Revision : $Revision$ +// Revision_date : $Date$ +// +// Author(s) : Eric Berberich +// +// ============================================================================ + +/*! \file SoX/GAPS/Z_stack_helpers.h + * \brief definition of Z_stack helpers + */ + +#ifndef SoX_GAPS_Z_STACK_HELPERS_H +#define SoX_GAPS_Z_STACK_HELPERS_H 1 + +#include + +namespace SoX { + +namespace Intern { + +/*!\brief + * Functor to determine the order and overlaps of two realrootisolators + */ +template < class RealRootIsolator > +class Refineable_interval_helper { +public: + //! this instance's template parameter + typedef RealRootIsolator Real_root_isolator; + + /*!\brief + * Compute order of two refineable intervals given by + * pair \c isolator1 and \c i1 and the pair + * \c isolator2 and \c i2 + * + * Returns \c CGAL::SMALLER if \c i1 -th interval of \c isolator1 + * is strictly less then the other, CGAL::LARGER, if strictly + * greater, and CGAL::EQUAL in case of an overlap of the two + * given intervals. + */ + ::CGAL::Comparison_result overlap_or_order( + const Real_root_isolator& isolator1, int i1, + const Real_root_isolator& isolator2, int i2) const { + + if (isolator1.right_boundary(i1) < isolator2.left_boundary(i2)) { + return CGAL::SMALLER; + } + if (isolator2.right_boundary(i2) < isolator1.left_boundary(i1)) { + return CGAL::LARGER; + } + return CGAL::EQUAL; + } + + /*!\brief + * Return whether interval defined by the pair \c isolator1 and \c i1 + * is included in intervals defined by pair \c isolator2 and \c i2 + * + * \c strict = true indicates a inclusion in interior + */ + bool is_included( + const Real_root_isolator& isolator1, int i1, + const Real_root_isolator& isolator2, int i2, + bool strict = false) const { + + if (strict) { + return + (isolator1.left_boundary(i1) > isolator2.left_boundary(i2)) && + (isolator1.right_boundary(i1) < isolator2.right_boundary(i2)); + } + // else + return + (isolator1.left_boundary(i1) >= isolator2.left_boundary(i2)) && + (isolator1.right_boundary(i1) <= isolator2.right_boundary(i2)); + } + + + /*!\brief + * computs all pairwise refined overlapping id + * + * value_type of OutputIterator is std::pair< int, int > + */ + template < class OutputIterator > + OutputIterator refined_overlaps( + const Real_root_isolator& isolator1, + const Real_root_isolator& isolator2, + OutputIterator oi + ) const { + + // obtain minimal number of overlaps, i.e., all each interval + // of one isolator overlaps with at most one interval of the + // other isolator + + const int num1 = isolator1.number_of_real_roots(); + const int num2 = isolator2.number_of_real_roots(); + CGAL_precondition(num1 > 0); + CGAL_precondition(num2 > 0); + + // TASK cache result for two isolators? + + // compute all overlaps + + std::list< int > empty; + + std::vector< std::list< int > > ovl1(num1, empty); + std::vector< std::list< int > > ovl2(num2, empty); + + int i1 = 0; + int i2 = 0; + while (true) { + if (i1 == num1) { + // no further overlaps possible + break; + } + if (i2 == num2) { + // no further overlaps possible + break; + } + CGAL::Comparison_result res = + this->overlap_or_order(isolator1, i1, isolator2, i2); + + switch (res) { + case CGAL::SMALLER: + //std::cout << "SM" << std::endl; + i1++; + break; + case CGAL::EQUAL: + //std::cout << "EQUAL " << i1 << " " << i2 << std::endl; + ovl1[i1].push_back(i2); + ovl2[i2].push_back(i1); + if (isolator1.right_boundary(i1) <= + isolator2.right_boundary(i2)) { + i1++; + } else if (isolator1.right_boundary(i1) >= + isolator2.right_boundary(i2)) { + i2++; + } + break; + case CGAL::LARGER: + //std::cout << "LA" << std::endl; + i2++; + break; + } + } + + for (i1 = 0; i1 < num1; i1++) { + switch (ovl1[i1].size()) { + case 0: + // nothing to do + case 1: + // is also fine + break; + default: + // size is greater than 1, refine against all found interval + do { + std::list< typename std::list< int >::iterator > erase; + for (typename std::list< int >::iterator it = + ovl1[i1].begin(); + it != ovl1[i1].end(); it++) { + isolator1.refine_interval(i1); + isolator2.refine_interval(*it); + if (this->overlap_or_order(isolator1, i1, + isolator2, *it) != + CGAL::EQUAL) { + erase.push_back(it); + } + } + for (typename + std::list< + typename std::list< int >::iterator >::iterator + eit = erase.begin(); eit != erase.end(); eit++) { + ovl2[*(*eit)].remove(i1); + ovl1[i1].erase(*eit); + } + } while (static_cast< int >(ovl1[i1].size()) > 1); + } + } + + for (i2 = 0; i2 < num2; i2++) { + switch (ovl2[i2].size()) { + case 0: + // nothing to do + case 1: + // is also fine + break; + default: + // size is greater than 1, refine against all found interval + do { + std::list< typename std::list< int >::iterator > erase; + for (typename std::list< int >::iterator it = + ovl2[i2].begin(); + it != ovl2[i2].end(); it++) { + isolator1.refine_interval(*it); + isolator2.refine_interval(i2); + if (this->overlap_or_order(isolator1, *it, + isolator2, i2) != + CGAL::EQUAL) { + erase.push_back(it); + } + } + for (typename + std::list< typename std::list< int >::iterator >:: + iterator + eit = erase.begin(); eit != erase.end(); eit++) { + ovl1[*(*eit)].remove(i2); + ovl2[i2].erase(*eit); + } + } while (static_cast< int >(ovl2[i2].size()) > 1); + } + } + + for (i1 = 0; i1 < num1; i1++) { + switch (ovl1[i1].size()) { + case 0: + // nothing to do + break; + case 1: { + // check whether it exist in other + int i2 = ovl1[i1].front(); + CGAL_assertion(static_cast< int >(ovl2[i2].size()) == 1); + CGAL_assertion(ovl2[i2].front() == i1); + *oi++ = std::make_pair(i1,i2); + break; + } + default: + CGAL_error_msg("should not occur"); + } + } + + return oi; + } + + //! returns the unique overlapping pair of intervals of two isolators + template< class InputIterator > + std::pair< int, int > unique_overlap( + const Real_root_isolator& isolator1, + const Real_root_isolator& isolator2, + InputIterator begin, + InputIterator end + ) const { + + // obtain the unique overlap by further refinement, as we know + // from the external that is is exactly one overlap + std::list< std::pair< int, int > > overlaps; + std::copy(begin,end, std::back_inserter(overlaps)); + + CGAL_precondition(static_cast< int >(overlaps.size()) > 0); + + while (static_cast< int >(overlaps.size()) > 1) { + std::list< typename std::list< std::pair< int, int > >::iterator > + erase; + + // for each overlapping pair + for (typename std::list< std::pair< int, int > >::iterator + it = overlaps.begin(); it != overlaps.end(); it++) { + // refine overlapping intervals + isolator1.refine_interval(it->first); + isolator2.refine_interval(it->second); + + // and remove non-overlapping pairs + if (this->overlap_or_order( + isolator1, it->first, isolator2, it->second + ) != CGAL::EQUAL) { + erase.push_back(it); + } + } + for (typename + std::list< typename std::list< + std::pair< int, int > >::iterator >::iterator it = + erase.begin(); it != erase.end(); it++) { + overlaps.erase(*it); + } + } + CGAL_assertion(static_cast< int >(overlaps.size()) == 1); + + CGAL_postcondition(overlaps.begin()->first >= 0); + CGAL_postcondition(overlaps.begin()->second >= 0); + return std::make_pair( + overlaps.begin()->first, overlaps.begin()->second + ); + } + +}; // Refineable_interval_helper + + +} // namespace Intern + +} // namespace SoX + +#endif // SoX_GAPS_Z_STACK_HELPERS_H +// EOF diff --git a/Curved_kernel_via_analysis_2/include/CGAL/Arrangement_2l/macros.h b/Curved_kernel_via_analysis_2/include/CGAL/Arrangement_2l/macros.h new file mode 100644 index 00000000000..7d9c7485e3a --- /dev/null +++ b/Curved_kernel_via_analysis_2/include/CGAL/Arrangement_2l/macros.h @@ -0,0 +1,52 @@ +// ============================================================================ +// +// Copyright (c) 2001-2008 Max-Planck-Institut Saarbruecken (Germany). +// All rights reserved. +// +// This file is part of EXACUS (http://www.mpi-inf.mpg.de/projects/EXACUS/); +// you may redistribute it under the terms of the Q Public License version 1.0. +// See the file LICENSE.QPL distributed with EXACUS. +// +// Licensees holding a valid commercial license may use this file in +// accordance with the commercial license agreement provided with the software. +// +// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE +// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. +// +// ---------------------------------------------------------------------------- +// +// Library : SoX +// File : include/SoX/GAPS/macros.h +// SoX_release : $Name: $ +// Revision : $Revision$ +// Revision_date : $Date$ +// +// Author(s) : Eric Berberich +// +// ============================================================================ + +/*! \file SoX/GAPS/macros.h + \brief some basic macro declarations +*/ + +#ifndef SoX_GAPS_MACROS_H +#define SoX_GAPS_MACROS_H 1 + +#include + +// TASK document +#define SoX_SURFACE_Z_AT_XY_ISOLATOR_TRAITS_SNAP_TYPEDEFS(Traits) \ + typedef typename Traits::Surface_3 Surface_3; \ + typedef typename Traits::Arrangement_traits_2 Arrangement_traits_2; \ + typedef typename Arrangement_traits_2::Curve_kernel_2 Curve_kernel_2; \ + typedef typename Traits::Z_at_xy_isolator Z_at_xy_isolator;\ + typedef typename Traits::Polynomial_3 Polynomial_3; \ + typedef typename Traits::Polynomial_2 Polynomial_2; \ + typedef typename Curve_kernel_2::Curve_analysis_2 Curve_analysis_2; \ + typedef typename Traits::Point_2 Point_2; \ + typedef typename Traits::X_monotone_curve_2 X_monotone_curve_2; \ +// end #define SoX_SURFACE_Z_AT_XY_ISOLATOR_TRAITS_SNAP_TYPEDEFS(Traits) + + +#endif // SoX_GAPS_MACROS_H +// EOF diff --git a/Curved_kernel_via_analysis_2/include/CGAL/Arrangement_2l/z_stack_predeclarations.h b/Curved_kernel_via_analysis_2/include/CGAL/Arrangement_2l/z_stack_predeclarations.h new file mode 100644 index 00000000000..2607f4b1515 --- /dev/null +++ b/Curved_kernel_via_analysis_2/include/CGAL/Arrangement_2l/z_stack_predeclarations.h @@ -0,0 +1,82 @@ +// ============================================================================ +// +// Copyright (c) 2001-2008 Max-Planck-Institut Saarbruecken (Germany). +// All rights reserved. +// +// This file is part of EXACUS (http://www.mpi-inf.mpg.de/projects/EXACUS/); +// you may redistribute it under the terms of the Q Public License version 1.0. +// See the file LICENSE.QPL distributed with EXACUS. +// +// Licensees holding a valid commercial license may use this file in +// accordance with the commercial license agreement provided with the software. +// +// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE +// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. +// +// ---------------------------------------------------------------------------- +// +// Library : SoX +// File : include/SoX/GAPS/z_stack_predclarations.h +// SoX_release : $Name: $ +// Revision : $Revision$ +// Revision_date : $Date$ +// +// Author(s) : Eric Berberich +// +// ============================================================================ + +/*! \file SoX/GAPS/z_stack_predeclarations.h + \brief contains all required pre-declarations +*/ + +#ifndef SoX_GAPS_Z_STACK_PREDECLARATIONS_H +#define SoX_GAPS_Z_STACK_PREDECLARATIONS_H 1 + +#include + +namespace SoX { + +// pre-declaration +template < class SurfaceZAtXyIsolatorTraits > +class Restricted_cad_3; + +// pre-declaration +template < class SurfaceZAtXyIsolatorTraits > +class P_dcel_info; + +// pre-declaration +template < class SurfaceZAtXyIsolatorTraits, class DcelData > +class Z_stack; + +// pre-declaration +template < class SurfaceZAtXyIsolatorTraits > +class Create_restricted_cad_3; + +// pre-declaration +template < class SurfaceZAtXyIsolatorTraits > +class Overlay_restricted_cad_3; + +// pre-declaration +template < class Arrangement_2_ > +class Arr_p_dcel_info_overlay_traits; + +// pre-declaration +template < class SurfaceZAtXyIsolatorTraits, class Rep_ > +class Surface_pair_3; + +namespace Intern { + +// pre-declaration +template < class SurfaceZAtXyIsolatorTraits > +class Restricted_cad_3_cache; + +// pre-declaration +template < class SurfaceZAtXyIsolatorTraits, class DcelData > +class Z_cell; + +} // namespace Intern + +} // namespace SoX + +#endif // SoX_GAPS_Z_STACK_PREDECLARATIONS_H +// EOF