Merge pull request #2425 from MaelRL/Interpolation-Fix_gradient_fitting-GF

Installation/changes.html

@@ -325,6 +325,15 @@ and <code>src/</code> directories).
       refines <code>DelaunayTriangulationTraits_2</code>.
     </li>
   </ul>
+  <h3> Interpolation </h3>
+  <ul>
+    <li><b>Breaking change</b>: The concept <code>GradientFittingTraits</code>
+    now additionally requests a weighted point type <code>Weighted_point_d</code>
+    and a functor <code>Construct_point_d</code>. The model
+    <code>CGAL::Interpolation_gradient_fitting_traits_2</code> has been appropriately
+    modified to still be a model of the concept <code>GradientFittingTraits</code>.
+    </li>
+  </ul>
   <h3> 2D and 3D Triangulations </h3>
   <ul>
     <li><b>Breaking change</b>: Added a new functor

Interpolation/doc/Interpolation/Concepts/GradientFittingTraits.h

@@ -36,6 +36,11 @@ which the function is defined and interpolated.
 */ 
 typedef unspecified_type Point_d; 
 
+/*!
+The weighted point type.
+*/
+typedef unspecified_type Weighted_point_d;
+
 /*!
 The corresponding vector type. 
 */ 
@@ -55,6 +60,16 @@ multiplication of `tr` with `v`.
 */ 
 typedef unspecified_type Aff_transformation_d; 
 
