propagate_winding_numbers.cpp 13 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319
  1. // This file is part of libigl, a simple c++ geometry processing library.
  2. //
  3. // Copyright (C) 2015 Qingnan Zhou <qnzhou@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. //
  9. #include "propagate_winding_numbers.h"
  10. #include "../../extract_manifold_patches.h"
  11. #include "../../extract_non_manifold_edge_curves.h"
  12. #include "../../facet_components.h"
  13. #include "../../unique_edge_map.h"
  14. #include "../../piecewise_constant_winding_number.h"
  15. #include "../../writeOBJ.h"
  16. #include "../../writePLY.h"
  17. #include "../../get_seconds.h"
  18. #include "../../LinSpaced.h"
  19. #include "order_facets_around_edge.h"
  20. #include "outer_facet.h"
  21. #include "closest_facet.h"
  22. #include "assign.h"
  23. #include "extract_cells.h"
  24. #include "cell_adjacency.h"
  25. #include <stdexcept>
  26. #include <limits>
  27. #include <vector>
  28. #include <tuple>
  29. #include <queue>
  30. //#define PROPAGATE_WINDING_NUMBER_TIMING
  31. template<
  32. typename DerivedV,
  33. typename DerivedF,
  34. typename DerivedL,
  35. typename DerivedW>
  36. IGL_INLINE bool igl::copyleft::cgal::propagate_winding_numbers(
  37. const Eigen::PlainObjectBase<DerivedV>& V,
  38. const Eigen::PlainObjectBase<DerivedF>& F,
  39. const Eigen::PlainObjectBase<DerivedL>& labels,
  40. Eigen::PlainObjectBase<DerivedW>& W)
  41. {
  42. #ifdef PROPAGATE_WINDING_NUMBER_TIMING
  43. const auto & tictoc = []() -> double
  44. {
  45. static double t_start = igl::get_seconds();
  46. double diff = igl::get_seconds()-t_start;
  47. t_start += diff;
  48. return diff;
  49. };
  50. const auto log_time = [&](const std::string& label) -> void {
  51. std::cout << "propagate_winding_num." << label << ": "
  52. << tictoc() << std::endl;
  53. };
  54. tictoc();
  55. #endif
  56. Eigen::MatrixXi E, uE;
  57. Eigen::VectorXi EMAP;
  58. std::vector<std::vector<size_t> > uE2E;
  59. igl::unique_edge_map(F, E, uE, EMAP, uE2E);
  60. Eigen::VectorXi P;
  61. const size_t num_patches = igl::extract_manifold_patches(F, EMAP, uE2E, P);
  62. DerivedW per_patch_cells;
  63. const size_t num_cells =
  64. igl::copyleft::cgal::extract_cells(
  65. V, F, P, E, uE, uE2E, EMAP, per_patch_cells);
  66. #ifdef PROPAGATE_WINDING_NUMBER_TIMING
  67. log_time("cell_extraction");
  68. #endif
  69. return propagate_winding_numbers(V, F,
  70. uE, uE2E,
  71. num_patches, P,
  72. num_cells, per_patch_cells,
  73. labels, W);
  74. }
  75. template<
  76. typename DerivedV,
  77. typename DerivedF,
  78. typename DeriveduE,
  79. typename uE2EType,
  80. typename DerivedP,
  81. typename DerivedC,
  82. typename DerivedL,
  83. typename DerivedW>
  84. IGL_INLINE bool igl::copyleft::cgal::propagate_winding_numbers(
  85. const Eigen::PlainObjectBase<DerivedV>& V,
  86. const Eigen::PlainObjectBase<DerivedF>& F,
  87. const Eigen::PlainObjectBase<DeriveduE>& uE,
  88. const std::vector<std::vector<uE2EType> >& uE2E,
  89. const size_t num_patches,
  90. const Eigen::PlainObjectBase<DerivedP>& P,
  91. const size_t num_cells,
  92. const Eigen::PlainObjectBase<DerivedC>& C,
  93. const Eigen::PlainObjectBase<DerivedL>& labels,
  94. Eigen::PlainObjectBase<DerivedW>& W)
  95. {
  96. #ifdef PROPAGATE_WINDING_NUMBER_TIMING
  97. const auto & tictoc = []() -> double
  98. {
  99. static double t_start = igl::get_seconds();
  100. double diff = igl::get_seconds()-t_start;
  101. t_start += diff;
  102. return diff;
  103. };
  104. const auto log_time = [&](const std::string& label) -> void {
  105. std::cout << "propagate_winding_num." << label << ": "
  106. << tictoc() << std::endl;
  107. };
  108. tictoc();
  109. #endif
  110. bool valid = true;
  111. if (!piecewise_constant_winding_number(F, uE, uE2E))
  112. {
  113. assert(false && "Input mesh is not orientable");
  114. std::cerr << "Input mesh is not orientable!" << std::endl;
  115. valid = false;
  116. }
  117. const size_t num_faces = F.rows();
  118. typedef std::tuple<typename DerivedC::Scalar, bool, size_t> CellConnection;
  119. std::vector<std::set<CellConnection> > cell_adj;
  120. igl::copyleft::cgal::cell_adjacency(C, num_cells, cell_adj);
  121. #ifdef PROPAGATE_WINDING_NUMBER_TIMING
  122. log_time("cell_connectivity");
  123. #endif
  124. auto save_cell = [&](const std::string& filename, size_t cell_id) -> void{
  125. std::vector<size_t> faces;
  126. for (size_t i=0; i<num_patches; i++) {
  127. if ((C.row(i).array() == cell_id).any()) {
  128. for (size_t j=0; j<num_faces; j++) {
  129. if ((size_t)P[j] == i) {
  130. faces.push_back(j);
  131. }
  132. }
  133. }
  134. }
  135. Eigen::MatrixXi cell_faces(faces.size(), 3);
  136. for (size_t i=0; i<faces.size(); i++) {
  137. cell_faces.row(i) = F.row(faces[i]);
  138. }
  139. Eigen::MatrixXd vertices;
  140. assign(V,vertices);
  141. writePLY(filename, vertices, cell_faces);
  142. };
  143. #ifndef NDEBUG
  144. {
  145. // Check for odd cycle.
  146. Eigen::VectorXi cell_labels(num_cells);
  147. cell_labels.setZero();
  148. Eigen::VectorXi parents(num_cells);
  149. parents.setConstant(-1);
  150. auto trace_parents = [&](size_t idx) -> std::list<size_t> {
  151. std::list<size_t> path;
  152. path.push_back(idx);
  153. while ((size_t)parents[path.back()] != path.back()) {
  154. path.push_back(parents[path.back()]);
  155. }
  156. return path;
  157. };
  158. for (size_t i=0; i<num_cells; i++) {
  159. if (cell_labels[i] == 0) {
  160. cell_labels[i] = 1;
  161. std::queue<size_t> Q;
  162. Q.push(i);
  163. parents[i] = i;
  164. while (!Q.empty()) {
  165. size_t curr_idx = Q.front();
  166. Q.pop();
  167. int curr_label = cell_labels[curr_idx];
  168. for (const auto& neighbor : cell_adj[curr_idx]) {
  169. if (cell_labels[std::get<0>(neighbor)] == 0)
  170. {
  171. cell_labels[std::get<0>(neighbor)] = curr_label * -1;
  172. Q.push(std::get<0>(neighbor));
  173. parents[std::get<0>(neighbor)] = curr_idx;
  174. } else
  175. {
  176. if (cell_labels[std::get<0>(neighbor)] != curr_label * -1)
  177. {
  178. std::cerr << "Odd cell cycle detected!" << std::endl;
  179. auto path = trace_parents(curr_idx);
  180. path.reverse();
  181. auto path2 = trace_parents(std::get<0>(neighbor));
  182. path.insert(path.end(), path2.begin(), path2.end());
  183. for (auto cell_id : path)
  184. {
  185. std::cout << cell_id << " ";
  186. std::stringstream filename;
  187. filename << "cell_" << cell_id << ".ply";
  188. save_cell(filename.str(), cell_id);
  189. }
  190. std::cout << std::endl;
  191. valid = false;
  192. }
  193. // Do not fail when odd cycle is detected because the resulting
  194. // integer winding number field, although inconsistent, may still
  195. // be used if the problem region is local and embedded within a
  196. // valid volume.
  197. //assert(cell_labels[std::get<0>(neighbor)] == curr_label * -1);
  198. }
  199. }
  200. }
  201. }
  202. }
  203. #ifdef PROPAGATE_WINDING_NUMBER_TIMING
  204. log_time("odd_cycle_check");
  205. #endif
  206. }
  207. #endif
  208. size_t outer_facet;
  209. bool flipped;
  210. Eigen::VectorXi I = igl::LinSpaced<Eigen::VectorXi>(num_faces, 0, num_faces-1);
  211. igl::copyleft::cgal::outer_facet(V, F, I, outer_facet, flipped);
  212. #ifdef PROPAGATE_WINDING_NUMBER_TIMING
  213. log_time("outer_facet");
  214. #endif
  215. const size_t outer_patch = P[outer_facet];
  216. const size_t infinity_cell = C(outer_patch, flipped?1:0);
  217. Eigen::VectorXi patch_labels(num_patches);
  218. const int INVALID = std::numeric_limits<int>::max();
  219. patch_labels.setConstant(INVALID);
  220. for (size_t i=0; i<num_faces; i++) {
  221. if (patch_labels[P[i]] == INVALID) {
  222. patch_labels[P[i]] = labels[i];
  223. } else {
  224. assert(patch_labels[P[i]] == labels[i]);
  225. }
  226. }
  227. assert((patch_labels.array() != INVALID).all());
  228. const size_t num_labels = patch_labels.maxCoeff()+1;
  229. Eigen::MatrixXi per_cell_W(num_cells, num_labels);
  230. per_cell_W.setConstant(INVALID);
  231. per_cell_W.row(infinity_cell).setZero();
  232. std::queue<size_t> Q;
  233. Q.push(infinity_cell);
  234. while (!Q.empty()) {
  235. size_t curr_cell = Q.front();
  236. Q.pop();
  237. for (const auto& neighbor : cell_adj[curr_cell]) {
  238. size_t neighbor_cell, patch_idx;
  239. bool direction;
  240. std::tie(neighbor_cell, direction, patch_idx) = neighbor;
  241. if ((per_cell_W.row(neighbor_cell).array() == INVALID).any()) {
  242. per_cell_W.row(neighbor_cell) = per_cell_W.row(curr_cell);
  243. for (size_t i=0; i<num_labels; i++) {
  244. int inc = (patch_labels[patch_idx] == (int)i) ?
  245. (direction ? -1:1) :0;
  246. per_cell_W(neighbor_cell, i) =
  247. per_cell_W(curr_cell, i) + inc;
  248. }
  249. Q.push(neighbor_cell);
  250. } else {
  251. #ifndef NDEBUG
  252. // Checking for winding number consistency.
  253. // This check would inevitably fail for meshes that contain open
  254. // boundary or non-orientable. However, the inconsistent winding number
  255. // field would still be useful in some cases such as when problem region
  256. // is local and embedded within the volume. This, unfortunately, is the
  257. // best we can do because the problem of computing integer winding
  258. // number is ill-defined for open and non-orientable surfaces.
  259. for (size_t i=0; i<num_labels; i++) {
  260. if ((int)i == patch_labels[patch_idx]) {
  261. int inc = direction ? -1:1;
  262. //assert(per_cell_W(neighbor_cell, i) ==
  263. // per_cell_W(curr_cell, i) + inc);
  264. } else {
  265. //assert(per_cell_W(neighbor_cell, i) ==
  266. // per_cell_W(curr_cell, i));
  267. }
  268. }
  269. #endif
  270. }
  271. }
  272. }
  273. #ifdef PROPAGATE_WINDING_NUMBER_TIMING
  274. log_time("propagate_winding_number");
  275. #endif
  276. W.resize(num_faces, num_labels*2);
  277. for (size_t i=0; i<num_faces; i++)
  278. {
  279. const size_t patch = P[i];
  280. const size_t positive_cell = C(patch, 0);
  281. const size_t negative_cell = C(patch, 1);
  282. for (size_t j=0; j<num_labels; j++) {
  283. W(i,j*2 ) = per_cell_W(positive_cell, j);
  284. W(i,j*2+1) = per_cell_W(negative_cell, j);
  285. }
  286. }
  287. #ifdef PROPAGATE_WINDING_NUMBER_TIMING
  288. log_time("store_result");
  289. #endif
  290. return valid;
  291. }
  292. #ifdef IGL_STATIC_LIBRARY
  293. // Explicit template instantiation
  294. // generated by autoexplicit.sh
  295. template bool igl::copyleft::cgal::propagate_winding_numbers<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, -1, 1, -1, -1>, Eigen::Matrix<int, -1, 3, 1, -1, 3>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, unsigned long, 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> >(Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, -1, 1, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 1, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, std::vector<std::vector<unsigned long, std::allocator<unsigned long> >, std::allocator<std::vector<unsigned long, std::allocator<unsigned long> > > > const&, unsigned long, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, unsigned long, Eigen::PlainObjectBase<Eigen::Matrix<int, -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> >&);
  296. // generated by autoexplicit.sh
  297. template bool igl::copyleft::cgal::propagate_winding_numbers<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, -1, 1, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, unsigned long, 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> >(Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, -1, 1, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, std::vector<std::vector<unsigned long, std::allocator<unsigned long> >, std::allocator<std::vector<unsigned long, std::allocator<unsigned long> > > > const&, unsigned long, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, unsigned long, Eigen::PlainObjectBase<Eigen::Matrix<int, -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> >&);
  298. template bool igl::copyleft::cgal::propagate_winding_numbers<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -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::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -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> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
  299. #endif