add some functionality to prevent flipping triangles

Algorithms/BaseRemesherT.hh

@@ -86,10 +86,11 @@ public:
     FACE_SELECTION
   };
 
-  typedef typename Mesh::Scalar        Scalar;
-  typedef typename Mesh::Point         Point;
-  typedef typename Mesh::EdgeHandle    EdgeHandle;
-  typedef typename Mesh::VertexHandle  VertexHandle;
+  typedef typename Mesh::Scalar         Scalar;
+  typedef typename Mesh::Point          Point;
+  typedef typename Mesh::EdgeHandle     EdgeHandle;
+  typedef typename Mesh::HalfedgeHandle HalfedgeHandle;
+  typedef typename Mesh::VertexHandle   VertexHandle;
 
 
   BaseRemesherT(Mesh& _mesh, ProgressEmitter* _progress = NULL);
@@ -125,6 +126,8 @@ protected:
   virtual bool is_too_long  (VertexHandle _v0, VertexHandle _v1) const = 0;
   virtual bool is_too_short (VertexHandle _v0, VertexHandle _v1) const = 0;
 
+  bool edge_flip_flips_normal(EdgeHandle _eh);
+  bool collapse_flips_normal(HalfedgeHandle _heh);
 
   
 protected:

Algorithms/BaseRemesherT_impl.hh

@@ -723,6 +723,12 @@ collapse_short_edges()
               hcol01 = false;
           }
 
+          // check if collapse flips triangles
+          if (collapse_flips_normal(h01))
+            hcol01 = false;
+          if (collapse_flips_normal(h10))
+            hcol10 = false;
+
           // try v1 -> v0
           if (hcol10) {
             // don't create too long edges
@@ -850,7 +856,7 @@ flip_edges()
 
           ve_after = ve0 + ve1 + ve2 + ve3;
 
-          if (ve_before > ve_after && mesh_.is_flip_ok(*e_it)) {
+          if (ve_before > ve_after && mesh_.is_flip_ok(*e_it) && !edge_flip_flips_normal(*e_it)) {
             mesh_.flip(*e_it);
             --mesh_.property(valences_, v0);
             --mesh_.property(valences_, v1);
@@ -888,6 +894,9 @@ tangential_smoothing(bool _use_projection)
 				   !mesh_.status(*v_it).feature() &&
 				   !mesh_.is_boundary(*v_it) );
 
+  mesh_.update_face_normals();
+
+  double damping = 1.0;
 
   // smooth
   for (int iters=0; iters<10; ++iters)
@@ -896,31 +905,51 @@ tangential_smoothing(bool _use_projection)
     {
       if (mesh_.status(*v_it).tagged())
       {
-	u.vectorize(0.0);
-	valence = 0;
-
-	for (vv_it=mesh_.cvv_iter(*v_it); vv_it.is_valid(); ++vv_it)
-	{
-	  u += mesh_.point(*vv_it);
-	  ++valence;
-	}
-
-	if (valence)
-	{
-	  u *= (1.0/valence);
-	  u -= mesh_.point(*v_it);
-	  n  = mesh_.normal(*v_it);
-	  n *= (u | n);
-	  u -= n;
-	}
-
-	mesh_.property(update_, *v_it) = u;
+        u.vectorize(0.0);
+        valence = 0;
+
+        for (vv_it=mesh_.cvv_iter(*v_it); vv_it.is_valid(); ++vv_it)
+        {
+          u += mesh_.point(*vv_it);
+          ++valence;
+        }
+
+        if (valence)
+        {
+          u *= (1.0/valence);
+          u -= mesh_.point(*v_it);
+          n  = mesh_.normal(*v_it);
+          n *= (u | n);
+          u -= n;
+        }
+
+        mesh_.property(update_, *v_it) = u;
       }
     }
 
     for (v_it=mesh_.vertices_begin(); v_it!=v_end; ++v_it)
       if (mesh_.status(*v_it).tagged())
-	mesh_.point(*v_it) += mesh_.property(update_, *v_it);
+        mesh_.point(*v_it) += damping*mesh_.property(update_, *v_it);
+
+    // check if normals changed
+    bool normals_changed = false;
+    for (auto fh : mesh_.faces())
+    {
+      if ((mesh_.calc_face_normal(fh) | mesh_.normal(fh)) < 0)
+      {
+        normals_changed = true;
+        break;
+      }
+    }
+
+    if (normals_changed)
+    {
+      // revert update and try again with smaller step
+      for (v_it=mesh_.vertices_begin(); v_it!=v_end; ++v_it)
+        if (mesh_.status(*v_it).tagged())
+          mesh_.point(*v_it) -= damping*mesh_.property(update_, *v_it);
+      damping *= 0.5;
+    }
   }
 
 
@@ -1124,6 +1153,83 @@ remove_caps()
 }
 
 
