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

Merge pull request #84 from mglisse/NewKernel_d-insphere-glisse 

Add functors to Epick_d

https://cgal.geometryfactory.com/CGAL/Members/wiki/Features/Small_Features/Side_of_bounded_diametral_sphere
This commit is contained in:
Laurent Rineau 2015-05-18 11:21:01 +02:00
parent a2ffb7d6f8 3a322b76e5
commit 81b7171da2
6 changed files with 211 additions and 35 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

31
Kernel_d/doc/Kernel_d/CGAL/Epick_d.h
View File

@ -73,7 +73,36 @@ Cartesian_const_iterator_d cartesian_begin()const;
/*! returns an iterator pointing beyond the last Cartesian coordinate. */
Cartesian_const_iterator_d cartesian_end()const;
};
/// @}
/*! \cgalModels `Kernel_d::Center_of_sphere_d`
*/
struct Construct_circumcenter_d {
/*! returns the center of the sphere defined by `A=tuple[first,last)`. The sphere is centered in the affine hull of A and passes through all the points of A. The order of the points of A does not matter.
\pre A is affinely independant.
\cgalRequires The value type of `ForwardIterator` is `Epick_d::Point_d`.
*/
template<typename ForwardIterator>
Point_d operator()(ForwardIterator first, ForwardIterator last);
};
struct Compute_squared_radius_d {
/*! returns the radius of the sphere defined by `A=tuple[first,last)`. The sphere is centered in the affine hull of A and passes through all the points of A. The order of the points of A does not matter.
\pre A is affinely independant.
\cgalRequires The value type of `ForwardIterator` is `Epick_d::Point_d`.
*/
template<class ForwardIterator>
Point_d operator()(ForwardIterator first, ForwardIterator last);
};
/*! \cgalModels `Kernel_d::Side_of_bounded_sphere_d`
*/
struct Side_of_bounded_sphere_d {
/*! returns the relative position of point p to the sphere defined by `A=tuple[first,last)`. The sphere is centered in the affine hull of A and passes through all the points of A. The order of the points of A does not matter.
\pre A is affinely independant.
\cgalRequires The value type of `ForwardIterator` is `Epick_d::Point_d`.
*/
template<class ForwardIterator>
Bounded_side operator()(ForwardIterator first, ForwardIterator last, const Point_d&p);
};
Construct_circumcenter_d construct_circumcenter_d_object();
Compute_squared_radius_d compute_squared_radius_d_object();
}; /* end Epick_d */
} /* end namespace CGAL */

16
NewKernel_d/include/CGAL/NewKernel_d/Kernel_d_interface.h
View File

