doublearea.cpp 13 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250
  1. // This file is part of libigl, a simple c++ geometry processing library.
  2. //
  3. // Copyright (C) 2013 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 "doublearea.h"
  9. #include "edge_lengths.h"
  10. #include "parallel_for.h"
  11. #include "sort.h"
  12. #include <cassert>
  13. #include <iostream>
  14. #include <limits>
  15. template <typename DerivedV, typename DerivedF, typename DeriveddblA>
  16. IGL_INLINE void igl::doublearea(
  17. const Eigen::MatrixBase<DerivedV> & V,
  18. const Eigen::MatrixBase<DerivedF> & F,
  19. Eigen::PlainObjectBase<DeriveddblA> & dblA)
  20. {
  21. // quads are handled by a specialized function
  22. if (F.cols() == 4) return doublearea_quad(V,F,dblA);
  23. const int dim = V.cols();
  24. // Only support triangles
  25. assert(F.cols() == 3);
  26. const size_t m = F.rows();
  27. // Compute edge lengths
  28. Eigen::Matrix<typename DerivedV::Scalar, Eigen::Dynamic, 3> l;
  29. // Projected area helper
  30. const auto & proj_doublearea =
  31. [&V,&F](const int x, const int y, const int f)->double
  32. {
  33. auto rx = V(F(f,0),x)-V(F(f,2),x);
  34. auto sx = V(F(f,1),x)-V(F(f,2),x);
  35. auto ry = V(F(f,0),y)-V(F(f,2),y);
  36. auto sy = V(F(f,1),y)-V(F(f,2),y);
  37. return rx*sy - ry*sx;
  38. };
  39. switch(dim)
  40. {
  41. case 3:
  42. {
  43. dblA = DeriveddblA::Zero(m,1);
  44. for(size_t f = 0;f<m;f++)
  45. {
  46. for(int d = 0;d<3;d++)
  47. {
  48. double dblAd = proj_doublearea(d,(d+1)%3,f);
  49. dblA(f) += dblAd*dblAd;
  50. }
  51. }
  52. dblA = dblA.array().sqrt().eval();
  53. break;
  54. }
  55. case 2:
  56. {
  57. dblA.resize(m,1);
  58. for(size_t f = 0;f<m;f++)
  59. {
  60. dblA(f) = proj_doublearea(0,1,f);
  61. }
  62. break;
  63. }
  64. default:
  65. {
  66. edge_lengths(V,F,l);
  67. return doublearea(l,0.,dblA);
  68. }
  69. }
  70. }
  71. template <
  72. typename DerivedA,
  73. typename DerivedB,
  74. typename DerivedC,
  75. typename DerivedD>
  76. IGL_INLINE void igl::doublearea(
  77. const Eigen::MatrixBase<DerivedA> & A,
  78. const Eigen::MatrixBase<DerivedB> & B,
  79. const Eigen::MatrixBase<DerivedC> & C,
  80. Eigen::PlainObjectBase<DerivedD> & D)
  81. {
  82. assert((B.cols() == A.cols()) && "dimensions of A and B should match");
  83. assert((C.cols() == A.cols()) && "dimensions of A and C should match");
  84. assert(A.rows() == B.rows() && "corners should have same length");
  85. assert(A.rows() == C.rows() && "corners should have same length");
  86. switch(A.cols())
  87. {
  88. case 2:
  89. {
  90. // For 2d compute signed area
  91. const auto & R = A-C;
  92. const auto & S = B-C;
  93. D = R.col(0).array()*S.col(1).array() - R.col(1).array()*S.col(0).array();
  94. break;
  95. }
  96. default:
  97. {
  98. Eigen::Matrix<typename DerivedD::Scalar,DerivedD::RowsAtCompileTime,3>
  99. uL(A.rows(),3);
  100. uL.col(0) = (B-C).rowwise().norm();
  101. uL.col(1) = (C-A).rowwise().norm();
  102. uL.col(2) = (A-B).rowwise().norm();
  103. doublearea(uL,D);
  104. }
  105. }
  106. }
  107. template <
  108. typename DerivedA,
  109. typename DerivedB,
  110. typename DerivedC>
  111. IGL_INLINE typename DerivedA::Scalar igl::doublearea_single(
  112. const Eigen::MatrixBase<DerivedA> & A,
  113. const Eigen::MatrixBase<DerivedB> & B,
  114. const Eigen::MatrixBase<DerivedC> & C)
  115. {
  116. assert(A.size() == 2 && "Vertices should be 2D");
  117. assert(B.size() == 2 && "Vertices should be 2D");
  118. assert(C.size() == 2 && "Vertices should be 2D");
  119. auto r = A-C;
  120. auto s = B-C;
  121. return r(0)*s(1) - r(1)*s(0);
  122. }
  123. template <typename Derivedl, typename DeriveddblA>
  124. IGL_INLINE void igl::doublearea(
  125. const Eigen::MatrixBase<Derivedl> & ul,
  126. Eigen::PlainObjectBase<DeriveddblA> & dblA)
  127. {
  128. // Default is to leave NaNs and fire asserts in debug mode
  129. return doublearea(
  130. ul,std::numeric_limits<typename Derivedl::Scalar>::quiet_NaN(),dblA);
  131. }
  132. template <typename Derivedl, typename DeriveddblA>
  133. IGL_INLINE void igl::doublearea(
  134. const Eigen::MatrixBase<Derivedl> & ul,
  135. const typename Derivedl::Scalar nan_replacement,
  136. Eigen::PlainObjectBase<DeriveddblA> & dblA)
  137. {
  138. using namespace Eigen;
  139. using namespace std;
  140. typedef typename Derivedl::Index Index;
  141. // Only support triangles
  142. assert(ul.cols() == 3);
  143. // Number of triangles
  144. const Index m = ul.rows();
  145. Eigen::Matrix<typename Derivedl::Scalar, Eigen::Dynamic, 3> l;
  146. MatrixXi _;
  147. //
  148. // "Miscalculating Area and Angles of a Needle-like Triangle"
  149. // https://people.eecs.berkeley.edu/~wkahan/Triangle.pdf
  150. igl::sort(ul,2,false,l,_);
  151. // semiperimeters
  152. //Matrix<typename Derivedl::Scalar,Dynamic,1> s = l.rowwise().sum()*0.5;
  153. //assert((Index)s.rows() == m);
  154. // resize output
  155. dblA.resize(l.rows(),1);
  156. parallel_for(
  157. m,
  158. [&l,&dblA,&nan_replacement](const int i)
  159. {
  160. // Kahan's Heron's formula
  161. typedef typename Derivedl::Scalar Scalar;
  162. const Scalar arg =
  163. (l(i,0)+(l(i,1)+l(i,2)))*
  164. (l(i,2)-(l(i,0)-l(i,1)))*
  165. (l(i,2)+(l(i,0)-l(i,1)))*
  166. (l(i,0)+(l(i,1)-l(i,2)));
  167. dblA(i) = 2.0*0.25*sqrt(arg);
  168. // Alec: If the input edge lengths were computed from floating point
  169. // vertex positions then there's no guarantee that they fulfill the
  170. // triangle inequality (in their floating point approximations). For
  171. // nearly degenerate triangles the round-off error during side-length
  172. // computation may be larger than (or rather smaller) than the height of
  173. // the triangle. In "Lecture Notes on Geometric Robustness" Shewchuck 09,
  174. // Section 3.1 http://www.cs.berkeley.edu/~jrs/meshpapers/robnotes.pdf,
  175. // he recommends computing the triangle areas for 2D and 3D using 2D
  176. // signed areas computed with determinants.
  177. assert(
  178. (nan_replacement == nan_replacement ||
  179. (l(i,2) - (l(i,0)-l(i,1)))>=0)
  180. && "Side lengths do not obey the triangle inequality.");
  181. if(dblA(i) != dblA(i))
  182. {
  183. dblA(i) = nan_replacement;
  184. }
  185. assert(dblA(i) == dblA(i) && "DOUBLEAREA() PRODUCED NaN");
  186. },
  187. 1000l);
  188. }
  189. template <typename DerivedV, typename DerivedF, typename DeriveddblA>
  190. IGL_INLINE void igl::doublearea_quad(
  191. const Eigen::MatrixBase<DerivedV> & V,
  192. const Eigen::MatrixBase<DerivedF> & F,
  193. Eigen::PlainObjectBase<DeriveddblA> & dblA)
  194. {
  195. assert(V.cols() == 3); // Only supports points in 3D
  196. assert(F.cols() == 4); // Only support quads
  197. const size_t m = F.rows();
  198. // Split the quads into triangles
  199. Eigen::MatrixXi Ft(F.rows()*2,3);
  200. for(size_t i=0; i<m;++i)
  201. {
  202. Ft.row(i*2 ) << F(i,0), F(i,1), F(i,2);
  203. Ft.row(i*2 + 1) << F(i,2), F(i,3), F(i,0);
  204. }
  205. // Compute areas
  206. Eigen::VectorXd doublearea_tri;
  207. igl::doublearea(V,Ft,doublearea_tri);
  208. dblA.resize(F.rows(),1);
  209. for(unsigned i=0; i<F.rows();++i)
  210. {
  211. dblA(i) = doublearea_tri(i*2) + doublearea_tri(i*2 + 1);
  212. }
  213. }
  214. #ifdef IGL_STATIC_LIBRARY
  215. // Explicit template instantiation
  216. // generated by autoexplicit.sh
  217. template void igl::doublearea<Eigen::Matrix<double, -1, -1, 1, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, 1, 0, -1, 1> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 1, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&);
  218. template void igl::doublearea<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, 1, 0, -1, 1> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&);
  219. template void igl::doublearea<Eigen::Matrix<float, -1, 3, 1, -1, 3>, Eigen::Matrix<unsigned int, -1, 3, 1, -1, 3>, Eigen::Matrix<float, -1, 1, 0, -1, 1> >(Eigen::MatrixBase<Eigen::Matrix<float, -1, 3, 1, -1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<unsigned int, -1, 3, 1, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<float, -1, 1, 0, -1, 1> >&);
  220. template void igl::doublearea<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<double, -1, 1, 0, -1, 1> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&);
  221. template void igl::doublearea<Eigen::Matrix<float, -1, 3, 1, -1, 3>, Eigen::Matrix<unsigned int, -1, -1, 1, -1, -1>, Eigen::Matrix<double, -1, 1, 0, -1, 1> >(Eigen::MatrixBase<Eigen::Matrix<float, -1, 3, 1, -1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<unsigned int, -1, -1, 1, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&);
  222. template void igl::doublearea<Eigen::Matrix<double, -1, 3, 1, -1, 3>, Eigen::Matrix<unsigned int, -1, -1, 1, -1, -1>, Eigen::Matrix<double, -1, 1, 0, -1, 1> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, 3, 1, -1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<unsigned int, -1, -1, 1, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&);
  223. template void igl::doublearea<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -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<double, -1, 1, 0, -1, 1> >&);
  224. template void igl::doublearea<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, 1, 0, -1, 1> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&);
  225. template void igl::doublearea<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -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<int, -1, 1, 0, -1, 1> >&);
  226. template void igl::doublearea<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<double, -1, 1, 0, -1, 1> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&);
  227. template void igl::doublearea<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
  228. template void igl::doublearea<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -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<double, -1, -1, 0, -1, -1> >&);
  229. template void igl::doublearea<Eigen::Matrix<double, 1, 3, 1, 1, 3>, Eigen::Matrix<double, 1, 3, 1, 1, 3>, Eigen::Matrix<double, 1, 3, 1, 1, 3>, Eigen::Matrix<double, 1, 1, 0, 1, 1> >(Eigen::MatrixBase<Eigen::Matrix<double, 1, 3, 1, 1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<double, 1, 3, 1, 1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<double, 1, 3, 1, 1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, 1, 1, 0, 1, 1> >&);
  230. template void igl::doublearea<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
  231. template void igl::doublearea<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
  232. template void igl::doublearea<Eigen::Matrix<double, -1, 3, 1, -1, 3>, Eigen::Matrix<int, -1, 3, 1, -1, 3>, Eigen::Matrix<double, -1, 1, 0, -1, 1> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, 3, 1, -1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 3, 1, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&);
  233. template Eigen::Matrix<double, -1, -1, 0, -1, -1>::Scalar igl::doublearea_single<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&);
  234. template Eigen::Matrix<double, 2, 1, 0, 2, 1>::Scalar igl::doublearea_single<Eigen::Matrix<double, 2, 1, 0, 2, 1>, Eigen::Matrix<double, 2, 1, 0, 2, 1>, Eigen::Matrix<double, 2, 1, 0, 2, 1> >(Eigen::MatrixBase<Eigen::Matrix<double, 2, 1, 0, 2, 1> > const&, Eigen::MatrixBase<Eigen::Matrix<double, 2, 1, 0, 2, 1> > const&, Eigen::MatrixBase<Eigen::Matrix<double, 2, 1, 0, 2, 1> > const&);
  235. #endif