+template<class Mesh>
+bool
+BaseRemesherT<Mesh>::
+edge_flip_flips_normal(EdgeHandle _eh)
+{
+  if (mesh_.is_boundary(_eh))
+    return true;
+
+  auto heh = mesh_.halfedge_handle(_eh, 0);
+
+  // get the four points of the two incident faces in ccw order
+  auto p0 = mesh_.point(mesh_.to_vertex_handle(heh));
+  auto p1 = mesh_.point(mesh_.to_vertex_handle(mesh_.next_halfedge_handle(heh)));
+  auto p2 = mesh_.point(mesh_.from_vertex_handle(heh));
+  auto p3 = mesh_.point(mesh_.to_vertex_handle(mesh_.next_halfedge_handle(mesh_.opposite_halfedge_handle(heh))));
+
+  // compute normals before flip
+  auto n0_before = (p1-p0) % (p2-p0);
+  auto n1_before = (p2-p0) % (p3-p0);
+
+  // compute normals after flip
+  auto n0_after = (p1-p0) % (p3-p0);
+  auto n1_after = (p3-p2) % (p1-p2);
+
+  // compare all pairs of before and after normals
+  if ((n0_before | n0_after) < 0 || (n1_before | n1_after) < 0 ||
+      (n0_before | n1_after) < 0 || (n1_before | n0_after) < 0)
+    return true;
+
+  return false;
+}
+
+
+template<class Mesh>
+bool
+BaseRemesherT<Mesh>::
+collapse_flips_normal(HalfedgeHandle _heh)
+{
+  if (!mesh_.is_collapse_ok(_heh))
+    return true;
+
+  // get faces that are removed by the collapse
+  auto fh0 = mesh_.face_handle(_heh);
+  auto fh1 = mesh_.face_handle(mesh_.opposite_halfedge_handle(_heh));
+
+  auto collapsing_vertex = mesh_.from_vertex_handle(_heh);
+  auto point_before = mesh_.point(collapsing_vertex);
+
+  // compute normals before collapse
+  std::vector<decltype (mesh_.calc_face_normal(OpenMesh::FaceHandle(0)))> normals_before;
+  for (auto fh : mesh_.vf_range(collapsing_vertex))
+    normals_before.push_back(mesh_.calc_face_normal(fh));
+
+  // move point
+  mesh_.point(collapsing_vertex) = mesh_.point(mesh_.to_vertex_handle(_heh));
+
+  bool collapse_ok = true;
+
+  int i = 0;
+  for (auto fh : mesh_.vf_range(collapsing_vertex))
+  {
+    if (fh != fh0 && fh != fh1) // we don't care about faces that are removec by the collapse
+    {
+      auto normal_after = mesh_.calc_face_normal(fh);
+      if ((normal_after | normals_before[i]) <= 0)
+        collapse_ok = false;
+    }
+    ++i;
+  }
+
+  // move point back
+  mesh_.point(collapsing_vertex) = point_before;
+
+  return !collapse_ok;
+}
+
+
 //=============================================================================
 } // namespace Remeshing
 //=============================================================================