Merge pull request #4082 from mglisse/NewKernel_d-alphaweight-glisse

Add functors for weighted alpha-complexes

Installation/CHANGES.md

@@ -51,6 +51,11 @@ Release date: June 2020
         k-NN search to interrupt some distance computations before its end,
         saving precious milliseconds, in particular in medium-to-high dimension.
 
+### dD Geometry Kernel
+-   Epick\_d and Epeck\_d gain 2 new functors: `Power_side_of_bounded_power_sphere_d` and
+    `Compute_squared_radius_smallest_orthogonal_sphere_d`. Those are
+    essential for the computation of weighted alpha-complexes.
+
 [Release 5.0](https://github.com/CGAL/cgal/releases/tag/releases%2FCGAL-5.0)
 -----------
 

Kernel_d/doc/Kernel_d/CGAL/Epeck_d.h

@@ -120,6 +120,15 @@ public:
 template<class ForwardIterator>
 FT operator()(ForwardIterator first, ForwardIterator last);
 };
+class Compute_squared_radius_smallest_orthogonal_sphere_d {
+public:
+/*! returns the radius of the sphere defined by `A=tuple[first,last)`. The sphere is centered in the affine hull of A and orthogonal to all the spheres of A. The order of the points of A does not matter.
+    \pre A is affinely independent.
+    \tparam ForwardIterator has `Epeck_d::Weighted_point_d` as value type.
+    */
+template<class ForwardIterator>
+FT operator()(ForwardIterator first, ForwardIterator last);
+};
 /*! \cgalModels `Kernel_d::Side_of_bounded_sphere_d`
  */
 class Side_of_bounded_sphere_d {
@@ -131,7 +140,18 @@ public:
 template<class ForwardIterator>
 Bounded_side operator()(ForwardIterator first, ForwardIterator last, const Point_d&p);
 };
+class Power_side_of_bounded_power_sphere_d {
+public:
+/*! returns the relative position of weighted point p to the sphere defined by `A=tuple[first,last)`. The sphere is centered in the affine hull of A and orthogonal to all the spheres of A. The order of the points of A does not matter.
+    \pre A is affinely independent.
+    \tparam ForwardIterator has `Epeck_d::Weighted_point_d` as value type.
+    */
+template<class ForwardIterator>
+Bounded_side operator()(ForwardIterator first, ForwardIterator last, const Weighted_point_d&p);
+};
 Construct_circumcenter_d construct_circumcenter_d_object();
 Compute_squared_radius_d compute_squared_radius_d_object();
+Compute_squared_radius_smallest_orthogonal_sphere_d compute_squared_radius_smallest_orthogonal_sphere_d_object();
+Power_side_of_bounded_power_sphere_d power_side_of_bounded_power_sphere_d_object();
 }; /* end Epeck_d */
 } /* end namespace CGAL */

Kernel_d/doc/Kernel_d/CGAL/Epick_d.h

@@ -109,6 +109,15 @@ public:
 template<class ForwardIterator>
 FT operator()(ForwardIterator first, ForwardIterator last);
 };
+class Compute_squared_radius_smallest_orthogonal_sphere_d {
+public:
+/*! returns the radius of the sphere defined by `A=tuple[first,last)`. The sphere is centered in the affine hull of A and orthogonal to all the spheres of A. The order of the points of A does not matter.
+    \pre A is affinely independent.
+    \tparam ForwardIterator has `Epick_d::Weighted_point_d` as value type.
+    */
+template<class ForwardIterator>
+FT operator()(ForwardIterator first, ForwardIterator last);
+};
 /*! \cgalModels `Kernel_d::Side_of_bounded_sphere_d`
  */
 class Side_of_bounded_sphere_d {
@@ -120,7 +129,18 @@ public:
 template<class ForwardIterator>
 Bounded_side operator()(ForwardIterator first, ForwardIterator last, const Point_d&p);
 };
+class Power_side_of_bounded_power_sphere_d {
+public:
+/*! returns the relative position of weighted point p to the sphere defined by `A=tuple[first,last)`. The sphere is centered in the affine hull of A and orthogonal to all the spheres of A. The order of the points of A does not matter.
+    \pre A is affinely independent.
+    \tparam ForwardIterator has `Epick_d::Weighted_point_d` as value type.
+    */
+template<class ForwardIterator>
+Bounded_side operator()(ForwardIterator first, ForwardIterator last, const Weighted_point_d&p);
+};
 Construct_circumcenter_d construct_circumcenter_d_object();
 Compute_squared_radius_d compute_squared_radius_d_object();
