Merge branch 'Triangulation-add_regular_tri-cjamin_mglisse-old' into Triangulation-add_regular_tri-cjamin_mglisse

This commit is contained in:
Clement Jamin 2014-07-30 16:42:33 +02:00
commit 9d8d255307
41 changed files with 4555 additions and 2142 deletions

View File

@ -26,6 +26,7 @@
#include <CGAL/NewKernel_d/Kernel_d_interface.h>
#include <CGAL/internal/Exact_type_selector.h>
#include <CGAL/Interval_nt.h>
#include <CGAL/NewKernel_d/Types/Weighted_point.h>
namespace CGAL {

View File

@ -87,6 +87,7 @@ struct Cartesian_LA_base_d : public Dimension_base<Dim_>
::add<Segment_tag>::type
::add<Hyperplane_tag>::type
::add<Sphere_tag>::type
::add<Weighted_point_tag>::type
Object_list;
typedef typeset< Point_cartesian_const_iterator_tag>::type

View File

@ -180,7 +180,7 @@ template<class R_> struct Compute_cartesian_coordinate {
#ifdef CGAL_CXX11
typedef decltype(std::declval<const first_argument_type>()[0]) result_type;
#else
typedef RT const& result_type;
typedef RT result_type;
// RT const& doesn't work with some LA (Eigen2 for instance) so we
// should use plain RT or find a way to detect this.
#endif

View File

@ -24,6 +24,7 @@
#include <CGAL/NewKernel_d/KernelD_converter.h>
#include <CGAL/NewKernel_d/Filtered_predicate2.h>
#include <boost/mpl/if.hpp>
#include <boost/mpl/and.hpp>
namespace CGAL {
@ -45,6 +46,12 @@ struct Cartesian_filter_K : public Base_,
typedef typename Store_kernel2<EK_>::reference2_type EK_rt;
EK_rt exact_kernel()const{return this->kernel2();}
// MSVC is too dumb to perform the empty base optimization.
typedef boost::mpl::and_<
internal::Do_not_store_kernel<Kernel_base>,
internal::Do_not_store_kernel<AK>,
internal::Do_not_store_kernel<EK> > Do_not_store_kernel;
//TODO: C2A/C2E could be able to convert *this into this->kernel() or this->kernel2().
typedef KernelD_converter<Kernel_base,AK> C2A;
typedef KernelD_converter<Kernel_base,EK> C2E;

View File

@ -78,6 +78,10 @@ template<class R_> struct Construct_flat_orientation : private Store_kernel<R_>
// the points are affinely independent.
template<class Iter>
result_type operator()(Iter f, Iter e)const{
/*std::cerr << "Kernel flat orientation - points: " ; // CJTODO DEBUG
for (Iter it = f ; it != e ; ++it)
std::cerr << (*it)[0] << " ";
std::cerr << std::endl; // CJTODO DEBUG*/
Iter f_save = f;
PD pd (this->kernel());
CCC ccc (this->kernel());
@ -89,6 +93,7 @@ template<class R_> struct Construct_flat_orientation : private Store_kernel<R_>
std::vector<int>& rest=o.rest; rest.reserve(dim+1);
for(int i=0; i<dim+1; ++i) rest.push_back(i);
for( ; f != e ; ++col, ++f ) {
//std::cerr << "(*f)[0]=" << (*f)[0] << std::endl;
Point const&p=*f;
// use a coordinate iterator instead?
for(int i=0; i<dim; ++i) coord(col, i) = ccc(p, i);
@ -268,11 +273,61 @@ template<class R_> struct In_flat_side_of_oriented_sphere : private Store_kernel
}
};
template<class R_> struct In_flat_power_test_raw : private Store_kernel<R_> {
CGAL_FUNCTOR_INIT_STORE(In_flat_power_test_raw)
typedef R_ R;
typedef typename Get_type<R, FT_tag>::type FT;
typedef typename Get_type<R, Point_tag>::type Point;
typedef typename Get_type<R, Orientation_tag>::type result_type;
typedef typename Increment_dimension<typename R::Default_ambient_dimension,2>::type D1;
typedef typename Increment_dimension<typename R::Max_ambient_dimension,2>::type D2;
typedef typename R::LA::template Rebind_dimension<D1,D2>::Other LA;
typedef typename LA::Square_matrix Matrix;
template<class Iter, class IterW, class Wt>
result_type operator()(Flat_orientation const&o, Iter f, Iter e, IterW fw, Point const&x, Wt const&w) const {
// TODO: can't work in the projection, but we should at least remove the row of 1s.
typename Get_functor<R, Compute_point_cartesian_coordinate_tag>::type c(this->kernel());
typename Get_functor<R, Point_dimension_tag>::type pd(this->kernel());
int d=pd(*f);
Matrix m(d+2,d+2);
int i=0;
for(;f!=e;++f,++fw,++i) {
Point const& p=*f;
m(i,0)=1;
m(i,d+1)=-*fw;
for(int j=0;j<d;++j){
m(i,j+1)=c(p,j);
m(i,d+1)+=CGAL_NTS square(m(i,j+1));
}
}
for(std::vector<int>::const_iterator it = o.rest.begin(); it != o.rest.end() /* i<d+1 */; ++i, ++it) {
m(i,0)=1;
for(int j=0;j<d;++j){
m(i,j+1)=0; // unneeded if the matrix is initialized to 0
}
if(*it != d) m(i,d+1)=m(i,1+*it)=1;
else m(i,d+1)=0;
}
m(d+1,0)=1;
m(d+1,d+1)=-w;
for(int j=0;j<d;++j){
m(d+1,j+1)=c(x,j);
m(d+1,d+1)+=CGAL_NTS square(m(d+1,j+1));
}
result_type ret = -LA::sign_of_determinant(CGAL_MOVE(m));
if(o.reverse) ret=-ret;
return ret;
}
};
}
CGAL_KD_DEFAULT_TYPE(Flat_orientation_tag,(CGAL::CartesianDKernelFunctors::Flat_orientation),(),());
CGAL_KD_DEFAULT_FUNCTOR(In_flat_orientation_tag,(CartesianDKernelFunctors::In_flat_orientation<K>),(Point_tag),(Compute_point_cartesian_coordinate_tag,Point_dimension_tag));
CGAL_KD_DEFAULT_FUNCTOR(In_flat_side_of_oriented_sphere_tag,(CartesianDKernelFunctors::In_flat_side_of_oriented_sphere<K>),(Point_tag),(Compute_point_cartesian_coordinate_tag,Point_dimension_tag));
CGAL_KD_DEFAULT_FUNCTOR(In_flat_power_test_raw_tag,(CartesianDKernelFunctors::In_flat_power_test_raw<K>),(Point_tag),(Compute_point_cartesian_coordinate_tag,Point_dimension_tag));
CGAL_KD_DEFAULT_FUNCTOR(Construct_flat_orientation_tag,(CartesianDKernelFunctors::Construct_flat_orientation<K>),(Point_tag),(Compute_point_cartesian_coordinate_tag,Point_dimension_tag,In_flat_orientation_tag));
CGAL_KD_DEFAULT_FUNCTOR(Contained_in_affine_hull_tag,(CartesianDKernelFunctors::Contained_in_affine_hull<K>),(Point_tag),(Compute_point_cartesian_coordinate_tag,Point_dimension_tag));
}

View File

@ -93,7 +93,7 @@ class KernelD_converter_
//typedef typename KOC::argument_type K1_Obj;
//typedef typename KOC::result_type K2_Obj;
public:
using Base::operator(); // don't use directly, just make it accessible to the next level
using Base::operator(); // don't use directly, just make it accessible to the next level
K2_Obj helper(K1_Obj const& o,CGAL_BOOSTD true_type)const{
return KOC()(this->myself().kernel(),this->myself().kernel2(),this->myself(),o);
}

View File

@ -54,6 +54,7 @@ template <class Base_> struct Kernel_d_interface : public Base_ {
typedef typename Get_type<Base, Ray_tag>::type Ray_d;
typedef typename Get_type<Base, Iso_box_tag>::type Iso_box_d;
typedef typename Get_type<Base, Aff_transformation_tag>::type Aff_transformation_d;
typedef typename Get_type<Base, Weighted_point_tag>::type Weighted_point_d;
typedef typename Get_functor<Base, Compute_point_cartesian_coordinate_tag>::type Compute_coordinate_d;
typedef typename Get_functor<Base, Compare_lexicographically_tag>::type Compare_lexicographically_d;
typedef typename Get_functor<Base, Equal_points_tag>::type Equal_d;
@ -64,10 +65,12 @@ template <class Base_> struct Kernel_d_interface : public Base_ {
typedef typename Get_functor<Base, Less_point_cartesian_coordinate_tag>::type Less_coordinate_d;
typedef typename Get_functor<Base, Point_dimension_tag>::type Point_dimension_d;
typedef typename Get_functor<Base, Side_of_oriented_sphere_tag>::type Side_of_oriented_sphere_d;
typedef typename Get_functor<Base, Power_test_tag>::type Power_test_d;
typedef typename Get_functor<Base, Contained_in_affine_hull_tag>::type Contained_in_affine_hull_d;
typedef typename Get_functor<Base, Construct_flat_orientation_tag>::type Construct_flat_orientation_d;
typedef typename Get_functor<Base, In_flat_orientation_tag>::type In_flat_orientation_d;
typedef typename Get_functor<Base, In_flat_side_of_oriented_sphere_tag>::type In_flat_side_of_oriented_sphere_d;
typedef typename Get_functor<Base, In_flat_power_test_tag>::type In_flat_power_test_d;
typedef typename Get_functor<Base, Point_to_vector_tag>::type Point_to_vector_d;
typedef typename Get_functor<Base, Vector_to_point_tag>::type Vector_to_point_d;
typedef typename Get_functor<Base, Construct_ttag<Point_tag> >::type Construct_point_d;
@ -80,6 +83,7 @@ template <class Base_> struct Kernel_d_interface : public Base_ {
typedef typename Get_functor<Base, Construct_ttag<Ray_tag> >::type Construct_ray_d;
typedef typename Get_functor<Base, Construct_ttag<Iso_box_tag> >::type Construct_iso_box_d;
typedef typename Get_functor<Base, Construct_ttag<Aff_transformation_tag> >::type Construct_aff_transformation_d;
typedef typename Get_functor<Base, Construct_ttag<Weighted_point_tag> >::type Construct_weighted_point_d;
typedef typename Get_functor<Base, Midpoint_tag>::type Midpoint_d;
struct Component_accessor_d : private Store_kernel<Kernel> {
typedef Kernel R_; // for the macro
@ -145,6 +149,9 @@ template <class Base_> struct Kernel_d_interface : public Base_ {
typedef typename Get_functor<Base, Orthogonal_vector_tag>::type Orthogonal_vector_d;
typedef typename Get_functor<Base, Linear_base_tag>::type Linear_base_d;
typedef typename Get_functor<Base, Point_weight_tag>::type Point_weight_d;
typedef typename Get_functor<Base, Point_drop_weight_tag>::type Point_drop_weight_d;
//TODO:
//typedef ??? Intersect_d;
@ -161,6 +168,7 @@ template <class Base_> struct Kernel_d_interface : public Base_ {
Point_dimension_d point_dimension_d_object()const{ return Point_dimension_d(*this); }
Point_of_sphere_d point_of_sphere_d_object()const{ return Point_of_sphere_d(*this); }
Side_of_oriented_sphere_d side_of_oriented_sphere_d_object()const{ return Side_of_oriented_sphere_d(*this); }
Power_test_d power_test_d_object()const{ return Power_test_d(*this); }
Side_of_bounded_sphere_d side_of_bounded_sphere_d_object()const{ return Side_of_bounded_sphere_d(*this); }
Contained_in_affine_hull_d contained_in_affine_hull_d_object()const{ return Contained_in_affine_hull_d(*this); }
Contained_in_linear_hull_d contained_in_linear_hull_d_object()const{ return Contained_in_linear_hull_d(*this); }
@ -168,6 +176,7 @@ template <class Base_> struct Kernel_d_interface : public Base_ {
Construct_flat_orientation_d construct_flat_orientation_d_object()const{ return Construct_flat_orientation_d(*this); }
In_flat_orientation_d in_flat_orientation_d_object()const{ return In_flat_orientation_d(*this); }
In_flat_side_of_oriented_sphere_d in_flat_side_of_oriented_sphere_d_object()const{ return In_flat_side_of_oriented_sphere_d(*this); }
In_flat_power_test_d in_flat_power_test_d_object()const{ return In_flat_power_test_d(*this); }
Point_to_vector_d point_to_vector_d_object()const{ return Point_to_vector_d(*this); }
Vector_to_point_d vector_to_point_d_object()const{ return Vector_to_point_d(*this); }
Affine_rank_d affine_rank_d_object()const{ return Affine_rank_d(*this); }
@ -193,6 +202,10 @@ template <class Base_> struct Kernel_d_interface : public Base_ {
Construct_ray_d construct_ray_d_object()const{ return Construct_ray_d(*this); }
Construct_iso_box_d construct_iso_box_d_object()const{ return Construct_iso_box_d(*this); }
Construct_aff_transformation_d construct_aff_transformation_d_object()const{ return Construct_aff_transformation_d(*this); }
Construct_weighted_point_d construct_weighted_point_d_object()const{ return Construct_weighted_point_d(*this); }
Point_weight_d point_weight_d_object()const{ return Point_weight_d(*this); }
Point_drop_weight_d point_drop_weight_d_object()const{ return Point_drop_weight_d(*this); }
// Dummies for those required functors missing a concept.
typedef Null_functor Position_on_line_d;

View File

@ -117,5 +117,17 @@ template <class K1, class K2> struct KO_converter<Sphere_tag,K1,K2>{
}
};
template <class K1, class K2> struct KO_converter<Weighted_point_tag,K1,K2>{
typedef typename Get_type<K1, Weighted_point_tag>::type argument_type;
typedef typename Get_type<K2, Weighted_point_tag>::type result_type;
template <class C>
result_type operator()(K1 const& k1, K2 const& k2, C const& conv, argument_type const& s) const {
typename Get_functor<K1, Point_drop_weight_tag>::type pdw(k1);
typename Get_functor<K1, Point_weight_tag>::type pw(k1);
typename Get_functor<K2, Construct_ttag<Weighted_point_tag> >::type cwp(k2);
return cwp(conv(pdw(s)),conv(pw(s)));
}
};
}
#endif

View File

@ -0,0 +1,125 @@
// Copyright (c) 2014
// INRIA Saclay-Ile de France (France)
//
// This file is part of CGAL (www.cgal.org); you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public License as
// published by the Free Software Foundation; either version 3 of the License,
// or (at your option) any later version.
//
// 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) : Marc Glisse
#ifndef CGAL_KD_TYPE_WP_H
#define CGAL_KD_TYPE_WP_H
#include <CGAL/NewKernel_d/store_kernel.h>
#include <boost/iterator/counting_iterator.hpp>
namespace CGAL {
namespace KerD {
template <class R_> class Weighted_point {
typedef typename Get_type<R_, FT_tag>::type FT_;
typedef typename Get_type<R_, Point_tag>::type Point_;
Point_ c_;
FT_ w_;
public:
Weighted_point(Point_ const&p, FT_ const&w): c_(p), w_(w) {}
// TODO: Add a piecewise constructor?
Point_ const& point()const{return c_;}
FT_ const& weight()const{return w_;}
};
}
namespace CartesianDKernelFunctors {
template <class R_> struct Construct_weighted_point : Store_kernel<R_> {
CGAL_FUNCTOR_INIT_STORE(Construct_weighted_point)
typedef typename Get_type<R_, Weighted_point_tag>::type result_type;
typedef typename Get_type<R_, Point_tag>::type Point;
typedef typename Get_type<R_, FT_tag>::type FT;
result_type operator()(Point const&a, FT const&b)const{
return result_type(a,b);
}
// Not really needed
result_type operator()()const{
typename Get_functor<R_, Construct_ttag<Point_tag> >::type cp(this->kernel());
return result_type(cp(),0);
}
};
template <class R_> struct Point_drop_weight {
CGAL_FUNCTOR_INIT_IGNORE(Point_drop_weight)
typedef typename Get_type<R_, Weighted_point_tag>::type argument_type;
typedef typename Get_type<R_, Point_tag>::type result_type;
// Returning a reference would be too fragile
result_type operator()(argument_type const&s)const{
return s.point();
}
};
template <class R_> struct Point_weight {
CGAL_FUNCTOR_INIT_IGNORE(Point_weight)
typedef typename Get_type<R_, Weighted_point_tag>::type argument_type;
typedef typename Get_type<R_, FT_tag>::type result_type;
result_type operator()(argument_type const&s)const{
return s.weight();
}
};
template<class R_> struct Power_test : private Store_kernel<R_> {
CGAL_FUNCTOR_INIT_STORE(Power_test)
typedef R_ R;
typedef typename Get_type<R, Oriented_side_tag>::type result_type;
template<class Iter, class Pt>
result_type operator()(Iter const& f, Iter const& e, Pt const& p0) const {
typename Get_functor<R, Power_test_raw_tag>::type ptr(this->kernel());
typename Get_functor<R, Point_drop_weight_tag>::type pdw(this->kernel());
typename Get_functor<R, Point_weight_tag>::type pw(this->kernel());
return ptr (
make_transforming_iterator (f, pdw),
make_transforming_iterator (e, pdw),
make_transforming_iterator (f, pw),
pdw (p0),
pw (p0));
}
};
template<class R_> struct In_flat_power_test : private Store_kernel<R_> {
CGAL_FUNCTOR_INIT_STORE(In_flat_power_test)
typedef R_ R;
typedef typename Get_type<R, Oriented_side_tag>::type result_type;
template<class Fo, class Iter, class Pt>
result_type operator()(Fo const& fo, Iter const& f, Iter const& e, Pt const& p0) const {
typename Get_functor<R, In_flat_power_test_raw_tag>::type ptr(this->kernel());
typename Get_functor<R, Point_drop_weight_tag>::type pdw(this->kernel());
typename Get_functor<R, Point_weight_tag>::type pw(this->kernel());
return ptr (
fo,
make_transforming_iterator (f, pdw),
make_transforming_iterator (e, pdw),
make_transforming_iterator (f, pw),
pdw (p0),
pw (p0));
}
};
}
CGAL_KD_DEFAULT_TYPE(Weighted_point_tag,(CGAL::KerD::Weighted_point<K>),(Point_tag),());
CGAL_KD_DEFAULT_FUNCTOR(Construct_ttag<Weighted_point_tag>,(CartesianDKernelFunctors::Construct_weighted_point<K>),(Weighted_point_tag,Point_tag),());
CGAL_KD_DEFAULT_FUNCTOR(Point_drop_weight_tag,(CartesianDKernelFunctors::Point_drop_weight<K>),(Weighted_point_tag,Point_tag),());
CGAL_KD_DEFAULT_FUNCTOR(Point_weight_tag,(CartesianDKernelFunctors::Point_weight<K>),(Weighted_point_tag,Point_tag),());
CGAL_KD_DEFAULT_FUNCTOR(Power_test_tag,(CartesianDKernelFunctors::Power_test<K>),(Weighted_point_tag),(Power_test_raw_tag,Point_drop_weight_tag,Point_weight_tag));
CGAL_KD_DEFAULT_FUNCTOR(In_flat_power_test_tag,(CartesianDKernelFunctors::In_flat_power_test<K>),(Weighted_point_tag),(In_flat_power_test_raw_tag,Point_drop_weight_tag,Point_weight_tag));
} // namespace CGAL
#endif

View File

@ -28,6 +28,7 @@
#include <CGAL/NewKernel_d/Wrapper/Segment_d.h>
#include <CGAL/NewKernel_d/Wrapper/Sphere_d.h>
#include <CGAL/NewKernel_d/Wrapper/Hyperplane_d.h>
#include <CGAL/NewKernel_d/Wrapper/Weighted_point_d.h>
#include <CGAL/NewKernel_d/Wrapper/Ref_count_obj.h>
@ -78,16 +79,17 @@ struct Forward_rep {
//};
//#else
template <class T,bool=Is_wrapper<T>::value,bool=Is_wrapper_iterator<T>::value> struct result_;
template <class T> struct result_<T,false,false>{typedef T const& type;};
template <class T> struct result_<T,true,false>{typedef typename decay<T>::type::Rep const& type;};
template <class T> struct result_<T,false,false>{typedef T type;}; // const&
template <class T> struct result_<T,true,false>{typedef typename decay<T>::type::Rep type;}; // const&
template <class T> struct result_<T,false,true>{typedef transforming_iterator<Forward_rep,typename decay<T>::type> type;};
template<class> struct result;
template<class T> struct result<Forward_rep(T)> : result_<T> {};
template <class T> typename boost::disable_if<boost::mpl::or_<Is_wrapper<T>,Is_wrapper_iterator<T> >,T>::type const& operator()(T const& t) const {return t;}
template <class T> typename boost::disable_if<boost::mpl::or_<Is_wrapper<T>,Is_wrapper_iterator<T> >,T>::type operator()(T const& t) const {return t;} // const&
template <class T> typename boost::disable_if<boost::mpl::or_<Is_wrapper<T>,Is_wrapper_iterator<T> >,T>::type& operator()(T& t) const {return t;}
template <class T> typename T::Rep const& operator()(T const& t, typename boost::enable_if<Is_wrapper<T> >::type* = 0) const {return t.rep();}
// FIXME: We should return const&, but it causes trouble inside a transform_iterator of an iterator that returns a prvalue :-(
template <class T> typename T::Rep operator()(T const& t, typename boost::enable_if<Is_wrapper<T> >::type* = 0) const {return t.rep();}
template <class T> transforming_iterator<Forward_rep,typename boost::enable_if<Is_wrapper_iterator<T>,T>::type> operator()(T const& t) const {return make_transforming_iterator(t,Forward_rep());}
//#endif
@ -106,6 +108,7 @@ CGAL_REGISTER_OBJECT_WRAPPER(Vector);
CGAL_REGISTER_OBJECT_WRAPPER(Segment);
CGAL_REGISTER_OBJECT_WRAPPER(Sphere);
CGAL_REGISTER_OBJECT_WRAPPER(Hyperplane);
CGAL_REGISTER_OBJECT_WRAPPER(Weighted_point);
#undef CGAL_REGISTER_OBJECT_WRAPPER
// Note: this tends to be an all or nothing thing currently, wrapping

View File

@ -0,0 +1,129 @@
// Copyright (c) 2014
// INRIA Saclay-Ile de France (France)
//
// This file is part of CGAL (www.cgal.org); you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public License as
// published by the Free Software Foundation; either version 3 of the License,
// or (at your option) any later version.
//
// 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) : Marc Glisse
#ifndef CGAL_WRAPPER_WEIGHTED_POINT_D_H
#define CGAL_WRAPPER_WEIGHTED_POINT_D_H
#include <CGAL/representation_tags.h>
#include <boost/static_assert.hpp>
#include <boost/type_traits.hpp>
#include <CGAL/Kernel/Return_base_tag.h>
#include <CGAL/Dimension.h>
#ifndef CGAL_CXX11
#include <boost/preprocessor/repetition.hpp>
#endif
#include <boost/utility/result_of.hpp>
namespace CGAL {
namespace Wrap {
template <class R_>
class Weighted_point_d : public Get_type<typename R_::Kernel_base, Weighted_point_tag>::type
{
typedef typename Get_type<R_, FT_tag>::type FT_;
typedef typename R_::Kernel_base Kbase;
typedef typename Get_type<R_, Point_tag>::type Point_;
typedef typename Get_functor<Kbase, Construct_ttag<Weighted_point_tag> >::type CWPBase;
typedef typename Get_functor<Kbase, Point_drop_weight_tag>::type PDWBase;
typedef typename Get_functor<Kbase, Point_weight_tag>::type PWBase;
typedef Weighted_point_d Self;
BOOST_STATIC_ASSERT((boost::is_same<Self, typename Get_type<R_, Weighted_point_tag>::type>::value));
public:
typedef Tag_true Is_wrapper;
typedef typename R_::Default_ambient_dimension Ambient_dimension;
typedef Dimension_tag<0> Feature_dimension;
typedef typename Get_type<Kbase, Weighted_point_tag>::type Rep;
const Rep& rep() const
{
return *this;
}
Rep& rep()
{
return *this;
}
typedef R_ R;
#ifdef CGAL_CXX11
template<class...U,class=typename std::enable_if<!std::is_same<std::tuple<typename std::decay<U>::type...>,std::tuple<Weighted_point_d> >::value>::type> explicit Weighted_point_d(U&&...u)
: Rep(CWPBase()(std::forward<U>(u)...)){}
// // called from Construct_point_d
// template<class...U> explicit Point_d(Eval_functor&&,U&&...u)
// : Rep(Eval_functor(), std::forward<U>(u)...){}
template<class F,class...U> explicit Weighted_point_d(Eval_functor&&,F&&f,U&&...u)
: Rep(std::forward<F>(f)(std::forward<U>(u)...)){}
#if 0
// the new standard may make this necessary
Point_d(Point_d const&)=default;
Point_d(Point_d &);//=default;
Point_d(Point_d &&)=default;
#endif
// try not to use these
Weighted_point_d(Rep const& v) : Rep(v) {}
Weighted_point_d(Rep& v) : Rep(static_cast<Rep const&>(v)) {}
Weighted_point_d(Rep&& v) : Rep(std::move(v)) {}
#else
Weighted_point_d() : Rep(CWPBase()()) {}
Weighted_point_d(Rep const& v) : Rep(v) {} // try not to use it
#define CGAL_CODE(Z,N,_) template<BOOST_PP_ENUM_PARAMS(N,class T)> \
explicit Weighted_point_d(BOOST_PP_ENUM_BINARY_PARAMS(N,T,const&t)) \
: Rep(CWPBase()( \
BOOST_PP_ENUM_PARAMS(N,t))) {} \
\
template<class F,BOOST_PP_ENUM_PARAMS(N,class T)> \
Weighted_point_d(Eval_functor,F const& f,BOOST_PP_ENUM_BINARY_PARAMS(N,T,const&t)) \
: Rep(f(BOOST_PP_ENUM_PARAMS(N,t))) {}
/*
template<BOOST_PP_ENUM_PARAMS(N,class T)> \
Point_d(Eval_functor,BOOST_PP_ENUM_BINARY_PARAMS(N,T,const&t)) \
: Rep(Eval_functor(), BOOST_PP_ENUM_PARAMS(N,t)) {}
*/
BOOST_PP_REPEAT_FROM_TO(1,11,CGAL_CODE,_)
#undef CGAL_CODE
#endif
//TODO: use references?
Point_ point()const{
return Point_(Eval_functor(),PDWBase(),rep());
}
FT_ weight()const{
return PWBase()(rep());
}
};
} //namespace Wrap
} //namespace CGAL
#endif // CGAL_WRAPPER_SPHERE_D_H

View File

@ -529,6 +529,60 @@ template<class R_> struct Orientation<R_,false> : private Store_kernel<R_> {
}
#endif
namespace CartesianDKernelFunctors {
template<class R_> struct Power_test_raw : private Store_kernel<R_> {
CGAL_FUNCTOR_INIT_STORE(Power_test_raw)
typedef R_ R;
typedef typename Get_type<R, RT_tag>::type RT;
typedef typename Get_type<R, FT_tag>::type FT;
typedef typename Get_type<R, Point_tag>::type Point;
typedef typename Get_type<R, Oriented_side_tag>::type result_type;
typedef typename Increment_dimension<typename R::Default_ambient_dimension>::type D1;
typedef typename Increment_dimension<typename R::Max_ambient_dimension>::type D2;
typedef typename R::LA::template Rebind_dimension<D1,D2>::Other LA;
typedef typename LA::Square_matrix Matrix;
template<class IterP, class IterW, class Pt, class Wt>
result_type operator()(IterP f, IterP const& e, IterW fw, Pt const& p0, Wt const& w0) const {
typedef typename Get_functor<R, Squared_distance_to_origin_tag>::type Sqdo;
typename Get_functor<R, Compute_point_cartesian_coordinate_tag>::type c(this->kernel());
typename Get_functor<R, Point_dimension_tag>::type pd(this->kernel());
int d=pd(p0);
Matrix m(d+1,d+1);
if(CGAL::Is_stored<Sqdo>::value) {
Sqdo sqdo(this->kernel());
FT const& h0 = sqdo(p0) - w0;
for(int i=0;f!=e;++f,++fw,++i) {
Point const& p=*f;
for(int j=0;j<d;++j){
RT const& x=c(p,j);
m(i,j)=x-c(p0,j);
}
m(i,d) = sqdo(p) - *fw - h0;
}
} else {
for(int i=0;f!=e;++f,++fw,++i) {
Point const& p=*f;
m(i,d) = w0 - *fw;
for(int j=0;j<d;++j){
RT const& x=c(p,j);
m(i,j)=x-c(p0,j);
m(i,d)+=CGAL::square(m(i,j));
}
}
}
if(d%2)
return -LA::sign_of_determinant(CGAL_MOVE(m));
else
return LA::sign_of_determinant(CGAL_MOVE(m));
}
};
}
CGAL_KD_DEFAULT_FUNCTOR(Power_test_raw_tag,(CartesianDKernelFunctors::Power_test_raw<K>),(Point_tag),(Point_dimension_tag,Squared_distance_to_origin_tag,Compute_point_cartesian_coordinate_tag));
// TODO: make Side_of_oriented_sphere call Power_test_raw
namespace CartesianDKernelFunctors {
template<class R_> struct Side_of_oriented_sphere : private Store_kernel<R_> {
CGAL_FUNCTOR_INIT_STORE(Side_of_oriented_sphere)

View File

@ -175,6 +175,7 @@ namespace CGAL {
CGAL_DECL_OBJ(Iso_box, Object);
CGAL_DECL_OBJ(Bbox, Object);
CGAL_DECL_OBJ(Aff_transformation, Object);
CGAL_DECL_OBJ(Weighted_point, Object);
#undef CGAL_DECL_OBJ_
#undef CGAL_DECL_OBJ
@ -216,6 +217,7 @@ namespace CGAL {
CGAL_DECL_COMPUTE(Scalar_product);
CGAL_DECL_COMPUTE(Hyperplane_translation);
CGAL_DECL_COMPUTE(Value_at);
CGAL_DECL_COMPUTE(Point_weight);
#undef CGAL_DECL_COMPUTE
#define CGAL_DECL_ITER_OBJ(X,Y,Z,C) struct X##_tag {}; \
@ -265,6 +267,7 @@ namespace CGAL {
CGAL_DECL_CONSTRUCT(Translated_point,Point);
CGAL_DECL_CONSTRUCT(Point_to_vector,Vector);
CGAL_DECL_CONSTRUCT(Vector_to_point,Point);
CGAL_DECL_CONSTRUCT(Point_drop_weight,Point);
#undef CGAL_DECL_CONSTRUCT
#if 0
#define CGAL_DECL_ITER_CONSTRUCT(X,Y) struct X##_tag {}; \
@ -304,6 +307,10 @@ namespace CGAL {
CGAL_DECL_PREDICATE(Affinely_independent);
CGAL_DECL_PREDICATE(Contained_in_linear_hull);
CGAL_DECL_PREDICATE(Contained_in_simplex);
CGAL_DECL_PREDICATE(Power_test_raw);
CGAL_DECL_PREDICATE(Power_test);
CGAL_DECL_PREDICATE(In_flat_power_test_raw);
CGAL_DECL_PREDICATE(In_flat_power_test);
#undef CGAL_DECL_PREDICATE
#define CGAL_DECL_MISC(X) struct X##_tag {}; \

View File

@ -26,17 +26,21 @@
namespace CGAL {
namespace internal {
BOOST_MPL_HAS_XXX_TRAIT_DEF(Do_not_store_kernel)
template<class T,bool=has_Do_not_store_kernel<T>::value> struct Do_not_store_kernel {
template<class T,bool=boost::is_empty<T>::value,bool=has_Do_not_store_kernel<T>::value> struct Do_not_store_kernel {
enum { value=false };
typedef Tag_false type;
};
template<class T> struct Do_not_store_kernel<T,true> {
template<class T> struct Do_not_store_kernel<T,true,false> {
enum { value=true };
typedef Tag_true type;
};
template<class T,bool b> struct Do_not_store_kernel<T,b,true> {
typedef typename T::Do_not_store_kernel type;
enum { value=type::value };
};
}
template<class R_,bool=boost::is_empty<R_>::value||internal::Do_not_store_kernel<R_>::value>
template<class R_,bool=internal::Do_not_store_kernel<R_>::value>
struct Store_kernel {
Store_kernel(){}
Store_kernel(R_ const&){}
@ -63,7 +67,7 @@ struct Store_kernel<R_,false> {
};
//For a second kernel. TODO: find something more elegant
template<class R_,bool=boost::is_empty<R_>::value||internal::Do_not_store_kernel<R_>::value>
template<class R_,bool=internal::Do_not_store_kernel<R_>::value>
struct Store_kernel2 {
Store_kernel2(){}
Store_kernel2(R_ const&){}

View File

@ -12,6 +12,7 @@
#include <CGAL/Interval_nt.h>
#include <CGAL/use.h>
#include <iostream>
#include <CGAL/NewKernel_d/Types/Weighted_point.h>
//typedef CGAL::Cartesian_base_d<double,CGAL::Dimension_tag<2> > K0;
//typedef CGAL::Cartesian_base_d<CGAL::Interval_nt_advanced,CGAL::Dimension_tag<2> > KA;
@ -466,6 +467,7 @@ void test3(){
P x4=cp(0,0,1);
P x5=cp(0,0,0);
P x6=cp(0,0,-1);
assert(!ed(x1,x2));
P tab2[]={x1,x2,x3,x4,x5};
assert(cis(tab2+0,tab2+4,x5));
assert(po(tab2+0,tab2+4)==CGAL::POSITIVE);
@ -513,6 +515,26 @@ void test3(){
assert(ifsos(fozn, tz+0, tz+3, tz[4]) == CGAL::ON_NEGATIVE_SIDE);
assert(ifsos(fozp, tz+0, tz+3, tz[5]) == CGAL::ON_NEGATIVE_SIDE);
assert(ifsos(fozn, tz+0, tz+3, tz[5]) == CGAL::ON_POSITIVE_SIDE);
typedef typename K1::Weighted_point_d WP;
typedef typename K1::Construct_weighted_point_d CWP;
typedef typename K1::Point_drop_weight_d PDW;
typedef typename K1::Point_weight_d PW;
typedef typename K1::Power_test_d PT;
typedef typename K1::In_flat_power_test_d IFPT;
CWP cwp Kinit(construct_weighted_point_d_object);
PDW pdw Kinit(point_drop_weight_d_object);
PW pw Kinit(point_weight_d_object);
PT pt Kinit(power_test_d_object);
IFPT ifpt Kinit(in_flat_power_test_d_object);
WP wp;
wp = cwp (x1, 2);
WP xw6 = cwp (x6, 0);
assert (pw(wp) == 2);
assert (ed(pdw(wp), x1));
WP tabw[]={cwp(x1,0),cwp(x2,0),cwp(x3,0),cwp(x4,0),cwp(x5,0)};
assert(pt(tabw+0,tabw+4,tabw[4])==CGAL::ON_POSITIVE_SIDE);
assert(ifpt(fo4,tabw+0,tabw+3,xw6)==CGAL::ON_POSITIVE_SIDE);
}
template struct CGAL::Epick_d<CGAL::Dimension_tag<2> >;
template struct CGAL::Epick_d<CGAL::Dimension_tag<3> >;

View File

@ -1,15 +1,17 @@
all: normal cxx11
CXX = g++
CGAL_INC = -I. -I../../include -I../../../Algebraic_foundations/include -I../../../STL_Extension/include -I../../../Number_types/include -I../../../Kernel_23/include -I../../../Installation/include -DCGAL_DISABLE_ROUNDING_MATH_CHECK
EIGEN_INC = `pkg-config --cflags eigen3|sed -e 's/-I/-isystem/g'` -DCGAL_EIGEN3_ENABLED
normal:
g++ Epick_d.cpp -O2 -lCGAL -lboost_thread -frounding-math -Wall -Wextra -lmpfr -lgmp ${CGAL_INC} ${EIGEN_INC}
${CXX} Epick_d.cpp -O2 -lCGAL -lboost_thread -frounding-math -Wall -Wextra -lmpfr -lgmp ${CGAL_INC} ${EIGEN_INC}
./a.out
cxx11:
g++ -std=c++0x -O2 Epick_d.cpp -lCGAL -lboost_thread -frounding-math -Wall -Wextra -lmpfr -lgmp ${CGAL_INC} ${EIGEN_INC} -o b.out
${CXX} -std=c++0x -O2 Epick_d.cpp -lCGAL -lboost_thread -frounding-math -Wall -Wextra -lmpfr -lgmp ${CGAL_INC} ${EIGEN_INC} -o b.out
./b.out
#-DBOOST_RESULT_OF_USE_DECLTYPE

View File

@ -155,6 +155,17 @@ namespace Eigen {
MulCost = 100
};
};
namespace internal {
template<>
struct significant_decimals_impl<CGAL::Gmpq>
{
static inline int run()
{
return 0;
}
};
}
}
//since types are included by Gmp_coercion_traits.h:

View File

@ -1283,6 +1283,13 @@ namespace Eigen {
MulCost = 10
};
};
namespace internal {
template<class> struct significant_decimals_impl;
template<bool b>
struct significant_decimals_impl<CGAL::Interval_nt<b> >
: significant_decimals_impl<typename CGAL::Interval_nt<b>::value_type> { };
}
}
#endif // CGAL_INTERVAL_NT_H

View File

@ -0,0 +1,11 @@
2
0.0071 1.6899 0
0.3272 1.3694 0.05
1.3697 1.8296 0.1
0.6722 0.3012 0.15
1.1726 0.1899 0.2
0.4374 2.8541 0.25
2.5923 0.1904 0.3
1.3083 2.5462 0.35
1.4981 1.3929 0.4
2.1304 2.055 0.45

View File

@ -0,0 +1,11 @@
3
0.0071 1.6899 2.521 0
0.3272 1.3694 3.15 0.05
1.3697 1.8296 2.654 0.1
0.6722 0.3012 0.1548 0.15
1.1726 0.1899 0.3658 0.2
0.4374 2.8541 1.45894 0.25
2.5923 0.1904 0.6971 0.3
1.3083 2.5462 1.3658 0.35
1.4981 1.3929 2.949 0.4
2.1304 2.055 0.6597455 0.45

View File

@ -0,0 +1,43 @@
#include <CGAL/Epick_d.h>
#include <CGAL/Regular_triangulation_euclidean_traits.h>
#include <CGAL/Regular_triangulation.h>
#include <CGAL/IO/Triangulation_off_ostream.h>
#include <fstream>
typedef CGAL::Epick_d<CGAL::Dynamic_dimension_tag> K;
typedef CGAL::Regular_triangulation_euclidean_traits<K> Traits;
typedef CGAL::Regular_triangulation<Traits> RT;
void test(int dim)
{
std::stringstream input_filename;
input_filename << "data/points_" << dim << ".cin";
std::ifstream in(input_filename.str());
RT::Weighted_point wp;
std::vector<RT::Weighted_point> wpoints;
int dim_from_file;
in >> dim_from_file;
while(in >> wp)
wpoints.push_back(wp);
// Build the Regular Triangulation
RT rt(dim_from_file);
rt.insert(wpoints.begin(), wpoints.end());
CGAL_assertion(rt.is_valid(true));
// Export
std::stringstream output_filename;
output_filename << "data/rt_dim" << dim << ".off";
std::ofstream off_stream(output_filename.str());
CGAL::export_triangulation_to_off(off_stream, rt);
}
int main()
{
test(2);
test(3);
return 0;
}

View File

@ -65,6 +65,5 @@ int main(int argc, char **argv)
go<7>(N);
go<8>(N);
return 0;
}

View File

@ -1,3 +1,2 @@
TODO
include/CGAL/Convex_hull.h
include/CGAL/Regular_triangulation.h

View File

@ -55,7 +55,7 @@ class Delaunay_triangulation
public: // PUBLIC NESTED TYPES
typedef DCTraits Geom_traits;
typedef typename Base::Triangulation_ds Triangulation_ds;
typedef typename Base::Triangulation_ds Triangulation_ds;
typedef typename Base::Vertex Vertex;
typedef typename Base::Full_cell Full_cell;
@ -71,21 +71,25 @@ public: // PUBLIC NESTED TYPES
typedef typename Base::Vertex_const_handle Vertex_const_handle;
typedef typename Base::Vertex_const_iterator Vertex_const_iterator;
typedef typename Base::Full_cell_handle Full_cell_handle;
typedef typename Base::Full_cell_iterator Full_cell_iterator;
typedef typename Base::Full_cell_const_handle Full_cell_const_handle;
typedef typename Base::Full_cell_const_iterator Full_cell_const_iterator;
typedef typename Base::Full_cell_handle Full_cell_handle;
typedef typename Base::Full_cell_iterator Full_cell_iterator;
typedef typename Base::Full_cell_const_handle Full_cell_const_handle;
typedef typename Base::Full_cell_const_iterator Full_cell_const_iterator;
typedef typename Base::size_type size_type;
typedef typename Base::difference_type difference_type;
typedef typename Base::Locate_type Locate_type;
//Tag to distinguish triangulations with weighted_points
typedef Tag_false Weighted_tag;
protected: // DATA MEMBERS
public:
using typename Base::Rotor;
using Base::maximal_dimension;
using Base::are_incident_full_cells_valid;
using Base::coaffine_orientation_predicate;
@ -95,11 +99,12 @@ public:
//using Base::incident_full_cells;
using Base::geom_traits;
using Base::index_of_covertex;
using Base::index_of_second_covertex;
using Base::rotate_rotor;
using Base::infinite_vertex;
using Base::insert_in_hole;
using Base::insert_outside_convex_hull_1;
using Base::is_infinite;
using Base::is_valid;
using Base::locate;
using Base::points_begin;
using Base::set_neighbors;
@ -142,32 +147,6 @@ private:
}
};
public:
// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - UTILITIES
// A co-dimension 2 sub-simplex. called a Rotor because we can rotate
// the two "covertices" around the sub-simplex. Useful for traversing the
// boundary of a hole. NOT DOCUMENTED
typedef cpp11::tuple<Full_cell_handle, int, int> Rotor;
Full_cell_handle full_cell(const Rotor & r) const // NOT DOCUMENTED
{
return cpp11::get<0>(r);
}
int index_of_covertex(const Rotor & r) const // NOT DOCUMENTED
{
return cpp11::get<1>(r);
}
int index_of_second_covertex(const Rotor & r) const // NOT DOCUMENTED
{
return cpp11::get<2>(r);
}
Rotor rotate_rotor(Rotor & r) // NOT DOCUMENTED...
{
int opposite = full_cell(r)->mirror_index(index_of_covertex(r));
Full_cell_handle s = full_cell(r)->neighbor(index_of_covertex(r));
int new_second = s->index(full_cell(r)->vertex(index_of_second_covertex(r)));
return Rotor(s, new_second, opposite);
}
// - - - - - - - - - - - - - - - - - - - - - - - - - - CREATION / CONSTRUCTORS
@ -339,6 +318,10 @@ public:
return pred_(dc_.full_cell(f)->neighbor(dc_.index_of_covertex(f)));
}
};
// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - VALIDITY
bool is_valid(bool verbose = false, int level = 0) const;
private:
// Some internal types to shorten notation
@ -352,27 +335,6 @@ private:
Conflict_traversal_pred_in_subspace;
typedef Conflict_traversal_predicate<Conflict_pred_in_fullspace>
Conflict_traversal_pred_in_fullspace;
// This is used in the |remove(v)| member function to manage sets of Full_cell_handles
template< typename FCH >
struct Full_cell_set : public std::vector<FCH>
{
typedef std::vector<FCH> Base_set;
using Base_set::begin;
using Base_set::end;
void make_searchable()
{ // sort the full cell handles
std::sort(begin(), end());
}
bool contains(const FCH & fch) const
{
return std::binary_search(begin(), end(), fch);
}
bool contains_1st_and_not_2nd(const FCH & fst, const FCH & snd) const
{
return ( ! contains(snd) ) && ( contains(fst) );
}
};
};
// = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
@ -427,7 +389,7 @@ Delaunay_triangulation<DCTraits, TDS>
// THE CASE cur_dim >= 2
// Gather the finite vertices sharing an edge with |v|
typedef Full_cell_set<Full_cell_handle> Simplices;
typedef typename Base::template Full_cell_set<Full_cell_handle> Simplices;
Simplices simps;
std::back_insert_iterator<Simplices> out(simps);
tds().incident_full_cells(v, out);
@ -563,7 +525,7 @@ Delaunay_triangulation<DCTraits, TDS>
Dark_s_handle dark_ret_s = dark_s;
Full_cell_handle ret_s;
typedef Full_cell_set<Dark_s_handle> Dark_full_cells;
typedef typename Base::template Full_cell_set<Dark_s_handle> Dark_full_cells;
Dark_full_cells conflict_zone;
std::back_insert_iterator<Dark_full_cells> dark_out(conflict_zone);
@ -777,7 +739,6 @@ Delaunay_triangulation<DCTraits, TDS>
::insert_in_conflicting_cell(const Point & p, const Full_cell_handle s)
{
typedef std::vector<Full_cell_handle> Full_cell_h_vector;
typedef typename Full_cell_h_vector::iterator SHV_iterator;
static Full_cell_h_vector cs; // for storing conflicting full_cells.
cs.clear();
// cs.reserve(64);
@ -888,6 +849,48 @@ Delaunay_triangulation<DCTraits, TDS>
}
}
// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - VALIDITY
template< typename DCTraits, typename TDS >
bool
Delaunay_triangulation<DCTraits, TDS>
::is_valid(bool verbose, int level) const
{
if (!Base::is_valid(verbose, level))
return false;
int dim = current_dimension();
if (dim == maximal_dimension())
{
for (Finite_full_cell_const_iterator cit = finite_full_cells_begin() ;
cit != finite_full_cells_end() ; ++cit )
{
Full_cell_const_handle ch = cit.base();
for(int i = 0; i < dim+1 ; ++i )
{
// If the i-th neighbor is not an infinite cell
Vertex_handle opposite_vh =
ch->neighbor(i)->vertex(ch->neighbor(i)->index(ch));
if (!is_infinite(opposite_vh))
{
Side_of_oriented_sphere_d side =
geom_traits().side_of_oriented_sphere_d_object();
if (side(Point_const_iterator(ch->vertices_begin()),
Point_const_iterator(ch->vertices_end()),
opposite_vh->point()) == ON_BOUNDED_SIDE)
{
if (verbose)
CGAL_warning_msg(false, "Non-empty sphere");
return false;
}
}
}
}
}
return true;
}
} //namespace CGAL
#endif // CGAL_DELAUNAY_COMPLEX_H

View File

@ -0,0 +1,255 @@
// Copyright (c) 2014 INRIA Sophia-Antipolis (France).
// All rights reserved.
//
// This file is part of CGAL (www.cgal.org); you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public License as
// published by the Free Software Foundation; either version 3 of the License,
// or (at your option) any later version.
//
// 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) : Clement Jamin
#ifndef CGAL_TRIANGULATION_IO_H
#define CGAL_TRIANGULATION_IO_H
#include <CGAL/Epick_d.h>
#include <CGAL/Triangulation.h>
#include <sstream>
#include <iostream>
namespace CGAL {
namespace Triangulation_IO
{
// TODO: test if the stream is binary or text?
template<typename Traits, typename P>
void
output_point(std::ostream & os, const Traits &traits, const P & p)
{
typedef typename Traits::Compute_coordinate_d Ccd;
const Ccd ccd = traits.compute_coordinate_d_object();
const int dim = traits.point_dimension_d_object()(p);
if (dim > 0)
{
os << ccd(p, 0);
for (int i = 1 ; i < dim ; ++i)
os << " " << CGAL::to_double(ccd(p, i));
}
}
// TODO: test if the stream is binary or text?
/*template<typename Traits, typename P>
void
input_point(std::istream & is, const Traits &traits, P & p)
{
typedef typename Traits::FT FT;
std::vector<FT> coords;
std::string line;
for(;;)
{
if (!std::getline(is, line))
return is;
if (line != "")
break;
}
std::stringstream line_sstr(line);
FT temp;
while (line_sstr >> temp)
coords.push_back(temp);
p = traits.construct_point_d_object()(coords.begin(), coords.end());
}*/
} // namespace Triangulation_IO
///////////////////////////////////////////////////////////////
// TODO: replace these operator>> by an "input_point" function
///////////////////////////////////////////////////////////////
// TODO: test if the stream is binary or text?
template<typename K>
std::istream &
operator>>(std::istream &is, typename Wrap::Point_d<K> & p)
{
typedef typename Wrap::Point_d<K> P;
typedef typename K::FT FT;
std::vector<FT> coords;
std::string line;
for(;;)
{
if (!std::getline(is, line))
return is;
if (line != "")
break;
}
std::stringstream line_sstr(line);
FT temp;
while (line_sstr >> temp)
coords.push_back(temp);
p = P(coords.begin(), coords.end());
return is;
}
// TODO: test if the stream is binary or text?
template<typename K>
std::istream &
operator>>(std::istream &is, typename Wrap::Weighted_point_d<K> & wp)
{
typedef typename Wrap::Point_d<K> P;
typedef typename Wrap::Weighted_point_d<K> WP;
typedef typename K::FT FT;
std::string line;
for(;;)
{
if (!std::getline(is, line))
return is;
if (line != "")
break;
}
std::stringstream line_sstr(line);
FT temp;
std::vector<FT> coords;
while (line_sstr >> temp)
coords.push_back(temp);
std::vector<FT>::iterator last = coords.end() - 1;
P p = P(coords.begin(), last);
wp = WP(p, *last);
return is;
}
template < class GT, class TDS >
std::ostream &
export_triangulation_to_off(std::ostream & os,
const Triangulation<GT,TDS> & tr,
bool in_3D_export_surface_only = false)
{
typedef Triangulation<GT,TDS> Tr;
typedef typename Tr::Vertex_const_handle Vertex_handle;
typedef typename Tr::Vertex_const_iterator Vertex_iterator;
typedef typename Tr::Finite_vertex_const_iterator Finite_vertex_iterator;
typedef typename Tr::Full_cell_const_handle Full_cell_handle;
typedef typename Tr::Finite_full_cell_const_iterator Finite_full_cell_iterator;
typedef typename Tr::Full_cell_const_iterator Full_cell_iterator;
typedef typename Tr::Full_cell Full_cell;
typedef typename Full_cell::Vertex_handle_const_iterator Full_cell_vertex_iterator;
if (tr.maximal_dimension() < 2 || tr.maximal_dimension() > 3)
{
std::cerr << "Warning: export_tds_to_off => dimension should be 2 or 3.";
os << "Warning: export_tds_to_off => dimension should be 2 or 3.";
return os;
}
size_t n = tr.number_of_vertices();
std::stringstream output;
// write the vertices
std::map<Vertex_handle, int> index_of_vertex;
int i = 0;
for(Finite_vertex_iterator it = tr.finite_vertices_begin();
it != tr.finite_vertices_end(); ++it, ++i)
{
Triangulation_IO::output_point(output, tr.geom_traits(), it->point());
if (tr.maximal_dimension() == 2)
output << " 0";
output << std::endl;
index_of_vertex[it.base()] = i;
}
CGAL_assertion( i == n );
size_t number_of_triangles = 0;
if (tr.maximal_dimension() == 2)
{
for (Finite_full_cell_iterator fch = tr.finite_full_cells_begin() ;
fch != tr.finite_full_cells_end() ; ++fch)
{
output << "3 ";
for (Full_cell_vertex_iterator vit = fch->vertices_begin() ;
vit != fch->vertices_end() ; ++vit)
{
output << index_of_vertex[*vit] << " ";
}
output << std::endl;
++number_of_triangles;
}
}
else if (tr.maximal_dimension() == 3)
{
if (in_3D_export_surface_only)
{
// Parse boundary facets
for (Full_cell_iterator fch = tr.full_cells_begin() ;
fch != tr.full_cells_end() ; ++fch)
{
if (tr.is_infinite(fch))
{
output << "3 ";
for (Full_cell_vertex_iterator vit = fch->vertices_begin() ;
vit != fch->vertices_end() ; ++vit)
{
if (!tr.is_infinite(*vit))
output << index_of_vertex[*vit] << " ";
}
output << std::endl;
++number_of_triangles;
}
}
}
else
{
// Parse finite cells
for (Finite_full_cell_iterator fch = tr.finite_full_cells_begin() ;
fch != tr.finite_full_cells_end() ; ++fch)
{
output << "3 "
<< index_of_vertex[fch->vertex(0)] << " "
<< index_of_vertex[fch->vertex(1)] << " "
<< index_of_vertex[fch->vertex(2)]
<< std::endl;
output << "3 "
<< index_of_vertex[fch->vertex(0)] << " "
<< index_of_vertex[fch->vertex(2)] << " "
<< index_of_vertex[fch->vertex(3)]
<< std::endl;
output << "3 "
<< index_of_vertex[fch->vertex(1)] << " "
<< index_of_vertex[fch->vertex(2)] << " "
<< index_of_vertex[fch->vertex(3)]
<< std::endl;
output << "3 "
<< index_of_vertex[fch->vertex(0)] << " "
<< index_of_vertex[fch->vertex(1)] << " "
<< index_of_vertex[fch->vertex(3)]
<< std::endl;
number_of_triangles += 4;
}
}
}
os << "OFF \n"
<< n << " "
<< number_of_triangles << " 0\n"
<< output.str();
return os;
}
} //namespace CGAL
#endif // CGAL_TRIANGULATION_IO_H

View File

@ -1,4 +1,4 @@
// Copyright (c) 2009-2014 INRIA Sophia-Antipolis (France).
// Copyright (c) 2014 INRIA Sophia-Antipolis (France).
// All rights reserved.
//
// This file is part of CGAL (www.cgal.org).
@ -15,7 +15,7 @@
// $URL$
// $Id$
//
// Author(s) : Samuel Hornus
// Author(s) : Clement Jamin
#ifndef CGAL_REGULAR_TRIANGULATION_H
#define CGAL_REGULAR_TRIANGULATION_H
@ -28,27 +28,909 @@ namespace CGAL {
template< typename RTTraits, typename TDS_ = Default >
class Regular_triangulation
: public Triangulation<RTTraits,
typename Default::Get<TDS_, Triangulation_data_structure<
typename Maximal_dimension<typename RTTraits::Point_d>::type,
Triangulation_vertex<RTTraits>,
Triangulation_full_cell<RTTraits> >
>::type >
: public Triangulation<
RTTraits,
typename Default::Get<TDS_, Triangulation_data_structure<
typename RTTraits::Dimension,
Triangulation_vertex<RTTraits>,
Triangulation_full_cell<RTTraits> >
>::type >
{
typedef typename Maximal_dimension<typename RTTraits::Point_d>::type
Maximal_dimension_;
typedef typename Default::Get<TDS_, Triangulation_data_structure<
Maximal_dimension_,
Triangulation_vertex<RTTraits>,
Triangulation_full_cell<RTTraits> >
>::type TDS;
typedef Triangulation<RTTraits, TDS> Base;
typedef Regular_triangulation<RTTraits, TDS_> Self;
typedef typename RTTraits::Dimension Maximal_dimension_;
typedef typename Default::Get<
TDS_,
Triangulation_data_structure<
Maximal_dimension_,
Triangulation_vertex<RTTraits>,
Triangulation_full_cell<RTTraits>
> >::type TDS;
typedef Triangulation<RTTraits, TDS> Base;
typedef Regular_triangulation<RTTraits, TDS_> Self;
typedef typename RTTraits::Orientation_d Orientation_d;
typedef typename RTTraits::Construct_weighted_point_d Construct_weighted_point_d;
typedef typename RTTraits::Power_test_d Power_test_d;
typedef typename RTTraits::In_flat_power_test_d In_flat_power_test_d;
typedef typename RTTraits::Flat_orientation_d Flat_orientation_d;
typedef typename RTTraits::Construct_flat_orientation_d Construct_flat_orientation_d;
typedef typename RTTraits::In_flat_orientation_d In_flat_orientation_d;
public: // PUBLIC NESTED TYPES
typedef RTTraits Geom_traits;
typedef typename Base::Triangulation_ds Triangulation_ds;
typedef typename Base::Vertex Vertex;
typedef typename Base::Full_cell Full_cell;
typedef typename Base::Facet Facet;
typedef typename Base::Face Face;
typedef Maximal_dimension_ Maximal_dimension;
typedef typename RTTraits::Bare_point Bare_point;
typedef typename RTTraits::Weighted_point Weighted_point;
typedef typename Base::Point_const_iterator Point_const_iterator;
typedef typename Base::Vertex_handle Vertex_handle;
typedef typename Base::Vertex_iterator Vertex_iterator;
typedef typename Base::Vertex_const_handle Vertex_const_handle;
typedef typename Base::Vertex_const_iterator Vertex_const_iterator;
typedef typename Base::Full_cell_handle Full_cell_handle;
typedef typename Base::Full_cell_iterator Full_cell_iterator;
typedef typename Base::Full_cell_const_handle Full_cell_const_handle;
typedef typename Base::Full_cell_const_iterator Full_cell_const_iterator;
typedef typename Base::Finite_full_cell_const_iterator
Finite_full_cell_const_iterator;
typedef typename Base::size_type size_type;
typedef typename Base::difference_type difference_type;
typedef typename Base::Locate_type Locate_type;
//Tag to distinguish Delaunay from Regular triangulations
typedef Tag_true Weighted_tag;
protected: // DATA MEMBERS
public:
typedef Maximal_dimension_ Maximal_dimension;
};
using typename Base::Rotor;
using Base::maximal_dimension;
using Base::are_incident_full_cells_valid;
using Base::coaffine_orientation_predicate;
using Base::reset_flat_orientation;
using Base::current_dimension;
using Base::geom_traits;
using Base::index_of_covertex;
using Base::index_of_second_covertex;
using Base::rotate_rotor;
using Base::infinite_vertex;
using Base::insert_in_hole;
using Base::insert_outside_convex_hull_1;
using Base::is_infinite;
using Base::locate;
using Base::points_begin;
using Base::set_neighbors;
using Base::new_full_cell;
using Base::number_of_vertices;
using Base::orientation;
using Base::tds;
using Base::reorient_full_cells;
using Base::full_cell;
using Base::full_cells_begin;
using Base::full_cells_end;
using Base::finite_full_cells_begin;
using Base::finite_full_cells_end;
using Base::vertices_begin;
using Base::vertices_end;
private:
//*** Power_test_in_flat_d *** CJTODO: better name?
// Wrapper
struct Power_test_in_flat_d
{
boost::optional<Flat_orientation_d>* fop;
Construct_flat_orientation_d cfo;
In_flat_power_test_d ifpt;
Power_test_in_flat_d(
boost::optional<Flat_orientation_d>& x,
Construct_flat_orientation_d const&y,
In_flat_power_test_d const&z)
: fop(&x), cfo(y), ifpt(z) {}
template<class Iter>
CGAL::Orientation operator()(Iter a, Iter b, const Weighted_point & p)const
{
if(!*fop)
*fop=cfo(a,b);
return ifpt(fop->get(),a,b,p);
}
};
public:
// - - - - - - - - - - - - - - - - - - - - - - - - - - CREATION / CONSTRUCTORS
Regular_triangulation(int dim, const Geom_traits k = Geom_traits())
: Base(dim, k)
{
}
// With this constructor,
// the user can specify a Flat_orientation_d object to be used for
// orienting simplices of a specific dimension
// (= preset_flat_orientation_.first)
// It it used by the dark triangulations created by DT::remove
Regular_triangulation(
int dim,
const std::pair<int, const Flat_orientation_d *> &preset_flat_orientation,
const Geom_traits k = Geom_traits())
: Base(dim, preset_flat_orientation, k)
{
}
~Regular_triangulation() {}
// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - ACCESS
// Not Documented
Power_test_in_flat_d power_test_in_flat_predicate() const
{
return Power_test_in_flat_d (
flat_orientation_,
geom_traits().construct_flat_orientation_d_object(),
geom_traits().in_flat_power_test_d_object()
);
}
// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - REMOVALS
Full_cell_handle remove(Vertex_handle);
Full_cell_handle remove(const Weighted_point & p, Full_cell_handle hint = Full_cell_handle())
{
Locate_type lt;
Face f(maximal_dimension());
Facet ft;
Full_cell_handle s = locate(p, lt, f, ft, hint);
if( Base::ON_VERTEX == lt )
{
return remove(s->vertex(f.index(0)));
}
return Full_cell_handle();
}
template< typename ForwardIterator >
void remove(ForwardIterator start, ForwardIterator end)
{
while( start != end )
remove(*start++);
}
// Not documented
void remove_decrease_dimension(Vertex_handle);
// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - INSERTIONS
template< typename ForwardIterator >
size_type insert(ForwardIterator start, ForwardIterator end)
{
size_type n = number_of_vertices();
typedef std::vector<Weighted_point> WP_vec;
WP_vec points(start, end);
typename Geom_traits::Point_drop_weight_d pdw =
geom_traits().point_drop_weight_d_object();
spatial_sort(
boost::make_transform_iterator(points.begin(), pdw),
boost::make_transform_iterator(points.end(), pdw),
typename Geom_traits::Base());
//spatial_sort(points.begin(), points.end(), geom_traits()); // CJTODO TEMP A REMETTRE
//spatial_sort(points.begin(), points.end(), Geom_traits::Base());
Full_cell_handle hint;
for(typename WP_vec::const_iterator p = points.begin(); p != points.end(); ++p )
{
Locate_type lt;
Face f(maximal_dimension());
Facet ft;
Full_cell_handle c = locate (*p, lt, f, ft, hint);
Vertex_handle v = insert (*p, lt, f, ft, c);
hint = v == Vertex_handle() ? c : v->full_cell();
}
return number_of_vertices() - n;
}
Vertex_handle insert(const Weighted_point &,
const Locate_type,
const Face &,
const Facet &,
const Full_cell_handle);
Vertex_handle insert(const Weighted_point & p,
const Full_cell_handle start = Full_cell_handle())
{
Locate_type lt;
Face f(maximal_dimension());
Facet ft;
Full_cell_handle s = locate(p, lt, f, ft, start);
return insert(p, lt, f, ft, s);
}
Vertex_handle insert(const Weighted_point & p, const Vertex_handle hint)
{
CGAL_assertion( Vertex_handle() != hint );
return insert(p, hint->full_cell());
}
Vertex_handle insert_outside_affine_hull(const Weighted_point &);
Vertex_handle insert_in_conflicting_cell(const Weighted_point &, const Full_cell_handle);
// - - - - - - - - - - - - - - - - - - - - - - - - - GATHERING CONFLICTING SIMPLICES
bool is_in_conflict(const Weighted_point &, Full_cell_const_handle) const;
template< class OrientationPredicate >
Oriented_side perturbed_power_test(const Weighted_point &,
Full_cell_const_handle, const OrientationPredicate &) const;
template< typename OutputIterator >
Facet compute_conflict_zone(const Weighted_point &, const Full_cell_handle, OutputIterator) const;
template < typename OrientationPredicate, typename PowerTestPredicate >
class Conflict_predicate
{
const Self & rt_;
const Weighted_point & p_;
OrientationPredicate ori_;
PowerTestPredicate power_test_;
int cur_dim_;
public:
Conflict_predicate(
const Self & rt,
const Weighted_point & p,
const OrientationPredicate & ori,
const PowerTestPredicate & power_test)
: rt_(rt), p_(p), ori_(ori), power_test_(power_test), cur_dim_(rt.current_dimension()) {}
inline
bool operator()(Full_cell_const_handle s) const
{
bool ok;
if( ! rt_.is_infinite(s) )
{
Oriented_side power_test = power_test_(rt_.points_begin(s), rt_.points_begin(s) + cur_dim_ + 1, p_);
if( ON_POSITIVE_SIDE == power_test )
ok = true;
else if( ON_NEGATIVE_SIDE == power_test )
ok = false;
else
ok = ON_POSITIVE_SIDE == rt_.perturbed_power_test<OrientationPredicate>(p_, s, ori_);
}
else
{
typedef typename Full_cell::Vertex_handle_const_iterator VHCI;
typedef Substitute_point_in_vertex_iterator<VHCI> F;
F spivi(rt_.infinite_vertex(), &p_);
Orientation o = ori_(
boost::make_transform_iterator(s->vertices_begin(), spivi),
boost::make_transform_iterator(s->vertices_begin() + cur_dim_ + 1,
spivi));
if( POSITIVE == o )
ok = true;
else if( o == NEGATIVE )
ok = false;
else
ok = (*this)(s->neighbor( s->index( rt_.infinite_vertex() ) ));
}
return ok;
}
};
template < typename ConflictPredicate >
class Conflict_traversal_predicate
{
const Self & rt_;
const ConflictPredicate & pred_;
public:
Conflict_traversal_predicate(const Self & rt, const ConflictPredicate & pred)
: rt_(rt), pred_(pred)
{}
inline
bool operator()(const Facet & f) const
{
return pred_(rt_.full_cell(f)->neighbor(rt_.index_of_covertex(f)));
}
};
// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - VALIDITY
bool is_valid(bool verbose = false, int level = 0) const;
private:
// Some internal types to shorten notation
using typename Base::Coaffine_orientation_d;
using Base::flat_orientation_;
typedef Conflict_predicate<Coaffine_orientation_d, Power_test_in_flat_d>
Conflict_pred_in_subspace;
typedef Conflict_predicate<Orientation_d, Power_test_d>
Conflict_pred_in_fullspace;
typedef Conflict_traversal_predicate<Conflict_pred_in_subspace>
Conflict_traversal_pred_in_subspace;
typedef Conflict_traversal_predicate<Conflict_pred_in_fullspace>
Conflict_traversal_pred_in_fullspace;
}; // class Regular_triangulation
// = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
// FUNCTIONS THAT ARE MEMBER METHODS:
// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - REMOVALS
template< typename RTTraits, typename TDS >
typename Regular_triangulation<RTTraits, TDS>::Full_cell_handle
Regular_triangulation<RTTraits, TDS>
::remove( Vertex_handle v )
{
CGAL_precondition( ! is_infinite(v) );
CGAL_expensive_precondition( is_vertex(v) );
// THE CASE cur_dim == 0
if( 0 == current_dimension() )
{
remove_decrease_dimension(v);
return Full_cell_handle();
}
else if( 1 == current_dimension() )
{ // THE CASE cur_dim == 1
if( 2 == number_of_vertices() )
{
remove_decrease_dimension(v);
return Full_cell_handle();
}
Full_cell_handle left = v->full_cell();
if( is_infinite(left) && left->neighbor(0)->index(left) == 0 ) // we are on the infinite right.
left = left->neighbor(0);
if( 0 == left->index(v) )
left = left->neighbor(1);
CGAL_assertion( 1 == left->index(v) );
Full_cell_handle right = left->neighbor(0);
if( ! is_infinite(right) )
{
tds().associate_vertex_with_full_cell(left, 1, right->vertex(1));
set_neighbors(left, 0, right->neighbor(0), right->mirror_index(0));
}
else
{
tds().associate_vertex_with_full_cell(left, 1, left->vertex(0));
tds().associate_vertex_with_full_cell(left, 0, infinite_vertex());
set_neighbors(left, 0, left->neighbor(1), left->mirror_index(1));
set_neighbors(left, 1, right->neighbor(1), right->mirror_index(1));
}
tds().delete_vertex(v);
tds().delete_full_cell(right);
return left;
}
// THE CASE cur_dim >= 2
// Gather the finite vertices sharing an edge with |v|
typedef typename Base::template Full_cell_set<Full_cell_handle> Simplices;
Simplices simps;
std::back_insert_iterator<Simplices> out(simps);
tds().incident_full_cells(v, out);
typedef std::set<Vertex_handle> Vertex_set;
Vertex_set verts;
Vertex_handle vh;
for( typename Simplices::iterator it = simps.begin(); it != simps.end(); ++it )
for( int i = 0; i <= current_dimension(); ++i )
{
vh = (*it)->vertex(i);
if( is_infinite(vh) )
continue;
if( vh == v )
continue;
verts.insert(vh);
}
// After gathering finite neighboring vertices, create their Dark Delaunay triangulation
typedef Triangulation_vertex<Geom_traits, Vertex_handle> Dark_vertex_base;
typedef Triangulation_full_cell<
Geom_traits,
internal::Triangulation::Dark_full_cell_data<Self> > Dark_full_cell_base;
typedef Triangulation_data_structure<Maximal_dimension,
Dark_vertex_base,
Dark_full_cell_base
> Dark_tds;
typedef Regular_triangulation<RTTraits, Dark_tds> Dark_triangulation;
typedef typename Dark_triangulation::Face Dark_face;
typedef typename Dark_triangulation::Facet Dark_facet;
typedef typename Dark_triangulation::Vertex_handle Dark_v_handle;
typedef typename Dark_triangulation::Full_cell_handle Dark_s_handle;
// If flat_orientation_ is defined, we give it the Dark triangulation
// so that the orientation it uses for "current_dimension()"-simplices is
// coherent with the global triangulation
Dark_triangulation dark_side(
maximal_dimension(),
flat_orientation_ ?
std::pair<int, const Flat_orientation_d *>(current_dimension(), flat_orientation_.get_ptr())
: std::pair<int, const Flat_orientation_d *>(std::numeric_limits<int>::max(), NULL) );
Dark_s_handle dark_s;
Dark_v_handle dark_v;
typedef std::map<Vertex_handle, Dark_v_handle> Vertex_map;
Vertex_map light_to_dark;
typename Vertex_set::iterator vit = verts.begin();
while( vit != verts.end() )
{
dark_v = dark_side.insert((*vit)->point(), dark_s);
dark_s = dark_v->full_cell();
dark_v->data() = *vit;
light_to_dark[*vit] = dark_v;
++vit;
}
if( dark_side.current_dimension() != current_dimension() )
{
CGAL_assertion( dark_side.current_dimension() + 1 == current_dimension() );
// Here, the finite neighbors of |v| span a affine subspace of
// dimension one less than the current dimension. Two cases are possible:
if( (size_type)(verts.size() + 1) == number_of_vertices() )
{
remove_decrease_dimension(v);
return Full_cell_handle();
}
else
{ // |v| is strictly outside the convex hull of the rest of the points. This is an
// easy case: first, modify the finite full_cells, then, delete the infinite ones.
// We don't even need the Dark triangulation.
Simplices infinite_simps;
{
Simplices finite_simps;
for( typename Simplices::iterator it = simps.begin(); it != simps.end(); ++it )
if( is_infinite(*it) )
infinite_simps.push_back(*it);
else
finite_simps.push_back(*it);
simps.swap(finite_simps);
} // now, simps only contains finite simplices
// First, modify the finite full_cells:
for( typename Simplices::iterator it = simps.begin(); it != simps.end(); ++it )
{
int v_idx = (*it)->index(v);
tds().associate_vertex_with_full_cell(*it, v_idx, infinite_vertex());
if( v_idx != 0 )
{
// we must put the infinite vertex at index 0.
// OK, now with the new convention that the infinite vertex
// does not have to be at index 0, this is not necessary,
// but still, I prefer to keep this piece of code here. [-- Samuel Hornus]
(*it)->swap_vertices(0, v_idx);
// Now, we preserve the positive orientation of the full_cell
(*it)->swap_vertices(current_dimension() - 1, current_dimension());
}
}
// Make the handles to infinite full cells searchable
infinite_simps.make_searchable();
// Then, modify the neighboring relation
for( typename Simplices::iterator it = simps.begin(); it != simps.end(); ++it )
{
for( int i = 1; i <= current_dimension(); ++i )
{
(*it)->vertex(i)->set_full_cell(*it);
Full_cell_handle n = (*it)->neighbor(i);
// Was |n| a finite full cell prior to removing |v| ?
if( ! infinite_simps.contains(n) )
continue;
int n_idx = n->index(v);
set_neighbors(*it, i, n->neighbor(n_idx), n->neighbor(n_idx)->index(n));
}
}
Full_cell_handle ret_s;
// Then, we delete the infinite full_cells
for( typename Simplices::iterator it = infinite_simps.begin(); it != infinite_simps.end(); ++it )
tds().delete_full_cell(*it);
tds().delete_vertex(v);
return simps.front();
}
}
else // From here on, dark_side.current_dimension() == current_dimension()
{
dark_side.infinite_vertex()->data() = infinite_vertex();
light_to_dark[infinite_vertex()] = dark_side.infinite_vertex();
}
// Now, compute the conflict zone of v->point() in
// the dark side. This is precisely the set of full_cells
// that we have to glue back into the light side.
Dark_face dark_f(dark_side.maximal_dimension());
Dark_facet dark_ft;
typename Dark_triangulation::Locate_type lt;
dark_s = dark_side.locate(v->point(), lt, dark_f, dark_ft);
CGAL_assertion( lt != Dark_triangulation::ON_VERTEX
&& lt != Dark_triangulation::OUTSIDE_AFFINE_HULL );
// |ret_s| is the full_cell that we return
Dark_s_handle dark_ret_s = dark_s;
Full_cell_handle ret_s;
typedef typename Base::template Full_cell_set<Dark_s_handle> Dark_full_cells;
Dark_full_cells conflict_zone;
std::back_insert_iterator<Dark_full_cells> dark_out(conflict_zone);
dark_ft = dark_side.compute_conflict_zone(v->point(), dark_s, dark_out);
// Make the dark simplices in the conflict zone searchable
conflict_zone.make_searchable();
// THE FOLLOWING SHOULD MAYBE GO IN TDS.
// Here is the plan:
// 1. Pick any Facet from boundary of the light zone
// 2. Find corresponding Facet on boundary of dark zone
// 3. stitch.
// 1. Build a facet on the boudary of the light zone:
Full_cell_handle light_s = *simps.begin();
Facet light_ft(light_s, light_s->index(v));
// 2. Find corresponding Dark_facet on boundary of the dark zone
Dark_full_cells dark_incident_s;
for( int i = 0; i <= current_dimension(); ++i )
{
if( index_of_covertex(light_ft) == i )
continue;
Dark_v_handle dark_v = light_to_dark[full_cell(light_ft)->vertex(i)];
dark_incident_s.clear();
dark_out = std::back_inserter(dark_incident_s);
dark_side.tds().incident_full_cells(dark_v, dark_out);
for(typename Dark_full_cells::iterator it = dark_incident_s.begin();
it != dark_incident_s.end();
++it)
{
(*it)->data().count_ += 1;
}
}
for( typename Dark_full_cells::iterator it = dark_incident_s.begin(); it != dark_incident_s.end(); ++it )
{
if( current_dimension() != (*it)->data().count_ )
continue;
if( ! conflict_zone.contains(*it) )
continue;
// We found a full_cell incident to the dark facet corresponding to the light facet |light_ft|
int ft_idx = 0;
while( light_s->has_vertex( (*it)->vertex(ft_idx)->data() ) )
++ft_idx;
dark_ft = Dark_facet(*it, ft_idx);
break;
}
// Pre-3. Now, we are ready to traverse both boundary and do the stiching.
// But first, we create the new full_cells in the light triangulation,
// with as much adjacency information as possible.
// Create new full_cells with vertices
for( typename Dark_full_cells::iterator it = conflict_zone.begin(); it != conflict_zone.end(); ++it )
{
Full_cell_handle new_s = new_full_cell();
(*it)->data().light_copy_ = new_s;
for( int i = 0; i <= current_dimension(); ++i )
tds().associate_vertex_with_full_cell(new_s, i, (*it)->vertex(i)->data());
if( dark_ret_s == *it )
ret_s = new_s;
}
// Setup adjacencies inside the hole
for( typename Dark_full_cells::iterator it = conflict_zone.begin(); it != conflict_zone.end(); ++it )
{
Full_cell_handle new_s = (*it)->data().light_copy_;
for( int i = 0; i <= current_dimension(); ++i )
if( conflict_zone.contains((*it)->neighbor(i)) )
tds().set_neighbors(new_s, i, (*it)->neighbor(i)->data().light_copy_, (*it)->mirror_index(i));
}
// 3. Stitch
simps.make_searchable();
typedef std::queue<std::pair<Facet, Dark_facet> > Queue;
Queue q;
q.push(std::make_pair(light_ft, dark_ft));
dark_s = dark_side.full_cell(dark_ft);
int dark_i = dark_side.index_of_covertex(dark_ft);
// mark dark_ft as visited:
// TODO try by marking with Dark_v_handle (vertex)
dark_s->neighbor(dark_i)->set_neighbor(dark_s->mirror_index(dark_i), Dark_s_handle());
while( ! q.empty() )
{
std::pair<Facet, Dark_facet> p = q.front();
q.pop();
light_ft = p.first;
dark_ft = p.second;
light_s = full_cell(light_ft);
int light_i = index_of_covertex(light_ft);
dark_s = dark_side.full_cell(dark_ft);
int dark_i = dark_side.index_of_covertex(dark_ft);
Full_cell_handle light_n = light_s->neighbor(light_i);
set_neighbors(dark_s->data().light_copy_, dark_i, light_n, light_s->mirror_index(light_i));
for( int di = 0; di <= current_dimension(); ++di )
{
if( di == dark_i )
continue;
int li = light_s->index(dark_s->vertex(di)->data());
Rotor light_r(light_s, li, light_i);
typename Dark_triangulation::Rotor dark_r(dark_s, di, dark_i);
while( simps.contains(full_cell(light_r)->neighbor(index_of_covertex(light_r))) )
light_r = rotate_rotor(light_r);
while( conflict_zone.contains(dark_side.full_cell(dark_r)->neighbor(dark_side.index_of_covertex(dark_r))) )
dark_r = dark_side.rotate_rotor(dark_r);
Dark_s_handle dark_ns = dark_side.full_cell(dark_r);
int dark_ni = dark_side.index_of_covertex(dark_r);
Full_cell_handle light_ns = full_cell(light_r);
int light_ni = index_of_covertex(light_r);
// mark dark_r as visited:
// TODO try by marking with Dark_v_handle (vertex)
Dark_s_handle outside = dark_ns->neighbor(dark_ni);
Dark_v_handle mirror = dark_ns->mirror_vertex(dark_ni, current_dimension());
int dn = outside->index(mirror);
if( Dark_s_handle() == outside->neighbor(dn) )
continue;
outside->set_neighbor(dn, Dark_s_handle());
q.push(std::make_pair(Facet(light_ns, light_ni), Dark_facet(dark_ns, dark_ni)));
}
}
tds().delete_full_cells(simps.begin(), simps.end());
tds().delete_vertex(v);
return ret_s;
}
template< typename RTTraits, typename TDS >
void
Regular_triangulation<RTTraits, TDS>
::remove_decrease_dimension(Vertex_handle v)
{
CGAL_precondition( current_dimension() >= 0 );
tds().remove_decrease_dimension(v, infinite_vertex());
// reset the predicates:
reset_flat_orientation();
if( 1 <= current_dimension() )
{
// FIXME: infinite vertex is NOT at index 0 a priori.
Full_cell_handle s = infinite_vertex()->full_cell()->neighbor(0);
Orientation o = orientation(s);
CGAL_assertion( ZERO != o );
if( NEGATIVE == o )
reorient_full_cells();
}
}
// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - INSERTIONS
template< typename RTTraits, typename TDS >
typename Regular_triangulation<RTTraits, TDS>::Vertex_handle
Regular_triangulation<RTTraits, TDS>
::insert(const Weighted_point & p, const Locate_type lt, const Face & f, const Facet & ft, const Full_cell_handle s)
{
switch( lt )
{
case Base::OUTSIDE_AFFINE_HULL:
return insert_outside_affine_hull(p);
break;
case Base::ON_VERTEX:
{
Vertex_handle v = s->vertex(f.index(0));
v->set_point(p);
return v;
break;
}
default:
if( 1 == current_dimension() )
{
if( Base::OUTSIDE_CONVEX_HULL == lt )
{
return insert_outside_convex_hull_1(p, s);
}
Vertex_handle v = tds().insert_in_full_cell(s);
v->set_point(p);
return v;
}
else
return insert_in_conflicting_cell(p, s);
break;
}
}
template< typename RTTraits, typename TDS >
typename Regular_triangulation<RTTraits, TDS>::Vertex_handle
Regular_triangulation<RTTraits, TDS>
::insert_outside_affine_hull(const Weighted_point & p)
{
// we don't use Base::insert_outside_affine_hull(...) because here, we
// also need to reset the side_of_oriented_subsphere functor.
CGAL_precondition( current_dimension() < maximal_dimension() );
Vertex_handle v = tds().insert_increase_dimension(infinite_vertex());
// reset the predicates:
reset_flat_orientation();
v->set_point(p);
if( current_dimension() >= 1 )
{
// FIXME: infinite vertex is NOT at index 0 a priori.
Full_cell_handle s = infinite_vertex()->full_cell()->neighbor(0);
Orientation o = orientation(s);
CGAL_assertion( ZERO != o );
if( NEGATIVE == o )
reorient_full_cells();
}
return v;
}
template< typename RTTraits, typename TDS >
typename Regular_triangulation<RTTraits, TDS>::Vertex_handle
Regular_triangulation<RTTraits, TDS>
::insert_in_conflicting_cell(const Weighted_point & p, const Full_cell_handle s)
{
typedef std::vector<Full_cell_handle> Full_cell_h_vector;
static Full_cell_h_vector cs; // for storing conflicting full_cells.
cs.clear();
// cs.reserve(64);
std::back_insert_iterator<Full_cell_h_vector> out(cs);
Facet ft = compute_conflict_zone(p, s, out);
return insert_in_hole(p, cs.begin(), cs.end(), ft);
}
// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - GATHERING CONFLICTING SIMPLICES
// NOT DOCUMENTED
template< typename RTTraits, typename TDS >
template< typename OrientationPred >
Oriented_side
Regular_triangulation<RTTraits, TDS>
::perturbed_power_test(const Weighted_point & p, Full_cell_const_handle s,
const OrientationPred & ori) const
{
CGAL_precondition_msg( ! is_infinite(s), "full cell must be finite");
CGAL_expensive_precondition( POSITIVE == orientation(s) );
typedef std::vector<const Weighted_point *> Points;
Points points(current_dimension() + 2);
int i(0);
for( ; i <= current_dimension(); ++i )
points[i] = &(s->vertex(i)->point());
points[i] = &p;
std::sort(points.begin(), points.end(),
internal::Triangulation::Compare_points_for_perturbation<Self>(*this));
typename Points::const_reverse_iterator cut_pt = points.rbegin();
Points test_points;
while( cut_pt != points.rend() )
{
if( &p == *cut_pt )
// because the full_cell "s" is assumed to be positively oriented
return ON_NEGATIVE_SIDE; // we consider |p| to lie outside the sphere
test_points.clear();
Point_const_iterator spit = points_begin(s);
int adjust_sign = -1;
for( i = 0; i < current_dimension(); ++i )
{
if( &(*spit) == *cut_pt )
{
++spit;
adjust_sign = (((current_dimension() + i) % 2) == 0) ? -1 : +1;
}
test_points.push_back(&(*spit));
++spit;
}
test_points.push_back(&p);
typedef typename CGAL::Iterator_project<
typename Points::iterator,
internal::Triangulation::Point_from_pointer<Self>,
const Weighted_point &, const Weighted_point *
> Point_pointer_iterator;
Orientation ori_value = ori(
Point_pointer_iterator(test_points.begin()),
Point_pointer_iterator(test_points.end()));
if( ZERO != ori_value )
return Oriented_side( - adjust_sign * ori_value );
++cut_pt;
}
CGAL_assertion(false); // we should never reach here
return ON_NEGATIVE_SIDE;
}
template< typename RTTraits, typename TDS >
bool
Regular_triangulation<RTTraits, TDS>
::is_in_conflict(const Weighted_point & p, Full_cell_const_handle s) const
{
CGAL_precondition( 2 <= current_dimension() );
if( current_dimension() < maximal_dimension() )
{
Conflict_pred_in_subspace c(
*this, p,
coaffine_orientation_predicate(),
power_test_in_flat_predicate());
return c(s);
}
else
{
Orientation_d ori = geom_traits().orientation_d_object();
Power_test_d side = geom_traits().power_test_d_object();
Conflict_pred_in_fullspace c(*this, p, ori, side);
return c(s);
}
}
template< typename RTTraits, typename TDS >
template< typename OutputIterator >
typename Regular_triangulation<RTTraits, TDS>::Facet
Regular_triangulation<RTTraits, TDS>
::compute_conflict_zone(const Weighted_point & p, const Full_cell_handle s, OutputIterator out) const
{
CGAL_precondition( 2 <= current_dimension() );
if( current_dimension() < maximal_dimension() )
{
Conflict_pred_in_subspace c(
*this, p,
coaffine_orientation_predicate(),
power_test_in_flat_predicate());
Conflict_traversal_pred_in_subspace tp(*this, c);
return tds().gather_full_cells(s, tp, out);
}
else
{
Orientation_d ori = geom_traits().orientation_d_object();
Power_test_d side = geom_traits().power_test_d_object();
Conflict_pred_in_fullspace c(*this, p, ori, side);
Conflict_traversal_pred_in_fullspace tp(*this, c);
return tds().gather_full_cells(s, tp, out);
}
}
// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - VALIDITY
template< typename RTTraits, typename TDS >
bool
Regular_triangulation<RTTraits, TDS>
::is_valid(bool verbose, int level) const
{
if (!Base::is_valid(verbose, level))
return false;
int dim = current_dimension();
if (dim == maximal_dimension())
{
for (Finite_full_cell_const_iterator cit = finite_full_cells_begin() ;
cit != finite_full_cells_end() ; ++cit )
{
Full_cell_const_handle ch = cit.base();
for(int i = 0; i < dim+1 ; ++i )
{
// If the i-th neighbor is not an infinite cell
Vertex_handle opposite_vh =
ch->neighbor(i)->vertex(ch->neighbor(i)->index(ch));
if (!is_infinite(opposite_vh))
{
Power_test_d side =
geom_traits().power_test_d_object();
if (side(Point_const_iterator(ch->vertices_begin()),
Point_const_iterator(ch->vertices_end()),
opposite_vh->point()) == ON_POSITIVE_SIDE)
{
if (verbose)
CGAL_warning_msg(false, "Non-empty sphere");
return false;
}
}
}
}
}
return true;
}
} //namespace CGAL
#endif CGAL_REGULAR_TRIANGULATION_H
#endif //CGAL_REGULAR_TRIANGULATION_H

View File

@ -0,0 +1,256 @@
// Copyright (c) 2014 INRIA Sophia-Antipolis (France).
// All rights reserved.
//
// This file is part of CGAL (www.cgal.org).
// You can redistribute it and/or modify it under the terms of the GNU
// General Public License as published by the Free Software Foundation,
// either version 3 of the License, or (at your option) any later version.
//
// 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) : Clement Jamin
#ifndef CGAL_REGULAR_TRIANGULATION_EUCLIDEAN_TRAITS_H
#define CGAL_REGULAR_TRIANGULATION_EUCLIDEAN_TRAITS_H
#include <CGAL/basic.h>
#include <CGAL/triangulation_assertions.h>
#include <CGAL/Weighted_point.h>
#include <CGAL/representation_tags.h>
#include <CGAL/Kernel_traits.h>
#include <boost/iterator/transform_iterator.hpp>
namespace CGAL {
template < class K, class Weight = typename K::RT >
class Regular_triangulation_euclidean_traits
: public K
{
public:
typedef K Base;
typedef Regular_triangulation_euclidean_traits<K, Weight> Self;
// Types from K
typedef K Kernel;
typedef typename K::Dimension Dimension;
typedef typename K::FT FT;
typedef typename K::Point_d Bare_point;
typedef typename K::Weighted_point_d Weighted_point;
typedef Weighted_point Weighted_point_d;
typedef Weighted_point Point_d;
typedef typename K::Construct_weighted_point_d Construct_weighted_point_d;
typedef typename K::Power_test_d Power_test_d;
typedef typename K::In_flat_power_test_d In_flat_power_test_d;
typedef typename K::Flat_orientation_d Flat_orientation_d;
typedef typename K::Point_drop_weight_d Point_drop_weight_d;
//=============================================================================
// Custom types
//=============================================================================
class Orientation_d
{
const K &m_kernel;
public:
typedef Orientation result_type;
Orientation_d(const K &kernel)
: m_kernel(kernel) {}
template <typename ForwardIterator>
result_type operator()(ForwardIterator start, ForwardIterator end) const
{
Point_drop_weight_d pdw = m_kernel.point_drop_weight_d_object();
return m_kernel.orientation_d_object() (
boost::make_transform_iterator(start, pdw),
boost::make_transform_iterator(end, pdw)
);
}
};
//=============================================================================
class Construct_flat_orientation_d
{
const K &m_kernel;
public:
typedef Flat_orientation_d result_type;
Construct_flat_orientation_d(const K &kernel)
: m_kernel(kernel) {}
template <typename ForwardIterator>
result_type operator()(ForwardIterator start, ForwardIterator end) const
{
Point_drop_weight_d pdw = m_kernel.point_drop_weight_d_object();
return m_kernel.construct_flat_orientation_d_object() (
boost::make_transform_iterator(start, pdw),
boost::make_transform_iterator(end, pdw)
);
}
};
//=============================================================================
class In_flat_orientation_d
{
const K &m_kernel;
public:
typedef Orientation result_type;
In_flat_orientation_d(const K &kernel)
: m_kernel(kernel) {}
template <typename ForwardIterator>
result_type operator()(Flat_orientation_d orient,
ForwardIterator start, ForwardIterator end) const
{
Point_drop_weight_d pdw = m_kernel.point_drop_weight_d_object();
return m_kernel.in_flat_orientation_d_object() (
orient,
boost::make_transform_iterator(start, pdw),
boost::make_transform_iterator(end, pdw)
);
}
};
//=============================================================================
class Contained_in_affine_hull_d
{
const K &m_kernel;
public:
typedef bool result_type;
Contained_in_affine_hull_d(const K &kernel)
: m_kernel(kernel) {}
template <typename ForwardIterator>
result_type operator()(ForwardIterator start, ForwardIterator end,
const Weighted_point_d & p) const
{
Point_drop_weight_d pdw = m_kernel.point_drop_weight_d_object();
return m_kernel.contained_in_affine_hull_d_object() (
boost::make_transform_iterator(start, pdw),
boost::make_transform_iterator(end, pdw),
pdw(p)
);
}
};
//=============================================================================
class Compare_lexicographically_d
{
const K &m_kernel;
public:
typedef Comparison_result result_type;
Compare_lexicographically_d(const K &kernel)
: m_kernel(kernel) {}
result_type operator()(
const Weighted_point_d & p, const Weighted_point_d & q) const
{
Point_drop_weight_d pdw = m_kernel.point_drop_weight_d_object();
return m_kernel.compare_lexicographically_d_object()(pdw(p), pdw(q));
}
};
//=============================================================================
class Compute_coordinate_d
{
const K &m_kernel;
public:
typedef FT result_type;
Compute_coordinate_d(const K &kernel)
: m_kernel(kernel) {}
result_type operator()(
const Weighted_point_d & p, const int i) const
{
Point_drop_weight_d pdw = m_kernel.point_drop_weight_d_object();
auto pp = pdw(p);
auto ddsd = m_kernel.compute_coordinate_d_object();
ddsd(pp, i);
return m_kernel.compute_coordinate_d_object()(pdw(p), i);
}
};
//=============================================================================
class Point_dimension_d
{
const K &m_kernel;
public:
typedef int result_type;
Point_dimension_d(const K &kernel)
: m_kernel(kernel) {}
result_type operator()(
const Weighted_point_d & p) const
{
Point_drop_weight_d pdw = m_kernel.point_drop_weight_d_object();
return m_kernel.point_dimension_d_object()(pdw(p));
}
};
//=============================================================================
// Object creation
//=============================================================================
Contained_in_affine_hull_d contained_in_affine_hull_d_object() const
{
return Contained_in_affine_hull_d(*this);
}
Orientation_d orientation_d_object() const
{
return Orientation_d(*this);
}
Construct_flat_orientation_d construct_flat_orientation_d_object() const
{
return Construct_flat_orientation_d(*this);
}
In_flat_orientation_d in_flat_orientation_d_object() const
{
return In_flat_orientation_d(*this);
}
Compare_lexicographically_d compare_lexicographically_d_object() const
{
return Compare_lexicographically_d(*this);
}
Compute_coordinate_d compute_coordinate_d_object() const
{
return Compute_coordinate_d(*this);
}
Point_dimension_d point_dimension_d_object() const
{
return Point_dimension_d(*this);
}
};
} //namespace CGAL
#endif // CGAL_REGULAR_TRIANGULATION_EUCLIDEAN_TRAITS_H

File diff suppressed because it is too large Load Diff

File diff suppressed because it is too large Load Diff

View File

@ -128,7 +128,7 @@ operator>>(std::istream & is, No_full_cell_data &)
}
std::ostream &
operator<<(std::ostream & os, const No_full_cell_data & nd)
operator<<(std::ostream & os, const No_full_cell_data &)
{
return os;
}

View File

@ -100,7 +100,7 @@ operator>>(std::istream & is, No_vertex_data &)
}
std::ostream &
operator<<(std::ostream & os, const No_vertex_data & nd)
operator<<(std::ostream & os, const No_vertex_data &)
{
return os;
}

View File

@ -98,7 +98,7 @@ public:
template< class T >
struct Compare_points_for_perturbation
{
typedef typename T::Point_d Point;
typedef typename T::Geom_traits::Point_d Point;
const T & t_;
@ -119,8 +119,8 @@ public:
template< class T >
struct Point_from_pointer
{
typedef const typename T::Point_d * argument_type;
typedef const typename T::Point_d result_type;
typedef const typename T::Geom_traits::Point_d * argument_type;
typedef const typename T::Geom_traits::Point_d result_type;
result_type & operator()(argument_type & x) const
{
return (*x);

View File

@ -0,0 +1,133 @@
#include <CGAL/Epick_d.h>
#include <CGAL/point_generators_d.h>
#include <CGAL/Regular_triangulation.h>
#include <CGAL/Regular_triangulation_euclidean_traits.h>
#include <CGAL/algorithm.h>
#include <tilted_grid.h>
#include <vector>
#include <string>
#include <fstream>
#include <cstdlib>
#include <algorithm>
using namespace std;
template<typename RTri>
void test(const int d, const string & type, const int N)
{
typedef typename RTri::Full_cell_handle Full_cell_handle;
typedef typename RTri::Face Face;
typedef typename RTri::Point Point;
typedef typename RTri::Bare_point Bare_point;
typedef typename RTri::Finite_full_cell_const_iterator Finite_full_cell_const_iterator;
typedef typename RTri::Finite_vertex_iterator Finite_vertex_iterator;
typedef CGAL::Random_points_in_cube_d<Bare_point> Random_points_iterator;
RTri rt(d);
cerr << "\nBuilding Regular triangulation of (" << type << d << ") dimension with " << N << " points";
assert(rt.empty());
vector<Point> points;
//CGAL::Random rng;
//Random_points_iterator rand_it(d, 2.0, rng); // CJTODO: unused
srand(10);
for( int i = 0; i < N; ++i )
{
vector<double> coords(d);
for( int j = 0; j < d; ++j )
coords[j] = static_cast<double>(rand() % 100000)/10000;
points.push_back(Point(
Bare_point(d, coords.begin(), coords.end()),
/*static_cast<double>(rand() % 100000)/100000*/static_cast<double>(i)/20
));
}
rt.insert(points.begin(), points.end());
cerr << "\nChecking topology and geometry...";
assert( rt.is_valid(true) );
cerr << "\nTraversing finite full_cells... ";
size_t nbfs(0), nbis(0);
Finite_full_cell_const_iterator fsit = rt.finite_full_cells_begin();
while( fsit != rt.finite_full_cells_end() )
++fsit, ++nbfs;
cerr << nbfs << " + ";
vector<Full_cell_handle> infinite_full_cells;
rt.tds().incident_full_cells(rt.infinite_vertex(), back_inserter(infinite_full_cells));
nbis = infinite_full_cells.size();
cerr << nbis << " = " << (nbis+nbfs)
<< " = " << rt.number_of_full_cells();
cerr << "\nThe triangulation has current dimension " << rt.current_dimension();
CGAL_assertion( rt.number_of_full_cells() == nbis+nbfs);
cerr << "\nTraversing finite vertices... ";
size_t nbfv(0);
Finite_vertex_iterator fvit = rt.finite_vertices_begin();
while( fvit != rt.finite_vertices_end() )
++fvit, ++nbfv;
cerr << nbfv <<endl;
// Count convex hull vertices:
if( rt.maximal_dimension() > 1 )
{
typedef vector<Face> Faces;
Faces edges;
back_insert_iterator<Faces> out(edges);
rt.tds().incident_faces(rt.infinite_vertex(), 1, out);
cout << "\nThere are " << edges.size() << " vertices on the convex hull.";
edges.clear();
}
else // rt.maximal_dimension() == 1
{
typedef vector<Full_cell_handle> Cells;
Cells cells;
back_insert_iterator<Cells> out(cells);
rt.tds().incident_full_cells(rt.infinite_vertex(), out);
cout << "\nThere are " << cells.size() << " vertices on the convex hull.";
cells.clear();
}
// Remove all !
cerr << "\nBefore removal: " << rt.number_of_vertices() << " vertices. After: ";
random_shuffle(points.begin(), points.end());
rt.remove(points.begin(), points.end());
assert( rt.is_valid() );
//std::cerr << ((rt.is_valid(true)) ? "VALID!" : "NOT VALID :(") << std::endl;
cerr << rt.number_of_vertices() << " vertices.";
// assert( rt.empty() ); NOT YET !
// CLEAR
rt.clear();
assert( -1 == rt.current_dimension() );
assert( rt.empty() );
assert( rt.is_valid() );
//std::cerr << ((rt.is_valid(true)) ? "VALID!" : "NOT VALID :(") << std::endl;
}
template< int D >
void go(const int N)
{
//typedef CGAL::Epick_d<CGAL::Dynamic_dimension_tag> FK;
typedef CGAL::Epick_d<CGAL::Dimension_tag<D> > FK;
typedef CGAL::Regular_triangulation_euclidean_traits<FK> Traits;
typedef CGAL::Regular_triangulation<Traits> Triangulation;
//test<Triangulation>(D, "dynamic", N);
test<Triangulation>(D, "static", N);
}
int main(int argc, char **argv)
{
srand(static_cast<unsigned int>(time(NULL)));
int N = 10;
if( argc > 1 )
N = atoi(argv[1]);
go<5>(N);
go<4>(N);
go<3>(N);
go<2>(N);
go<1>(N);
cerr << endl;
return 0;
}

View File

@ -19,92 +19,88 @@ void test(const int d, const string & type, const int N)
// we must write 'typename' below, because we are in a template-function,
// so the parser has no way to know that DC contains sub-types, before
// instanciating the function.
typedef typename DC::Vertex Vertex;
typedef typename DC::Vertex_handle Vertex_handle;
typedef typename DC::Full_cell Full_cell;
typedef typename DC::Full_cell_handle Full_cell_handle;
typedef typename DC::Facet Facet;
typedef typename DC::Face Face;
typedef typename DC::Point Point;
typedef typename DC::Geom_traits::RT RT;
typedef typename DC::Finite_full_cell_const_iterator Finite_full_cell_const_iterator;
typedef typename DC::Finite_vertex_iterator Finite_vertex_iterator;
typedef CGAL::Random_points_in_cube_d<Point> Random_points_iterator;
DC pc(d);
DC dt(d);
cerr << "\nBuilding Delaunay triangulation of (" << type << d << ") dimension with " << N << " points";
assert(pc.empty());
assert(dt.empty());
vector<Point> points;
CGAL::Random rng;
Random_points_iterator rand_it(d, 2.0, rng);
//CGAL::Random rng;
//Random_points_iterator rand_it(d, 2.0, rng);
//CGAL::cpp11::copy_n(rand_it, N, back_inserter(points));
vector<int> coords(d);
srand(10);
for( int i = 0; i < N; ++i )
{
vector<double> coords(d);
for( int j = 0; j < d; ++j )
coords[j] = rand() % 100000;
coords[j] = static_cast<double>(rand() % 100000)/10000;
points.push_back(Point(d, coords.begin(), coords.end()));
}
pc.insert(points.begin(), points.end());
dt.insert(points.begin(), points.end());
cerr << "\nChecking topology and geometry...";
assert( pc.is_valid() );
assert( dt.is_valid() );
cerr << "\nTraversing finite full_cells... ";
size_t nbfs(0), nbis(0);
Finite_full_cell_const_iterator fsit = pc.finite_full_cells_begin();
while( fsit != pc.finite_full_cells_end() )
Finite_full_cell_const_iterator fsit = dt.finite_full_cells_begin();
while( fsit != dt.finite_full_cells_end() )
++fsit, ++nbfs;
cerr << nbfs << " + ";
vector<Full_cell_handle> infinite_full_cells;
pc.tds().incident_full_cells(pc.infinite_vertex(), back_inserter(infinite_full_cells));
dt.tds().incident_full_cells(dt.infinite_vertex(), back_inserter(infinite_full_cells));
nbis = infinite_full_cells.size();
cerr << nbis << " = " << (nbis+nbfs)
<< " = " << pc.number_of_full_cells();
cerr << "\nThe triangulation has current dimension " << pc.current_dimension();
CGAL_assertion( pc.number_of_full_cells() == nbis+nbfs);
<< " = " << dt.number_of_full_cells();
cerr << "\nThe triangulation has current dimension " << dt.current_dimension();
CGAL_assertion( dt.number_of_full_cells() == nbis+nbfs);
cerr << "\nTraversing finite vertices... ";
size_t nbfv(0);
Finite_vertex_iterator fvit = pc.finite_vertices_begin();
while( fvit != pc.finite_vertices_end() )
Finite_vertex_iterator fvit = dt.finite_vertices_begin();
while( fvit != dt.finite_vertices_end() )
++fvit, ++nbfv;
cerr << nbfv <<endl;
// Count convex hull vertices:
if( pc.maximal_dimension() > 1 )
if( dt.maximal_dimension() > 1 )
{
typedef vector<Face> Faces;
Faces edges;
back_insert_iterator<Faces> out(edges);
pc.tds().incident_faces(pc.infinite_vertex(), 1, out);
dt.tds().incident_faces(dt.infinite_vertex(), 1, out);
cout << "\nThere are " << edges.size() << " vertices on the convex hull.";
edges.clear();
}
else // pc.maximal_dimension() == 1
else // dt.maximal_dimension() == 1
{
typedef vector<Full_cell_handle> Cells;
Cells cells;
back_insert_iterator<Cells> out(cells);
pc.tds().incident_full_cells(pc.infinite_vertex(), out);
dt.tds().incident_full_cells(dt.infinite_vertex(), out);
cout << "\nThere are " << cells.size() << " vertices on the convex hull.";
cells.clear();
}
// Remove all !
cerr << "\nBefore removal: " << pc.number_of_vertices() << " vertices. After: ";
cerr << "\nBefore removal: " << dt.number_of_vertices() << " vertices. After: ";
random_shuffle(points.begin(), points.end());
pc.remove(points.begin(), points.end());
assert( pc.is_valid() );
cerr << pc.number_of_vertices() << " vertices.";
// assert( pc.empty() ); NOT YET !
dt.remove(points.begin(), points.end());
assert( dt.is_valid() );
cerr << dt.number_of_vertices() << " vertices.";
// assert( dt.empty() ); NOT YET !
// CLEAR
pc.clear();
assert( -1 == pc.current_dimension() );
assert( pc.empty() );
assert( pc.is_valid() );
dt.clear();
assert( -1 == dt.current_dimension() );
assert( dt.empty() );
assert( dt.is_valid() );
}
template< int D >
@ -120,14 +116,14 @@ void go(const int N)
int main(int argc, char **argv)
{
srand(static_cast<unsigned int>(time(NULL)));
int N = 100;
int N = 10;
if( argc > 1 )
N = atoi(argv[1]);
go<5>(N);
go<4>(N);
go<3>(N);
go<2>(N);
go<1>(N);
//go<5>(N);
//go<4>(N);
//go<3>(N);
go<2>(N);
//go<1>(N);
cerr << endl;
return 0;

View File

@ -119,10 +119,10 @@ int main(int argc, char **argv)
int N = 1000;
if( argc > 1 )
N = atoi(argv[1]);
go<5>(N);
go<3>(N);
//go<5>(N);
//go<3>(N);
go<2>(N);
go<1>(N);
//go<1>(N);
cerr << std::endl;
return 0;

View File

@ -0,0 +1,10 @@
0.0071 1.6899 0
0.3272 1.3694 0.05
1.3697 1.8296 0.1
0.6722 0.3012 0.15
1.1726 0.1899 0.2
0.4374 2.8541 0.25
2.5923 0.1904 0.3
1.3083 2.5462 0.35
1.4981 1.3929 0.4
2.1304 2.055 0.45

View File

@ -0,0 +1,27 @@
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Regular_triangulation_euclidean_traits_2.h>
#include <CGAL/Regular_triangulation_filtered_traits_2.h>
#include <CGAL/Regular_triangulation_2.h>
#include <CGAL/IO/Triangulation_off_ostream_2.h>
#include <fstream>
typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
typedef CGAL::Regular_triangulation_filtered_traits_2<K> Traits;
typedef CGAL::Regular_triangulation_2<Traits> Regular_triangulation;
int main()
{
std::ifstream in("data/points.cin");
Regular_triangulation::Weighted_point wp;
std::vector<Regular_triangulation::Weighted_point> wpoints;
while(in >> wp)
wpoints.push_back(wp);
Regular_triangulation rt(wpoints.begin(), wpoints.end());
std::ofstream off_stream("data/rt2.off");
CGAL::export_triangulation_2_to_off(off_stream, rt);
return 0;
}

View File

@ -0,0 +1,79 @@
// Copyright (c) 2014 INRIA Sophia-Antipolis (France).
// All rights reserved.
//
// This file is part of CGAL (www.cgal.org); you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public License as
// published by the Free Software Foundation; either version 3 of the License,
// or (at your option) any later version.
//
// 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) : Clement Jamin
#ifndef CGAL_TRIANGULATION_OFF_OSTREAM_2_H
#define CGAL_TRIANGULATION_OFF_OSTREAM_2_H
#include <CGAL/Triangulation_2.h>
#include <sstream>
#include <iostream>
namespace CGAL {
template < class GT, class TDS >
std::ostream &
export_triangulation_2_to_off(std::ostream & os,
const Triangulation_2<GT,TDS> & tr)
{
typedef Triangulation_2<GT,TDS> Tr;
typedef typename Tr::Vertex_handle Vertex_handle;
typedef typename Tr::Vertex_iterator Vertex_iterator;
typedef typename Tr::Finite_vertices_iterator Finite_vertex_iterator;
typedef typename Tr::Finite_faces_iterator Finite_faces_iterator;
size_t n = tr.number_of_vertices();
std::stringstream output;
// write the vertices
std::map<Vertex_handle, int> index_of_vertex;
int i = 0;
for(Finite_vertex_iterator it = tr.finite_vertices_begin();
it != tr.finite_vertices_end(); ++it, ++i)
{
output << it->point().x() << " " << it->point().y() << " 0" << std::endl;
index_of_vertex[it.base()] = i;
}
CGAL_assertion( i == n );
size_t number_of_triangles = 0;
for (Finite_faces_iterator fit = tr.finite_faces_begin() ;
fit != tr.finite_faces_end() ; ++fit)
{
output << "3 "
<< index_of_vertex[fit->vertex(0)] << " "
<< index_of_vertex[fit->vertex(1)] << " "
<< index_of_vertex[fit->vertex(2)]
<< std::endl;
++number_of_triangles;
}
os << "OFF \n"
<< n << " "
<< number_of_triangles << " 0\n"
<< output.str();
return os;
}
} //namespace CGAL
#endif // CGAL_TRIANGULATION_OFF_OSTREAM_2_H

View File

@ -0,0 +1,10 @@
0.0071 1.6899 2.521 0
0.3272 1.3694 3.15 0.05
1.3697 1.8296 2.654 0.1
0.6722 0.3012 0.1548 0.15
1.1726 0.1899 0.3658 0.2
0.4374 2.8541 1.45894 0.25
2.5923 0.1904 0.6971 0.3
1.3083 2.5462 1.3658 0.35
1.4981 1.3929 2.949 0.4
2.1304 2.055 0.6597455 0.45

View File

@ -0,0 +1,26 @@
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Regular_triangulation_3.h>
#include <CGAL/Regular_triangulation_euclidean_traits_3.h>
#include <CGAL/IO/Triangulation_off_ostream_3.h>
#include <fstream>
typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
typedef CGAL::Regular_triangulation_euclidean_traits_3<K> Traits;
typedef CGAL::Regular_triangulation_3<Traits> Regular_triangulation;
int main()
{
std::ifstream in("data/points.cin");
Regular_triangulation::Weighted_point wp;
std::vector<Regular_triangulation::Weighted_point> wpoints;
while(in >> wp)
wpoints.push_back(wp);
Regular_triangulation rt(wpoints.begin(), wpoints.end());
std::ofstream off_stream("data/rt3.off");
CGAL::export_triangulation_3_to_off(off_stream, rt);
return 0;
}

View File

@ -0,0 +1,119 @@
// Copyright (c) 2014 INRIA Sophia-Antipolis (France).
// All rights reserved.
//
// This file is part of CGAL (www.cgal.org); you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public License as
// published by the Free Software Foundation; either version 3 of the License,
// or (at your option) any later version.
//
// 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) : Clement Jamin
#ifndef CGAL_TRIANGULATION_OFF_OSTREAM_3_H
#define CGAL_TRIANGULATION_OFF_OSTREAM_3_H
#include <CGAL/Triangulation_3.h>
#include <sstream>
#include <iostream>
namespace CGAL {
template < class GT, class TDS >
std::ostream &
export_triangulation_3_to_off(std::ostream & os,
const Triangulation_3<GT,TDS> & tr,
bool export_surface_only = false)
{
typedef Triangulation_3<GT,TDS> Tr;
typedef typename Tr::Vertex_handle Vertex_handle;
typedef typename Tr::Vertex_iterator Vertex_iterator;
typedef typename Tr::Finite_vertices_iterator Finite_vertex_iterator;
typedef typename Tr::All_cells_iterator Cells_iterator;
typedef typename Tr::Finite_cells_iterator Finite_cells_iterator;
size_t n = tr.number_of_vertices();
std::stringstream output;
// write the vertices
std::map<Vertex_handle, int> index_of_vertex;
int i = 0;
for(Finite_vertex_iterator it = tr.finite_vertices_begin();
it != tr.finite_vertices_end(); ++it, ++i)
{
output << it->point().x() << " "
<< it->point().y() << " "
<< it->point().z() << std::endl;
index_of_vertex[it.base()] = i;
}
CGAL_assertion( i == n );
size_t number_of_triangles = 0;
if (export_surface_only)
{
for (Cells_iterator cit = tr.cells_begin() ;
cit != tr.cells_end() ; ++cit)
{
if (tr.is_infinite(cit))
{
output << "3 ";
for (int i = 0 ; i < 4 ; ++i)
{
if (!tr.is_infinite(cit->vertex(i)))
output << index_of_vertex[cit->vertex(i)] << " ";
}
output << std::endl;
++number_of_triangles;
}
}
}
else
{
for (Finite_cells_iterator cit = tr.finite_cells_begin() ;
cit != tr.finite_cells_end() ; ++cit)
{
output << "3 "
<< index_of_vertex[cit->vertex(0)] << " "
<< index_of_vertex[cit->vertex(1)] << " "
<< index_of_vertex[cit->vertex(2)]
<< std::endl;
output << "3 "
<< index_of_vertex[cit->vertex(0)] << " "
<< index_of_vertex[cit->vertex(2)] << " "
<< index_of_vertex[cit->vertex(3)]
<< std::endl;
output << "3 "
<< index_of_vertex[cit->vertex(1)] << " "
<< index_of_vertex[cit->vertex(2)] << " "
<< index_of_vertex[cit->vertex(3)]
<< std::endl;
output << "3 "
<< index_of_vertex[cit->vertex(0)] << " "
<< index_of_vertex[cit->vertex(1)] << " "
<< index_of_vertex[cit->vertex(3)]
<< std::endl;
number_of_triangles += 4;
}
}
os << "OFF \n"
<< n << " "
<< number_of_triangles << " 0\n"
<< output.str();
return os;
}
} //namespace CGAL
#endif // CGAL_TRIANGULATION_OFF_OSTREAM_3_H