From 2df0fc716d4366a76b974f67638389d14d38f888 Mon Sep 17 00:00:00 2001 From: Lingjie Zhu Date: Mon, 25 Feb 2019 11:34:14 +0800 Subject: [PATCH] handle degenerate faces with zero area use boost::math::isnormal and CGAL::collinear predicates --- .../L21_metric_plane_proxy.h | 13 ++++-- .../CGAL/Variational_shape_approximation.h | 40 +++++++++++++------ 2 files changed, 38 insertions(+), 15 deletions(-) diff --git a/Surface_mesh_approximation/include/CGAL/Surface_mesh_approximation/L21_metric_plane_proxy.h b/Surface_mesh_approximation/include/CGAL/Surface_mesh_approximation/L21_metric_plane_proxy.h index 71a446118bf..43bc45ed46d 100644 --- a/Surface_mesh_approximation/include/CGAL/Surface_mesh_approximation/L21_metric_plane_proxy.h +++ b/Surface_mesh_approximation/include/CGAL/Surface_mesh_approximation/L21_metric_plane_proxy.h @@ -29,6 +29,7 @@ #include #include +#include #include #include @@ -97,7 +98,12 @@ public: const Point_3 &p0 = vpmap[source(he, tm)]; const Point_3 &p1 = vpmap[target(he, tm)]; const Point_3 &p2 = vpmap[target(next(he, tm), tm)]; - put(m_fnmap, f, CGAL::unit_normal(p0, p1, p2)); + if (CGAL::collinear(p0, p1, p2)) { + std::cout << CGAL::unit_normal(p0, p1, p2) << std::endl; + put(m_fnmap, f, CGAL::NULL_VECTOR); + } + else + put(m_fnmap, f, CGAL::unit_normal(p0, p1, p2)); put(m_famap, f, std::sqrt(CGAL::to_double(CGAL::squared_area(p0, p1, p2)))); } } @@ -134,8 +140,9 @@ public: norm = m_sum_functor(norm, m_scale_functor(get(m_fnmap, f), get(m_famap, f))); } - norm = m_scale_functor(norm, - FT(1.0 / std::sqrt(CGAL::to_double(norm.squared_length())))); + const FT inv_len = FT(1.0 / std::sqrt(CGAL::to_double(norm.squared_length()))); + if ((boost::math::isnormal)(inv_len)) + norm = m_scale_functor(norm, inv_len); return norm; } diff --git a/Surface_mesh_approximation/include/CGAL/Variational_shape_approximation.h b/Surface_mesh_approximation/include/CGAL/Variational_shape_approximation.h index c4815d69ad4..7228c83c52f 100644 --- a/Surface_mesh_approximation/include/CGAL/Variational_shape_approximation.h +++ b/Surface_mesh_approximation/include/CGAL/Variational_shape_approximation.h @@ -38,6 +38,7 @@ #include +#include #include #include #include @@ -240,7 +241,6 @@ private: const Error_metric *m_metric; Construct_vector_3 vector_functor; - Construct_point_3 point_functor; Construct_scaled_vector_3 scale_functor; Construct_sum_of_vectors_3 sum_functor; Compute_scalar_product_3 scalar_product_functor; @@ -290,7 +290,6 @@ public: Geom_traits traits; vector_functor = traits.construct_vector_3_object(); - point_functor = traits.construct_point_3_object(); scale_functor = traits.construct_scaled_vector_3_object(); sum_functor = traits.construct_sum_of_vectors_3_object(); scalar_product_functor = traits.compute_scalar_product_3_object(); @@ -1430,12 +1429,19 @@ private: const Point_3 &p1 = m_vpoint_map[target(he, *m_ptm)]; const Point_3 &p2 = m_vpoint_map[target(next(he, *m_ptm), *m_ptm)]; const FT farea(std::sqrt(CGAL::to_double(CGAL::squared_area(p0, p1, p2)))); - const Vector_3 fnorm = CGAL::unit_normal(p0, p1, p2); + const Vector_3 fnorm = + CGAL::collinear(p0, p1, p2) ? CGAL::NULL_VECTOR : CGAL::unit_normal(p0, p1, p2); norm = sum_functor(norm, scale_functor(fnorm, farea)); area += farea; } - norm = scale_functor(norm, FT(1.0 / std::sqrt(CGAL::to_double(norm.squared_length())))); + const FT inv_len = FT(1.0 / std::sqrt(CGAL::to_double(norm.squared_length()))); + if ((boost::math::isnormal)(inv_len)) + norm = scale_functor(norm, inv_len); + else + norm = Vector_3(FT(0.0), FT(0.0), FT(1.0)); + + std::cout << inv_len << " # " << norm << " # " << area << std::endl; m_px_planes.push_back(Proxy_plane(fit_plane, norm, area)); } @@ -1910,7 +1916,7 @@ private: // projection FT sum_area(0.0); - Vector_3 vec = vector_functor(CGAL::ORIGIN, point_functor(CGAL::ORIGIN)); + Vector_3 vec = CGAL::NULL_VECTOR; const Point_3 vtx_pt = m_vpoint_map[v]; BOOST_FOREACH(const std::size_t px_idx, px_set) { const Vector_3 proj = vector_functor( @@ -1919,9 +1925,13 @@ private: vec = sum_functor(vec, scale_functor(proj, area)); sum_area += area; } - vec = scale_functor(vec, FT(1.0) / sum_area); - - a.pos = translate_point_functor(CGAL::ORIGIN, vec); + const FT inv_sum_area = FT(1.0) / sum_area; + if ((boost::math::isnormal)(inv_sum_area)) + a.pos = translate_point_functor( + CGAL::ORIGIN, + scale_functor(vec, inv_sum_area)); + else + a.pos = vtx_pt; } } @@ -1965,15 +1975,21 @@ private: Vector_3 vec = vector_functor(CGAL::ORIGIN, CGAL::centroid(p0, p1, p2)); FT farea(std::sqrt(CGAL::to_double(CGAL::squared_area(p0, p1, p2)))); - Vector_3 fnorm = CGAL::unit_normal(p0, p1, p2); + Vector_3 fnorm = + CGAL::collinear(p0, p1, p2) ? CGAL::NULL_VECTOR : CGAL::unit_normal(p0, p1, p2); norm = sum_functor(norm, scale_functor(fnorm, farea)); cent = sum_functor(cent, scale_functor(vec, farea)); sum_area += farea; } - norm = scale_functor(norm, - FT(1.0 / std::sqrt(CGAL::to_double(norm.squared_length())))); - cent = scale_functor(cent, FT(1.0) / sum_area); + const FT inv_len = FT(1.0 / std::sqrt(CGAL::to_double(norm.squared_length()))); + if ((boost::math::isnormal)(inv_len)) + norm = scale_functor(norm, inv_len); + else + norm = Vector_3(FT(0.0), FT(0.0), FT(1.0)); + const FT inv_sum_area = FT(1.0) / sum_area; + if ((boost::math::isnormal)(inv_sum_area)) + cent = scale_functor(cent, inv_sum_area); return Plane_3(CGAL::ORIGIN + cent, norm); }