include/CGAL/Convex_hull_d.h

// Copyright (c) 1997-2000  Max-Planck-Institute Saarbruecken (Germany).
// All rights reserved.
//
// This file is part of CGAL (www.cgal.org).
//
// $URL: https://github.com/CGAL/cgal/blob/v5.3/Convex_hull_d/include/CGAL/Convex_hull_d.h $
// $Id: Convex_hull_d.h 0d66e19 2020-07-24T17:05:10+02:00 Mael Rouxel-Labbé
// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial
//
//
// Author(s)     : Michael Seel <seel@mpi-sb.mpg.de>
//---------------------------------------------------------------------
// file generated by notangle from Convex_hull_d.lw
// please debug or modify web file
// mails and bugs: Michael.Seel@mpi-sb.mpg.de
// based on LEDA architecture by S. Naeher, C. Uhrig
// coding: K. Mehlhorn, M. Seel
// debugging and templatization: M. Seel
//---------------------------------------------------------------------

#ifndef CGAL_CONVEX_HULL_D_H
#define CGAL_CONVEX_HULL_D_H

#include <CGAL/license/Convex_hull_d.h>

#define CGAL_DEPRECATED_HEADER "<CGAL/Convex_hull_d.h>"
#define CGAL_DEPRECATED_MESSAGE_DETAILS \
  "The Triangulation package (see https://doc.cgal.org/latest/Triangulation) should be used instead."
#include <CGAL/internal/deprecation_warning.h>

/*{\Manpage {Convex_hull_d}{R}{Convex Hulls}{C}}*/

/*{\Mdefinition An instance |\Mvar| of type |\Mname| is the convex
hull of a multi-set |S| of points in $d$-dimensional space. We call
|S| the underlying point set and $d$ or |dim| the dimension of the
underlying space. We use |dcur| to denote the affine dimension of |S|.
The data type supports incremental construction of hulls.

The closure of the hull is maintained as a simplicial complex, i.e.,
as a collection of simplices. The intersection of any two is a face of
both\cgalFootnote{The empty set if a facet of every simplex.}. In the
sequel we reserve the word simplex for the simplices of dimension
|dcur|. For each simplex there is a handle of type |Simplex_handlex|
and for each vertex there is a handle of type |Vertex_handle|.  Each
simplex has $1 + |dcur|$ vertices indexed from $0$ to |dcur|; for a
simplex $s$ and an index $i$, |C.vertex(s,i)| returns the $i$-th
vertex of $s$. For any simplex $s$ and any index $i$ of $s$ there may
or may not be a simplex $t$ opposite to the $i$-th vertex of $s$.  The
function |C.opposite_simplex(s,i)| returns $t$ if it exists and
returns |Simplex_handle()| (the undefined handle) otherwise. If $t$
exists then $s$ and $t$ share |dcur| vertices, namely all but the
vertex with index $i$ of $s$ and the vertex with index
|C.index_of_vertex_in_opposite_simplex(s,i)| of $t$. Assume that $t$
exists and let |j = C.index_of_vertex_in_opposite_simplex(s,i)|.  Then
$s =$ |C.opposite_simplex(t,j)| and $i =$
|C.index_of_vertex_in_opposite_simplex(t,j)|.

The boundary of the hull is also a simplicial complex. All simplices
in this complex have dimension $|dcur| - 1$.  For each boundary simplex
there is a handle of type |Facet_handle|.  Each facet has |dcur| vertices
indexed from $0$ to $|dcur| - 1$. If |dcur > 1| then for each facet $f$
and each index $i$, $0 \le i < |dcur|$, there is a facet $g$ opposite
to the $i$-th vertex of $f$.  The function |C.opposite_facet(f,i)|
returns $g$.  Two neighboring facets $f$ and $g$ share |dcur - 1|
vertices, namely all but the vertex with index $i$ of $f$ and the
vertex with index |C.index_of_vertex_in_opposite_facet(f,i)| of $g$.
Let |j = C.index_of_vertex_in_opposite_facet(f,i)|. Then
|f = C.opposite_facet(g,j)| and
|i = C.index_of_vertex_in_opposite_facet(g,j)|.}*/

#include <CGAL/basic.h>
#include <CGAL/Origin.h>
#include <CGAL/Unique_hash_map.h>
#include <CGAL/Regular_complex_d.h>
#include <CGAL/Handle_for.h>
#include <list>
#include <vector>

#include <CGAL/Kernel_d/debug.h>

namespace CGAL {

template <typename HP, typename H> class Facet_iterator_rep_;
template <typename HP, typename H> class Facet_iterator_;

template <typename Hull_pointer, typename Handle>
class Facet_iterator_rep_
{
  CGAL::Unique_hash_map<Handle,bool>* pvisited_;
  std::list<Handle>* pcandidates_;
  Hull_pointer hull_;
  Handle current_;
  friend class Facet_iterator_<Hull_pointer,Handle>;

  void add_candidates()
  { CGAL_assertion(pvisited_ && pcandidates_ && hull_);
    for(int i = 1; i <= hull_->current_dimension(); ++i) {
      Handle f = hull_->opposite_simplex(current_,i);
      if ( !(*pvisited_)[f] ) {
        pcandidates_->push_front(f);
        (*pvisited_)[f] = true;
      }
    }
  }
public:
  Facet_iterator_rep_() :
    pvisited_(0), pcandidates_(0), hull_(0), current_() {}
  Facet_iterator_rep_(Hull_pointer p, Handle h) :
    pvisited_(0), pcandidates_(0), hull_(p), current_(h) {}
  ~Facet_iterator_rep_()
  { if (pvisited_) delete pvisited_;
    if (pcandidates_) delete pcandidates_; }
};

template <typename Hull_pointer, typename Handle>
class Facet_iterator_ : private
  Handle_for<Facet_iterator_rep_<Hull_pointer,Handle> >
{ typedef Facet_iterator_<Hull_pointer,Handle> Self;
  typedef Facet_iterator_rep_<Hull_pointer,Handle> Rep;
  typedef Handle_for< Facet_iterator_rep_<Hull_pointer,Handle> > Base;

  using Base::ptr;

public:
  typedef typename Handle::value_type value_type;
  typedef typename Handle::pointer    pointer;
  typedef typename Handle::reference  reference;
  typedef typename Handle::difference_type difference_type;
  typedef std::forward_iterator_tag   iterator_category;

  Facet_iterator_() : Base( Rep() ) {}
  Facet_iterator_(Hull_pointer p, Handle h) : Base( Rep(p,h) )
  { ptr()->pvisited_ = new Unique_hash_map<Handle,bool>(false);
    ptr()->pcandidates_ = new std::list<Handle>;
    (*(ptr()->pvisited_))[ptr()->current_]=true;
    ptr()->add_candidates();
  }

  reference operator*() const
  { return ptr()->current_.operator*(); }
  pointer operator->() const
  { return ptr()->current_.operator->(); }

  Self& operator++()
  { if ( ptr()->current_ == Handle() ) return *this;
    if ( ptr()->pcandidates_->empty() ) ptr()->current_ = Handle();
    else {
      ptr()->current_ = ptr()->pcandidates_->back();
      ptr()->pcandidates_->pop_back();
      ptr()->add_candidates();
    }
    return *this;
  }
  Self  operator++(int)
  { Self tmp = *this; ++*this; return tmp; }
  bool operator==(const Self& x) const
  { return ptr()->current_ == x.ptr()->current_; }
  bool operator!=(const Self& x) const
  { return !operator==(x); }

  operator Handle() { return ptr()->current_; }

};


template <typename HP, typename VH, typename FH>
class Hull_vertex_iterator_rep_;
template <typename HP, typename VH, typename FH>
class Hull_vertex_iterator_;

template <typename Hull_pointer, typename VHandle, typename FHandle>
class Hull_vertex_iterator_rep_
{
  CGAL::Unique_hash_map<VHandle,bool>* pvisited_;
  Hull_pointer hull_;
  VHandle v_; FHandle f_; int i_;
  friend class Hull_vertex_iterator_<Hull_pointer,VHandle,FHandle>;