+Compute_squared_radius_smallest_orthogonal_sphere_d compute_squared_radius_smallest_orthogonal_sphere_d_object();
+Power_side_of_bounded_power_sphere_d power_side_of_bounded_power_sphere_d_object();
 }; /* end Epick_d */
 } /* end namespace CGAL */

NewKernel_d/include/CGAL/NewKernel_d/Kernel_d_interface.h

@@ -58,6 +58,7 @@ template <class Base_> struct Kernel_d_interface : public Base_ {
 	typedef typename Get_functor<Base, Point_dimension_tag>::type Point_dimension_d;
 	typedef typename Get_functor<Base, Side_of_oriented_sphere_tag>::type Side_of_oriented_sphere_d;
 	typedef typename Get_functor<Base, Power_side_of_power_sphere_tag>::type Power_side_of_power_sphere_d;
+        typedef typename Get_functor<Base, Power_side_of_bounded_power_circumsphere_tag>::type Power_side_of_bounded_power_sphere_d;
 	typedef typename Get_functor<Base, Power_center_tag>::type Power_center_d;
 	typedef typename Get_functor<Base, Power_distance_tag>::type Power_distance_d;
 	typedef typename Get_functor<Base, Contained_in_affine_hull_tag>::type Contained_in_affine_hull_d;
@@ -166,6 +167,16 @@ template <class Base_> struct Kernel_d_interface : public Base_ {
 	    return typename Get_functor<Base, Squared_circumradius_tag>::type(this->kernel())(b,e);
 	  }
 	};
+        struct Compute_squared_radius_smallest_orthogonal_sphere_d : private Store_kernel<Kernel> {
+          typedef Kernel R_; // for the macro
+          CGAL_FUNCTOR_INIT_STORE(Compute_squared_radius_smallest_orthogonal_sphere_d)
+          typedef FT result_type;
+          template<class I> FT operator()(I b, I e)const{
+            typename Get_functor<Base, Point_weight_tag>::type pw(this->kernel());
+            typename Get_functor<Base, Power_center_tag>::type pc(this->kernel());
+            return pw(pc(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;
@@ -210,6 +221,7 @@ template <class Base_> struct Kernel_d_interface : public Base_ {
 	Point_of_sphere_d point_of_sphere_d_object()const{ return Point_of_sphere_d(*this); }
 	Side_of_oriented_sphere_d side_of_oriented_sphere_d_object()const{ return Side_of_oriented_sphere_d(*this); }
 	Power_side_of_power_sphere_d power_side_of_power_sphere_d_object()const{ return Power_side_of_power_sphere_d(*this); }
+        Power_side_of_bounded_power_sphere_d power_side_of_bounded_power_sphere_d_object()const{ return Power_side_of_bounded_power_sphere_d(*this); }
 	Power_center_d power_center_d_object()const{ return Power_center_d(*this); }
 	Power_distance_d power_distance_d_object()const{ return Power_distance_d(*this); }
 	Side_of_bounded_sphere_d side_of_bounded_sphere_d_object()const{ return Side_of_bounded_sphere_d(*this); }
@@ -243,6 +255,7 @@ template <class Base_> struct Kernel_d_interface : public Base_ {
 	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); }
+        Compute_squared_radius_smallest_orthogonal_sphere_d compute_squared_radius_smallest_orthogonal_sphere_d_object()const{ return Compute_squared_radius_smallest_orthogonal_sphere_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); }

NewKernel_d/include/CGAL/NewKernel_d/Types/Weighted_point.h

@@ -169,29 +169,107 @@ template <class R_> struct Power_center : Store_kernel<R_> {
     WPoint const& wp0 = *f;
     Point const& p0 = pdw(wp0);
     int d = pd(p0);
-    FT const& n0 = sdo(p0) - pw(wp0);
-    Matrix m(d,d);
-    Vec b = typename CVec::Dimension()(d);
-    // Write the point coordinates in lines.
-    int i;
-    for(i=0; ++f!=e; ++i) {
-      WPoint const& wp=*f;
-      Point const& p=pdw(wp);
-      for(int j=0;j<d;++j) {
-	m(i,j)=2*(c(p,j)-c(p0,j));
+    int k = static_cast<int>(std::distance(f,e));
+    if (d+1 == k)
+    {
+      FT const& n0 = sdo(p0) - pw(wp0);
+      Matrix m(d,d);
+      Vec b = typename CVec::Dimension()(d);
+      // Write the point coordinates in lines.
+      int i;
+      for(i=0; ++f!=e; ++i) {
+        WPoint const& wp=*f;
+        Point const& p=pdw(wp);
+        for(int j=0;j<d;++j) {
+          m(i,j)=2*(c(p,j)-c(p0,j));
+        }
+        b[i] = sdo(p) - pw(wp) - n0;
       }
-      b[i] = sdo(p) - pw(wp) - 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;
+      Point center = cp(d,LA::vector_begin(res),LA::vector_end(res));
+      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;
+      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;
+        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));
+        }
+      }
+      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){
+        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));
     }
