diff --git a/NewKernel_d/include/CGAL/NewKernel_d/Cartesian_LA_functors.h b/NewKernel_d/include/CGAL/NewKernel_d/Cartesian_LA_functors.h index a379cee7ed9..b717a80bac4 100644 --- a/NewKernel_d/include/CGAL/NewKernel_d/Cartesian_LA_functors.h +++ b/NewKernel_d/include/CGAL/NewKernel_d/Cartesian_LA_functors.h @@ -224,6 +224,7 @@ template struct Scalar_product { }; template 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; diff --git a/NewKernel_d/include/CGAL/NewKernel_d/Types/Weighted_point.h b/NewKernel_d/include/CGAL/NewKernel_d/Types/Weighted_point.h index 8dfe875b8c8..6e970a5f13a 100644 --- a/NewKernel_d/include/CGAL/NewKernel_d/Types/Weighted_point.h +++ b/NewKernel_d/include/CGAL/NewKernel_d/Types/Weighted_point.h @@ -149,6 +149,7 @@ template struct Power_center : Store_kernel { typedef typename Get_type::type WPoint; typedef WPoint result_type; typedef typename Get_type::type Point; + typedef typename Get_type::type Vector; typedef typename Get_type::type FT; template result_type operator()(Iter f, Iter e)const{ @@ -168,11 +169,26 @@ template struct Power_center : Store_kernel { WPoint const& wp0 = *f; Point const& p0 = pdw(wp0); + FT const& w0 = pw(wp0); int d = pd(p0); int k = static_cast(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::type dp(this->kernel()); + typename Get_functor::type sl(this->kernel()); + typename Get_functor::type tp(this->kernel()); + typename Get_functor::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 struct Power_center : Store_kernel { 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::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::type D2; + typedef typename R_::LA::template Rebind_dimension::Other LA; + typedef typename LA::Square_matrix Matrix; + typedef typename LA::Vector Vec; + typedef typename LA::Construct_vector CVec; + typename Get_functor::type tp(this->kernel()); + typename Get_functor::type sv(this->kernel()); + typename Get_functor::type dp(this->kernel()); typename Get_functor::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::type pv(this->kernel()); + typename Get_functor::type sl(this->kernel()); + + Matrix m(k-1,k-1); + Vec b = typename CVec::Dimension()(k-1); + std::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),(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),(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),(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),(Weighted_point_tag,Point_tag),(Compute_point_cartesian_coordinate_tag,Construct_ttag,Construct_ttag,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),(Weighted_point_tag,Point_tag,Vector_tag),(Compute_point_cartesian_coordinate_tag,Construct_ttag,Construct_ttag,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),(Weighted_point_tag),(Power_distance_tag,Power_center_tag)); } // namespace CGAL #endif diff --git a/NewKernel_d/include/CGAL/NewKernel_d/function_objects_cartesian.h b/NewKernel_d/include/CGAL/NewKernel_d/function_objects_cartesian.h index 5f534070db2..3ab117131dc 100644 --- a/NewKernel_d/include/CGAL/NewKernel_d/function_objects_cartesian.h +++ b/NewKernel_d/include/CGAL/NewKernel_d/function_objects_cartesian.h @@ -27,6 +27,7 @@ #include #include #include +#include namespace CGAL { namespace CartesianDKernelFunctors { @@ -603,25 +604,57 @@ namespace CartesianDKernelFunctors { template struct Construct_circumcenter : Store_kernel { CGAL_FUNCTOR_INIT_STORE(Construct_circumcenter) typedef typename Get_type::type Point; + typedef typename Get_type::type Vector; typedef Point result_type; typedef typename Get_type::type FT; template result_type operator()(Iter f, Iter e)const{ - typedef typename Get_type::type Point; - typedef typename R_::LA LA; - typename Get_functor::type c(this->kernel()); - typename Get_functor >::type cp(this->kernel()); typename Get_functor::type pd(this->kernel()); - typename Get_functor::type sdo(this->kernel()); Point const& p0=*f; int d = pd(p0); - if (d+1 == std::distance(f,e)) + int k = static_cast(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::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::type sl(this->kernel()); + typename Get_functor::type sp(this->kernel()); + typename Get_functor::type sv(this->kernel()); + typename Get_functor::type dp(this->kernel()); + typename Get_functor::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::type c(this->kernel()); + typename Get_functor >::type cp(this->kernel()); + typename Get_functor::type sdo(this->kernel()); FT const& n0 = sdo(p0); Matrix m(d,d); Vec b = typename CVec::Dimension()(d); @@ -635,72 +668,59 @@ template struct Construct_circumcenter : Store_kernel { 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::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::type D2; + typedef typename R_::LA::template Rebind_dimension::Other LA; + typedef typename LA::Square_matrix Matrix; + typedef typename LA::Vector Vec; + typedef typename LA::Construct_vector CVec; + typename Get_functor::type tp(this->kernel()); + typename Get_functor::type sv(this->kernel()); + typename Get_functor::type dp(this->kernel()); typename Get_functor::type sp(this->kernel()); - int k=static_cast(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 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),(Point_tag),(Construct_ttag,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),(Point_tag,Vector_tag),(Construct_ttag,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 struct Squared_circumradius : Store_kernel { diff --git a/NewKernel_d/test/NewKernel_d/Epick_d.cpp b/NewKernel_d/test/NewKernel_d/Epick_d.cpp index 21c64288380..3a3b150432f 100644 --- a/NewKernel_d/test/NewKernel_d/Epick_d.cpp +++ b/NewKernel_d/test/NewKernel_d/Epick_d.cpp @@ -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 +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 >; template struct CGAL::Epick_d >; template struct CGAL::Epick_d; @@ -733,6 +765,7 @@ int main(){ test2>>(); test3>>(); test3>(); + test4>(); #endif } diff --git a/Number_types/include/CGAL/Interval_nt.h b/Number_types/include/CGAL/Interval_nt.h index ca77127d6e0..eb09e9ae3fd 100644 --- a/Number_types/include/CGAL/Interval_nt.h +++ b/Number_types/include/CGAL/Interval_nt.h @@ -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::max)(), std::numeric_limits::infinity()); } + static inline Real lowest() { return Real(-std::numeric_limits::infinity(), std::numeric_limits::lowest()); } // Costs could depend on b. enum {