From 4003fd8520453ab62fa82e485b3bb13cd5ea81a8 Mon Sep 17 00:00:00 2001 From: Laurent Saboret Date: Fri, 19 Jun 2009 16:00:59 +0000 Subject: [PATCH] Fixed infinite loop in Delaunay refinement with "teeth" model (BPoint_MM): use filtered robust version of circumcenter computation. (thanks to Stephane Tayeb for providing Robust_circumcenter_filtered_traits_3 class) --- .../Point_set_demo/CMakeLists.txt | 1 + .../poisson/CMakeLists.txt | 1 + .../CMakeLists.txt | 1 + .../include/CGAL/Mesh_3/Refine_tets.h | 6 + .../CGAL/Poisson_reconstruction_function.h | 7 +- .../Robust_circumcenter_filtered_traits_3.h | 227 ++++++++++++++++++ .../CGAL/poisson_refine_triangulation.h | 7 +- .../CMakeLists.txt | 1 + 8 files changed, 246 insertions(+), 5 deletions(-) create mode 100644 Surface_reconstruction_points_3/include/CGAL/Robust_circumcenter_filtered_traits_3.h diff --git a/Surface_reconstruction_points_3/demo/Surface_reconstruction_points_3/Point_set_demo/CMakeLists.txt b/Surface_reconstruction_points_3/demo/Surface_reconstruction_points_3/Point_set_demo/CMakeLists.txt index 39a1f85fd30..c375df831cc 100644 --- a/Surface_reconstruction_points_3/demo/Surface_reconstruction_points_3/Point_set_demo/CMakeLists.txt +++ b/Surface_reconstruction_points_3/demo/Surface_reconstruction_points_3/Point_set_demo/CMakeLists.txt @@ -115,6 +115,7 @@ if(CGAL_Qt4_FOUND AND QT4_FOUND AND OPENGL_FOUND AND QGLVIEWER_FOUND) # Temporary debugging stuff ADD_DEFINITIONS( "-DDEBUG_TRACE" ) + ADD_DEFINITIONS( "-DCGAL_PROFILE" ) qt4_wrap_ui( UI_FILES MainWindow.ui) diff --git a/Surface_reconstruction_points_3/demo/Surface_reconstruction_points_3/poisson/CMakeLists.txt b/Surface_reconstruction_points_3/demo/Surface_reconstruction_points_3/poisson/CMakeLists.txt index 578d395cc93..73609d4cef1 100644 --- a/Surface_reconstruction_points_3/demo/Surface_reconstruction_points_3/poisson/CMakeLists.txt +++ b/Surface_reconstruction_points_3/demo/Surface_reconstruction_points_3/poisson/CMakeLists.txt @@ -82,6 +82,7 @@ if(CGAL_FOUND AND MSVC AND OPENGL_FOUND AND TAUCS_FOUND) # Temporary debugging stuff ADD_DEFINITIONS( "-DDEBUG_TRACE" ) + ADD_DEFINITIONS( "-DCGAL_PROFILE" ) # Creates Poisson executable ADD_EXECUTABLE(Poisson WIN32 ChildFrm.cpp DialogOptions.cpp director.cpp MainFrm.cpp Poisson.cpp PoissonDoc.cpp PoissonView.cpp Poisson.rc) diff --git a/Surface_reconstruction_points_3/examples/Surface_reconstruction_points_3/CMakeLists.txt b/Surface_reconstruction_points_3/examples/Surface_reconstruction_points_3/CMakeLists.txt index 97801ef1b86..eb0d88ebba1 100644 --- a/Surface_reconstruction_points_3/examples/Surface_reconstruction_points_3/CMakeLists.txt +++ b/Surface_reconstruction_points_3/examples/Surface_reconstruction_points_3/CMakeLists.txt @@ -54,6 +54,7 @@ if ( CGAL_FOUND ) # Temporary debugging stuff ADD_DEFINITIONS( "-DDEBUG_TRACE" ) + ADD_DEFINITIONS( "-DCGAL_PROFILE" ) # Executables that do *not* require BLAS, LAPACK nor TAUCS create_single_source_cgal_program( "APSS_reconstruction.cpp" ) diff --git a/Surface_reconstruction_points_3/include/CGAL/Mesh_3/Refine_tets.h b/Surface_reconstruction_points_3/include/CGAL/Mesh_3/Refine_tets.h index c3289265f3c..698e13c47a5 100644 --- a/Surface_reconstruction_points_3/include/CGAL/Mesh_3/Refine_tets.h +++ b/Surface_reconstruction_points_3/include/CGAL/Mesh_3/Refine_tets.h @@ -195,6 +195,12 @@ public: { Zone zone; +#ifdef CGAL_MESHES_DEBUG_REFINEMENT_POINTS + // Check if triangulation's geometric traits provides a robust circumcenter computation + if (triangulation_ref_impl().side_of_sphere(c, p, true) != ON_BOUNDED_SIDE) + std::cerr << "Refine_tets_with_oracle_base::conflicts_zone_impl: ERROR: circumcenter out of sphere!\n"; +#endif + zone.cell = c; zone.locate_type = Tr::CELL; diff --git a/Surface_reconstruction_points_3/include/CGAL/Poisson_reconstruction_function.h b/Surface_reconstruction_points_3/include/CGAL/Poisson_reconstruction_function.h index 29a7d855312..370be9e5f80 100644 --- a/Surface_reconstruction_points_3/include/CGAL/Poisson_reconstruction_function.h +++ b/Surface_reconstruction_points_3/include/CGAL/Poisson_reconstruction_function.h @@ -33,6 +33,7 @@ #include #include #include +#include #include @@ -63,8 +64,10 @@ public: typedef Gt Geom_traits; ///< Geometric traits class - /// Typedef to Reconstruction_triangulation_3 - typedef Reconstruction_triangulation_3 Triangulation; + /// Internal 3D triangulation, of type Reconstruction_triangulation_3. + // Note: poisson_refine_triangulation() requires a robust circumcenter computation. + typedef Reconstruction_triangulation_3 > + Triangulation; // Geometric types typedef typename Geom_traits::FT FT; ///< typedef to Geom_traits::FT diff --git a/Surface_reconstruction_points_3/include/CGAL/Robust_circumcenter_filtered_traits_3.h b/Surface_reconstruction_points_3/include/CGAL/Robust_circumcenter_filtered_traits_3.h new file mode 100644 index 00000000000..f337c405a3b --- /dev/null +++ b/Surface_reconstruction_points_3/include/CGAL/Robust_circumcenter_filtered_traits_3.h @@ -0,0 +1,227 @@ +// Copyright (c) 2009 INRIA Sophia-Antipolis (France). +// All rights reserved. +// +// This file is part of CGAL (www.cgal.org); you may redistribute it under +// the terms of the Q Public License version 1.0. +// See the file LICENSE.QPL distributed with CGAL. +// +// Licensees holding a valid commercial license may use this file in +// accordance with the commercial license agreement provided with the software. +// +// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE +// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. +// +// $URL$ +// $Id$ +// +// +// Author(s) : Stéphane Tayeb +// +//****************************************************************************** +// File Description : +// +//****************************************************************************** + +#ifndef ROBUST_CIRCUMCENTER_FILTERED_TRAITS_3_H_ +#define ROBUST_CIRCUMCENTER_FILTERED_TRAITS_3_H_ + + +#include +#include +#include +#include +#include + + +namespace CGAL { + + +template < typename K > +class Robust_filtered_construct_circumcenter_3 +{ +public: + typedef typename K::Point_3 Point_3; + typedef typename K::FT FT; + typedef typename K::Sphere_3 Sphere_3; + typedef Point_3 result_type; + + typedef Exact_predicates_exact_constructions_kernel EK2; + typedef Regular_triangulation_euclidean_traits_3 EK; + typedef Cartesian_converter To_exact; + typedef Cartesian_converter Back_from_exact; + + + Point_3 operator() ( const Point_3 & p, + const Point_3 & q, + const Point_3 & r, + const Point_3 & s ) const + { + typename K::Construct_circumcenter_3 circumcenter = + K().construct_circumcenter_3_object(); + typename K::Has_on_bounded_side_3 on_bounded_side = + K().has_on_bounded_side_3_object(); + + // Compute denominator to swith to exact if it is 0. + // TODO: replace hard coded comparison with 1E-14 by static filter. + const FT denom = compute_denom(p,q,r,s); + if (denom < -1E-14 || denom > 1E-14) + { + result_type point = circumcenter(p,q,r,s); + + // Fast output + if ( on_bounded_side(Sphere_3(p,q,r,s),point) ) + return point; + } + + // Switch to exact + To_exact to_exact; + Back_from_exact back_from_exact; + EK::Construct_circumcenter_3 exact_circumcenter = + EK().construct_circumcenter_3_object(); + + return back_from_exact(exact_circumcenter(to_exact(p), + to_exact(q), + to_exact(r), + to_exact(s))); + } + + Point_3 operator() ( const Point_3 & p, + const Point_3 & q, + const Point_3 & r ) const + { + typename K::Construct_circumcenter_3 circumcenter = + K().construct_circumcenter_3_object(); + typename K::Has_on_bounded_side_3 on_bounded_side = + K().has_on_bounded_side_3_object(); + + // Compute denominator to swith to exact if it is 0. + // TODO: replace hard coded comparison with 1E-14 by static filter. + const FT denom = compute_denom(p,q,r); + if (denom < -1E-14 || denom > 1E-14) + { + result_type point = circumcenter(p,q,r); + + // Fast output + if ( on_bounded_side(Sphere_3(p,q,r),point) ) + return point; + } + + // Switch to exact + To_exact to_exact; + Back_from_exact back_from_exact; + EK::Construct_circumcenter_3 exact_circumcenter = + EK().construct_circumcenter_3_object(); + + return back_from_exact(exact_circumcenter(to_exact(p), + to_exact(q), + to_exact(r))); + } + + Point_3 operator() ( const Point_3 & p, + const Point_3 & q ) const + { + typename K::Construct_circumcenter_3 circumcenter = + K().construct_circumcenter_3_object(); + typename K::Has_on_bounded_side_3 on_bounded_side = + K().has_on_bounded_side_3_object(); + + // No division here + result_type point = circumcenter(p,q); + + // Fast output + if ( on_bounded_side(Sphere_3(p,q),point) ) + return point; + + // Switch to exact + To_exact to_exact; + Back_from_exact back_from_exact; + EK::Construct_circumcenter_3 exact_circumcenter = + EK().construct_circumcenter_3_object(); + + return back_from_exact(exact_circumcenter(to_exact(p), + to_exact(q))); + } + +private: + + FT compute_denom(const Point_3 & p, + const Point_3 & q, + const Point_3 & r, + const Point_3 & s) const + { + return compute_denom(p.x(),p.y(),p.z(), + q.x(),q.y(),q.z(), + r.x(),r.y(),r.z(), + s.x(),s.y(),s.z()); + } + + FT compute_denom(const Point_3 & p, + const Point_3 & q, + const Point_3 & r) const + { + return compute_denom(p.x(),p.y(),p.z(), + q.x(),q.y(),q.z(), + r.x(),r.y(),r.z()); + } + + FT compute_denom(const FT &px, const FT &py, const FT &pz, + const FT &qx, const FT &qy, const FT &qz, + const FT &rx, const FT &ry, const FT &rz, + const FT &sx, const FT &sy, const FT &sz) const + { + const FT qpx = qx-px; + const FT qpy = qy-py; + const FT qpz = qz-pz; + const FT rpx = rx-px; + const FT rpy = ry-py; + const FT rpz = rz-pz; + const FT spx = sx-px; + const FT spy = sy-py; + const FT spz = sz-pz; + + return determinant(qpx,qpy,qpz, + rpx,rpy,rpz, + spx,spy,spz); + } + + FT compute_denom(const FT &px, const FT &py, const FT &pz, + const FT &qx, const FT &qy, const FT &qz, + const FT &rx, const FT &ry, const FT &rz) const + { + const FT qpx = qx-px; + const FT qpy = qy-py; + const FT qpz = qz-pz; + const FT rpx = rx-px; + const FT rpy = ry-py; + const FT rpz = rz-pz; + const FT sx = qpy*rpz-qpz*rpy; + const FT sy = qpz*rpx-qpx*rpz; + const FT sz = qpx*rpy-qpy*rpx; + + return determinant(qpx,qpy,qpz, + rpx,rpy,rpz, + sx,sy,sz); + } + +}; + +/** + * @class Robust_circumcenter_filtered_traits_3 + */ +template +struct Robust_circumcenter_filtered_traits_3 +: public K +{ + typedef CGAL::Robust_filtered_construct_circumcenter_3 + Construct_circumcenter_3; + + Construct_circumcenter_3 + construct_circumcenter_3_object() const + { return Construct_circumcenter_3(); } + +}; // end class Robust_circumcenter_filtered_traits_3 + + +} // end namespace CGAL + +#endif // ROBUST_CIRCUMCENTER_FILTERED_TRAITS_3_H_ diff --git a/Surface_reconstruction_points_3/include/CGAL/poisson_refine_triangulation.h b/Surface_reconstruction_points_3/include/CGAL/poisson_refine_triangulation.h index 5eaa5a7df41..013f699eb26 100644 --- a/Surface_reconstruction_points_3/include/CGAL/poisson_refine_triangulation.h +++ b/Surface_reconstruction_points_3/include/CGAL/poisson_refine_triangulation.h @@ -47,12 +47,12 @@ class Poisson_mesher_level_impl_base : { typedef Mesh_3::Refine_tets_with_oracle_base Base; +public: // Inherited methods and fields used below using Base::triangulation_ref_impl; using Base::oracle; using Base::surface; -public: typedef typename Tr::Geom_traits Geom_traits; typedef typename Tr::Vertex_handle Vertex_handle; typedef typename Tr::Cell_handle Cell_handle; @@ -192,8 +192,9 @@ public: /// bad means badly shaped or too big). /// @return the number of vertices inserted. /// -/// @commentheading Precondition: -/// convergence is guaranteed if radius_edge_ratio_bound >= 1.0. +/// @commentheading Preconditions: +/// - Tr must use a geometric traits with robust circumcenter computation. +/// - convergence is guaranteed if radius_edge_ratio_bound >= 1.0. /// /// @commentheading Template Parameters: /// @param Tr 3D Delaunay triangulation. diff --git a/Surface_reconstruction_points_3/test/Surface_reconstruction_points_3/CMakeLists.txt b/Surface_reconstruction_points_3/test/Surface_reconstruction_points_3/CMakeLists.txt index 5548331a82b..9ffb6b1904a 100644 --- a/Surface_reconstruction_points_3/test/Surface_reconstruction_points_3/CMakeLists.txt +++ b/Surface_reconstruction_points_3/test/Surface_reconstruction_points_3/CMakeLists.txt @@ -54,6 +54,7 @@ if ( CGAL_FOUND ) # Temporary debugging stuff ADD_DEFINITIONS( "-DDEBUG_TRACE" ) + ADD_DEFINITIONS( "-DCGAL_PROFILE" ) # Executables that do *not* require BLAS, LAPACK nor TAUCS create_single_source_cgal_program( "APSS_reconstruction_test.cpp" )