-    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;
-    Point center = cp(d,LA::vector_begin(res),LA::vector_end(res));
-    FT const& r2 = pdp (wp0, center);
-    return cwp(std::move(center), r2);
   }
 };
+template<class R_> struct Power_side_of_bounded_power_circumsphere : private Store_kernel<R_> {
+        CGAL_FUNCTOR_INIT_STORE(Power_side_of_bounded_power_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_, Power_center_tag>::type pc(this->kernel());
+          typename Get_functor<R_, Power_distance_tag>::type pd(this->kernel());
+
+          // ON_UNBOUNDED_SIDE = -1
+          return enum_cast<Bounded_side>(-CGAL::sign(pd(pc(f, e), p0)));
+        }
+};
+
 }
 CGAL_KD_DEFAULT_TYPE(Weighted_point_tag,(CGAL::KerD::Weighted_point<K>),(Point_tag),());
 CGAL_KD_DEFAULT_FUNCTOR(Construct_ttag<Weighted_point_tag>,(CartesianDKernelFunctors::Construct_weighted_point<K>),(Weighted_point_tag,Point_tag),());
@@ -202,5 +280,6 @@ CGAL_KD_DEFAULT_FUNCTOR(In_flat_power_side_of_power_sphere_tag,(CartesianDKernel
 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_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

NewKernel_d/include/CGAL/NewKernel_d/functor_tags.h

@@ -303,6 +303,7 @@ namespace CGAL {
 	CGAL_DECL_PREDICATE(Contained_in_simplex);
 	CGAL_DECL_PREDICATE(Power_side_of_power_sphere_raw);
 	CGAL_DECL_PREDICATE(Power_side_of_power_sphere);
+        CGAL_DECL_PREDICATE(Power_side_of_bounded_power_circumsphere);
 	CGAL_DECL_PREDICATE(In_flat_power_side_of_power_sphere_raw);
 	CGAL_DECL_PREDICATE(In_flat_power_side_of_power_sphere);
 #undef CGAL_DECL_PREDICATE

NewKernel_d/test/NewKernel_d/Epick_d.cpp

@@ -143,8 +143,9 @@ void test2(){
   typedef typename K1::Translated_point_d TP;
   typedef typename K1::Power_center_d PC;
   typedef typename K1::Power_distance_d PoD;
-  typedef typename K1::Weighted_point_d WP;
   typedef typename K1::Construct_weighted_point_d CWP;
+  typedef typename K1::Power_side_of_bounded_power_sphere_d PSBPS;
+  typedef typename K1::Compute_squared_radius_smallest_orthogonal_sphere_d CSRSOS;
   //typedef typename K1::Point_drop_weight_d PDW;
   typedef CP PDW;
   typedef typename K1::Compute_weight_d PW;
@@ -215,6 +216,8 @@ void test2(){
   PDW const& pdw = cp;
   PW pw Kinit(compute_weight_d_object);
   PoD pod Kinit(power_distance_d_object);
+  PSBPS psbps Kinit(power_side_of_bounded_power_sphere_d_object);
+  CSRSOS csrsos Kinit(compute_squared_radius_smallest_orthogonal_sphere_d_object);
 
   CGAL_USE(bc);
   CGAL_USE(pol);
@@ -378,7 +381,13 @@ void test2(){
   assert(pdw(xw)[0]>2.5);
   assert(pw(xw)<2.95);
   assert(pw(xw)>2.5);
-
+  assert(psbps(tw+0,tw+3,cwp(cp(5,0),1.499)) == CGAL::ON_UNBOUNDED_SIDE);
+  assert(psbps(tw+0,tw+3,cwp(cp(5,0),1.500)) == CGAL::ON_BOUNDARY);
+  assert(psbps(tw+0,tw+3,cwp(cp(5,0),1.501)) == CGAL::ON_BOUNDED_SIDE);
+  WP tw2[]={cwp(cp(-3,2),1),cwp(cp(5,2),2),cwp(cp(1,6),1)};
+  assert(psbps(tw2+0,tw2+2,tw2[2]) == CGAL::ON_UNBOUNDED_SIDE);
+  assert(abs(csrsos(tw2+0,tw2+2)-14.5039)<.0001);
+  assert(abs(csrsos(tw2+0,tw2+1)+1)<.0001);
 
   P tl=cp(2,5);
   P br=cp(4,-1);