diff --git a/Interpolation/examples/Interpolation/interpolation_2.cpp b/Interpolation/examples/Interpolation/interpolation_2.cpp
index 96497c203f5..a6f480313ef 100644
--- a/Interpolation/examples/Interpolation/interpolation_2.cpp
+++ b/Interpolation/examples/Interpolation/interpolation_2.cpp
@@ -49,7 +49,7 @@ int main()
   // Create n+m-4 points within a disc of radius 2
   double r_d = 3;
   CGAL::Random rng(1513114263);
-  
+
   CGAL::Random_points_in_disc_2<Point> g(r_d,rng );
   CGAL::cpp11::copy_n( g, n+m, std::back_inserter(points));
 
@@ -114,7 +114,7 @@ int main()
     l_value = CGAL::linear_interpolation(coords.begin(), coords.end(),
                                          norm,
                                          CGAL::Data_access<Point_value_map>(values));
-    
+
     error = CGAL_NTS abs(l_value - exact_value);
     l_total += error;
     if (error > l_max) l_max = error;
@@ -208,6 +208,5 @@ int main()
             << CGAL::to_double(ssquare_total)/n << "  max "
             << CGAL::to_double(ssquare_max)  << std::endl;
 
-  std::cout << "done" << std::endl;
-  return 0;
+  return EXIT_SUCCESS;
 }
diff --git a/Interpolation/examples/Interpolation/interpolation_vertex_with_info_2.cpp b/Interpolation/examples/Interpolation/interpolation_vertex_with_info_2.cpp
index 9e29f538dfe..387fa8f135f 100644
--- a/Interpolation/examples/Interpolation/interpolation_vertex_with_info_2.cpp
+++ b/Interpolation/examples/Interpolation/interpolation_vertex_with_info_2.cpp
@@ -16,72 +16,58 @@
 #include <cassert>
 
 typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
-typedef K::FT                                         Coord_type;
-typedef K::Vector_2                                   Vector;
-typedef K::Point_2                                    Point;
-
+typedef K::FT                                               Coord_type;
+typedef K::Vector_2                                         Vector;
+typedef K::Point_2                                          Point;
 
 template <typename V, typename G>
