Loading...
Searching...
No Matches
Implicit_manifold_intersection_oracle.h
1/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT.
2 * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details.
3 * Author(s): Siargey Kachanovich
4 *
5 * Copyright (C) 2019 Inria
6 *
7 * Modification(s):
8 * - YYYY/MM Author: Description of the modification
9 */
10
11#ifndef IMPLICIT_MANIFOLD_INTERSECTION_ORACLE_H_
12#define IMPLICIT_MANIFOLD_INTERSECTION_ORACLE_H_
13
14#include <Eigen/Dense>
15
16#include <gudhi/Permutahedral_representation/face_from_indices.h>
17#include <gudhi/Functions/Constant_function.h>
18#include <gudhi/Functions/PL_approximation.h>
19#include <gudhi/Coxeter_triangulation/Query_result.h>
20#include <gudhi/Debug_utils.h> // for GUDHI_CHECK
21
22#include <vector>
23#include <limits> // for std::numeric_limits<>
24#include <cmath> // for std::fabs
25
26namespace Gudhi {
27
28namespace coxeter_triangulation {
29
40template <class Function_, class Domain_function_ = Constant_function>
42 /* Computes the affine coordinates of the intersection point of the implicit manifold
43 * and the affine hull of the simplex. */
44 template <class Simplex_handle, class Triangulation>
45 Eigen::VectorXd compute_lambda(const Simplex_handle& simplex, const Triangulation& triangulation) const {
46 std::size_t cod_d = this->cod_d();
47 Eigen::MatrixXd matrix(cod_d + 1, cod_d + 1);
48 for (std::size_t i = 0; i < cod_d + 1; ++i) matrix(0, i) = 1;
49 std::size_t j = 0;
50 for (auto v : simplex.vertex_range()) {
51 Eigen::VectorXd v_coords = fun_(triangulation.cartesian_coordinates(v));
52 for (std::size_t i = 1; i < cod_d + 1; ++i) matrix(i, j) = v_coords(i - 1);
53 j++;
54 }
55 Eigen::VectorXd z(cod_d + 1);
56 z(0) = 1;
57 for (std::size_t i = 1; i < cod_d + 1; ++i) z(i) = 0;
58 Eigen::VectorXd lambda = matrix.colPivHouseholderQr().solve(z);
59 if (!z.isApprox(matrix*lambda)) {
60 // NaN non valid results
61 for (std::size_t i = 0; i < (std::size_t)lambda.size(); ++i) lambda(i) =
62 std::numeric_limits<double>::quiet_NaN();
63 }
64 return lambda;
65 }
66
67 /* Computes the affine coordinates of the intersection point of the boundary
68 * of the implicit manifold and the affine hull of the simplex. */
69 template <class Simplex_handle, class Triangulation>
70 Eigen::VectorXd compute_boundary_lambda(const Simplex_handle& simplex, const Triangulation& triangulation) const {
71 std::size_t cod_d = this->cod_d();
72 Eigen::MatrixXd matrix(cod_d + 2, cod_d + 2);
73 for (std::size_t i = 0; i < cod_d + 2; ++i) matrix(0, i) = 1;
74 std::size_t j = 0;
75 for (auto v : simplex.vertex_range()) {
76 Eigen::VectorXd v_coords = fun_(triangulation.cartesian_coordinates(v));
77 for (std::size_t i = 1; i < cod_d + 1; ++i) matrix(i, j) = v_coords(i - 1);
78 Eigen::VectorXd bv_coords = domain_fun_(triangulation.cartesian_coordinates(v));
79 matrix(cod_d + 1, j) = bv_coords(0);
80 j++;
81 }
82 Eigen::VectorXd z(cod_d + 2);
83 z(0) = 1;
84 for (std::size_t i = 1; i < cod_d + 2; ++i) z(i) = 0;
85 Eigen::VectorXd lambda = matrix.colPivHouseholderQr().solve(z);
86 if (!z.isApprox(matrix*lambda)) {
87 // NaN non valid results
88 for (std::size_t i = 0; i < (std::size_t)lambda.size(); ++i) lambda(i) =
89 std::numeric_limits<double>::quiet_NaN();
90 }
91 return lambda;
92 }
93
94 /* Computes the intersection result for a given simplex in a triangulation. */
95 template <class Simplex_handle, class Triangulation>
96 Query_result<Simplex_handle> intersection_result(const Eigen::VectorXd& lambda, const Simplex_handle& simplex,
97 const Triangulation& triangulation) const {
99 std::size_t amb_d = triangulation.dimension();
100 std::size_t cod_d = simplex.dimension();
101 for (std::size_t i = 0; i < (std::size_t)lambda.size(); ++i) {
102 if (std::isnan(lambda(i))) return QR({Eigen::VectorXd(), false});
103 GUDHI_CHECK((std::fabs(lambda(i) - 1.) > std::numeric_limits<double>::epsilon() &&
104 std::fabs(lambda(i) - 0.) > std::numeric_limits<double>::epsilon()),
105 std::invalid_argument("A vertex of the triangulation lies exactly on the manifold"));
106 if (lambda(i) < 0. || lambda(i) > 1.) return QR({Eigen::VectorXd(), false});
107 }
108 Eigen::MatrixXd vertex_matrix(cod_d + 1, amb_d);
109 auto v_range = simplex.vertex_range();
110 auto v_it = v_range.begin();
111 for (std::size_t i = 0; i < cod_d + 1 && v_it != v_range.end(); ++v_it, ++i) {
112 Eigen::VectorXd v_coords = triangulation.cartesian_coordinates(*v_it);
113 for (std::size_t j = 0; j < amb_d; ++j) vertex_matrix(i, j) = v_coords(j);
114 }
115 Eigen::VectorXd intersection = lambda.transpose() * vertex_matrix;
116 return QR({intersection, true});
117 }
118
119 public:
121 std::size_t amb_d() const { return fun_.amb_d(); }
122
124 std::size_t cod_d() const { return fun_.cod_d(); }
125
127 Eigen::VectorXd seed() const { return fun_.seed(); }
128
149 template <class Simplex_handle, class Triangulation>
150 Query_result<Simplex_handle> intersects(const Simplex_handle& simplex, const Triangulation& triangulation) const {
151 Eigen::VectorXd lambda = compute_lambda(simplex, triangulation);
152 return intersection_result(lambda, simplex, triangulation);
153 }
154
175 template <class Simplex_handle, class Triangulation>
177 const Triangulation& triangulation) const {
178 //std::cout << "intersects_boundary" << std::endl;
179 Eigen::VectorXd lambda = compute_boundary_lambda(simplex, triangulation);
180 return intersection_result(lambda, simplex, triangulation);
181 }
182
193 template <class Triangulation>
194 bool lies_in_domain(const Eigen::VectorXd& p, const Triangulation& triangulation) const {
195 Eigen::VectorXd pl_p = make_pl_approximation(domain_fun_, triangulation)(p);
196 return pl_p(0) < 0;
197 }
198
200 const Function_& function() const { return fun_; }
201
210 Implicit_manifold_intersection_oracle(const Function_& function, const Domain_function_& domain_function)
211 : fun_(function), domain_fun_(domain_function) {}
212
223 : fun_(function), domain_fun_(function.amb_d(), 1, Eigen::VectorXd::Constant(1, -1)) {}
224
225 private:
226 Function_ fun_;
227 Domain_function_ domain_fun_;
228};
229
239template <class Function_, class Domain_function_>
241 const Function_& function, const Domain_function_& domain_function) {
243}
244
252template <class Function_>
255}
256
257} // namespace coxeter_triangulation
258
259} // namespace Gudhi
260
261#endif
An oracle that supports the intersection query on an implicit manifold.
Definition: Implicit_manifold_intersection_oracle.h:41
Eigen::VectorXd seed() const
The seed point of the implicit manifold.
Definition: Implicit_manifold_intersection_oracle.h:127
Implicit_manifold_intersection_oracle(const Function_ &function, const Domain_function_ &domain_function)
Constructs an intersection oracle for an implicit manifold potentially with boundary from given funct...
Definition: Implicit_manifold_intersection_oracle.h:210
Query_result< Simplex_handle > intersects(const Simplex_handle &simplex, const Triangulation &triangulation) const
Intersection query with the relative interior of the manifold.
Definition: Implicit_manifold_intersection_oracle.h:150
Implicit_manifold_intersection_oracle(const Function_ &function)
Constructs an intersection oracle for an implicit manifold without boundary from a given function.
Definition: Implicit_manifold_intersection_oracle.h:222
const Function_ & function() const
Returns the function that defines the interior of the manifold.
Definition: Implicit_manifold_intersection_oracle.h:200
std::size_t amb_d() const
Ambient dimension of the implicit manifold.
Definition: Implicit_manifold_intersection_oracle.h:121
bool lies_in_domain(const Eigen::VectorXd &p, const Triangulation &triangulation) const
Returns true if the input point lies inside the piecewise-linear domain induced by the given ambient ...
Definition: Implicit_manifold_intersection_oracle.h:194
std::size_t cod_d() const
Codimension of the implicit manifold.
Definition: Implicit_manifold_intersection_oracle.h:124
Query_result< Simplex_handle > intersects_boundary(const Simplex_handle &simplex, const Triangulation &triangulation) const
Intersection query with the boundary of the manifold.
Definition: Implicit_manifold_intersection_oracle.h:176
PL_approximation< Function_, Triangulation_ > make_pl_approximation(const Function_ &function, const Triangulation_ &triangulation)
Static constructor of the piecewise-linear approximation of a function induced by an ambient triangul...
Definition: PL_approximation.h:102
Implicit_manifold_intersection_oracle< Function_, Domain_function_ > make_oracle(const Function_ &function, const Domain_function_ &domain_function)
Static constructor of an intersection oracle from a function with a domain.
Definition: Implicit_manifold_intersection_oracle.h:240
The result of a query by an oracle such as Implicit_manifold_intersection_oracle.
Definition: Query_result.h:26