  void advance()
  { CGAL_assertion(pvisited_ &&  hull_);
    if ( f_ == FHandle() ) return;
    bool search_next = true; ++i_;
    while ( search_next ) {
      while(i_ < hull_->current_dimension()) {
        v_ = hull_->vertex_of_facet(f_,i_);
        if ( !(*pvisited_)[v_] ) { search_next=false; break; }
        ++i_;
      }
      if ( search_next ) { i_=0; ++f_; }
      if ( f_ == FHandle() )
      { search_next=false; v_ = VHandle(); }
    }
    (*pvisited_)[v_] = true;
  }
public:
  Hull_vertex_iterator_rep_() :
    pvisited_(0), hull_(0), i_(0) {}
  Hull_vertex_iterator_rep_(Hull_pointer p, FHandle f) :
    pvisited_(0), hull_(p), f_(f), i_(-1) {}
  ~Hull_vertex_iterator_rep_()
  { if (pvisited_) delete pvisited_; }
};

template <typename Hull_pointer, typename VHandle, typename FHandle>
class Hull_vertex_iterator_ : private
  Handle_for< Hull_vertex_iterator_rep_<Hull_pointer,VHandle,FHandle> >
{ typedef Hull_vertex_iterator_<Hull_pointer,VHandle,FHandle> Self;
  typedef Hull_vertex_iterator_rep_<Hull_pointer,VHandle,FHandle> Rep;
  typedef Handle_for< Rep > Base;

  using Base::ptr;

public:
  typedef typename VHandle::value_type value_type;
  typedef typename VHandle::pointer    pointer;
  typedef typename VHandle::reference  reference;
  typedef typename VHandle::difference_type difference_type;
  typedef std::forward_iterator_tag   iterator_category;

  Hull_vertex_iterator_() : Base( Rep() ) {}
  Hull_vertex_iterator_(Hull_pointer p, FHandle f) : Base( Rep(p,f) )
  { ptr()->pvisited_ = new Unique_hash_map<VHandle,bool>(false);
    ptr()->advance();
  }

  reference operator*() const
  { return ptr()->v_.operator*(); }
  pointer operator->() const
  { return ptr()->v_.operator->(); }

  Self& operator++()
  { if ( ptr()->v_ == VHandle() ) return *this;
    ptr()->advance();
    return *this;
  }
  Self  operator++(int)
  { Self tmp = *this; ++*this; return tmp; }
  bool operator==(const Self& x) const
  { return ptr()->v_ == x.ptr()->v_; }
  bool operator!=(const Self& x) const
  { return !operator==(x); }

  operator VHandle() { return ptr()->v_; }

};

template <class Vertex, class Point>
struct Point_from_Vertex {
  typedef Vertex argument_type;
  typedef Point  result_type;
  result_type& operator()(argument_type& x) const
  { return x.point(); }
  const result_type& operator()(const argument_type& x) const
  { return x.point(); }
};


template <typename H1, typename H2>
struct list_collector {
  std::list<H1>& L_;
  list_collector(std::list<H1>& L) : L_(L) {}
  void operator()(H2 f) const
  { L_.push_back(static_cast<H1>(f)); }
};


template <class R_>
class Convex_hull_d : public Regular_complex_d<R_>
{
typedef Regular_complex_d<R_> Base;
typedef Convex_hull_d<R_>     Self;

public:

  using Base::new_simplex;
  using Base::new_vertex;
  using Base::associate_vertex_with_simplex;
  using Base::associate_point_with_vertex;
  using Base::set_neighbor;

  using Base::kernel;
  using Base::dcur;

/*{\Xgeneralization Regular_complex_d<R>}*/
/*{\Mtypes 6.5}*/
typedef R_ R;
/*{\Mtypemember the representation class.}*/
typedef typename R::Point_d Point_d;
/*{\Mtypemember the point type.}*/
typedef typename R::Hyperplane_d Hyperplane_d;
/*{\Mtypemember the hyperplane type.}*/

typedef typename R::Vector_d Vector_d;
typedef typename R::Ray_d Ray_d;
typedef typename R::RT RT;
typedef std::list<Point_d> Point_list;
// make traits types locally available

typedef typename Base::Simplex_handle Simplex_handle;
/*{\Mtypemember handle for simplices.}*/

typedef typename Base::Simplex_handle Facet_handle;
/*{\Mtypemember handle for facets.}*/

typedef typename Base::Vertex_handle Vertex_handle;
/*{\Mtypemember handle for vertices.}*/

typedef typename Base::Simplex_iterator Simplex_iterator;
/*{\Mtypemember iterator for simplices.}*/

typedef typename Base::Vertex_iterator Vertex_iterator;
/*{\Mtypemember iterator for vertices.}*/

typedef Facet_iterator_<Self*,Facet_handle> Facet_iterator;
/*{\Mtypemember iterator for facets.}*/

typedef Hull_vertex_iterator_<Self*,Vertex_handle,Facet_iterator>
  Hull_vertex_iterator;
/*{\Mtypemember iterator for vertices that are part of the convex hull.}*/

/*{\Mtext Note that each iterator fits the handle concept, i.e. iterators
can be used as handles. Note also that all iterator and handle types
come also in a const flavor, e.g., |Vertex_const_iterator| is the
constant version of |Vertex_iterator|. Thus use the const version
whenever the convex hull object is referenced as constant.}*/

#define CGAL_USING(t) typedef typename Base::t t
CGAL_USING(Simplex_const_iterator);CGAL_USING(Vertex_const_iterator);
CGAL_USING(Simplex_const_handle);CGAL_USING(Vertex_const_handle);
#undef CGAL_USING
typedef Simplex_const_handle Facet_const_handle;
typedef Facet_iterator_<const Self*,Facet_const_handle> Facet_const_iterator;
typedef Hull_vertex_iterator_<const Self*,Vertex_const_handle,
  Facet_const_iterator> Hull_vertex_const_iterator;

typedef typename Point_list::const_iterator Point_const_iterator;
/*{\Mtypemember const iterator for all inserted points.}*/

typedef typename Vertex_handle::value_type Vertex;

typedef CGAL::Iterator_project<
  Hull_vertex_const_iterator, Point_from_Vertex<Vertex,Point_d>,
  const Point_d&, const Point_d*> Hull_point_const_iterator;
/*{\Mtypemember const iterator for all points of the hull.}*/

protected:
  Point_list         all_pnts_;
  Vector_d           quasi_center_;   // sum of the vertices of origin simplex
  Simplex_handle     origin_simplex_; // pointer to the origin simplex
  Facet_handle       start_facet_;    // pointer to some facet on the surface
  Vertex_handle      anti_origin_;

public: // until there are template friend functions possible
  Point_d center() const // compute center from quasi center
  { typename R::Vector_to_point_d to_point =
      kernel().vector_to_point_d_object();
    return to_point(quasi_center_/RT(dcur + 1)); }
  const Vector_d& quasi_center() const
  { return quasi_center_; }
  Simplex_const_handle origin_simplex() const
  { return origin_simplex_; }
  Hyperplane_d base_facet_plane(Simplex_handle s) const
  { return s -> hyperplane_of_base_facet(); }
  Hyperplane_d base_facet_plane(Simplex_const_handle s) const
  { return s -> hyperplane_of_base_facet(); }
  bool& visited_mark(Simplex_handle s) const
  { return s->visited(); }

protected:
  std::size_t num_of_bounded_simplices;
  std::size_t num_of_unbounded_simplices;
  std::size_t num_of_visibility_tests;
  std::size_t num_of_vertices;

