One more (also untested).

This commit is contained in:
Marc Glisse 2014-02-15 16:12:20 +01:00
parent 6fea55e611
commit f8a07a51da
3 changed files with 56 additions and 2 deletions

View File

@ -109,7 +109,7 @@ template <class Base_> struct Kernel_d_interface : public Base_ {
typedef typename Get_functor<Base, Affine_rank_tag>::type Affine_rank_d;
typedef typename Get_functor<Base, Affinely_independent_tag>::type Affinely_independent_d;
typedef typename Get_functor<Base, Contained_in_linear_hull_tag>::type Contained_in_linear_hull_d;
//typedef typename Get_functor<Base, Contained_in_simplex_tag>::type Contained_in_simplex_d;
typedef typename Get_functor<Base, Contained_in_simplex_tag>::type Contained_in_simplex_d;
typedef typename Get_functor<Base, Has_on_positive_side_tag>::type Has_on_positive_side_d;
typedef typename Get_functor<Base, Linear_rank_tag>::type Linear_rank_d;
typedef typename Get_functor<Base, Linearly_independent_tag>::type Linearly_independent_d;
@ -137,7 +137,7 @@ template <class Base_> struct Kernel_d_interface : public Base_ {
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); }
//Contained_in_simplex_d contained_in_simplex_d_object()const{ return Contained_in_simplex_d(*this); }
Contained_in_simplex_d contained_in_simplex_d_object()const{ return Contained_in_simplex_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); }
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); }

View File

@ -345,6 +345,54 @@ template<class R_> struct Affinely_independent : private Store_kernel<R_> {
CGAL_KD_DEFAULT_FUNCTOR(Affinely_independent_tag,(CartesianDKernelFunctors::Affinely_independent<K>),(Point_tag),(Point_dimension_tag,Affine_rank_tag));
namespace CartesianDKernelFunctors {
template<class R_> struct Contained_in_simplex : private Store_kernel<R_> {
CGAL_FUNCTOR_INIT_STORE(Contained_in_simplex)
typedef R_ R;
typedef typename Get_type<R, Point_tag>::type Point;
// Computing a sensible Uncertain<*> is not worth it
typedef typename Get_type<R, Bounded_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::Dynamic_matrix Matrix;
typedef typename LA::Dynamic_vector DynVec;
typedef typename LA::Vector Vec;
template<class Iter, class P>
result_type operator()(Iter f, Iter e, P const&q)const{
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 n=std::distance(f,e);
if (n==0) return ON_UNBOUNDED_SIDE;
int d=pd(q);
Matrix m(d+1,n);
DynVec a(n);
// FIXME: Should use the proper vector constructor (Iterator_and_last)
Vec b(d+1);
for(int j=0;j<d;++j) b[j]=c(q,j);
b[d]=1;
for(int i=0; ++f!=e; ++i){
Point const& p = *f;
for(int j=0;j<d;++j){
m(i,j)=c(p,j);
}
m(i,d)=1;
}
if (!LA::solve(a,CGAL_MOVE(m),CGAL_MOVE(b))) return false;
result_type res = ON_BOUNDED_SIDE;
for(int i=0;i<n;++i){
if (a[i]<0) return ON_UNBOUNDED_SIDE;
if (a[i]==0) res = ON_BOUNDARY;
}
return res;
}
};
}
CGAL_KD_DEFAULT_FUNCTOR(Contained_in_simplex_tag,(CartesianDKernelFunctors::Contained_in_simplex<K>),(Point_tag),(Point_dimension_tag,Compute_point_cartesian_coordinate_tag));
#if 0
namespace CartesianDKernelFunctors {
template<class R_,bool=boost::is_same<typename R_::Point,typename R_::Vector>::value> struct Orientation : private Store_kernel<R_> {

View File

@ -103,6 +103,12 @@ template<class NT_,class Dim_,class Max_dim_=Dim_> struct LA_eigen {
return decomp.rank();
}
// m*a==b
static bool solve(Dynamic_vector&a, Dynamic_matrix const&m, Vector const& b){
a = m.colPivHouseholderQr().solve(b);
return b.isApprox(m*a);
}
template<class Vec1,class Vec2> static Vector homogeneous_add(Vec1 const&a,Vec2 const&b){
//TODO: use compile-time size when available
int d=a.size();