diff --git a/Algorithms/BaseRemesherT.hh b/Algorithms/BaseRemesherT.hh index f75fa3a6721316ec63aed5ec35536ec821bc237b..6ce4e945108ce553148778c44d544130e4f00a60 100644 --- a/Algorithms/BaseRemesherT.hh +++ b/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: diff --git a/Algorithms/BaseRemesherT_impl.hh b/Algorithms/BaseRemesherT_impl.hh index 5a54a14593f775ebd2eafe3292150f8dfc5a0726..176ed624b5edf71f5809498cc164f7e7f6506724 100644 --- a/Algorithms/BaseRemesherT_impl.hh +++ b/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); @@ -723,6 +738,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 +871,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); @@ -878,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; @@ -888,6 +908,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 +919,79 @@ 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); + int valence = 0; + int feature_valence = 0; + + for (auto heh : mesh_.voh_range(*v_it)) + { + u += mesh_.point(mesh_.to_vertex_handle(heh)); + if (mesh_.status(mesh_.edge_handle(heh)).feature()) + ++feature_valence; + ++valence; + } + + if (valence) + { + u *= (1.0/valence); + u -= mesh_.point(*v_it); + n = mesh_.normal(*v_it); + n *= (u | n); + 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; } } 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 +1195,83 @@ remove_caps() } +template +bool +BaseRemesherT:: +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 +bool +BaseRemesherT:: +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 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 //=============================================================================