doublearea.cpp 5.7 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126
  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 "sort.h"
  11. #include <cassert>
  12. #include <iostream>
  13. template <typename DerivedV, typename DerivedF, typename DeriveddblA>
  14. IGL_INLINE void igl::doublearea(
  15. const Eigen::PlainObjectBase<DerivedV> & V,
  16. const Eigen::PlainObjectBase<DerivedF> & F,
  17. Eigen::PlainObjectBase<DeriveddblA> & dblA)
  18. {
  19. const int dim = V.cols();
  20. // Only support triangles
  21. assert(F.cols() == 3);
  22. const int m = F.rows();
  23. // Compute edge lengths
  24. Eigen::PlainObjectBase<DerivedV> l;
  25. // "Lecture Notes on Geometric Robustness" Shewchuck 09, Section 3.1
  26. // http://www.cs.berkeley.edu/~jrs/meshpapers/robnotes.pdf
  27. switch(dim)
  28. {
  29. case 3:
  30. {
  31. dblA = Eigen::PlainObjectBase<DeriveddblA>::Zero(m,1);
  32. for(int d = 0;d<3;d++)
  33. {
  34. Eigen::Matrix<typename DerivedV::Scalar, DerivedV::RowsAtCompileTime, 2> Vd(V.rows(),2);
  35. Vd << V.col(d), V.col((d+1)%3);
  36. Eigen::PlainObjectBase<DeriveddblA> dblAd;
  37. doublearea(Vd,F,dblAd);
  38. dblA.array() += dblAd.array().square();
  39. }
  40. dblA = dblA.array().sqrt().eval();
  41. break;
  42. }
  43. case 2:
  44. {
  45. dblA.resize(m,1);
  46. for(int f = 0;f<m;f++)
  47. {
  48. auto r = V.row(F(f,0))-V.row(F(f,2));
  49. auto s = V.row(F(f,1))-V.row(F(f,2));
  50. dblA(f) = r(0)*s(1) - r(1)*s(0);
  51. }
  52. break;
  53. }
  54. default:
  55. {
  56. edge_lengths(V,F,l);
  57. return doublearea(l,dblA);
  58. }
  59. }
  60. }
  61. template <
  62. typename DerivedA,
  63. typename DerivedB,
  64. typename DerivedC>
  65. IGL_INLINE typename DerivedA::Scalar igl::doublearea_single(
  66. const Eigen::PlainObjectBase<DerivedA> & A,
  67. const Eigen::PlainObjectBase<DerivedB> & B,
  68. const Eigen::PlainObjectBase<DerivedC> & C)
  69. {
  70. auto r = A-C;
  71. auto s = B-C;
  72. return r(0)*s(1) - r(1)*s(0);
  73. }
  74. template <typename Derivedl, typename DeriveddblA>
  75. IGL_INLINE void igl::doublearea(
  76. const Eigen::PlainObjectBase<Derivedl> & ul,
  77. Eigen::PlainObjectBase<DeriveddblA> & dblA)
  78. {
  79. using namespace Eigen;
  80. using namespace std;
  81. // Only support triangles
  82. assert(ul.cols() == 3);
  83. // Number of triangles
  84. const int m = ul.rows();
  85. Eigen::PlainObjectBase<Derivedl> l;
  86. MatrixXi _;
  87. sort(ul,2,false,l,_);
  88. // semiperimeters
  89. Matrix<typename Derivedl::Scalar,Dynamic,1> s = l.rowwise().sum()*0.5;
  90. assert(s.rows() == m);
  91. // resize output
  92. dblA.resize(l.rows(),1);
  93. for(int i = 0;i<m;i++)
  94. {
  95. //// Heron's formula for area
  96. //const typename Derivedl::Scalar arg =
  97. // s(i)*(s(i)-l(i,0))*(s(i)-l(i,1))*(s(i)-l(i,2));
  98. //assert(arg>=0);
  99. //dblA(i) = 2.0*sqrt(arg);
  100. // Kahan's Heron's formula
  101. const typename Derivedl::Scalar arg =
  102. (l(i,0)+(l(i,1)+l(i,2)))*
  103. (l(i,2)-(l(i,0)-l(i,1)))*
  104. (l(i,2)+(l(i,0)-l(i,1)))*
  105. (l(i,0)+(l(i,1)-l(i,2)));
  106. dblA(i) = 2.0*0.25*sqrt(arg);
  107. assert( l(i,2) - (l(i,0)-l(i,1)) && "FAILED KAHAN'S ASSERTION");
  108. assert(dblA(i) == dblA(i) && "DOUBLEAREA() PRODUCED NaN");
  109. }
  110. }
  111. #ifndef IGL_HEADER_ONLY
  112. // Explicit template specialization
  113. // generated by autoexplicit.sh
  114. template void igl::doublearea<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<double, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&);
  115. // generated by autoexplicit.sh
  116. 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::PlainObjectBase<Eigen::Matrix<float, -1, 3, 1, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<unsigned int, -1, -1, 1, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&);
  117. 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::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> >&);
  118. template void igl::doublearea<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&);
  119. 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::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<int, -1, 1, 0, -1, 1> >&);
  120. 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::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, 1, 0, -1, 1> >&);
  121. 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::PlainObjectBase<Eigen::Matrix<double, 2, 1, 0, 2, 1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, 2, 1, 0, 2, 1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, 2, 1, 0, 2, 1> > const&);
  122. #endif