diff --git a/Interpolation/examples/Interpolation/linear_interpolation_of_vector_3.cpp b/Interpolation/examples/Interpolation/linear_interpolation_of_vector_3.cpp
new file mode 100644
index 00000000000..590dac81f34
--- /dev/null
+++ b/Interpolation/examples/Interpolation/linear_interpolation_of_vector_3.cpp
@@ -0,0 +1,86 @@
+#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>
+
+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::Vector_3                                         Vector_3;
+typedef K::Point_2                                          Point_2;
+
+
+namespace boost {
+  
+template<class K>
+class value_initialized<CGAL::Vector_3<K> >
+{
+  private :
+  
+    typedef CGAL::Vector_3<K> T;
+    T m_data;
+
+  public :
+    
+    value_initialized()
+    :
+      m_data(0,0,0)
+    { }
+
+    T const & data() const
+    {
+      return m_data;
+    }
+  
+    T& data()
+    {
+      return m_data;
+    }
+
+    operator T const &() const
+    {
+      return m_data;
+    }
+
+    operator T&()
+    {
+      return m_data;
+    }
+} ;
+
+} // nammespace boost
+
+
+int main()
+{
+  Delaunay_triangulation T;
+
+  typedef std::map<Point_2, Vector_3, K::Less_xy_2>           Coord_map;
+  typedef CGAL::Data_access<Coord_map>                      Value_access;
+
+  Coord_map value_function;
+  double a(0.25), bx(1.3), by(-0.7);
+
+  for (int y=0 ; y<255 ; y++){
+    for (int x=0 ; x<255 ; x++){
+      K::Point_2_2 p(x,y);
+      T.insert(p);
+      value_function.insert(std::make_pair(p, Vector_3(x,y,1)));
+    }
+  }
+
+  //coordinate computation
+  K::Point_2_2 p(1.3, 0.34);
+  std::vector<std::pair<Point_2, double> > coords;
+
+  double norm = CGAL::natural_neighbor_coordinates_2(T, p, std::back_inserter(coords)).second;
+  Vector_3 res =  CGAL::linear_interpolation(coords.begin(), coords.end(), norm,
+                                           Value_access(value_function));
+
+  std::cout << "Tested interpolation on " << p << " interpolation: "
+            << res << std::endl;
+
+  return EXIT_SUCCESS;
+}
diff --git a/Interpolation/include/CGAL/interpolation_functions.h b/Interpolation/include/CGAL/interpolation_functions.h
index cac67d5aa40..8d5e4565425 100644
--- a/Interpolation/include/CGAL/interpolation_functions.h
+++ b/Interpolation/include/CGAL/interpolation_functions.h
@@ -28,7 +28,7 @@
 #include <CGAL/use.h>
 
 #include <boost/utility/result_of.hpp>
-
+#include <boost/utility/value_init.hpp> 
 #include <iterator>
 #include <utility>
 #include <vector>
@@ -59,12 +59,13 @@ struct Data_access
 };
 
 //the interpolation functions:
-template < class ForwardIterator, class ValueFunctor >
+  template < class ForwardIterator, class ValueFunctor>
 typename boost::result_of<
            ValueFunctor(typename std::iterator_traits<ForwardIterator>::value_type::first_type)>
              ::type::first_type
 linear_interpolation(ForwardIterator first, ForwardIterator beyond,
                      const typename std::iterator_traits<ForwardIterator>::value_type::second_type& norm,
+                     //const Norm& norm,
                      ValueFunctor value_function)
 {
   CGAL_precondition(norm > 0);
@@ -73,13 +74,13 @@ linear_interpolation(ForwardIterator first, ForwardIterator beyond,
   typedef typename boost::result_of<ValueFunctor(arg_type)>::type                 result_type;
   typedef typename result_type::first_type                                        Value_type;
 
-  Value_type result(0);
+  boost::value_initialized<Value_type> result;
   result_type val;
   for(; first!=beyond; ++first)
   {
     val = value_function(first->first);
     CGAL_assertion(val.second);
-    result += (first->second / norm) * val.first;
+    get(result) += (first->second / norm) * val.first;
   }
   return result;
 }