  void compute_equation_of_base_facet(Simplex_handle s);
  /*{\Xop computes the equation of the base facet of $s$ and sets the
  |base_facet| member of $s$. The equation is normalized such
  that the origin simplex lies in the negative halfspace.}*/

  bool is_base_facet_visible(Simplex_handle s, const Point_d& x) const
  { typename R::Has_on_positive_side_d has_on_positive_side =
      kernel().has_on_positive_side_d_object();
    return has_on_positive_side(s->hyperplane_of_base_facet(),x); }
  /*{\Xop returns true if $x$ sees the base facet of $s$, i.e., lies in
  its positive halfspace.}*/

  bool contains_in_base_facet(Simplex_handle s, const Point_d& x) const;
  /*{\Xop returns true if $x$ lies in the closure of the base facet of
  $s$.}*/


  void visibility_search(Simplex_handle S, const Point_d& x,
                         std::list<Simplex_handle>& visible_simplices,
                         std::size_t& num_of_visited_simplices,
                         int& location, Facet_handle& f) const;
  /*{\Xop adds all unmarked unbounded simplices with $x$-visible base
          facet to |visible_simplices| and marks them. In |location| the
          procedure returns the position of |x| with respect to the
          current hull: $-1$ for inside, $0$ for on the boundary,
          and $+1$ for outside; the initial value of |location| for the
          outermost call must be $-1$. If $x$ is contained in the
          boundary of |\Mvar| then a facet incident to $x$ is returned
          in $f$.}*/

  void clear_visited_marks(Simplex_handle s) const;
 /*{\Xop removes the mark bits from all marked simplices reachable from $s$.}*/

  void dimension_jump(Simplex_handle S, Vertex_handle x);
  /*{\Xop Adds a new vertex $x$ to the triangulation. The point associated
  with $x$ lies outside the affine hull of the current point set.  }*/

  void visible_facets_search(Simplex_handle S, const Point_d& x,
                             std::list< Facet_handle >& VisibleFacets,
                             std::size_t& num_of_visited_facets) const;

public:

  /*{\Mcreation 3}*/

  Convex_hull_d(int d, const R& Kernel = R());
  /*{\Mcreate creates an instance |\Mvar| of type |\Mtype|. The
  dimension of the underlying space is $d$ and |S| is initialized to the
  empty point set. The traits class |R| specifies the models of
  all types and the implementations of all geometric primitives used by
  the convex hull class. The default model is one of the $d$-dimensional
  representation classes (e.g., |Homogeneous_d|).}*/

  protected:
  /*{\Mtext The data type |\Mtype| offers neither copy constructor nor
  assignment operator.}*/
    Convex_hull_d(const Self&);
    Convex_hull_d& operator=(const Self&);
  public:


  /*{\Moperations 3}*/
  /*{\Mtext All operations below that take a point |x| as argument have
  the common precondition that |x| is a point of ambient space.}*/

  int dimension() const { return Base::dimension(); }
  /*{\Mop returns the dimension of ambient space}*/
  int current_dimension() const { return Base::current_dimension(); }
  /*{\Mop returns the affine dimension |dcur| of $S$.}*/

  Point_d  associated_point(Vertex_handle v) const
  { return Base::associated_point(v); }
  /*{\Mop returns the point associated with vertex $v$.}*/
  Point_d  associated_point(Vertex_const_handle v) const
  { return Base::associated_point(v); }

  Vertex_handle  vertex_of_simplex(Simplex_handle s, int i) const
  { return Base::vertex(s,i); }
  /*{\Mop returns the vertex corresponding to the $i$-th vertex of $s$.\\
  \precond $0 \leq i \leq |dcur|$. }*/
  Vertex_const_handle  vertex_of_simplex(Simplex_const_handle s, int i) const
  { return Base::vertex(s,i); }

  Point_d  point_of_simplex(Simplex_handle s,int i) const
  { return associated_point(vertex_of_simplex(s,i)); }
  /*{\Mop same as |C.associated_point(C.vertex_of_simplex(s,i))|. }*/
  Point_d  point_of_simplex(Simplex_const_handle s,int i) const
  { return associated_point(vertex_of_simplex(s,i)); }


  Simplex_handle opposite_simplex(Simplex_handle s,int i) const
  { return Base::opposite_simplex(s,i); }
  /*{\Mop returns the simplex opposite to the $i$-th vertex of $s$
  (|Simplex_handle()| if there is no such simplex).
  \precond $0 \leq i \leq |dcur|$. }*/
  Simplex_const_handle opposite_simplex(Simplex_const_handle s,int i) const
  { return Base::opposite_simplex(s,i); }


  int  index_of_vertex_in_opposite_simplex(Simplex_handle s,int i) const
  { return Base::index_of_opposite_vertex(s,i); }
  /*{\Mop returns the index of the vertex opposite to the $i$-th vertex
  of $s$. \precond $0 \leq i \leq |dcur|$ and there is a
  simplex opposite to the $i$-th vertex of $s$. }*/
  int  index_of_vertex_in_opposite_simplex(Simplex_const_handle s,int i) const
  { return Base::index_of_opposite_vertex(s,i); }

  Simplex_handle  simplex(Vertex_handle v) const
  { return Base::simplex(v); }
  /*{\Mop returns a simplex of which $v$ is a node. Note that this
  simplex is not unique. }*/
  Simplex_const_handle  simplex(Vertex_const_handle v) const
  { return Base::simplex(v); }

  int index(Vertex_handle v) const { return Base::index(v); }
  /*{\Mop returns the index of $v$ in |simplex(v)|.}*/
  int index(Vertex_const_handle v) const { return Base::index(v); }

  Vertex_handle  vertex_of_facet(Facet_handle f, int i) const
  { return vertex_of_simplex(f,i+1); }
  /*{\Mop returns the vertex corresponding to the $i$-th vertex of $f$.
  \precond $0 \leq i < |dcur|$. }*/
  Vertex_const_handle  vertex_of_facet(Facet_const_handle f, int i) const
  { return vertex_of_simplex(f,i+1); }

  Point_d point_of_facet(Facet_handle f, int i) const
  { return point_of_simplex(f,i+1); }
  /*{\Mop same as |C.associated_point(C.vertex_of_facet(f,i))|.}*/
  Point_d point_of_facet(Facet_const_handle f, int i) const
  { return point_of_simplex(f,i+1); }

  Facet_handle opposite_facet(Facet_handle f, int i) const
  { return opposite_simplex(f,i+1); }
  /*{\Mop returns the facet opposite to the $i$-th vertex of $f$
  (|Facet_handle()| if there is no such facet). \precond $0 \leq i <
  |dcur|$ and |dcur > 1|. }*/
  Facet_const_handle opposite_facet(Facet_const_handle f, int i) const
  { return opposite_simplex(f,i+1); }

  int index_of_vertex_in_opposite_facet(Facet_handle f, int i) const
  { return index_of_vertex_in_opposite_simplex(f,i+1) - 1; }
  /*{\Mop returns the index of the vertex opposite to the $i$-th vertex of
  $f$. \precond $0 \leq i < |dcur|$ and |dcur > 1|.}*/
  int index_of_vertex_in_opposite_facet(Facet_const_handle f, int i) const
  { return index_of_vertex_in_opposite_simplex(f,i+1) - 1; }

  Hyperplane_d hyperplane_supporting(Facet_handle f) const
  { return f -> hyperplane_of_base_facet(); }
  /*{\Mop returns a hyperplane supporting facet |f|. The hyperplane is
  oriented such that the interior of |\Mvar| is on the negative
  side of it. \precond |f| is a facet of |\Mvar| and |dcur > 1|.}*/
  Hyperplane_d hyperplane_supporting(Facet_const_handle f) const
  { return f -> hyperplane_of_base_facet(); }


