diff --git a/Packages/Interpolation/doc_tex/Interpolation_ref/Data_access.tex b/Packages/Interpolation/doc_tex/Interpolation_ref/Data_access.tex
new file mode 100644
index 00000000000..2bf9d2d0029
--- /dev/null
+++ b/Packages/Interpolation/doc_tex/Interpolation_ref/Data_access.tex
@@ -0,0 +1,64 @@
+% +------------------------------------------------------------------------+
+% | Reference manual page: Data_access.tex
+% +------------------------------------------------------------------------+
+% | 
+% | August 2003  Julia Floetotto
+% | Package:   Interpolation
+% | 
+% |
+% +------------------------------------------------------------------------+
+\begin{ccRefClass}{Data_access<Map>}  %% add template arg's if necessary
+\ccIndexMainItemBegin{Data_access<Map>}
+
+\ccDefinition
+  
+The struct \ccRefName\ implements a functor that allows to retrieve
+data from a associative container. The functor keeps a reference to
+the container. Given an instance of the container's key type, it
+returns a pair of the container's value type and a boolean indicating
+whether the retrieval was successful.\\
+
+This class can be used to provide the function\_values and
+function\_gradients of the interpolation functions.
+
+\ccInclude{CGAL/interpolation_functions.h} 
+
+\ccParameters The class
+\ccRefName\ has the container type \ccc{Map} as template parameter.
+
+\ccTypes 
+\ccTypedef{typedef Map::mapped_type Data_type;}{}
+\ccGlue
+\ccTypedef{typedef Map::key_type Key_type;}{}
+\ccCreation
+\ccCreationVariable{data_access}  %% choose variable name
+\ccConstructor{Data_acces(const Map& map);}
+{Introduces a \ccc{Data_access} to the container \ccc{map}.}
+
+
+\ccMethod{std::pair< Data_type, bool> operator()(const Key_type& p);} {If
+  there is an entry for \ccc{p} in the container \ccc{map}, then the
+  pair of \ccc{map.find(p)} and \ccc{true} is returned. Otherwise, the
+  boolean value of the pair is \ccc{false}.}
+
+\ccSeeAlso
+\ccRefIdfierPage{CGAL::linear_interpolation}\\
+\ccRefIdfierPage{CGAL::sibson_c1_interpolation} \\
+\ccRefIdfierPage{CGAL::sibson_gradient_fitting} \\
+\ccRefIdfierPage{CGAL::farin_c1_interpolation} \\
+\ccRefIdfierPage{CGAL::quadratic_interpolation} \\
+%\ccRefIdfierPage{CGAL::natural_neighbor_coordinates_2}\\
+%\ccRefIdfierPage{CGAL::regular_neighbor_coordinates_2} \\
+%\ccRefIdfierPage{CGAL::surface_neighbor_coordinates_2_3}\\
+%\ccImplementation 
+
+\ccIndexSubitemEnd{Interpolation}{Data_access}
+\end{ccRefClass}
+%\renewcommand{\ccRefPageBegin}{}
+%\renewcommand{\ccRefPageEnd}{}
+
+% +------------------------------------------------------------------------+
+%%RefPage: end of main body, begin of footer
+% EOF
+% +------------------------------------------------------------------------+
+
diff --git a/Packages/Interpolation/doc_tex/Interpolation_ref/InterpolationTraits.tex b/Packages/Interpolation/doc_tex/Interpolation_ref/InterpolationTraits.tex
index bb45b83d100..72475e48ec5 100644
--- a/Packages/Interpolation/doc_tex/Interpolation_ref/InterpolationTraits.tex
+++ b/Packages/Interpolation/doc_tex/Interpolation_ref/InterpolationTraits.tex
@@ -1,5 +1,5 @@
 % +------------------------------------------------------------------------+
-% | Reference manual page: Interpolation_traits.tex
+% | Reference manual page: InterpolationTraits.tex
 % +------------------------------------------------------------------------+
 % | 06.04.2000   Author
 % | Package: Package
@@ -24,7 +24,7 @@ defines the primitives used in the algorithms.  The concept
   \ccc{FieldNumberType}.}  \ccNestedType{Point_d}{The point type on
   which the function is defined and interpolated.}
 \ccNestedType{Vector_d}{The corresponding vector type.}
-\ccNestedType{Construct_vector} {A constructor object for
+\ccNestedType{Construct_vector_d} {A constructor object for
   \ccc{Vector_d}.
   Provides~: \\
   \ccc{Vector_d operator() (Point_d a, Point_d b)} which produces the
@@ -32,40 +32,28 @@ defines the primitives used in the algorithms.  The concept
     a} and \\
   \ccc{Vector_d operator() (Null_vector NULL_VECTOR)} which introduces
   the null vector.}
