diff --git a/Surface_mesh_approximation/examples/Surface_mesh_approximation/vsa_example.cpp b/Surface_mesh_approximation/examples/Surface_mesh_approximation/vsa_example.cpp
index 63e5e22acc4..472c2896a50 100644
--- a/Surface_mesh_approximation/examples/Surface_mesh_approximation/vsa_example.cpp
+++ b/Surface_mesh_approximation/examples/Surface_mesh_approximation/vsa_example.cpp
@@ -19,11 +19,13 @@ typedef Polyhedron::Halfedge_const_handle Halfedge_const_handle;
 typedef Polyhedron::Facet_const_iterator Facet_const_iterator;
 typedef boost::associative_property_map<std::map<Facet_const_handle, Vector> > FacetNormalMap;
 typedef boost::associative_property_map<std::map<Facet_const_handle, FT> > FacetAreaMap;
+typedef boost::property_map<Polyhedron, boost::vertex_point_t>::type VertexPointMap;
 
 typedef CGAL::PlaneProxy<Polyhedron> PlaneProxy;
 typedef CGAL::L21Metric<PlaneProxy, FacetNormalMap, FacetAreaMap> L21Metric;
 typedef CGAL::L21ProxyFitting<PlaneProxy, L21Metric, FacetNormalMap, FacetAreaMap> L21ProxyFitting;
-typedef CGAL::L21ApproximationTrait<PlaneProxy, L21Metric, L21ProxyFitting, FacetNormalMap, FacetAreaMap> L21ApproximationTrait;
+typedef CGAL::PlaneFitting<Polyhedron, FacetAreaMap, FacetNormalMap, VertexPointMap> PlaneFitting;
+typedef CGAL::L21ApproximationTrait<PlaneProxy, Polyhedron, L21Metric, L21ProxyFitting, PlaneFitting, VertexPointMap, FacetNormalMap, FacetAreaMap> L21ApproximationTrait;
 
 int main(int argc, char *argv[])
 {
@@ -53,6 +55,7 @@ int main(int argc, char *argv[])
   }
   FacetNormalMap normal_pmap(facet_normals);
   FacetAreaMap area_pmap(facet_areas);
+  VertexPointMap point_pmap = get(boost::vertex_point, const_cast<Polyhedron &>(mesh));
 
   // create a property-map for segment-ids
   typedef std::map<Facet_const_handle, std::size_t> Facet_id_map;
@@ -74,13 +77,13 @@ int main(int argc, char *argv[])
     num_proxies,
     num_iterations,
     proxy_patch_map,
-    get(boost::vertex_point, const_cast<Polyhedron &>(mesh)),
+    point_pmap,
     area_pmap,
     tris,
     anchor_pos,
     anchor_vtx,
     bdrs,
-    L21ApproximationTrait(normal_pmap, area_pmap),
+    L21ApproximationTrait(mesh, point_pmap, normal_pmap, area_pmap),
     Kernel());
 
   return EXIT_SUCCESS;
diff --git a/Surface_mesh_approximation/examples/Surface_mesh_approximation/vsa_metric_example.cpp b/Surface_mesh_approximation/examples/Surface_mesh_approximation/vsa_metric_example.cpp
index c75ea245d90..a3be094e3a7 100644
--- a/Surface_mesh_approximation/examples/Surface_mesh_approximation/vsa_metric_example.cpp
+++ b/Surface_mesh_approximation/examples/Surface_mesh_approximation/vsa_metric_example.cpp
@@ -21,6 +21,7 @@ typedef Polyhedron::Facet_const_iterator Facet_const_iterator;
 typedef boost::associative_property_map<std::map<Facet_const_handle, Vector> > FacetNormalMap;
 typedef boost::associative_property_map<std::map<Facet_const_handle, FT> > FacetAreaMap;
 typedef boost::associative_property_map<std::map<Facet_const_handle, Point> > FacetCenterMap;