  Vertex_handle insert(const Point_d& x);
  /*{\Mop adds point |x| to the underlying set of points.  If $x$ is
  equal to (the point associated with) a vertex of the current hull this
  vertex is returned and its associated point is changed to $x$. If $x$
  lies outside the current hull, a new vertex |v| with associated point
  $x$ is added to the hull and returned. In all other cases, i.e., if
  $x$ lies in the interior of the hull or on the boundary but not on a
  vertex, the current hull is not changed and |Vertex_handle()| is
  returned. If |CGAL_CHECK_EXPENSIVE| is defined then the validity
  check |is_valid(true)| is executed as a post condition.}*/

  template <typename Forward_iterator>
  void insert(Forward_iterator first, Forward_iterator last)
  { while (first != last) insert(*first++); }
  /*{\Mop adds |S = set [first,last)| to the underlying set of
  points. If any point |S[i]| is equal to (the point associated with) a
  vertex of the current hull its associated point is changed to |S[i]|.}*/


  bool is_dimension_jump(const Point_d& x) const
  /*{\Mop returns true if $x$ is not contained in the affine hull of |S|.}*/
  {
    if (current_dimension() == dimension()) return false;
    typename R::Contained_in_affine_hull_d contained_in_affine_hull =
      kernel().contained_in_affine_hull_d_object();
    return ( !contained_in_affine_hull(origin_simplex_->points_begin(),
      origin_simplex_->points_begin()+current_dimension()+1,x) );
  }


  std::list<Facet_handle>  facets_visible_from(const Point_d& x);
  /*{\Mop returns the list of all facets that are visible from |x|.\\
  \precond |x| is contained in the affine hull of |S|.}*/

  Bounded_side bounded_side(const Point_d& x);
  /*{\Mop returns |ON_BOUNDED_SIDE| (|ON_BOUNDARY|,|ON_UNBOUNDED_SIDE|)
  if |x| is contained in the interior (lies on the boundary, is contained
  in the exterior) of |\Mvar|. \precond |x| is contained in the affine
  hull of |S|.}*/


  bool is_unbounded_simplex(Simplex_handle S) const
  { return vertex_of_simplex(S,0) == anti_origin_; }
  bool is_unbounded_simplex(Simplex_const_handle S) const
  { return vertex_of_simplex(S,0) == anti_origin_; }
  bool is_bounded_simplex(Simplex_handle S) const
  { return vertex_of_simplex(S,0) != anti_origin_; }
  bool is_bounded_simplex(Simplex_const_handle S) const
  { return vertex_of_simplex(S,0) != anti_origin_; }

  void clear(int d)
  /*{\Mop reinitializes |\Mvar| to an empty hull in $d$-dimensional space.}*/
  {
    typename R::Construct_vector_d create =
      kernel().construct_vector_d_object();
    quasi_center_ = create(d,NULL_VECTOR);
    anti_origin_ = Vertex_handle();
    origin_simplex_ = Simplex_handle();
    all_pnts_.clear();
    Base::clear(d);
    num_of_vertices = 0;
    num_of_unbounded_simplices = num_of_bounded_simplices = 0;
    num_of_visibility_tests = 0;
  }

  std::size_t number_of_vertices() const
  /*{\Mop returns the number of vertices of |\Mvar|.}*/
  { return num_of_vertices; }

  std::size_t number_of_facets() const
  /*{\Mop returns the number of facets of |\Mvar|.}*/
  { return num_of_unbounded_simplices; }

  std::size_t number_of_simplices() const
  /*{\Mop returns the number of bounded simplices of |\Mvar|.}*/
  { return num_of_bounded_simplices; }

  void print_statistics()
  /*{\Mop gives information about the size of the current hull and the
  number of visibility tests performed.}*/
  {
    std::cout << "Convex_hull_d ("
              << current_dimension() << "/" << dimension()
              << ") - statistic" << std::endl;
    std::cout<<" # points = " << all_pnts_.size() << std::endl;
    std::cout<<" # vertices = " << num_of_vertices << std::endl;
    std::cout<<" # unbounded simplices = " << num_of_unbounded_simplices
      << std::endl;
    std::cout<<" # bounded simplices = " << num_of_bounded_simplices
      << std::endl;
    std::cout<<" # visibility tests = " << num_of_visibility_tests
      << std::endl;
  }

  class chull_has_double_coverage {};
  class chull_has_local_non_convexity {};
  class chull_has_center_on_wrong_side_of_hull_facet {};

  bool is_valid(bool throw_exceptions = false) const;
  /*{\Mop checks the validity of the data structure.
  If |throw_exceptions == thrue| then the program throws
  the following exceptions to inform about the problem.\\
  [[chull_has_center_on_wrong_side_of_hull_facet]] the hyperplane
  supporting a facet has the wrong orientation.\\
  [[chull_has_local_non_convexity]] a ridge is locally non convex.\\
  [[chull_has_double_coverage]] the hull has a winding number larger
  than 1.
  }*/