-\ccNestedType{Construct_scaled_vector}{Constructor object for
+\ccNestedType{Construct_scaled_vector_d}{Constructor object for
   \ccc{Vector_d}.
   Provides~: \\
   \ccc{Vector_d operator() (Vector_d v,FT scale)} which produces the
   vector \ccc{v} scaled by a factor \ccc{scale}.}
-\ccNestedType{Compute_squared_distance} {Constructor
+\ccNestedType{Compute_squared_distance_d} {Constructor
   object for \ccc{FT}. Provides the operator:\\
   \ccc{FT operator() (Point_d a, Point_d b)} returning the squared
   distance between \ccc{a} and \ccc{b}.}
 
 %\ccCreation
 \ccCreationVariable{traits}  %% choose variable name
-%A default constructor and a copy constructor
-% and an assignement operator are required. 
-%Note that further constructors
-%can be provided. 
-%%\ccThree{InterpolationTraits}{traits=gtrxx  }{}
 \ccConstructor{InterpolationTraits();}{default constructor.}
-%\ccGlue
-%\ccConstructor{InterpolationTraits(InterpolationTraits gtr);}
-%{Copy constructor}
-%\ccMethod{InterpolationTraits operator=(InterpolationTraits gtr);}
-%{Assignment operator.}
 
-%\ccOperations 
 
 \ccHeading{Construction objects} The following functions that
 create instances of the above constructor object types must exist.
-
-\ccMethod{Construct_vector construct_vector_object();}{}
+\ccMethod{Construct_vector_d construct_vector_d_object();}{}
 \ccGlue
-\ccMethod{Construct_scaled_vector construct_scaled_vector_object();}{}
+\ccMethod{Construct_scaled_vector_d construct_scaled_vector_d_object();}{}
 \ccGlue
-\ccMethod{Compute_squared_distance compute_squared_distance_object();}{}
+\ccMethod{Compute_squared_distance_d compute_squared_distance_d_object();}{}
 
 \ccHasModels
 \ccc{CGAL::Interpolation_traits_2<K>} \\
@@ -78,7 +66,7 @@ create instances of the above constructor object types must exist.
 \ccRefIdfierPage{CGAL::quadratic_interpolation} \\
 \ccRefIdfierPage{CGAL::natural_neighbor_coordinates_2}\\
 \ccRefIdfierPage{CGAL::regular_neighbor_coordinates_2} \\
-\ccRefIdfierPage{CGAL::surface_neighbor_coordinates_2_3}\\
+\ccRefIdfierPage{CGAL::surface_neighbor_coordinates_3}\\
 \end{ccRefConcept}
 
 % +------------------------------------------------------------------------+
diff --git a/Packages/Interpolation/doc_tex/Interpolation_ref/Interpolation_gradient_fitting_traits_2.tex b/Packages/Interpolation/doc_tex/Interpolation_ref/Interpolation_gradient_fitting_traits_2.tex
index 6a7a7099c90..aeab07aabe9 100644
--- a/Packages/Interpolation/doc_tex/Interpolation_ref/Interpolation_gradient_fitting_traits_2.tex
+++ b/Packages/Interpolation/doc_tex/Interpolation_ref/Interpolation_gradient_fitting_traits_2.tex
@@ -32,22 +32,22 @@
   \ccIndexSubitem[c]{InterpolationTraits}{model} \\
   \ccIndexSubitem[c]{GradientFittingTraits}{model} \\
    \ccTypes
-  \ccSetThreeColumns{typedef R::Less_signed_distance_to_line_2 }{Less_signed_distance_to_line_2}{}
+  \ccSetThreeColumns{typedef K::Less_signed_distance_to_line_2 }{Less_signed_distance_to_line_2}{}
   \ccThreeToTwo 
-  \ccTypedef{typedef R::FT                            FT;}{}
+  \ccTypedef{typedef K::FT                            FT;}{}
   \ccGlue
-  \ccTypedef{typedef R::Point_2                        Point_d;}{}
+  \ccTypedef{typedef K::Point_2                        Point_d;}{}
   \ccGlue
-  \ccTypedef{typedef R::Vector_2                       Vector_d;}{}
+  \ccTypedef{typedef K::Vector_2                       Vector_d;}{}
   \ccGlue
-\ccTypedef{typedef R::Aff_transformation_2         Aff_transformation_d;}{}
+\ccTypedef{typedef K::Aff_transformation_2         Aff_transformation_d;}{}
 \ccGlue
-  \ccTypedef{typedef R::Construct_vector_2     Construct_vector_d;}{}
+  \ccTypedef{typedef K::Construct_vector_2     Construct_vector_d;}{}
   \ccGlue
-  \ccTypedef{typedef R::Construct_scaled_vector_2 
+  \ccTypedef{typedef K::Construct_scaled_vector_2 
     Construct_scaled_vector_d;}{}
   \ccGlue
-  \ccTypedef{typedef R::Compute_squared_distance_2 
+  \ccTypedef{typedef K::Compute_squared_distance_2 
     Compute_squared_distance_d;}{}
   \ccGlue
   \ccTypedef{typedef Construct_null_matrix_2<Aff_transformation_d>   
@@ -58,7 +58,7 @@
   \ccGlue
   \ccTypedef{typedef Construct_sum_matrix_2<Aff_transformation_d> Construct_sum_matrix_d;}{}
   \ccGlue
-  \ccTypedef{typedef Construct_outer_product_2<Rep>             Construct_outer_product_d;}{}
+  \ccTypedef{typedef Construct_outer_product_2<K>             Construct_outer_product_d;}{}
                       
 \ccOperations
 \ccAutoIndexingOff
diff --git a/Packages/Interpolation/doc_tex/Interpolation_ref/Interpolation_traits_2.tex b/Packages/Interpolation/doc_tex/Interpolation_ref/Interpolation_traits_2.tex
index 6058eaeb091..622c0327f3c 100644
--- a/Packages/Interpolation/doc_tex/Interpolation_ref/Interpolation_traits_2.tex
+++ b/Packages/Interpolation/doc_tex/Interpolation_ref/Interpolation_traits_2.tex
@@ -11,7 +11,7 @@
 % +------------------------------------------------------------------------+
 
 
-\begin{ccRefClass}{Interpolation_traits_2<R>}
+\begin{ccRefClass}{Interpolation_traits_2<K>}
   
 
 
@@ -24,7 +24,7 @@
   \ccc{InterpolationTraits} and can be used to instantiate the
   geometric traits class of interpolation methods applied on a
   bivariate function over a two-dimensional domain. The traits class
-  is templated by a kernel class \ccc{R}.
+  is templated by a kernel class \ccc{K}.
   
   \ccInclude{CGAL/Interpolation_traits_2.h}
 
@@ -33,20 +33,19 @@
   \ccIndexSubitem[c]{InterpolationTraits}{model} \\
   
   \ccTypes
-  \ccSetThreeColumns{typedef R::Less_signed_distance_to_line_2 }{Less_signed_distance_to_line_2}{}
   \ccThreeToTwo 
-  \ccTypedef{typedef R::FT                            FT;}{}
+  \ccTypedef{typedef K::FT                            FT;}{}
   \ccGlue
-  \ccTypedef{typedef R::Point_2                        Point_d;}{}
+  \ccTypedef{typedef K::Point_2                        Point_d;}{}
   \ccGlue
-  \ccTypedef{typedef R::Vector_2                       Vector_d;}{}
+  \ccTypedef{typedef K::Vector_2                       Vector_d;}{}
   \ccGlue
-  \ccTypedef{typedef R::Construct_vector_2     Construct_vector_d;}{}
+  \ccTypedef{typedef K::Construct_vector_2     Construct_vector_d;}{}
   \ccGlue
-  \ccTypedef{typedef R::Construct_scaled_vector_2 
+  \ccTypedef{typedef K::Construct_scaled_vector_2 
     Construct_scaled_vector_d;}{}
   \ccGlue
-  \ccTypedef{typedef R::Compute_squared_distance_2 
+  \ccTypedef{typedef K::Compute_squared_distance_2 
     Compute_squared_distance_d;}{}
                                            
 \ccOperations
diff --git a/Packages/Interpolation/doc_tex/Interpolation_ref/farin_c1_interpolation.tex b/Packages/Interpolation/doc_tex/Interpolation_ref/farin_c1_interpolation.tex
index bbcb967ed9f..15dde9d6e64 100644
--- a/Packages/Interpolation/doc_tex/Interpolation_ref/farin_c1_interpolation.tex
+++ b/Packages/Interpolation/doc_tex/Interpolation_ref/farin_c1_interpolation.tex
@@ -55,7 +55,7 @@ to provide the square root operation.
 \ccRefIdfierPage{CGAL::Interpolation_traits_2} \\
 \ccRefIdfierPage{CGAL::natural_neighbor_coordinates_2}\\
 \ccRefIdfierPage{CGAL::regular_neighbor_coordinates_2} \\
-\ccRefIdfierPage{CGAL::surface_neighbor_coordinates_2_3}\\
+\ccRefIdfierPage{CGAL::surface_neighbor_coordinates_3}\\
 %\ccImplementation
 
 \ccIndexSubitemEnd{Interpolation}{farin_c1_interpolation}
diff --git a/Packages/Interpolation/doc_tex/Interpolation_ref/intro.tex b/Packages/Interpolation/doc_tex/Interpolation_ref/intro.tex
index 9de70b86e9d..cf6b5c97e65 100644
--- a/Packages/Interpolation/doc_tex/Interpolation_ref/intro.tex
+++ b/Packages/Interpolation/doc_tex/Interpolation_ref/intro.tex
@@ -73,7 +73,11 @@ For a more thorough introduction see the user manual.
 \ccHeading{Natural neighbor coordinate computation}
 \ccRefIdfierPage{CGAL::natural_neighbor_coordinates_2} \\
 \ccRefIdfierPage{CGAL::regular_neighbor_coordinates_2} \\
-\ccRefIdfierPage{CGAL::surface_neighbor_coordinates_2_3} \\
+
+\ccHeading{Surface neighbor coordinate computation}
+\ccRefIdfierPage{CGAL::Voronoi_intersection_2_traits_3<K>}\\
+\ccRefIdfierPage{CGAL::surface_neighbor_coordinates_3} \\
+
 \clearpage
 
 \lcHtml{\ccHeading{Alphabetical Listing of Reference Pages}}
diff --git a/Packages/Interpolation/doc_tex/Interpolation_ref/linear_interpolation.tex b/Packages/Interpolation/doc_tex/Interpolation_ref/linear_interpolation.tex
index 383a405afd9..3dd2828ebeb 100644
--- a/Packages/Interpolation/doc_tex/Interpolation_ref/linear_interpolation.tex
+++ b/Packages/Interpolation/doc_tex/Interpolation_ref/linear_interpolation.tex
@@ -65,7 +65,7 @@ values provided via a functor.
   \ccRefIdfierPage{CGAL::quadratic_interpolation} \\
   \ccRefIdfierPage{CGAL::natural_neighbor_coordinates_2}\\
   \ccRefIdfierPage{CGAL::regular_neighbor_coordinates_2} \\
-  \ccRefIdfierPage{CGAL::surface_neighbor_coordinates_2_3}\\
+  \ccRefIdfierPage{CGAL::surface_neighbor_coordinates_3}\\
                                 %\ccImplementation 
   
   \ccIndexSubitemEnd{Interpolation}{linear_interpolation}
diff --git a/Packages/Interpolation/doc_tex/Interpolation_ref/main.tex b/Packages/Interpolation/doc_tex/Interpolation_ref/main.tex
index 04434d04a3e..e17c8b8c09d 100644
--- a/Packages/Interpolation/doc_tex/Interpolation_ref/main.tex
+++ b/Packages/Interpolation/doc_tex/Interpolation_ref/main.tex
@@ -7,6 +7,12 @@
 \input{Interpolation_ref/InterpolationTraits}
 \input{Interpolation_ref/Interpolation_traits_2}
 \input{Interpolation_ref/natural_neighbor_coordinates_2}
+\input{Interpolation_ref/regular_neighbor_coordinates_2}
+\input{Interpolation_ref/Regular_neighbor_coordinates_traits_2}
 \input{Interpolation_ref/sibson_gradient_fitting}
 \input{Interpolation_ref/GradientFittingTraits}
-\input{Interpolation_ref/Interpolation_gradient_fitting_traits_2}
\ No newline at end of file
+\input{Interpolation_ref/Interpolation_gradient_fitting_traits_2}
+%
+%Surface coordinate stuff:
+\input{Interpolation_ref/surface_neighbor_coordinates_3}
+\input{Interpolation_ref/Voronoi_intersection_2_traits_3}
\ No newline at end of file
diff --git a/Packages/Interpolation/doc_tex/Interpolation_ref/natural_neighbor_coordinates_2.tex b/Packages/Interpolation/doc_tex/Interpolation_ref/natural_neighbor_coordinates_2.tex
index 3429354fd70..b213ab627b2 100644
--- a/Packages/Interpolation/doc_tex/Interpolation_ref/natural_neighbor_coordinates_2.tex
+++ b/Packages/Interpolation/doc_tex/Interpolation_ref/natural_neighbor_coordinates_2.tex
@@ -31,7 +31,7 @@ of the triangulation.
   computes the natural neighbor coordinates for \ccc{p} with respect
   to the points in the two-dimensional Delaunay triangulation \ccc{dt}.
   The template class \ccc{Dt} should be of type
-  \ccc{Delaunay_triangulation_2}.  The value type of the
+  \ccc{Delaunay_triangulation_2<Traits, Tds>}.  The value type of the
   \ccc{OutputIterator} is a pair of \ccc{Dt::Point_2} and the
   coordinate value of type \ccc{Dt::Geom_traits::FT}. The sequence of
   point-coordinate pairs that is computed by the function is placed
@@ -63,7 +63,7 @@ of the triangulation.
   over the boundary edges of the conflict zone of \ccc{p} in the
   triangulation. It is the result of the function
   \ccc{T.get_boundary_of_conflicts(p,std::back_inserter(hole),
-    start);}, see \ccc{Delaunay_triangulation_2}.}
+    start)}, see \ccc{Delaunay_triangulation_2<Traits, Tds>}.}
 
 \ccFunction{template <class Dt, class OutputIterator> std::pair<
   OutputIterator, typename Dt::Geom_traits::FT >
@@ -74,10 +74,17 @@ of the triangulation.
   \ccc{Dt} excluding \ccc{vh->point()}. The same as above for the
   remaining parameters.}
 
+\ccFunction{template <class Dt, class OutputIterator> 
+  std::pair< OutputIterator, typename Dt::Geom_traits::FT > 
+  regular_neighbor_coordinates_2(const Dt& rt,
+  typename Dt::Vertex_handle vh, OutputIterator out);}{The template
+  parameter \ccc{Traits} of the function above is instantiated to
+  \ccc{Dt::Geom_traits}.}
+
 \ccHeading{Requirements}
 \begin{enumerate}
 \item \ccc{Dt} are equivalent to the class
-  \ccc{Delaunay_triangulation_2}.
+  \ccc{Delaunay_triangulation_2<Traits, Tds>}.
 \item \ccc{OutputIterator::value_type} is equivalent to
   \ccc{std::pair<Dt::Point_2, Dt::Geom_traits::FT>}, i.e.\ a pair
   asscociating a point and its natural neighbor coordinate.
@@ -94,10 +101,10 @@ of the triangulation.
 \end{enumerate}
 
 \ccSeeAlso
-\ccRefIdfierPage{CGAL::linear_interpolation.h} \\
-\ccRefIdfierPage{CGAL::Sibson_c1_interpolation.h}\\ 
-\ccRefIdfierPage{CGAL::surface_coordinates_3.h} \\
-\ccRefIdfierPage{CGAL::regular_neighbor_coordinates_2.h} 
+\ccRefIdfierPage{CGAL::linear_interpolation} \\
+\ccRefIdfierPage{CGAL::sibson_c1_interpolation}\\ 
+\ccRefIdfierPage{CGAL::surface_neighbor_coordinates_3} \\
+\ccRefIdfierPage{CGAL::regular_neighbor_coordinates_2} 
 
 \ccImplementation This function computes the areas stolen from the
 Voronoi cells of points in \ccc{dt} by the insertion of \ccc{p}. The
diff --git a/Packages/Interpolation/doc_tex/Interpolation_ref/sibson_c1_interpolation.tex b/Packages/Interpolation/doc_tex/Interpolation_ref/sibson_c1_interpolation.tex
index 63ecfe41cd2..f335f72b6ac 100644
--- a/Packages/Interpolation/doc_tex/Interpolation_ref/sibson_c1_interpolation.tex
+++ b/Packages/Interpolation/doc_tex/Interpolation_ref/sibson_c1_interpolation.tex
@@ -70,11 +70,11 @@ function gradient given a point.
   the type of the gradient of the function value and a boolean. The
   function gradient type must provide a multiplication operation with
   \ccc{Traits::Vector_d}.
-\item A model of the functor types \ccc{Functor} and \ccc{GradFunctor}
-  resp.\ is provided by the struct \ccc{Data_access}. It must be
-  instantiated accordingly with an associative container (e.g.\ \stl\ 
-  \ccc{std::map}) having the point type as \ccc{key_type} and the
-  function value type and function gradient type resp.\ as
+\item A model of the functor types \ccc{Functor} (resp.\ 
+  \ccc{GradFunctor}) is provided by the struct \ccc{Data_access}. It
+  must be instantiated accordingly with an associative container
+  (e.g.\ \stl\ \ccc{std::map}) having the point type as \ccc{key_type}
+  and the function value type (resp.\ function gradient type) as
   \ccc{mapped_type}.
 \item The number type \ccc{FT} provided by \ccc{Traits} must support
   the square root operation \ccc{sqrt()}.
@@ -82,7 +82,7 @@ function gradient given a point.
 
 \ccFunction{template < class ForwardIterator, class Functor, class
   GradFunctor, class Traits> typename Functor::result_type
-  Sibson_c1_interpolation_square(ForwardIterator first,
+  sibson_c1_interpolation_square(ForwardIterator first,
   ForwardIterator beyond, const typename
   std::iterator_traits<ForwardIterator>::value_type::second_type&
   norm, Functor function_value, GradFunctor function_gradient,const
@@ -101,7 +101,7 @@ function gradient given a point.
 \ccRefIdfierPage{CGAL::Interpolation_gradient_fitting_traits_2<K>} \\
 \ccRefIdfierPage{CGAL::natural_neighbor_coordinates_2}\\
 \ccRefIdfierPage{CGAL::regular_neighbor_coordinates_2} \\
-\ccRefIdfierPage{CGAL::surface_neighbor_coordinates_2_3}\\
+\ccRefIdfierPage{CGAL::surface_neighbor_coordinates_3}\\
 %\ccImplementation
 
 \ccIndexSubitemEnd{Interpolation}{sibson_c1_interpolation}
diff --git a/Packages/Interpolation/doc_tex/Interpolation_ref/sibson_gradient_fitting.tex b/Packages/Interpolation/doc_tex/Interpolation_ref/sibson_gradient_fitting.tex
index 6e804778bf9..ae408b636a0 100644
--- a/Packages/Interpolation/doc_tex/Interpolation_ref/sibson_gradient_fitting.tex
+++ b/Packages/Interpolation/doc_tex/Interpolation_ref/sibson_gradient_fitting.tex
@@ -80,10 +80,10 @@ coordinates.
 \ccRefIdfierPage{CGAL::sibson_c1_interpolation} \\
 \ccRefIdfierPage{CGAL::farin_c1_interpolation} \\
 \ccRefIdfierPage{CGAL::quadratic_interpolation} \\
-\ccRefIdfierPage{CGAL::Interpolation_gradient_fitting_traits_2} \\
+\ccRefIdfierPage{CGAL::Interpolation_gradient_fitting_traits_2<K>} \\
 \ccRefIdfierPage{CGAL::natural_neighbor_coordinates_2}\\
 \ccRefIdfierPage{CGAL::regular_neighbor_coordinates_2} \\
-\ccRefIdfierPage{CGAL::surface_neighbor_coordinates_2_3}\\
+\ccRefIdfierPage{CGAL::surface_neighbor_coordinates_3}\\
 
 \ccImplementation This function implements Sibson's gradient
 estimation method based on natural neighbor coordinates
diff --git a/Packages/Interpolation/doc_tex/basic/Interpolation_ref/Data_access.tex b/Packages/Interpolation/doc_tex/basic/Interpolation_ref/Data_access.tex
new file mode 100644
index 00000000000..2bf9d2d0029
--- /dev/null
+++ b/Packages/Interpolation/doc_tex/basic/Interpolation_ref/Data_access.tex
@@ -0,0 +1,64 @@
+% +------------------------------------------------------------------------+
+% | Reference manual page: Data_access.tex
+% +------------------------------------------------------------------------+
+% | 
+% | August 2003  Julia Floetotto
+% | Package:   Interpolation
+% | 
+% |
+% +------------------------------------------------------------------------+
+\begin{ccRefClass}{Data_access<Map>}  %% add template arg's if necessary
+\ccIndexMainItemBegin{Data_access<Map>}
+
+\ccDefinition
+  
+The struct \ccRefName\ implements a functor that allows to retrieve
+data from a associative container. The functor keeps a reference to
+the container. Given an instance of the container's key type, it
+returns a pair of the container's value type and a boolean indicating
+whether the retrieval was successful.\\
+
+This class can be used to provide the function\_values and
+function\_gradients of the interpolation functions.
+
+\ccInclude{CGAL/interpolation_functions.h} 
+
+\ccParameters The class
+\ccRefName\ has the container type \ccc{Map} as template parameter.
+
+\ccTypes 
+\ccTypedef{typedef Map::mapped_type Data_type;}{}
+\ccGlue
+\ccTypedef{typedef Map::key_type Key_type;}{}
+\ccCreation
+\ccCreationVariable{data_access}  %% choose variable name
+\ccConstructor{Data_acces(const Map& map);}
+{Introduces a \ccc{Data_access} to the container \ccc{map}.}
+
+
+\ccMethod{std::pair< Data_type, bool> operator()(const Key_type& p);} {If
+  there is an entry for \ccc{p} in the container \ccc{map}, then the
+  pair of \ccc{map.find(p)} and \ccc{true} is returned. Otherwise, the
+  boolean value of the pair is \ccc{false}.}
+
+\ccSeeAlso
+\ccRefIdfierPage{CGAL::linear_interpolation}\\
+\ccRefIdfierPage{CGAL::sibson_c1_interpolation} \\
+\ccRefIdfierPage{CGAL::sibson_gradient_fitting} \\
+\ccRefIdfierPage{CGAL::farin_c1_interpolation} \\
+\ccRefIdfierPage{CGAL::quadratic_interpolation} \\
+%\ccRefIdfierPage{CGAL::natural_neighbor_coordinates_2}\\
+%\ccRefIdfierPage{CGAL::regular_neighbor_coordinates_2} \\
+%\ccRefIdfierPage{CGAL::surface_neighbor_coordinates_2_3}\\
+%\ccImplementation 
+
+\ccIndexSubitemEnd{Interpolation}{Data_access}
+\end{ccRefClass}
+%\renewcommand{\ccRefPageBegin}{}
+%\renewcommand{\ccRefPageEnd}{}
+
+% +------------------------------------------------------------------------+
+%%RefPage: end of main body, begin of footer
+% EOF
+% +------------------------------------------------------------------------+
+
diff --git a/Packages/Interpolation/doc_tex/basic/Interpolation_ref/InterpolationTraits.tex b/Packages/Interpolation/doc_tex/basic/Interpolation_ref/InterpolationTraits.tex
index bb45b83d100..72475e48ec5 100644
--- a/Packages/Interpolation/doc_tex/basic/Interpolation_ref/InterpolationTraits.tex
+++ b/Packages/Interpolation/doc_tex/basic/Interpolation_ref/InterpolationTraits.tex
@@ -1,5 +1,5 @@
 % +------------------------------------------------------------------------+
-% | Reference manual page: Interpolation_traits.tex
+% | Reference manual page: InterpolationTraits.tex
 % +------------------------------------------------------------------------+
 % | 06.04.2000   Author
 % | Package: Package
@@ -24,7 +24,7 @@ defines the primitives used in the algorithms.  The concept
   \ccc{FieldNumberType}.}  \ccNestedType{Point_d}{The point type on
   which the function is defined and interpolated.}
 \ccNestedType{Vector_d}{The corresponding vector type.}
-\ccNestedType{Construct_vector} {A constructor object for
+\ccNestedType{Construct_vector_d} {A constructor object for
   \ccc{Vector_d}.
   Provides~: \\
   \ccc{Vector_d operator() (Point_d a, Point_d b)} which produces the
@@ -32,40 +32,28 @@ defines the primitives used in the algorithms.  The concept
     a} and \\
   \ccc{Vector_d operator() (Null_vector NULL_VECTOR)} which introduces
   the null vector.}
-\ccNestedType{Construct_scaled_vector}{Constructor object for
+\ccNestedType{Construct_scaled_vector_d}{Constructor object for
   \ccc{Vector_d}.
   Provides~: \\
   \ccc{Vector_d operator() (Vector_d v,FT scale)} which produces the
   vector \ccc{v} scaled by a factor \ccc{scale}.}
-\ccNestedType{Compute_squared_distance} {Constructor
+\ccNestedType{Compute_squared_distance_d} {Constructor
   object for \ccc{FT}. Provides the operator:\\
   \ccc{FT operator() (Point_d a, Point_d b)} returning the squared
   distance between \ccc{a} and \ccc{b}.}
 
 %\ccCreation
 \ccCreationVariable{traits}  %% choose variable name
-%A default constructor and a copy constructor
-% and an assignement operator are required. 
-%Note that further constructors
-%can be provided. 
-%%\ccThree{InterpolationTraits}{traits=gtrxx  }{}
 \ccConstructor{InterpolationTraits();}{default constructor.}
-%\ccGlue
-%\ccConstructor{InterpolationTraits(InterpolationTraits gtr);}
-%{Copy constructor}
-%\ccMethod{InterpolationTraits operator=(InterpolationTraits gtr);}
-%{Assignment operator.}
 
-%\ccOperations 
 
 \ccHeading{Construction objects} The following functions that
 create instances of the above constructor object types must exist.
-
-\ccMethod{Construct_vector construct_vector_object();}{}
+\ccMethod{Construct_vector_d construct_vector_d_object();}{}
 \ccGlue
-\ccMethod{Construct_scaled_vector construct_scaled_vector_object();}{}
+\ccMethod{Construct_scaled_vector_d construct_scaled_vector_d_object();}{}
 \ccGlue
-\ccMethod{Compute_squared_distance compute_squared_distance_object();}{}
+\ccMethod{Compute_squared_distance_d compute_squared_distance_d_object();}{}
 
 \ccHasModels
 \ccc{CGAL::Interpolation_traits_2<K>} \\
@@ -78,7 +66,7 @@ create instances of the above constructor object types must exist.
 \ccRefIdfierPage{CGAL::quadratic_interpolation} \\
 \ccRefIdfierPage{CGAL::natural_neighbor_coordinates_2}\\
 \ccRefIdfierPage{CGAL::regular_neighbor_coordinates_2} \\
-\ccRefIdfierPage{CGAL::surface_neighbor_coordinates_2_3}\\
+\ccRefIdfierPage{CGAL::surface_neighbor_coordinates_3}\\
 \end{ccRefConcept}
 
 % +------------------------------------------------------------------------+
diff --git a/Packages/Interpolation/doc_tex/basic/Interpolation_ref/Interpolation_gradient_fitting_traits_2.tex b/Packages/Interpolation/doc_tex/basic/Interpolation_ref/Interpolation_gradient_fitting_traits_2.tex
index 6a7a7099c90..aeab07aabe9 100644
--- a/Packages/Interpolation/doc_tex/basic/Interpolation_ref/Interpolation_gradient_fitting_traits_2.tex
+++ b/Packages/Interpolation/doc_tex/basic/Interpolation_ref/Interpolation_gradient_fitting_traits_2.tex
@@ -32,22 +32,22 @@
   \ccIndexSubitem[c]{InterpolationTraits}{model} \\
   \ccIndexSubitem[c]{GradientFittingTraits}{model} \\
    \ccTypes
-  \ccSetThreeColumns{typedef R::Less_signed_distance_to_line_2 }{Less_signed_distance_to_line_2}{}
+  \ccSetThreeColumns{typedef K::Less_signed_distance_to_line_2 }{Less_signed_distance_to_line_2}{}
   \ccThreeToTwo 
-  \ccTypedef{typedef R::FT                            FT;}{}
+  \ccTypedef{typedef K::FT                            FT;}{}
   \ccGlue
-  \ccTypedef{typedef R::Point_2                        Point_d;}{}
+  \ccTypedef{typedef K::Point_2                        Point_d;}{}
   \ccGlue
-  \ccTypedef{typedef R::Vector_2                       Vector_d;}{}
+  \ccTypedef{typedef K::Vector_2                       Vector_d;}{}
   \ccGlue
-\ccTypedef{typedef R::Aff_transformation_2         Aff_transformation_d;}{}
+\ccTypedef{typedef K::Aff_transformation_2         Aff_transformation_d;}{}
 \ccGlue
-  \ccTypedef{typedef R::Construct_vector_2     Construct_vector_d;}{}
+  \ccTypedef{typedef K::Construct_vector_2     Construct_vector_d;}{}
   \ccGlue
-  \ccTypedef{typedef R::Construct_scaled_vector_2 
+  \ccTypedef{typedef K::Construct_scaled_vector_2 
     Construct_scaled_vector_d;}{}
   \ccGlue
-  \ccTypedef{typedef R::Compute_squared_distance_2 
+  \ccTypedef{typedef K::Compute_squared_distance_2 
     Compute_squared_distance_d;}{}
   \ccGlue
   \ccTypedef{typedef Construct_null_matrix_2<Aff_transformation_d>   
@@ -58,7 +58,7 @@
   \ccGlue
   \ccTypedef{typedef Construct_sum_matrix_2<Aff_transformation_d> Construct_sum_matrix_d;}{}
   \ccGlue
-  \ccTypedef{typedef Construct_outer_product_2<Rep>             Construct_outer_product_d;}{}
+  \ccTypedef{typedef Construct_outer_product_2<K>             Construct_outer_product_d;}{}
                       
 \ccOperations
 \ccAutoIndexingOff
diff --git a/Packages/Interpolation/doc_tex/basic/Interpolation_ref/Interpolation_traits_2.tex b/Packages/Interpolation/doc_tex/basic/Interpolation_ref/Interpolation_traits_2.tex
index 6058eaeb091..622c0327f3c 100644
--- a/Packages/Interpolation/doc_tex/basic/Interpolation_ref/Interpolation_traits_2.tex
+++ b/Packages/Interpolation/doc_tex/basic/Interpolation_ref/Interpolation_traits_2.tex
@@ -11,7 +11,7 @@
 % +------------------------------------------------------------------------+
 
 
-\begin{ccRefClass}{Interpolation_traits_2<R>}
+\begin{ccRefClass}{Interpolation_traits_2<K>}
   
 
 
@@ -24,7 +24,7 @@
   \ccc{InterpolationTraits} and can be used to instantiate the
   geometric traits class of interpolation methods applied on a
   bivariate function over a two-dimensional domain. The traits class
-  is templated by a kernel class \ccc{R}.
+  is templated by a kernel class \ccc{K}.
   
   \ccInclude{CGAL/Interpolation_traits_2.h}
 
@@ -33,20 +33,19 @@
   \ccIndexSubitem[c]{InterpolationTraits}{model} \\
   
   \ccTypes
-  \ccSetThreeColumns{typedef R::Less_signed_distance_to_line_2 }{Less_signed_distance_to_line_2}{}
   \ccThreeToTwo 
-  \ccTypedef{typedef R::FT                            FT;}{}
+  \ccTypedef{typedef K::FT                            FT;}{}
   \ccGlue
-  \ccTypedef{typedef R::Point_2                        Point_d;}{}
+  \ccTypedef{typedef K::Point_2                        Point_d;}{}
   \ccGlue
-  \ccTypedef{typedef R::Vector_2                       Vector_d;}{}
+  \ccTypedef{typedef K::Vector_2                       Vector_d;}{}
   \ccGlue
-  \ccTypedef{typedef R::Construct_vector_2     Construct_vector_d;}{}
+  \ccTypedef{typedef K::Construct_vector_2     Construct_vector_d;}{}
   \ccGlue
-  \ccTypedef{typedef R::Construct_scaled_vector_2 
+  \ccTypedef{typedef K::Construct_scaled_vector_2 
     Construct_scaled_vector_d;}{}
   \ccGlue
-  \ccTypedef{typedef R::Compute_squared_distance_2 
+  \ccTypedef{typedef K::Compute_squared_distance_2 
     Compute_squared_distance_d;}{}
                                            
 \ccOperations
diff --git a/Packages/Interpolation/doc_tex/basic/Interpolation_ref/farin_c1_interpolation.tex b/Packages/Interpolation/doc_tex/basic/Interpolation_ref/farin_c1_interpolation.tex
index bbcb967ed9f..15dde9d6e64 100644
--- a/Packages/Interpolation/doc_tex/basic/Interpolation_ref/farin_c1_interpolation.tex
+++ b/Packages/Interpolation/doc_tex/basic/Interpolation_ref/farin_c1_interpolation.tex
@@ -55,7 +55,7 @@ to provide the square root operation.
 \ccRefIdfierPage{CGAL::Interpolation_traits_2} \\
 \ccRefIdfierPage{CGAL::natural_neighbor_coordinates_2}\\
 \ccRefIdfierPage{CGAL::regular_neighbor_coordinates_2} \\
-\ccRefIdfierPage{CGAL::surface_neighbor_coordinates_2_3}\\
+\ccRefIdfierPage{CGAL::surface_neighbor_coordinates_3}\\
 %\ccImplementation
 
 \ccIndexSubitemEnd{Interpolation}{farin_c1_interpolation}
diff --git a/Packages/Interpolation/doc_tex/basic/Interpolation_ref/intro.tex b/Packages/Interpolation/doc_tex/basic/Interpolation_ref/intro.tex
index 9de70b86e9d..cf6b5c97e65 100644
--- a/Packages/Interpolation/doc_tex/basic/Interpolation_ref/intro.tex
+++ b/Packages/Interpolation/doc_tex/basic/Interpolation_ref/intro.tex
@@ -73,7 +73,11 @@ For a more thorough introduction see the user manual.
 \ccHeading{Natural neighbor coordinate computation}
 \ccRefIdfierPage{CGAL::natural_neighbor_coordinates_2} \\
 \ccRefIdfierPage{CGAL::regular_neighbor_coordinates_2} \\
-\ccRefIdfierPage{CGAL::surface_neighbor_coordinates_2_3} \\
+
+\ccHeading{Surface neighbor coordinate computation}
+\ccRefIdfierPage{CGAL::Voronoi_intersection_2_traits_3<K>}\\
+\ccRefIdfierPage{CGAL::surface_neighbor_coordinates_3} \\
+
 \clearpage
 
 \lcHtml{\ccHeading{Alphabetical Listing of Reference Pages}}
diff --git a/Packages/Interpolation/doc_tex/basic/Interpolation_ref/linear_interpolation.tex b/Packages/Interpolation/doc_tex/basic/Interpolation_ref/linear_interpolation.tex
index 383a405afd9..3dd2828ebeb 100644
--- a/Packages/Interpolation/doc_tex/basic/Interpolation_ref/linear_interpolation.tex
+++ b/Packages/Interpolation/doc_tex/basic/Interpolation_ref/linear_interpolation.tex
@@ -65,7 +65,7 @@ values provided via a functor.
   \ccRefIdfierPage{CGAL::quadratic_interpolation} \\
   \ccRefIdfierPage{CGAL::natural_neighbor_coordinates_2}\\
   \ccRefIdfierPage{CGAL::regular_neighbor_coordinates_2} \\
-  \ccRefIdfierPage{CGAL::surface_neighbor_coordinates_2_3}\\
+  \ccRefIdfierPage{CGAL::surface_neighbor_coordinates_3}\\
                                 %\ccImplementation 
   
   \ccIndexSubitemEnd{Interpolation}{linear_interpolation}
diff --git a/Packages/Interpolation/doc_tex/basic/Interpolation_ref/main.tex b/Packages/Interpolation/doc_tex/basic/Interpolation_ref/main.tex
index 04434d04a3e..e17c8b8c09d 100644
--- a/Packages/Interpolation/doc_tex/basic/Interpolation_ref/main.tex
+++ b/Packages/Interpolation/doc_tex/basic/Interpolation_ref/main.tex
@@ -7,6 +7,12 @@
 \input{Interpolation_ref/InterpolationTraits}
 \input{Interpolation_ref/Interpolation_traits_2}
 \input{Interpolation_ref/natural_neighbor_coordinates_2}
+\input{Interpolation_ref/regular_neighbor_coordinates_2}
+\input{Interpolation_ref/Regular_neighbor_coordinates_traits_2}
 \input{Interpolation_ref/sibson_gradient_fitting}
 \input{Interpolation_ref/GradientFittingTraits}
-\input{Interpolation_ref/Interpolation_gradient_fitting_traits_2}
\ No newline at end of file
+\input{Interpolation_ref/Interpolation_gradient_fitting_traits_2}
+%
+%Surface coordinate stuff:
+\input{Interpolation_ref/surface_neighbor_coordinates_3}
+\input{Interpolation_ref/Voronoi_intersection_2_traits_3}
\ No newline at end of file
diff --git a/Packages/Interpolation/doc_tex/basic/Interpolation_ref/natural_neighbor_coordinates_2.tex b/Packages/Interpolation/doc_tex/basic/Interpolation_ref/natural_neighbor_coordinates_2.tex
index 3429354fd70..b213ab627b2 100644
--- a/Packages/Interpolation/doc_tex/basic/Interpolation_ref/natural_neighbor_coordinates_2.tex
+++ b/Packages/Interpolation/doc_tex/basic/Interpolation_ref/natural_neighbor_coordinates_2.tex
@@ -31,7 +31,7 @@ of the triangulation.
   computes the natural neighbor coordinates for \ccc{p} with respect
   to the points in the two-dimensional Delaunay triangulation \ccc{dt}.
   The template class \ccc{Dt} should be of type
-  \ccc{Delaunay_triangulation_2}.  The value type of the
+  \ccc{Delaunay_triangulation_2<Traits, Tds>}.  The value type of the
   \ccc{OutputIterator} is a pair of \ccc{Dt::Point_2} and the
   coordinate value of type \ccc{Dt::Geom_traits::FT}. The sequence of
   point-coordinate pairs that is computed by the function is placed
@@ -63,7 +63,7 @@ of the triangulation.
   over the boundary edges of the conflict zone of \ccc{p} in the
   triangulation. It is the result of the function
   \ccc{T.get_boundary_of_conflicts(p,std::back_inserter(hole),
-    start);}, see \ccc{Delaunay_triangulation_2}.}
+    start)}, see \ccc{Delaunay_triangulation_2<Traits, Tds>}.}
 
 \ccFunction{template <class Dt, class OutputIterator> std::pair<
   OutputIterator, typename Dt::Geom_traits::FT >
@@ -74,10 +74,17 @@ of the triangulation.
   \ccc{Dt} excluding \ccc{vh->point()}. The same as above for the
   remaining parameters.}
 
+\ccFunction{template <class Dt, class OutputIterator> 
+  std::pair< OutputIterator, typename Dt::Geom_traits::FT > 
+  regular_neighbor_coordinates_2(const Dt& rt,
+  typename Dt::Vertex_handle vh, OutputIterator out);}{The template
+  parameter \ccc{Traits} of the function above is instantiated to
+  \ccc{Dt::Geom_traits}.}
+
 \ccHeading{Requirements}
 \begin{enumerate}
 \item \ccc{Dt} are equivalent to the class
-  \ccc{Delaunay_triangulation_2}.
+  \ccc{Delaunay_triangulation_2<Traits, Tds>}.
 \item \ccc{OutputIterator::value_type} is equivalent to
   \ccc{std::pair<Dt::Point_2, Dt::Geom_traits::FT>}, i.e.\ a pair
   asscociating a point and its natural neighbor coordinate.
@@ -94,10 +101,10 @@ of the triangulation.
 \end{enumerate}
 
 \ccSeeAlso
-\ccRefIdfierPage{CGAL::linear_interpolation.h} \\
-\ccRefIdfierPage{CGAL::Sibson_c1_interpolation.h}\\ 
-\ccRefIdfierPage{CGAL::surface_coordinates_3.h} \\
-\ccRefIdfierPage{CGAL::regular_neighbor_coordinates_2.h} 
+\ccRefIdfierPage{CGAL::linear_interpolation} \\
+\ccRefIdfierPage{CGAL::sibson_c1_interpolation}\\ 
+\ccRefIdfierPage{CGAL::surface_neighbor_coordinates_3} \\
+\ccRefIdfierPage{CGAL::regular_neighbor_coordinates_2} 
 
 \ccImplementation This function computes the areas stolen from the
 Voronoi cells of points in \ccc{dt} by the insertion of \ccc{p}. The
diff --git a/Packages/Interpolation/doc_tex/basic/Interpolation_ref/sibson_c1_interpolation.tex b/Packages/Interpolation/doc_tex/basic/Interpolation_ref/sibson_c1_interpolation.tex
index 63ecfe41cd2..f335f72b6ac 100644
--- a/Packages/Interpolation/doc_tex/basic/Interpolation_ref/sibson_c1_interpolation.tex
+++ b/Packages/Interpolation/doc_tex/basic/Interpolation_ref/sibson_c1_interpolation.tex
@@ -70,11 +70,11 @@ function gradient given a point.
   the type of the gradient of the function value and a boolean. The
   function gradient type must provide a multiplication operation with
   \ccc{Traits::Vector_d}.
-\item A model of the functor types \ccc{Functor} and \ccc{GradFunctor}
-  resp.\ is provided by the struct \ccc{Data_access}. It must be
-  instantiated accordingly with an associative container (e.g.\ \stl\ 
-  \ccc{std::map}) having the point type as \ccc{key_type} and the
-  function value type and function gradient type resp.\ as
+\item A model of the functor types \ccc{Functor} (resp.\ 
+  \ccc{GradFunctor}) is provided by the struct \ccc{Data_access}. It
+  must be instantiated accordingly with an associative container
+  (e.g.\ \stl\ \ccc{std::map}) having the point type as \ccc{key_type}
+  and the function value type (resp.\ function gradient type) as
   \ccc{mapped_type}.
 \item The number type \ccc{FT} provided by \ccc{Traits} must support
   the square root operation \ccc{sqrt()}.
@@ -82,7 +82,7 @@ function gradient given a point.
 
 \ccFunction{template < class ForwardIterator, class Functor, class
   GradFunctor, class Traits> typename Functor::result_type
-  Sibson_c1_interpolation_square(ForwardIterator first,
+  sibson_c1_interpolation_square(ForwardIterator first,
   ForwardIterator beyond, const typename
   std::iterator_traits<ForwardIterator>::value_type::second_type&
   norm, Functor function_value, GradFunctor function_gradient,const
@@ -101,7 +101,7 @@ function gradient given a point.
 \ccRefIdfierPage{CGAL::Interpolation_gradient_fitting_traits_2<K>} \\
 \ccRefIdfierPage{CGAL::natural_neighbor_coordinates_2}\\
 \ccRefIdfierPage{CGAL::regular_neighbor_coordinates_2} \\
-\ccRefIdfierPage{CGAL::surface_neighbor_coordinates_2_3}\\
+\ccRefIdfierPage{CGAL::surface_neighbor_coordinates_3}\\
 %\ccImplementation
 
 \ccIndexSubitemEnd{Interpolation}{sibson_c1_interpolation}
diff --git a/Packages/Interpolation/doc_tex/basic/Interpolation_ref/sibson_gradient_fitting.tex b/Packages/Interpolation/doc_tex/basic/Interpolation_ref/sibson_gradient_fitting.tex
index 6e804778bf9..ae408b636a0 100644
--- a/Packages/Interpolation/doc_tex/basic/Interpolation_ref/sibson_gradient_fitting.tex
+++ b/Packages/Interpolation/doc_tex/basic/Interpolation_ref/sibson_gradient_fitting.tex
@@ -80,10 +80,10 @@ coordinates.
 \ccRefIdfierPage{CGAL::sibson_c1_interpolation} \\
 \ccRefIdfierPage{CGAL::farin_c1_interpolation} \\
 \ccRefIdfierPage{CGAL::quadratic_interpolation} \\
-\ccRefIdfierPage{CGAL::Interpolation_gradient_fitting_traits_2} \\
+\ccRefIdfierPage{CGAL::Interpolation_gradient_fitting_traits_2<K>} \\
 \ccRefIdfierPage{CGAL::natural_neighbor_coordinates_2}\\
 \ccRefIdfierPage{CGAL::regular_neighbor_coordinates_2} \\
-\ccRefIdfierPage{CGAL::surface_neighbor_coordinates_2_3}\\
+\ccRefIdfierPage{CGAL::surface_neighbor_coordinates_3}\\
 
 \ccImplementation This function implements Sibson's gradient
 estimation method based on natural neighbor coordinates