From 1a0fb4f0bd71489d479a4b88d0849a684f9e438d Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Julia=20Fl=C3=B6totto?= <Julia.Flototto@sophia.inria.fr>
Date: Fri, 28 Nov 2003 10:37:33 +0000
Subject: [PATCH] - first addition of example and demo files for interpolation
 and 2D coordinate computation.

---
 .../Interpolation/demo/Interpolation/README   |  30 ++
 .../demo/Interpolation/demo_interpolation_2.C | 348 ++++++++++++++++++
 .../demo/Interpolation/demo_numerical_2.C     | 239 ++++++++++++
 .../Interpolation/demo/Interpolation/makefile |  56 +++
 .../demo/Interpolation/makefile_local         |  59 +++
 .../examples/Interpolation/README             |  17 +
 .../Interpolation/linear_interpolation_2.C    |  67 ++++
 .../examples/Interpolation/nn_coordinates_2.C |  29 ++
 .../examples/Interpolation/rn_coordinates_2.C |  34 ++
 .../Interpolation/sibson_interpolation_2.C    |  81 ++++
 10 files changed, 960 insertions(+)
 create mode 100644 Packages/Interpolation/demo/Interpolation/README
 create mode 100644 Packages/Interpolation/demo/Interpolation/demo_interpolation_2.C
 create mode 100644 Packages/Interpolation/demo/Interpolation/demo_numerical_2.C
 create mode 100644 Packages/Interpolation/demo/Interpolation/makefile
 create mode 100644 Packages/Interpolation/demo/Interpolation/makefile_local
 create mode 100644 Packages/Interpolation/examples/Interpolation/README
 create mode 100644 Packages/Interpolation/examples/Interpolation/linear_interpolation_2.C
 create mode 100644 Packages/Interpolation/examples/Interpolation/nn_coordinates_2.C
 create mode 100644 Packages/Interpolation/examples/Interpolation/rn_coordinates_2.C
 create mode 100644 Packages/Interpolation/examples/Interpolation/sibson_interpolation_2.C

