Medusa  1.1
Coordinate Free Mehless Method implementation
compute_normal.hpp
Go to the documentation of this file.
1 #ifndef MEDUSA_BITS_DOMAINS_COMPUTE_NORMAL_HPP_
2 #define MEDUSA_BITS_DOMAINS_COMPUTE_NORMAL_HPP_
3 
9 #include "compute_normal_fwd.hpp"
10 #include <medusa/bits/utils/assert.hpp>
11 #include <Eigen/Core>
12 #include <Eigen/SVD>
13 #include <Eigen/LU>
14 
15 namespace mm {
16 namespace surface_fill_internal {
17 
19 template <typename scalar_t, int dim_from, int dim_to>
20 Vec<scalar_t, dim_to> compute_normal(Eigen::Matrix<scalar_t, dim_to, dim_from> J) {
21  static_assert(dim_from < dim_to, "At least one free dimension must be present.");
22  Eigen::JacobiSVD<decltype(J)> svd(J, Eigen::ComputeFullU);
23  assert_msg(svd.rank() == dim_from, "Jacobi matrix does not have full rank.");
24  Vec<scalar_t, dim_to> normal = svd.matrixU().col(dim_from);
25  // Find correct orientation
26  Eigen::Matrix<double, dim_from + 1, dim_from + 1> M;
27  M.template topLeftCorner<dim_from + 1, dim_from>() =
28  J.template topLeftCorner<dim_from + 1, dim_from>();
29  M.template rightCols<1>() = normal.template tail<dim_from + 1>();
30  normal *= (M.determinant() > 0) ? -1.0 : 1.0;
31  return normal;
32 }
34 
35 } // namespace surface_fill_internal
36 } // namespace mm
37 
38 #endif // MEDUSA_BITS_DOMAINS_COMPUTE_NORMAL_HPP_
mm
Root namespace for the whole library.
Definition: Gaussian.hpp:14
mm::Vec
Eigen::Matrix< scalar_t, dim, 1, Eigen::ColMajor|Eigen::AutoAlign, dim, 1 > Vec
Fixed size vector type, representing a mathematical 1d/2d/3d vector.
Definition: Vec_fwd.hpp:31
assert_msg
#define assert_msg(cond,...)
Assert with better error reporting.
Definition: assert.hpp:75
compute_normal_fwd.hpp
assert.hpp
mm::surface_fill_internal::compute_normal
Vec< scalar_t, dim_to > compute_normal(Eigen::Matrix< scalar_t, dim_to, dim_from > J)
Compute normal from given Jacobian.