  template <typename Forward_iterator>
  void initialize(Forward_iterator first, Forward_iterator last)
  /*{\Xop initializes the complex with the set |S = set [first,last)|.
  The initialization uses the quickhull approach.}*/
  { CGAL_assertion(current_dimension()==-1);
    Vertex_handle z;
    Forward_iterator it(first);
    std::list< Point_d > OtherPoints;
    typename std::list< Point_d >::iterator pit, pred;
    while ( it != last ) {
      Point_d x = *it++;
      if ( current_dimension() == -1 ) {

        Simplex_handle outer_simplex; // a pointer to the outer simplex
        dcur = 0;  // we jump from dimension - 1 to dimension 0
        origin_simplex_ = new_simplex();  num_of_bounded_simplices ++;
        outer_simplex  = new_simplex();   num_of_unbounded_simplices ++;
        start_facet_ = origin_simplex_;
        z = new_vertex(x);                num_of_vertices ++;
        associate_vertex_with_simplex(origin_simplex_,0,z);
          // z is the only point and the peak
        associate_vertex_with_simplex(outer_simplex,0,anti_origin_);
        set_neighbor(origin_simplex_,0,outer_simplex,0);
        typename R::Point_to_vector_d to_vector =
          kernel().point_to_vector_d_object();
        quasi_center_ = to_vector(x);


      }
      else if ( is_dimension_jump(x) ) {
        dcur++;
        z = new_vertex(x); num_of_vertices++;
        typename R::Point_to_vector_d to_vector =
          kernel().point_to_vector_d_object();
        quasi_center_ = quasi_center_ + to_vector(x);
        dimension_jump(origin_simplex_, z);
        clear_visited_marks(origin_simplex_);
        Simplex_iterator S;
        forall_rc_simplices(S,*this) compute_equation_of_base_facet(S);
        num_of_unbounded_simplices += num_of_bounded_simplices;
        if (dcur > 1) {
          start_facet_ = opposite_simplex(origin_simplex_,dcur);
          CGAL_assertion(vertex_of_simplex(start_facet_,0)==Vertex_handle());
        }


      } else {
        OtherPoints.push_back(x);
      }
    }
    all_pnts_.insert(all_pnts_.end(),first,last);

    // what is the point of this line? ...
    //int dcur = current_dimension();
    Unique_hash_map<Facet_handle, std::list<Point_d> > PointsOf;
    std::list<Facet_handle> FacetCandidates;
    typename R::Oriented_side_d side_of =
      kernel().oriented_side_d_object();
    for (int i=0; i<=dcur; ++i) {
      Simplex_handle f = opposite_simplex(origin_simplex_,i);
      Hyperplane_d h = f->hyperplane_of_base_facet();
      std::list<Point_d>& L = PointsOf[f];
      pit = OtherPoints.begin();
      while ( pit != OtherPoints.end() ) {
        if ( side_of(h,*pit) == ON_POSITIVE_SIDE ) {
          L.push_back(*pit); pred=pit; ++pit; OtherPoints.erase(pred);
        } else ++pit; // not above h
      }
      if ( !L.empty() ) FacetCandidates.push_back(f);
    }
    OtherPoints.clear();

    while ( !FacetCandidates.empty() ) {
      Facet_handle f = *FacetCandidates.begin();
      FacetCandidates.pop_front();
      Hyperplane_d h = f->hyperplane_of_base_facet();
      std::list<Point_d>& L = PointsOf[f];
      if (L.empty()) { CGAL_assertion( is_bounded_simplex(f) ); continue; }
      Point_d p = *L.begin();
      typename R::Value_at_d value_at = kernel().value_at_d_object();
      RT maxdist = value_at(h,p), dist;
      for (pit = ++L.begin(); pit != L.end(); ++pit) {
        dist = value_at(h,*pit);
        if ( dist > maxdist ) { maxdist = dist; p = *pit; }
      }
      num_of_visibility_tests += L.size();

      std::size_t num_of_visited_facets = 0;
      std::list<Facet_handle> VisibleFacets;
      VisibleFacets.push_back(f);
      visible_facets_search(f, p, VisibleFacets, num_of_visited_facets);
      num_of_visibility_tests += num_of_visited_facets;
      num_of_bounded_simplices += VisibleFacets.size();
      clear_visited_marks(f);

      ++num_of_vertices;
      Vertex_handle z = new_vertex(p);
      std::list<Simplex_handle> NewSimplices;
      typename std::list<Facet_handle>::iterator it;
      CGAL_assertion(OtherPoints.empty());
      for (it = VisibleFacets.begin(); it != VisibleFacets.end(); ++it) {
        OtherPoints.splice(OtherPoints.end(),PointsOf[*it]);
        Facet_handle S = *it;
        associate_vertex_with_simplex(S,0,z);
        for (int k = 1; k <= dcur; ++k) {
          if (!is_base_facet_visible(opposite_simplex(S,k),p)) {
            Simplex_handle T = new_simplex();
            NewSimplices.push_back(T);

            /* set the vertices of T as described above */
            for (int i = 1; i < dcur; i++) {
              if ( i != k )
                associate_vertex_with_simplex(T,i,vertex_of_simplex(S,i));
            }
            if (k != dcur)
              associate_vertex_with_simplex(T,k,vertex_of_simplex(S,dcur));
            associate_vertex_with_simplex(T,dcur,z);
            associate_vertex_with_simplex(T,0,anti_origin_);
            /* in the above, it is tempting to drop the tests ( i != k ) and
               ( k != dcur ) since the subsequent lines after will correct the
               erroneous assignment.  This reasoning is fallacious as the
               procedure assoc_vertex_with_simplex also the internal data of
               the third argument. */

            /* compute the equation of its base facet */
            compute_equation_of_base_facet(T);

            /* record adjacency information for the two known neighbors */
            set_neighbor(T,dcur,opposite_simplex(S,k),
                         index_of_vertex_in_opposite_simplex(S,k));
            set_neighbor(T,0,S,k);


          }
        }
      }
      num_of_unbounded_simplices -= VisibleFacets.size();
      if ( vertex_of_simplex(start_facet_,0) != Vertex_handle() )
        start_facet_ = *(--NewSimplices.end());
      CGAL_assertion( vertex_of_simplex(start_facet_,0)==Vertex_handle() );
      for (it = NewSimplices.begin(); it != NewSimplices.end(); ++it) {
        Simplex_handle Af = *it;
        for (int k = 1; k < dcur ; k++) {
          // neighbors 0 and dcur are already known
          if (opposite_simplex(Af,k) == Simplex_handle()) {
            // we have not performed the walk in the opposite direction yet
            Simplex_handle T = opposite_simplex(Af,0);
            int y1 = 0;
            while ( vertex_of_simplex(T,y1) != vertex_of_simplex(Af,k) )
              y1++;
            // exercise: show that we can also start with y1 = 1
            int y2 = index_of_vertex_in_opposite_simplex(Af,0);

            while ( vertex_of_simplex(T,0) == z ) {
              // while T has peak x do find new y_1 */
              int new_y1 = 0;
              while (vertex_of_simplex(opposite_simplex(T,y1),new_y1) !=
                     vertex_of_simplex(T,y2))
                new_y1++;
                // exercise: show that we can also start with new_y1 = 1
              y2 = index_of_vertex_in_opposite_simplex(T,y1);
              T =  opposite_simplex(T,y1);
              y1 = new_y1;
            }
            set_neighbor(Af,k,T,y1); // update adjacency information
          }
        }
      }


      for (it = NewSimplices.begin(); it != NewSimplices.end(); ++it) {
        Facet_handle f = *it;
        CGAL_assertion( is_unbounded_simplex(f) );
        Hyperplane_d h = f->hyperplane_of_base_facet();
        std::list<Point_d>& L = PointsOf[f];
        pit = OtherPoints.begin();
        while ( pit != OtherPoints.end() ) {
          if ( side_of(h,*pit) == ON_POSITIVE_SIDE ) {
            L.push_back(*pit); pred=pit; ++pit; OtherPoints.erase(pred);
          } else ++pit; // not above h
        }
        if ( !L.empty() ) FacetCandidates.push_back(f);
      }
      OtherPoints.clear();

    }

  #ifdef CGAL_CHECK_EXPENSIVE
    CGAL_assertion(is_valid(true));
  #endif
  }


  /*{\Mtext \headerline{Lists and Iterators}
  \setopdims{3.5cm}{3.5cm}}*/

  Vertex_iterator vertices_begin() { return Base::vertices_begin(); }
  /*{\Mop returns an iterator to start iteration over all vertices
  of |\Mvar|.}*/
  Vertex_iterator vertices_end()   { return Base::vertices_end(); }
  /*{\Mop the past the end iterator for vertices.}*/
  Simplex_iterator simplices_begin() { return Base::simplices_begin(); }
  /*{\Mop returns an iterator to start iteration over all simplices
  of |\Mvar|.}*/
  Simplex_iterator simplices_end()   { return Base::simplices_end(); }
  /*{\Mop the past the end iterator for simplices.}*/
  Facet_iterator facets_begin() { return Facet_iterator(this,start_facet_); }
  /*{\Mop returns an iterator to start iteration over all facets of |\Mvar|.}*/
  Facet_iterator facets_end() { return Facet_iterator(); }
  /*{\Mop the past the end iterator for facets.}*/
  Hull_vertex_iterator hull_vertices_begin()
  { return Hull_vertex_iterator(this,facets_begin()); }
  /*{\Mop returns an iterator to start iteration over all hull vertex
  of |\Mvar|. Remember that the hull is a simplicial complex.}*/
  Hull_vertex_iterator hull_vertices_end()
  { return Hull_vertex_iterator(); }
  /*{\Mop the past the end iterator for hull vertices.}*/

  Point_const_iterator points_begin() const { return all_pnts_.begin(); }
  /*{\Mop returns the start iterator for all points that have been
  inserted to construct |\Mvar|.}*/
  Point_const_iterator points_end() const { return all_pnts_.end(); }
  /*{\Mop returns the past the end iterator for all points.}*/

  Hull_point_const_iterator hull_points_begin() const
  { return Hull_point_const_iterator(hull_vertices_begin()); }
  /*{\Mop returns an iterator to start iteration over all inserted
  points that are part of the convex hull |\Mvar|. Remember that the
  hull is a simplicial complex.}*/
  Hull_point_const_iterator hull_points_end() const
  { return Hull_point_const_iterator(hull_vertices_end()); }
  /*{\Mop returns the past the end iterator for points of the hull.}*/

