diff options
Diffstat (limited to 'src/Coxeter_triangulation/include/gudhi/Implicit_manifold_intersection_oracle.h')
| -rw-r--r-- | src/Coxeter_triangulation/include/gudhi/Implicit_manifold_intersection_oracle.h | 261 |
1 files changed, 261 insertions, 0 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 new file mode 100644 index 00000000..277f8b6c --- /dev/null +++ b/src/Coxeter_triangulation/include/gudhi/Implicit_manifold_intersection_oracle.h @@ -0,0 +1,261 @@ +/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT. + * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. + * Author(s): Siargey Kachanovich + * + * Copyright (C) 2019 Inria + * + * Modification(s): + * - YYYY/MM Author: Description of the modification + */ + +#ifndef IMPLICIT_MANIFOLD_INTERSECTION_ORACLE_H_ +#define IMPLICIT_MANIFOLD_INTERSECTION_ORACLE_H_ + +#include <Eigen/Dense> + +#include <gudhi/Permutahedral_representation/face_from_indices.h> +#include <gudhi/Functions/Constant_function.h> +#include <gudhi/Functions/PL_approximation.h> +#include <gudhi/Coxeter_triangulation/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 { + +namespace coxeter_triangulation { + +/** \class Implicit_manifold_intersection_oracle + * \brief An oracle that supports the intersection query on an implicit manifold. + * + * \tparam Function_ The function template parameter. Should be a model of + * the concept FunctionForImplicitManifold. + * \tparam Domain_function_ The domain function template parameter. Should be a model of + * the concept FunctionForImplicitManifold. + * + * \ingroup coxeter_triangulation + */ +template <class Function_, class Domain_function_ = Constant_function> +class Implicit_manifold_intersection_oracle { + /* Computes the affine coordinates of the intersection point of the implicit manifold + * and the affine hull of the simplex. */ + template <class Simplex_handle, class Triangulation> + Eigen::VectorXd compute_lambda(const Simplex_handle& simplex, const Triangulation& triangulation) const { + std::size_t cod_d = this->cod_d(); + Eigen::MatrixXd matrix(cod_d + 1, cod_d + 1); + for (std::size_t i = 0; i < cod_d + 1; ++i) matrix(0, i) = 1; + std::size_t j = 0; + for (auto v : simplex.vertex_range()) { + Eigen::VectorXd v_coords = fun_(triangulation.cartesian_coordinates(v)); + for (std::size_t i = 1; i < cod_d + 1; ++i) matrix(i, j) = v_coords(i - 1); + j++; + } + Eigen::VectorXd z(cod_d + 1); + 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; + } + + /* Computes the affine coordinates of the intersection point of the boundary + * of the implicit manifold and the affine hull of the simplex. */ + template <class Simplex_handle, class Triangulation> + Eigen::VectorXd compute_boundary_lambda(const Simplex_handle& simplex, const Triangulation& triangulation) const { + std::size_t cod_d = this->cod_d(); + Eigen::MatrixXd matrix(cod_d + 2, cod_d + 2); + for (std::size_t i = 0; i < cod_d + 2; ++i) matrix(0, i) = 1; + std::size_t j = 0; + for (auto v : simplex.vertex_range()) { + Eigen::VectorXd v_coords = fun_(triangulation.cartesian_coordinates(v)); + for (std::size_t i = 1; i < cod_d + 1; ++i) matrix(i, j) = v_coords(i - 1); + Eigen::VectorXd bv_coords = domain_fun_(triangulation.cartesian_coordinates(v)); + matrix(cod_d + 1, j) = bv_coords(0); + j++; + } + Eigen::VectorXd z(cod_d + 2); + 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; + } + + /* Computes the intersection result for a given simplex in a triangulation. */ + template <class Simplex_handle, class Triangulation> + Query_result<Simplex_handle> intersection_result(const Eigen::VectorXd& lambda, const Simplex_handle& simplex, + const Triangulation& triangulation) const { + 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 (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(); + for (std::size_t i = 0; i < cod_d + 1 && v_it != v_range.end(); ++v_it, ++i) { + Eigen::VectorXd v_coords = triangulation.cartesian_coordinates(*v_it); + for (std::size_t j = 0; j < amb_d; ++j) vertex_matrix(i, j) = v_coords(j); + } + Eigen::VectorXd intersection = lambda.transpose() * vertex_matrix; + return QR({intersection, true}); + } + + public: + /** \brief Ambient dimension of the implicit manifold. */ + std::size_t amb_d() const { return fun_.amb_d(); } + + /** \brief Codimension of the implicit manifold. */ + std::size_t cod_d() const { return fun_.cod_d(); } + + /** \brief The seed point of the implicit manifold. */ + Eigen::VectorXd seed() const { return fun_.seed(); } + + /** \brief Intersection query with the relative interior of the manifold. + * + * \details The returned structure Query_result contains the boolean value + * that is true only if the intersection point of the query simplex and + * the relative interior of the manifold exists, the intersection point + * and the face of the query simplex that contains + * the intersection point. + * + * \tparam Simplex_handle The class of the query simplex. + * Needs to be a model of the concept SimplexInCoxeterTriangulation. + * \tparam Triangulation The class of the triangulation. + * Needs to be a model of the concept TriangulationForManifoldTracing. + * + * @param[in] simplex The query simplex. The dimension of the simplex + * should be the same as the codimension of the manifold + * (the codomain dimension of the function). + * @param[in] triangulation The ambient triangulation. The dimension of + * the triangulation should be the same as the ambient dimension of the manifold + * (the domain dimension of the function). + */ + template <class Simplex_handle, class Triangulation> + Query_result<Simplex_handle> intersects(const Simplex_handle& simplex, const Triangulation& triangulation) const { + Eigen::VectorXd lambda = compute_lambda(simplex, triangulation); + return intersection_result(lambda, simplex, triangulation); + } + + /** \brief Intersection query with the boundary of the manifold. + * + * \details The returned structure Query_result contains the boolean value + * that is true only if the intersection point of the query simplex and + * the boundary of the manifold exists, the intersection point + * and the face of the query simplex that contains + * the intersection point. + * + * \tparam Simplex_handle The class of the query simplex. + * Needs to be a model of the concept SimplexInCoxeterTriangulation. + * \tparam Triangulation The class of the triangulation. + * Needs to be a model of the concept TriangulationForManifoldTracing. + * + * @param[in] simplex The query simplex. The dimension of the simplex + * should be the same as the codimension of the boundary of the manifold + * (the codomain dimension of the function + 1). + * @param[in] triangulation The ambient triangulation. The dimension of + * the triangulation should be the same as the ambient dimension of the manifold + * (the domain dimension of the function). + */ + 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); + } + + /** \brief Returns true if the input point lies inside the piecewise-linear + * domain induced by the given ambient triangulation that defines the relative + * interior of the piecewise-linear approximation of the manifold. + * + * @param p The input point. Needs to have the same dimension as the ambient + * dimension of the manifold (the domain dimension of the function). + * @param triangulation The ambient triangulation. Needs to have the same + * dimension as the ambient dimension of the manifold + * (the domain dimension of the function). + */ + template <class Triangulation> + bool lies_in_domain(const Eigen::VectorXd& p, const Triangulation& triangulation) const { + Eigen::VectorXd pl_p = make_pl_approximation(domain_fun_, triangulation)(p); + return pl_p(0) < 0; + } + + /** \brief Returns the function that defines the interior of the manifold */ + const Function_& function() const { return fun_; } + + /** \brief Constructs an intersection oracle for an implicit manifold potentially + * with boundary from given function and domain. + * + * @param function The input function that represents the implicit manifold + * before the restriction with the domain. + * @param domain_function The input domain function that can be used to define an implicit + * manifold with boundary. + */ + Implicit_manifold_intersection_oracle(const Function_& function, const Domain_function_& domain_function) + : fun_(function), domain_fun_(domain_function) {} + + /** \brief Constructs an intersection oracle for an implicit manifold + * without boundary from a given function. + * + * \details To use this constructor, the template Domain_function_ needs to be left + * at its default value (Gudhi::coxeter_triangulation::Constant_function). + * + * @param function The input function that represents the implicit manifold + * without boundary. + */ + Implicit_manifold_intersection_oracle(const Function_& function) + : fun_(function), domain_fun_(function.amb_d(), 1, Eigen::VectorXd::Constant(1, -1)) {} + + private: + Function_ fun_; + Domain_function_ domain_fun_; +}; + +/** \brief Static constructor of an intersection oracle from a function with a domain. + * + * @param function The input function that represents the implicit manifold + * before the restriction with the domain. + * @param domain_function The input domain function that can be used to define an implicit + * manifold with boundary. + * + * \ingroup coxeter_triangulation + */ +template <class Function_, class Domain_function_> +Implicit_manifold_intersection_oracle<Function_, Domain_function_> make_oracle( + const Function_& function, const Domain_function_& domain_function) { + return Implicit_manifold_intersection_oracle<Function_, Domain_function_>(function, domain_function); +} + +/** \brief Static constructor of an intersection oracle from a function without a domain. + * + * @param function The input function that represents the implicit manifold + * without boundary. + * + * \ingroup coxeter_triangulation + */ +template <class Function_> +Implicit_manifold_intersection_oracle<Function_> make_oracle(const Function_& function) { + return Implicit_manifold_intersection_oracle<Function_>(function); +} + +} // namespace coxeter_triangulation + +} // namespace Gudhi + +#endif |