-struct Value_and_gradient {
-  Value_and_gradient()
-    : value(), gradient(CGAL::NULL_VECTOR)
-  {}
-  
+struct Value_and_gradient
+{
+  Value_and_gradient() : value(), gradient(CGAL::NULL_VECTOR) {}
+
   V value;
   G gradient;
 };
 
-typedef CGAL::Triangulation_vertex_base_with_info_2<Value_and_gradient<Coord_type,Vector>, K>    Vb;
-typedef CGAL::Triangulation_data_structure_2<Vb>      Tds;
-typedef CGAL::Delaunay_triangulation_2<K,Tds>         Delaunay_triangulation;
-typedef Delaunay_triangulation::Vertex_handle         Vertex_handle;
-typedef CGAL::Interpolation_traits_2<K>               Traits;
+typedef CGAL::Triangulation_vertex_base_with_info_2<
+                Value_and_gradient<Coord_type, Vector>, K>  Vb;
+typedef CGAL::Triangulation_data_structure_2<Vb>            Tds;
+typedef CGAL::Delaunay_triangulation_2<K, Tds>              Delaunay_triangulation;
+typedef Delaunay_triangulation::Vertex_handle               Vertex_handle;
+typedef CGAL::Interpolation_traits_2<K>                     Traits;
 
-typedef std::vector< std::pair<Point, Coord_type> >   Coordinate_vector;
+typedef std::vector< std::pair<Point, Coord_type> >         Coordinate_vector;
 
 template <typename V, typename T>
-struct Function_value {
-  typedef V argument_type;
-  typedef std::pair<T, bool> result_type;
-  
-  result_type operator()(const argument_type& a)const
-  {
+struct Function_value
+{
+  typedef V                                                 argument_type;
+  typedef std::pair<T, bool>                                result_type;
+
+  result_type operator()(const argument_type& a) const {
     return result_type(a->info().value, true);
   }
-  
 };
 
-
 template <typename V, typename G>
 struct Function_gradient
-  : public std::iterator<std::output_iterator_tag,void,void,void,void> {
-  
-  typedef V argument_type;
-  typedef std::pair<G,bool> result_type;
-  
-  
-  result_type
-  operator()(const argument_type& a)const
-  {
-    return std::make_pair(a->info().gradient,a->info().gradient != CGAL::NULL_VECTOR) ;
+  : public std::iterator<std::output_iterator_tag, void, void, void, void>
+{
+
+  typedef V                                                 argument_type;
+  typedef std::pair<G, bool>                                result_type;
+
+  result_type operator()(const argument_type& a) const {
+    return std::make_pair(a->info().gradient, a->info().gradient != CGAL::NULL_VECTOR);
   }
 
-
-  const Function_gradient& operator=(const std::pair<V, G>& p) const
-    {
-      p.first->info().gradient = p.second;
-      return *this;
-    }
-
-  const Function_gradient& operator++(int) const
-  {
-    return *this;
-  }
-  
-  const Function_gradient& operator*() const
-  {
+  const Function_gradient& operator=(const std::pair<V, G>& p) const {
+    p.first->info().gradient = p.second;
     return *this;
   }
+
+  const Function_gradient& operator++(int) const { return *this; }
+  const Function_gradient& operator*() const { return *this; }
 };
 
 int main()
@@ -103,24 +89,25 @@ int main()
   // Create n+m-4 points within a disc of radius 2
   double r_d = 3;
   CGAL::Random rng(1513114263);
-  CGAL::Random_points_in_disc_2<Point> g(r_d,rng );
+  CGAL::Random_points_in_disc_2<Point> g(r_d, rng);
   CGAL::cpp11::copy_n( g, n+m, std::back_inserter(points));
 
   Delaunay_triangulation T;
 
-  Function_value<Vertex_handle,Coord_type> function_value;
-  Function_gradient<Vertex_handle,Vector> function_gradient;
-  
+  Function_value<Vertex_handle, Coord_type> function_value;
+  Function_gradient<Vertex_handle, Vector> function_gradient;
+
   //parameters for quadratic function:
   Coord_type alpha = Coord_type(1.0),
-      beta1 = Coord_type(2.0),
-      beta2 = Coord_type(1.0),
-      gamma1 = Coord_type(0.3),
-      gamma2 = Coord_type(0.0),
-      gamma3 = Coord_type(0.0),
-      gamma4 = Coord_type(0.3);
+             beta1 = Coord_type(2.0),
+             beta2 = Coord_type(1.0),
+             gamma1 = Coord_type(0.3),
+             gamma2 = Coord_type(0.0),
+             gamma3 = Coord_type(0.0),
+             gamma4 = Coord_type(0.3);
 
-  for(int j=0; j<n ; j++){
+  for(int j=0; j<n ; j++)
+  {
     Vertex_handle vh = T.insert(points[j]);
 
     //determine function value/gradient:
@@ -129,8 +116,8 @@ int main()
 
     Coord_type value = alpha + beta1*x + beta2*y + gamma1*(x*x) +
                        gamma4*(y*y) + (gamma2+ gamma3) *(x*y);
-    Vector gradient(beta1+ (gamma2+ gamma3)*y + Coord_type(2)*(gamma1*x),
-                    beta2+ (gamma2+ gamma3)*x + Coord_type(2)*(gamma4*y));
+    Vector gradient(beta1+ (gamma2 + gamma3)*y + Coord_type(2)*(gamma1*x),
+                    beta2+ (gamma2 + gamma3)*x + Coord_type(2)*(gamma4*y));
     vh->info().value = value;
     vh->info().gradient = gradient;
   }
@@ -138,15 +125,16 @@ int main()
   //variables for statistics:
   std::pair<Coord_type, bool> res;
   Coord_type error, l_total = Coord_type(0),
-      q_total(l_total), f_total(l_total), s_total(l_total),
-      ssquare_total(l_total), l_max(l_total),
-      q_max(l_total), f_max(l_total), s_max(l_total),
-      ssquare_max(l_total),
-      total_value(l_total), l_value(l_total);
+                    q_total(l_total), f_total(l_total), s_total(l_total),
+                    ssquare_total(l_total), l_max(l_total),
+                    q_max(l_total), f_max(l_total), s_max(l_total),
+                    ssquare_max(l_total),
+                    total_value(l_total), l_value(l_total);
   int failure(0);
 
   //interpolation + error statistics
-  for(int i=n;i<n+m;i++){
+  for(int i=n; i<n+m; i++)
+  {
     Coord_type x(points[i].x());
     Coord_type y(points[i].y());
 
@@ -155,32 +143,23 @@ int main()
 
     total_value += exact_value;
 
-
-    
     //Coordinate_vector:
-    std::vector< std::pair< Vertex_handle, Coord_type > > coords;
-    typedef CGAL::Identity<std::pair< Vertex_handle, Coord_type > > Identity;
-    Coord_type norm =
-        CGAL::natural_neighbor_coordinates_2(T,
-                                             points[i],
-                                             std::back_inserter(coords),
-                                             Identity()).second;
-
-    assert(norm>0);
+    std::vector< std::pair<Vertex_handle, Coord_type> > coords;
+    typedef CGAL::Identity<std::pair< Vertex_handle, Coord_type> > Identity;
+    Coord_type norm = CGAL::natural_neighbor_coordinates_2(T,
+                                                           points[i],
+                                                           std::back_inserter(coords),
+                                                           Identity()).second;
+    assert(norm > 0);
 
     //linear interpolant:
-    l_value =
-      CGAL::linear_interpolation(coords.begin(), coords.end(),
-                                 norm,
-                                 function_value);
-
-
+    l_value = CGAL::linear_interpolation(coords.begin(), coords.end(),
+                                         norm, function_value);
 
     error = CGAL_NTS abs(l_value - exact_value);
     l_total += error;
     if (error > l_max) l_max = error;
 
-    
     //Farin interpolant:
     res = CGAL::farin_c1_interpolation(coords.begin(),
                                        coords.end(), norm,points[i],
@@ -189,59 +168,79 @@ int main()
                                        Traits());
 
 
-    if(res.second){
+    if(res.second)
+    {
       error = CGAL_NTS abs(res.first - exact_value);
       f_total += error;
-      if (error > f_max) f_max = error;
-    } else ++failure;
+      if (error > f_max)
+        f_max = error;
+    }
+    else
+    {
+      ++failure;
+    }
 
     //quadratic interpolant:
-    res =  CGAL::quadratic_interpolation(coords.begin(), coords.end(),
-                                         norm,points[i],
-                                         function_value,
-                                         function_gradient,
-                                         Traits());
-    if(res.second){
+    res = CGAL::quadratic_interpolation(coords.begin(), coords.end(),
+                                        norm, points[i],
+                                        function_value,
+                                        function_gradient,
+                                        Traits());
+
+    if(res.second)
+    {
       error = CGAL_NTS abs(res.first - exact_value);
       q_total += error;
-      if (error > q_max) q_max = error;
-    } else ++failure;
-
-
+      if (error > q_max)
+        q_max = error;
+    }
+    else
+    {
+      ++failure;
+    }
 
     //Sibson interpolant: version without sqrt:
-    res =  CGAL::sibson_c1_interpolation_square(coords.begin(),
-                                                coords.end(), norm,
-                                                points[i],
-                                                function_value,
-                                                function_gradient,
-                                                Traits());
+    res = CGAL::sibson_c1_interpolation_square(coords.begin(),
+                                               coords.end(), norm,
+                                               points[i],
+                                               function_value,
+                                               function_gradient,
+                                               Traits());
     //error statistics
-    if(res.second){
+    if(res.second)
+    {
       error = CGAL_NTS abs(res.first - exact_value);
       ssquare_total += error;
-      if (error > ssquare_max) ssquare_max = error;
-    } else ++failure;
+      if (error > ssquare_max)
+        ssquare_max = error;
+    }
+    else
+    {
+      ++failure;
+    }
 
     //with sqrt(the traditional):
-    res =  CGAL::sibson_c1_interpolation(coords.begin(),
-                                         coords.end(), norm,
-                                         points[i],
-                                         function_value,
-                                         function_gradient,
-                                         Traits());
+    res = CGAL::sibson_c1_interpolation(coords.begin(),
+                                        coords.end(), norm,
+                                        points[i],
+                                        function_value,
+                                        function_gradient,
+                                        Traits());
 
     //error statistics
-    if(res.second){
+    if(res.second)
+    {
       error = CGAL_NTS abs(res.first - exact_value);
       s_total += error;
-      if (error > s_max) s_max = error;
-    } else ++failure;
-
+      if (error > s_max)
+        s_max = error;
+    }
+    else
+    {
+      ++failure;
+    }
   }
 
-
-
   /************** end of Interpolation: dump statistics **************/
   std::cout << "Result: -----------------------------------" << std::endl;
   std::cout <<  "Interpolation of '" << alpha <<" + "
@@ -270,7 +269,5 @@ int main()
             << CGAL::to_double(ssquare_total)/n << "  max "
             << CGAL::to_double(ssquare_max)  << std::endl;
 
-
-  std::cout << "done" << std::endl;
-  return 0;
+  return EXIT_SUCCESS;
 }
diff --git a/Interpolation/examples/Interpolation/linear_interpolation_2.cpp b/Interpolation/examples/Interpolation/linear_interpolation_2.cpp
index 701c59c5ec6..d399e2b6c11 100644
--- a/Interpolation/examples/Interpolation/linear_interpolation_2.cpp
+++ b/Interpolation/examples/Interpolation/linear_interpolation_2.cpp
@@ -6,42 +6,39 @@
 #include <CGAL/interpolation_functions.h>
 
 typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
-typedef CGAL::Delaunay_triangulation_2<K>             Delaunay_triangulation;
-typedef CGAL::Interpolation_traits_2<K>               Traits;
-typedef K::FT                                         Coord_type;
-typedef K::Point_2                                    Point;
+typedef CGAL::Delaunay_triangulation_2<K>                   Delaunay_triangulation;
+typedef CGAL::Interpolation_traits_2<K>                     Traits;
+typedef K::FT                                               Coord_type;
+typedef K::Point_2                                          Point;
 
 int main()
 {
   Delaunay_triangulation T;
-  std::map<Point, Coord_type, K::Less_xy_2> function_values;
-  typedef CGAL::Data_access< std::map<Point, Coord_type, K::Less_xy_2 > >
-      Value_access;
 
+  typedef std::map<Point, Coord_type, K::Less_xy_2>         Coord_map;
+  typedef CGAL::Data_access<Coord_map>                      Value_access;
+
+  Coord_map function_values;
   Coord_type a(0.25), bx(1.3), by(-0.7);
 
   for (int y=0 ; y<3 ; y++){
     for (int x=0 ; x<3 ; x++){
       K::Point_2 p(x,y);
       T.insert(p);
-      function_values.insert(std::make_pair(p,a + bx* x+ by*y));
+      function_values.insert(std::make_pair(p, a + bx*x + by*y));
     }
   }
 
   //coordinate computation
-  K::Point_2 p(1.3,0.34);
-  std::vector< std::pair< Point, Coord_type > > coords;
-  Coord_type norm =
-      CGAL::natural_neighbor_coordinates_2
-      (T, p,std::back_inserter(coords)).second;
+  K::Point_2 p(1.3, 0.34);
+  std::vector<std::pair<Point, Coord_type> > coords;
 
-  Coord_type res =  CGAL::linear_interpolation(coords.begin(), coords.end(),
-                                               norm,
+  Coord_type norm = CGAL::natural_neighbor_coordinates_2(T, p, std::back_inserter(coords)).second;
+  Coord_type res =  CGAL::linear_interpolation(coords.begin(), coords.end(), norm,
                                                Value_access(function_values));
 
   std::cout << "Tested interpolation on " << p << " interpolation: "
-            << res << " exact: " << a + bx* p.x()+ by* p.y()<< std::endl;
+            << res << " exact: " << a + bx*p.x() + by*p.y() << std::endl;
 
-  std::cout << "done" << std::endl;
-  return 0;
+  return EXIT_SUCCESS;
 }
diff --git a/Interpolation/examples/Interpolation/nn_coordinates_2.cpp b/Interpolation/examples/Interpolation/nn_coordinates_2.cpp
index 96e3d6567be..897d815f268 100644
--- a/Interpolation/examples/Interpolation/nn_coordinates_2.cpp
+++ b/Interpolation/examples/Interpolation/nn_coordinates_2.cpp
@@ -1,36 +1,40 @@
 #include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
+
 #include <CGAL/Delaunay_triangulation_2.h>
 #include <CGAL/natural_neighbor_coordinates_2.h>
 
+#include <iostream>
+#include <iterator>
+#include <utility>
+#include <vector>
+
 typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
-typedef CGAL::Delaunay_triangulation_2<K>              Delaunay_triangulation;
-typedef std::vector< std::pair< K::Point_2, K::FT  > > Point_coordinate_vector;
+typedef CGAL::Delaunay_triangulation_2<K>                   Delaunay_triangulation;
+typedef std::vector< std::pair<K::Point_2, K::FT> >         Point_coordinate_vector;
 
 int main()
 {
   Delaunay_triangulation dt;
 
-  for (int y=0 ; y<3 ; y++)
-    for (int x=0 ; x<3 ; x++)
-      dt.insert(K::Point_2(x,y));
+  for (int y=0; y<3; y++)
+    for (int x=0; x<3; x++)
+      dt.insert(K::Point_2(x, y));
 
   //coordinate computation
   K::Point_2 p(1.2, 0.7);
   Point_coordinate_vector coords;
-  CGAL::Triple<std::back_insert_iterator<Point_coordinate_vector>,
-               K::FT, bool> result =
+  CGAL::Triple<std::back_insert_iterator<Point_coordinate_vector>, K::FT, bool> result =
       CGAL::natural_neighbor_coordinates_2(dt, p, std::back_inserter(coords));
 
-  if(!result.third){
-    std::cout << "The coordinate computation was not successful."
-              << std::endl;
-    std::cout << "The point (" <<p << ") lies outside the convex hull."
-              << std::endl;
+  if(!result.third)
+  {
+    std::cout << "The coordinate computation was not successful." << std::endl;
+    std::cout << "The point (" << p << ") lies outside the convex hull." << std::endl;
   }
 
-  K::FT  norm = result.second;
+  K::FT norm = result.second;
   std::cout << "Coordinate computation successful." << std::endl;
-  std::cout << "Normalization factor: " <<norm << std::endl;
-  std::cout << "done" << std::endl;
-  return 0;
+  std::cout << "Normalization factor: " << norm << std::endl;
+
+  return EXIT_SUCCESS;
 }
diff --git a/Interpolation/examples/Interpolation/nn_coordinates_3.cpp b/Interpolation/examples/Interpolation/nn_coordinates_3.cpp
index 1ff94563f89..8fc91c7937d 100644
--- a/Interpolation/examples/Interpolation/nn_coordinates_3.cpp
+++ b/Interpolation/examples/Interpolation/nn_coordinates_3.cpp
@@ -1,27 +1,33 @@
 #include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
+
 #include <CGAL/Delaunay_triangulation_3.h>
 #include <CGAL/natural_neighbor_coordinates_3.h>
+
 #include <fstream>
+#include <iostream>
+#include <iterator>
+#include <utility>
+#include <vector>
 
 typedef double NT; //Number Type
 
-typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
+typedef CGAL::Exact_predicates_inexact_constructions_kernel    K;
 
-typedef K::Point_3 Point3;
-typedef K::Vector_3 Vector3;
-typedef K::Sphere_3 Sphere_3;
+typedef K::Point_3                                             Point3;
+typedef K::Vector_3                                            Vector3;
+typedef K::Sphere_3                                            Sphere_3;
 
 typedef CGAL::Delaunay_triangulation_3<K, CGAL::Fast_location> Dh;
 
-typedef Dh::Facet Facet;
-typedef Dh::Vertex_handle Vertex_handle;
-typedef Dh::Cell_handle Cell_handle;
-typedef Dh::Finite_vertices_iterator Finite_vertices_iterator;
-typedef Dh::Vertex_iterator Vertex_iterator;
-typedef Dh::Facet_circulator Facet_circulator;
-typedef Dh::Cell_iterator Cell_iterator;
+typedef Dh::Facet                                              Facet;
+typedef Dh::Vertex_handle                                      Vertex_handle;
+typedef Dh::Cell_handle                                        Cell_handle;
+typedef Dh::Finite_vertices_iterator                           Finite_vertices_iterator;
+typedef Dh::Vertex_iterator                                    Vertex_iterator;
+typedef Dh::Facet_circulator                                   Facet_circulator;
+typedef Dh::Cell_iterator                                      Cell_iterator;
 
-typedef K::Construct_circumcenter_3 Construct_circumcenter_3;
+typedef K::Construct_circumcenter_3                            Construct_circumcenter_3;
 
 int main()
 {
@@ -35,11 +41,11 @@ int main()
 
   Point3 pp[3];
   std::cout << "Consider the natural coordinates of P1, P2 and P3 with regard to the triangulation of data/points3 " << std::endl;
-  pp[0]=Point3(106.55,112.57,110.4); //inside data/points3 convex hull
+  pp[0] = Point3(106.55,112.57,110.4); //inside data/points3 convex hull
   std::cout << "P1 is inside the convex hull" << std::endl;
-  pp[1]=Point3(250,100,140); //on data/points3 convex hull
+  pp[1] = Point3(250,100,140); //on data/points3 convex hull
   std::cout << "P2 is on a vertex of the triangulation" << std::endl;
-  pp[2]=Point3(0,0,0); //outside data/points3 convex hull
+  pp[2] = Point3(0,0,0); //outside data/points3 convex hull
   std::cout << "P2 is outside the convex hull" << std::endl;
 
   for(int ii=0; ii<3; ++ii)
@@ -48,9 +54,9 @@ int main()
     std::vector< std::pair< Vertex_handle,NT> > coor_sibson;
     NT norm_coeff_laplace, norm_coeff_sibson;
 
-    std::cout << "Point P"<< ii+1 << " : "<<pp[ii].x() << " "
-              << pp[ii].y() << " "
-              << pp[ii].z() << std::endl;
+    std::cout << "Point P" << ii+1 << " : " << pp[ii].x() << " "
+                                            << pp[ii].y() << " "
+                                            << pp[ii].z() << std::endl;
 
     laplace_natural_neighbor_coordinates_3(triangulation,pp[ii],
                                            std::back_inserter(coor_laplace),
@@ -67,6 +73,7 @@ int main()
     sibson_natural_neighbor_coordinates_3(triangulation,pp[ii],
                                           std::back_inserter(coor_sibson),
                                           norm_coeff_sibson);
+
     std::cout << "Linear combination of natural neighbors with Sibson natural coordinates" << std::endl;
     std::cout << " + correctness (ensured only with an exact number type)" << std::endl;
     std::cout << is_correct_natural_neighborhood(triangulation,pp[ii],
@@ -76,6 +83,5 @@ int main()
               << std::endl;
   }
 
-  std::cout << "done" << std::endl;
-  return 0;
+  return EXIT_SUCCESS;
 }
diff --git a/Interpolation/examples/Interpolation/rn_coordinates_2.cpp b/Interpolation/examples/Interpolation/rn_coordinates_2.cpp
index 7402eae3be9..dded72129fe 100644
--- a/Interpolation/examples/Interpolation/rn_coordinates_2.cpp
+++ b/Interpolation/examples/Interpolation/rn_coordinates_2.cpp
@@ -1,40 +1,42 @@
 #include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
+
 #include <CGAL/Regular_triangulation_2.h>
 #include <CGAL/regular_neighbor_coordinates_2.h>
 
+#include <iostream>
+#include <iterator>
+#include <vector>
+#include <utility>
+
 typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
 
-typedef CGAL::Regular_triangulation_2<K>               Regular_triangulation;
-typedef Regular_triangulation::Bare_point              Bare_point;
-typedef Regular_triangulation::Weighted_point          Weighted_point;
-typedef std::vector< std::pair< Weighted_point, K::FT  > > Point_coordinate_vector;
+typedef CGAL::Regular_triangulation_2<K>                    Regular_triangulation;
+typedef Regular_triangulation::Bare_point                   Bare_point;
+typedef Regular_triangulation::Weighted_point               Weighted_point;
+typedef std::vector<std::pair<Weighted_point, K::FT> >      Point_coordinate_vector;
 
 int main()
 {
   Regular_triangulation rt;
 
-  for (int y=0 ; y<3 ; y++)
-    for (int x=0 ; x<3 ; x++)
+  for (int y=0; y<3; y++)
+    for (int x=0; x<3; x++)
       rt.insert(Weighted_point(Bare_point(x,y), 0));
 
   //coordinate computation
-  Weighted_point wp(Bare_point(1.2, 0.7),2);
-  Point_coordinate_vector  coords;
-  CGAL::Triple<std::back_insert_iterator<Point_coordinate_vector>,
-               K::FT, bool> result =
+  Weighted_point wp(Bare_point(1.2, 0.7), 2);
+  Point_coordinate_vector coords;
+  CGAL::Triple<std::back_insert_iterator<Point_coordinate_vector>, K::FT, bool> result =
       CGAL::regular_neighbor_coordinates_2(rt, wp, std::back_inserter(coords));
 
   if(!result.third){
-    std::cout << "The coordinate computation was not successful."
-              << std::endl;
-    std::cout << "The point (" <<wp.point() << ") lies outside the convex hull."
-              << std::endl;
+    std::cout << "The coordinate computation was not successful." << std::endl;
+    std::cout << "The point (" <<wp.point() << ") lies outside the convex hull." << std::endl;
   }
 
   K::FT  norm = result.second;
   std::cout << "Coordinate computation successful." << std::endl;
-  std::cout << "Normalization factor: " <<norm << std::endl;
+  std::cout << "Normalization factor: " << norm << std::endl;
 
-  std::cout << "done" << std::endl;
-  return 0;
+  return EXIT_SUCCESS;
 }
diff --git a/Interpolation/examples/Interpolation/sibson_interpolation_2.cpp b/Interpolation/examples/Interpolation/sibson_interpolation_2.cpp
index 449e65e0f83..1499e977df4 100644
--- a/Interpolation/examples/Interpolation/sibson_interpolation_2.cpp
+++ b/Interpolation/examples/Interpolation/sibson_interpolation_2.cpp
@@ -6,14 +6,20 @@
 #include <CGAL/sibson_gradient_fitting.h>
 #include <CGAL/interpolation_functions.h>
 
-typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
-typedef CGAL::Delaunay_triangulation_2<K>               Delaunay_triangulation;
-typedef CGAL::Interpolation_gradient_fitting_traits_2<K> Traits;
+#include <iostream>
+#include <iterator>
+#include <map>
+#include <utility>
+#include <vector>
 
-typedef K::FT                                            Coord_type;
-typedef K::Point_2                                       Point;
-typedef std::map<Point, Coord_type, K::Less_xy_2>        Point_value_map ;
-typedef std::map<Point, K::Vector_2 , K::Less_xy_2 >     Point_vector_map;
+typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
+typedef CGAL::Delaunay_triangulation_2<K>                   Delaunay_triangulation;
+typedef CGAL::Interpolation_gradient_fitting_traits_2<K>    Traits;
+
+typedef K::FT                                               Coord_type;
+typedef K::Point_2                                          Point;
+typedef std::map<Point, Coord_type, K::Less_xy_2>           Point_value_map ;
+typedef std::map<Point, K::Vector_2 , K::Less_xy_2>         Point_vector_map;
 
 int main()
 {
@@ -24,16 +30,16 @@ int main()
 
   //parameters for spherical function:
   Coord_type a(0.25), bx(1.3), by(-0.7), c(0.2);
-  for (int y=0 ; y<4 ; y++){
-    for (int x=0 ; x<4 ; x++){
+  for (int y=0; y<4; y++) {
+    for (int x=0; x<4; x++) {
       K::Point_2 p(x,y);
       T.insert(p);
       function_values.insert(std::make_pair(p,a + bx* x+ by*y + c*(x*x+y*y)));
     }
   }
 
-  sibson_gradient_fitting_nn_2(T,std::inserter(function_gradients,
-                                               function_gradients.begin()),
+  sibson_gradient_fitting_nn_2(T, std::inserter(function_gradients,
+                                                function_gradients.begin()),
                                CGAL::Data_access<Point_value_map>(function_values),
                                Traits());
 
@@ -41,8 +47,8 @@ int main()
   {
     std::cout << it->first << "  "  << it->second << std::endl;
   }
-  //coordinate computation
-  K::Point_2 p(1.6,1.4);
+  // coordinate computation
+  K::Point_2 p(1.6, 1.4);
   std::vector< std::pair< Point, Coord_type > > coords;
   Coord_type norm = CGAL::natural_neighbor_coordinates_2(T, p, std::back_inserter
                                                          (coords)).second;
@@ -50,23 +56,21 @@ int main()
 
   //Sibson interpolant: version without sqrt:
   std::pair<Coord_type, bool> res =
-      CGAL::sibson_c1_interpolation_square(
-        coords.begin(),
-        coords.end(),norm,p,
-        CGAL::Data_access<Point_value_map>(function_values),
-        CGAL::Data_access<Point_vector_map>(function_gradients),
-        Traits());
+    CGAL::sibson_c1_interpolation_square(coords.begin(),
+                                         coords.end(),norm,p,
+                                         CGAL::Data_access<Point_value_map>(function_values),
+                                         CGAL::Data_access<Point_vector_map>(function_gradients),
+                                         Traits());
 
   if(res.second)
     std::cout << "Tested interpolation on " << p
               << " interpolation: " << res.first << " exact: "
-              << a + bx * p.x()+ by * p.y()+ c*(p.x()*p.x()+p.y()*p.y())
+              << a + bx*p.x() + by*p.y() + c*(p.x()*p.x()+p.y()*p.y())
               << std::endl;
   else
     std::cout << "C^1 Interpolation not successful." << std::endl
-              << " not all function_gradients are provided."  << std::endl
+              << " not all function_gradients are provided." << std::endl
               << " You may resort to linear interpolation." << std::endl;
 
-  std::cout << "done" << std::endl;
-  return 0;
+  return EXIT_SUCCESS;
 }
diff --git a/Interpolation/examples/Interpolation/sibson_interpolation_rn_2.cpp b/Interpolation/examples/Interpolation/sibson_interpolation_rn_2.cpp
index 171964fdffb..b72cd85996a 100644
--- a/Interpolation/examples/Interpolation/sibson_interpolation_rn_2.cpp
+++ b/Interpolation/examples/Interpolation/sibson_interpolation_rn_2.cpp
@@ -1,28 +1,34 @@
 #include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
-#include <CGAL/Regular_triangulation_2.h>
 
 #include <CGAL/natural_neighbor_coordinates_2.h>
 #include <CGAL/Interpolation_gradient_fitting_traits_2.h>
 #include <CGAL/sibson_gradient_fitting.h>
 #include <CGAL/interpolation_functions.h>
 
+#include <CGAL/Regular_triangulation_2.h>
+
+#include <iostream>
+#include <iterator>
+#include <map>
+#include <utility>
+#include <vector>
+
 typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
-typedef CGAL::Regular_triangulation_2<K>               Regular_triangulation;
-typedef CGAL::Interpolation_gradient_fitting_traits_2<K> Traits;
+typedef CGAL::Regular_triangulation_2<K>                    Regular_triangulation;
+typedef CGAL::Interpolation_gradient_fitting_traits_2<K>    Traits;
 
-typedef K::FT                                            Coord_type;
-typedef K::Weighted_point_2                              Point;
+typedef K::FT                                               Coord_type;
+typedef K::Weighted_point_2                                 Point;
 
-struct Less {
-  bool operator()(const Point& p, const Point& q) const
-  {
+struct Less
+{
+  bool operator()(const Point& p, const Point& q) const {
     return K::Less_xy_2()(p.point(), q.point());
   }
 };
 
-
-typedef std::map<Point, Coord_type, Less>        Point_value_map ;
-typedef std::map<Point, K::Vector_2 , Less >     Point_vector_map;
+typedef std::map<Point, Coord_type, Less>                   Point_value_map ;
+typedef std::map<Point, K::Vector_2 , Less>                 Point_vector_map;
 
 int main()
 {
@@ -33,11 +39,11 @@ int main()
 
   //parameters for spherical function:
   Coord_type a(0.25), bx(1.3), by(-0.7), c(0.2);
-  for (int y=0 ; y<4 ; y++){
-    for (int x=0 ; x<4 ; x++){
+  for (int y=0; y<4; y++) {
+    for (int x=0; x<4; x++) {
       Point p(x,y);
       T.insert(p);
-      function_values.insert(std::make_pair(p,a + bx* x+ by*y + c*(x*x+y*y)));
+      function_values.insert(std::make_pair(p, a + bx*x + by*y + c*(x*x+y*y)));
     }
   }
 
@@ -46,36 +52,34 @@ int main()
                                CGAL::Data_access<Point_value_map>(function_values),
                                Traits());
 
-  for(Point_vector_map::iterator it = function_gradients.begin(); it != function_gradients.end(); ++it)
-  {
+  for(Point_vector_map::iterator it = function_gradients.begin();
+                                 it != function_gradients.end(); ++it) {
     std::cout << it->first << "  "  << it->second << std::endl;
   }
-  //coordinate computation
-  Point p(1.6,1.4);
-  std::vector< std::pair< Point, Coord_type > > coords;
-  Coord_type norm = CGAL::regular_neighbor_coordinates_2(T, p, std::back_inserter
-                                                         (coords)).second;
 
+  //coordinate computation
+  Point p(1.6, 1.4);
+  std::vector<std::pair<Point, Coord_type> > coords;
+  Coord_type norm = CGAL::regular_neighbor_coordinates_2(T, p, std::back_inserter(coords)).second;
 
   //Sibson interpolant: version without sqrt:
-  std::pair<Coord_type, bool> res =
-      CGAL::sibson_c1_interpolation_square(
-        coords.begin(),
-        coords.end(),norm,p,
-        CGAL::Data_access<Point_value_map>(function_values),
-        CGAL::Data_access<Point_vector_map>(function_gradients),
-        Traits());
+  std::pair<Coord_type, bool> res = CGAL::sibson_c1_interpolation_square(coords.begin(),
+                                                                         coords.end(),
+                                                                         norm,
+                                                                         p,
+                                                                         CGAL::Data_access<Point_value_map>(function_values),
+                                                                         CGAL::Data_access<Point_vector_map>(function_gradients),
+                                                                         Traits());
 
   if(res.second)
     std::cout << "Tested interpolation on " << p
               << " interpolation: " << res.first << " exact: "
-              << a + bx * p.x()+ by * p.y()+ c*(p.x()*p.x()+p.y()*p.y())
+              << a + bx*p.x() + by*p.y()+ c*(p.x()*p.x()+p.y()*p.y())
               << std::endl;
   else
     std::cout << "C^1 Interpolation not successful." << std::endl
               << " not all function_gradients are provided."  << std::endl
               << " You may resort to linear interpolation." << std::endl;
 
-  std::cout << "done" << std::endl;
   return 0;
 }
diff --git a/Interpolation/examples/Interpolation/sibson_interpolation_rn_vertex_with_info_2.cpp b/Interpolation/examples/Interpolation/sibson_interpolation_rn_vertex_with_info_2.cpp
index b50a1f06bb7..176c35416a7 100644
--- a/Interpolation/examples/Interpolation/sibson_interpolation_rn_vertex_with_info_2.cpp
+++ b/Interpolation/examples/Interpolation/sibson_interpolation_rn_vertex_with_info_2.cpp
@@ -1,95 +1,87 @@
 #include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
-#include <CGAL/Triangulation_vertex_base_with_info_2.h>
-#include <CGAL/Regular_triangulation_2.h>
 
 #include <CGAL/natural_neighbor_coordinates_2.h>
 #include <CGAL/Interpolation_gradient_fitting_traits_2.h>
 #include <CGAL/sibson_gradient_fitting.h>
 #include <CGAL/interpolation_functions.h>
 
-typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
-typedef CGAL::Interpolation_gradient_fitting_traits_2<K> Traits;
+#include <CGAL/Triangulation_vertex_base_with_info_2.h>
+#include <CGAL/Regular_triangulation_2.h>
 
-typedef K::FT                                            Coord_type;
-typedef K::Weighted_point_2                              Point;
-typedef K::Vector_2                                      Vector;
+#include <iostream>
+#include <iterator>
+#include <utility>
+#include <vector>
+
+typedef CGAL::Exact_predicates_inexact_constructions_kernel    K;
+typedef CGAL::Interpolation_gradient_fitting_traits_2<K>       Traits;
+
+typedef K::FT                                                  Coord_type;
+typedef K::Weighted_point_2                                    Point;
+typedef K::Vector_2                                            Vector;
 
 template <typename V, typename G>
-struct Value_and_gradient {
-  Value_and_gradient()
-    : value(), gradient(CGAL::NULL_VECTOR)
-  {}
-  
+struct Value_and_gradient
+{
+  Value_and_gradient() : value(), gradient(CGAL::NULL_VECTOR) {}
+
   V value;
   G gradient;
 };
 
-typedef CGAL::Triangulation_vertex_base_with_info_2<Value_and_gradient<Coord_type,Vector>, K, CGAL::Regular_triangulation_vertex_base_2<K> >    Vb;
+typedef CGAL::Triangulation_vertex_base_with_info_2<
+                Value_and_gradient<Coord_type, Vector>, K,
+                CGAL::Regular_triangulation_vertex_base_2<K> > Vb;
 typedef CGAL::Regular_triangulation_face_base_2<K> Fb;
-typedef CGAL::Triangulation_data_structure_2<Vb,Fb>     Tds;
-typedef CGAL::Regular_triangulation_2<K,Tds>            Regular_triangulation;
-typedef Regular_triangulation::Vertex_handle            Vertex_handle;
-
+typedef CGAL::Triangulation_data_structure_2<Vb, Fb>           Tds;
+typedef CGAL::Regular_triangulation_2<K, Tds>                  Regular_triangulation;
+typedef Regular_triangulation::Vertex_handle                   Vertex_handle;
 
 template <typename V, typename T>
-struct Function_value {
-  typedef V argument_type;
-  typedef std::pair<T, bool> result_type;
-  
-  result_type operator()(const argument_type& a)const
-  {
+struct Function_value
+{
+  typedef V                                                    argument_type;
+  typedef std::pair<T, bool>                                   result_type;
+
+  result_type operator()(const argument_type& a) const {
     return result_type(a->info().value, true);
   }
-  
 };
 
-
 template <typename V, typename G>
 struct Function_gradient
-  : public std::iterator<std::output_iterator_tag,void,void,void,void> {
-  
-  typedef V argument_type;
-  typedef std::pair<G,bool> result_type;
-  
-  
-  result_type
-  operator()(const argument_type& a)const
-  {
-    return std::make_pair(a->info().gradient,a->info().gradient != CGAL::NULL_VECTOR) ;
+  : public std::iterator<std::output_iterator_tag, void, void, void, void>
+{
+  typedef V                                                    argument_type;
+  typedef std::pair<G, bool>                                   result_type;
+
+  result_type operator()(const argument_type& a) const {
+    return std::make_pair(a->info().gradient, a->info().gradient != CGAL::NULL_VECTOR);
   }
 
-
-  const Function_gradient& operator=(const std::pair<V, G>& p) const
-    {
-      p.first->info().gradient = p.second;
-      return *this;
-    }
-
-  const Function_gradient& operator++(int) const
-  {
-    return *this;
-  }
-  
-  const Function_gradient& operator*() const
-  {
+  const Function_gradient& operator=(const std::pair<V, G>& p) const {
+    p.first->info().gradient = p.second;
     return *this;
   }
+
+  const Function_gradient& operator++(int) const { return *this; }
+  const Function_gradient& operator*() const { return *this; }
 };
 
 int main()
 {
   Regular_triangulation rt;
 
-  Function_value<Vertex_handle,Coord_type> function_value;
-  Function_gradient<Vertex_handle,Vector> function_gradient;
-  
+  Function_value<Vertex_handle, Coord_type> function_value;
+  Function_gradient<Vertex_handle, Vector> function_gradient;
+
   //parameters for spherical function:
   Coord_type a(0.25), bx(1.3), by(-0.7), c(0.2);
-  for (int y=0 ; y<4 ; y++){
-    for (int x=0 ; x<4 ; x++){
+  for (int y=0; y<4; y++) {
+    for (int x=0; x<4; x++) {
       Point p(x,y);
       Vertex_handle vh = rt.insert(p);
-      Coord_type value = a + bx* x+ by*y + c*(x*x+y*y);
+      Coord_type value = a + bx*x + by*y + c*(x*x+y*y);
       vh->info().value = value;
     }
   }
@@ -99,16 +91,16 @@ int main()
                                function_value,
                                CGAL::Identity<std::pair<Vertex_handle, Vector> >(),
                                Traits());
+
   //coordinate computation
   Point p(1.6,1.4);
-  std::vector< std::pair< Vertex_handle, Coord_type > > coords;
-  typedef CGAL::Identity<std::pair< Vertex_handle, Coord_type > > Identity;
+  std::vector<std::pair<Vertex_handle, Coord_type> > coords;
+  typedef CGAL::Identity<std::pair<Vertex_handle, Coord_type> > Identity;
   Coord_type norm = CGAL::regular_neighbor_coordinates_2(rt,
                                                          p,
                                                          std::back_inserter(coords),
                                                          Identity()).second;
 
-
   //Sibson interpolant: version without sqrt:
   std::pair<Coord_type, bool> res = CGAL::sibson_c1_interpolation_square(coords.begin(),
                                                                          coords.end(),
@@ -117,7 +109,7 @@ int main()
                                                                          function_value,
                                                                          function_gradient,
                                                                          Traits());
-  
+
   if(res.second)
     std::cout << "Tested interpolation on " << p
               << " interpolation: " << res.first << " exact: "
@@ -128,6 +120,5 @@ int main()
               << " not all function_gradients are provided."  << std::endl
               << " You may resort to linear interpolation." << std::endl;
 
-  std::cout << "done" << std::endl;
-  return 0;
+  return EXIT_SUCCESS;
 }
diff --git a/Interpolation/examples/Interpolation/sibson_interpolation_vertex_with_info_2.cpp b/Interpolation/examples/Interpolation/sibson_interpolation_vertex_with_info_2.cpp
index b18876e719f..4bfde90b3e7 100644
--- a/Interpolation/examples/Interpolation/sibson_interpolation_vertex_with_info_2.cpp
+++ b/Interpolation/examples/Interpolation/sibson_interpolation_vertex_with_info_2.cpp
@@ -7,85 +7,75 @@
 #include <CGAL/sibson_gradient_fitting.h>
 #include <CGAL/interpolation_functions.h>
 
-typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
-typedef CGAL::Interpolation_gradient_fitting_traits_2<K> Traits;
+#include <iostream>
+#include <iterator>
+#include <utility>
+#include <vector>
 
-typedef K::FT                                            Coord_type;
-typedef K::Point_2                                       Point;
-typedef K::Vector_2                                      Vector;
+typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
+typedef CGAL::Interpolation_gradient_fitting_traits_2<K>    Traits;
+
+typedef K::FT                                               Coord_type;
+typedef K::Point_2                                          Point;
+typedef K::Vector_2                                         Vector;
 
 template <typename V, typename G>
-struct Value_and_gradient {
-  Value_and_gradient()
-    : value(), gradient(CGAL::NULL_VECTOR)
-  {}
-  
+struct Value_and_gradient
+{
+  Value_and_gradient() : value(), gradient(CGAL::NULL_VECTOR) {}
+
   V value;
   G gradient;
 };
 
-typedef CGAL::Triangulation_vertex_base_with_info_2<Value_and_gradient<Coord_type,Vector>, K>    Vb;
-typedef CGAL::Triangulation_data_structure_2<Vb>         Tds;
-typedef CGAL::Delaunay_triangulation_2<K,Tds>            Delaunay_triangulation;
-typedef Delaunay_triangulation::Vertex_handle            Vertex_handle;
-
+typedef CGAL::Triangulation_vertex_base_with_info_2<
+                Value_and_gradient<Coord_type, Vector>, K>  Vb;
+typedef CGAL::Triangulation_data_structure_2<Vb>            Tds;
+typedef CGAL::Delaunay_triangulation_2<K,Tds>               Delaunay_triangulation;
+typedef Delaunay_triangulation::Vertex_handle               Vertex_handle;
 
 template <typename V, typename T>
-struct Function_value {
-  typedef V argument_type;
-  typedef std::pair<T, bool> result_type;
-  
-  result_type operator()(const argument_type& a)const
-  {
+struct Function_value
+{
+  typedef V                                                 argument_type;
+  typedef std::pair<T, bool>                                result_type;
+
+  result_type operator()(const argument_type& a) const {
     return result_type(a->info().value, true);
   }
-  
 };
 
-
 template <typename V, typename G>
 struct Function_gradient
-  : public std::iterator<std::output_iterator_tag,void,void,void,void> {
-  
-  typedef V argument_type;
-  typedef std::pair<G,bool> result_type;
-  
-  
-  result_type
-  operator()(const argument_type& a)const
-  {
-    return std::make_pair(a->info().gradient,a->info().gradient != CGAL::NULL_VECTOR) ;
+  : public std::iterator<std::output_iterator_tag, void, void, void, void>
+{
+  typedef V                                                argument_type;
+  typedef std::pair<G, bool>                               result_type;
+
+  result_type operator()(const argument_type& a) const {
+    return std::make_pair(a->info().gradient, a->info().gradient != CGAL::NULL_VECTOR);
   }
 
-
-  const Function_gradient& operator=(const std::pair<V, G>& p) const
-    {
-      p.first->info().gradient = p.second;
-      return *this;
-    }
-
-  const Function_gradient& operator++(int) const
-  {
-    return *this;
-  }
-  
-  const Function_gradient& operator*() const
-  {
+  const Function_gradient& operator=(const std::pair<V, G>& p) const {
+    p.first->info().gradient = p.second;
     return *this;
   }
+
+  const Function_gradient& operator++(int) const { return *this; }
+  const Function_gradient& operator*() const { return *this; }
 };
 
 int main()
 {
   Delaunay_triangulation dt;
 
-  Function_value<Vertex_handle,Coord_type> function_value;
-  Function_gradient<Vertex_handle,Vector> function_gradient;
-  
+  Function_value<Vertex_handle, Coord_type> function_value;
+  Function_gradient<Vertex_handle, Vector> function_gradient;
+
   //parameters for spherical function:
   Coord_type a(0.25), bx(1.3), by(-0.7), c(0.2);
-  for (int y=0 ; y<4 ; y++){
-    for (int x=0 ; x<4 ; x++){
+  for (int y=0 ; y<4 ; y++) {
+    for (int x=0 ; x<4 ; x++) {
       K::Point_2 p(x,y);
       Vertex_handle vh = dt.insert(p);
       Coord_type value = a + bx* x+ by*y + c*(x*x+y*y);
@@ -100,15 +90,14 @@ int main()
                                Traits());
 
   //coordinate computation
-  K::Point_2 p(1.6,1.4);
-  std::vector< std::pair< Vertex_handle, Coord_type > > coords;
-  typedef CGAL::Identity<std::pair< Vertex_handle, Coord_type > > Identity;
+  K::Point_2 p(1.6, 1.4);
+  std::vector<std::pair<Vertex_handle, Coord_type> > coords;
+  typedef CGAL::Identity<std::pair< Vertex_handle, Coord_type> > Identity;
   Coord_type norm = CGAL::natural_neighbor_coordinates_2(dt,
                                                          p,
                                                          std::back_inserter(coords),
                                                          Identity()).second;
 
-
   //Sibson interpolant: version without sqrt:
   std::pair<Coord_type, bool> res = CGAL::sibson_c1_interpolation_square(coords.begin(),
                                                                          coords.end(),
@@ -117,17 +106,16 @@ int main()
                                                                          function_value,
                                                                          function_gradient,
                                                                          Traits());
-  
+
   if(res.second)
     std::cout << "Tested interpolation on " << p
               << " interpolation: " << res.first << " exact: "
-              << a + bx * p.x()+ by * p.y()+ c*(p.x()*p.x()+p.y()*p.y())
+              << a + bx*p.x() + by*p.y() + c*(p.x()*p.x() + p.y()*p.y())
               << std::endl;
   else
     std::cout << "C^1 Interpolation not successful." << std::endl
-              << " not all function_gradients are provided."  << std::endl
+              << " not all function_gradients are provided." << std::endl
               << " You may resort to linear interpolation." << std::endl;
 
-  std::cout << "done" << std::endl;
-  return 0;
+  return EXIT_SUCCESS;
 }
diff --git a/Interpolation/examples/Interpolation/surface_neighbor_coordinates_3.cpp b/Interpolation/examples/Interpolation/surface_neighbor_coordinates_3.cpp
index 0b533a66431..5cd6ba9d563 100644
--- a/Interpolation/examples/Interpolation/surface_neighbor_coordinates_3.cpp
+++ b/Interpolation/examples/Interpolation/surface_neighbor_coordinates_3.cpp
@@ -20,7 +20,7 @@ typedef std::vector< std::pair< Point_3, K::FT > >           Point_coordinate_ve
 
 int main()
 {
-  int n=100;
+  int n = 100;
   std::vector< Point_3> points;
   points.reserve(n);
 
@@ -28,32 +28,32 @@ int main()
   CGAL::Random_points_on_sphere_3<Point_3> g(1);
   CGAL::cpp11::copy_n(g, n, std::back_inserter(points));
 
-  Point_3 p(1, 0,0);
+  Point_3 p(1, 0, 0);
   Vector_3 normal(p - CGAL::ORIGIN);
   std::cout << "Compute surface neighbor coordinates for " << p << std::endl;
   Point_coordinate_vector coords;
-  CGAL::Triple<std::back_insert_iterator<Point_coordinate_vector>,
-               K::FT, bool> result =
+  CGAL::Triple<std::back_insert_iterator<Point_coordinate_vector>, K::FT, bool> result =
     CGAL::surface_neighbor_coordinates_3(points.begin(), points.end(),
                                          p, normal,
                                          std::back_inserter(coords),
                                          K());
-  if(!result.third){
+  if(!result.third)
+  {
     //Undersampling:
     std::cout << "The coordinate computation was not successful." << std::endl;
     return 0;
   }
+
   K::FT norm = result.second;
 
   std::cout << "Testing the barycentric property " << std::endl;
   Point_3 b(0, 0, 0);
   for(std::vector< std::pair< Point_3, Coord_type > >::const_iterator
       it = coords.begin(); it!=coords.end(); ++it)
-    b = b + (it->second/norm)* (it->first - CGAL::ORIGIN);
+    b = b + (it->second/norm) * (it->first - CGAL::ORIGIN);
 
   std::cout << "    weighted barycenter: " << b <<std::endl;
   std::cout << "    squared distance: " << CGAL::squared_distance(p,b) << std::endl;
 
-  std::cout << "done" << std::endl;
-  return 0;
+  return EXIT_SUCCESS;
 }
diff --git a/Interpolation/include/CGAL/Interpolation/internal/helpers.h b/Interpolation/include/CGAL/Interpolation/internal/helpers.h
index 6e5857cb440..4bfc0094045 100644
--- a/Interpolation/include/CGAL/Interpolation/internal/helpers.h
+++ b/Interpolation/include/CGAL/Interpolation/internal/helpers.h
@@ -1,4 +1,4 @@
-// Copyright (c) 2003,2017  INRIA Sophia-Antipolis (France).
+// Copyright (c) 2003, 2017 INRIA Sophia-Antipolis (France).
 // All rights reserved.
 //
 // This file is part of CGAL (www.cgal.org).
@@ -23,134 +23,137 @@
 
 #include <CGAL/license/Interpolation.h>
 
-
 namespace CGAL {
 namespace Interpolation {
 namespace internal {
-  
-  template < class InterpolationTraits>
-  struct V2P
-  {    
-    typedef typename InterpolationTraits::Point_d Point;
-    typedef typename InterpolationTraits::Weighted_point_d Weighted_point;
 
-    V2P(const InterpolationTraits& traits)
-      : traits(traits)
-    {}
-    
-    template <typename VH>
-    const Point& operator()(const VH& vh) const
-    {
-      return traits.construct_point_d_object()(vh->point());
-    }
-   
-    const Point& operator()(const Point& p) const
-    {
-      return p;
-    }
-    
-    Point operator()(const Weighted_point& wp) const
-    {
-      return traits.construct_point_d_object()(wp);
-    }
+template < class InterpolationTraits >
+struct V2P
+{
+  typedef typename InterpolationTraits::Point_d             Point;
+  typedef typename InterpolationTraits::Weighted_point_d    Weighted_point;
 
-  private:
-    InterpolationTraits traits;
-  };
+  V2P(const InterpolationTraits& traits)
+    : traits(traits)
+  {}
 
-
-  template < typename Dt, typename T2>
-  struct Vertex2Point {
-    typedef typename Dt::Vertex_handle Vertex_handle;
-    typedef typename Dt::Point Point;
-    
-    typedef std::pair<Vertex_handle, T2> argument_type;
-    typedef std::pair<Point, T2> result_type;
-    
-    result_type operator()(const argument_type& vp) const
-    {
-      return std::make_pair(vp.first->point(), vp.second);
-    }
-  };
-
-  template < typename Dt, typename T2>
-  struct Vertex2WPoint {
-    typedef typename Dt::Vertex_handle Vertex_handle;
-    typedef typename Dt::Weighted_point Point;
-    
-    typedef std::pair<Vertex_handle, T2> argument_type;
-    typedef std::pair<Point, T2> result_type;
-    
-    result_type operator()(const argument_type& vp) const
-    {
-      return std::make_pair(vp.first->point(), vp.second);
-    }
-  };
-
-  
-  template <typename Dt, typename Map>
-  struct Vertex2Vertex {
-    typedef typename Dt::Vertex_handle Vertex_handle;
-    typedef typename Dt::Geom_traits::FT FT;
-    typedef std::pair<Vertex_handle,FT> argument_type;
-    typedef std::pair<Vertex_handle,FT> result_type;
-    
-    const Map& map;
-    const Dt& dt;
-
-    Vertex2Vertex(const Map& map, const Dt& dt)
-      : map(map), dt(dt)
-    {}
-
-    result_type operator()(const argument_type& vp) const
-    {
-      typename Map::const_iterator it = map.find(vp.first);
-      CGAL_assertion(it != map.end());
-      CGAL_assertion(dt.tds().is_vertex(it->second));
-      return std::make_pair(it->second, vp.second);
-    }
-  };
-  
-  // the struct "Project_vertex_output_iterator"
-  // is used in the (next two) functions
-  // as well as in regular_neighbor_coordinates_2 and
-  // in surface_neighbor_coordinates_3
-  //
-  //projection of iterator over std::pair <Vertex_handle, T>
-  //to iterator over std::pair< Point, T>
-  template < class OutputIterator, class Fct = void>
-  struct Project_vertex_output_iterator
+  template <typename VH>
+  const Point& operator()(const VH& vh) const
   {
-    // this class wraps the OutputIterator with value type
-    // std::pair<Vertex_handle,T>
-    // into an output iterator with value type std::pair<Point, T>
-    // Conditions: OutputIterator has value type std::pair<Vertex_handle, T>
-    //             and Vertex_handle has a function ->point()
-    //             with return type const Point&
+    return traits.construct_point_d_object()(vh->point());
+  }
 
-    OutputIterator _base;
-    Fct fct;
+  const Point& operator()(const Point& p) const
+  {
+    return p;
+  }
 
-    //creation:
-    Project_vertex_output_iterator(OutputIterator o, Fct fct)
-      : _base(o), fct(fct)
-    {}
+  Point operator()(const Weighted_point& wp) const
+  {
+    return traits.construct_point_d_object()(wp);
+  }
 
-    OutputIterator base() {return _base;}
+private:
+  InterpolationTraits traits;
+};
 
-    Project_vertex_output_iterator& operator++(){_base++; return *this;}
-    Project_vertex_output_iterator& operator++(int){_base++; return *this;}
-    Project_vertex_output_iterator& operator*(){return *this;}
 
-    template<class Vertex_pair>
-    Project_vertex_output_iterator&
-    operator=(const Vertex_pair& vp){
-      *_base = fct(vp);
-      return *this;
-    }
-  };
+template < typename Dt, typename T2 >
+struct Vertex2Point
+{
+  typedef typename Dt::Vertex_handle                  Vertex_handle;
+  typedef typename Dt::Point                          Point;
+
+  typedef std::pair<Vertex_handle, T2>                argument_type;
+  typedef std::pair<Point, T2>                        result_type;
+
+  result_type operator()(const argument_type& vp) const
+  {
+    return std::make_pair(vp.first->point(), vp.second);
+  }
+};
+
+
+template < typename Dt, typename T2 >
+struct Vertex2WPoint
+{
+  typedef typename Dt::Vertex_handle                  Vertex_handle;
+  typedef typename Dt::Weighted_point                 Point;
+
+  typedef std::pair<Vertex_handle, T2>                argument_type;
+  typedef std::pair<Point, T2>                        result_type;
+
+  result_type operator()(const argument_type& vp) const
+  {
+    return std::make_pair(vp.first->point(), vp.second);
+  }
+};
+
+
+template < typename Dt, typename Map >
+struct Vertex2Vertex
+{
+  typedef typename Dt::Vertex_handle                 Vertex_handle;
+  typedef typename Dt::Geom_traits::FT               FT;
+  typedef std::pair<Vertex_handle, FT>               argument_type;
+  typedef std::pair<Vertex_handle, FT>               result_type;
+
+  const Map& map;
+  const Dt& dt;
+
+  Vertex2Vertex(const Map& map, const Dt& dt)
+    : map(map), dt(dt)
+  {}
+
+  result_type operator()(const argument_type& vp) const
+  {
+    typename Map::const_iterator it = map.find(vp.first);
+    CGAL_assertion(it != map.end());
+    CGAL_assertion(dt.tds().is_vertex(it->second));
+    return std::make_pair(it->second, vp.second);
+  }
+};
+
+
+// the struct "Project_vertex_output_iterator"
+// is used in the (next two) functions
+// as well as in regular_neighbor_coordinates_2 and
+// in surface_neighbor_coordinates_3
+//
+//projection of iterator over std::pair <Vertex_handle, T>
+//to iterator over std::pair< Point, T>
+template < class OutputIterator, class Fct = void >
+struct Project_vertex_output_iterator
+{
+  // this class wraps the OutputIterator with value type
+  // std::pair<Vertex_handle,T>
+  // into an output iterator with value type std::pair<Point, T>
+  // Conditions: OutputIterator has value type std::pair<Vertex_handle, T>
+  //             and Vertex_handle has a function ->point()
+  //             with return type const Point&
+
+  OutputIterator _base;
+  Fct fct;
+
+  //creation:
+  Project_vertex_output_iterator(OutputIterator o, Fct fct)
+    : _base(o), fct(fct)
+  {}
+
+  OutputIterator base() {return _base;}
+
+  Project_vertex_output_iterator& operator++(){_base++; return *this;}
+  Project_vertex_output_iterator& operator++(int){_base++; return *this;}
+  Project_vertex_output_iterator& operator*(){return *this;}
+
+  template<class Vertex_pair>
+  Project_vertex_output_iterator& operator=(const Vertex_pair& vp)
+  {
+    *_base = fct(vp);
+    return *this;
+  }
+};
 
-  
 } // namespace internal
 } // namespace Interpolation
 } // namespace CGAL
diff --git a/Interpolation/include/CGAL/Voronoi_intersection_2_traits_3.h b/Interpolation/include/CGAL/Voronoi_intersection_2_traits_3.h
index 3d1ecae2cab..57fdfb1ff47 100644
--- a/Interpolation/include/CGAL/Voronoi_intersection_2_traits_3.h
+++ b/Interpolation/include/CGAL/Voronoi_intersection_2_traits_3.h
@@ -163,16 +163,16 @@ public:
 
   Comparison_result operator()(const Point& p, const Point& q) const
   {
-    if(normal.x()!=Coord_type(0))
+    if(normal.x() != Coord_type(0))
       return (Comparison_result) CGAL_NTS
-          sign(Vector(normal.y(),-normal.x(),Coord_type(0))*(p-q));
-    if(normal.y()!= Coord_type(0))
+          sign(Vector(normal.y(), -normal.x(), Coord_type(0))*(p-q));
+    if(normal.y() != Coord_type(0))
       return (Comparison_result) CGAL_NTS
-          sign(Vector(-normal.y(),normal.x(),Coord_type(0))*(p-q));
+          sign(Vector(-normal.y(), normal.x(), Coord_type(0))*(p-q));
 
-    CGAL_assertion(normal.z()!= Coord_type(0));
+    CGAL_assertion(normal.z() != Coord_type(0));
     return (Comparison_result) CGAL_NTS
-        sign(Vector(-normal.z(),Coord_type(0),normal.x())*(p-q));
+        sign(Vector(-normal.z(), Coord_type(0), normal.x())*(p-q));
   }
 
 private:
@@ -195,16 +195,16 @@ public:
 
   Comparison_result operator()(const Point& p, const Point& q) const
   {
-    if(normal.x()!=Coord_type(0))
+    if(normal.x() != Coord_type(0))
       return (Comparison_result) CGAL_NTS
-          sign(Vector(normal.z(),Coord_type(0),-normal.x())*(p-q));
-    if(normal.y()!= Coord_type(0))
+          sign(Vector(normal.z(), Coord_type(0), -normal.x())*(p-q));
+    if(normal.y() != Coord_type(0))
       return (Comparison_result) CGAL_NTS
-          sign(Vector(Coord_type(0),normal.z(),-normal.y())*(p-q));
+          sign(Vector(Coord_type(0), normal.z(), -normal.y())*(p-q));
 
-    CGAL_assertion(normal.z()!= Coord_type(0));
+    CGAL_assertion(normal.z() != Coord_type(0));
     return (Comparison_result) CGAL_NTS
-        sign(Vector(Coord_type(0),-normal.z(),normal.y())*(p-q));
+        sign(Vector(Coord_type(0), -normal.z(), normal.y())*(p-q));
   }
 
 private:
diff --git a/Interpolation/include/CGAL/natural_neighbor_coordinates_2.h b/Interpolation/include/CGAL/natural_neighbor_coordinates_2.h
index 1edd26c347f..3ec4b39889f 100644
--- a/Interpolation/include/CGAL/natural_neighbor_coordinates_2.h
+++ b/Interpolation/include/CGAL/natural_neighbor_coordinates_2.h
@@ -22,19 +22,21 @@
 #define CGAL_NATURAL_NEIGHBOR_COORDINATES_2_H
 
 #include <CGAL/license/Interpolation.h>
+
 #include <CGAL/Interpolation/internal/helpers.h>
+
 #include <CGAL/Iterator_project.h>
 #include <CGAL/Polygon_2.h>
 #include <CGAL/number_utils_classes.h>
 #include <CGAL/utility.h>
 
+#include <iterator>
 #include <list>
 #include <utility>
+#include <vector>
 
 namespace CGAL {
 
-
-  
 // The following natural_neighbor_coordinate_2 functions fix the
 // traits class to be Dt::Geom_traits. The following signatures could
 // be used if one wants to pass a traits class as argument:
@@ -60,8 +62,8 @@ template < class Dt, class OutputIterator >
 Triple< OutputIterator, typename Dt::Geom_traits::FT, bool >
 natural_neighbors_2(const Dt& dt,
                     const typename Dt::Geom_traits::Point_2& p,
-                    OutputIterator out, typename Dt::Face_handle start
-                    = typename Dt::Face_handle())
+                    OutputIterator out,
+                    typename Dt::Face_handle start = typename Dt::Face_handle())
 {
   typedef typename Dt::Geom_traits       Traits;
   typedef typename Traits::FT            Coord_type;
@@ -96,13 +98,13 @@ natural_neighbors_2(const Dt& dt,
     Equal_x_2 equal_x_2;
     if(!equal_x_2(p1,p2))
     {
-      coef1 =  (p.x() - p2.x()) / (p1.x() - p2.x());
+      coef1 = (p.x() - p2.x()) / (p1.x() - p2.x());
       coef2 = 1 - coef1;
       *out++ = std::make_pair(v1,coef1);
       *out++ = std::make_pair(v2,coef2);
     } else {
       coef1 = (p.y() - p2.y()) / (p1.y() - p2.y());
-      coef2 = 1-coef1;
+      coef2 = 1 - coef1;
       *out++ = std::make_pair(v1,coef1);
       *out++ = std::make_pair(v2,coef2);
     }
@@ -122,7 +124,7 @@ natural_neighbors_2(const Dt& dt,
   return natural_neighbors_2(dt, p, out, hole.begin(), hole.end());
 }
 
-  
+
 //function call if the conflict zone is known:
 // OutputIterator has value type
 //        std::pair<Dt::Vertex_handle, Dt::Geom_traits::FT>
@@ -148,6 +150,7 @@ natural_neighbors_2(const Dt& dt,
   Coord_type area_sum(0);
   EdgeIterator hit = hole_end;
   --hit;
+
   //in the beginning: prev is the "last" vertex of the hole:
   // later: prev is the last vertex processed (previously)
   Vertex_handle prev = hit->first->vertex(dt.cw(hit->second));
@@ -190,7 +193,7 @@ natural_neighbors_2(const Dt& dt,
   return make_triple(out, area_sum, true);
 }
 
-  
+
 template < class Dt, class OutputIterator, class Fct >
 Triple< OutputIterator, typename Dt::Geom_traits::FT, bool >
 natural_neighbor_coordinates_2(const Dt& dt,
@@ -205,8 +208,7 @@ natural_neighbor_coordinates_2(const Dt& dt,
   Interpolation::internal::Project_vertex_output_iterator<OutputIterator, Fct> op(out, fct);
 
   Triple<Interpolation::internal::Project_vertex_output_iterator<OutputIterator,Fct>,
-         typename Dt::Geom_traits::FT, bool > result =
-      natural_neighbors_2(dt, p, op, start);
+         typename Dt::Geom_traits::FT, bool > result = natural_neighbors_2(dt, p, op, start);
 
   return make_triple(result.first.base(), result.second, result.third);
 }
@@ -255,6 +257,7 @@ natural_neighbor_coordinates_2(const Dt& dt,
   return make_triple(result.first.base(), result.second, result.third);
 }
 
+
 //OutputIterator has value type
 //   std::pair< Dt::Geom_traits::Point_2, Dt::Geom_traits::FT>
 //function call if the conflict zone is known:
@@ -269,7 +272,8 @@ natural_neighbor_coordinates_2(const Dt& dt,
   return natural_neighbor_coordinates_2(dt,p, out, hole_begin, hole_end,
                                         Interpolation::internal::Vertex2Point<Dt, typename Dt::Geom_traits::FT>());
 }
-  
+
+
 /**********************************************************/
 //compute the coordinates for a vertex of the triangulation
 // with respect to the other points in the triangulation
@@ -295,11 +299,11 @@ natural_neighbor_coordinates_2(const Dt& dt,
 
   typedef std::map<Vertex_handle /*t2*/, Vertex_handle/*dt*/> V2V;
   V2V correspondence_map;
-  
+
   do{
     CGAL_assertion(!dt.is_infinite(vc));
     Vertex_handle vh2 = t2.insert(vc->point());
-    correspondence_map[vh2] = vc; 
+    correspondence_map[vh2] = vc;
   }
   while(++vc!=done);
 
@@ -308,7 +312,7 @@ natural_neighbor_coordinates_2(const Dt& dt,
   return natural_neighbor_coordinates_2(t2, vh->point(), out, cfct);
 }
 
-  
+
 template <class Dt, class OutputIterator>
 Triple< OutputIterator, typename Dt::Geom_traits::FT, bool >
 natural_neighbor_coordinates_2(const Dt& dt,
@@ -340,7 +344,6 @@ public:
   }
 };
 
- 
 } //namespace CGAL
 
 #endif // CGAL_NATURAL_NEIGHBOR_COORDINATES_2_H
diff --git a/Interpolation/include/CGAL/natural_neighbor_coordinates_3.h b/Interpolation/include/CGAL/natural_neighbor_coordinates_3.h
index cacc1ec9f98..18b92a32c76 100644
--- a/Interpolation/include/CGAL/natural_neighbor_coordinates_3.h
+++ b/Interpolation/include/CGAL/natural_neighbor_coordinates_3.h
@@ -72,8 +72,8 @@ construct_circumcenter(const typename DT::Facet& f,
 
 template <class Dt, class OutputIterator>
 Triple< OutputIterator,  // iterator with value type std::pair<Dt::Vertex_handle, Dt::Geom_traits::FT>
-  typename Dt::Geom_traits::FT,  // Should provide 0 and 1
-  bool >
+        typename Dt::Geom_traits::FT,  // Should provide 0 and 1
+        bool >
 laplace_natural_neighbor_coordinates_3(const Dt& dt,
                                        const typename Dt::Geom_traits::Point_3& Q,
                                        OutputIterator nn_out,
@@ -108,23 +108,25 @@ laplace_natural_neighbor_coordinates_3(const Dt& dt,
   std::set<Cell_handle> cells;
   // To replace the forbidden access to the "in conflict" flag :
   // std::find operations on this set
-  std::vector<Facet> bound_facets; bound_facets.reserve(32);
-  typename std::vector<Facet>::iterator bound_it;
+  std::vector<Facet> bound_facets;
+  bound_facets.reserve(32);
+
   // Find the cells in conflict with Q
   dt.find_conflicts(Q, c,
                     std::back_inserter(bound_facets),
-                    std::inserter(cells,cells.begin()));
+                    std::inserter(cells, cells.begin()));
 
   std::map<Vertex_handle,Coord_type> coordinate;
   typename std::map<Vertex_handle,Coord_type>::iterator coor_it;
 
+  typename std::vector<Facet>::iterator bound_it;
   for (bound_it = bound_facets.begin(); bound_it != bound_facets.end(); ++bound_it)
   {
     //for each facet on the boundary
     Facet f1 = *bound_it;
     Cell_handle cc1 = f1.first;
     if (dt.is_infinite(cc1))
-      return make_triple(nn_out,norm_coeff=Coord_type(1), false);//point outside the convex-hull
+      return make_triple(nn_out, norm_coeff=Coord_type(1), false);//point outside the convex-hull
 
     CGAL_triangulation_assertion_code(Cell_handle cc2 = cc1->neighbor(f1.second);)
     CGAL_triangulation_assertion(std::find(cells.begin(),cells.end(),cc1) != cells.end());//TODO : Delete
@@ -172,8 +174,8 @@ laplace_natural_neighbor_coordinates_3(const Dt& dt,
 
 template <class Dt, class OutputIterator>
 Triple< OutputIterator,  // iterator with value type std::pair<Dt::Vertex_handle, Dt::Geom_traits::FT>
-  typename Dt::Geom_traits::FT,  // Should provide 0 and 1
-  bool >
+        typename Dt::Geom_traits::FT,  // Should provide 0 and 1
+        bool >
 sibson_natural_neighbor_coordinates_3(const Dt& dt,
                                       const typename Dt::Geom_traits::Point_3& Q,
                                       OutputIterator nn_out,
diff --git a/Interpolation/include/CGAL/regular_neighbor_coordinates_2.h b/Interpolation/include/CGAL/regular_neighbor_coordinates_2.h
index 33d332d31e6..3818525661e 100644
--- a/Interpolation/include/CGAL/regular_neighbor_coordinates_2.h
+++ b/Interpolation/include/CGAL/regular_neighbor_coordinates_2.h
@@ -120,8 +120,7 @@ regular_neighbor_coordinates_vertex_2(const Rt& rt,
 
   rt.get_boundary_of_conflicts_and_hidden_vertices(p,
                                                    std::back_inserter(hole),
-                                                   std::back_inserter
-                                                   (hidden_vertices),
+                                                   std::back_inserter(hidden_vertices),
                                                    fh);
   return regular_neighbor_coordinates_vertex_2(rt, p, out, vor_vertices,
                                                hole.begin(),hole.end(),
@@ -171,8 +170,8 @@ regular_neighbor_coordinates_vertex_2(const Rt& rt,
   typedef typename Traits::FT              Coord_type;
   typedef typename Rt::Bare_point          Bare_point;
 
-  typedef typename Rt::Vertex_handle     Vertex_handle;
-  typedef typename Rt::Face_circulator   Face_circulator;
+  typedef typename Rt::Vertex_handle       Vertex_handle;
+  typedef typename Rt::Face_circulator     Face_circulator;
 
   //no hole because only (exactly!) one vertex is hidden:
   if(hole_begin==hole_end){
@@ -193,54 +192,55 @@ regular_neighbor_coordinates_vertex_2(const Rt& rt,
   Vertex_handle prev = hit->first->vertex(rt.cw(hit->second));
   hit = hole_begin;
   while(hit != hole_end)
+  {
+    Coord_type area(0);
+    Vertex_handle current = hit->first->vertex(rt.cw(hit->second));
+
+    //a first Voronoi vertex of the cell of p:
+    vor[0] = rt.geom_traits().construct_weighted_circumcenter_2_object()
+             (current->point(),
+              hit->first->vertex(rt.ccw(hit->second))->point(), p);
+    *vor_vertices++= vor[0];
+
+    //triangulation of the Voronoi subcell:
+    //a second vertex as base
+    Face_circulator fc = rt.incident_faces(current, hit->first);
+    ++fc;
+    vor[1] = rt.dual(fc);
+    // iteration over all other "old" Voronoi vertices
+    while(!fc->has_vertex(prev))
     {
-      Coord_type area(0);
-      Vertex_handle current = hit->first->vertex(rt.cw(hit->second));
-
-      //a first Voronoi vertex of the cell of p:
-      vor[0] = rt.geom_traits().construct_weighted_circumcenter_2_object()
-               (current->point(),
-                hit->first->vertex(rt.ccw(hit->second))->point(), p);
-      *vor_vertices++= vor[0];
-
-      //triangulation of the Voronoi subcell:
-      //a second vertex as base
-      Face_circulator fc = rt.incident_faces(current, hit->first);
       ++fc;
-      vor[1] = rt.dual(fc);
-      // iteration over all other "old" Voronoi vertices
-      while(!fc->has_vertex(prev))
-      {
-        ++fc;
-        vor[2] = rt.dual(fc);
-
-        area += polygon_area_2(vor.begin(), vor.end(), rt.geom_traits());
-        vor[1] = vor[2];
-      }
-      //the second Voronoi vertex of the cell of p:
-      vor[2] =
-          rt.geom_traits().construct_weighted_circumcenter_2_object()
-          (prev->point(),current->point(),p);
-      *vor_vertices++= vor[2];
+      vor[2] = rt.dual(fc);
 
       area += polygon_area_2(vor.begin(), vor.end(), rt.geom_traits());
-      *out++= std::make_pair(current,area);
-
-      area_sum += area;
-
-      //update prev and hit:
-      prev= current;
-      ++hit;
+      vor[1] = vor[2];
     }
 
+    //the second Voronoi vertex of the cell of p:
+    vor[2] =
+        rt.geom_traits().construct_weighted_circumcenter_2_object()
+        (prev->point(),current->point(),p);
+    *vor_vertices++= vor[2];
+
+    area += polygon_area_2(vor.begin(), vor.end(), rt.geom_traits());
+    *out++= std::make_pair(current,area);
+
+    area_sum += area;
+
+    //update prev and hit:
+    prev= current;
+    ++hit;
+  }
+
   //get coordinate for hidden vertices
   //                   <=> the area of their Voronoi cell.
   //decomposition of the cell into triangles
   //        vor1: dual of first triangle
   //        vor2, vor 3: duals of two consecutive triangles
   Face_circulator fc, fc_begin;
-  for(; hidden_vertices_begin != hidden_vertices_end;
-      ++hidden_vertices_begin){
+  for(; hidden_vertices_begin != hidden_vertices_end; ++hidden_vertices_begin)
+  {
     Coord_type area(0);
     fc_begin = rt.incident_faces(*hidden_vertices_begin);
     vor[0] = rt.dual(fc_begin);
@@ -248,7 +248,8 @@ regular_neighbor_coordinates_vertex_2(const Rt& rt,
     ++fc;
     vor[1] = rt.dual(fc);
     ++fc;
-    while(fc != fc_begin){
+    while(fc != fc_begin)
+    {
       vor[2] = rt.dual(fc);
       area += polygon_area_2(vor.begin(), vor.end(), rt.geom_traits());
 
@@ -263,7 +264,7 @@ regular_neighbor_coordinates_vertex_2(const Rt& rt,
   return make_triple(out, area_sum, true);
 }
 
-  
+
 ////////////////////////////////////////////////////////////
 //the cast from vertex to point:
 // the following functions return an Output_iterator over
@@ -289,7 +290,7 @@ regular_neighbor_coordinates_2(const Rt& rt,
                                Fct fct,
                                typename boost::disable_if_c<
                                  boost::is_convertible<Fct,
-                                                       typename Rt::Face_handle>::value 
+                                                       typename Rt::Face_handle>::value
                                >::type* = 0)
 {
   return regular_neighbor_coordinates_2(rt, p, out, fct, typename Rt::Face_handle());
@@ -346,7 +347,7 @@ regular_neighbor_coordinates_2(const Rt& rt,
   return make_triple(result.first.base(), result.second, result.third);
 }
 
-  
+
 template <class Rt, class OutputIterator, class OutputIteratorVorVertices>
 Triple< OutputIterator, typename Rt::Geom_traits::FT, bool >
 regular_neighbor_coordinates_2(const Rt& rt,
@@ -361,7 +362,7 @@ regular_neighbor_coordinates_2(const Rt& rt,
                                  vor_vertices, start);
 }
 
-  
+
 //OutputIterator has value type
 //   std::pair< Rt::Geom_traits::Weighted_point_2, Rt::Geom_traits::FT>
 template <class Rt, class OutputIterator, class EdgeIterator,
@@ -410,7 +411,7 @@ regular_neighbor_coordinates_2(const Rt& rt,
   return make_triple(result.first.base(), result.second, result.third);
 }
 
-  
+
 template <class Rt, class OutputIterator, class EdgeIterator,
   class VertexIterator , class OutputIteratorVorVertices >
 Triple< OutputIterator, typename Rt::Geom_traits::FT, bool >
@@ -454,11 +455,11 @@ regular_neighbor_coordinates_2(const Rt& rt,
 
   typedef std::map<Vertex_handle /*t2*/, Vertex_handle/*dt*/> V2V;
   V2V correspondence_map;
-  
+
   do{
     CGAL_assertion(!rt.is_infinite(vc));
     Vertex_handle vh2 = t2.insert(vc->point());
-    correspondence_map[vh2] = vc; 
+    correspondence_map[vh2] = vc;
   }
   while(++vc!=done);
 
@@ -468,7 +469,7 @@ regular_neighbor_coordinates_2(const Rt& rt,
   return regular_neighbor_coordinates_2(t2, vh->point(), out, cfct);
 }
 
-  
+
 template <class Rt, class OutputIterator>
 Triple< OutputIterator, typename Rt::Geom_traits::FT, bool >
 regular_neighbor_coordinates_2(const Rt& rt,
@@ -481,7 +482,7 @@ regular_neighbor_coordinates_2(const Rt& rt,
                                         Interpolation::internal::Vertex2Point<Rt, typename Rt::Geom_traits::FT>());
 }
 
-  
+
 //class providing a function object:
 //OutputIterator has value type
 //   std::pair< Rt::Geom_traits::Weighted_point_2, Rt::Geom_traits::FT>
diff --git a/Interpolation/include/CGAL/sibson_gradient_fitting.h b/Interpolation/include/CGAL/sibson_gradient_fitting.h
index 8f2d60e80be..536b1853ce6 100644
--- a/Interpolation/include/CGAL/sibson_gradient_fitting.h
+++ b/Interpolation/include/CGAL/sibson_gradient_fitting.h
@@ -82,7 +82,7 @@ sibson_gradient_fitting(ForwardIterator first,
   return Hn.inverse().transform(pn);
 }
 
-  
+
 template < class Triangul, class ForwardIterator, class Functor, class Traits, class VH>
 typename Traits::Vector_d
 sibson_gradient_fitting(const Triangul& tr,
@@ -98,9 +98,9 @@ sibson_gradient_fitting(const Triangul& tr,
 {
   const typename Traits::Point_d& bare_p = traits.construct_point_d_object()(vh->point());
   typename Functor::result_type fn = function_value(bare_p);
- 
+
   CGAL_assertion(fn.second);
-                                                         
+
   return sibson_gradient_fitting(first,
                                  beyond,
                                  norm,
@@ -109,7 +109,7 @@ sibson_gradient_fitting(const Triangul& tr,
                                  function_value,
                                  traits);
 }
- 
+
 template < class Triangul, class ForwardIterator, class Functor, class Traits, class VH>
 typename Traits::Vector_d
 sibson_gradient_fitting(const Triangul& tr,
@@ -126,7 +126,7 @@ sibson_gradient_fitting(const Triangul& tr,
   const typename Traits::Point_d& bare_p = traits.construct_point_d_object()(vh->point());
   typename Functor::result_type fn = function_value(vh->point());
   CGAL_assertion(fn.second);
-                                                         
+
   return sibson_gradient_fitting(first,
                                  beyond,
                                  norm,
@@ -135,7 +135,7 @@ sibson_gradient_fitting(const Triangul& tr,
                                  function_value,
                                  traits);
 }
-  
+
 template < class Triangul, class ForwardIterator, class Functor, class Traits, class VH>
 typename Traits::Vector_d
 sibson_gradient_fitting(const Triangul& tr,
@@ -152,7 +152,7 @@ sibson_gradient_fitting(const Triangul& tr,
   const typename Traits::Point_d& bare_p = traits.construct_point_d_object()(vh->point());
   typename Functor::result_type fn = function_value(vh);
    CGAL_assertion(fn.second);
-                                                         
+
   return sibson_gradient_fitting(first,
                                  beyond,
                                  norm,
@@ -162,7 +162,7 @@ sibson_gradient_fitting(const Triangul& tr,
                                  traits);
 }
 
-  
+
 template < class Triangul, class OutputIterator, class Functor,
            class CoordFunctor, class OIF, class Traits>
 OutputIterator
@@ -197,14 +197,14 @@ sibson_gradient_fitting_internal(const Triangul& tr,
                                                           function_value,
                                                           traits,
                                                           typename Functor::argument_type())));
-      
+
       coords.clear();
     }
   }
   return out;
 }
 
-  
+
 //the following functions allow to fit the gradients for all points in
 // a triangulation except the convex hull points.
 // -> _nn2: natural_neighbor_coordinates_2
@@ -242,7 +242,7 @@ template < class Dt, class OutputIterator, class Functor, class Traits>
 OutputIterator
 sibson_gradient_fitting_nn_2(const Dt& dt,
                              OutputIterator out,
-                             Functor function_value,    
+                             Functor function_value,
                              const Traits& traits)
 {
   return sibson_gradient_fitting_nn_2(dt, out, function_value,