diff --git a/Surface_mesh_approximation/examples/Surface_mesh_approximation/vsa_class_interface_example.cpp b/Surface_mesh_approximation/examples/Surface_mesh_approximation/vsa_class_interface_example.cpp
index e5f9e3281c3..a4312888bf3 100644
--- a/Surface_mesh_approximation/examples/Surface_mesh_approximation/vsa_class_interface_example.cpp
+++ b/Surface_mesh_approximation/examples/Surface_mesh_approximation/vsa_class_interface_example.cpp
@@ -42,7 +42,7 @@ int main()
     l21_approx.run_one_step();
 
   // add proxies to the one with the maximum fitting error
-  l21_approx.add_proxies_furthest(3);
+  l21_approx.add_proxies_furthest(3, 5);
   for (std::size_t i = 0; i < 10; ++i)
     l21_approx.run_one_step();
 
diff --git a/Surface_mesh_approximation/include/CGAL/VSA_approximation.h b/Surface_mesh_approximation/include/CGAL/VSA_approximation.h
index 515a56905da..67c14e0d0ab 100644
--- a/Surface_mesh_approximation/include/CGAL/VSA_approximation.h
+++ b/Surface_mesh_approximation/include/CGAL/VSA_approximation.h
@@ -454,28 +454,24 @@ public:
   }
 
   /*!
-   * @brief Add proxies to the regions with the maximum fitting error.
-   * Except for the first one, a coarse re-fitting is performed before each proxy is inserted.
+   * @brief Add proxies to the worst regions one by one.
+   * The re-fitting is performed after each proxy is inserted.
    * @pre current facet proxy map is valid
-   * @note after the addition, the facet proxy map is invalid
+   * @note after the addition, the facet proxy map remains valid
    * @param num_proxies number of proxies to be added
-   * @param num_iterations the number of iterations of coarse re-fitting
-   * @return number of proxies successfully added
+   * @param num_iterations the number of iterations of re-fitting
+   * @return number of proxies added
    */
   std::size_t add_proxies_furthest(const std::size_t num_proxies,
-    const std::size_t num_iterations = 5) {
-    // when insert only one proxy, it has the same effect of add_proxy_furthest()
-    if (num_proxies == 0 || !add_proxy_furthest())
-      return 0;
-
-    std::size_t num_added = 1;
+    const std::size_t num_iterations) {
+    std::size_t num_added = 0;
     for (; num_added < num_proxies; ++num_added) {
+      if (!add_proxy_furthest())
+        break;
       for (std::size_t i = 0; i < num_iterations; ++i) {
         partition();
         fit();
       }
-      if (!add_proxy_furthest())
-        return num_added;
     }
     return num_added;
   }
@@ -1083,44 +1079,46 @@ private:
   }
 
   /*!
-   * @brief Inserts a proxy at the furthest facet of the region with the maximum fitting error.
-   * No re-fitting is performed.
-   * @pre current facet proxy map is valid
-   * @note after the addition, the facet proxy map is invalid
-   * @return true if insertion success, false otherwise
+   * @brief Add a proxy seed at the facet with the maximum fitting error.
+   * @pre current facet proxy map is valid, proxy error is computed
+   * @note No re-fitting is performed. After the operation, the facet proxy map remains valid.
+   * @return true add successfully, false otherwise
    */
   bool add_proxy_furthest() {
-    std::vector<FT> px_error(proxies.size(), FT(0.0));
-    std::vector<FT> max_facet_error(proxies.size(), FT(0.0));
-    std::vector<face_descriptor> max_facet(proxies.size());
+#ifdef CGAL_SURFACE_MESH_APPROXIMATION_DEBUG
+    std::cerr << "add furthest " << proxies.size() << std::endl;
+#endif
+    compute_fitting_error();
+    FT max_error = proxies.front().err;
+    std::size_t px_worst = 0;
+    for (std::size_t i = 0; i < proxies.size(); ++i) {
+      if (max_error < proxies[i].err) {
+        max_error = proxies[i].err;
+        px_worst = i;
+      }
+    }
 
+    face_descriptor fworst;
+    bool first = true;
     BOOST_FOREACH(face_descriptor f, faces(*m_pmesh)) {
       std::size_t px_idx = fproxy_map[f];
+      if (px_idx != px_worst || f == proxies[px_idx].seed)
+        continue;
+
       FT err = (*fit_error)(f, proxies[px_idx].px);
-      px_error[px_idx] += err;
-
-      if (err > max_facet_error[px_idx]) {
-        max_facet_error[px_idx] = err;
-        max_facet[px_idx] = f;
+      if (first || max_error < err) {
+        first = false;
+        max_error = err;
+        fworst = f;
       }
     }
 
-    FT max_px_error = px_error.front();
-    std::size_t max_px_idx = 0;
-    for (std::size_t i = 0; i < proxies.size(); ++i) {
-      if (px_error[i] > max_px_error) {
-        max_px_error = px_error[i];
-        max_px_idx = i;
-      }
-    }
-    if (max_facet[max_px_idx] == proxies[max_px_idx].seed) {
-#ifdef CGAL_SURFACE_MESH_APPROXIMATION_DEBUG
-      std::cerr << "add furthest failed" << std::endl;
-#endif
+    if (first)
       return false;
-    }
 
-    proxies.push_back(fit_new_proxy(max_facet[max_px_idx]));
+    fproxy_map[fworst] = proxies.size();
+    proxies.push_back(fit_new_proxy(fworst));
+
     return true;
   }
 
diff --git a/Surface_mesh_approximation/test/Surface_mesh_approximation/vsa_class_interface_test.cpp b/Surface_mesh_approximation/test/Surface_mesh_approximation/vsa_class_interface_test.cpp
index 8814f440bd1..34729596ab9 100644
--- a/Surface_mesh_approximation/test/Surface_mesh_approximation/vsa_class_interface_test.cpp
+++ b/Surface_mesh_approximation/test/Surface_mesh_approximation/vsa_class_interface_test.cpp
@@ -61,7 +61,7 @@ int main()
 
   // incremental add and run until convergence
   std::cout << "incremental add and run until convergence" << std::endl;
-  l2_approx.add_proxies_furthest(3);
+  l2_approx.add_proxies_furthest(3, 5);
   if (l2_approx.run_until_convergence(0.1))
     std::cout << "Converged." << std::endl;
   if (l2_approx.get_proxies_size() != 13)