is_intrinsic_delaunay.cpp 3.7 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103
  1. // This file is part of libigl, a simple c++ geometry processing library.
  2. //
  3. // Copyright (C) 2018 Alec Jacobson <alecjacobson@gmail.com>
  4. //
  5. // This Source Code Form is subject to the terms of the Mozilla Public License
  6. // v. 2.0. If a copy of the MPL was not distributed with this file, You can
  7. // obtain one at http://mozilla.org/MPL/2.0/.
  8. #include "is_intrinsic_delaunay.h"
  9. #include "unique_edge_map.h"
  10. #include "tan_half_angle.h"
  11. #include <cassert>
  12. template <
  13. typename Derivedl,
  14. typename DerivedF,
  15. typename DerivedD>
  16. IGL_INLINE void igl::is_intrinsic_delaunay(
  17. const Eigen::MatrixBase<Derivedl> & l,
  18. const Eigen::MatrixBase<DerivedF> & F,
  19. Eigen::PlainObjectBase<DerivedD> & D)
  20. {
  21. typedef Eigen::Matrix<typename DerivedF::Scalar,Eigen::Dynamic,2> MatrixX2I;
  22. typedef Eigen::Matrix<typename DerivedF::Scalar,Eigen::Dynamic,1> VectorXI;
  23. MatrixX2I E,uE;
  24. VectorXI EMAP;
  25. std::vector<std::vector<typename DerivedF::Scalar> > uE2E;
  26. igl::unique_edge_map(F, E, uE, EMAP, uE2E);
  27. const int num_faces = F.rows();
  28. D.setConstant(F.rows(),F.cols(),false);
  29. // loop over all unique edges
  30. for(int ue = 0;ue < uE2E.size(); ue++)
  31. {
  32. const bool ue_is_d = is_intrinsic_delaunay(l,F,uE2E,ue);
  33. // Set for all instances
  34. for(int e = 0;e<uE2E[ue].size();e++)
  35. {
  36. D( uE2E[ue][e]%num_faces, uE2E[ue][e]/num_faces) = ue_is_d;
  37. }
  38. }
  39. }
  40. template <
  41. typename Derivedl,
  42. typename DerivedF,
  43. typename uE2EType,
  44. typename ueiType>
  45. IGL_INLINE bool igl::is_intrinsic_delaunay(
  46. const Eigen::MatrixBase<Derivedl> & l,
  47. const Eigen::MatrixBase<DerivedF> & F,
  48. const std::vector<std::vector<uE2EType> > & uE2E,
  49. const ueiType uei)
  50. {
  51. if(uE2E[uei].size() == 1) return true;
  52. if(uE2E[uei].size() > 2) return false;
  53. typedef typename Derivedl::Scalar Scalar;
  54. typedef typename DerivedF::Scalar Index;
  55. const auto cot_alpha = [](
  56. const Scalar & a,
  57. const Scalar & b,
  58. const Scalar & c)->Scalar
  59. {
  60. // Fisher 2007
  61. const Scalar t = tan_half_angle(a,b,c);
  62. return (1.0-t*t)/(2*t);
  63. };
  64. // . //
  65. // /|\ //
  66. // a/ | \d //
  67. // / e \ //
  68. // / | \ //
  69. // .α---|-f-β. //
  70. // \ | / //
  71. // \ | / //
  72. // b\ | /c //
  73. // \|/ //
  74. // . //
  75. const Index num_faces = F.rows();
  76. assert(uE2E[uei].size() == 2 && "edge should have 2 incident faces");
  77. const Index he1 = uE2E[uei][0];
  78. const Index he2 = uE2E[uei][1];
  79. const Index f1 = he1%num_faces;
  80. const Index c1 = he1/num_faces;
  81. const Index f2 = he2%num_faces;
  82. const Index c2 = he2/num_faces;
  83. assert( std::abs(l(f1,c1)-l(f2,c2) < igl::EPS<Scalar>()) );
  84. const Scalar e = l(f1,c1);
  85. const Scalar a = l(f1,(c1+1)%3);
  86. const Scalar b = l(f1,(c1+2)%3);
  87. const Scalar c = l(f2,(c2+1)%3);
  88. const Scalar d = l(f2,(c2+2)%3);
  89. const Scalar w = cot_alpha(e,a,b) + cot_alpha(e,c,d);
  90. return w >= 0;
  91. }
  92. #ifdef IGL_STATIC_LIBRARY
  93. // Explicit template instantiation
  94. // generated by autoexplicit.sh
  95. template bool igl::is_intrinsic_delaunay<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, int, unsigned long>(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, std::vector<std::vector<int, std::allocator<int> >, std::allocator<std::vector<int, std::allocator<int> > > > const&, unsigned long);
  96. template void igl::is_intrinsic_delaunay<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<bool, -1, -1, 0, -1, -1> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<bool, -1, -1, 0, -1, -1> >&);
  97. #endif