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

Merge pull request #4828 from mglisse/circumcenter-lowdim 

Rewrite circumcenter computations in NewKernel_d
This commit is contained in:
Laurent Rineau 2020-09-11 15:07:52 +02:00
parent 6cf203f5ca 53ed991b5d
commit dc74ea8c81
5 changed files with 179 additions and 117 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
NewKernel_d/include/CGAL/NewKernel_d/Cartesian_LA_functors.h
View File

@ -224,6 +224,7 @@ template<class R_> struct Scalar_product {
};
template<class R_> struct Squared_distance_to_origin_stored {
// What about weighted points, should they store sdo-w?
CGAL_FUNCTOR_INIT_IGNORE(Squared_distance_to_origin_stored)
typedef R_ R;
typedef typename R::LA_vector LA;

118
NewKernel_d/include/CGAL/NewKernel_d/Types/Weighted_point.h
View File

@ -149,6 +149,7 @@ template <class R_> struct Power_center : Store_kernel<R_> {
typedef typename Get_type<R_, Weighted_point_tag>::type WPoint;
typedef WPoint result_type;
typedef typename Get_type<R_, Point_tag>::type Point;
typedef typename Get_type<R_, Vector_tag>::type Vector;
typedef typename Get_type<R_, FT_tag>::type FT;
template <class Iter>
result_type operator()(Iter f, Iter e)const{
@ -168,11 +169,26 @@ template <class R_> struct Power_center : Store_kernel<R_> {
WPoint const& wp0 = *f;
Point const& p0 = pdw(wp0);
FT const& w0 = pw(wp0);
int d = pd(p0);
int k = static_cast<int>(std::distance(f,e));
if (d+1 == k)
{
FT const& n0 = sdo(p0) - pw(wp0);
if (k == 1) return cwp(p0, -w0);
// TODO: check for degenerate cases?
if (k == 2) {
typename Get_functor<R_, Difference_of_points_tag>::type dp(this->kernel());
typename Get_functor<R_, Squared_length_tag>::type sl(this->kernel());
typename Get_functor<R_, Translated_point_tag>::type tp(this->kernel());
typename Get_functor<R_, Scaled_vector_tag>::type sv(this->kernel());
WPoint const& wp1 = *++f;
Point const& p1 = pdw(wp1);
FT const& w1 = pw(wp1);
Vector v01 = dp(p1, p0);
FT l01 = sl(v01);
FT coef = ((w0 - w1) / l01 + 1) / 2;
return cwp(tp(p0, sv(v01, coef)), CGAL::square(coef) * l01 - w0);
}
if (d+1 == k) {
FT const& n0 = sdo(p0) - w0;
Matrix m(d,d);
Vec b = typename CVec::Dimension()(d);
// Write the point coordinates in lines.
@ -194,64 +210,54 @@ template <class R_> struct Power_center : Store_kernel<R_> {
FT const& r2 = pdp (wp0, center);
return cwp(std::move(center), r2);
}
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, ...) plus weights
* 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;
// The general case. ui=p(i+1)-p0, center-p0=c=sum ai*ui, c.2ui=ui²+w-w, M*a=b with M symmetric
typedef typename Increment_dimension<typename R_::Max_ambient_dimension>::type D2;
typedef typename R_::LA::template Rebind_dimension<Dynamic_dimension_tag,D2>::Other LA;
typedef typename LA::Square_matrix Matrix;
typedef typename LA::Vector Vec;
typedef typename LA::Construct_vector CVec;
typename Get_functor<R_, Translated_point_tag>::type tp(this->kernel());
typename Get_functor<R_, Scaled_vector_tag>::type sv(this->kernel());
typename Get_functor<R_, Difference_of_points_tag>::type dp(this->kernel());
typename Get_functor<R_, Scalar_product_tag>::type sp(this->kernel());
Matrix m(k,k);
Vec b(k);
Vec l(k);
int j,i=0;
for(Iter f2=f; f2!=e; ++f2,++i){
WPoint const& wp = *f2;
typename Get_functor<R_, Sum_of_vectors_tag>::type pv(this->kernel());
typename Get_functor<R_, Squared_length_tag>::type sl(this->kernel());
Matrix m(k-1,k-1);
Vec b = typename CVec::Dimension()(k-1);
std::vector<Vector> vecs; vecs.reserve(k-1);
for(int i=0; ++f!=e; ++i) {
WPoint const& wp = *f;
Point const& p = pdw(wp);
b(i) = m(i,i) = sdo(p) - pw(wp);
j=0;
for(Iter f3=f; f3!=e; ++f3,++j){
// FIXME: scalar product of points ???
m(j,i) = m(i,j) = sp(p,pdw(*f3));
}
vecs.emplace_back(dp(p, p0));
b[i] = w0 - pw(wp);
}
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));
// Only need to fill the lower half
for(int i = 0; i < k-1; ++i){
for(int j = i; j < k-1; ++j){
m(j, i) = sp(vecs[i], vecs[j]);
#if ! EIGEN_VERSION_AT_LEAST(3, 3, 5)
m(i, j) = m(j, i);
#endif
}
b[i] += m(i, i);
b[i] /= 2;
}
for(j=0;j<k;++j) m(0,j)=1;
b(0)=1;
LAd::solve(l,std::move(m),std::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){
WPoint const& wp = *f2;
Point const& p = pdw(wp);
for(i=0;i<d;++i){
center(i)+=l(j)*c(p,i);
}
}
Point c = cp(d, LA::vector_begin(center), LA::vector_end(center));
FT r2 = pdp (wp0, c);
return cwp(std::move(c), std::move(r2));
// Assumes Eigen...
#if EIGEN_VERSION_AT_LEAST(3, 3, 5)
Vec res = m.ldlt().solve(b);
#else
// Older versions of Eigen use 1/highest as tolerance,
// which we have no way to set to 0 for exact types.
// Use something slow but that should work.
Vec res = m.fullPivLu().solve(b);
#endif
Vector to_center = sv(vecs[0], res[0]);
for(int i=1;i<k-1;++i)
to_center = pv(to_center, sv(vecs[i],res[i]));
return cwp(tp(p0, to_center), sl(to_center) - w0);
}
}
};
@ -279,7 +285,7 @@ CGAL_KD_DEFAULT_FUNCTOR(Power_side_of_power_sphere_tag,(CartesianDKernelFunctors
CGAL_KD_DEFAULT_FUNCTOR(In_flat_power_side_of_power_sphere_tag,(CartesianDKernelFunctors::In_flat_power_side_of_power_sphere<K>),(Weighted_point_tag),(In_flat_power_side_of_power_sphere_raw_tag,Point_drop_weight_tag,Point_weight_tag));
CGAL_KD_DEFAULT_FUNCTOR(Power_distance_tag,(CartesianDKernelFunctors::Power_distance<K>),(Weighted_point_tag,Point_tag),(Squared_distance_tag,Point_drop_weight_tag,Point_weight_tag));
CGAL_KD_DEFAULT_FUNCTOR(Power_distance_to_point_tag,(CartesianDKernelFunctors::Power_distance_to_point<K>),(Weighted_point_tag,Point_tag),(Squared_distance_tag,Point_drop_weight_tag,Point_weight_tag));
CGAL_KD_DEFAULT_FUNCTOR(Power_center_tag,(CartesianDKernelFunctors::Power_center<K>),(Weighted_point_tag,Point_tag),(Compute_point_cartesian_coordinate_tag,Construct_ttag<Point_tag>,Construct_ttag<Weighted_point_tag>,Point_dimension_tag,Squared_distance_to_origin_tag,Point_drop_weight_tag,Point_weight_tag,Power_distance_to_point_tag));
CGAL_KD_DEFAULT_FUNCTOR(Power_center_tag,(CartesianDKernelFunctors::Power_center<K>),(Weighted_point_tag,Point_tag,Vector_tag),(Compute_point_cartesian_coordinate_tag,Construct_ttag<Point_tag>,Construct_ttag<Weighted_point_tag>,Point_dimension_tag,Squared_distance_to_origin_tag,Point_drop_weight_tag,Point_weight_tag,Power_distance_to_point_tag,Translated_point_tag,Scaled_vector_tag,Difference_of_points_tag,Scalar_product_tag,Sum_of_vectors_tag,Squared_length_tag));
CGAL_KD_DEFAULT_FUNCTOR(Power_side_of_bounded_power_circumsphere_tag,(CartesianDKernelFunctors::Power_side_of_bounded_power_circumsphere<K>),(Weighted_point_tag),(Power_distance_tag,Power_center_tag));
} // namespace CGAL
#endif

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

@ -27,6 +27,7 @@
#include <CGAL/predicates/sign_of_determinant.h>
#include <functional>
#include <initializer_list>
#include <vector>
namespace CGAL {
namespace CartesianDKernelFunctors {
@ -603,25 +604,57 @@ 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 typename Get_type<R_, Vector_tag>::type Vector;
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))
int k = static_cast<int>(std::distance(f,e));
// Sorted from fastest to slowest, whether the dimension is 2, 3 or 4. It may be worth checking at some point.
CGAL_assume(k>=1);
if(k==1) return p0;
if(k==2){
typename Get_functor<R_, Midpoint_tag>::type mid(this->kernel());
return mid(p0, *++f);
}
// TODO: check for degenerate cases in all the following cases?
if(k==3){
// Same equations as in the general case, but solved by hand (Cramer)
// (c-r).(p-r)=(p-r)²/2
// (c-r).(q-r)=(q-r)²/2
typename Get_functor<R_, Squared_length_tag>::type sl(this->kernel());
typename Get_functor<R_, Scalar_product_tag>::type sp(this->kernel());
typename Get_functor<R_, Scaled_vector_tag>::type sv(this->kernel());
typename Get_functor<R_, Difference_of_points_tag>::type dp(this->kernel());
typename Get_functor<R_, Translated_point_tag>::type tp(this->kernel());
Iter f2=f;
Point const& q=*++f2;
Point const& r=*++f2;
Vector u = dp(p0, r);
Vector v = dp(q, r);
FT uv = sp(u, v);
FT u2 = sl(u);
FT v2 = sl(v);
FT den = 2 * (u2 * v2 - CGAL::square(uv));
FT a = (u2 - uv) * v2 / den;
FT b = (v2 - uv) * u2 / den;
// Wasteful if we only want the radius
return tp(tp(r, sv(u, a)), sv(v, b));
}
if (k == d+1)
{
// 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;
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_, Squared_distance_to_origin_tag>::type sdo(this->kernel());
FT const& n0 = sdo(p0);
Matrix m(d,d);
Vec b = typename CVec::Dimension()(d);
@ -635,72 +668,59 @@ template <class R_> struct Construct_circumcenter : Store_kernel<R_> {
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, std::move(m), std::move(b));
//std::cout << "Sol: " << res << std::endl;
//Vec res = typename CVec::Dimension()(d);;
//LA::solve(res, std::move(m), std::move(b));
// We already assume Eigen below...
Vec res=m.partialPivLu().solve(b);
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;
// The general case. ui=p(i+1)-p0, center-p0=c=sum ai*ui, c.ui=ui²/2, M*a=b with M symmetric
typedef typename Increment_dimension<typename R_::Max_ambient_dimension>::type D2;
typedef typename R_::LA::template Rebind_dimension<Dynamic_dimension_tag,D2>::Other LA;
typedef typename LA::Square_matrix Matrix;
typedef typename LA::Vector Vec;
typedef typename LA::Construct_vector CVec;
typename Get_functor<R_, Translated_point_tag>::type tp(this->kernel());
typename Get_functor<R_, Scaled_vector_tag>::type sv(this->kernel());
typename Get_functor<R_, Difference_of_points_tag>::type dp(this->kernel());
typename Get_functor<R_, Scalar_product_tag>::type sp(this->kernel());
int k=static_cast<int>(std::distance(f,e));
Matrix m(k,k);
Vec b(k);
Vec l(k);
int j,i=0;
// We are doing a quadratic number of *f, which can be costly with transforming_iterator.
for(Iter f2=f;f2!=e;++f2,++i){
Point const& p2=*f2;
b(i)=m(i,i)=sdo(p2);
j=0;
for(Iter f3=f;f3!=e;++f3,++j){
m(j,i)=m(i,j)=sp(p2,*f3);
Matrix m(k-1,k-1);
Vec b = typename CVec::Dimension()(k-1);
std::vector<Vector> vecs; vecs.reserve(k-1);
while(++f!=e)
vecs.emplace_back(dp(*f,p0));
// Only need to fill the lower half
for(int i=0;i<k-1;++i){
for(int j=i;j<k-1;++j) {
m(j,i)=sp(vecs[i],vecs[j]);
#if ! EIGEN_VERSION_AT_LEAST(3, 3, 5)
m(i,j)=m(j,i);
#endif
}
b[i]=m(i,i)/2;
}
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,std::move(m),std::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));
// Assumes Eigen...
#if EIGEN_VERSION_AT_LEAST(3, 3, 5)
Vec res=m.ldlt().solve(b);
#else
// Older versions of Eigen use 1/highest as tolerance,
// which we have no way to set to 0 for exact types.
// Use something slow but that should work.
Vec res=m.fullPivLu().solve(b);
#endif
Point center=p0;
// Wasteful if we only want the radius
for(int i=0;i<k-1;++i)
center=tp(center,sv(vecs[i],res[i]));
return 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));
CGAL_KD_DEFAULT_FUNCTOR(Construct_circumcenter_tag,(CartesianDKernelFunctors::Construct_circumcenter<K>),(Point_tag,Vector_tag),(Construct_ttag<Point_tag>,Compute_point_cartesian_coordinate_tag,Scalar_product_tag,Squared_distance_to_origin_tag,Point_dimension_tag,Translated_point_tag,Scaled_vector_tag,Difference_of_points_tag,Squared_length_tag));
namespace CartesianDKernelFunctors {
template <class R_> struct Squared_circumradius : Store_kernel<R_> {

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

@ -712,6 +712,38 @@ void test3(){
std::ostringstream sv1; sv1 << v1; assert(sv1.str()=="3 3 2 1");
std::istringstream sv2("3 4 5 6"); sv2 >> v1; assert(v1[0]==4&&v1[1]==5);
}
template<class Ker>
void test4(){
typedef typename Ker::Point_d P;
typedef typename Ker::Weighted_point_d WP;
typedef typename Ker::Construct_circumcenter_d CCc;
typedef typename Ker::Equal_d E;
typedef typename Ker::Power_center_d PC;
typedef typename Ker::Power_distance_d PoD;
typedef typename Ker::Affine_rank_d AR;
Ker k(4);
CCc ccc Kinit(construct_circumcenter_d_object);
E ed Kinit(equal_d_object);
PC pc Kinit(power_center_d_object);
PoD pod Kinit(power_distance_d_object);
AR ar Kinit(affine_rank_d_object);
auto mkpt=[](auto...x){double l[]{(double)x...};return P(std::begin(l), std::end(l));};
P tab1[]={mkpt(15,20,40,80),mkpt(10,23,36,80),mkpt(10,20,40,85),mkpt(10,15,40,80),mkpt(13,20,40,76)};
assert(ed(ccc(tab1+0, tab1+5),mkpt(10,20,40,80)));
P tab2[]={mkpt(15,20,40,80),mkpt(13,24,40,80),mkpt(10,25,40,80),mkpt(10,20,43,84)};
assert(ed(ccc(tab2+0, tab2+4),mkpt(10,20,40,80)));
P tab3[]={mkpt(15,20,35,80),mkpt(10,25,40,75),mkpt(13,24,37,76)};
assert(ed(ccc(tab3+0, tab3+3),mkpt(10,20,40,80)));
auto mkwpt=[](auto...x){double l[]{(double)x...};auto last=std::prev(std::end(l));return WP(P(std::begin(l), last),*last);};
WP tab4[]={mkwpt(89,17,29,97,14),mkwpt(86,99,64,26,44),mkwpt(40,9,13,91,20),mkwpt(41,30,93,13,10),mkwpt(45,6,98,9,0),mkwpt(0,0,0,0,0)};
for(int i=5;i>=1;--i){
tab4[i]=pc(tab4+0, tab4+i);
for(int j=0;j<i;++j)
assert(pod(tab4[i],tab4[j])==0);
auto drop=[](WP const&x){return x.point();};
assert(ar(CGAL::make_transforming_iterator(tab4+0,drop), CGAL::make_transforming_iterator(tab4+i+1,drop))==i-1);
}
}
template struct CGAL::Epick_d<CGAL::Dimension_tag<2> >;
template struct CGAL::Epick_d<CGAL::Dimension_tag<3> >;
template struct CGAL::Epick_d<CGAL::Dynamic_dimension_tag>;
@ -733,6 +765,7 @@ int main(){
test2<CGAL::Epeck_d<CGAL::Dimension_tag<2>>>();
test3<CGAL::Epeck_d<CGAL::Dimension_tag<3>>>();
test3<CGAL::Epeck_d<CGAL::Dynamic_dimension_tag>>();
test4<CGAL::Epeck_d<CGAL::Dynamic_dimension_tag>>();
#endif
}

2
Number_types/include/CGAL/Interval_nt.h
View File

@ -1573,6 +1573,8 @@ namespace Eigen {
static inline Real epsilon() { return 0; }
static inline Real dummy_precision() { return 0; }
static inline Real highest() { return Real((std::numeric_limits<double>::max)(), std::numeric_limits<double>::infinity()); }
static inline Real lowest() { return Real(-std::numeric_limits<double>::infinity(), std::numeric_limits<double>::lowest()); }
// Costs could depend on b.
enum {