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

1043 lines
39 KiB
C++
Raw Blame History

// ============================================================================
//
// Copyright (c) 1997-2000 The CGAL Consortium
//
// This software and related documentation is part of an INTERNAL release
// of the Computational Geometry Algorithms Library (CGAL). It is not
// intended for general use.
//
// ----------------------------------------------------------------------------
//
// release : $CGAL_Revision$
// release_date : $CGAL_Date$
//
// file : include/CGAL/Convex_hull_d.h
// package : Kernel_d
// chapter : Basic
//
// source : ddgeo/Convex_hull_d.lw
// revision : $Revision$
// revision_date : $Date$
//
// author(s) : Michael Seel <seel@mpi-sb.mpg.de>
// maintainer : Michael Seel <seel@mpi-sb.mpg.de>
// coordinator : Susan Hert <hert@mpi-sb.mpg.de>
//
// implementation: Higher dimensional geometry
// ============================================================================
//---------------------------------------------------------------------
// 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
/*{\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\footnote{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/Hash_map.h>
#include <CGAL/Regular_complex_d.h>
#include <list>
#include <vector>
#undef _DEBUG
#define _DEBUG 93
#include <CGAL/Kernel_d/debug.h>
CGAL_BEGIN_NAMESPACE
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:
/*{\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 typename Point_list::const_iterator Point_const_iterator;
/*{\Mtypemember const iterator for points.}*/
/*{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 the convex hull object is refereced as constant.}*/
#define USING(t) typedef typename Base::t t
USING(Simplex_const_iterator);USING(Vertex_const_iterator);
USING(Simplex_const_handle);USING(Vertex_const_handle);
#undef USING
typedef Simplex_const_handle Facet_const_handle;
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:
int num_of_bounded_simplices;
int num_of_unbounded_simplices;
int num_of_visibility_tests;
int num_of_points;
int 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,
int& 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 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. }*/
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.}*/
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_points = num_of_vertices = 0;
num_of_unbounded_simplices = num_of_bounded_simplices = 0;
num_of_visibility_tests = 0;
}
int number_of_vertices() const
/*{\Mop returns the number of vertices of |\Mvar|.}*/
{ return num_of_vertices; }
int number_of_facets() const
/*{\Mop returns the number of facets of |\Mvar|.}*/
{ return num_of_unbounded_simplices; }
int number_of_simplices() const
/*{\Mop returns the number of bounded simplices of |\Mvar|.}*/
{ return num_of_bounded_simplices; }
#define STATISTIC(t) std::cout << #t << " = " << t << std::endl
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 (" << dimension() << ") - statistic"
<< std::endl;
STATISTIC(num_of_points);
STATISTIC(num_of_vertices);
STATISTIC(num_of_unbounded_simplices);
STATISTIC(num_of_bounded_simplices);
STATISTIC(num_of_visibility_tests);
}
#undef STATISTIC
bool is_valid() const;
/*{\Mop checks the validity of the data structure.}*/
/*{\Mtext \headerline{Lists and Iterators}
\setopdims{3.5cm}{3.5cm}}*/
Point_const_iterator points_begin() const
/*{\Mop returns the start iterator for points in |\Mvar|.}*/
{ return all_pnts_.begin(); }
Point_const_iterator points_end() const
/*{\Mop returns the past the end iterator for points in |\Mvar|.}*/
{ return all_pnts_.end(); }
Vertex_iterator vertices_begin() { return Base::vertices_begin(); }
/*{\Mop the first vertex 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 the first simplex of |\Mvar|.}*/
Simplex_iterator simplices_end() { return Base::simplices_end(); }
/*{\Mop the past the end iterator for simplices.}*/
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(); }
/* 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) {
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(); }
template <typename H>
struct list_collector {
std::list<H>& L_;
list_collector(std::list<H>& L) : L_(L) {}
void operator()(Facet_handle f) const
{ L_.push_back(static_cast<H>(f)); }
};
std::list<Facet_handle> all_facets()
/*{\Mop returns a list of all facets of |\Mvar|.}*/
{ std::list<Facet_handle> L;
list_collector<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> 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)
/*{\Mtext
\headerline{Iteration Statements}
{\bf forall\_ch\_simplices}($s,C$)
$\{$ ``the simplices of $C$ are successively assigned to $s$'' $\}$
{\bf forall\_ch\_vertices}($v,C$)
$\{$ ``the vertices of $C$ are successively assigned to $v$'' $\}$
}*/
/*{\Mimplementation The implementation of type |\Mtype| is based on
\cite{CMS} and \cite{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_points = 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);
num_of_points++;
if (current_dimension() == -1) { // |x| is the first point to be inserted
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;
int 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 += NewSimplices.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());
#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,
int& 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
int 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
int 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_assertion(0); 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 stands 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() const
{
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) {
CGAL_assertion_msg(0,"check: center on wrong side of a hull facet");
}
}
/* 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) {
CGAL_assertion_msg(0,"check: detected local non-convexity.");
}
}
}
}
/* 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)) ) {
CGAL_assertion_msg( !contained_in_simplex(
S->points_begin()+1,S->points_begin()+1+current_dimension(),p),
"check: current hull has double coverage.");
}
}
}
return true;
}
CGAL_END_NAMESPACE
#endif // CGAL_CONVEX_HULL_D_H
Reference in New Issue View Git Blame Copy Permalink