@ -130,6 +130,17 @@ template <class Base_> struct Kernel_d_interface : public Base_ {
return operator()(v, End_tag());
}
};
struct Compute_squared_radius_d : private Store_kernel<Kernel> {
typedef Kernel R_; // for the macro
CGAL_FUNCTOR_INIT_STORE(Compute_squared_radius_d)
typedef FT result_type;
template<class S> FT operator()(CGAL_FORWARDABLE(S) s)const{
return typename Get_functor<Base, Squared_radius_tag>::type(this->kernel())(CGAL_FORWARD(S,s));
}
template<class I> FT operator()(I b, I e)const{
return typename Get_functor<Base, Squared_circumradius_tag>::type(this->kernel())(b,e);
}
};
typedef typename Construct_cartesian_const_iterator_d::result_type Cartesian_const_iterator_d;
typedef typename Get_functor<Base, Squared_distance_tag>::type Squared_distance_d;
typedef typename Get_functor<Base, Squared_length_tag>::type Squared_length_d;
@ -142,9 +153,10 @@ template <class Base_> struct Kernel_d_interface : public Base_ {
typedef typename Get_functor<Base, Linear_rank_tag>::type Linear_rank_d;
typedef typename Get_functor<Base, Linearly_independent_tag>::type Linearly_independent_d;
typedef typename Get_functor<Base, Oriented_side_tag>::type Oriented_side_d;
typedef typename Get_functor<Base, Side_of_bounded_sphere_tag>::type Side_of_bounded_sphere_d;
typedef typename Get_functor<Base, Side_of_bounded_circumsphere_tag>::type Side_of_bounded_sphere_d;
typedef typename Get_functor<Base, Center_of_sphere_tag>::type Center_of_sphere_d;
typedef typename Get_functor<Base, Construct_circumcenter_tag>::type Construct_circumcenter_d;
typedef typename Get_functor<Base, Value_at_tag>::type Value_at_d;
typedef typename Get_functor<Base, Point_of_sphere_tag>::type Point_of_sphere_d;
typedef typename Get_functor<Base, Orthogonal_vector_tag>::type Orthogonal_vector_d;
@ -196,10 +208,12 @@ template <class Base_> struct Kernel_d_interface : public Base_ {
Construct_segment_d construct_segment_d_object()const{ return Construct_segment_d(*this); }
Construct_sphere_d construct_sphere_d_object()const{ return Construct_sphere_d(*this); }
Construct_hyperplane_d construct_hyperplane_d_object()const{ return Construct_hyperplane_d(*this); }
Compute_squared_radius_d compute_squared_radius_d_object()const{ return Compute_squared_radius_d(*this); }
Squared_distance_d squared_distance_d_object()const{ return Squared_distance_d(*this); }
Squared_length_d squared_length_d_object()const{ return Squared_length_d(*this); }
Scalar_product_d scalar_product_d_object()const{ return Scalar_product_d(*this); }
Center_of_sphere_d center_of_sphere_d_object()const{ return Center_of_sphere_d(*this); }
Construct_circumcenter_d construct_circumcenter_d_object()const{ return Construct_circumcenter_d(*this); }
Construct_direction_d construct_direction_d_object()const{ return Construct_direction_d(*this); }
Construct_line_d construct_line_d_object()const{ return Construct_line_d(*this); }
Construct_ray_d construct_ray_d_object()const{ return Construct_ray_d(*this); }

34
NewKernel_d/include/CGAL/NewKernel_d/Types/Sphere.h
View File

@ -52,39 +52,13 @@ template <class R_> struct Construct_sphere : Store_kernel<R_> {
}
template <class Iter>
result_type operator()(Iter f, Iter e)const{
// 2*(x-y).c == x^2-y^2
typedef typename R_::LA LA;
typedef typename LA::Square_matrix Matrix;
typedef typename LA::Vector Vec;
typedef typename LA::Construct_vector CVec;
typedef typename Get_type<R_, Point_tag>::type Point;
typename Get_functor<R_, Compute_point_cartesian_coordinate_tag>::type c(this->kernel());
typename Get_functor<R_, Construct_ttag<Point_tag> >::type cp(this->kernel());
typename Get_functor<R_, Point_dimension_tag>::type pd(this->kernel());
typename Get_functor<R_, Squared_distance_to_origin_tag>::type sdo(this->kernel());
typename Get_functor<R_, Construct_circumcenter_tag>::type cc(this->kernel());
typename Get_functor<R_, Squared_distance_tag>::type sd(this->kernel());
Point const& p0=*f;
int d = pd(p0);
FT const& n0 = sdo(p0);
Matrix m(d,d);
Vec b = typename CVec::Dimension()(d);
// Write the point coordinates in lines.
int i;
for(i=0; ++f!=e; ++i) {
Point const& p=*f;
for(int j=0;j<d;++j) {
m(i,j)=2*(c(p,j)-c(p0,j));
b[i] = sdo(p) - n0;
}
}
CGAL_assertion (i == d);
Vec res = typename CVec::Dimension()(d);;
//std::cout << "Mat: " << m << "\n Vec: " << one << std::endl;
LA::solve(res, CGAL_MOVE(m), CGAL_MOVE(b));
//std::cout << "Sol: " << res << std::endl;
Point center = cp(d,LA::vector_begin(res),LA::vector_end(res));
FT const& r2 = sd (center, p0);
// It should be possible to avoid copying the center by moving this code to a constructor.
Point center = cc(f, e);
FT const& r2 = sd(center, *f);
return this->operator()(CGAL_MOVE(center), r2);
}
};

142
NewKernel_d/include/CGAL/NewKernel_d/function_objects_cartesian.h
View File

@ -620,6 +620,123 @@ template<class R_> struct Side_of_oriented_sphere : private Store_kernel<R_> {
CGAL_KD_DEFAULT_FUNCTOR(Side_of_oriented_sphere_tag,(CartesianDKernelFunctors::Side_of_oriented_sphere<K>),(Point_tag),(Point_dimension_tag,Squared_distance_to_origin_tag,Compute_point_cartesian_coordinate_tag));
namespace CartesianDKernelFunctors {
template <class R_> struct Construct_circumcenter : Store_kernel<R_> {
CGAL_FUNCTOR_INIT_STORE(Construct_circumcenter)
typedef typename Get_type<R_, Point_tag>::type Point;
typedef Point result_type;
typedef typename Get_type<R_, FT_tag>::type FT;
template <class Iter>
result_type operator()(Iter f, Iter e)const{
typedef typename Get_type<R_, Point_tag>::type Point;
typedef typename R_::LA LA;
typename Get_functor<R_, Compute_point_cartesian_coordinate_tag>::type c(this->kernel());
typename Get_functor<R_, Construct_ttag<Point_tag> >::type cp(this->kernel());
typename Get_functor<R_, Point_dimension_tag>::type pd(this->kernel());
typename Get_functor<R_, Squared_distance_to_origin_tag>::type sdo(this->kernel());
Point const& p0=*f;
int d = pd(p0);
if (d+1 == std::distance(f,e))
{
// 2*(x-y).c == x^2-y^2
typedef typename LA::Square_matrix Matrix;
typedef typename LA::Vector Vec;
typedef typename LA::Construct_vector CVec;
FT const& n0 = sdo(p0);
Matrix m(d,d);
Vec b = typename CVec::Dimension()(d);
// Write the point coordinates in lines.
int i;
for(i=0; ++f!=e; ++i) {
Point const& p=*f;
for(int j=0;j<d;++j) {
m(i,j)=2*(c(p,j)-c(p0,j));
b[i] = sdo(p) - n0;
}
}
CGAL_assertion (i == d);
Vec res = typename CVec::Dimension()(d);;
//std::cout << "Mat: " << m << "\n Vec: " << one << std::endl;
LA::solve(res, CGAL_MOVE(m), CGAL_MOVE(b));
//std::cout << "Sol: " << res << std::endl;
return cp(d,LA::vector_begin(res),LA::vector_end(res));
}
else
{
/*
* Matrix P=(p1, p2, ...) (each point as a column)
* Matrix Q=2*t(p2-p1,p3-p1, ...) (each vector as a line)
* Matrix M: QP, adding a line of 1 at the top
* Vector B: (1, p2^2-p1^2, p3^2-p1^2, ...)
* Solve ML=B, the center of the sphere is PL
*
* It would likely be faster to write P then transpose, multiply,
* etc instead of doing it by hand.
*/
// TODO: check for degenerate cases?
typedef typename R_::Max_ambient_dimension D2;
typedef typename R_::LA::template Rebind_dimension<Dynamic_dimension_tag,D2>::Other LAd;
typedef typename LAd::Square_matrix Matrix;
typedef typename LAd::Vector Vec;
typename Get_functor<R_, Scalar_product_tag>::type sp(this->kernel());
int k=static_cast<int>(std::distance(f,e));
int d=pd(p0);
Matrix m(k,k);
Vec b(k);
Vec l(k);
int j,i=0;
for(Iter f2=f;f2!=e;++f2,++i){
b(i)=m(i,i)=sdo(*f2);
j=0;
for(Iter f3=f;f3!=e;++f3,++j){
m(j,i)=m(i,j)=sp(*f2,*f3);
}
}
for(i=1;i<k;++i){
b(i)-=b(0);
for(j=0;j<k;++j){
m(i,j)=2*(m(i,j)-m(0,j));
}
}
for(j=0;j<k;++j) m(0,j)=1;
b(0)=1;
LAd::solve(l,CGAL_MOVE(m),CGAL_MOVE(b));
typename LA::Vector center=typename LA::Construct_vector::Dimension()(d);
for(i=0;i<d;++i) center(i)=0;
j=0;
for(Iter f2=f;f2!=e;++f2,++j){
for(i=0;i<d;++i){
center(i)+=l(j)*c(*f2,i);
}
}
return cp(LA::vector_begin(center),LA::vector_end(center));
}
}
};
}
CGAL_KD_DEFAULT_FUNCTOR(Construct_circumcenter_tag,(CartesianDKernelFunctors::Construct_circumcenter<K>),(Point_tag),(Construct_ttag<Point_tag>,Compute_point_cartesian_coordinate_tag,Scalar_product_tag,Squared_distance_to_origin_tag,Point_dimension_tag));
namespace CartesianDKernelFunctors {
template <class R_> struct Squared_circumradius : Store_kernel<R_> {
CGAL_FUNCTOR_INIT_STORE(Squared_circumradius)
typedef typename Get_type<R_, FT_tag>::type result_type;
template <class Iter>
result_type operator()(Iter f, Iter e)const{
typename Get_functor<R_, Construct_circumcenter_tag>::type cc(this->kernel());
typename Get_functor<R_, Squared_distance_tag>::type sd(this->kernel());
return sd(cc(f, e), *f);
}
};
}
CGAL_KD_DEFAULT_FUNCTOR(Squared_circumradius_tag,(CartesianDKernelFunctors::Squared_circumradius<K>),(Point_tag),(Construct_circumcenter_tag,Squared_distance_tag));
namespace CartesianDKernelFunctors {
// TODO: implement it directly, it should be at least as fast as Side_of_oriented_sphere.
template<class R_> struct Side_of_bounded_sphere : private Store_kernel<R_> {
@ -631,6 +748,9 @@ template<class R_> struct Side_of_bounded_sphere : private Store_kernel<R_> {
template<class Iter>
result_type operator()(Iter f, Iter const& e) const {
Point const& p0 = *f++; // *--e ?
typename Get_functor<R, Point_dimension_tag>::type pd(this->kernel());
//FIXME: Doesn't work for non-full dimension.
CGAL_assertion (std::distance(f,e) == pd(p0)+1);
return operator() (f, e, p0);
}
@ -660,6 +780,24 @@ template<class R_> struct Side_of_bounded_sphere : private Store_kernel<R_> {
CGAL_KD_DEFAULT_FUNCTOR(Side_of_bounded_sphere_tag,(CartesianDKernelFunctors::Side_of_bounded_sphere<K>),(Point_tag),(Side_of_oriented_sphere_tag,Orientation_of_points_tag));
namespace CartesianDKernelFunctors {
template<class R_> struct Side_of_bounded_circumsphere : private Store_kernel<R_> {
CGAL_FUNCTOR_INIT_STORE(Side_of_bounded_circumsphere)
typedef typename Get_type<R_, Bounded_side_tag>::type result_type;
template<class Iter, class P>
result_type operator()(Iter f, Iter const& e, P const& p0) const {
// TODO: Special case when the dimension is full.
typename Get_functor<R_, Construct_circumcenter_tag>::type cc(this->kernel());
typename Get_functor<R_, Compare_distance_tag>::type cd(this->kernel());
return enum_cast<Bounded_side>(cd(cc(f, e), *f, p0));
}
};
}
CGAL_KD_DEFAULT_FUNCTOR(Side_of_bounded_circumsphere_tag,(CartesianDKernelFunctors::Side_of_bounded_circumsphere<K>),(Point_tag),(Squared_distance_tag,Construct_circumcenter_tag));
namespace CartesianDKernelFunctors {
template<class R_> struct Point_to_vector : private Store_kernel<R_> {
CGAL_FUNCTOR_INIT_STORE(Point_to_vector)
@ -687,13 +825,13 @@ template<class R_> struct Vector_to_point : private Store_kernel<R_> {
typedef typename Get_type<R, RT_tag>::type RT;
typedef typename Get_type<R, Vector_tag>::type Vector;
typedef typename Get_type<R, Point_tag>::type Point;
typedef typename Get_functor<R, Construct_ttag<Point_tag> >::type CV;
typedef typename Get_functor<R, Construct_ttag<Point_tag> >::type CP;
typedef typename Get_functor<R, Construct_ttag<Vector_cartesian_const_iterator_tag> >::type CI;
typedef Point result_type;
typedef Vector argument_type;
result_type operator()(argument_type const&v)const{
CI ci(this->kernel());
return CV(this->kernel())(ci(v,Begin_tag()),ci(v,End_tag()));
return CP(this->kernel())(ci(v,Begin_tag()),ci(v,End_tag()));
}
};
}

3
NewKernel_d/include/CGAL/NewKernel_d/functor_tags.h
View File

@ -210,6 +210,7 @@ namespace CGAL {
CGAL_DECL_COMPUTE(Squared_distance_to_origin);
CGAL_DECL_COMPUTE(Squared_length);
CGAL_DECL_COMPUTE(Squared_radius);
CGAL_DECL_COMPUTE(Squared_circumradius);
CGAL_DECL_COMPUTE(Scalar_product);
CGAL_DECL_COMPUTE(Hyperplane_translation);
CGAL_DECL_COMPUTE(Value_at);
@ -264,6 +265,7 @@ namespace CGAL {
CGAL_DECL_CONSTRUCT(Vector_to_point,Point);
CGAL_DECL_CONSTRUCT(Construct_min_vertex,Point);
CGAL_DECL_CONSTRUCT(Construct_max_vertex,Point);
CGAL_DECL_CONSTRUCT(Construct_circumcenter,Point);
#undef CGAL_DECL_CONSTRUCT
#if 0
#define CGAL_DECL_ITER_CONSTRUCT(X,Y) struct X##_tag {}; \
@ -293,6 +295,7 @@ namespace CGAL {
CGAL_DECL_PREDICATE(Orientation_of_vectors);
CGAL_DECL_PREDICATE(Side_of_oriented_sphere);
CGAL_DECL_PREDICATE(Side_of_bounded_sphere);
CGAL_DECL_PREDICATE(Side_of_bounded_circumsphere);
CGAL_DECL_PREDICATE(Contained_in_affine_hull);
CGAL_DECL_PREDICATE(In_flat_orientation);
CGAL_DECL_PREDICATE(In_flat_side_of_oriented_sphere);

20
NewKernel_d/test/NewKernel_d/Epick_d.cpp
View File

@ -81,6 +81,7 @@ void test2(){
typedef typename K1::Construct_sphere_d CSp;
typedef typename CGAL::Get_functor<K1, CGAL::Segment_extremity_tag>::type CSE;
typedef typename K1::Construct_cartesian_const_iterator_d CCI;
typedef typename K1::Construct_circumcenter_d CCc;
typedef typename K1::Linear_base_d LB;
typedef typename K1::Orientation_d PO;
typedef typename K1::Side_of_oriented_sphere_d SOS;
@ -122,6 +123,7 @@ void test2(){
typedef typename K1::Difference_of_points_d DP;
typedef typename K1::Construct_min_vertex_d CmV;
typedef typename K1::Construct_max_vertex_d CMV;
typedef typename K1::Compute_squared_radius_d SR;
CGAL_USE_TYPE(AT);
CGAL_USE_TYPE(D);
@ -138,6 +140,7 @@ void test2(){
CV cv Kinit(construct_vector_d_object);
CCI ci Kinit(construct_cartesian_const_iterator_d_object);
CC cc Kinit(compute_coordinate_d_object);
CCc ccc Kinit(construct_circumcenter_d_object);
PO po Kinit(orientation_d_object);
CS cs Kinit(construct_segment_d_object);
CSp csp Kinit(construct_sphere_d_object);
@ -181,6 +184,7 @@ void test2(){
DP dp Kinit(difference_of_points_d_object);
CmV cmv Kinit(construct_min_vertex_d_object);
CMV cMv Kinit(construct_max_vertex_d_object);
SR sr Kinit(compute_squared_radius_d_object);
CGAL_USE(bc);
CGAL_USE(pol);
@ -229,6 +233,12 @@ void test2(){
assert(sos(tabn+0,tabn+3,P(3,3))==CGAL::ON_POSITIVE_SIDE);
assert(sbs(tabp+0,tabp+3,P(3,3))==CGAL::ON_UNBOUNDED_SIDE);
assert(sbs(tabn+0,tabn+3,P(3,3))==CGAL::ON_UNBOUNDED_SIDE);
assert(sbs(tabp+1,tabp+3,P(1,1))==CGAL::ON_BOUNDARY);
assert(ccc(tabp+1,tabp+2)==tabp[1]);
assert(ccc(tabn+0,tabn+2)==P(0,.5));
assert(sr(tabp+2,tabp+3)==0);
assert(sr(tabp+1,tabp+3)==.5);
assert(sbs(tabp+1,tabp+3,P(10,-1))==CGAL::ON_UNBOUNDED_SIDE);
assert(sos(tabp+0,tabp+3,P(.5,.5))==CGAL::ON_POSITIVE_SIDE);
assert(sos(tabn+0,tabn+3,P(.5,.5))==CGAL::ON_NEGATIVE_SIDE);
assert(sbs(tabp+0,tabp+3,P(.5,.5))==CGAL::ON_BOUNDED_SIDE);
@ -401,6 +411,7 @@ void test3(){
typedef typename K1::Point_dimension_d PD;
typedef typename K1::Affinely_independent_d AI;
typedef typename K1::Scaled_vector_d SV;
typedef typename K1::Side_of_bounded_sphere_d SBDS;
Ker k
#if 1
@ -432,6 +443,7 @@ void test3(){
SD sd Kinit(squared_distance_d_object);
PD pd Kinit(point_dimension_d_object);
AI ai Kinit(affinely_independent_d_object);
SBDS sbds Kinit(side_of_bounded_sphere_d_object);
P a; // Triangulation needs this :-(
a=cp(2,3,4);
assert(pd(a)==3);
@ -459,7 +471,7 @@ void test3(){
P tab[]={a,b,c,d,e};
std::cout << po (&tab[0],tab+4) << ' ';
std::cout << sos(&tab[0],tab+5) << ' ';
std::cout << sbs(&tab[0],tab+5) << std::endl;
std::cout << sbs(tab+1,tab+5,tab[0]) << std::endl;
FO fo=cfo(&tab[0],tab+3);
std::cout << fo;
P x[]={cp(2,2,3),cp(2,2,0),cp(1,2,1)};
@ -559,6 +571,12 @@ void test3(){
std::ostringstream output;
output << showit;
assert(output.str()=="3 1 2 4");
P t1[]={cp(1,2,3),cp(3,2,1),cp(2,4,2)};
assert(sbds(t1+0,t1+2,cp(2,2,3.414)) == CGAL::ON_BOUNDED_SIDE);
assert(sbds(t1+0,t1+2,cp(1,2,3)) == CGAL::ON_BOUNDARY);
assert(sbds(t1+0,t1+2,cp(2,2,3.415)) == CGAL::ON_UNBOUNDED_SIDE);
assert(sbds(t1+0,t1+3,cp(2.1,3.5,1.9)) == CGAL::ON_BOUNDED_SIDE);
assert(sbds(t1+0,t1+3,cp(10,10,10)) == CGAL::ON_UNBOUNDED_SIDE);
}
template struct CGAL::Epick_d<CGAL::Dimension_tag<2> >;
template struct CGAL::Epick_d<CGAL::Dimension_tag<3> >;