fix error rounding up

code and document

Surface_mesh_approximation/doc/Surface_mesh_approximation/surface_mesh_approximation.txt

@@ -60,10 +60,10 @@ Assuming a clustering partition of \f$n\f$ regions with errors \f$ \{E_k\}_{k=1\
 - <b>Hierarchical</b>. \f$m\f$ seed facets are dispatched within the current partition, where each partition is refined with a number of proxies chosen in accordance to their fitting error:
   - calculate total error \f$ E_{total} \f$, then average error \f$ E_{avg} = E_{total} / m \f$ (assuming that each new proxy shares the same amount of error)
   - sort errors \f$ \{E_{min},\cdots,E_{max}\} \f$
-  - from \f$ E_{min} \f$ to \f$ E_{max} \f$, we diffuse the error <em>hierarchically</em> one after another. More specifically, the number of proxies added to the \f$k\f$th region is proportional to its error given by the formula:
+  - from \f$ E_{min} \f$ to \f$ E_{max} \f$, we diffuse the error <em>hierarchically</em> one after another. More specifically, the number of proxies \f$N_k\f$ added to the \f$k\f$th region is proportional to its error:
-\f[ nb\_to\_add\_k = rounded\_to\_nearest\_integer(E_k / E_{avg}), \f]
+\f[ N_k = \lfloor E_k / E_{avg} + 0.5 \rfloor, \f]
 and the remaining error is added to the next proxy error in order to keep the total error unchanged:
-\f[ E'_{k+1} = (E_k - nb\_to\_add\_k * E_{avg}) + E_{k+1}. \f]
+\f[ E'_{k+1} = (E_k - N_k * E_{avg}) + E_{k+1}. \f]
 
 \cgalFigureBegin{seeding_method, seeding_method.png}
 Comparison of seeding methods on the sphere-cube model. From left to right: initial partition (\f$ \mathcal{L}^{2,1} \f$ metrics and 20 proxies), add 5 proxy seeds (red facets) with random, incremental and hierarchical methods respectively.

Surface_mesh_approximation/include/CGAL/vsa_approximation.h

@@ -545,59 +545,54 @@ public:
     std::cerr << "#px " << m_proxies.size() << std::endl;
 #endif
 
-    const FT sum_error = compute_total_error();
+    const double sum_error = CGAL::to_double(compute_total_error());
-    const FT avg_error = sum_error / FT(static_cast<double>(num_proxies));
+    const double avg_error = sum_error / static_cast<double>(num_proxies);
-
-    std::vector<Proxy_error> px_error;
-    for (std::size_t i = 0; i < m_proxies.size(); ++i)
-      px_error.push_back(Proxy_error(i, m_proxies[i].err));
-    // sort partition by error
-    std::sort(px_error.begin(), px_error.end());
 
     // number of proxies to be added to each region
     std::vector<std::size_t> num_to_add(m_proxies.size(), 0);
-    if (avg_error == FT(0.0)) {
+    if (avg_error <= 0.0) {
       // rare case on extremely regular geometry like a cube
 #ifdef CGAL_SURFACE_MESH_APPROXIMATION_DEBUG
       std::cerr << "zero error, diffuse w.r.t. number of facets" << std::endl;
 #endif
-      const FT avg_facet = FT(
+      const double avg_facet =
-        static_cast<double>(num_faces(*m_ptm)) / static_cast<double>(num_proxies));
+        static_cast<double>(num_faces(*m_ptm)) / static_cast<double>(num_proxies);
-      std::vector<FT> px_size(m_proxies.size(), FT(0.0));
+      std::vector<double> px_size(m_proxies.size(), 0.0);
       BOOST_FOREACH(face_descriptor f, faces(*m_ptm))
-        px_size[get(m_fproxy_map, f)] += FT(1.0);
+        px_size[get(m_fproxy_map, f)] += 1.0;
-      FT residual(0.0);
+      double residual = 0.0;
       for (std::size_t i = 0; i < m_proxies.size(); ++i) {
-        FT to_add = (residual + px_size[i]) / avg_facet;
+        const double to_add = (residual + px_size[i]) / avg_facet;
-        FT floor_to_add = FT(std::floor(CGAL::to_double(to_add)));
+        const double to_add_round_up = std::floor(to_add + 0.5);
-        FT q_to_add = FT(CGAL::to_double(
+        residual = (to_add - to_add_round_up) * avg_facet;
-          ((to_add - floor_to_add) > FT(0.5)) ? (floor_to_add + FT(1.0)) : floor_to_add));
+        num_to_add[i] = static_cast<std::size_t>(to_add_round_up);
-        residual = (to_add - q_to_add) * avg_facet;
-        num_to_add[i] = static_cast<std::size_t>(CGAL::to_double(q_to_add));
       }
     }
     else {
+      std::vector<Proxy_error> px_error;
+      for (std::size_t i = 0; i < m_proxies.size(); ++i)
+        px_error.push_back(Proxy_error(i, m_proxies[i].err));
+      // sort partition by error
+      std::sort(px_error.begin(), px_error.end());
+
       // residual from previous proxy in range (-0.5, 0.5] * avg_error
-      FT residual(0.0);
+      double residual = 0.0;
-      BOOST_FOREACH(const Proxy_error &pxe, px_error) {
+      for (std::size_t i = 0; i < m_proxies.size(); ++i) {
         // add error residual from previous proxy
         // to_add maybe negative but greater than -0.5
-        FT to_add = (residual + pxe.err) / avg_error;
+        const double to_add = (residual + CGAL::to_double(px_error[i].err)) / avg_error;
-        // floor_to_add maybe negative but no less than -1
+        const double to_add_round_up = std::floor(to_add + 0.5);
-        FT floor_to_add = FT(std::floor(CGAL::to_double(to_add)));
+        residual = (to_add - to_add_round_up) * avg_error;
-        FT q_to_add = FT(CGAL::to_double(
+        num_to_add[i] = static_cast<std::size_t>(to_add_round_up);
-          ((to_add - floor_to_add) > FT(0.5)) ? (floor_to_add + FT(1.0)) : floor_to_add));
-        residual = (to_add - q_to_add) * avg_error;
-        num_to_add[pxe.px] = static_cast<std::size_t>(CGAL::to_double(q_to_add));
       }
-    }
 
 #ifdef CGAL_SURFACE_MESH_APPROXIMATION_DEBUG
-    for (std::size_t i = 0; i < px_error.size(); ++i)
+      for (std::size_t i = 0; i < px_error.size(); ++i)
-      std::cerr << "#px " << px_error[i].px
+        std::cerr << "#px " << px_error[i].px
-        << ", #error " << px_error[i].err
+          << ", #error " << px_error[i].err
-        << ", #num_to_add " << num_to_add[px_error[i].px] << std::endl;
+          << ", #num_to_add " << num_to_add[px_error[i].px] << std::endl;
 #endif
+    }
 
     std::size_t num_added = 0;
     BOOST_FOREACH(face_descriptor f, faces(*m_ptm)) {