diff options
Diffstat (limited to 'src/Coxeter_triangulation')
-rw-r--r-- | src/Coxeter_triangulation/include/gudhi/Implicit_manifold_intersection_oracle.h | 25 | ||||
-rw-r--r-- | src/Coxeter_triangulation/test/manifold_tracing_test.cpp | 4 |
2 files changed, 23 insertions, 6 deletions
diff --git a/src/Coxeter_triangulation/include/gudhi/Implicit_manifold_intersection_oracle.h b/src/Coxeter_triangulation/include/gudhi/Implicit_manifold_intersection_oracle.h index 25f688de..dca97331 100644 --- a/src/Coxeter_triangulation/include/gudhi/Implicit_manifold_intersection_oracle.h +++ b/src/Coxeter_triangulation/include/gudhi/Implicit_manifold_intersection_oracle.h @@ -17,8 +17,11 @@ #include <gudhi/Functions/Constant_function.h> #include <gudhi/Functions/PL_approximation.h> #include <gudhi/Query_result.h> +#include <gudhi/Debug_utils.h> // for GUDHI_CHECK #include <vector> +#include <limits> // for std::numeric_limits<> +#include <cmath> // for std::fabs namespace Gudhi { @@ -53,6 +56,11 @@ class Implicit_manifold_intersection_oracle { z(0) = 1; for (std::size_t i = 1; i < cod_d + 1; ++i) z(i) = 0; Eigen::VectorXd lambda = matrix.colPivHouseholderQr().solve(z); + if (!z.isApprox(matrix*lambda)) { + // NaN non valid results + for (std::size_t i = 0; i < (std::size_t)lambda.size(); ++i) lambda(i) = + std::numeric_limits<double>::quiet_NaN(); + } return lambda; } @@ -75,6 +83,11 @@ class Implicit_manifold_intersection_oracle { z(0) = 1; for (std::size_t i = 1; i < cod_d + 2; ++i) z(i) = 0; Eigen::VectorXd lambda = matrix.colPivHouseholderQr().solve(z); + if (!z.isApprox(matrix*lambda)) { + // NaN non valid results + for (std::size_t i = 0; i < (std::size_t)lambda.size(); ++i) lambda(i) = + std::numeric_limits<double>::quiet_NaN(); + } return lambda; } @@ -85,10 +98,13 @@ class Implicit_manifold_intersection_oracle { using QR = Query_result<Simplex_handle>; std::size_t amb_d = triangulation.dimension(); std::size_t cod_d = simplex.dimension(); - - for (std::size_t i = 0; i < (std::size_t)lambda.size(); ++i) - if (lambda(i) < 0 || lambda(i) > 1) return QR({Eigen::VectorXd(), false}); - + for (std::size_t i = 0; i < (std::size_t)lambda.size(); ++i) { + if (std::isnan(lambda(i))) return QR({Eigen::VectorXd(), false}); + GUDHI_CHECK((std::fabs(lambda(i) - 1.) > std::numeric_limits<double>::epsilon() && + std::fabs(lambda(i) - 0.) > std::numeric_limits<double>::epsilon()), + std::invalid_argument("A vertex of the triangulation lies exactly on the manifold")); + if (lambda(i) < 0. || lambda(i) > 1.) return QR({Eigen::VectorXd(), false}); + } Eigen::MatrixXd vertex_matrix(cod_d + 1, amb_d); auto v_range = simplex.vertex_range(); auto v_it = v_range.begin(); @@ -159,6 +175,7 @@ class Implicit_manifold_intersection_oracle { template <class Simplex_handle, class Triangulation> Query_result<Simplex_handle> intersects_boundary(const Simplex_handle& simplex, const Triangulation& triangulation) const { + //std::cout << "intersects_boundary" << std::endl; Eigen::VectorXd lambda = compute_boundary_lambda(simplex, triangulation); return intersection_result(lambda, simplex, triangulation); } diff --git a/src/Coxeter_triangulation/test/manifold_tracing_test.cpp b/src/Coxeter_triangulation/test/manifold_tracing_test.cpp index 68b79487..2d7d6087 100644 --- a/src/Coxeter_triangulation/test/manifold_tracing_test.cpp +++ b/src/Coxeter_triangulation/test/manifold_tracing_test.cpp @@ -40,7 +40,7 @@ BOOST_AUTO_TEST_CASE(manifold_tracing) { BOOST_CHECK(si_pair.second.size() == (long int)oracle.function().amb_d()); } std::cout << "out_simplex_map.size() = " << out_simplex_map.size() << "\n"; - BOOST_CHECK(out_simplex_map.size() == 1140); + BOOST_CHECK(out_simplex_map.size() == 1118); // manifold with boundary Function_Sm_in_Rd fun_boundary(3.0, 2, fun_sph.seed()); @@ -58,5 +58,5 @@ BOOST_AUTO_TEST_CASE(manifold_tracing) { BOOST_CHECK(si_pair.second.size() == (long int)oracle.function().amb_d()); } std::cout << "boundary_simplex_map.size() = " << boundary_simplex_map.size() << "\n"; - BOOST_CHECK(boundary_simplex_map.size() == 55); + BOOST_CHECK(boundary_simplex_map.size() == 54); } |