+typedef boost::property_map<Polyhedron, boost::vertex_point_t>::type VertexPointMap;
 
 struct PointProxy {
   Facet_handle seed;
@@ -100,11 +101,14 @@ struct ApproxTrait {
   typedef PointProxy Proxy;
   typedef CompactMetric ErrorMetric;
   typedef PointProxyFitting ProxyFitting;
+  typedef CGAL::PlaneFitting<Polyhedron, FacetAreaMap, FacetNormalMap, VertexPointMap> PlaneFitting;
 
-  ApproxTrait(const FacetCenterMap _center_pmap,
+  ApproxTrait(const Polyhedron &_mesh,
+    const VertexPointMap &_point_pmap,
+    const FacetCenterMap &_center_pmap,
     const FacetAreaMap &_area_pmap,
     const FacetNormalMap &_normal_pmap)
-    : center_pmap(_center_pmap), area_pmap(_area_pmap), normal_pmap(_normal_pmap) {}
+    : mesh(_mesh), point_pmap(_point_pmap), center_pmap(_center_pmap), area_pmap(_area_pmap), normal_pmap(_normal_pmap) {}
 
   ErrorMetric construct_fit_error_functor() const {
     return ErrorMetric(center_pmap);
@@ -114,6 +118,12 @@ struct ApproxTrait {
     return ProxyFitting(center_pmap, area_pmap, normal_pmap);
   }
 
+  PlaneFitting construct_plane_fitting_functor() const {
+    return PlaneFitting(mesh, area_pmap, normal_pmap, point_pmap);
+  }
+
+  const Polyhedron &mesh;
+  const VertexPointMap point_pmap;
   const FacetCenterMap center_pmap;
   const FacetAreaMap area_pmap;
   const FacetNormalMap normal_pmap;
@@ -150,6 +160,7 @@ int main(int argc, char *argv[])
   FacetNormalMap normal_pmap(facet_normals);
   FacetAreaMap area_pmap(facet_areas);
   FacetCenterMap center_pmap(facet_centers);
+  VertexPointMap point_pmap = get(boost::vertex_point, const_cast<Polyhedron &>(mesh));
 
   // create a property-map for segment-ids
   typedef std::map<Facet_const_handle, std::size_t> Facet_id_map;
@@ -171,13 +182,13 @@ int main(int argc, char *argv[])
     num_proxies,
     num_iterations,
     proxy_patch_map,
-    get(boost::vertex_point, const_cast<Polyhedron &>(mesh)),
+    point_pmap,
     area_pmap,
     tris,
     anchor_pos,
     anchor_vtx,
     bdrs,
-    ApproxTrait(center_pmap, area_pmap, normal_pmap),
+    ApproxTrait(mesh, point_pmap, center_pmap, area_pmap, normal_pmap),
     Kernel());
 
   return EXIT_SUCCESS;
diff --git a/Surface_mesh_approximation/include/CGAL/internal/Surface_mesh_approximation/VSA.h b/Surface_mesh_approximation/include/CGAL/internal/Surface_mesh_approximation/VSA.h
index a56707acea8..7408271014d 100644
--- a/Surface_mesh_approximation/include/CGAL/internal/Surface_mesh_approximation/VSA.h
+++ b/Surface_mesh_approximation/include/CGAL/internal/Surface_mesh_approximation/VSA.h
@@ -46,6 +46,7 @@ private:
   typedef typename ApproximationTrait::Proxy Proxy;
   typedef typename ApproximationTrait::ErrorMetric ErrorMetric;
   typedef typename ApproximationTrait::ProxyFitting ProxyFitting;
+  typedef typename ApproximationTrait::PlaneFitting PlaneFitting;
 
   typedef typename GeomTraits::FT FT;
   typedef typename GeomTraits::Point_3 Point;
@@ -126,6 +127,9 @@ private:
   // The proxy fitting functor.
   ProxyFitting proxy_fitting;
 
+  // The proxy plane fitting functor.
+  PlaneFitting plane_fitting;
+
   //member functions
 public:
   /**
@@ -144,7 +148,8 @@ public:
     area_pmap(_facet_area_map),
     vertex_status_pmap(vertex_status_map),
     fit_error(_appx_trait.construct_fit_error_functor()),
-    proxy_fitting(_appx_trait.construct_proxy_fitting_functor()) {
+    proxy_fitting(_appx_trait.construct_proxy_fitting_functor()),
+    plane_fitting(_appx_trait.construct_plane_fitting_functor()) {
 
     GeomTraits traits;
     vector_functor = traits.construct_vector_3_object();
@@ -232,7 +237,13 @@ public:
 
     pseudo_CDT(seg_pmap, tris);
 
-    compute_anchor_position(seg_pmap);
+    std::vector<Plane> px_planes(proxies.size());
+    std::vector<std::list<face_descriptor> > px_facets(proxies.size());
+    BOOST_FOREACH(face_descriptor f, faces(mesh))
+      px_facets[seg_pmap[f]].push_back(f);
+    for (std::size_t i = 0; i < proxies.size(); ++i)
+      px_planes[i] = plane_fitting(px_facets[i].begin(), px_facets[i].end());
+    compute_anchor_position(seg_pmap, px_planes);
     std::vector<Point> vtx;
     BOOST_FOREACH(const Anchor &a, anchors)
       vtx.push_back(a.pos);
@@ -1073,11 +1084,11 @@ private:
    * @param 
    */
   template<typename FacetSegmentMap>
-  void compute_anchor_position(const FacetSegmentMap &seg_pmap) {
-    std::vector<Point> proxies_center; // The proxy center.
+  void compute_anchor_position(
+    const FacetSegmentMap &seg_pmap,
+    const std::vector<Plane> &px_planes) {
     std::vector<FT> proxies_area; // The proxy area.
     compute_proxy_area(seg_pmap, proxies_area);
-    compute_proxy_center(seg_pmap, proxies_center);
 
     anchors = std::vector<Anchor>(anchor_index);
     BOOST_FOREACH(vertex_descriptor v, vertices(mesh)) {
@@ -1095,8 +1106,7 @@ private:
         for (std::set<std::size_t>::iterator pxitr = px_set.begin();
           pxitr != px_set.end(); ++pxitr) {
           std::size_t px_idx = *pxitr;
-          Plane px_plane(proxies_center[px_idx], proxies[px_idx].normal);
-          Point proj = px_plane.projection(vtx_pt);
+          Point proj = px_planes[px_idx].projection(vtx_pt);
           FT area = proxies_area[px_idx];
           avgx += proj.x() * area;
           avgy += proj.y() * area;
diff --git a/Surface_mesh_approximation/include/CGAL/vsa_mesh_approximation_traits.h b/Surface_mesh_approximation/include/CGAL/vsa_mesh_approximation_traits.h
index 9d3a803db4c..16831d90e9d 100644
--- a/Surface_mesh_approximation/include/CGAL/vsa_mesh_approximation_traits.h
+++ b/Surface_mesh_approximation/include/CGAL/vsa_mesh_approximation_traits.h
@@ -164,21 +164,32 @@ template<typename PlaneProxy,
 
 template<typename TriangleMesh,
   typename FacetAreaMap,
+  typename FacetNormalMap,
   typename VertexPointMap,
   typename GeomTraits = typename TriangleMesh::Traits>
   struct PlaneFitting
 {
   PlaneFitting(const TriangleMesh &_mesh,
     const FacetAreaMap &_area_pmap,
+    const FacetNormalMap &_facet_normal_map,
     const VertexPointMap &_point_pmap)
     : mesh(_mesh),
     area_pmap(_area_pmap),
-    point_pmap(_point_pmap) {}
+    normal_pmap(_facet_normal_map),
+    point_pmap(_point_pmap) {
+      GeomTraits traits;
+      vector_functor = traits.construct_vector_3_object();
+      sum_functor = traits.construct_sum_of_vectors_3_object();
+      scale_functor = traits.construct_scaled_vector_3_object();
+    }
 
   typedef typename GeomTraits::FT FT;
   typedef typename GeomTraits::Point_3 Point_3;
   typedef typename GeomTraits::Vector_3 Vector_3;
   typedef typename GeomTraits::Plane_3 Plane_3;
+  typedef typename GeomTraits::Construct_vector_3 Construct_vector_3;
+  typedef typename GeomTraits::Construct_scaled_vector_3 Construct_scaled_vector_3;
+  typedef typename GeomTraits::Construct_sum_of_vectors_3 Construct_sum_of_vectors_3;
   typedef typename boost::graph_traits<TriangleMesh>::halfedge_descriptor halfedge_descriptor;
 
   template<typename FacetIterator>
@@ -197,24 +208,29 @@ template<typename TriangleMesh,
     FT sum_area(0);
     Vector_3 centroid = CGAL::NULL_VECTOR;
     for (FacetIterator fitr = beg; fitr != end; ++fitr) {
-      const halfedge_descriptor he = halfedge(f, mesh);
+      const halfedge_descriptor he = halfedge(*fitr, mesh);
       Point_3 pt = CGAL::centroid(
         point_pmap[source(he, mesh)],
         point_pmap[target(he, mesh)],
         point_pmap[target(next(he, mesh), mesh)]);
       Vector_3 vec = vector_functor(CGAL::ORIGIN, pt);
-      FT area = area_pmap[f];
+      FT area = area_pmap[*fitr];
       centroid = sum_functor(centroid, scale_functor(vec, area));
       sum_area += area;
     }
     centroid = scale_functor(centroid, FT(1) / sum_area);
 
-    return Plane_3(CGAL::ORIGIN + centroid, normal);
+    return Plane_3(CGAL::ORIGIN + centroid, norm);
   }
 
   const TriangleMesh &mesh;
   const FacetAreaMap area_pmap;
+  const FacetNormalMap normal_pmap;
   const VertexPointMap point_pmap;
+  Construct_vector_3 vector_functor;
+  Construct_scaled_vector_3 scale_functor;
+  Construct_sum_of_vectors_3 sum_functor;
+
 };
 
 template<typename PlaneProxy,
@@ -288,8 +304,11 @@ template<typename PlaneProxy,
 
 // Bundled l21 approximation traits
 template<typename PlaneProxy,
+  typename TriangleMesh,
   typename L21ErrorMetric,
   typename L21ProxyFitting,
+  typename L21PlaneFitting,
+  typename VertexPointMap,
   typename FacetNormalMap,
   typename FacetAreaMap>
   struct L21ApproximationTrait
@@ -298,11 +317,16 @@ template<typename PlaneProxy,
   typedef PlaneProxy Proxy;
   typedef L21ErrorMetric ErrorMetric;
   typedef L21ProxyFitting ProxyFitting;
+  typedef L21PlaneFitting PlaneFitting;
 
   L21ApproximationTrait(
+    const TriangleMesh &_mesh,
+    const VertexPointMap &_point_pmap,
     const FacetNormalMap &_facet_normal_map,
     const FacetAreaMap &_facet_area_map)
-    : normal_pmap(_facet_normal_map),
+    : mesh(_mesh),
+    point_pmap(_point_pmap),
+    normal_pmap(_facet_normal_map),
     area_pmap(_facet_area_map) {}
 
   // traits function object form
@@ -316,7 +340,14 @@ template<typename PlaneProxy,
     return ProxyFitting(normal_pmap, area_pmap);
   }
 
+  // construct plane fitting functor
+  L21PlaneFitting construct_plane_fitting_functor() const {
+    return L21PlaneFitting(mesh, area_pmap, normal_pmap, point_pmap);
+  }
+
 private:
+  const TriangleMesh &mesh;
+  const VertexPointMap point_pmap;
   const FacetNormalMap normal_pmap;
   const FacetAreaMap area_pmap;
 };