diff --git a/.gitattributes b/.gitattributes index 9fe0b687eb3..58e3416e03d 100644 --- a/.gitattributes +++ b/.gitattributes @@ -1417,9 +1417,16 @@ Convex_decomposition_3/test/Convex_decomposition_3/star.nef3 -text Convex_hull_2/demo/Convex_hull_2/help/index.html svneol=native#text/html Convex_hull_2/doc_tex/Convex_hull_2/convex_hull.png -text Convex_hull_2/doc_tex/Convex_hull_2/saarhull.png -text svneol=unset#image/png +Convex_hull_3/benchmark/Convex_hull_3/compare_different_approach.cpp -text +Convex_hull_3/benchmark/Convex_hull_3/is_on_positive_side.cpp -text Convex_hull_3/demo/Convex_hull_3/CMakeLists.txt -text Convex_hull_3/doc_tex/Convex_hull_3/bunny.png -text Convex_hull_3/doc_tex/Convex_hull_3/bunny.wrl.gz -text +Convex_hull_3/doc_tex/Convex_hull_3/chull_bimba.png -text +Convex_hull_3/doc_tex/Convex_hull_3_ref/convex_hull_3_to_polyhedron_3.tex -text +Convex_hull_3/examples/Convex_hull_3/incremental_hull_class_3.cpp -text +Convex_hull_3/include/CGAL/convex_hull_3_to_polyhedron_3.h -text +Convex_hull_3/test/Convex_hull_3/quick_hull_default_traits.cpp -text Developers_manual/doc_tex/Developers_manual/fig/Cartesian_ipoint.gif -text svneol=unset#image/gif Developers_manual/doc_tex/Developers_manual/fig/Cartesian_orientation.png -text svneol=unset#image/png Developers_manual/doc_tex/Developers_manual/fig/Object.gif -text svneol=unset#image/gif diff --git a/Convex_hull_3/benchmark/Convex_hull_3/compare_different_approach.cpp b/Convex_hull_3/benchmark/Convex_hull_3/compare_different_approach.cpp new file mode 100644 index 00000000000..145243337aa --- /dev/null +++ b/Convex_hull_3/benchmark/Convex_hull_3/compare_different_approach.cpp @@ -0,0 +1,95 @@ +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + + +typedef CGAL::Exact_predicates_inexact_constructions_kernel K; +typedef CGAL::Exact_predicates_exact_constructions_kernel EK; +typedef CGAL::Polyhedron_3 Polyhedron_3; +typedef CGAL::Polyhedron_3 EK_Polyhedron_3; +typedef K::Segment_3 Segment_3; +typedef CGAL::Delaunay_triangulation_3 Delaunay; + +// define point creator +typedef K::Point_3 Point_3; +typedef CGAL::Creator_uniform_3 PointCreator; + +void load_from_file(const char* path,std::vector& points) +{ + std::ifstream infile (path); + std::size_t nbpt; + infile >> nbpt; + points.reserve(nbpt); + Point_3 p; + while (--nbpt>0) + { + infile >> p; + points.push_back(p); + } +} + +int main(int argc,char** argv) +{ + std::vector points; + if (argc==1){ + CGAL::Random_points_in_sphere_3 gen(1.0); + int nbpt=1000000; + CGAL::copy_n( gen, nbpt, std::back_inserter(points) ); + std::cout << "Using " << 1000000 << " random points in the unit ball\n"; + } + else{ + load_from_file(argv[1],points); + std::cout << "Using a model with " << points.size() << " points.\n"; + } + + Polyhedron_3 poly; + + // compute convex hull + CGAL::Timer time; + time.start(); + CGAL::convex_hull_3(points.begin(), points.end(), poly); + time.stop(); + std::cout << "Static " << time.time() <<" "<< poly.size_of_vertices() << std::endl; + + poly.clear(); + + time.reset(); + time.start(); + Delaunay T(points.begin(), points.end()); + time.stop(); + std::cout << "Delaunay " << time.time() << std::endl; + time.start(); + CGAL::convex_hull_3_to_polyhedron_3(T,poly); + time.stop(); + std::cout << "Delaunay+to_poly " << time.time() <<" "<< poly.size_of_vertices() << std::endl; + poly.clear(); + + time.reset(); + time.start(); + CGAL::convex_hull_incremental_3( points.begin(), points.end(), poly, false); + time.stop(); + std::cout << "incremental EPIC " << time.time() <<" "<< poly.size_of_vertices() << std::endl; + + EK_Polyhedron_3 poly2; + std::vector ek_points; + ek_points.reserve(points.size()); + CGAL::Cartesian_converter convert; + for (std::vector::iterator it=points.begin();it!=points.end();++it){ + ek_points.push_back(convert(*it)); + } + time.reset(); + time.start(); + CGAL::convex_hull_incremental_3( ek_points.begin(), ek_points.end(), poly2, false); + time.stop(); + std::cout << "incremental EPEC " << time.time() <<" "<< poly2.size_of_vertices() << std::endl; + + return 0; +} diff --git a/Convex_hull_3/benchmark/Convex_hull_3/is_on_positive_side.cpp b/Convex_hull_3/benchmark/Convex_hull_3/is_on_positive_side.cpp new file mode 100644 index 00000000000..a9e37bcdd28 --- /dev/null +++ b/Convex_hull_3/benchmark/Convex_hull_3/is_on_positive_side.cpp @@ -0,0 +1,182 @@ +#include +#include +#include +#include +#include + + +#include + +namespace CGAL{ + +template +struct Is_on_positive_side_of_plane_3; + +template +struct Is_on_positive_side_of_plane_3{ + typedef typename Kernel::Point_3 Point_3; + Point_3 p,q,r; + typename Kernel::Orientation_3 orientation; +public: + Is_on_positive_side_of_plane_3(const Kernel& kernel,const Point_3& p_,const Point_3& q_,const Point_3& r_) + :p(p_),q(q_),r(r_),orientation(kernel.orientation_3_object()) {} + + bool operator() (const Point_3& s) const + { + return orientation(p,q,r,s) == CGAL::POSITIVE; + } +}; + +template +struct Is_on_positive_side_of_plane_3{ + typedef typename Kernel::Point_3 Point_3; + double m10,m20,m21,Maxx,Maxy,Maxz; + const Point_3& p_; + + int static_filtered(double psx,double psy, double psz) const{ + using std::fabs; + + // Then semi-static filter. + double maxx = Maxx; + if (maxx < fabs(psx)) maxx = fabs(psx); + double maxy = Maxy; + if (maxy < fabs(psy)) maxy = fabs(psy); + double maxz = Maxz; + if (maxz < fabs(psz)) maxz = fabs(psz); + double det = psx*m10 - m20*psy + m21*psz; + + // Sort maxx < maxy < maxz. + if (maxx > maxz) + std::swap(maxx, maxz); + if (maxy > maxz) + std::swap(maxy, maxz); + else if (maxy < maxx) + std::swap(maxx, maxy); + + // Protect against underflow in the computation of eps. + if (maxx < 1e-97) /* cbrt(min_double/eps) */ { + if (maxx == 0) + return 0; + } + // Protect against overflow in the computation of det. + else if (maxz < 1e102) /* cbrt(max_double [hadamard]/4) */ { + double eps = 5.1107127829973299e-15 * maxx * maxy * maxz; + if (det > eps) return 1; + if (det < -eps) return -1; + } + return 555; + } +public: + Is_on_positive_side_of_plane_3(const Kernel&,const Point_3& p,const Point_3& q,const Point_3& r):p_(p) + { + double pqx = q.x() - p.x(); + double pqy = q.y() - p.y(); + double pqz = q.z() - p.z(); + double prx = r.x() - p.x(); + double pry = r.y() - p.y(); + double prz = r.z() - p.z(); + + m10 = pqy*prz - pry*pqz; + m20 = pqx*prz - prx*pqz; + m21 = pqx*pry - prx*pqy; + + Maxx = fabs(pqx); + if (Maxx < fabs(prx)) Maxx = fabs(prx); + Maxy = fabs(pqy); + if (Maxy < fabs(pry)) Maxy = fabs(pry); + Maxz = fabs(pqz); + if (Maxz < fabs(prz)) Maxz = fabs(prz); + } + + bool operator() (const Point_3& s) const + { + double psx = s.x() - p_.x(); + double psy = s.y() - p_.y(); + double psz = s.z() - p_.z(); + + int static_res = static_filtered(psx,psy,psz); + if (static_res != 555) + return static_res == 1; + + std::cerr << "ERROR static predicate failure!!!\n"; + exit(EXIT_FAILURE); + } +}; + + +template +struct Is_on_positive_side_of_plane_3{ + typedef Simple_cartesian CK; + typedef typename Kernel::Point_3 Point_3; + + Cartesian_converter to_CK; + typename CK::Plane_3 ck_plane; + + Is_on_positive_side_of_plane_3(const Kernel&,const Point_3& p,const Point_3& q,const Point_3& r): + ck_plane(to_CK(p),to_CK(q),to_CK(r)) + {} + + bool operator() (const Point_3& s) const + { + try{ + return ck_plane.has_on_positive_side(to_CK(s)); + } + catch (Uncertain_conversion_exception){ + std::cerr << "ERROR Interval filtering failure\n"; + exit(EXIT_FAILURE); + } + } +}; + +}//namespace CGAL + + +typedef CGAL::Exact_predicates_inexact_constructions_kernel K; +typedef CGAL::Is_on_positive_side_of_plane_3 Pred_using_kernel; +typedef CGAL::Is_on_positive_side_of_plane_3 Pred_using_optimized_static; +typedef CGAL::Is_on_positive_side_of_plane_3 Pred_using_intervals; + + +template +void run(const K::Point_3& p,const K::Point_3& q,const K::Point_3& r,const std::vector& queries,Outputiterator out){ + CGAL::Timer time; time.start(); + Predicate predicate(K(),p,q,r); + for (std::vector::const_iterator it=queries.begin();it!=queries.end();++it) + *out++=predicate(*it); + time.stop(); + std::cout << time.time() << std::endl; +} + +int main() +{ + std::size_t nb_pts=20000000; + typedef CGAL::Creator_uniform_3 Creator; + CGAL::Random_points_in_sphere_3 gen(1); + + std::vector points; + points.reserve(nb_pts); + + K::Point_3 p=*gen++,q=*gen++,r=*gen++; + CGAL::copy_n(gen,nb_pts,std::back_inserter(points)); + + std::vector res0; res0.reserve(nb_pts); + std::vector res1; res1.reserve(nb_pts); + std::vector res2; res2.reserve(nb_pts); + + std::cout << "Running kernel predicates: "; + run(p,q,r,points,std::back_inserter(res0)); + std::cout << "Running static optimized predicates: "; + run(p,q,r,points,std::back_inserter(res1)); + std::cout << "Running predicates with intervals: "; + run(p,q,r,points,std::back_inserter(res2)); + + for (std::size_t k=0;k} should be used +(\ccc{R} being the input kernel). Note that the default traits class takes this into +account. \subsection{Convexity Checking} @@ -45,7 +48,7 @@ resulting hull is a segment or a polyhedron. The function \ccc{convex_hull_incremental_3} % \ccIndexMainItem[C]{convex_hull_incremental_3} provides an interface similar to \ccc{convex_hull_3} for the $d$-dimensional -incremental construction algorithm \cite{cms-frric-93}. +incremental construction algorithm \cite{cms-frric-93} implemented by the class \ccc{CGAL::Convex_hull_d} that is specialized to three dimensions. This function accepts an iterator range over a set of input points and returns a polyhedron, but it does not have a traits class @@ -53,9 +56,13 @@ in its interface. It uses the kernel class \ccc{Kernel} used in the polyhedron type to define an instance of the adapter traits class \ccc{CGAL::Convex_hull_d_traits_3}. -In most cases, the function \ccc{convex_hull_3} will be faster than -\ccc{convex_hull_incremental_3}. The latter is provided mainly -for comparison purposes. +In almost all cases, the static and the dynamic version will +be faster than the incremental convex hull algorithm (mainly +because of the lack of efficient filtering and the overhead +of the general d-dimension). The incremental version is provided for +completeness and educational purposes. You should use the dynamic +version when you need an efficient incremental convex hull algorithm. + To use the full functionality available with the $d$-dimensional class \ccc{CGAL::Convex_hull_d} in three dimensions (\textit{e.g.}, the ability @@ -66,7 +73,7 @@ example. \subsection{Example} -\ccIncludeExampleCode{Convex_hull_3/incremental_hull_3_demo.cpp} +\ccIncludeExampleCode{Convex_hull_3/incremental_hull_class_3.cpp} \section{Dynamic Convex Hull Construction} \ccIndexSubitem{convex hull, 3D}{dynamic} @@ -87,4 +94,33 @@ not all of them are vertices of the hull. \subsection{Example} \ccIncludeExampleCode{Convex_hull_3/dynamic_hull_3.cpp} +\section{Performance} +\begin{figure} +\begin{ccTexOnly} +\begin{center} +\includegraphics[width=12cm]{Convex_hull_3/chull_bimba.png} +\end{center} +\end{ccTexOnly} +\begin{ccHtmlOnly} +
+the convex hull of the bimba model +
+\end{ccHtmlOnly} +\caption{The convex hull of a model made of 192135 points. +\label{fig-ch-bimba}} +\end{figure} + +In the following, we compare the running times of the three approaches to compute 3D convex hulls. +For the static version (using \ccc{CGAL::convex_hull_3}) and the dynamic version +(using \ccc{CGAL::Delaunay_triangulation_3} and \ccc{CGAL::convex_hull_3_to_polyhedron_3}), the kernel +used was \ccc{CGAL::Exact_predicates_inexact_constructions_kernel}. For the incremental version +(using \ccc{CGAL::convex_hull_incremental_3}), the kernel used was \ccc{CGAL::Exact_predicates_exact_constructions_kernel}. + +To compute the convex hull of a million of random points in a unit ball the static approach needed 1.63s, while +the dynamic and incremental approaches needed 9.50s and 11.54s respectively. +To compute the convex hull of the model of Figure \ref{fig-ch-bimba} featuring 192135 points, +the static approach needed 0.18s, while the dynamic and incremental approaches needed 1.90s and 6.80s respectively. + +The measurements have been performed using \cgal\ 3.9, using the \gnu\ \CC\ compiler version 4.3.5, under Linux (Debian distribution), +with the compilation options \texttt{-O3 -DCGAL\_NDEBUG}. The computer used was equipped with a 64bit Intel Xeon 2.27GHz processor and 12GB of RAM. \ No newline at end of file diff --git a/Convex_hull_3/doc_tex/Convex_hull_3_ref/ConvexHullTraits_3.tex b/Convex_hull_3/doc_tex/Convex_hull_3_ref/ConvexHullTraits_3.tex index e32d7ea1830..efd62c4723c 100644 --- a/Convex_hull_3/doc_tex/Convex_hull_3_ref/ConvexHullTraits_3.tex +++ b/Convex_hull_3/doc_tex/Convex_hull_3_ref/ConvexHullTraits_3.tex @@ -63,26 +63,26 @@ provides \ccc{bool operator()(Plane_3 p, Point_3 q, Point_3 r)}, which returns true iff the signed distance from \ccc{q} to \ccc{p} is smaller than the signed distance from \ccc{r} to \ccc{p}} -To handle the degenerate case when all points are coplanar, the following -three types that are default-constructable are necessary: +%To handle the degenerate case when all points are coplanar, the following +%three types that are default-constructable are necessary: +% +%\ccNestedType{Traits_xy}{A model of \ccc{ConvexHullTraits_2} for points projected +% into the $xy$-plane} +%\ccNestedType{Traits_xz}{A model of \ccc{ConvexHullTraits_2} for points projected +% into the $xz$-plane} +%\ccNestedType{Traits_yz}{A model of \ccc{ConvexHullTraits_2} for points projected +% into the $yz$-plane} -\ccNestedType{Traits_xy}{A model of \ccc{ConvexHullTraits_2} for points projected - into the $xy$-plane} -\ccNestedType{Traits_xz}{A model of \ccc{ConvexHullTraits_2} for points projected - into the $xz$-plane} -\ccNestedType{Traits_yz}{A model of \ccc{ConvexHullTraits_2} for points projected - into the $yz$-plane} +%One also needs the following function object to help choose which of the above +%traits classes to use: +% +%\ccNestedType{Max_coordinate_3}{Function object type that provides +%\ccc{int operator()(Vector_3 v)}, which returns the index (0, 1, or 2 for +%$x$, $y$, or $z$, respectively) of the coordinate of $v$ with maximum absolute +%value.} -One also needs the following function object to help choose which of the above -traits classes to use: - -\ccNestedType{Max_coordinate_3}{Function object type that provides -\ccc{int operator()(Vector_3 v)}, which returns the index (0, 1, or 2 for -$x$, $y$, or $z$, respectively) of the coordinate of $v$ with maximum absolute -value.} - -These types need not be provided when it is known that the points are -not all coplanar. +%These types need not be provided when it is known that the points are +%not all coplanar. \ccCreation \ccCreationVariable{traits} @@ -103,7 +103,8 @@ predicate object type. For example: \ccAutoIndexingOn \ccHasModels -\ccRefIdfierPage{CGAL::Convex_hull_traits_3} +\ccRefIdfierPage{CGAL::Convex_hull_traits_3}\\ +All kernels of CGAL %\ccSeeAlso diff --git a/Convex_hull_3/doc_tex/Convex_hull_3_ref/Convex_hull_traits_3.tex b/Convex_hull_3/doc_tex/Convex_hull_3_ref/Convex_hull_traits_3.tex index 1a86a05ac8d..63af9466c22 100644 --- a/Convex_hull_3/doc_tex/Convex_hull_3_ref/Convex_hull_traits_3.tex +++ b/Convex_hull_3/doc_tex/Convex_hull_3_ref/Convex_hull_traits_3.tex @@ -16,7 +16,8 @@ The class \ccRefName\ serves as a traits class for the function \ccc{convex_hull_3}. This is the default traits class for this -function. +function when \ccc{R} is a kernel with exact predicates but inexact constructions +(note below that the type \ccc{Plane_3} is a triple of \ccc{Point_3} and not \ccc{R::Plane_3}). \ccInclude{CGAL/Convex_hull_traits_3.h} @@ -101,7 +102,7 @@ function. \ccCreation \ccCreationVariable{traits} %% choose variable name -\ccConstructor{Convex_hull_traits_3(Convex_hull_traits_2& t);}% +\ccConstructor{Convex_hull_traits_3(Convex_hull_traits_3& t);}% {copy constructor.} \ccOperations diff --git a/Convex_hull_3/doc_tex/Convex_hull_3_ref/IsStronglyConvexTraits_3.tex b/Convex_hull_3/doc_tex/Convex_hull_3_ref/IsStronglyConvexTraits_3.tex index ba2da15fa49..f9504e0244a 100644 --- a/Convex_hull_3/doc_tex/Convex_hull_3_ref/IsStronglyConvexTraits_3.tex +++ b/Convex_hull_3/doc_tex/Convex_hull_3_ref/IsStronglyConvexTraits_3.tex @@ -62,6 +62,7 @@ predicate object type. For example: \ccHasModels \ccRefIdfierPage{CGAL::Convex_hull_traits_3} \\ +All kernels of CGAL %\ccc{CGAL::Kernel_traits_3} \ccSeeAlso diff --git a/Convex_hull_3/doc_tex/Convex_hull_3_ref/convex_hull_3.tex b/Convex_hull_3/doc_tex/Convex_hull_3_ref/convex_hull_3.tex index d7839bcd4e3..11408b3bbcc 100644 --- a/Convex_hull_3/doc_tex/Convex_hull_3_ref/convex_hull_3.tex +++ b/Convex_hull_3/doc_tex/Convex_hull_3_ref/convex_hull_3.tex @@ -47,9 +47,10 @@ Both functions require the following: \begin{enumerate} \item \ccc{InputIterator::value_type} is equivalent to \ccc{Traits::Point_3}. \item \ccc{Traits} is a model of the concept \ccc{ConvexHullTraits_3} - \ccIndexMainItem[c]{ConvexHullTraits_3}. When it is known that - the input points are not all coplanar, the types \ccc{Traits_xy}, - \ccc{Traits_yx}, and \ccc{Traits_yz} need not be provided. + \ccIndexMainItem[c]{ConvexHullTraits_3}. + %When it is known that + %the input points are not all coplanar, the types \ccc{Traits_xy}, + %\ccc{Traits_yx}, and \ccc{Traits_yz} need not be provided. For the purposes of checking the postcondition that the convex hull is valid, \ccc{Traits} should also be a model of the concept \ccc{IsStronglyConvexTraits_3}. @@ -67,9 +68,11 @@ and for the second, it is required that \ccc{ConvexHullPolyhedron_3}. \end{itemize} -The default traits class for both versions of \ccc{convex_hull_3} is -\ccc{Convex_hull_traits_3},% -with the representation \ccc{R} determined by \ccc{InputIterator::value_type}. +For both versions, if the kernel \ccc{R} of the points determined by \ccc{InputIterator::value_type} +is a kernel with exact predicates but inexact constructions +(in practice we check \ccc{R::Has_filtered_predicates_tag} is \ccc{Tag_true} and \ccc{R::FT} is a floating point type), +then the default traits class of \ccc{convex_hull_3} is \ccc{Convex_hull_traits_3}, and \ccc{R} otherwise. + \ccSeeAlso diff --git a/Convex_hull_3/doc_tex/Convex_hull_3_ref/convex_hull_3_to_polyhedron_3.tex b/Convex_hull_3/doc_tex/Convex_hull_3_ref/convex_hull_3_to_polyhedron_3.tex new file mode 100644 index 00000000000..8a0e83f2fcd --- /dev/null +++ b/Convex_hull_3/doc_tex/Convex_hull_3_ref/convex_hull_3_to_polyhedron_3.tex @@ -0,0 +1,34 @@ +\begin{ccRefFunction}{convex_hull_3_to_polyhedron_3} + +\ccDefinition + +The function \ccRefName\ fills a polyhedron with the convex hull +of a set of 3D points contained into a 3D triangulation of \cgal. + +\ccInclude{CGAL/convex_hull_3_to_polyhedron_3.h} + +\ccFunction{ +template +void convex_hull_3_to_polyhedron_3(const Triangulation_3& T,Polyhedron_3& P); +} +{ +The polyhedron \ccc{P} is cleared and the convex hull of the set of 3D points is stored in \ccc{P}. +The plane equations of each face are not computed. +\ccPrecond{\ccc{T.dimension()}==3}. +} + + +\ccHeading{Requirements} +This function requires the following: +\begin{enumerate} + \item \ccc{Triangulation_3} is a CGAL 3D triangulation. + \item \ccc{Polyhedron_3} is an instantiation of \ccc{CGAL::Polyhedron_3}. +\end{enumerate} + + + +\ccSeeAlso + +\ccRefIdfierPage{CGAL::convex_hull_3} + +\end{ccRefFunction} diff --git a/Convex_hull_3/doc_tex/Convex_hull_3_ref/convex_hull_incremental_3.tex b/Convex_hull_3/doc_tex/Convex_hull_3_ref/convex_hull_incremental_3.tex index 2b4b5ee163e..613c7ff19d8 100644 --- a/Convex_hull_3/doc_tex/Convex_hull_3_ref/convex_hull_incremental_3.tex +++ b/Convex_hull_3/doc_tex/Convex_hull_3_ref/convex_hull_incremental_3.tex @@ -14,8 +14,12 @@ \ccDefinition The function \ccRefName\ computes the convex hull polyhedron from a set -of given three-dimensional points. +of given three-dimensional points. +This function is provided for completeness and educational +purposes. When an efficient incremental implementation is needed, +using \ccc{CGAL::Delaunay_triangulation_3} together with +\ccc{CGAL::convex_hull_3_to_polyhedron_3} is highly recommended. \ccInclude{CGAL/convex_hull_incremental_3.h} diff --git a/Convex_hull_3/doc_tex/Convex_hull_3_ref/intro.tex b/Convex_hull_3/doc_tex/Convex_hull_3_ref/intro.tex index c1086a63c9e..acee79a62a3 100644 --- a/Convex_hull_3/doc_tex/Convex_hull_3_ref/intro.tex +++ b/Convex_hull_3/doc_tex/Convex_hull_3_ref/intro.tex @@ -58,7 +58,8 @@ defining \ccc{CGAL_CH_CHECK_EXPENSIVE}% \ccHeading{Convex Hull Functions} \ccRefIdfierPage{CGAL::convex_hull_3} \\ -\ccRefIdfierPage{CGAL::convex_hull_incremental_3} +\ccRefIdfierPage{CGAL::convex_hull_incremental_3}\\ +\ccRefIdfierPage{CGAL::convex_hull_3_to_polyhedron_3} \ccHeading{Convexity Checking Function} diff --git a/Convex_hull_3/doc_tex/Convex_hull_3_ref/main.tex b/Convex_hull_3/doc_tex/Convex_hull_3_ref/main.tex index 77b9c7ea0c4..b9bb833e6c9 100644 --- a/Convex_hull_3/doc_tex/Convex_hull_3_ref/main.tex +++ b/Convex_hull_3/doc_tex/Convex_hull_3_ref/main.tex @@ -15,6 +15,7 @@ \input{Convex_hull_3_ref/ConvexHullPolyhedron_3.tex} \input{Convex_hull_3_ref/ConvexHullTraits_3.tex} \input{Convex_hull_3_ref/Convex_hull_traits_3.tex} +\input{Convex_hull_3_ref/convex_hull_3_to_polyhedron_3.tex} \input{Convex_hull_3_ref/is_strongly_convex_3.tex} \input{Convex_hull_3_ref/IsStronglyConvexTraits_3.tex} diff --git a/Convex_hull_3/examples/Convex_hull_3/dynamic_hull_3.cpp b/Convex_hull_3/examples/Convex_hull_3/dynamic_hull_3.cpp index 6da2c4e83f5..46585091a6d 100644 --- a/Convex_hull_3/examples/Convex_hull_3/dynamic_hull_3.cpp +++ b/Convex_hull_3/examples/Convex_hull_3/dynamic_hull_3.cpp @@ -1,14 +1,17 @@ #include #include #include +#include +#include #include #include -typedef CGAL::Exact_predicates_inexact_constructions_kernel K; -typedef K::Point_3 Point_3; -typedef CGAL::Delaunay_triangulation_3 Delaunay; -typedef Delaunay::Vertex_handle Vertex_handle; +typedef CGAL::Exact_predicates_inexact_constructions_kernel K; +typedef K::Point_3 Point_3; +typedef CGAL::Delaunay_triangulation_3 Delaunay; +typedef Delaunay::Vertex_handle Vertex_handle; +typedef CGAL::Polyhedron_3 Polyhedron_3; int main() { @@ -34,9 +37,15 @@ int main() v_set_it++; } - vertices.clear(); - T.incident_vertices(T.infinite_vertex(), std::back_inserter(vertices)); + //copy the convex hull of points into a polyhedron and use it + //to get the number of points on the convex hull + Polyhedron_3 chull; + CGAL::convex_hull_3_to_polyhedron_3(T,chull); + std::cout << "After removal of 25 points, there are " - << vertices.size() << " points on the convex hull." << std::endl; + << chull.size_of_vertices() << " points on the convex hull." << std::endl; + + + return 0; } diff --git a/Convex_hull_3/examples/Convex_hull_3/incremental_hull_class_3.cpp b/Convex_hull_3/examples/Convex_hull_3/incremental_hull_class_3.cpp new file mode 100644 index 00000000000..3e823bfe349 --- /dev/null +++ b/Convex_hull_3/examples/Convex_hull_3/incremental_hull_class_3.cpp @@ -0,0 +1,68 @@ +// Copyright (c) 2002-2011 Max Planck Institut fuer Informatik (Germany). +// All rights reserved. +// +// This file is part of CGAL (www.cgal.org); you may redistribute it under +// the terms of the Q Public License version 1.0. +// See the file LICENSE.QPL distributed with CGAL. +// +// Licensees holding a valid commercial license may use this file in +// accordance with the commercial license agreement provided with the software. +// +// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE +// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. +// +// $URL: svn+ssh://sloriot@scm.gforge.inria.fr/svn/cgal/branches/experimental-packages/Faster_static_convex_hull_3/demo/Convex_hull_3/incremental_hull_3_demo.cpp $ +// $Id: incremental_hull_3_demo.cpp 44910 2008-08-12 12:58:18Z spion $ +// +// +// Author(s) : Susan Hert +// + + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#ifdef CGAL_USE_GMP +#include +typedef CGAL::Gmpz RT; +#else +#include +typedef CGAL::MP_Float RT; +#endif + +typedef CGAL::Homogeneous K; +typedef K::Point_3 Point_3; +typedef CGAL::Polyhedron_3< K> Polyhedron_3; + +typedef CGAL::Convex_hull_d_traits_3 Hull_traits_3; +typedef CGAL::Convex_hull_d< Hull_traits_3 > Convex_hull_3; +typedef CGAL::Creator_uniform_3 Creator; + +int main () +{ + Convex_hull_3 CH(3); // create instance of the class with dimension == 3 + + // generate 250 points randomly on a sphere of radius 100 + // and insert them into the convex hull + CGAL::Random_points_in_sphere_3 gen(100); + + for (int i = 0; i < 250 ; i++, ++gen) + CH.insert(*gen); + + assert(CH.is_valid()); + + // define polyhedron to hold convex hull and create it + Polyhedron_3 P; + CGAL::convex_hull_d_to_polyhedron_3(CH,P); + + std::cout << "The convex hull has " << P.size_of_vertices() + << " vertices" << std::endl; + return 0; +} diff --git a/Convex_hull_3/examples/Convex_hull_3/quickhull_3.cpp b/Convex_hull_3/examples/Convex_hull_3/quickhull_3.cpp index a8946950c8b..3605845d243 100644 --- a/Convex_hull_3/examples/Convex_hull_3/quickhull_3.cpp +++ b/Convex_hull_3/examples/Convex_hull_3/quickhull_3.cpp @@ -1,14 +1,13 @@ #include #include #include -#include +#include #include #include typedef CGAL::Exact_predicates_inexact_constructions_kernel K; -typedef CGAL::Convex_hull_traits_3 Traits; -typedef Traits::Polyhedron_3 Polyhedron_3; +typedef CGAL::Polyhedron_3 Polyhedron_3; typedef K::Segment_3 Segment_3; // define point creator @@ -32,10 +31,12 @@ int main() CGAL::convex_hull_3(points.begin(), points.end(), ch_object); // determine what kind of object it is - if (CGAL::object_cast(&ch_object) ) - std::cout << "convex hull is a segment " << std::endl; - else if (CGAL::object_cast(&ch_object) ) - std::cout << "convex hull is a polyhedron " << std::endl; + if ( const Segment_3* segment=CGAL::object_cast(&ch_object) ){ + std::cout << "convex hull is the segment " << *segment << std::endl; + } + else if (const Polyhedron_3* poly = CGAL::object_cast(&ch_object) ) + std::cout << "convex hull is a polyhedron with " + << poly->size_of_vertices() << " vertices" << std::endl; else std::cout << "convex hull error!" << std::endl; diff --git a/Convex_hull_3/include/CGAL/Convex_hull_face_base_2.h b/Convex_hull_3/include/CGAL/Convex_hull_face_base_2.h new file mode 100644 index 00000000000..11de5e7a67f --- /dev/null +++ b/Convex_hull_3/include/CGAL/Convex_hull_face_base_2.h @@ -0,0 +1,70 @@ +// Copyright (c) 2003 INRIA Sophia-Antipolis (France). +// All rights reserved. +// +// This file is part of CGAL (www.cgal.org); you may redistribute it under +// the terms of the Q Public License version 1.0. +// See the file LICENSE.QPL distributed with CGAL. +// +// Licensees holding a valid commercial license may use this file in +// accordance with the commercial license agreement provided with the software. +// +// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE +// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. +// +// $URL$ +// $Id$ +// +// +// Author(s) : Mariette Yvinec,Sylvain Pion + +// face of a triangulation of any dimension <=3 + +#ifndef CGAL_CONVEX_HULL_FACE_BASE_2_H +#define CGAL_CONVEX_HULL_FACE_BASE_2_H + +#include + +namespace CGAL { + +template < typename Info_, typename GT, + typename Fb = Triangulation_face_base_2 > +class Convex_hull_face_base_2 + : public Fb +{ + Info_ _info; +public: + typedef typename Fb::Vertex_handle Vertex_handle; + typedef typename Fb::Face_handle Face_handle; + typedef Info_ Info; + + typename std::list::iterator it; + std::list points; + template < typename TDS2 > + struct Rebind_TDS { + typedef typename Fb::template Rebind_TDS::Other Fb2; + typedef Convex_hull_face_base_2 Other; + }; + + Convex_hull_face_base_2() + : Fb(), _info(0) {} + + Convex_hull_face_base_2(Vertex_handle v0, + Vertex_handle v1, + Vertex_handle v2) + : Fb(v0, v1, v2), _info(0) {} + + Convex_hull_face_base_2(Vertex_handle v0, + Vertex_handle v1, + Vertex_handle v2, + Face_handle n0, + Face_handle n1, + Face_handle n2 ) + : Fb(v0, v1, v2, n0, n1, n2), _info(0) {} + + const Info& info() const { return _info; } + Info& info() { return _info; } +}; + +} //namespace CGAL + +#endif // CGAL_CONVEX_HULL_FACE_BASE_2_H diff --git a/Convex_hull_3/include/CGAL/Convex_hull_traits_3.h b/Convex_hull_3/include/CGAL/Convex_hull_traits_3.h index c8dc3a07be1..0f90e68d2e7 100644 --- a/Convex_hull_3/include/CGAL/Convex_hull_traits_3.h +++ b/Convex_hull_3/include/CGAL/Convex_hull_traits_3.h @@ -21,9 +21,11 @@ #define CGAL_CONVEX_HULL_TRAITS_3_H #include +#include #include #include #include +#include namespace CGAL { template < class R_ > @@ -34,7 +36,7 @@ protected: typedef typename R_::Point_3 Point_3; typedef typename R_::Vector_3 Vector_3; - Point_3 p_, q_, r_; + Point_3 p_, q_, r_; public: typedef R_ R; @@ -115,9 +117,15 @@ public: const Point_3& hp = h.p(); const Point_3& hq = h.q(); const Point_3& hr = h.r(); - typename OldK::Less_signed_distance_to_plane_3 - less_signed_distance_to_plane_3; - return less_signed_distance_to_plane_3(hp, hq, hr, p, q); + //typename OldK::Less_signed_distance_to_plane_3 + // less_signed_distance_to_plane_3; + // return less_signed_distance_to_plane_3(hp, hq, hr, p, q); + return has_smaller_signed_dist_to_planeC3(hp.x(), hp.y(), hp.z(), + hq.x(), hq.y(), hq.z(), + hr.x(), hr.y(), hr.z(), + p.x(), p.y(), p.z(), + q.x(), q.y(), q.z()); + } }; @@ -143,6 +151,14 @@ public: } }; + +template +struct GT3 { + typedef typename GT::Point_3 Point_2; +}; + + + template class Convex_hull_traits_3 { @@ -157,7 +173,6 @@ class Convex_hull_traits_3 typedef CGAL::Polyhedron_3 Polyhedron_3; - typedef typename R::Construct_segment_3 Construct_segment_3; typedef typename R::Construct_ray_3 Construct_ray_3; @@ -186,9 +201,6 @@ class Convex_hull_traits_3 Less_signed_distance_to_plane_3; // required for degenerate case of all points coplanar - typedef Projection_traits_xy_3 Traits_xy; - typedef Projection_traits_xz_3 Traits_xz; - typedef Projection_traits_yz_3 Traits_yz; typedef CGAL::Max_coordinate_3 Max_coordinate_3; // for postcondition checking diff --git a/Convex_hull_3/include/CGAL/convex_hull_3.h b/Convex_hull_3/include/CGAL/convex_hull_3.h index f6e521297b3..fa9f2c3ef48 100644 --- a/Convex_hull_3/include/CGAL/convex_hull_3.h +++ b/Convex_hull_3/include/CGAL/convex_hull_3.h @@ -1,4 +1,4 @@ -// Copyright (c) 2001 Max-Planck-Institute Saarbruecken (Germany). +// Copyright (c) 2001-2011 Max-Planck-Institute Saarbruecken (Germany). // All rights reserved. // // This file is part of CGAL (www.cgal.org); you may redistribute it under @@ -17,23 +17,34 @@ // // Author(s) : Susan Hert // : Amol Prakash +// : Andreas Fabri #ifndef CGAL_CONVEX_HULL_3_H #define CGAL_CONVEX_HULL_3_H - #include -#include #include #include #include +#include +#include +#include #include #include +#include +#include +#include +#include #include #include #include #include +#include #include #include +#include +#include +#include + #ifndef CGAL_CH_NO_POSTCONDITIONS #include @@ -42,10 +53,9 @@ namespace CGAL { - -namespace internal{ namespace CH3 { - - struct Plane_equation { +namespace internal{ namespace Convex_hull_3{ + +struct Plane_equation { template typename Facet::Plane_3 operator()( Facet& f) { typename Facet::Halfedge_handle h = f.halfedge(); @@ -55,8 +65,184 @@ namespace internal{ namespace CH3 { h->next()->next()->vertex()->point()); } }; + +//struct to select the default traits class for computing convex hull +template< class Point_3, + class Is_floating_point=typename boost::is_floating_point::Kernel::FT>::type, + class Has_filtered_predicates_tag=typename Kernel_traits::Kernel::Has_filtered_predicates_tag > +struct Default_traits_for_Chull_3{ + typedef typename Kernel_traits::Kernel type; +}; + +//FT is a floating point type and Kernel is a filtered kernel +template +struct Default_traits_for_Chull_3{ + typedef Convex_hull_traits_3< typename Kernel_traits::Kernel > type; +}; + +template +struct Default_polyhedron_for_Chull_3{ + typedef CGAL::Polyhedron_3 type; +}; + +template +struct Default_polyhedron_for_Chull_3 >{ + typedef typename Convex_hull_traits_3::Polyhedron_3 type; +}; -} } //namespace internal::CH3 +//utility class to select the right version of internal predicate Is_on_positive_side_of_plane_3 +template ::Kernel::FT>::type, + class Has_filtered_predicates_tag=typename Kernel_traits::Kernel::Has_filtered_predicates_tag, + class Has_cartesian_tag=typename Kernel_traits::Kernel::Kernel_tag > +struct Use_advanced_filtering{ + typedef CGAL::Tag_false type; +}; + +template +struct Use_advanced_filtering{ + typedef typename Kernel_traits::Kernel K; + typedef CGAL::Boolean_tag type; +}; + +//Predicates internally used +template ::type > +class Is_on_positive_side_of_plane_3{ + typedef typename Traits::Point_3 Point_3; + typename Traits::Plane_3 plane; + typename Traits::Has_on_positive_side_3 has_on_positive_side; +public: + typedef Protect_FPU_rounding Protector; + + Is_on_positive_side_of_plane_3(const Traits& traits,const Point_3& p,const Point_3& q,const Point_3& r) + :plane(traits.construct_plane_3_object()(p,q,r)),has_on_positive_side(traits.has_on_positive_side_3_object()) {} + + bool operator() (const Point_3& s) const + { + return has_on_positive_side(plane,s); + } +}; + + +//This predicate uses copy of the code from the statically filtered version of +//Orientation_3. The rational is that the plane is a member of the functor +//so optimization are done to avoid doing several time operations on the plane. +//The main operator() first tries the static version of the predicate, then uses +//interval arithmetic (the protector must be created before using this predicate) +//and in case of failure, exact arithmetic is used. +template +class Is_on_positive_side_of_plane_3,Tag_true>{ + typedef Simple_cartesian::Type> PK; + typedef Simple_cartesian CK; + typedef Convex_hull_traits_3 Traits; + typedef typename Traits::Point_3 Point_3; + + Cartesian_converter to_CK; + Cartesian_converter to_PK; + + const Point_3& p,q,r; + mutable typename CK::Plane_3* ck_plane; + mutable typename PK::Plane_3* pk_plane; + + double m10,m20,m21,Maxx,Maxy,Maxz; + + static const int STATIC_FILTER_FAILURE = 555; + + //this function is a made from the statically filtered version of Orientation_3 + int static_filtered(double psx,double psy, double psz) const{ + + // Then semi-static filter. + double apsx = CGAL::abs(psx); + double apsy = CGAL::abs(psy); + double apsz = CGAL::abs(psz); + + double maxx = (Maxx < apsx)? apsx : Maxx; + double maxy = (Maxy < apsy)? apsy : Maxy; + double maxz = (Maxz < apsz)? apsz : Maxz; + + double det = psx*m10 - m20*psy + m21*psz; + + // Sort maxx < maxy < maxz. + if (maxx > maxz) + std::swap(maxx, maxz); + if (maxy > maxz) + std::swap(maxy, maxz); + else if (maxy < maxx) + std::swap(maxx, maxy); + + // Protect against underflow in the computation of eps. + if (maxx < 1e-97) /* cbrt(min_double/eps) */ { + if (maxx == 0) + return 0; + } + // Protect against overflow in the computation of det. + else if (maxz < 1e102) /* cbrt(max_double [hadamard]/4) */ { + double eps = 5.1107127829973299e-15 * maxx * maxy * maxz; + if (det > eps) return 1; + if (det < -eps) return -1; + } + return STATIC_FILTER_FAILURE; + } + +public: + typedef typename Interval_nt_advanced::Protector Protector; + + Is_on_positive_side_of_plane_3(const Traits&,const Point_3& p_,const Point_3& q_,const Point_3& r_) + :p(p_),q(q_),r(r_),ck_plane(NULL),pk_plane(NULL) + { + double pqx = q.x() - p.x(); + double pqy = q.y() - p.y(); + double pqz = q.z() - p.z(); + double prx = r.x() - p.x(); + double pry = r.y() - p.y(); + double prz = r.z() - p.z(); + + + m10 = pqy*prz - pry*pqz; + m20 = pqx*prz - prx*pqz; + m21 = pqx*pry - prx*pqy; + + double aprx = CGAL::abs(prx); + double apry = CGAL::abs(pry); + double aprz = CGAL::abs(prz); + + Maxx = CGAL::abs(pqx); + if (Maxx < aprx) Maxx = aprx; + Maxy = CGAL::abs(pqy); + if (Maxy < apry) Maxy = apry; + Maxz = CGAL::abs(pqz); + if (Maxz < aprz) Maxz = aprz; + } + + ~Is_on_positive_side_of_plane_3(){ + if (ck_plane!=NULL) delete ck_plane; + if (pk_plane!=NULL) delete pk_plane; + } + + bool operator() (const Point_3& s) const + { + double psx = s.x() - p.x(); + double psy = s.y() - p.y(); + double psz = s.z() - p.z(); + + int static_res = static_filtered(psx,psy,psz); + if (static_res != STATIC_FILTER_FAILURE) + return static_res == 1; + + try{ + if (ck_plane==NULL) + ck_plane=new typename CK::Plane_3(to_CK(p),to_CK(q),to_CK(r)); + return ck_plane->has_on_positive_side(to_CK(s)); + } + catch (Uncertain_conversion_exception){ + if (pk_plane==NULL) + pk_plane=new typename PK::Plane_3(to_PK(p),to_PK(q),to_PK(r)); + return pk_plane->has_on_positive_side(to_PK(s)); + } + } +}; + template class Build_coplanar_poly : public Modifier_base { @@ -98,27 +284,21 @@ class Build_coplanar_poly : public Modifier_base { ForwardIterator end; }; - template void coplanar_3_hull(InputIterator first, InputIterator beyond, Plane_3 plane, Polyhedron_3& P, const Traits& traits) { - typedef typename Traits::R R; typedef typename Traits::Point_3 Point_3; + typedef typename Kernel_traits::Kernel R; typedef typename Traits::Vector_3 Vector_3; - typedef typename Traits::Max_coordinate_3 Max_coordinate_3; - - typedef typename Traits::Traits_xy Traits_xy; - typedef typename Traits::Traits_xz Traits_xz; - typedef typename Traits::Traits_yz Traits_yz; - + typedef Max_coordinate_3 Max_coordinate_3; typedef Polyhedron_3 Polyhedron; std::list CH_2; typedef typename std::list::iterator CH_2_iterator; typedef typename Traits::Construct_orthogonal_vector_3 Construct_normal_vec; - Max_coordinate_3 max_coordinate = traits.max_coordinate_3_object(); + Max_coordinate_3 max_coordinate; Construct_normal_vec c_normal = traits.construct_orthogonal_vector_3_object(); @@ -129,19 +309,19 @@ void coplanar_3_hull(InputIterator first, InputIterator beyond, case 0: { convex_hull_points_2(first, beyond, std::back_inserter(CH_2), - Traits_yz()); + Projection_traits_yz_3()); break; } case 1: { convex_hull_points_2(first, beyond, std::back_inserter(CH_2), - Traits_xz()); + Projection_traits_xz_3()); break; } case 2: { convex_hull_points_2(first, beyond, std::back_inserter(CH_2), - Traits_xy()); + Projection_traits_xy_3()); break; } default: @@ -158,340 +338,341 @@ void coplanar_3_hull(InputIterator first, InputIterator beyond, // visible is the set of facets visible from point and reachable from // start_facet. // -template +template void -find_visible_set(const typename Traits::Point_3& point, - Facet_handle start_facet, - std::list& visible, +find_visible_set(TDS_2& tds, + const typename Traits::Point_3& point, + typename TDS_2::Face_handle start, + std::list& visible, + std::map& outside, const Traits& traits) { - typedef typename Facet_handle::value_type Facet; - typedef typename Facet::Halfedge_around_facet_circulator Halfedge_circ; - typedef typename Facet::Halfedge_handle Halfedge_handle; typedef typename Traits::Plane_3 Plane_3; - + typedef typename TDS_2::Face_handle Face_handle; + typedef typename TDS_2::Vertex_handle Vertex_handle; typename Traits::Has_on_positive_side_3 has_on_positive_side = traits.has_on_positive_side_3_object(); + std::vector vertices; + vertices.reserve(10); + int VISITED=1, BORDER=2; visible.clear(); - typename std::list::iterator vis_it; - CGAL::Unique_hash_map visited(false); - visible.push_back(start_facet); - visited[start_facet] = true; - Facet_handle current; + typename std::list::iterator vis_it; + visible.push_back(start); + start->info() = VISITED; + vertices.push_back(start->vertex(0)); + vertices.push_back(start->vertex(1)); + vertices.push_back(start->vertex(2)); + start->vertex(0)->info() = start->vertex(1)->info() = start->vertex(2)->info() = VISITED; + for (vis_it = visible.begin(); vis_it != visible.end(); vis_it++) { - // check all the neighbors of the current facet to see if they have + // check all the neighbors of the current face to see if they have // already been visited or not and if not whether they are visible // or not. - current = *vis_it; - Halfedge_circ hdl_init = (*current).facet_begin(); - Halfedge_circ hdl_curr = hdl_init; - do - { - // the facet on the other side of the current halfedge - Facet_handle f = (*(*hdl_curr).opposite()).facet(); - // if haven't already seen this facet - if ( !visited[f] ) - { - visited[f] = true; - Plane_3 plane; - get_plane(plane,f); - if ( has_on_positive_side(plane, point) ) // is visible - { - visible.push_back(f); - } + + for(int i=0; i < 3; i++) { + // the facet on the other side of the current halfedge + Face_handle f = (*vis_it)->neighbor(i); + // if haven't already seen this facet + if (f->info() == 0) { + f->info() = VISITED; + Plane_3 plane(f->vertex(0)->point(),f->vertex(1)->point(),f->vertex(2)->point()); + int ind = f->index(*vis_it); + if ( has_on_positive_side(plane, point) ){ // is visible + visible.push_back(f); + Vertex_handle vh = f->vertex(ind); + if(vh->info() == 0){ vertices.push_back(vh); vh->info() = VISITED;} + } else { + f->info() = BORDER; + f->vertex(TDS_2::cw(ind))->info() = BORDER; + f->vertex(TDS_2::ccw(ind))->info() = BORDER; + outside.insert(std::make_pair(f->vertex(TDS_2::cw(ind)), + typename TDS_2::Edge(f,ind))); } - hdl_curr++; + } else if(f->info() == BORDER) { + int ind = f->index(*vis_it); + f->vertex(TDS_2::cw(ind))->info() = BORDER; + f->vertex(TDS_2::ccw(ind))->info() = BORDER; + outside.insert(std::make_pair(f->vertex(TDS_2::cw(ind)), + typename TDS_2::Edge(f,ind))); + } } - while (hdl_curr != hdl_init); } + + for(typename std::vector::iterator vit = vertices.begin(); + vit != vertices.end(); + ++vit){ + if((*vit)->info() != BORDER){ + tds.delete_vertex(*vit); + } else { + (*vit)->info() = 0; + } + } + } // using a third template parameter for the point instead of getting it from // the traits class as it should be is required by M$VC6 -template -Point -farthest_outside_point(Facet_handle f_handle, std::list& outside_set, +template +typename std::list::iterator +farthest_outside_point(Face_handle f, std::list& outside_set, const Traits& traits) { - typedef typename std::list::iterator Outside_set_iterator; - typename Traits::Plane_3 plane; - get_plane(plane, f_handle); + typedef typename std::list::iterator Outside_set_iterator; CGAL_ch_assertion(!outside_set.empty()); + + typename Traits::Plane_3 plane(f->vertex(0)->point(),f->vertex(1)->point(),f->vertex(2)->point()); + typename Traits::Less_signed_distance_to_plane_3 less_dist_to_plane = traits.less_signed_distance_to_plane_3_object(); Outside_set_iterator farthest_it = std::max_element(outside_set.begin(), outside_set.end(), boost::bind(less_dist_to_plane, plane, _1, _2)); - - return *farthest_it; + return farthest_it; } -template -void compute_plane_equation(Facet_handle f) -{ - typedef typename Facet_handle::value_type Facet; - typedef typename Facet::Halfedge_handle Halfedge_handle; - typedef typename Facet::Plane_3 Plane_3; - - Halfedge_handle h = (*f).halfedge(); - (*f).plane() = Plane_3(h->opposite()->vertex()->point(), - h->vertex()->point(), - h->next()->vertex()->point()); -} - - -template -void get_plane(Plane& plane, Facet_handle f) -{ - typedef typename Facet_handle::value_type Facet; - typedef typename Facet::Halfedge_handle Halfedge_handle; - - Halfedge_handle h = (*f).halfedge(); - plane = Plane(h->opposite()->vertex()->point(), - h->vertex()->point(), - h->next()->vertex()->point()); -} - -// using a template for the Unique_hash_map is required by M$VC7 -// using a template for the Point type instead of getting it from -// the traits class as it should be is required by M$VC6 -template +template void -partition_outside_sets(const std::list& new_facets, - std::list& vis_outside_set, - UHM& outside_sets, - std::list& pending_facets, - const Traits& traits) +partition_outside_sets(const std::list& new_facets, + std::list& vis_outside_set, + std::list& pending_facets, + const Traits& traits) { - typedef typename UHM::Data Data; - typedef typename Traits::Plane_3 Plane_3; - typename std::list::const_iterator f_list_it; - typename std::list::iterator point_it; - - typename Traits::Has_on_positive_side_3 has_on_positive_side = - traits.has_on_positive_side_3_object(); - - // walk through all the new facets and check each unassigned outside point - // to see if it belongs to the outside set of this new facet. - for (f_list_it = new_facets.begin(); f_list_it != new_facets.end(); - f_list_it++) - { - Plane_3 plane; - get_plane(plane, *f_list_it); - Data& point_list = outside_sets[(*f_list_it)]; - for (point_it = vis_outside_set.begin(); - point_it != vis_outside_set.end();) - { - if ( has_on_positive_side(plane, *point_it) ) - { - point_list.push_back(*point_it); - point_it = vis_outside_set.erase(point_it); - } - else - point_it++; - } - } - // put all the new facets with non-empty outside sets in the pending facets - // list. - for (f_list_it = new_facets.begin(); f_list_it != new_facets.end(); - f_list_it++) - { - if (!outside_sets[*f_list_it].empty()) - pending_facets.push_back(*f_list_it); - } + typename std::list::const_iterator f_list_it; + typename std::list::iterator point_it, to_splice; + // walk through all the new facets and check each unassigned outside point + // to see if it belongs to the outside set of this new facet. + for (f_list_it = new_facets.begin(); (f_list_it != new_facets.end()) && (! vis_outside_set.empty()); + ++f_list_it) + { + Face_handle f = *f_list_it; + Is_on_positive_side_of_plane_3 is_on_positive_side( + traits,f->vertex(0)->point(),f->vertex(1)->point(),f->vertex(2)->point()); + std::list& point_list = f->points; + + for (point_it = vis_outside_set.begin();point_it != vis_outside_set.end();){ + if( is_on_positive_side(*point_it) ) { + to_splice = point_it; + ++point_it; + point_list.splice(point_list.end(), vis_outside_set, to_splice); + } else { + ++point_it; + } + } + if(! point_list.empty()){ + pending_facets.push_back(f); + f->it = boost::prior(pending_facets.end()); + } else { + f->it = pending_facets.end(); + } + } + + + for (; f_list_it != new_facets.end();++f_list_it) + (*f_list_it)->it = pending_facets.end(); } -template +template void -ch_quickhull_3_scan( - Polyhedron_3& P, - std::list& pending_facets, - CGAL::Unique_hash_map >& outside_sets, - const Traits& traits) +ch_quickhull_3_scan(TDS_2& tds, + std::list& pending_facets, + const Traits& traits) { - typedef Polyhedron_3 Polyhedron; - typedef typename Polyhedron::Halfedge_handle Halfedge_handle; - typedef typename Polyhedron::Halfedge_iterator Halfedge_iterator; - typedef typename Polyhedron::Facet_handle Facet_handle; + typedef typename TDS_2::Edge Edge; + typedef typename TDS_2::Face_circulator Face_circulator; + typedef typename TDS_2::Face_handle Face_handle; + typedef typename TDS_2::Vertex_handle Vertex_handle; typedef typename Traits::Point_3 Point_3; typedef std::list Outside_set; typedef typename std::list::iterator Outside_set_iterator; + typedef std::map Border_edges; - typedef typename CGAL::Unique_hash_map >::Data Data; - - std::list visible_set; - typename std::list::iterator vis_set_it; - Outside_set vis_outside_set; - Halfedge_iterator hole_halfedge; - Halfedge_handle new_pt_halfedge; - + std::list visible_set; + typename std::list::iterator vis_set_it; + Outside_set vis_outside_set; + Border_edges border; while (!pending_facets.empty()) { vis_outside_set.clear(); - Facet_handle f_handle = pending_facets.back(); - pending_facets.pop_back(); - Point_3 farthest_pt = - farthest_outside_point(f_handle, outside_sets[f_handle], traits); -#ifdef CGAL_CH_3_WINDOW_DEBUG - window << CGAL::RED; - window << farthest_pt; - cout << "farthest point is in red" << endl; - char ch; - cin >> ch; - CGAL_ch_assertion(P.is_valid(true)); - window.clear(); -#endif + Face_handle f_handle = pending_facets.front(); + + Outside_set_iterator farthest_pt_it = farthest_outside_point(f_handle, f_handle->points, traits); + Point_3 farthest_pt = *farthest_pt_it; + f_handle->points.erase(farthest_pt_it); + find_visible_set(tds, farthest_pt, f_handle, visible_set, border, traits); - find_visible_set(farthest_pt, f_handle, visible_set, traits); // for each visible facet for (vis_set_it = visible_set.begin(); vis_set_it != visible_set.end(); vis_set_it++) { + // add its outside set to the global outside set list - Data& point_list = outside_sets[*vis_set_it]; - vis_outside_set.splice(vis_outside_set.end(), point_list, point_list.begin(), point_list.end()); - // delete this visible facet - P.erase_facet((*(*vis_set_it)).halfedge()); - } -#ifdef CGAL_CH_3_WINDOW_DEBUG - window << CGAL::RED; - window << farthest_pt; - window << CGAL::BLUE; - window << P; - cout << "farthest point is in red" << endl; - cout << "after erasing visibile facets"; - cin >> ch; -#endif - for (hole_halfedge = P.halfedges_begin(); - hole_halfedge != P.halfedges_end() && !(*hole_halfedge).is_border(); - hole_halfedge++) - {} - CGAL_ch_assertion(hole_halfedge->is_border()); - CGAL_ch_assertion(hole_halfedge->next()->is_border()); - // add a new facet and vertex to the surface. This is the first - // new facet to be added. - new_pt_halfedge = P.add_vertex_and_facet_to_border(hole_halfedge, - (*hole_halfedge).next()); - // associate the farthest point with the new vertex. - (*(*new_pt_halfedge).vertex()).point() = farthest_pt; - CGAL_ch_assertion( !new_pt_halfedge->is_border() ); - CGAL_ch_assertion( new_pt_halfedge->opposite()->is_border() ); - - std::list new_facets; - new_facets.push_back(new_pt_halfedge->facet()); - Halfedge_handle start_hole_halfedge = new_pt_halfedge->opposite()->prev(); - CGAL_ch_assertion( start_hole_halfedge->is_border() ); - CGAL_ch_assertion( start_hole_halfedge->vertex()->point() == farthest_pt); - - // need to move to second next halfedge to get to a point where a - // triangular facet can be created - Halfedge_handle curr_halfedge = start_hole_halfedge->next()->next(); - CGAL_ch_assertion( curr_halfedge->is_border() ); - - Halfedge_handle new_halfedge; - - // now walk around all the border halfedges and add a facet incident to - // each one to connect it to the farthest point - while (curr_halfedge->next() != start_hole_halfedge) - { - new_halfedge = - P.add_facet_to_border(start_hole_halfedge, curr_halfedge); - CGAL_ch_assertion( !new_halfedge->is_border() ); - CGAL_ch_assertion( new_halfedge->opposite()->is_border() ); - new_facets.push_back(new_halfedge->facet()); - - // once the new facet is added curr->next() will be the next halfedge - // on this facet (i.e., not a border halfedge), so the next border - // halfedge will be the one after the opposite of the new halfedge - curr_halfedge = new_halfedge->opposite()->next(); - CGAL_ch_assertion( curr_halfedge->is_border() ); - } - // fill in the last triangular hole with a facet - new_halfedge = P.fill_hole(curr_halfedge); - new_facets.push_back(new_halfedge->facet()); -#ifdef CGAL_CH_3_WINDOW_DEBUG - window << CGAL::BLUE; - window << P; - cout << "after filling hole" << endl; - char c; - cin >> c; - CGAL_ch_assertion(P.is_valid(true)); -#endif - - // now partition the set of outside set points among the new facets. - partition_outside_sets(new_facets, vis_outside_set, outside_sets, - pending_facets, traits); - } - -} - -template -void non_coplanar_quickhull_3(std::list& points, - Polyhedron_3& P, const Traits& traits) -{ - typedef Polyhedron_3 Polyhedron; - - typedef typename Polyhedron::Facet_handle Facet_handle; - typedef typename Polyhedron::Facet_iterator Facet_iterator; - - typedef typename Traits::Point_3 Point_3; - typedef typename Traits::Plane_3 Plane_3; - typedef CGAL::Unique_hash_map > - Outside_set_map; - typedef typename std::list::iterator P3_iterator; - - std::list pending_facets; - - Facet_iterator f_it; - - typename Traits::Has_on_positive_side_3 has_on_positive_side = - traits.has_on_positive_side_3_object(); - - Outside_set_map outside_sets; - - // for each facet, look at each unassigned point and decide if it belongs - // to the outside set of this facet. - - for (f_it = P.facets_begin(); f_it != P.facets_end(); f_it++) - { - Plane_3 plane; - get_plane(plane, f_it); - for (P3_iterator point_it = points.begin() ; point_it != points.end(); ) - { - if ( has_on_positive_side(plane, *point_it) ){ - outside_sets[f_it].push_back(*point_it); - point_it = points.erase(point_it); - } else { - ++point_it; + std::list& point_list = (*vis_set_it)->points; + if(! point_list.empty()){ + vis_outside_set.splice(vis_outside_set.end(), point_list, point_list.begin(), point_list.end()); } - + if((*vis_set_it)->it != pending_facets.end()){ + pending_facets.erase((*vis_set_it)->it); + } + (*vis_set_it)->info() = 0; } + + std::vector edges; + edges.reserve(border.size()); + typename Border_edges::iterator it = border.begin(); + Edge e = it->second; + e.first->info() = 0; + edges.push_back(e); + border.erase(it); + while(! border.empty()){ + it = border.find(e.first->vertex(TDS_2::ccw(e.second))); + assert(it != border.end()); + e = it->second; + e.first->info() = 0; + edges.push_back(e); + border.erase(it); + } + + // If we want to reuse the faces we must only pass |edges| many, and call delete_face for the others. + // Also create facets if necessary + std::ptrdiff_t diff = visible_set.size() - edges.size(); + if(diff < 0){ + for(int i = 0; i<-diff;i++){ + visible_set.push_back(tds.create_face()); + } + } else { + for(int i = 0; ipoint() = farthest_pt; + vh->info() = 0; + + // now partition the set of outside set points among the new facets. + + partition_outside_sets(visible_set, vis_outside_set, + pending_facets, traits); + + } +} + +template +void non_coplanar_quickhull_3(std::list& points, + TDS_2& tds, const Traits& traits) +{ + typedef typename Traits::Point_3 Point_3; + typedef typename Traits::Plane_3 Plane_3; + + typedef typename TDS_2::Vertex_handle Vertex_handle; + typedef typename TDS_2::Face_handle Face_handle; + typedef typename TDS_2::Face_iterator Face_iterator; + typedef typename std::list::iterator P3_iterator; + + std::list pending_facets; + + typename Is_on_positive_side_of_plane_3::Protector p; + + // for each facet, look at each unassigned point and decide if it belongs + // to the outside set of this facet. + for(Face_iterator fit = tds.faces_begin(); fit != tds.faces_end(); ++fit){ + Is_on_positive_side_of_plane_3 is_on_positive_side( + traits,fit->vertex(0)->point(),fit->vertex(1)->point(),fit->vertex(2)->point() ); + for (P3_iterator point_it = points.begin() ; point_it != points.end(); ) + { + if( is_on_positive_side(*point_it) ) { + P3_iterator to_splice = point_it; + ++point_it; + fit->points.splice(fit->points.end(), points, to_splice); + } else { + ++point_it; + } + } } // add all the facets with non-empty outside sets to the set of facets for // further consideration - for (f_it = P.facets_begin(); f_it != P.facets_end(); f_it++) - if (!outside_sets[f_it].empty()) - pending_facets.push_back(f_it); - ch_quickhull_3_scan(P, pending_facets, outside_sets, traits); + for(Face_iterator fit = tds.faces_begin(); fit != tds.faces_end(); ++fit){ + if (! fit->points.empty()){ + pending_facets.push_back(fit); + fit->it = boost::prior(pending_facets.end()); + } else { + fit->it = pending_facets.end(); + } + } - CGAL_ch_expensive_postcondition(all_points_inside(points.begin(), - points.end(),P,traits)); - CGAL_ch_postcondition(is_strongly_convex_3(P, traits)); + + ch_quickhull_3_scan(tds, pending_facets, traits); + + //std::cout << "|V(tds)| = " << tds.number_of_vertices() << std::endl; +// CGAL_ch_expensive_postcondition(all_points_inside(points.begin(), +// points.end(),P,traits)); +// CGAL_ch_postcondition(is_strongly_convex_3(P, traits)); } +namespace internal{ + +template +class Build_convex_hull_from_TDS_2 : public CGAL::Modifier_base { + typedef std::map Vertex_map; + + const TDS& t; + template + static unsigned get_vertex_index( Vertex_map& vertex_map, + typename TDS::Vertex_handle vh, + Builder& builder, + unsigned& vindex) + { + std::pair + res=vertex_map.insert(std::make_pair(vh,vindex)); + if (res.second){ + builder.add_vertex(vh->point()); + ++vindex; + } + return res.first->second; + } + +public: + Build_convex_hull_from_TDS_2(const TDS& t_):t(t_) + { + CGAL_assertion(t.dimension()==2); + } + void operator()( HDS& hds) { + // Postcondition: `hds' is a valid polyhedral surface. + typedef typename HDS::Vertex Vertex; + typedef typename Vertex::Point Point; + + CGAL::Polyhedron_incremental_builder_3 B( hds, true); + Vertex_map vertex_map; + //start the surface + B.begin_surface( t.number_of_vertices(), t.number_of_faces()); + unsigned vindex=0; + for (typename TDS::Face_iterator it=t.faces_begin();it!=t.faces_end();++it) + { + unsigned i0=get_vertex_index(vertex_map,it->vertex(0),B,vindex); + unsigned i1=get_vertex_index(vertex_map,it->vertex(1),B,vindex); + unsigned i2=get_vertex_index(vertex_map,it->vertex(2),B,vindex); + B.begin_facet(); + B.add_vertex_to_facet( i0 ); + B.add_vertex_to_facet( i1 ); + B.add_vertex_to_facet( i2 ); + B.end_facet(); + } + B.end_surface(); + } +}; + +} //namespace internal template void @@ -504,6 +685,13 @@ ch_quickhull_polyhedron_3(std::list& points, typedef typename Traits::Plane_3 Plane_3; typedef typename std::list::iterator P3_iterator; + typedef typename Kernel_traits::Kernel R; + typedef Triangulation_data_structure_2< + Triangulation_vertex_base_with_info_2 >, + Convex_hull_face_base_2 > Tds; + typedef typename Tds::Vertex_handle Vertex_handle; + typedef typename Tds::Face_handle Face_handle; + // found three points that are not collinear, so construct the plane defined // by these points and then find a point that has maximum distance from this // plane. @@ -532,34 +720,45 @@ ch_quickhull_polyhedron_3(std::list& points, } else max_it = min_max.second; -#ifdef CGAL_CH_3_WINDOW_DEBUG - window << CGAL::GREEN; - window << *point1_it; - window << *point2_it; - window << *point3_it; - window << CGAL::RED; - window << *max_it; - char ch; - cin >> ch; -#endif // if the maximum distance point is on the plane then all are coplanar if (coplanar(*point1_it, *point2_it, *point3_it, *max_it)) { coplanar_3_hull(points.begin(), points.end(), plane, P, traits); } else { - P.make_tetrahedron(*point1_it, *point2_it, *point3_it, *max_it); - points.erase(point1_it); - points.erase(point2_it); - points.erase(point3_it); - points.erase(max_it); - if (!points.empty()) - non_coplanar_quickhull_3(points, P, traits); + Tds tds; + Vertex_handle v0 = tds.create_vertex(); v0->set_point(*point1_it); + Vertex_handle v1 = tds.create_vertex(); v1->set_point(*point2_it); + Vertex_handle v2 = tds.create_vertex(); v2->set_point(*point3_it); + Vertex_handle v3 = tds.create_vertex(); v3->set_point(*max_it); + + v0->info() = v1->info() = v2->info() = v3->info() = 0; + Face_handle f0 = tds.create_face(v0,v1,v2); + Face_handle f1 = tds.create_face(v3,v1,v0); + Face_handle f2 = tds.create_face(v3,v2,v1); + Face_handle f3 = tds.create_face(v3,v0,v2); + tds.set_dimension(2); + f0->set_neighbors(f2, f3, f1); + f1->set_neighbors(f0, f3, f2); + f2->set_neighbors(f0, f1, f3); + f3->set_neighbors(f0, f2, f1); + + points.erase(point1_it); + points.erase(point2_it); + points.erase(point3_it); + points.erase(max_it); + if (!points.empty()){ + non_coplanar_quickhull_3(points, tds, traits); + internal::Build_convex_hull_from_TDS_2 builder(tds); + P.delegate(builder); + } } - std::transform( P.facets_begin(), P.facets_end(), P.planes_begin(),internal::CH3::Plane_equation()); + std::transform( P.facets_begin(), P.facets_end(), P.planes_begin(),Plane_equation()); } +} } //namespace internal::Convex_hull_3 + template void convex_hull_3(InputIterator first, InputIterator beyond, @@ -643,21 +842,35 @@ convex_hull_3(InputIterator first, InputIterator beyond, return; } - typename Traits::Polyhedron_3 P; - // result will be a polyhedron - ch_quickhull_polyhedron_3(points, point1_it, point2_it, point3_it, - P, traits); + typename internal::Convex_hull_3::Default_polyhedron_for_Chull_3::type P; + + P3_iterator minx, maxx, miny, it; + minx = maxx = miny = it = points.begin(); + ++it; + for(; it != points.end(); ++it){ + if(it->x() < minx->x()) minx = it; + if(it->x() > maxx->x()) maxx = it; + if(it->y() < miny->y()) miny = it; + } + if(! collinear(*minx, *maxx, *miny) ){ + internal::Convex_hull_3::ch_quickhull_polyhedron_3(points, minx, maxx, miny, P, traits); + } else { + internal::Convex_hull_3::ch_quickhull_polyhedron_3(points, point1_it, point2_it, point3_it, P, traits); + } ch_object = make_object(P); + //std::cout << "|V| = " << P.size_of_vertices() << std::endl; } + template void convex_hull_3(InputIterator first, InputIterator beyond, Object& ch_object) { typedef typename std::iterator_traits::value_type Point_3; typedef typename Kernel_traits::Kernel K; - convex_hull_3(first, beyond, ch_object, Convex_hull_traits_3()); + typedef typename internal::Convex_hull_3::Default_traits_for_Chull_3::type Traits; + convex_hull_3(first, beyond, ch_object, Traits()); } @@ -703,8 +916,8 @@ void convex_hull_3(InputIterator first, InputIterator beyond, polyhedron.clear(); // result will be a polyhedron - ch_quickhull_polyhedron_3(points, point1_it, point2_it, point3_it, - polyhedron, traits); + internal::Convex_hull_3::ch_quickhull_polyhedron_3(points, point1_it, point2_it, point3_it, + polyhedron, traits); } @@ -714,8 +927,8 @@ void convex_hull_3(InputIterator first, InputIterator beyond, Polyhedron_3& polyhedron) { typedef typename std::iterator_traits::value_type Point_3; - typedef typename Kernel_traits::Kernel K; - convex_hull_3(first, beyond, polyhedron, Convex_hull_traits_3()); + typedef typename internal::Convex_hull_3::Default_traits_for_Chull_3::type Traits; + convex_hull_3(first, beyond, polyhedron, Traits()); } } // namespace CGAL diff --git a/Convex_hull_3/include/CGAL/convex_hull_3_to_polyhedron_3.h b/Convex_hull_3/include/CGAL/convex_hull_3_to_polyhedron_3.h new file mode 100644 index 00000000000..2bcc402e200 --- /dev/null +++ b/Convex_hull_3/include/CGAL/convex_hull_3_to_polyhedron_3.h @@ -0,0 +1,93 @@ +// Copyright (c) 2011 GeometryFactory (France). +// All rights reserved. +// +// This file is part of CGAL (www.cgal.org); you may redistribute it under +// the terms of the Q Public License version 1.0. +// See the file LICENSE.QPL distributed with CGAL. +// +// Licensees holding a valid commercial license may use this file in +// accordance with the commercial license agreement provided with the software. +// +// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE +// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. +// +// $URL$ +// $Id$ +// +// +// Author(s) : Sebastien Loriot +// + +#include + +#ifndef CGAL_CONVEX_HULL_3_TO_POLYHEDRON_3_H +#define CGAL_CONVEX_HULL_3_TO_POLYHEDRON_3_H + +namespace CGAL { + +template +class Convex_hull_modifier_from_triangulation_3 : public CGAL::Modifier_base { + typedef std::map Vertex_map; + + const Triangulation& t; + template + static unsigned get_vertex_index( Vertex_map& vertex_map, + typename Triangulation::Vertex_handle vh, + Builder& builder, + unsigned& vindex) + { + std::pair + res=vertex_map.insert(std::make_pair(vh,vindex)); + if (res.second){ + builder.add_vertex(vh->point()); + ++vindex; + } + return res.first->second; + } + +public: + Convex_hull_modifier_from_triangulation_3(const Triangulation& t_):t(t_) + { + CGAL_assertion(t.dimension()==3); + } + void operator()( HDS& hds) { + // Postcondition: `hds' is a valid polyhedral surface. + typedef typename HDS::Vertex Vertex; + typedef typename Vertex::Point Point; + + CGAL::Polyhedron_incremental_builder_3 B( hds, true); + std::vector ch_facets; + Vertex_map vertex_map; + t.incident_cells(t.infinite_vertex(),std::back_inserter(ch_facets)); + std::size_t nb_facets=ch_facets.size(); + //start the surface + B.begin_surface( nb_facets, nb_facets); + unsigned vindex=0; + for (typename std::vector::const_iterator it= + ch_facets.begin();it!=ch_facets.end();++it) + { + unsigned ifv_index= (*it)->index(t.infinite_vertex()); + bool is_odd=ifv_index%2==0; + unsigned i0=get_vertex_index(vertex_map,(*it)->vertex((ifv_index + (is_odd?3:1) )%4),B,vindex); + unsigned i1=get_vertex_index(vertex_map,(*it)->vertex((ifv_index + 2 )%4),B,vindex); + unsigned i2=get_vertex_index(vertex_map,(*it)->vertex((ifv_index + (is_odd?1:3) )%4),B,vindex); + B.begin_facet(); + B.add_vertex_to_facet( i0 ); + B.add_vertex_to_facet( i1 ); + B.add_vertex_to_facet( i2 ); + B.end_facet(); + } + B.end_surface(); + } +}; + +template +void convex_hull_3_to_polyhedron_3(const Triangulation_3& T,Polyhedron_3& P){ + P.clear(); + Convex_hull_modifier_from_triangulation_3 modifier(T); + P.delegate(modifier); +} + +} //namespace CGAL + +#endif //CGAL_CONVEX_HULL_3_TO_POLYHEDRON_3_H diff --git a/Convex_hull_3/test/Convex_hull_3/quick_hull_default_traits.cpp b/Convex_hull_3/test/Convex_hull_3/quick_hull_default_traits.cpp new file mode 100644 index 00000000000..69cb843ce32 --- /dev/null +++ b/Convex_hull_3/test/Convex_hull_3/quick_hull_default_traits.cpp @@ -0,0 +1,39 @@ +#include +#include +#include +#include +#include + +#include +#include + +#ifdef CGAL_USE_LEDA +# include +typedef leda_rational Precise_rational; +#elif defined CGAL_USE_GMP +# include +# include +typedef CGAL::Quotient Precise_rational; +#else +# include +# include +typedef CGAL::Quotient Precise_rational; +#endif + +typedef CGAL::Exact_predicates_exact_constructions_kernel EPEC; +typedef CGAL::Exact_predicates_inexact_constructions_kernel EPIC; +typedef CGAL::Simple_cartesian SCD; +typedef CGAL::Simple_cartesian SCR; + +using namespace CGAL::internal::Convex_hull_3; + +int main() +{ + BOOST_STATIC_ASSERT( (boost::is_same::type>::value) ); + BOOST_STATIC_ASSERT( (boost::is_same::type>::value) ); + BOOST_STATIC_ASSERT( (boost::is_same::type>::value) ); + BOOST_STATIC_ASSERT( (boost::is_same::type>::value) ); + BOOST_STATIC_ASSERT( (boost::is_same,Default_traits_for_Chull_3::type>::value) ); + BOOST_STATIC_ASSERT( (boost::is_same >::Protector,CGAL::Protect_FPU_rounding >::value) ); + return 0; +} \ No newline at end of file diff --git a/Convex_hull_3/test/Convex_hull_3/quickhull_test_3.cpp b/Convex_hull_3/test/Convex_hull_3/quickhull_test_3.cpp index 1aaa4b8130d..ed6fbaedf50 100644 --- a/Convex_hull_3/test/Convex_hull_3/quickhull_test_3.cpp +++ b/Convex_hull_3/test/Convex_hull_3/quickhull_test_3.cpp @@ -33,6 +33,33 @@ typedef CGAL::Random_points_in_sphere_3 Generator; const unsigned int num = 40; +template +void compute_plane_equation(Facet_handle f) +{ + typedef typename Facet_handle::value_type Facet; + typedef typename Facet::Halfedge_handle Halfedge_handle; + typedef typename Facet::Plane_3 Plane_3; + + Halfedge_handle h = (*f).halfedge(); + (*f).plane() = Plane_3(h->opposite()->vertex()->point(), + h->vertex()->point(), + h->next()->vertex()->point()); +} + + +template +void get_plane(Plane& plane, Facet_handle f) +{ + typedef typename Facet_handle::value_type Facet; + typedef typename Facet::Halfedge_handle Halfedge_handle; + + Halfedge_handle h = (*f).halfedge(); + plane = Plane(h->opposite()->vertex()->point(), + h->vertex()->point(), + h->next()->vertex()->point()); +} + + void test_tetrahedron_convexity() { Polyhedron_3 P; @@ -45,7 +72,7 @@ void test_tetrahedron_convexity() for( Polyhedron_3::Facet_iterator f = P.facets_begin(); f != P.facets_end(); ++f ) { - CGAL::compute_plane_equation(f); + compute_plane_equation(f); } assert( CGAL::is_strongly_convex_3(P) ); @@ -62,7 +89,7 @@ void test_triangle_convexity() for( Polyhedron_3::Facet_iterator f = P.facets_begin(); f != P.facets_end(); ++f ) { - CGAL::compute_plane_equation(f); + compute_plane_equation(f); } assert( CGAL::is_strongly_convex_3(P) );