better feature handling

Algorithms/BaseRemesherT_impl.hh

@@ -593,7 +593,6 @@ split_long_edges()
           eh = (is_boundary ? mesh_.edge_handle(mesh_.n_edges() - 2) : mesh_.edge_handle(mesh_.n_edges() - 3));
 
           mesh_.status(eh).set_feature(true);
-          mesh_.status(vh).set_feature(true);
         } else {
           project_to_reference(vh);
         }
@@ -658,9 +657,6 @@ collapse_short_edges()
           f1 = mesh_.status(v1).feature();
           hcol01 = hcol10 = true;
 
-          if (mesh_.status(*e_it).feature() && !f0 && !f1)
-            std::cerr << "Bad luck" << std::endl;
-
           // boundary rules
           if (b0 && b1) {
             if (!mesh_.is_boundary(*e_it))
@@ -678,7 +674,45 @@ collapse_short_edges()
           else if (l1)
             hcol10 = false;
 
-          // feature rules
+          // feature vertex rules
+          if (f0 && f1)
+            continue;
+          else if (f0)
+            hcol01 = false;
+          else if (f1)
+            hcol10 = false;
+
+          // feature edge rules
+          // if vertex is incident to feature edge but collapse edge is not feature dont collapse
+
+          // hcol01
+          {
+            int feature_valence_0 = 0;
+            for (auto eh : mesh_.ve_range(v0))
+              if (mesh_.status(eh).feature())
+                ++feature_valence_0;
+            if (feature_valence_0 == 2)
+            {
+              if (!mesh_.status(*e_it).feature())
+                hcol01 = false;
+            }
+            else if (feature_valence_0 != 0)
+              hcol01 = false;
+          }
+          // hcol10
+          {
+            int feature_valence_1 = 0;
+            for (auto eh : mesh_.ve_range(v1))
+              if (mesh_.status(eh).feature())
+                ++feature_valence_1;
+            if (feature_valence_1 == 2)
+            {
+              if (!mesh_.status(*e_it).feature())
+                hcol10 = false;
+            }
+            else if (feature_valence_1 != 0)
+              hcol10 = false;
+          }
 
           // the other two edges removed by collapse must not be features
           h0 = mesh_.prev_halfedge_handle(h01);
@@ -690,25 +724,6 @@ collapse_short_edges()
           if (mesh_.status(mesh_.edge_handle(h0)).feature() || mesh_.status(mesh_.edge_handle(h1)).feature())
             hcol10 = false;
 
-          if (f0 && f1) {
-            // edge must be feature
-            if (!mesh_.status(*e_it).feature())
-              continue;
-
-            // the other two edges removed by collapse must not be features
-            h0 = mesh_.prev_halfedge_handle(h01);
-            h1 = mesh_.next_halfedge_handle(h10);
-            if (mesh_.status(mesh_.edge_handle(h0)).feature() || mesh_.status(mesh_.edge_handle(h1)).feature())
-              hcol01 = false;
-            h0 = mesh_.prev_halfedge_handle(h10);
-            h1 = mesh_.next_halfedge_handle(h01);
-            if (mesh_.status(mesh_.edge_handle(h0)).feature() || mesh_.status(mesh_.edge_handle(h1)).feature())
-              hcol10 = false;
-          } else if (f0)
-            hcol01 = false;
-          else if (f1)
-            hcol10 = false;
-
           // topological rules
           if (hcol01)
             hcol01 = mesh_.is_collapse_ok(h01);
@@ -884,7 +899,6 @@ tangential_smoothing(bool _use_projection)
 {
   typename Mesh::VIter     v_it, v_end(mesh_.vertices_end());
   typename Mesh::CVVIter   vv_it;
-  typename Mesh::Scalar    valence;
   typename Mesh::Point     u, n;
 
 
@@ -906,11 +920,14 @@ tangential_smoothing(bool _use_projection)
       if (mesh_.status(*v_it).tagged())
       {
         u.vectorize(0.0);
-        valence = 0;
+        int valence = 0;
+        int feature_valence = 0;
 
-        for (vv_it=mesh_.cvv_iter(*v_it); vv_it.is_valid(); ++vv_it)
+        for (auto heh : mesh_.voh_range(*v_it))
         {
-          u += mesh_.point(*vv_it);
+          u += mesh_.point(mesh_.to_vertex_handle(heh));
+          if (mesh_.status(mesh_.edge_handle(heh)).feature())
+            ++feature_valence;
           ++valence;
         }
 
@@ -923,6 +940,31 @@ tangential_smoothing(bool _use_projection)
           u -= n;
         }
 
+        if (feature_valence == 2)
+        {
+          // project update onto feature directions
+          typename Mesh::Point feature_dir(0.0,0.0,0.0);
+          for (auto heh : mesh_.voh_range(*v_it))
+          {
+            if (mesh_.status(mesh_.edge_handle(heh)).feature())
+            {
+              auto dir = mesh_.point(mesh_.to_vertex_handle(heh)) - mesh_.point(mesh_.from_vertex_handle(heh));
+              if ((dir | feature_dir) > 0)
+                feature_dir += dir;
+              else
+                feature_dir -= dir;
+            }
+          }
+          if (feature_dir.sqrnorm() > 0)
+            feature_dir.normalize();
+          u = (feature_dir | u) * feature_dir;
+        }
+        else if (feature_valence != 0)
+        {
+          // more than two or exactly one incident feature edge means vertex should be preserved
+          u.vectorize(0.0);
+        }
+
         mesh_.property(update_, *v_it) = u;
       }
     }