main.cpp 23 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720721
  1. #include <igl/readOFF.h>
  2. #include <igl/readOBJ.h>
  3. #include <igl/n_polyvector.h>
  4. #include <igl/integrable_polyvector_fields.h>
  5. #include <igl/viewer/Viewer.h>
  6. #include <igl/local_basis.h>
  7. #include <igl/avg_edge_length.h>
  8. #include <igl/is_border_vertex.h>
  9. #include <igl/adjacency_list.h>
  10. #include <igl/vertex_triangle_adjacency.h>
  11. #include <igl/triangle_triangle_adjacency.h>
  12. #include <igl/edge_topology.h>
  13. #include <igl/jet.h>
  14. #include <igl/barycenter.h>
  15. #include <igl/polyvector_field_matchings.h>
  16. #include <igl/polyvector_field_singularities_from_matchings.h>
  17. #include <igl/polyvector_field_cut_mesh_with_singularities.h>
  18. #include <igl/polyvector_field_comb_from_matchings_and_cuts.h>
  19. #include <igl/polyvector_field_poisson_reconstruction.h>
  20. #include <igl/cut_mesh.h>
  21. #include <igl/slice.h>
  22. #include <igl/false_barycentric_subdivision.h>
  23. #include <iostream>
  24. #include <fstream>
  25. #include <igl/matlab_format.h>
  26. #include "tutorial_shared_path.h"
  27. using namespace std;
  28. // Input mesh
  29. Eigen::MatrixXd V;
  30. Eigen::MatrixXi F;
  31. std::vector<bool> V_border;
  32. std::vector<std::vector<int> > VF, VFi;
  33. std::vector<std::vector<int> > VV;
  34. Eigen::MatrixXi TT, TTi;
  35. Eigen::MatrixXi E, E2F, F2E;
  36. // Per face bases (only needed to generate constraints)
  37. Eigen::MatrixXd B1,B2,B3;
  38. // "Subdivided" mesh obtained by splitting each triangle into 3 (only needed for display)
  39. Eigen::MatrixXd Vbs;
  40. Eigen::MatrixXi Fbs;
  41. // Scale for visualizing the fields
  42. double global_scale;
  43. // Scale for visualizing textures
  44. double uv_scale;
  45. // Data for original PolyVector field
  46. Eigen::MatrixXd two_pv_ori; // field
  47. Eigen::VectorXi singularities_ori; // singularities
  48. Eigen::VectorXd curl_ori; // curl per edge
  49. Eigen::MatrixXi cuts_ori; // cut edges
  50. Eigen::MatrixXd two_pv_poisson_ori; // field after poisson integration
  51. Eigen::VectorXf poisson_error_ori; // poisson integration error
  52. Eigen::MatrixXd scalars_ori;
  53. Eigen::MatrixXd Vcut_ori;
  54. Eigen::MatrixXi Fcut_ori;
  55. // Data for curl-free PolyVector field
  56. Eigen::MatrixXd two_pv; // field
  57. Eigen::VectorXi singularities; // singularities
  58. Eigen::VectorXd curl; // curl per edge
  59. Eigen::MatrixXi cuts; // cut edges
  60. Eigen::MatrixXd two_pv_poisson; // field after poisson integration
  61. Eigen::VectorXf poisson_error; // poisson integration error
  62. Eigen::MatrixXd scalars;
  63. Eigen::MatrixXd Vcut;
  64. Eigen::MatrixXi Fcut;
  65. // Vector of constrained faces
  66. Eigen::VectorXi b;
  67. // Matrix of constraints
  68. Eigen::MatrixXd bc;
  69. // "constraint level" flag (level=2 indicates that both directions are constrained,
  70. // level = 1 indicates a partially constrained face, i.e. only the first vector will
  71. // be constrained)
  72. Eigen::VectorXi blevel;
  73. // Face Barycenters (only needed for display)
  74. Eigen::MatrixXd B;
  75. // percentage of constrained faces
  76. double constraint_percentage = 0.002;
  77. // Random length factor
  78. double rand_factor = 5;
  79. // The set of parameters for calculating the curl-free fields
  80. igl::integrable_polyvector_fields_parameters params;
  81. // Solver data (needed for precomputation)
  82. igl::IntegrableFieldSolverData<Eigen::MatrixXd, Eigen::MatrixXi, Eigen::MatrixXd, Eigen::MatrixXd> ipfdata;
  83. //texture image
  84. Eigen::Matrix<unsigned char,Eigen::Dynamic,Eigen::Dynamic> texture_R, texture_G, texture_B;
  85. int display_mode = 1;
  86. int iter = 0;
  87. // Create a texture that hides the integer translation in the parametrization
  88. void line_texture(Eigen::Matrix<unsigned char,Eigen::Dynamic,Eigen::Dynamic> &texture_R,
  89. Eigen::Matrix<unsigned char,Eigen::Dynamic,Eigen::Dynamic> &texture_G,
  90. Eigen::Matrix<unsigned char,Eigen::Dynamic,Eigen::Dynamic> &texture_B)
  91. {
  92. unsigned size = 128;
  93. unsigned size2 = size/2;
  94. unsigned lineWidth = 3;
  95. texture_R.setConstant(size, size, 255);
  96. for (unsigned i=0; i<size; ++i)
  97. for (unsigned j=size2-lineWidth; j<=size2+lineWidth; ++j)
  98. texture_R(i,j) = 0;
  99. for (unsigned i=size2-lineWidth; i<=size2+lineWidth; ++i)
  100. for (unsigned j=0; j<size; ++j)
  101. texture_R(i,j) = 0;
  102. texture_G = texture_R;
  103. texture_B = texture_R;
  104. }
  105. // Create a random set of tangent vectors
  106. void generate_constraints()
  107. {
  108. b.resize(42,1);
  109. b<<
  110. 663,
  111. 513,
  112. 3872,
  113. 2601,
  114. 3549,
  115. 2721,
  116. 3796,
  117. 594,
  118. 868,
  119. 1730,
  120. 1581,
  121. 3081,
  122. 1471,
  123. 1650,
  124. 454,
  125. 2740,
  126. 2945,
  127. 3808,
  128. 3679,
  129. 3589,
  130. 450,
  131. 2656,
  132. 1791,
  133. 1792,
  134. 2917,
  135. 3744,
  136. 1536,
  137. 2809,
  138. 3866,
  139. 3658,
  140. 1665,
  141. 2670,
  142. 1601,
  143. 1793,
  144. 3614,
  145. 524,
  146. 2877,
  147. 449,
  148. 455,
  149. 3867,
  150. 3871,
  151. 2592;
  152. blevel.setOnes(b.rows(),1);
  153. bc.resize(b.rows(),6);
  154. bc<<
  155. -0.88046298335147721,0.27309862654264377,0.38755912468723769,-0.350632259447135,-0.92528970817792766,-0.14455440005410564,
  156. 0.91471470003012889,0.392936119054695,-0.094330397492144599,0.32234487030777614,-0.85027369799342767,-0.41608703787410195,
  157. 0.94335566040683105,0.073867667925654024,-0.32345581709658111,0.19950360079371404,-0.90525435056476755,0.37511714710727789,
  158. -0.92054671613540229,0.15077598183983737,0.36036141124232496,-0.27998315313687211,-0.89796618385425386,-0.33950871360506074,
  159. 0.88944399663239937,0.23035525634795684,-0.39474780902172396,0.27297422303141039,-0.96047177712172194,0.054580572670497061,
  160. -0.83112706922096102,-0.55599928943162547,0.0096221078617792517,0.52546831822766438,-0.79091522174894457,-0.31358596675362826,
  161. 0.90724658517664569,-0.41046292080872998,-0.091781394228251156,-0.34252813327252363,-0.84767620917196618,0.40511667741613094,
  162. -0.8932101465465786,0.23975524191487588,0.38038540729184012,-0.33645713296414853,-0.91759364410558497,-0.21170380718016926,
  163. -0.87839308390284521,0.27039404931731387,0.39409725734320344,-0.29712518405497651,-0.95481177255558192,-0.0071487054467167244,
  164. -0.91448048788760505,-0.17055891298176448,0.36692655188106316,0.29811257890714044,-0.89715315396744022,0.32595261714489804,
  165. 0.82798126471567035,-0.56074230404745851,0.003885065171440813,-0.53510484459763941,-0.78801608401899037,0.30445600111594384,
  166. -0.87831929581593793,0.25312706437601257,0.40556368658667746,-0.26531767440854004,-0.9637845762158106,0.026941089342378936,
  167. -0.87482003689209031,-0.27011021313654948,0.4021571531272935,0.32303198334357713,-0.94388894366288889,0.0687313594225408,
  168. 0.87408456883093666,-0.48487387939766946,-0.029554823793924323,-0.43846604347950752,-0.81368808449189478,0.38165328489525413,
  169. 0.94988212941972827,-0.041936960956176939,-0.30978255521381903,-0.16637246118261648,-0.90677959514398765,-0.3873899456497869,
  170. 0.87516493768857484,-0.27181042881473483,-0.40025669591913515,-0.36755520380602424,-0.91147911093961553,-0.18468622708756641,
  171. -0.87064244687577641,0.27922257819020818,0.40498948323008854,-0.32176729617260508,-0.94599403842079244,-0.039510585747255161,
  172. -0.91274615133859638,-0.1832657459904497,0.36511385835536858,0.29782083933521708,-0.91026141603074595,0.28762284704690655,
  173. 0.875611546674125,0.28258715176515403,-0.39172556846369444,0.36000975242683186,-0.92250014843287764,0.13923524804764778,
  174. 0.76763693171870195,-0.64088483679642994,0.00040868803559811201,-0.63058113310112196,-0.75518119878562417,0.17907761327874747,
  175. 0.963064265211517,0.17044712473620016,-0.20845862597111031,0.061832174999749308,-0.89345471128813481,-0.44487690546019126,
  176. -0.88228998376691692,-0.46837234310148523,0.046815945597605227,0.41604986062280985,-0.82249303168905052,-0.38782434980116298,
  177. -0.96608602970701829,0.11121907649833783,0.23304098400879364,0.010641270548624404,-0.88457418950525291,0.46627810008860171,
  178. -0.96329451047686887,0.055809409647239343,0.26258140810033831,0.07182051046944142,-0.88891411988025926,0.45240855623364673,
  179. -0.71244584326772997,-0.70122065397026967,-0.026655484539588895,0.70046172163981768,-0.70836773631021255,-0.086997279682342638,
  180. 0.88646445996853696,0.2549240118236365,-0.38625705094979518,0.35132981358631576,-0.91395520354543514,0.20310895597591658,
  181. -0.86109327343809683,-0.30822574449366841,0.40437020769461601,0.37896596246993836,-0.91928725525816557,0.10628142645421024,
  182. 0.86443027504389158,-0.29669958642983363,-0.40586901212079213,-0.37200560813855077,-0.92052106924988175,-0.11938504337027039,
  183. 0.95370728000967508,-0.24689991217686594,-0.17170572915195079,-0.14736898596800915,-0.88138264597997584,0.4488284898935197,
  184. -0.81439393313167019,0.57995723960933832,0.020300785774083896,-0.55494919604589421,-0.78855235001798585,0.26498411478639117,
  185. 0.89527216270596455,0.22395367264061938,-0.38513959442592183,0.33680943342191538,-0.90609008785063272,0.25604717974787594,
  186. -0.9003647006267198,0.20802946062196581,0.38218732236782926,-0.32431023000528064,-0.90640636884236436,-0.27064805418831556,
  187. -0.87050937437709508,-0.28614105672408718,0.40042068475344922,0.37746788793940733,-0.91025870352880611,0.17013843253251126,
  188. -0.95715079751532439,0.0030851788865496879,0.28957353554324744,0.12908381923211401,-0.89056292562302997,0.43615942397041058,
  189. -0.87324677619075319,-0.28591869514051466,0.39457644080913162,0.3438918663433696,-0.93530416305293196,0.083333707698197687,
  190. 0.91999856277124803,-0.1621255206103257,-0.35681642348085474,-0.27672206872177485,-0.91342693749618353,-0.2984562389005877,
  191. -0.8669467282645521,0.29036174243712859,0.40508447128995645,-0.34873789125620602,-0.93406639205959408,-0.07682355385522964,
  192. -0.9570365266718136,-0.22821899053183647,0.17887755302603078,0.12590409644663494,-0.88275887883510706,-0.45264215483728532,
  193. -0.94033215083998489,0.087395510869996196,0.32884262311388451,-0.2131320783418921,-0.90465024471116184,-0.36902933748646671,
  194. -0.96131014054749453,0.18866284908038999,0.20072155603578434,-0.08260532909072589,-0.89255302833360861,0.44331191188407898,
  195. -0.95240414686152941,-0.02752900142620229,0.30359264668538755,0.15128346580527452,-0.9073021943457209,0.39232134929083828,
  196. -0.94070423353276911,-0.31552769387286655,0.12457053990729766,0.22959741970407915,-0.86253407908715607,-0.45091017650802745;
  197. }
  198. void drawCuts(igl::viewer::Viewer& viewer,
  199. const Eigen::MatrixXi &cuts)
  200. {
  201. int maxCutNum = cuts.sum();
  202. Eigen::MatrixXd start(maxCutNum,3);
  203. Eigen::MatrixXd end(maxCutNum,3);
  204. int ind = 0;
  205. for (unsigned int i=0;i<F.rows();i++)
  206. for (int j=0;j<3;j++)
  207. if (cuts(i,j))
  208. {
  209. start.row(ind) = V.row(F(i,j));
  210. end.row(ind) = V.row(F(i,(j+1)%3));
  211. ind++;
  212. }
  213. viewer.data.add_edges(start, end , Eigen::RowVector3d(1.,0,1.));
  214. }
  215. void drawField(igl::viewer::Viewer &viewer,
  216. const Eigen::MatrixXd &field,
  217. const Eigen::RowVector3d &color)
  218. {
  219. for (int n=0; n<2; ++n)
  220. {
  221. Eigen::MatrixXd VF = field.block(0,n*3,F.rows(),3);
  222. Eigen::VectorXd c = VF.rowwise().norm();
  223. viewer.data.add_edges(B - global_scale*VF, B + global_scale*VF , color);
  224. }
  225. }
  226. void drawConstraints(igl::viewer::Viewer &viewer)
  227. {
  228. for (int n=0; n<2; ++n)
  229. {
  230. Eigen::MatrixXd Bc = igl::slice(B, b, 1);
  231. Eigen::MatrixXd color;
  232. color.setZero(b.rows(),3);
  233. color.col(2).setOnes();
  234. for (int i =0; i<b.rows(); ++i)
  235. if (blevel[i] ==1 && n>0)
  236. color.row(i)<<0.7,0.7,0.7;
  237. // Eigen::RowVector3d color; color<<0.5,0.5,0.5;
  238. viewer.data.add_edges(Bc - global_scale*bc.block(0,n*3,bc.rows(),3), Bc + global_scale*bc.block(0,n*3,bc.rows(),3) , color);
  239. }
  240. }
  241. void colorEdgeMeshFaces(const Eigen::VectorXd &values,
  242. const double &minimum,
  243. const double &maximum,
  244. Eigen::MatrixXd &C)
  245. {
  246. C.setConstant(Fbs.rows(),3,1);
  247. Eigen::MatrixXd colors;
  248. igl::jet(values, minimum, maximum, colors);
  249. for (int ei = 0; ei<E.rows(); ++ei)
  250. {
  251. const Eigen::RowVector3d &this_color = colors.row(ei);
  252. int f0 = E2F(ei,0);
  253. int f1 = E2F(ei,1);
  254. if(f0 != -1)
  255. {
  256. int i0 = -1;
  257. for (int k = 0; k<3; k++)
  258. if (F2E(f0,k)== ei)
  259. {
  260. i0 = k;
  261. break;
  262. }
  263. C.row(3*f0+i0) = this_color;
  264. }
  265. if(f1 != -1)
  266. {
  267. int i1 = -1;
  268. for (int k = 0; k<3; k++)
  269. if (F2E(f1,k)== ei)
  270. {
  271. i1 = k;
  272. break;
  273. }
  274. C.row(3*f1+i1) = this_color;
  275. }
  276. }
  277. }
  278. void update_display(igl::viewer::Viewer& viewer)
  279. {
  280. using namespace std;
  281. using namespace Eigen;
  282. viewer.data.clear();
  283. viewer.data.lines.resize(0,9);
  284. viewer.data.points.resize(0,6);
  285. viewer.core.show_texture = false;
  286. if (display_mode == 1)
  287. {
  288. cerr<< "Displaying original field, its singularities and its cuts" <<endl;
  289. viewer.data.set_mesh(V, F);
  290. // Highlight in red the constrained faces
  291. MatrixXd C = MatrixXd::Constant(F.rows(),3,1);
  292. for (unsigned i=0; i<b.size();++i)
  293. C.row(b(i)) << 1, 0, 0;
  294. viewer.data.set_colors(C);
  295. //Draw constraints
  296. drawConstraints(viewer);
  297. // Draw Field
  298. Eigen::RowVector3d color; color<<0,0,1;
  299. drawField(viewer,two_pv_ori,color);
  300. // Draw Cuts
  301. drawCuts(viewer,cuts_ori);
  302. //Draw Singularities
  303. Eigen::MatrixXd singular_points = igl::slice(V, singularities_ori, 1);
  304. viewer.data.add_points(singular_points,Eigen::RowVector3d(239./255.,205./255.,57./255.));
  305. }
  306. if (display_mode == 2)
  307. {
  308. cerr<< "Displaying current field, its singularities and its cuts" <<endl;
  309. viewer.data.set_mesh(V, F);
  310. // Highlight in red the constrained faces
  311. MatrixXd C = MatrixXd::Constant(F.rows(),3,1);
  312. for (unsigned i=0; i<b.size();++i)
  313. C.row(b(i)) << 1, 0, 0;
  314. viewer.data.set_colors(C);
  315. //Draw constraints
  316. drawConstraints(viewer);
  317. // Draw Field
  318. Eigen::RowVector3d color; color<<0,0,1;
  319. drawField(viewer,two_pv,color);
  320. // Draw Cuts
  321. drawCuts(viewer,cuts);
  322. //Draw Singularities
  323. Eigen::MatrixXd singular_points = igl::slice(V, singularities, 1);
  324. viewer.data.add_points(singular_points,Eigen::RowVector3d(239./255.,205./255.,57./255.));
  325. }
  326. if (display_mode == 3)
  327. {
  328. cerr<< "Displaying original field and its curl" <<endl;
  329. viewer.data.set_mesh(Vbs, Fbs);
  330. Eigen::MatrixXd C;
  331. colorEdgeMeshFaces(curl_ori, 0, 0.2, C);
  332. viewer.data.set_colors(C);
  333. // Draw Field
  334. Eigen::RowVector3d color; color<<1,1,1;
  335. drawField(viewer,two_pv_ori,color);
  336. }
  337. if (display_mode == 4)
  338. {
  339. cerr<< "Displaying current field and its curl" <<endl;
  340. viewer.data.set_mesh(Vbs, Fbs);
  341. Eigen::MatrixXd C;
  342. colorEdgeMeshFaces(curl, 0, 0.2, C);
  343. viewer.data.set_colors(C);
  344. // Draw Field
  345. Eigen::RowVector3d color; color<<1,1,1;
  346. drawField(viewer,two_pv,color);
  347. }
  348. if (display_mode == 5)
  349. {
  350. cerr<< "Displaying original poisson-integrated field and original poisson error" <<endl;
  351. viewer.data.set_mesh(V, F);
  352. Eigen::MatrixXd C;
  353. igl::jet(poisson_error_ori, 0, 0.5, C);
  354. viewer.data.set_colors(C);
  355. // Draw Field
  356. Eigen::RowVector3d color; color<<1,1,1;
  357. drawField(viewer,two_pv_poisson_ori,color);
  358. }
  359. if (display_mode == 6)
  360. {
  361. cerr<< "Displaying current poisson-integrated field and current poisson error" <<endl;
  362. viewer.data.set_mesh(V, F);
  363. Eigen::MatrixXd C;
  364. igl::jet(poisson_error, 0, 0.5, C);
  365. viewer.data.set_colors(C);
  366. // Draw Field
  367. Eigen::RowVector3d color; color<<1,1,1;
  368. drawField(viewer,two_pv_poisson,color);
  369. }
  370. if (display_mode == 7)
  371. {
  372. cerr<< "Displaying original texture with cuts and singularities" <<endl;
  373. viewer.data.set_mesh(V, F);
  374. MatrixXd C = MatrixXd::Constant(F.rows(),3,1);
  375. viewer.data.set_colors(C);
  376. viewer.data.set_uv(uv_scale*scalars_ori, Fcut_ori);
  377. viewer.data.set_texture(texture_R, texture_B, texture_G);
  378. viewer.core.show_texture = true;
  379. // Draw Cuts
  380. drawCuts(viewer,cuts_ori);
  381. //Draw Singularities
  382. Eigen::MatrixXd singular_points = igl::slice(V, singularities_ori, 1);
  383. viewer.data.add_points(singular_points,Eigen::RowVector3d(239./255.,205./255.,57./255.));
  384. }
  385. if (display_mode == 8)
  386. {
  387. cerr<< "Displaying current texture with cuts and singularities" <<endl;
  388. viewer.data.set_mesh(V, F);
  389. MatrixXd C = MatrixXd::Constant(F.rows(),3,1);
  390. viewer.data.set_colors(C);
  391. viewer.data.set_uv(uv_scale*scalars, Fcut);
  392. viewer.data.set_texture(texture_R, texture_B, texture_G);
  393. viewer.core.show_texture = true;
  394. // Draw Cuts
  395. drawCuts(viewer,cuts);
  396. //Draw Singularities
  397. Eigen::MatrixXd singular_points = igl::slice(V, singularities, 1);
  398. viewer.data.add_points(singular_points,Eigen::RowVector3d(239./255.,205./255.,57./255.));
  399. }
  400. if (display_mode == 9)
  401. {
  402. cerr<< "Displaying original field overlayed onto the current integrated field" <<endl;
  403. viewer.data.set_mesh(V, F);
  404. // Highlight in red the constrained faces
  405. MatrixXd C = MatrixXd::Constant(F.rows(),3,1);
  406. for (unsigned i=0; i<b.size();++i)
  407. C.row(b(i)) << 1, 0, 0;
  408. viewer.data.set_colors(C);
  409. // Draw Field
  410. Eigen::RowVector3d color; color<<0,0,1;
  411. drawField(viewer,two_pv_ori,color);
  412. // Draw Integrated Field
  413. color<<.2,.2,.2;
  414. drawField(viewer,two_pv_poisson_ori,color);
  415. }
  416. if (display_mode == 0)
  417. {
  418. cerr<< "Displaying current field overlayed onto the current integrated field" <<endl;
  419. viewer.data.set_mesh(V, F);
  420. // Highlight in red the constrained faces
  421. MatrixXd C = MatrixXd::Constant(F.rows(),3,1);
  422. for (unsigned i=0; i<b.size();++i)
  423. C.row(b(i)) << 1, 0, 0;
  424. viewer.data.set_colors(C);
  425. // Draw Field
  426. Eigen::RowVector3d color; color<<0,0,1;
  427. drawField(viewer,two_pv,color);
  428. // Draw Integrated Field
  429. color<<.2,.2,.2;
  430. drawField(viewer,two_pv_poisson,color);
  431. }
  432. }
  433. bool key_down(igl::viewer::Viewer& viewer, unsigned char key, int modifier)
  434. {
  435. if (key == '1')
  436. {
  437. display_mode = 1;
  438. update_display(viewer);
  439. }
  440. if (key == '2')
  441. {
  442. display_mode = 2;
  443. update_display(viewer);
  444. }
  445. if (key == '3')
  446. {
  447. display_mode = 3;
  448. update_display(viewer);
  449. }
  450. if (key == '4')
  451. {
  452. display_mode = 4;
  453. update_display(viewer);
  454. }
  455. if (key == '5')
  456. {
  457. display_mode = 5;
  458. update_display(viewer);
  459. }
  460. if (key == '6')
  461. {
  462. display_mode = 6;
  463. update_display(viewer);
  464. }
  465. if (key == '7')
  466. {
  467. display_mode = 7;
  468. update_display(viewer);
  469. }
  470. if (key == '8')
  471. {
  472. display_mode = 8;
  473. update_display(viewer);
  474. }
  475. if (key == '9')
  476. {
  477. display_mode = 9;
  478. update_display(viewer);
  479. }
  480. if (key == '0')
  481. {
  482. display_mode = 0;
  483. update_display(viewer);
  484. }
  485. if (key == 'A')
  486. {
  487. //do a batch of iterations
  488. printf("--Improving Curl--\n");
  489. for (int bi = 0; bi<5; ++bi)
  490. {
  491. printf("\n\n **** Batch %d ****\n", iter);
  492. igl::integrable_polyvector_fields_solve(ipfdata, params, two_pv, iter ==0);
  493. iter++;
  494. params.wSmooth *= params.redFactor_wsmooth;
  495. }
  496. // Post process current field
  497. // Compute curl_minimizing matchings and curl
  498. printf("--Matchings and curl--\n");
  499. Eigen::MatrixXi match_ab, match_ba; // matchings across interior edges
  500. double avgCurl = igl::polyvector_field_matchings(two_pv, V, F, true, match_ab, match_ba, curl);
  501. double maxCurl = curl.maxCoeff();
  502. printf("curl -- max: %.5g, avg: %.5g\n", maxCurl, avgCurl);
  503. // Compute singularities
  504. printf("--Singularities--\n");
  505. igl::polyvector_field_singularities_from_matchings(V, F, match_ab, match_ba, singularities);
  506. printf("#singularities: %ld\n", singularities.rows());
  507. // Get mesh cuts based on singularities
  508. printf("--Cuts--\n");
  509. igl::polyvector_field_cut_mesh_with_singularities(V, F, singularities, cuts);
  510. // Comb field
  511. printf("--Combing--\n");
  512. Eigen::MatrixXd combed;
  513. igl::polyvector_field_comb_from_matchings_and_cuts(V, F, two_pv, match_ab, match_ba, cuts, combed);
  514. // Reconstruct integrable vector fields from combed field
  515. printf("--Cut mesh--\n");
  516. igl::cut_mesh(V, F, cuts, Vcut, Fcut);
  517. printf("--Poisson--\n");
  518. double avgPoisson = igl::polyvector_field_poisson_reconstruction(Vcut, Fcut, combed, scalars, two_pv_poisson, poisson_error);
  519. double maxPoisson = poisson_error.maxCoeff();
  520. printf("poisson error -- max: %.5g, avg: %.5g\n", maxPoisson, avgPoisson);
  521. update_display(viewer);
  522. }
  523. return false;
  524. }
  525. int main(int argc, char *argv[])
  526. {
  527. // Load a mesh
  528. igl::readOBJ(TUTORIAL_SHARED_PATH "/inspired_mesh.obj", V, F);
  529. printf("--Initialization--\n");
  530. V_border = igl::is_border_vertex(V,F);
  531. igl::adjacency_list(F, VV);
  532. igl::vertex_triangle_adjacency(V,F,VF,VFi);
  533. igl::triangle_triangle_adjacency(V,F,TT,TTi);
  534. igl::edge_topology(V,F,E,F2E,E2F);
  535. // Generate "subdivided" mesh for visualization of curl terms
  536. igl::false_barycentric_subdivision(V, F, Vbs, Fbs);
  537. // Compute scale for visualizing fields
  538. global_scale = .2*igl::avg_edge_length(V, F);
  539. //Compute scale for visualizing texture
  540. uv_scale = 0.6/igl::avg_edge_length(V, F);
  541. // Compute face barycenters
  542. igl::barycenter(V, F, B);
  543. // Compute local basis for faces
  544. igl::local_basis(V,F,B1,B2,B3);
  545. //Generate random vectors for constraints
  546. generate_constraints();
  547. // Interpolate a 2-PolyVector field to be used as the original field
  548. printf("--Initial solution--\n");
  549. igl::n_polyvector(V, F, b, bc, two_pv_ori);
  550. // Post process original field
  551. // Compute curl_minimizing matchings and curl
  552. Eigen::MatrixXi match_ab, match_ba; // matchings across interior edges
  553. printf("--Matchings and curl--\n");
  554. double avgCurl = igl::polyvector_field_matchings(two_pv_ori, V, F, true, match_ab, match_ba, curl_ori);
  555. double maxCurl = curl_ori.maxCoeff();
  556. printf("original curl -- max: %.5g, avg: %.5g\n", maxCurl, avgCurl);
  557. printf("--Singularities--\n");
  558. // Compute singularities
  559. igl::polyvector_field_singularities_from_matchings(V, F, V_border, VF, TT, E2F, F2E, match_ab, match_ba, singularities_ori);
  560. printf("original #singularities: %ld\n", singularities.rows());
  561. printf("--Cuts--\n");
  562. // Get mesh cuts based on singularities
  563. igl::polyvector_field_cut_mesh_with_singularities(V, F, VF, VV, TT, TTi, singularities_ori, cuts_ori);
  564. printf("--Combing--\n");
  565. // Comb field
  566. Eigen::MatrixXd combed;
  567. igl::polyvector_field_comb_from_matchings_and_cuts(V, F, TT, E2F, F2E, two_pv_ori, match_ab, match_ba, cuts_ori, combed);
  568. printf("--Cut mesh--\n");
  569. // Reconstruct integrable vector fields from combed field
  570. igl::cut_mesh(V, F, VF, VFi, TT, TTi, V_border, cuts_ori, Vcut_ori, Fcut_ori);
  571. printf("--Poisson--\n");
  572. double avgPoisson = igl::polyvector_field_poisson_reconstruction(Vcut_ori, Fcut_ori, combed, scalars_ori, two_pv_poisson_ori, poisson_error_ori);
  573. double maxPoisson = poisson_error_ori.maxCoeff();
  574. printf("poisson error -- max: %.5g, avg: %.5g\n", maxPoisson, avgPoisson);
  575. // Set the curl-free 2-PolyVector to equal the original field
  576. two_pv = two_pv_ori;
  577. singularities = singularities_ori;
  578. curl = curl_ori;
  579. cuts = cuts_ori;
  580. two_pv_poisson = two_pv_poisson_ori;
  581. poisson_error = poisson_error_ori;
  582. Vcut = Vcut_ori;
  583. Fcut = Fcut_ori;
  584. scalars = scalars_ori;
  585. printf("--Integrable - Precomputation--\n");
  586. // Precompute stuff for solver
  587. igl::integrable_polyvector_fields_precompute(V, F, b, bc, blevel, two_pv_ori, ipfdata);
  588. cerr<<"Done. Press keys 1-0 for various visualizations, 'A' to improve integrability." <<endl;
  589. igl::viewer::Viewer viewer;
  590. viewer.callback_key_down = &key_down;
  591. viewer.core.show_lines = false;
  592. key_down(viewer,'2',0);
  593. // Replace the standard texture with an integer shift invariant texture
  594. line_texture(texture_R, texture_G, texture_B);
  595. viewer.launch();
  596. return 0;
  597. }