intrinsic_delaunay_triangulation.cpp 11 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286
  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 "intrinsic_delaunay_triangulation.h"
  9. #include "is_intrinsic_delaunay.h"
  10. #include "tan_half_angle.h"
  11. #include "unique_edge_map.h"
  12. #include "flip_edge.h"
  13. #include "EPS.h"
  14. #include "matlab_format.h"
  15. #include <iostream>
  16. #include <queue>
  17. #include <map>
  18. template <
  19. typename Derivedl_in,
  20. typename DerivedF_in,
  21. typename Derivedl,
  22. typename DerivedF>
  23. IGL_INLINE void igl::intrinsic_delaunay_triangulation(
  24. const Eigen::MatrixBase<Derivedl_in> & l_in,
  25. const Eigen::MatrixBase<DerivedF_in> & F_in,
  26. Eigen::PlainObjectBase<Derivedl> & l,
  27. Eigen::PlainObjectBase<DerivedF> & F)
  28. {
  29. typedef Eigen::Matrix<typename DerivedF::Scalar,Eigen::Dynamic,2> MatrixX2I;
  30. typedef Eigen::Matrix<typename DerivedF::Scalar,Eigen::Dynamic,1> VectorXI;
  31. MatrixX2I E,uE;
  32. VectorXI EMAP;
  33. std::vector<std::vector<typename DerivedF::Scalar> > uE2E;
  34. return intrinsic_delaunay_triangulation(l_in,F_in,l,F,E,uE,EMAP,uE2E);
  35. }
  36. template <
  37. typename Derivedl_in,
  38. typename DerivedF_in,
  39. typename Derivedl,
  40. typename DerivedF,
  41. typename DerivedE,
  42. typename DeriveduE,
  43. typename DerivedEMAP,
  44. typename uE2EType>
  45. IGL_INLINE void igl::intrinsic_delaunay_triangulation(
  46. const Eigen::MatrixBase<Derivedl_in> & l_in,
  47. const Eigen::MatrixBase<DerivedF_in> & F_in,
  48. Eigen::PlainObjectBase<Derivedl> & l,
  49. Eigen::PlainObjectBase<DerivedF> & F,
  50. Eigen::PlainObjectBase<DerivedE> & E,
  51. Eigen::PlainObjectBase<DeriveduE> & uE,
  52. Eigen::PlainObjectBase<DerivedEMAP> & EMAP,
  53. std::vector<std::vector<uE2EType> > & uE2E)
  54. {
  55. igl::unique_edge_map(F_in, E, uE, EMAP, uE2E);
  56. // We're going to work in place
  57. l = l_in;
  58. F = F_in;
  59. typedef typename DerivedF::Scalar Index;
  60. typedef typename Derivedl::Scalar Scalar;
  61. const Index num_faces = F.rows();
  62. bool any_flips = true;
  63. while(any_flips)
  64. {
  65. any_flips = false;
  66. std::vector<Index> face_queue;
  67. face_queue.reserve(32);
  68. std::vector<Index> pushed;
  69. // 32 is faster than 8
  70. pushed.reserve(32);
  71. // Does edge (a,b) exist in the edges of all faces incident on
  72. // existing unique edge uei.
  73. //
  74. // Inputs:
  75. // a 1st end-point of query edge
  76. // b 2nd end-point of query edge
  77. // uei index into uE/uE2E of unique edge
  78. // uE2E map from unique edges to half-edges (see unique_edge_map)
  79. // E #F*3 by 2 list of half-edges
  80. //
  81. const auto edge_exists_near =
  82. [&](const Index & a,const Index & b,const Index & uei)->bool
  83. {
  84. face_queue.clear();
  85. pushed.clear();
  86. assert(a!=b);
  87. // Not handling case where (a,b) is edge of face incident on uei
  88. // since this can't happen for edge-flipping.
  89. assert(a!=uE(uei,0));
  90. assert(a!=uE(uei,1));
  91. assert(b!=uE(uei,0));
  92. assert(b!=uE(uei,1));
  93. // starting with the (2) faces incident on e, consider all faces
  94. // incident on edges containing either a or b.
  95. //
  96. // face_queue Queue containing faces incident on exactly one of a/b
  97. // Using a vector seems mildly faster
  98. const Index f1 = uE2E[uei][0]%num_faces;
  99. const Index f2 = uE2E[uei][1]%num_faces;
  100. // map is faster than unordered_map here, and vector + brute force
  101. // is_member check is even faster
  102. face_queue.push_back(f1);
  103. pushed.push_back(f1);
  104. face_queue.push_back(f2);
  105. pushed.push_back(f2);
  106. while(!face_queue.empty())
  107. {
  108. const Index f = face_queue.back();
  109. face_queue.pop_back();
  110. // consider each edge of this face
  111. for(int c = 0;c<3;c++)
  112. {
  113. // Unique edge id
  114. const Index uec = EMAP(c*num_faces+f);
  115. const Index s = uE(uec,0);
  116. const Index d = uE(uec,1);
  117. const bool ona = s == a || d == a;
  118. const bool onb = s == b || d == b;
  119. // Is this the edge we're looking for?
  120. if(ona && onb)
  121. {
  122. return true;
  123. }
  124. // not incident on either?
  125. if(!ona && !onb)
  126. {
  127. continue;
  128. }
  129. // loop over all incident half-edges
  130. for(const auto & he : uE2E[uec])
  131. {
  132. // face of this he
  133. const Index fhe = he%num_faces;
  134. bool already_pushed = false;
  135. for(const auto & fp : pushed)
  136. {
  137. if(fp == fhe)
  138. {
  139. already_pushed = true;
  140. break;
  141. }
  142. }
  143. if(!already_pushed)
  144. {
  145. pushed.push_back(fhe);
  146. face_queue.push_back(fhe);
  147. }
  148. }
  149. }
  150. }
  151. return false;
  152. };
  153. // Vector is faster than queue...
  154. std::vector<Index> Q;
  155. Q.reserve(uE2E.size());
  156. for (size_t uei=0; uei<uE2E.size(); uei++)
  157. {
  158. Q.push_back(uei);
  159. }
  160. // I tried using a "delaunay_since = iter" flag to avoid duplicates, but there
  161. // was no speed up.
  162. while(!Q.empty())
  163. {
  164. const Index uei = Q.back();
  165. Q.pop_back();
  166. if (uE2E[uei].size() == 2)
  167. {
  168. if(!is_intrinsic_delaunay(l,uE2E,num_faces,uei))
  169. {
  170. // update l just before flipping edge
  171. // . //
  172. // /|\ //
  173. // a/ | \d //
  174. // / e \ //
  175. // / | \ //
  176. // .----|-f--. //
  177. // \ | / //
  178. // \ | / //
  179. // b\α|δ/c //
  180. // \|/ //
  181. // . //
  182. // Annotated from flip_edge:
  183. // Edge to flip [v1,v2] --> [v3,v4]
  184. // Before:
  185. // F(f1,:) = [v1,v2,v4] // in some cyclic order
  186. // F(f2,:) = [v1,v3,v2] // in some cyclic order
  187. // After:
  188. // F(f1,:) = [v1,v3,v4] // in *this* order
  189. // F(f2,:) = [v2,v4,v3] // in *this* order
  190. //
  191. // v1 v1
  192. // /|\ / \
  193. // c/ | \b c/f1 \b
  194. // v3 /f2|f1\ v4 => v3 /__f__\ v4
  195. // \ e / \ f2 /
  196. // d\ | /a d\ /a
  197. // \|/ \ /
  198. // v2 v2
  199. //
  200. // Compute intrinsic length of oppposite edge
  201. assert(uE2E[uei].size() == 2 && "edge should have 2 incident faces");
  202. const Index f1 = uE2E[uei][0]%num_faces;
  203. const Index f2 = uE2E[uei][1]%num_faces;
  204. const Index c1 = uE2E[uei][0]/num_faces;
  205. const Index c2 = uE2E[uei][1]/num_faces;
  206. assert(c1 < 3);
  207. assert(c2 < 3);
  208. assert(f1 != f2);
  209. const Index v1 = F(f1, (c1+1)%3);
  210. const Index v2 = F(f1, (c1+2)%3);
  211. const Index v4 = F(f1, c1);
  212. const Index v3 = F(f2, c2);
  213. assert(F(f2, (c2+2)%3) == v1);
  214. assert(F(f2, (c2+1)%3) == v2);
  215. // From gptoolbox/mesh/flip_edge.m
  216. // "If edge-after-flip already exists then this will create a non-manifold
  217. // edge"
  218. // Yes, this can happen: e.g., an edge of a tetrahedron."
  219. // "If two edges will be the same edge after flip then this will create a
  220. // non-manifold edge."
  221. // I dont' think this can happen if we flip one at a time. gptoolbox
  222. // flips in parallel.
  223. //// Over 50% of the time is spent doing this check...
  224. //bool flippable = !edge_exists_near(v3,v4,uei);
  225. //if(flippable)
  226. if(true)
  227. {
  228. any_flips = true;
  229. assert( std::abs(l(f1,c1)-l(f2,c2)) < igl::EPS<Scalar>() );
  230. const Scalar e = l(f1,c1);
  231. const Scalar a = l(f1,(c1+1)%3);
  232. const Scalar b = l(f1,(c1+2)%3);
  233. const Scalar c = l(f2,(c2+1)%3);
  234. const Scalar d = l(f2,(c2+2)%3);
  235. // tan(α/2)
  236. const Scalar tan_a_2= tan_half_angle(a,b,e);
  237. // tan(δ/2)
  238. const Scalar tan_d_2 = tan_half_angle(d,e,c);
  239. // tan((α+δ)/2)
  240. const Scalar tan_a_d_2 = (tan_a_2 + tan_d_2)/(1.0-tan_a_2*tan_d_2);
  241. // cos(α+δ)
  242. const Scalar cos_a_d =
  243. (1.0 - tan_a_d_2*tan_a_d_2)/(1.0+tan_a_d_2*tan_a_d_2);
  244. const Scalar f = sqrt(b*b + c*c - 2.0*b*c*cos_a_d);
  245. l(f1,0) = f;
  246. l(f1,1) = b;
  247. l(f1,2) = c;
  248. l(f2,0) = f;
  249. l(f2,1) = d;
  250. l(f2,2) = a;
  251. flip_edge(F, E, uE, EMAP, uE2E, uei);
  252. // append neighbors to back
  253. const size_t e_24 = f1 + ((c1 + 1) % 3) * num_faces;
  254. const size_t e_41 = f1 + ((c1 + 2) % 3) * num_faces;
  255. const size_t e_13 = f2 + ((c2 + 1) % 3) * num_faces;
  256. const size_t e_32 = f2 + ((c2 + 2) % 3) * num_faces;
  257. const size_t ue_24 = EMAP(e_24);
  258. const size_t ue_41 = EMAP(e_41);
  259. const size_t ue_13 = EMAP(e_13);
  260. const size_t ue_32 = EMAP(e_32);
  261. Q.push_back(ue_24);
  262. Q.push_back(ue_41);
  263. Q.push_back(ue_13);
  264. Q.push_back(ue_32);
  265. }
  266. }
  267. }
  268. }
  269. }
  270. }
  271. #ifdef IGL_STATIC_LIBRARY
  272. // Explicit template instantiation
  273. // generated by autoexplicit.sh
  274. template void igl::intrinsic_delaunay_triangulation<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::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, int>(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> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, std::vector<std::vector<int, std::allocator<int> >, std::allocator<std::vector<int, std::allocator<int> > > >&);
  275. // generated by autoexplicit.sh
  276. template void igl::intrinsic_delaunay_triangulation<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::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<double, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
  277. #endif