+/*!
+A constructor object for `Point_d`.
+Provides :
+
+`Point_d operator() (Point_d p)` which simply returns `p`
+
+`Point_d operator() (Weighted_point wp)` which returns the bare point contained in `wp`.
+*/
+typedef unspecified_type Construct_point_d;
+
 /*!
 A constructor object for 
 `Vector_d`. 

Interpolation/include/CGAL/Interpolation_gradient_fitting_traits_2.h

@@ -100,8 +100,10 @@ public:
 
   typedef typename Rep::FT                           FT;
   typedef typename Rep::Point_2                      Point_d;
+  typedef typename Rep::Weighted_point_2             Weighted_point_d;
   typedef typename Rep::Vector_2                     Vector_d;
 
+  typedef typename Rep::Construct_point_2            Construct_point_d;
   typedef typename Rep::Construct_vector_2           Construct_vector_d;
   typedef typename Rep::Construct_scaled_vector_2    Construct_scaled_vector_d;
   //only one not needed by gradient fitting:
@@ -140,6 +142,10 @@ public:
   construct_scaled_vector_d_object()const
   {return Construct_scaled_vector_d();}
 
+  Construct_point_d
+  construct_point_d_object()const
+  {return Construct_point_d();}
+
   Construct_vector_d
   construct_vector_d_object()const
   {return Construct_vector_d();}

Interpolation/include/CGAL/sibson_gradient_fitting.h

@@ -48,7 +48,9 @@ sibson_gradient_fitting(ForwardIterator first, ForwardIterator beyond,
   typedef typename Traits::Aff_transformation_d Aff_transformation;
   typedef typename Traits::FT                   Coord_type;
 
-  typename Functor::result_type fn =  function_value(p);
+  const typename Traits::Point_d& bare_p = traits.construct_point_d_object()(p);
+
+  typename Functor::result_type fn = function_value(bare_p);
   CGAL_assertion(fn.second); //function value of p is valid
 
   typename Traits::Vector_d pn =
@@ -56,16 +58,17 @@ sibson_gradient_fitting(ForwardIterator first, ForwardIterator beyond,
   Aff_transformation scaling, m,
       Hn(traits.construct_null_matrix_d_object()());
 
-  for(;first!=beyond; ++first){
+  for(;first!=beyond; ++first)
-    Coord_type square_dist = traits.compute_squared_distance_d_object()
+  {
-                             (first->first, p);
+    const typename Traits::Point_d& bare_f = traits.construct_point_d_object()(first->first);
+
+    Coord_type square_dist = traits.compute_squared_distance_d_object()(bare_f, bare_p);
     CGAL_assertion(square_dist != 0);
     Coord_type scale = first->second / (norm*square_dist);
-    typename Traits::Vector_d d =
+    typename Traits::Vector_d d = traits.construct_vector_d_object()(bare_p, bare_f);
-        traits.construct_vector_d_object()(p, first->first);
 
     //compute the vector pn:
-    typename Functor::result_type f = function_value(first->first);
+    typename Functor::result_type f = function_value(bare_f);
     CGAL_assertion(f.second);//function value of first->first is valid
     pn = pn + traits.construct_scaled_vector_d_object()
          (d,scale * (f.first - fn.first));
@@ -89,7 +92,7 @@ sibson_gradient_fitting(const Triangul& tr,
                         CoordFunctor compute_coordinates,
                         const Traits& traits)
 {
-  typedef typename Traits::Point_d                        Point;
+  typedef typename Triangul::Point                        Point;
   typedef typename Traits::FT                             Coord_type;
 
   std::vector< std::pair< Point, Coord_type > > coords;
@@ -101,7 +104,7 @@ sibson_gradient_fitting(const Triangul& tr,
     if (!tr.is_edge(vit, tr.infinite_vertex()))
     {
       norm = compute_coordinates(tr, vit, std::back_inserter(coords)).second;
-      *out++ = std::make_pair(vit->point(),
+      *out++ = std::make_pair(traits.construct_point_d_object()(vit->point()),
                               sibson_gradient_fitting(coords.begin(),
                                                       coords.end(),
                                                       norm, vit->point(),
@@ -145,8 +148,8 @@ sibson_gradient_fitting_rn_2(const Rt& rt,
 {
   typedef typename std::back_insert_iterator<
                      std::vector<
-                       std::pair< typename Traits::Point_d,
+                       std::pair< typename Traits::Weighted_point_d,
-                                  typename  Traits::FT > > >   CoordInserter;
+                                  typename Traits::FT > > >   CoordInserter;
 
   return sibson_gradient_fitting(rt, out, function_value,
                                  regular_neighbor_coordinates_2_object< Rt,

Interpolation/test/Interpolation/include/CGAL/_test_regular_neighbors_2.cpp

@@ -28,6 +28,8 @@
 #include <CGAL/Origin.h>
 
 #include <CGAL/regular_neighbor_coordinates_2.h>
+#include <CGAL/Interpolation_gradient_fitting_traits_2.h>
+#include <CGAL/sibson_gradient_fitting.h>
 
 template < class ForwardIterator >
 bool test_norm(ForwardIterator first, ForwardIterator beyond,
@@ -55,6 +57,34 @@ bool test_barycenter(ForwardIterator first, ForwardIterator beyond,
   return p==b;
 }
 
+template <class Point, class FT>
+struct Functor
+{
+  typedef std::pair<FT, bool>                  result_type;
+
+  result_type operator()(const Point&) { return std::make_pair(0., true); }
+};
+
+template <class Triangul>
+void test_gradient_fitting(const Triangul& rt)
+{
+  typedef typename Triangul::Geom_traits::Kernel    K;
+  typedef typename Triangul::Geom_traits::FT        FT;
+  typedef typename Triangul::Geom_traits::Vector_2  Vector_2;
+  typedef typename Triangul::Geom_traits::Point_2   Point;
+
+  typedef typename std::back_insert_iterator<
+                     std::vector<
+                       std::pair< Point, Vector_2 > > >   OutputIterator;
+
+  std::vector<std::pair< Point, Vector_2 > > v;
+  OutputIterator out = std::back_inserter(v);
+  Functor<Point, FT> f;
+
+  CGAL::Interpolation_gradient_fitting_traits_2<K> traits;
+
+  sibson_gradient_fitting_rn_2(rt, out, f, traits);
+}
 
 template <class Triangul>
 void
@@ -113,6 +143,8 @@ _test_regular_neighbors_2( const Triangul & )
     coords.clear();
   }
 
+  test_gradient_fitting(T);
+
   //TESTING a GRID POINT SET
   std::cout << "RN2: Testing grid points." << std::endl;