comb_cross_field.cpp 5.4 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148
  1. // This file is part of libigl, a simple c++ geometry processing library.
  2. //
  3. // Copyright (C) 2014 Daniele Panozzo <daniele.panozzo@gmail.com>, Olga Diamanti <olga.diam@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 "comb_cross_field.h"
  9. #include <vector>
  10. #include <deque>
  11. #include <Eigen/Geometry>
  12. #include "per_face_normals.h"
  13. #include "is_border_vertex.h"
  14. #include "rotation_matrix_from_directions.h"
  15. #include "triangle_triangle_adjacency.h"
  16. namespace igl {
  17. template <typename DerivedV, typename DerivedF>
  18. class Comb
  19. {
  20. public:
  21. const Eigen::PlainObjectBase<DerivedV> &V;
  22. const Eigen::PlainObjectBase<DerivedF> &F;
  23. const Eigen::PlainObjectBase<DerivedV> &PD1;
  24. const Eigen::PlainObjectBase<DerivedV> &PD2;
  25. DerivedV N;
  26. private:
  27. // internal
  28. DerivedF TT;
  29. DerivedF TTi;
  30. private:
  31. static inline double Sign(double a){return (double)((a>0)?+1:-1);}
  32. private:
  33. // returns the 90 deg rotation of a (around n) most similar to target b
  34. /// a and b should be in the same plane orthogonal to N
  35. static inline Eigen::Matrix<typename DerivedV::Scalar, 3, 1> K_PI_new(const Eigen::Matrix<typename DerivedV::Scalar, 3, 1>& a,
  36. const Eigen::Matrix<typename DerivedV::Scalar, 3, 1>& b,
  37. const Eigen::Matrix<typename DerivedV::Scalar, 3, 1>& n)
  38. {
  39. Eigen::Matrix<typename DerivedV::Scalar, 3, 1> c = (a.cross(n)).normalized();
  40. typename DerivedV::Scalar scorea = a.dot(b);
  41. typename DerivedV::Scalar scorec = c.dot(b);
  42. if (fabs(scorea)>=fabs(scorec))
  43. return a*Sign(scorea);
  44. else
  45. return c*Sign(scorec);
  46. }
  47. public:
  48. inline Comb(const Eigen::PlainObjectBase<DerivedV> &_V,
  49. const Eigen::PlainObjectBase<DerivedF> &_F,
  50. const Eigen::PlainObjectBase<DerivedV> &_PD1,
  51. const Eigen::PlainObjectBase<DerivedV> &_PD2
  52. ):
  53. V(_V),
  54. F(_F),
  55. PD1(_PD1),
  56. PD2(_PD2)
  57. {
  58. igl::per_face_normals(V,F,N);
  59. igl::triangle_triangle_adjacency(F,TT,TTi);
  60. }
  61. inline void comb(Eigen::PlainObjectBase<DerivedV> &PD1out,
  62. Eigen::PlainObjectBase<DerivedV> &PD2out)
  63. {
  64. // PD1out = PD1;
  65. // PD2out = PD2;
  66. PD1out.setZero(F.rows(),3);PD1out<<PD1;
  67. PD2out.setZero(F.rows(),3);PD2out<<PD2;
  68. Eigen::VectorXi mark = Eigen::VectorXi::Constant(F.rows(),false);
  69. std::deque<int> d;
  70. d.push_back(0);
  71. mark(0) = true;
  72. while (!d.empty())
  73. {
  74. int f0 = d.at(0);
  75. d.pop_front();
  76. for (int k=0; k<3; k++)
  77. {
  78. int f1 = TT(f0,k);
  79. if (f1==-1) continue;
  80. if (mark(f1)) continue;
  81. Eigen::Matrix<typename DerivedV::Scalar, 3, 1> dir0 = PD1out.row(f0);
  82. Eigen::Matrix<typename DerivedV::Scalar, 3, 1> dir1 = PD1out.row(f1);
  83. Eigen::Matrix<typename DerivedV::Scalar, 3, 1> n0 = N.row(f0);
  84. Eigen::Matrix<typename DerivedV::Scalar, 3, 1> n1 = N.row(f1);
  85. Eigen::Matrix<typename DerivedV::Scalar, 3, 1> dir0Rot = igl::rotation_matrix_from_directions(n0, n1)*dir0;
  86. dir0Rot.normalize();
  87. Eigen::Matrix<typename DerivedV::Scalar, 3, 1> targD = K_PI_new(dir1,dir0Rot,n1);
  88. PD1out.row(f1) = targD;
  89. PD2out.row(f1) = n1.cross(targD).normalized();
  90. mark(f1) = true;
  91. d.push_back(f1);
  92. }
  93. }
  94. // everything should be marked
  95. for (int i=0; i<F.rows(); i++)
  96. {
  97. assert(mark(i));
  98. }
  99. }
  100. };
  101. }
  102. template <typename DerivedV, typename DerivedF>
  103. IGL_INLINE void igl::comb_cross_field(const Eigen::PlainObjectBase<DerivedV> &V,
  104. const Eigen::PlainObjectBase<DerivedF> &F,
  105. const Eigen::PlainObjectBase<DerivedV> &PD1,
  106. const Eigen::PlainObjectBase<DerivedV> &PD2,
  107. Eigen::PlainObjectBase<DerivedV> &PD1out,
  108. Eigen::PlainObjectBase<DerivedV> &PD2out)
  109. {
  110. igl::Comb<DerivedV, DerivedF> cmb(V, F, PD1, PD2);
  111. cmb.comb(PD1out, PD2out);
  112. }
  113. #ifdef IGL_STATIC_LIBRARY
  114. // Explicit template instantiation
  115. template void igl::comb_cross_field<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> >&);
  116. template void igl::comb_cross_field<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
  117. #endif