Test Contained_in_simplex

This commit is contained in:
Marc Glisse 2014-05-06 16:04:51 +02:00
parent de3314319f
commit d218be48be
3 changed files with 11 additions and 9 deletions

View File

@ -117,7 +117,8 @@ template<class NT_,class Dim_,class Max_dim_=Dim_> struct LA_eigen {
// m*a==b // m*a==b
template<class DV, class DM, class V> template<class DV, class DM, class V>
static bool solve(DV&a, DM const&m, V const& b){ static bool solve(DV&a, DM const&m, V const& b){
a = m.colPivHouseholderQr().solve(b); //a = m.colPivHouseholderQr().solve(b);
a = m.fullPivLu().solve(b);
return b.isApprox(m*a); return b.isApprox(m*a);
} }

View File

@ -353,7 +353,8 @@ template<class R_> struct Contained_in_simplex : private Store_kernel<R_> {
typedef R_ R; typedef R_ R;
typedef typename Get_type<R, Point_tag>::type Point; typedef typename Get_type<R, Point_tag>::type Point;
// Computing a sensible Uncertain<*> is not worth it // Computing a sensible Uncertain<*> is not worth it
typedef typename Get_type<R, Bounded_side_tag>::type result_type; // typedef typename Get_type<R, Boolean_tag>::type result_type;
typedef bool result_type;
typedef typename Increment_dimension<typename R::Default_ambient_dimension>::type D1; typedef typename Increment_dimension<typename R::Default_ambient_dimension>::type D1;
typedef typename Increment_dimension<typename R::Max_ambient_dimension>::type D2; typedef typename Increment_dimension<typename R::Max_ambient_dimension>::type D2;
typedef typename R::LA::template Rebind_dimension<D1,D2>::Other LA; typedef typename R::LA::template Rebind_dimension<D1,D2>::Other LA;
@ -375,20 +376,18 @@ template<class R_> struct Contained_in_simplex : private Store_kernel<R_> {
for(int j=0;j<d;++j) b[j]=c(q,j); for(int j=0;j<d;++j) b[j]=c(q,j);
b[d]=1; b[d]=1;
for(int i=0; ++f!=e; ++i){ for(int i=0; f!=e; ++i,++f){
Point const& p = *f; Point const& p = *f;
for(int j=0;j<d;++j){ for(int j=0;j<d;++j){
m(i,j)=c(p,j); m(j,i)=c(p,j);
} }
m(i,d)=1; m(d,i)=1;
} }
if (!LA::solve(a,CGAL_MOVE(m),CGAL_MOVE(b))) return false; 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){ for(int i=0;i<n;++i){
if (a[i]<0) return ON_UNBOUNDED_SIDE; if (a[i]<0) return false;
if (a[i]==0) res = ON_BOUNDARY;
} }
return res; return true;
} }
}; };
} }

View File

@ -409,6 +409,7 @@ void test3(){
FO fo2=cfo(&x[0],x+3); FO fo2=cfo(&x[0],x+3);
std::cout << fo2; std::cout << fo2;
P y[]={cp(0,2,4),cp(3,1,2),cp(3,3,6),cp(0,4,8)}; P y[]={cp(0,2,4),cp(3,1,2),cp(3,3,6),cp(0,4,8)};
assert(!cis(x+0,x+3,y[0]));
V yv[]={cv(0,2,4),cv(3,1,2),cv(3,3,6),cv(0,4,8)}; V yv[]={cv(0,2,4),cv(3,1,2),cv(3,3,6),cv(0,4,8)};
assert( clh(yv+0,yv+1,yv[3])); assert( clh(yv+0,yv+1,yv[3]));
assert( clh(yv+0,yv+2,yv[2])); assert( clh(yv+0,yv+2,yv[2]));
@ -450,6 +451,7 @@ void test3(){
P x5=cp(0,0,0); P x5=cp(0,0,0);
P x6=cp(0,0,-1); P x6=cp(0,0,-1);
P tab2[]={x1,x2,x3,x4,x5}; P tab2[]={x1,x2,x3,x4,x5};
assert(cis(tab2+0,tab2+4,x5));
assert(po(tab2+0,tab2+4)==CGAL::POSITIVE); assert(po(tab2+0,tab2+4)==CGAL::POSITIVE);
assert(sos(tab2+0,tab2+4,x5)==CGAL::ON_POSITIVE_SIDE); assert(sos(tab2+0,tab2+4,x5)==CGAL::ON_POSITIVE_SIDE);
assert(sbs(tab2+0,tab2+4,x5)==CGAL::ON_BOUNDED_SIDE); assert(sbs(tab2+0,tab2+4,x5)==CGAL::ON_BOUNDED_SIDE);