diff --git a/Packages/Interpolation/demo/Interpolation/README b/Packages/Interpolation/demo/Interpolation/README
new file mode 100644
index 00000000000..c716ba46741
--- /dev/null
+++ b/Packages/Interpolation/demo/Interpolation/README
@@ -0,0 +1,30 @@
+Some demos use Geomview 
+[see the chapter Geomview in the cgal manual - support library:
+Geomview 1.8.1 is required.  The geomview command must be in the user's $PATH,
+otherwise the program will not be able to execute.]
+ 
+------- demo_interpolation_2 -------------------------------------------------
+using Geomview
+
+This demo program dumps the plot of the different interpolation functions
+using a very simple data set
+shows:	the differentiablity (or not) at the data points
+	the closeness to the gradient at the data points	
+
+1) Construction of a Delaunay triangulation of the data points.
+ 
+2) Interpolation on a grid using a user chosen interpolation function.
+
+3) Dumps the plot of the interpolation functions (.off file)
+
+4) Displays the data points in geomview.
+--------------------------------------------------------------
+
+------- demo_numerical_2 -----------------------------------------------
+This demo program computes the interpolation of a quadratic function 
+known on a small set of random sample points. The gradients of the function are given.
+
+- verifies the barycentric property
+- computes the maximum and mean error of the different interpolation function 
+  compared to the exact function value.
+--------------------------------------------------------------
diff --git a/Packages/Interpolation/demo/Interpolation/demo_interpolation_2.C b/Packages/Interpolation/demo/Interpolation/demo_interpolation_2.C
new file mode 100644
index 00000000000..db82aa47e75
--- /dev/null
+++ b/Packages/Interpolation/demo/Interpolation/demo_interpolation_2.C
@@ -0,0 +1,348 @@
+// ============================================================================
+//
+// Copyright (c) 1998-1999 The CGAL Consortium
+//
+// This software and related documentation is part of an INTERNAL release
+// of the Computational Geometry Algorithms Library (CGAL). It is not
+// intended for general use.
+//
+// ----------------------------------------------------------------------------
+//
+// release       :
+// release_date  :
+//
+// file          : demo_2/
+// revision      : $Revision$
+// author(s)     : Julia Floetotto
+// coordinator   : INRIA Sophia Antipolis (Mariette Yvinec)
+//
+// ============================================================================
+// Geomview doesn't work on M$ at the moment, so we don't compile this
+// file.
+//**********************
+//demo 2D Interpolation over the plane - using 2D natural neighbors
+//**********************
+
+
+#if defined(__BORLANDC__) || defined(_MSC_VER)
+#include <iostream>
+int main()
+{
+  std::cerr << "Geomview doesn't work on Windows, so this demo doesn't work"
+            << std::endl;
+  return 0;
+}
+#else
+
+#include <CGAL/basic.h>
+#include <fstream>
+#include <iostream>
+#include <utility>
+
+
+#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
+
+//#include <CGAL/double.h>
+//#include <CGAL/Random.h>
+
+#include <CGAL/Delaunay_triangulation_2.h>
+#include <CGAL/natural_neighbor_coordinates_2.h>
+#include <CGAL/Interpolation_traits_2.h>
+#include <CGAL/interpolation_functions.h>
+
+#include <CGAL/point_generators_2.h> 
+#include <CGAL/copy_n.h>
+
+#include <CGAL/IO/Geomview_stream.h>
+
+
+struct K : CGAL::Exact_predicates_inexact_constructions_kernel {};
+
+typedef CGAL::Delaunay_triangulation_2<K>   Dt;
+typedef CGAL::Interpolation_traits_2<K>     ITraits;
+
+typedef K::FT                               Coord_type;
+typedef K::Point_2                          Point_2;
+typedef K::Vector_2                         Vector_2;
+
+typedef K::Point_3                          Point_3;
+typedef K::Vector_3                         Vector_3;
+typedef K::Segment_3                        Segment_3;
+
+typedef std::map<Point_2, Coord_type, K::Less_xy_2>     Point_value_map ;
+typedef std::map<Point_2, Vector_2, K::Less_xy_2 >      Point_vector_map;
+
+typedef std::vector<Point_3>                   Point_vector_3;
+typedef std::vector<Point_2>                     Point_vector_2;
+
+////////////////////// 
+// VISU GEOMVIEW
+////////////////////// 
+
+template<class Point_vector>
+void visu_points(CGAL::Geomview_stream & os, const Point_vector & points)
+{
+  int n = points.size();
+  for(int i=0; i<n ; i++)
+    os << points[i];
+}
+
+//////////////////////////////// 
+template< class Map >  
+struct DataAccess : public std::unary_function< typename Map::key_type,
+		    typename Map::mapped_type> {
+  typedef typename Map::mapped_type Data_type;
+  typedef typename Map::key_type  Point;
+  
+  //CONSTRUCTOR:
+  DataAccess< Map >(const Map& m): map(m){};
+  
+  //Functor
+  Data_type operator()(const Point& p) { 
+    
+    typename Map::const_iterator mit = map.find(p);
+    if(mit!= map.end())
+      return mit->second;
+    return Data_type();
+  };
+  
+  const Map& map;
+};
+
+//////////////////////////////////////////////////////////////////////////
+////POINT GENERATION:
+void generate_grid_points(Point_vector_2& points, int m, float h) 
+{
+  
+  //int n = (m+1)*(m+1);
+  int n = m*m;
+  points.clear();
+  points.reserve(n);
+
+  std::cout <<" generate " <<n   << " grid points in square of height " 
+	    << h <<std::endl;
+ 
+   
+  // Create n points from a 16 x 16 grid. Note that the double
+  // arithmetic _is_ sufficient to produce exact integer grid points.
+  // The distance between neighbors is 34 pixel = 510 / 15.
+  CGAL::points_on_square_grid_2((double) h, n,
+				std::back_inserter(points),
+				CGAL::Creator_uniform_2<double,Point_2>());
+  
+}
+template < class Point_vector>
+bool
+dump_off_file_quadrilateral(const Point_vector& points, int m) 
+{
+  //open file stream 
+  char fname[11];
+  CGAL_CLIB_STD::strcpy(fname, "result.off");
+
+  //dump reconstruction file:
+  std::ofstream os(fname,std::ios::out);
+  if (! os) {
+    std::cout <<"cannot open output file: result.off. "<< std::endl;
+    return false;
+  }
+
+
+
+  int n= points.size();
+  std::cout << n << " points " << (m-1)*(m-1) <<std::endl;
+  
+  //first line: OFF
+  // 2. line : number of vertices and quad_indices:
+  os << "OFF" << std::endl;
+  os << n << " " <<(m-1)*(m-1) << " 0" <<
+    std::endl;
+  os << std::endl;
+  
+  for(int i=0; i<n ; i++)
+    os << points[i]<< std::endl;
+  
+  //indices
+  for(int i=0; i< m-1;i++)
+    for(int j=0; j< m-1; j++)
+      os << "4  " << i*m + j  << "  "
+	 <<  j+1 + i*m 
+	 << "  " <<(j+1) + (i+1)*m << "  " <<j + (i+1)*m 
+	 << std::endl;
+  return true;  
+}
+
+
+////////////////////// 
+
+int main(int argc,  char* argv[])
+{
+
+  //parameters:
+  int m = 78;
+  double h = 10;
+  double g = 2;
+  Coord_type w(4);
+  Point_vector_2 sample;
+  //3D function graph: 2D points + function value in z-direction:
+  Point_vector_3 sample_3;
+
+  Dt T;
+  Point_value_map values;
+  Point_vector_map gradients;
+  
+  sample.push_back( Point_2(h,h));
+  sample.push_back( Point_2( -h,h));
+  sample.push_back( Point_2( h,-h ));
+  sample.push_back( Point_2( -h,-h));
+  sample.push_back( Point_2( 0.0, h/2));
+  sample.push_back( Point_2( 0.0,-h/2));
+  sample.push_back( Point_2( h/2, 0.0 ));
+  sample.push_back( Point_2( -h/2, 0.0));
+  sample.push_back( Point_2(0.0,0.0));
+  int s = sample.size();
+  for(int j=0; j<s ; j++){
+    T.insert(sample[j]);
+    
+    values.insert(std::make_pair(sample[j], Coord_type(0)));
+    gradients.insert(std::make_pair(sample[j], CGAL::NULL_VECTOR));
+    sample_3.push_back( Point_3(sample[j].x(),sample[j].y(),0));
+  }
+
+  Point_2 p = Point_2(h/3,h/3);
+  T.insert(p);
+  values.insert(std::make_pair(p, w));
+  gradients.insert(std::make_pair(p, Vector_2(-g,-g) ));
+   
+  p = Point_2(-h/3,-h/3);
+  T.insert(p);
+  values.insert(std::make_pair(p,w));
+  gradients.insert(std::make_pair(p, Vector_2(g,g) ));
+   
+  p = Point_2(-h/3,h/3);
+  T.insert(p);
+  values.insert(std::make_pair(p, w));
+  gradients.insert(std::make_pair(p, Vector_2(g,-g) ));
+   
+  p = Point_2(h/3,-h/3);
+  T.insert(p);
+  values.insert(std::make_pair(p, w));
+  gradients.insert(std::make_pair(p, Vector_2(-g,g) ));
+   
+  sample_3.push_back( Point_3(h/3,h/3,w)); 
+  sample_3.push_back( Point_3(-h/3,-h/3,w));
+  sample_3.push_back( Point_3(-h/3,h/3,w));
+  sample_3.push_back( Point_3(h/3,-h/3,w));
+ 
+  //Interpolation grid:
+  Point_vector_2 points;
+  generate_grid_points(points, m, 0.999* h); 
+  Point_vector_3 points_3;
+  
+  int method;
+  std::cout << "Enter the choice of the interpolation function: "<<std::endl;
+  std::cout << " 0 -- linear interpolation" <<std::endl;  
+  std::cout << " 1 -- simple quadratic interpolant"  <<std::endl; 
+  std::cout << " 2 -- Sibson's C1 interpolant "<<std::endl;
+  std::cout << " 3 -- Sibson's C1 interpolant without square-root"<<std::endl;
+  std::cout << " 4 -- Farin's C1 interpolant" << std::endl;
+  std::cin >> method;
+  
+
+  //INTERPOLATION:
+  Coord_type exact_value, res, norm;
+  std::vector< std::pair< Point_2, Coord_type > > coords;
+  int n = points.size();
+  ITraits traits;
+  
+  std::cout << "Interpolation at  "<<n  <<" grid points " << std::endl;
+  for(int i=0;i<n;i++){
+    norm = 
+      CGAL::natural_neighbor_coordinates_2(T, points[i],
+					   std::back_inserter(coords)).second;
+    assert(norm>0);  
+    
+    switch(method){
+    case 0: 
+      res = CGAL::linear_interpolation(coords.begin(),coords.end(),norm,
+				       DataAccess<Point_value_map>(values)); 
+      break;
+    case 1: 
+      res = CGAL::quadratic_interpolation(coords.begin(),coords.end(),
+					  norm, points[i], 
+					  DataAccess< Point_value_map>
+					  (values),
+					  DataAccess< Point_vector_map>
+					  (gradients),traits); break;
+    case 2:  
+      res = CGAL::sibson_c1_interpolation(coords.begin(),coords.end(),
+					  norm, points[i], 
+					  DataAccess<Point_value_map>
+					  (values),
+					  DataAccess<Point_vector_map>
+					  (gradients), traits); break;
+    case 3:  
+      res = CGAL::sibson_c1_interpolation_square(coords.begin(),coords.end(),
+						 norm, points[i], 
+						 DataAccess<Point_value_map>
+						 (values),
+						 DataAccess<Point_vector_map>
+						 (gradients), traits); break;
+    case 4:  
+      res = CGAL::farin_c1_interpolation(coords.begin(),coords.end(),
+					 norm, points[i], 
+					 DataAccess<Point_value_map>
+					 (values),
+					 DataAccess<Point_vector_map>
+					 (gradients), traits); break;
+      
+    default: std::cout <<"No valid choice of interpolant." <<
+		  std::endl; break;
+    }
+    points_3.push_back(Point_3(points[i].x(), points[i].y(),res));
+    coords.clear();
+  }
+  /************** end of Coordinate computation **************/
+  dump_off_file_quadrilateral(points_3, m);
+
+  char ch;
+  
+  //viewer
+  CGAL::Geomview_stream gv(CGAL::Bbox_3(0,0,0, 2, 2, 2));
+  gv.set_bg_color(CGAL::Color(0, 200, 200));
+  gv.clear();
+  gv.set_line_width(2);
+  
+  std::cout << "The data points are displayed in blue in the geomview"
+	    << " application." << std::endl;
+  gv << CGAL::BLUE;
+  visu_points(gv,sample_3);
+  
+  //show the gradients
+  if(method>0){
+    std::cout << "The function gradients are displayed by red lines "
+	      <<" in the geomview application." << std::endl;
+    gv <<CGAL::RED;
+    gv << Segment_3(Point_3(h/3,h/3,w),Point_3(h/3,h/3,w)+ Vector_3(-g,-g,0));
+    gv << Segment_3(Point_3(-h/3,h/3,w),Point_3(-h/3,h/3,w)+Vector_3(g,-g,0));
+    gv << Segment_3(Point_3(-h/3,-h/3,w),Point_3(-h/3,-h/3,w)+Vector_3(g,g,0));
+    gv << Segment_3(Point_3(h/3,-h/3,w),Point_3(h/3,-h/3,w)+Vector_3(-g,g,0));
+  }
+  std::cout << "You should open the file 'result.off' with geomview."
+	    << std::endl;
+  std::cout << "'result.off' contains the graph of the interpolation"
+	    << " function over a grid using ";
+  switch(method){
+  case 0: std::cout << " linear interpolation"; break;
+  case 1: std::cout << " simple quadratic interpolant"; break;
+  case 2: std::cout << " Sibson's C1 interpolant";break;
+  case 3: std::cout << " Sibson's C1 interpolant without square-root"; break;
+  case 4: std::cout << " Farin's C1 interpolant"; break;
+  }
+  std::cout << std::endl;
+  std::cout << "Enter any character to quit." << std::endl;
+  std::cin >> ch;
+  
+  return 1; 
+}
+#endif // if defined(__BORLANDC__) || defined(_MSC_VER)
+ 
diff --git a/Packages/Interpolation/demo/Interpolation/demo_numerical_2.C b/Packages/Interpolation/demo/Interpolation/demo_numerical_2.C
new file mode 100644
index 00000000000..c210769100e
--- /dev/null
+++ b/Packages/Interpolation/demo/Interpolation/demo_numerical_2.C
@@ -0,0 +1,239 @@
+//file: demo/Interpolation/demo_numerical_2.C
+// compares the result of several interpolation methods  
+// Author: Julia
+#include <CGAL/basic.h>
+#include <cstdio>
+#include <utility>
+
+#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
+
+#include <CGAL/Random.h>
+#include <CGAL/squared_distance_2.h>
+
+#include <CGAL/Delaunay_triangulation_2.h>
+#include <CGAL/Interpolation_traits_2.h>
+#include <CGAL/natural_neighbor_coordinates_2.h>
+#include <CGAL/interpolation_functions.h>
+
+#include <CGAL/point_generators_2.h> 
+#include <CGAL/copy_n.h>
+#include <CGAL/Origin.h>
+
+
+struct K : CGAL::Exact_predicates_exact_constructions_kernel {};
+typedef K::FT                                         Coord_type;
+typedef K::Vector_2                                   Vector;
+typedef K::Point_2                                    Point;
+
+typedef CGAL::Delaunay_triangulation_2<K>             Delaunay_triangulation;
+typedef CGAL::Interpolation_traits_2<K>               Traits;
+
+typedef std::vector< std::pair<Point, Coord_type> >   Coordinate_vector;
+typedef std::map<Point, Coord_type, K::Less_xy_2>     Point_value_map;
+typedef std::map<Point,  Vector, K::Less_xy_2 >       Point_vector_map;
+
+//////////////////////////////// 
+// functor class that provides access to the function values/ function
+//   gradients given a data point (using "map"):
+template< class Map >  
+struct DataAccess : public std::unary_function< typename Map::key_type,
+		    typename Map::mapped_type> {
+  typedef typename Map::mapped_type Data_type;
+  typedef typename Map::key_type  Point;
+  
+  //CONSTRUCTOR:
+  DataAccess< Map >(const Map& m): map(m){};
+  
+  //Functor
+  Data_type operator()(const Point& p) { 
+    
+    typename Map::const_iterator mit = map.find(p);
+    if(mit!= map.end())
+      return mit->second;
+    return Data_type();
+  };
+  
+  const Map& map;
+};
+
+
+//////////////test function:///////////////////////////////
+template < class Iterator > 
+void test_barycenter(Iterator first, Iterator beyond,
+		     typename std::iterator_traits<Iterator>::value_type::
+		     second_type norm,
+		     const typename std::iterator_traits<Iterator>::
+		     value_type::first_type& p)
+{
+  typedef typename
+    std::iterator_traits<Iterator>::value_type::first_type Point_2;
+  Point_2 b(0,0);
+  
+  for(; first !=beyond; first++)
+    b = b+ (first->second/norm) * (first->first - CGAL::ORIGIN);
+    
+  std::cout <<"Distance between test point and barycenter " <<
+    CGAL::squared_distance(p, b ) << " egality: "<< (p==b)<< std::endl;
+}
+
+////////////////////////////////
+
+int main()
+{
+  //number of sample points:
+  int n = 24;
+  //number of interpolation points:
+  int m = 20;
+
+  std::vector<Point> points;
+  points.reserve(n+m);
+
+  //put four bounding box points:
+  points.push_back(Point(-3,-3));
+  points.push_back(Point(3,-3));
+  points.push_back(Point(-3,3));
+  points.push_back(Point(3,3));
+
+  // Create n+m-4 points within a disc of radius 2
+  double r_d = 3;
+  CGAL::Random_points_in_disc_2<Point> g(r_d );
+  CGAL::copy_n( g, n+m, std::back_inserter(points));
+  
+  Delaunay_triangulation T;
+
+
+  Point_value_map values;
+  Point_vector_map gradients;
+
+
+  //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.2),
+    gamma2 = Coord_type(0.7),  
+    gamma3 = Coord_type(0.4),  
+    gamma4 = Coord_type(0.6);
+
+  for(int j=0; j<n ; j++){
+    T.insert(points[j]);
+  
+    //determine function value/gradient:
+    Coord_type x(points[j].x());
+    Coord_type y(points[j].y());
+      
+    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));
+    values.insert(std::make_pair(points[j], value));
+    gradients.insert(std::make_pair(points[j], gradient));
+  }
+
+  //variables for statistics:
+  Coord_type res, error, l_total = Coord_type(0),
+    q_total(l_total), f_total(l_total), s_total(l_total),
+    l_max(l_total), 
+    q_max(l_total), f_max(l_total), s_max(l_total), 
+    total_value(l_total);
+  int failure(0);
+
+  //interpolation + error statistics
+  for(int i=n;i<n+m;i++){
+    Coord_type x(points[i].x());
+    Coord_type y(points[i].y());
+      
+    Coord_type exact_value = alpha + beta1*x + beta2*y + gamma1*(x*x) + 
+      gamma4*(y*y) + (gamma2+ gamma3) *(x*y);
+  
+    total_value += exact_value;
+    
+    //Coordinate_vector:
+    std::vector< std::pair< Point, Coord_type > > coords;
+    Coord_type norm = 
+      CGAL::natural_neighbor_coordinates_2(T, points[i],
+					   std::back_inserter(coords)).second;
+    
+    if(norm>0){
+      //test barycenter
+      test_barycenter( coords.begin(), coords.end(),norm, points[i]);
+      
+      
+      //linear interpolant:
+      res =  CGAL::linear_interpolation(coords.begin(), coords.end(),
+					norm,
+					DataAccess<Point_value_map>(values));
+      
+      error = CGAL_NTS abs(res - 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], 
+					  DataAccess<Point_value_map>(values), 
+					  DataAccess<Point_vector_map>
+					  (gradients), 
+					  Traits());
+      error = CGAL_NTS abs(res - exact_value);
+      f_total += error;
+      if (error > f_max) f_max = error;
+      
+      //quadratic interpolant:
+      res =  CGAL::quadratic_interpolation(coords.begin(), coords.end(),
+					   norm,points[i], 
+					   DataAccess<Point_value_map>
+					   (values), 
+					   DataAccess<Point_vector_map>
+					   (gradients), 
+					   Traits());
+      error = CGAL_NTS abs(res - exact_value);
+      q_total += error;
+      if (error > q_max) q_max = error;
+    
+      //Sibson interpolant: version without sqrt:
+      res =  CGAL::sibson_c1_interpolation_square(coords.begin(),
+						  coords.end(), norm, 
+						  points[i], 
+						  DataAccess<Point_value_map>
+						  (values), 
+						  DataAccess<Point_vector_map>
+						  (gradients),
+						  Traits());
+      //error statistics
+      error = CGAL_NTS abs(res - exact_value);
+      s_total += error;
+      if (error > s_max) s_max = error;
+      
+    }else  failure++;
+  }
+  
+  /************** end of Interpolation: dump statistics **************/
+  std::cout << "Result: -----------------------------------" << std::endl;
+  std::cout <<  "Interpolation of '" << alpha <<" + " 
+	    << beta1<<" x + " 
+	    << beta2 << " y + " << gamma1 <<" x^2 + "  << gamma2+ gamma3
+	    <<" xy + "  << gamma4 << " y^2'" << std::endl; 
+  std::cout << "Knowing " << m << " sample points. Interpolation on "
+	    << n <<" random points. "<< std::endl; 
+  std::cout <<"Average function value "
+	    << (1.0/n) * CGAL::to_double(total_value)
+	    << ", nb of failures "<< failure << std::endl;
+  
+  std::cout << "linear interpolant mean error  " 
+	    << CGAL::to_double(l_total)/n << "  max " 
+	    << CGAL::to_double(l_max) <<std::endl;
+  std::cout << "quadratic interpolant  mean error  " 
+	    << CGAL::to_double(q_total)/n << "  max " 
+	    << CGAL::to_double(q_max) << std::endl;
+  std::cout << "Farin interpolant  mean error  " 
+	    << CGAL::to_double(f_total)/n << "  max " 
+	    << CGAL::to_double(f_max)  << std::endl;
+  std::cout << "Sibson interpolant  mean error  " 
+	    << CGAL::to_double(s_total)/n << "  max " 
+	    << CGAL::to_double(s_max)  << std::endl;  
+ 
+  return 1; 
+}
+
+//end of file
diff --git a/Packages/Interpolation/demo/Interpolation/makefile b/Packages/Interpolation/demo/Interpolation/makefile
new file mode 100644
index 00000000000..d155c1dfb34
--- /dev/null
+++ b/Packages/Interpolation/demo/Interpolation/makefile
@@ -0,0 +1,56 @@
+# Created by the script create_makefile
+# This is the makefile for compiling a CGAL application.
+
+#---------------------------------------------------------------------#
+#                    include platform specific settings
+#---------------------------------------------------------------------#
+# Choose the right include file from the <cgalroot>/make directory.
+
+# CGAL_MAKEFILE = ENTER_YOUR_INCLUDE_MAKEFILE_HERE
+include $(CGAL_MAKEFILE)
+
+#---------------------------------------------------------------------#
+#                    compiler flags
+#---------------------------------------------------------------------#
+
+CXXFLAGS = \
+           $(CGAL_CXXFLAGS) \
+           $(LONG_NAME_PROBLEM_CXXFLAGS) \
+           $(DEBUG_OPT)
+
+#---------------------------------------------------------------------#
+#                    linker flags
+#---------------------------------------------------------------------#
+
+LIBPATH = \
+          $(CGAL_LIBPATH)
+
+LDFLAGS = \
+          $(LONG_NAME_PROBLEM_LDFLAGS) \
+          $(CGAL_GEOWIN_LDFLAGS)
+
+#---------------------------------------------------------------------#
+#                    target entries
+#---------------------------------------------------------------------#
+
+all:            \
+                demo_interpolation_2$(EXE_EXT) \
+                demo_numerical_2$(EXE_EXT) 
+
+demo_interpolation_2$(EXE_EXT): demo_interpolation_2$(OBJ_EXT)
+	$(CGAL_CXX) $(LIBPATH) $(EXE_OPT)demo_interpolation_2 demo_interpolation_2$(OBJ_EXT) $(LDFLAGS)
+
+demo_numerical_2$(EXE_EXT): demo_numerical_2$(OBJ_EXT)
+	$(CGAL_CXX) $(LIBPATH) $(EXE_OPT)demo_numerical_2 demo_numerical_2$(OBJ_EXT) $(LDFLAGS)
+
+clean: \
+                   demo_interpolation_2.clean \
+                   demo_numerical_2.clean 
+
+#---------------------------------------------------------------------#
+#                    suffix rules
+#---------------------------------------------------------------------#
+
+.C$(OBJ_EXT):
+	$(CGAL_CXX) $(CXXFLAGS) $(OBJ_OPT) $<
+
diff --git a/Packages/Interpolation/demo/Interpolation/makefile_local b/Packages/Interpolation/demo/Interpolation/makefile_local
new file mode 100644
index 00000000000..2fea82f1bc9
--- /dev/null
+++ b/Packages/Interpolation/demo/Interpolation/makefile_local
@@ -0,0 +1,59 @@
+# Created by the script create_makefile
+# This is the makefile for compiling a CGAL application.
+
+#---------------------------------------------------------------------#
+#                    include platform specific settings
+#---------------------------------------------------------------------#
+# Choose the right include file from the <cgalroot>/make directory.
+
+# CGAL_MAKEFILE = ENTER_YOUR_INCLUDE_MAKEFILE_HERE
+CGAL_MAKEFILE = /0/prisme_util/CGAL/CGAL-I/make/makefile_i686_Linux-2.2.18_g++-3.3.0
+
+include $(CGAL_MAKEFILE)
+
+#---------------------------------------------------------------------#
+#                    compiler flags
+#---------------------------------------------------------------------#
+
+CXXFLAGS = \
+	   -I ../../include \
+           $(CGAL_CXXFLAGS) \
+           $(LONG_NAME_PROBLEM_CXXFLAGS) \
+           $(DEBUG_OPT)
+
+#---------------------------------------------------------------------#
+#                    linker flags
+#---------------------------------------------------------------------#
+
+LIBPATH = \
+          $(CGAL_LIBPATH)
+
+LDFLAGS = \
+          $(LONG_NAME_PROBLEM_LDFLAGS) \
+          $(CGAL_GEOWIN_LDFLAGS)
+
+#---------------------------------------------------------------------#
+#                    target entries
+#---------------------------------------------------------------------#
+
+all:            \
+                demo_interpolation_2$(EXE_EXT) \
+                demo_numerical_2$(EXE_EXT) 
+
+demo_interpolation_2$(EXE_EXT): demo_interpolation_2$(OBJ_EXT)
+	$(CGAL_CXX) $(LIBPATH) $(EXE_OPT)demo_interpolation_2 demo_interpolation_2$(OBJ_EXT) $(LDFLAGS)
+
+demo_numerical_2$(EXE_EXT): demo_numerical_2$(OBJ_EXT)
+	$(CGAL_CXX) $(LIBPATH) $(EXE_OPT)demo_numerical_2 demo_numerical_2$(OBJ_EXT) $(LDFLAGS)
+
+clean: \
+                   demo_interpolation_2.clean \
+                   demo_numerical_2.clean 
+
+#---------------------------------------------------------------------#
+#                    suffix rules
+#---------------------------------------------------------------------#
+
+.C$(OBJ_EXT):
+	$(CGAL_CXX) $(CXXFLAGS) $(OBJ_OPT) $<
+
diff --git a/Packages/Interpolation/examples/Interpolation/README b/Packages/Interpolation/examples/Interpolation/README
new file mode 100644
index 00000000000..bae2eb3f208
--- /dev/null
+++ b/Packages/Interpolation/examples/Interpolation/README
@@ -0,0 +1,17 @@
+To compule and run all these examples type :
+./cgal_test
+To compute and run only some of them type
+./cgal_test name-of_wanted_example
+ 
+nn_coordinates_2: shows how to compute 2D natural coordinates 
+	given a 2D Delaunay triangulation  and a query point.
+
+rn_coordinates_2: shows how to compute 2D regular (natural) coordinates 
+	given a 2D Regular triangulation  and a (weighted) query point.
+
+linear_interpolation_2: interpolates a linear function using the 
+	linear_interpolation function and 2D natural neighbor coordinates.
+
+sibson_interpolation_2: interpolates a spherical function using the 
+	sibson_gradient_fitting and sibson_c1_interpolation function with 
+	2D natural neighbor coordinates.
\ No newline at end of file
diff --git a/Packages/Interpolation/examples/Interpolation/linear_interpolation_2.C b/Packages/Interpolation/examples/Interpolation/linear_interpolation_2.C
new file mode 100644
index 00000000000..445b06c6909
--- /dev/null
+++ b/Packages/Interpolation/examples/Interpolation/linear_interpolation_2.C
@@ -0,0 +1,67 @@
+//
+//file: examples/Interpolation/linear_interoplation_2.C 
+//
+#include <CGAL/basic.h>
+#include <utility>
+
+#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
+#include <CGAL/Delaunay_triangulation_2.h>
+
+#include <CGAL/Interpolation_traits_2.h>
+#include <CGAL/natural_neighbor_coordinates_2.h>
+#include <CGAL/interpolation_functions.h>
+
+struct K : CGAL::Exact_predicates_inexact_constructions_kernel {};
+typedef CGAL::Delaunay_triangulation_2<K>             Delaunay_triangulation;
+typedef CGAL::Interpolation_traits_2<K>               Traits;
+typedef K::FT                                         Coord_type;
+
+//Functor for accessing the function values
+template< class Map >  
+struct DataAccess : public std::unary_function< typename Map::key_type,
+		    typename Map::mapped_type> {
+  typedef typename Map::mapped_type Data_type;
+  typedef typename Map::key_type  Point;
+  
+  DataAccess< Map >(const Map& m): map(m){};
+ 
+  Data_type operator()(const Point& p) { 
+    
+    typename Map::const_iterator mit = map.find(p);
+    if(mit!= map.end())
+      return mit->second;
+    return Data_type();
+  };
+  
+  const Map& map;
+};
+
+
+int main()
+{
+  Delaunay_triangulation T;
+  std::map<Point, Coord_type, K::Less_xy_2> values;
+  typedef DataAccess< std::map<Point, Coord_type, K::Less_xy_2 > > 
+                                                        Value_access;
+  
+  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);
+      values.insert(std::make_pair(p,a + bx* x+ by*y));
+    }
+  //coordiante 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;  
+ 
+  Coord_type res =  CGAL::linear_interpolation(coords.begin(), coords.end(), 
+				               norm,Value_access(values));
+  
+  std::cout << "   Tested interpolation on " << p << " interpolation: " << res 
+	    << " exact: " << a + bx* p.x()+ by* p.y()<< std::endl;
+  return 1; 
+}
diff --git a/Packages/Interpolation/examples/Interpolation/nn_coordinates_2.C b/Packages/Interpolation/examples/Interpolation/nn_coordinates_2.C
new file mode 100644
index 00000000000..9f51f6095d3
--- /dev/null
+++ b/Packages/Interpolation/examples/Interpolation/nn_coordinates_2.C
@@ -0,0 +1,29 @@
+//
+//file: examples/Interpolation/nn_coordinates_2.C 
+//
+#include <CGAL/basic.h>
+#include <utility>
+
+#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
+#include <CGAL/Delaunay_triangulation_2.h>
+#include <CGAL/natural_neighbor_coordinates_2.h>
+
+struct K : CGAL::Exact_predicates_inexact_constructions_kernel {};
+typedef CGAL::Delaunay_triangulation_2<K>             Delaunay_triangulation;
+
+int main()
+{
+  Delaunay_triangulation T;
+  
+  for (int y=0 ; y<3 ; y++)
+    for (int x=0 ; x<3 ; x++) 
+      T.insert(K::Point_2(x,y));
+   
+  //coordinate computation
+  K::Point_2 p(1.2, 0.7);
+  std::vector< std::pair< K::Point_2, K::FT  > > coords;
+  K::FT  norm = 
+    CGAL::natural_neighbor_coordinates_2(T, p,std::back_inserter(coords)).second;
+  
+  return 1; 
+}
diff --git a/Packages/Interpolation/examples/Interpolation/rn_coordinates_2.C b/Packages/Interpolation/examples/Interpolation/rn_coordinates_2.C
new file mode 100644
index 00000000000..ba42b66fd0f
--- /dev/null
+++ b/Packages/Interpolation/examples/Interpolation/rn_coordinates_2.C
@@ -0,0 +1,34 @@
+//
+//file: examples/Interpolation/rn_coordinates_2.C 
+//
+#include <CGAL/basic.h>
+#include <utility>
+
+#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
+#include <CGAL/Regular_triangulation_2.h>
+#include <CGAL/Regular_neighbor_coordinates_traits_2.h>
+#include <CGAL/regular_neighbor_coordinates_2.h>
+
+struct K : CGAL::Exact_predicates_inexact_constructions_kernel {};
+
+typedef CGAL::Regular_neighbor_coordinates_traits_2<K>     Gt;
+typedef CGAL::Regular_triangulation_2<Gt>     Regular_triangulation;
+typedef Regular_triangulation::Weighted_point Weighted_point;
+
+int main()
+{
+  Regular_triangulation rt;
+  
+  for (int y=0 ; y<3 ; y++)
+    for (int x=0 ; x<3 ; x++) 
+      rt.insert(Weighted_point(K::Point_2(x,y), 0));
+   
+  //coordinate computation
+  Weighted_point wp(K::Point_2(1.2, 0.7),2);
+  
+  std::vector< std::pair< K::Point_2, K::FT  > > coords;
+  K::FT  norm = 
+    CGAL::regular_neighbor_coordinates_2(rt, wp,std::back_inserter(coords)).second;
+  
+  return 1; 
+}
diff --git a/Packages/Interpolation/examples/Interpolation/sibson_interpolation_2.C b/Packages/Interpolation/examples/Interpolation/sibson_interpolation_2.C
new file mode 100644
index 00000000000..414681977aa
--- /dev/null
+++ b/Packages/Interpolation/examples/Interpolation/sibson_interpolation_2.C
@@ -0,0 +1,81 @@
+//
+//file:  examples/Interpolation/sibson_interoplation_2.C 
+//
+#include <CGAL/basic.h>
+#include <utility>
+
+#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
+#include <CGAL/Delaunay_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>
+
+struct K : CGAL::Exact_predicates_inexact_constructions_kernel {};
+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;
+
+//Functor class for accessing the function values/gradients
+template< class Map >  
+struct DataAccess : public std::unary_function< typename Map::key_type,
+		    typename Map::mapped_type> {
+  typedef typename Map::mapped_type Data_type;
+  typedef typename Map::key_type  Point;
+  
+  DataAccess< Map >(const Map& m): map(m){};
+  
+  Data_type operator()(const Point& p) { 
+    typename Map::const_iterator mit = map.find(p);
+    if(mit!= map.end())
+      return mit->second;
+    return Data_type();
+  };
+  
+  const Map& map;
+};
+
+int main()
+{
+  
+  Delaunay_triangulation T;
+
+  Point_value_map values;
+  Point_vector_map gradients;
+
+  //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++){ 
+      K::Point_2 p(x,y);
+      T.insert(p);
+      values.insert(std::make_pair(p,a + bx* x+ by*y + c*(x*x+y*y)));
+    }
+  sibson_gradient_fitting_nn_2(T,std::inserter(gradients,gradients.begin()),
+			       DataAccess<Point_value_map>(values), Traits());
+  
+ 
+  //coordiante 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;
+  
+  
+  //Sibson interpolant: version without sqrt:
+  Coord_type res =  CGAL::sibson_c1_interpolation_square(coords.begin(),
+							 coords.end(),norm,p, 
+							 DataAccess<Point_value_map>(values),
+							 DataAccess<Point_vector_map>(gradients),
+							 Traits());
+  std::cout << "   Tested interpolation on " << p << " interpolation: " << res 
+	    << " exact: " 
+	    << a + bx * p.x()+ by * p.y()+ c*(p.x()*p.x()+p.y()*p.y()) 
+	    << std::endl;
+  return 1; 
+}