  Vertex_const_iterator vertices_begin() const
  { return Base::vertices_begin(); }
  Vertex_const_iterator vertices_end()   const
  { return Base::vertices_end(); }
  Simplex_const_iterator simplices_begin() const
  { return Base::simplices_begin(); }
  Simplex_const_iterator simplices_end()   const
  { return Base::simplices_end(); }
  Facet_const_iterator facets_begin() const
  { return Facet_const_iterator(this,start_facet_); }
  Facet_const_iterator facets_end() const
  { return Facet_const_iterator(); }
  Hull_vertex_const_iterator hull_vertices_begin() const
  { return Hull_vertex_const_iterator(this,facets_begin()); }
  Hull_vertex_const_iterator hull_vertices_end() const
  { return Hull_vertex_const_iterator(); }

  /* We distinguish cases according to the current dimension.  If the
     dimension is less than one then the hull has no facets, if the
     dimension is one then the hull has two facets which we extract by a
     scan through the set of all simplices, and if the hull has
     dimension at least two the boundary is a connected set and we
     construct the list of facets by depth first search starting at
     |start_facet_|.*/

  /*{\Mtext\setopdims{5.5cm}{3.5cm}}*/
  template <typename Visitor>
  void visit_all_facets(const Visitor& V) const
  /*{\Mop each facet of |\Mvar| is visited by the visitor object |V|.
  |V| has to have a function call operator:\\
  |void operator()(Facet_handle) const|}*/
  {
    if (current_dimension() > 1) {
      Unique_hash_map<Facet_handle,bool> visited(false);
      std::list<Facet_handle> candidates;
      candidates.push_back(start_facet_);
      visited[start_facet_] = true;
      Facet_handle current;
      while ( !candidates.empty() ) {
        current = *(candidates.begin()); candidates.pop_front();
        CGAL_assertion(vertex_of_simplex(current,0)==Vertex_handle());
        V(current);
        for(int i = 1; i <= dcur; ++i) {
          Facet_handle f = opposite_simplex(current,i);
          if ( !visited[f] ) {
            candidates.push_back(f);
            visited[f] = true;
          }
        }
      }
    }
    else if ( current_dimension() == 1 ) { V(start_facet_); }
  }

  const std::list<Point_d>& all_points() const
  /*{\Mop returns a list of all points in |\Mvar|.}*/
  { return all_pnts_; }

  std::list<Vertex_handle> all_vertices()
  /*{\Mop returns a list of all vertices of |\Mvar|
  (also interior ones).}*/
  { return Base::all_vertices(); }
  std::list<Vertex_const_handle> all_vertices() const
  { return Base::all_vertices(); }

  std::list<Simplex_handle> all_simplices()
  /*{\Mop returns a list of all simplices in |\Mvar|.}*/
  { return Base::all_simplices(); }
  std::list<Simplex_const_handle> all_simplices() const
  { return Base::all_simplices(); }


  std::list<Facet_handle>  all_facets()
  /*{\Mop returns a list of all facets of |\Mvar|.}*/
  { std::list<Facet_handle> L;
    list_collector<Facet_handle,Facet_handle> visitor(L);
    visit_all_facets(visitor);
    return L;
  }

  std::list<Facet_const_handle>  all_facets() const
  { std::list<Facet_const_handle> L;
    list_collector<Facet_const_handle,Facet_handle> visitor(L);
    visit_all_facets(visitor);
    return L;
  }

  #define forall_ch_vertices(x,CH)\
  for(x = (CH).vertices_begin(); x != (CH).vertices_end(); ++x)
  #define forall_ch_simplices(x,CH)\
  for(x = (CH).simplices_begin(); x != (CH).simplices_end(); ++x)
  #define forall_ch_facets(x,CH)\
  for(x = (CH).facets_begin(); x != (CH).facets_end(); ++x)

  /*{\Mtext
  \headerline{Iteration Statements}

  {\bf forall\_ch\_vertices}($v,C$)
  $\{$ ``the vertices of $C$ are successively assigned to $v$'' $\}$

  {\bf forall\_ch\_simplices}($s,C$)
  $\{$ ``the simplices of $C$ are successively assigned to $s$'' $\}$

  {\bf forall\_ch\_facets}($f,C$)
  $\{$ ``the facets of $C$ are successively assigned to $f$'' $\}$

  }*/


  /*{\Mimplementation The implementation of type |\Mtype| is based on
  \cgalCite{cms:fourresults} and \cgalCite{BMS:degeneracy}.  The details of the
  implementation can be found in the implementation document available
  at the download site of this package.

  The time and space requirements are input dependent.  Let $C_1$, $C_2$,
  $C_3$, \ldots be the sequence of hulls constructed and for a point $x$
  let $k_i$ be the number of facets of $C_i$ that are visible from $x$
  and that are not already facets of $C_{i-1}$. Then the time for
  inserting $x$ is $O(|dim| \sum_i k_i)$ and the number of new simplices
  constructed during the insertion of $x$ is the number of facets of the
  hull which were not already facets of the hull before the insertion.

  The data type |\Mtype| is derived from |Regular_complex_d|. The space
  requirement of regular complexes is essentially $12(|dim| +2)$ bytes
  times the number of simplices plus the space for the points. |\Mtype|
  needs an additional $8 + (4 + x)|dim|$ bytes per simplex where $x$ is
  the space requirement of the underlying number type and an additional
  $12$ bytes per point. The total is therefore $(16 + x)|dim| + 32$
  bytes times the number of simplices plus $28 + x \cdot |dim|$ bytes
  times the number of points.}*/

