cotmatrix_entries.cpp 4.8 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154
  1. #include <test_common.h>
  2. #include <igl/PI.h>
  3. #include <igl/cotmatrix_entries.h>
  4. #include <igl/edge_lengths.h>
  5. #include <igl/EPS.h>
  6. TEST_CASE("cotmatrix_entries: simple", "[igl]")
  7. {
  8. //The allowed error for this test
  9. const double epsilon = 1e-15;
  10. Eigen::MatrixXd V;
  11. Eigen::MatrixXi F;
  12. //This is a cube of dimensions 1.0x1.0x1.0
  13. test_common::load_mesh("cube.obj", V, F);
  14. //Prepare another mesh with triangles along side diagonals of the cube
  15. //These triangles are form a regular tetrahedron of side sqrt(2)
  16. Eigen::MatrixXi F_tet(4,3);
  17. F_tet << 4,6,1,
  18. 6,4,3,
  19. 4,1,3,
  20. 1,6,3;
  21. //1. Check cotmatrix_entries
  22. Eigen::MatrixXd C1;
  23. igl::cotmatrix_entries(V,F,C1);
  24. REQUIRE (C1.rows() == F.rows());
  25. REQUIRE (C1.cols() == 3);
  26. //All angles in unit cube measure 45 or 90 degrees
  27. //Their (half)cotangent must value 0.5 or 0.0
  28. for(int f = 0;f<C1.rows();f++)
  29. {
  30. #ifdef IGL_EDGE_LENGTHS_SQUARED_H
  31. //Hard assert if we have edge_length_squared
  32. for(int v = 0;v<3;v++)
  33. if (C1(f,v) > 0.1)
  34. REQUIRE (C1(f,v) == 0.5);
  35. else
  36. REQUIRE (C1(f,v) == 0.0);
  37. //All cotangents sum 1.0 for those triangles
  38. REQUIRE (C1.row(f).sum() == 1.0);
  39. #else
  40. //Soft assert if we have not edge_length_squared
  41. for(int v = 0;v<3;v++)
  42. if (C1(f,v) > 0.1)
  43. REQUIRE (C1(f,v) == Approx (0.5).margin( epsilon));
  44. else
  45. REQUIRE (C1(f,v) == Approx (0.0).margin( epsilon));
  46. //All cotangents sum 1.0 for those triangles
  47. REQUIRE (C1.row(f).sum() == Approx (1.0).margin( epsilon));
  48. #endif
  49. }
  50. //Check the regular tetrahedron
  51. Eigen::MatrixXd C2;
  52. igl::cotmatrix_entries(V,F_tet,C2);
  53. REQUIRE (C2.rows() == F_tet.rows());
  54. REQUIRE (C2.cols() == 3);
  55. for(int f = 0;f<C2.rows();f++)
  56. {
  57. //Their (half)cotangent must value 0.5 / tan(igl::PI / 3.0)
  58. for(int v = 0;v<3;v++)
  59. REQUIRE (C2(f,v) == Approx (0.5 / tan(igl::PI / 3.0)).margin( epsilon));
  60. }
  61. //Scale the cube to have huge sides
  62. Eigen::MatrixXd V_huge = V * 1.0e8;
  63. igl::cotmatrix_entries(V_huge,F,C1);
  64. REQUIRE (C1.rows() == F.rows());
  65. REQUIRE (C1.cols() == 3);
  66. //All angles still measure 45 or 90 degrees
  67. //Their (half)cotangent must value 0.5 or 0.0
  68. for(int f = 0;f<C1.rows();f++)
  69. {
  70. #ifdef IGL_EDGE_LENGTHS_SQUARED_H
  71. //Hard assert if we have edge_length_squared
  72. for(int v = 0;v<3;v++)
  73. if (C1(f,v) > 0.1)
  74. REQUIRE (C1(f,v) == 0.5);
  75. else
  76. REQUIRE (C1(f,v) == 0.0);
  77. //All cotangents sum 1.0 for those triangles
  78. REQUIRE (C1.row(f).sum() == 1.0);
  79. #else
  80. //Soft assert if we have not edge_length_squared
  81. for(int v = 0;v<3;v++)
  82. if (C1(f,v) > 0.1)
  83. REQUIRE (C1(f,v) == Approx (0.5).margin( epsilon));
  84. else
  85. REQUIRE (C1(f,v) == Approx (0.0).margin( epsilon));
  86. //All cotangents sum 1.0 for those triangles
  87. REQUIRE (C1.row(f).sum() == Approx (1.0).margin( epsilon));
  88. #endif
  89. }
  90. //Check the huge regular tetrahedron
  91. igl::cotmatrix_entries(V_huge,F_tet,C2);
  92. REQUIRE (C2.rows() == F_tet.rows());
  93. REQUIRE (C2.cols() == 3);
  94. for(int f = 0;f<C2.rows();f++)
  95. {
  96. //Their (half)cotangent must value 0.5 / tan(igl::PI / 3.0)
  97. for(int v = 0;v<3;v++)
  98. REQUIRE (C2(f,v) == Approx (0.5 / tan(igl::PI / 3.0)).margin( epsilon));
  99. }
  100. //Scale the cube to have tiny sides
  101. Eigen::MatrixXd V_tiny = V * 1.0e-8;
  102. igl::cotmatrix_entries(V_tiny,F,C1);
  103. REQUIRE (C1.rows() == F.rows());
  104. REQUIRE (C1.cols() == 3);
  105. //All angles still measure 45 or 90 degrees
  106. //Their (half)cotangent must value 0.5 or 0.0
  107. for(int f = 0;f<C1.rows();f++)
  108. {
  109. for(int v = 0;v<3;v++)
  110. if (C1(f,v) > 0.1)
  111. REQUIRE (C1(f,v) == Approx (0.5).margin( epsilon));
  112. else
  113. REQUIRE (C1(f,v) == Approx (0.0).margin( epsilon));
  114. //All cotangents sum 1.0 for those triangles
  115. REQUIRE (C1.row(f).sum() == Approx (1.0).margin( epsilon));
  116. }
  117. //Check the tiny regular tetrahedron
  118. igl::cotmatrix_entries(V_tiny,F_tet,C2);
  119. REQUIRE (C2.rows() == F_tet.rows());
  120. REQUIRE (C2.cols() == 3);
  121. for(int f = 0;f<C2.rows();f++)
  122. {
  123. //Their (half)cotangent must value 0.5 / tan(igl::PI / 3.0)
  124. for(int v = 0;v<3;v++)
  125. REQUIRE (C2(f,v) == Approx (0.5 / tan(igl::PI / 3.0)).margin( epsilon));
  126. }
  127. }// TEST_CASE("cotmatrix_entries: simple", "[igl]")
  128. TEST_CASE("cotmatrix_entries: intrinsic", "[igl]")
  129. {
  130. Eigen::MatrixXd V;
  131. Eigen::MatrixXi F;
  132. //This is a cube of dimensions 1.0x1.0x1.0
  133. test_common::load_mesh("cube.obj", V, F);
  134. Eigen::MatrixXd Cext,Cint;
  135. // compute C extrinsically
  136. igl::cotmatrix_entries(V,F,Cext);
  137. // compute C intrinsically
  138. Eigen::MatrixXd l;
  139. igl::edge_lengths(V,F,l);
  140. igl::cotmatrix_entries(l,Cint);
  141. test_common::assert_near(Cext,Cint,igl::EPS<double>());
  142. }