add tests for the dihedral angle

Data/data/meshes/regular_tetrahedron.off

@@ -0,0 +1,11 @@
+OFF
+4 4 0
+-1 0 -0.707107
+1 0 -0.707107
+0 -1 0.707107
+0 1 0.707107
+3 0 1 2
+3 0 3 1
+3 0 2 3
+3 1 3 2
+

Kernel_23/include/CGAL/Kernel/function_objects.h

@@ -932,22 +932,35 @@ namespace CommonKernelFunctors {
      const Vector_3 abad = cross_product(ab,ad);
      const Vector_3 abac = cross_product(ab,ac);
 
-     // The dihedral angle we are interested in is the angle around the edge ab which is
-     // the same as the angle between the vectors abac and abad
+     // The dihedral angle we are interested in is the angle around the oriented
+     // edge ab which is the same (in absolute value) as the angle between the
+     // vectors ab^ac and ab^ad (cross-products).
      // (abac points inside the tetra abcd if its orientation is positive and outside otherwise)
-     // In order to increase the numerical precision of the computation, we consider the
-     // vector abad in the basis defined by (ab, abac, ab^abac).
-     // In this basis, adab=(ab * abad, abac * abad, [ab^abac] * abad), as all basis vectors are orthogonal.
-     // We have ab * abad = 0, so in the plane (yz) of the new basis
-     // the dihedral angle is the angle between the y axis and abad
-     // which is the arctan of y/z (up to normalization)
-     // (Note that ab^abac is in the plane abc, pointing outside the tetra if its orientation is positive and inside otherwise).
+     //
+     // We consider the vector abad in the basis defined by the three vectors
+     //    (<ab>, <abac>, <ab^abac>)
+     // where <u> denote the normalized vector u/|u|.
+     //
+     // In this orthonormal basis, the vector adab has the coordinates
+     //    x = <ab>      * abad
+     //    y = <abac>    * abad
+     //    z = <ab^abac> * abad
+     // We have x == 0, because abad and ab are orthogonal, and then abad is in
+     // the plane (yz) of the new basis.
+     //
+     // In that basic, the dihedral angle is the angle between the y axis and abad
+     // which is the arctan of y/z, or atan2(z, y).
+     //
+     // (Note that ab^abac is in the plane abc, pointing outside the tetra if
+     //  its orientation is positive and inside otherwise).
+     //
      // For the normalization, abad appears in both scalar products
      // in the quotient so we can ignore its norm. For the second
      // terms of the scalar products, we are left with ab^abac and abac.
-     // Since ab and abac are orthogonal the sinus of the angle between the vector is 1
-     // so the norms are ||ab||.||abac|| vs ||abac||, which is why we have a multiplication by ||ab||
-     // in y below.
+     // Since ab and abac are orthogonal, the sinus of the angle between the
+     // two vectors is 1.
+     // So the norms are |ab|.|abac| vs |abac|, which is why we have a
+     // multiplication by |ab| in y below.
      const double l_ab = CGAL::sqrt(CGAL::to_double(sq_distance(a,b)));
      const double y = l_ab * CGAL::to_double(scalar_product(abac, abad));
      const double z = CGAL::to_double(scalar_product(cross_product(ab,abac),abad));

Kernel_23/test/Kernel_23/test_approximate_dihedral_angle_3.cpp

@@ -29,29 +29,65 @@ void sign_test()
   assert( CGAL::approximate_dihedral_angle(d, c, a, b) < 0);
 }
 
+auto almost_equal_angle(double a, double b) {
+  return std::min(std::abs(a - b), std::abs(a + 360 - b)) < 0.1;
+}
+
+void test_regular_tetrahedron()
+{
+  auto half_root_of_2 = std::sqrt(2) / 2;
+
+  // Regular tetrahedron
+  Point_3 a{ -1,  0, -half_root_of_2};
+  Point_3 b{  1,  0, -half_root_of_2};
+  Point_3 c{  0,  1,  half_root_of_2};
+  Point_3 d{  0, -1,  half_root_of_2};
+  assert(orientation(a, b, c, d) == CGAL::POSITIVE);
+  assert(almost_equal_angle(CGAL::approximate_dihedral_angle(a, b, c, d), 70.5288));
+}
+
 int main() {
+  std::cout.precision(17);
   sign_test();
+  test_regular_tetrahedron();
 
   Point_3 a = {0, 0, 0};
-  Point_3 b = {0, 1, 0};
-  Point_3 c = {1, 0, 0};
+  Point_3 b = {0, -1, 0}; // ab is oriented so that it sees the plan xz positively.
+  [[maybe_unused]] Point_3 c = {1, 0, 0};
+  // c can be any point in the half-plane xy, with x>0
 
   const query queries[] = {
     { {  1, 0,  0},    0.},
-    { {  1, 0,  1},   -45.},
-    { {  0, 0,  1},   -90.},
-    { { -1, 0,  1},  -135.},
-    //{ { -1, 0,  0},  -180.},
-    { { -1, 0, -1}, 135.},
-    { {  0, 0, -1},  90.},
-    { {  1, 0, -1},  45.},
+    { {  1, 0,  1},   45.},
+    { {  0, 0,  1},   90.},
+    { { -1, 0,  1},  135.},
+    { { -1, 0,  0},  180.},
+    { { -1, 0, -1}, -135.},
+    { {  0, 0, -1},  -90.},
+    { {  1, 0, -1},  -45.},
   };
 
-  for(auto query: queries) {
-    const auto& expected = query.expected_angle;
-    const auto& p = query.p;
-    auto approx = CGAL::approximate_dihedral_angle(a, b, c, p);
-    std::cout << approx << "  -- " << expected << '\n';
-    assert( std::abs(approx - expected) < 0.1 );
+  auto cnt = 0u;
+  for(double yc = -10; yc < 10; yc += 0.1) {
+    Point_3 c{1, yc, 0};
+    // std::cout << "c = " << c << '\n';
+    for(const auto& query : queries) {
+      for(double yp = -10; yp < 10; yp += 0.3) {
+        const auto& expected = query.expected_angle;
+        const Point_3 p{query.p.x(), yp, query.p.z()};
+        // std::cout << "p = " << p << '\n';
+        auto approx = CGAL::approximate_dihedral_angle(a, b, c, p);
+        // std::cout << approx << "  -- " << expected << '\n';
+        if(!almost_equal_angle(approx, expected)) {
+          std::cout << "ERROR:\n";
+          std::cout << "CGAL::approximate_dihedral_angle(" << a << ", " << b << ", " << c << ", " << p << ") = " << approx << '\n';
+          std::cout << "expected: " << expected << '\n';
+          return 1;
+        }
+        ++cnt;
+      }
+    }
   }
+  std::cout << "OK (" << cnt << " tests)\n";
+  assert(cnt > 10000);
 }