  /*{\Mtext\headerline{Traits requirements}

  |\Mname| requires the following types from the kernel traits |R|:
  \begin{Mverb}
    Point_d Vector_d Ray_d Hyperplane_d
  \end{Mverb}
  and uses the following function objects from the kernel traits:
  \begin{Mverb}
    Construct_vector_d
    Construct_hyperplane_d
    Vector_to_point_d / Point_to_vector_d
    Orientation_d
    Orthogonal_vector_d
    Oriented_side_d / Has_on_positive_side_d
    Affinely_independent_d
    Contained_in_simplex_d
    Contained_in_affine_hull_d
    Intersect_d
  \end{Mverb}
  }*/

}; // Convex_hull_d<R>


template <class R>
Convex_hull_d<R>::Convex_hull_d(int d, const R& Kernel) : Base(d,Kernel)
{
  origin_simplex_ = Simplex_handle();
  start_facet_ = Facet_handle();
  anti_origin_ = Vertex_handle();
  num_of_vertices = 0;
  num_of_unbounded_simplices = num_of_bounded_simplices = 0;
  num_of_visibility_tests = 0;
  typename R::Construct_vector_d create =
    kernel().construct_vector_d_object();
  quasi_center_ = create(d,NULL_VECTOR);
}

template <class R>
bool Convex_hull_d<R>::
contains_in_base_facet(Simplex_handle s, const Point_d& x) const
{
  typename R::Contained_in_simplex_d contained_in_simplex =
    kernel().contained_in_simplex_d_object();
  return contained_in_simplex(s->points_begin()+1,
    s->points_begin()+current_dimension()+1,x);
}

template <class R>
void Convex_hull_d<R>::
compute_equation_of_base_facet(Simplex_handle S)
{
  typename R::Construct_hyperplane_d hyperplane_through_points =
    kernel().construct_hyperplane_d_object();
  S->set_hyperplane_of_base_facet( hyperplane_through_points(
    S->points_begin()+1, S->points_begin()+1+current_dimension(),
    center(), ON_NEGATIVE_SIDE)); // skip the first point !

  #ifdef CGAL_CHECK_EXPENSIVE
  { /* Let us check */
    typename R::Oriented_side_d side = kernel().oriented_side_d_object();
    for (int i = 1; i <= current_dimension(); i++)
      CGAL_assertion_msg(side(S->hyperplane_of_base_facet(),
                              point_of_simplex(S,i)) == ON_ORIENTED_BOUNDARY,
        " hyperplane does not support base ");
    CGAL_assertion_msg(side(S->hyperplane_of_base_facet(),center()) ==
                       ON_NEGATIVE_SIDE,
      " hyperplane has quasi center on wrong side ");
  }
  #endif


}

template <class R>
typename Convex_hull_d<R>::Vertex_handle
Convex_hull_d<R>::insert(const Point_d& x)
{
  Vertex_handle z = Vertex_handle();
  all_pnts_.push_back(x);
  if (current_dimension() == -1) {

    Simplex_handle outer_simplex; // a pointer to the outer simplex
    dcur = 0;  // we jump from dimension - 1 to dimension 0
    origin_simplex_ = new_simplex();  num_of_bounded_simplices ++;
    outer_simplex  = new_simplex();   num_of_unbounded_simplices ++;
    start_facet_ = origin_simplex_;
    z = new_vertex(x);                num_of_vertices ++;
    associate_vertex_with_simplex(origin_simplex_,0,z);
      // z is the only point and the peak
    associate_vertex_with_simplex(outer_simplex,0,anti_origin_);
    set_neighbor(origin_simplex_,0,outer_simplex,0);
    typename R::Point_to_vector_d to_vector =
      kernel().point_to_vector_d_object();
    quasi_center_ = to_vector(x);


  } else if ( is_dimension_jump(x) ) {
    dcur++;
    z = new_vertex(x); num_of_vertices++;
    typename R::Point_to_vector_d to_vector =
      kernel().point_to_vector_d_object();
    quasi_center_ = quasi_center_ + to_vector(x);
    dimension_jump(origin_simplex_, z);
    clear_visited_marks(origin_simplex_);
    Simplex_iterator S;
    forall_rc_simplices(S,*this) compute_equation_of_base_facet(S);
    num_of_unbounded_simplices += num_of_bounded_simplices;
    if (dcur > 1) {
      start_facet_ = opposite_simplex(origin_simplex_,dcur);
      CGAL_assertion(vertex_of_simplex(start_facet_,0)==Vertex_handle());
    }


  } else {
    if ( current_dimension() == 0 ) {
      z = vertex_of_simplex(origin_simplex_,0);
      associate_point_with_vertex(z,x);
      return z;
    }
    std::list<Simplex_handle> visible_simplices;
    int location = -1;
    Facet_handle f;

    std::size_t num_of_visited_simplices = 0;

    visibility_search(origin_simplex_, x, visible_simplices,
                      num_of_visited_simplices, location, f);

    num_of_visibility_tests += num_of_visited_simplices;

    #ifdef COUNTS
      cout << "\nthe number of visited simplices in this iteration is ";
      cout << num_of_visited_simplices << endl;
    #endif

    clear_visited_marks(origin_simplex_);


    #ifdef COUNTS
      cout << "\nthe number of bounded simplices constructed ";
      cout << " in this iteration is  " << visible_simplices.size() << endl;
    #endif

    num_of_bounded_simplices += visible_simplices.size();

    switch (location) {
      case -1:
        return Vertex_handle();
      case 0:
        { for (int i = 0; i < current_dimension(); i++) {
            if ( x == point_of_facet(f,i) ) {
              z = vertex_of_facet(f,i);
              associate_point_with_vertex(z,x);
              return z;
            }
          }
          return Vertex_handle();
        }
      case 1:
        { num_of_vertices++;
          z = new_vertex(x);
          std::list<Simplex_handle> NewSimplices; // list of new simplices
          typename std::list<Simplex_handle>::iterator it;

          for (it = visible_simplices.begin();
               it != visible_simplices.end(); ++it) {
            Simplex_handle S = *it;
            associate_vertex_with_simplex(S,0,z);
            for (int k = 1; k <= dcur; k++) {
              if (!is_base_facet_visible(opposite_simplex(S,k),x)) {
                Simplex_handle T = new_simplex();
                NewSimplices.push_back(T);

                /* set the vertices of T as described above */
                for (int i = 1; i < dcur; i++) {
                  if ( i != k )
                    associate_vertex_with_simplex(T,i,vertex_of_simplex(S,i));
                }
                if (k != dcur)
                  associate_vertex_with_simplex(T,k,vertex_of_simplex(S,dcur));
                associate_vertex_with_simplex(T,dcur,z);
                associate_vertex_with_simplex(T,0,anti_origin_);
                /* in the above, it is tempting to drop the tests ( i != k )
                   and ( k != dcur ) since the subsequent lines after will
                   correct the erroneous assignment.  This reasoning is
                   fallacious as the procedure assoc_vertex_with_simplex also
                   the internal data of the third argument. */

                /* compute the equation of its base facet */
                compute_equation_of_base_facet(T);

                /* record adjacency information for the two known neighbors */
                set_neighbor(T,dcur,opposite_simplex(S,k),
                             index_of_vertex_in_opposite_simplex(S,k));
                set_neighbor(T,0,S,k);


              }
            }
          }
          num_of_unbounded_simplices -= visible_simplices.size();
          if ( vertex_of_simplex(start_facet_,0) != Vertex_handle() )
            start_facet_ = *(--NewSimplices.end());
          CGAL_assertion( vertex_of_simplex(start_facet_,0)==Vertex_handle() );
          for (it = NewSimplices.begin(); it != NewSimplices.end(); ++it) {
            Simplex_handle Af = *it;
            for (int k = 1; k < dcur ; k++) {
              // neighbors 0 and dcur are already known
              if (opposite_simplex(Af,k) == Simplex_handle()) {
                // we have not performed the walk in the opposite direction yet
                Simplex_handle T = opposite_simplex(Af,0);
                int y1 = 0;
                while ( vertex_of_simplex(T,y1) != vertex_of_simplex(Af,k) )
                  y1++;
                // exercise: show that we can also start with y1 = 1
                int y2 = index_of_vertex_in_opposite_simplex(Af,0);

                while ( vertex_of_simplex(T,0) == z ) {
                  // while T has peak x do find new y_1 */
                  int new_y1 = 0;
                  while (vertex_of_simplex(opposite_simplex(T,y1),new_y1) !=
                         vertex_of_simplex(T,y2))
                    new_y1++;
                    // exercise: show that we can also start with new_y1 = 1
                  y2 = index_of_vertex_in_opposite_simplex(T,y1);
                  T =  opposite_simplex(T,y1);
                  y1 = new_y1;
                }
                set_neighbor(Af,k,T,y1); // update adjacency information
              }
            }
          }


        }
    }

  }
#ifdef CGAL_CHECK_EXPENSIVE
  CGAL_assertion(is_valid(true));
#endif
  return z;
}


template <class R>
void Convex_hull_d<R>::
visibility_search(Simplex_handle S, const Point_d& x,
                  std::list< Simplex_handle >& visible_simplices,
                  std::size_t& num_of_visited_simplices, int& location,
                  Simplex_handle& f) const
{
  num_of_visited_simplices ++;
  S->visited() = true; // we have visited S and never come back ...
  for(int i = 0; i <= current_dimension(); ++i) {
    Simplex_handle T = opposite_simplex(S,i); // for all neighbors T of S
    if ( !T->visited() ) {
      typename R::Oriented_side_d side_of =
        kernel().oriented_side_d_object();
      int side = side_of(T->hyperplane_of_base_facet(),x);
      if ( is_unbounded_simplex(T) ) {
        if ( side == ON_POSITIVE_SIDE ) {
          // T is an unbounded simplex with x-visible base facet
          visible_simplices.push_back(T);
          location = 1;
        }
        if ( side == ON_ORIENTED_BOUNDARY &&
             location == -1 && contains_in_base_facet(T,x) ) {
          location = 0; f = T;
          return;
        }
      }
      if ( side == ON_POSITIVE_SIDE ||
           (side == ON_ORIENTED_BOUNDARY && location == -1) ) {
        visibility_search(T,x,visible_simplices,
                          num_of_visited_simplices,location,f);
        // do the recursive search
      }
    } // end visited
    else {
    }
  } // end for
}

template <class R>
void Convex_hull_d<R>::clear_visited_marks(Simplex_handle S) const
{
  S->visited() = false; // clear the visit - bit
  for(int i = 0; i <= current_dimension(); i++) // for all neighbors of S
    if (opposite_simplex(S,i)->visited())
      // if the i - th neighbor has been visited
      clear_visited_marks(opposite_simplex(S,i));
      // clear its bit recursively
}

template <class R>
std::list< typename Convex_hull_d<R>::Simplex_handle >
Convex_hull_d<R>::facets_visible_from(const Point_d& x)
{
  std::list<Simplex_handle> visible_simplices;
  int location = -1;                       // intialization is important
  std::size_t num_of_visited_simplices = 0;     // irrelevant
  Facet_handle f;                          // irrelevant

  visibility_search(origin_simplex_, x, visible_simplices,
                    num_of_visited_simplices, location, f);
  clear_visited_marks(origin_simplex_);
  return visible_simplices;
}

template <class R>
Bounded_side Convex_hull_d<R>::bounded_side(const Point_d& x)
{
  if ( is_dimension_jump(x) ) return ON_UNBOUNDED_SIDE;
  std::list<Simplex_handle> visible_simplices;
  int location = -1;                       // intialization is important
  std::size_t num_of_visited_simplices = 0;     // irrelevant
  Facet_handle f;

  visibility_search(origin_simplex_, x, visible_simplices,
                    num_of_visited_simplices, location, f);
  clear_visited_marks(origin_simplex_);
  switch (location) {
    case -1: return ON_BOUNDED_SIDE;
    case  0: return ON_BOUNDARY;
    case  1: return ON_UNBOUNDED_SIDE;
  }
  CGAL_error(); return ON_BOUNDARY; // never come here
}


template <class R>
void Convex_hull_d<R>::
dimension_jump(Simplex_handle S, Vertex_handle x)
{
  Simplex_handle S_new;
  S->visited() = true;
  associate_vertex_with_simplex(S,dcur,x);
  if ( !is_unbounded_simplex(S) ) { // S is bounded
    S_new = new_simplex();
    set_neighbor(S,dcur,S_new,0);
    associate_vertex_with_simplex(S_new,0,anti_origin_);
    for (int k = 1; k <= dcur; k++) {
      associate_vertex_with_simplex(S_new,k,vertex_of_simplex(S,k-1));
    }

  }
  /* visit unvisited neighbors 0 to dcur - 1 */
  for (int k = 0; k <= dcur - 1; k++) {
    if ( !opposite_simplex(S,k) -> visited() ) {
      dimension_jump(opposite_simplex(S,k), x);
    }
  }
  /* the recursive calls ensure that all neighbors exist */
  if ( is_unbounded_simplex(S) ) {
    set_neighbor(S,dcur,opposite_simplex(opposite_simplex(S,0),dcur),
                 index_of_vertex_in_opposite_simplex(S,0) + 1);

  } else {
    for (int k = 0; k < dcur; k++) {
      if ( is_unbounded_simplex(opposite_simplex(S,k)) ) {
        // if F' is unbounded
        set_neighbor(S_new,k + 1,opposite_simplex(S,k),dcur);
        // the neighbor of S_new opposite to v is S' and
        // x is in position dcur
      } else { // F' is bounded
        set_neighbor(S_new,k + 1,opposite_simplex(opposite_simplex(S,k),dcur),
                     index_of_vertex_in_opposite_simplex(S,k) + 1);
        // neighbor of S_new opposite to v is S_new'
        // the vertex opposite to v remains the same but ...
        // remember the shifting of the vertices one step to the right
      }
    }

  }
}

template <class R>
bool Convex_hull_d<R>::is_valid(bool throw_exceptions) const
{
  this->check_topology();
  if (current_dimension() < 1) return true;

  /* Recall that center() gives us the center-point of the origin
     simplex. We check whether it is locally inside with respect to
     all hull facets.  */

  typename R::Oriented_side_d side =
    kernel().oriented_side_d_object();
  Point_d centerpoint = center();
  Simplex_const_iterator S;
  forall_rc_simplices(S,*this) {
    if ( is_unbounded_simplex(S) &&
         side(S->hyperplane_of_base_facet(),centerpoint) !=
         ON_NEGATIVE_SIDE) {
      if (throw_exceptions)
        throw chull_has_center_on_wrong_side_of_hull_facet();
      return false;
    }
  }
  /* next we check convexity at every ridge. Let |S| be any hull
     simplex and let |v| be any vertex of its base facet. The vertex
     opposite to |v| must not be on the positive side of the base
     facet.*/

  forall_rc_simplices(S,*this) {
    if ( is_unbounded_simplex(S) ) {
      for (int i = 1; i <= dcur; i++) {
        int k = index_of_vertex_in_opposite_simplex(S,i);
        if (side(S->hyperplane_of_base_facet(),
            point_of_simplex(opposite_simplex(S,i),k)) ==
            ON_POSITIVE_SIDE) {
          if (throw_exceptions)
            throw chull_has_local_non_convexity();
          return false;
        }
      }
    }
  }

  /* next we select one hull facet */

  Simplex_const_handle selected_hull_simplex;
  forall_rc_simplices(S,*this) {
    if ( is_unbounded_simplex(S) ) { selected_hull_simplex = S; break; }
  }

  /* we compute the center of gravity of the base facet of the
     hull simplex */

  typename R::Point_to_vector_d to_vector =
    kernel().point_to_vector_d_object();
  typename R::Vector_to_point_d to_point =
    kernel().vector_to_point_d_object();
  typename R::Construct_vector_d create =
    kernel().construct_vector_d_object();
  Vector_d center_of_hull_facet = create(dimension(),NULL_VECTOR);
  for (int i = 1; i <= current_dimension(); i++) {
    center_of_hull_facet = center_of_hull_facet +
    to_vector(point_of_simplex(selected_hull_simplex,i));
  }
  Point_d center_of_hull_facet_point =
    to_point(center_of_hull_facet/RT(dcur));

  /* we set up the ray from the center to the center of the hull facet */

  Ray_d l(centerpoint, center_of_hull_facet_point);

  /* and check whether it intersects the interior of any hull facet */

  typename R::Contained_in_simplex_d contained_in_simplex =
    kernel().contained_in_simplex_d_object();
  typename R::Intersect_d intersect =
    kernel().intersect_d_object();

  forall_rc_simplices(S,*this) {
    if ( is_unbounded_simplex(S) && S != selected_hull_simplex ) {
      Point_d p; Object op;
      Hyperplane_d  h = S->hyperplane_of_base_facet();
      if ( (op = intersect(l,h), assign(p,op)) ) {
        if ( contained_in_simplex(S->points_begin()+1,
               S->points_begin()+1+current_dimension(),p) ) {
          if ( throw_exceptions )
            throw chull_has_double_coverage();
          return false;
        }
      }
    }
  }
  return true;
}

template <class R>
void Convex_hull_d<R>::
visible_facets_search(Simplex_handle S, const Point_d& x,
                      std::list< Facet_handle >& VisibleFacets,
                      std::size_t& num_of_visited_facets) const
{
  ++num_of_visited_facets;
  S->visited() = true; // we have visited S and never come back ...
  for(int i = 1; i <= current_dimension(); ++i) {
    Simplex_handle T = opposite_simplex(S,i); // for all neighbors T of S
    if ( !T->visited() ) {
      typename R::Oriented_side_d side_of =
        kernel().oriented_side_d_object();
      int side = side_of(T->hyperplane_of_base_facet(),x);
      CGAL_assertion( is_unbounded_simplex(T) );
      if ( side == ON_POSITIVE_SIDE ) {
        VisibleFacets.push_back(T);
        visible_facets_search(T,x,VisibleFacets,num_of_visited_facets);
        // do the recursive search
      }
    } // end visited
  } // end for
}


} //namespace CGAL
#endif // CGAL_CONVEX_HULL_D_H