fix dihedral angle computation

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

@@ -930,11 +930,42 @@ namespace CommonKernelFunctors {
      const Vector_3 ad = vector(a,d);
 
      const Vector_3 abad = cross_product(ab,ad);
-     const double x = CGAL::to_double(scalar_product(cross_product(ab,ac), abad));
+     const Vector_3 abac = cross_product(ab,ac);
+
+     // 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)
+     //
+     // 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 thus abad is in
+     // the plane (yz) of the new basis.
+     //
+     // In that basis, 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
+     // 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(ac,abad));
+     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));
 
-     return FT(std::atan2(y, x) * 180 / CGAL_PI );
+     return FT(std::atan2(z, y) * 180 / CGAL_PI );
    }
  };
 

Kernel_23/test/Kernel_23/test_approximate_dihedral_angle_3.cpp

@@ -10,10 +10,51 @@ struct query {
   double expected_angle;
 };
 
+void sign_test()
+{
+  K::Point_3 a(0,0,0), b(1,0,0), c(0,1, 0), d(0,0,1);
+
+  assert( CGAL::approximate_dihedral_angle(a, b, c, d) > 0);
+  assert( CGAL::approximate_dihedral_angle(c, a, b, d) > 0);
+  assert( CGAL::approximate_dihedral_angle(a, d, b, c) > 0);
+  assert( CGAL::approximate_dihedral_angle(c, b, d, a) > 0);
+  assert( CGAL::approximate_dihedral_angle(d, b, a, c) > 0);
+  assert( CGAL::approximate_dihedral_angle(d, c, b, a) > 0);
+
+  assert( CGAL::approximate_dihedral_angle(a, b, d, c) < 0);
+  assert( CGAL::approximate_dihedral_angle(c, a, d, b) < 0);
+  assert( CGAL::approximate_dihedral_angle(a, d, c, b) < 0);
+  assert( CGAL::approximate_dihedral_angle(c, b, a, d) < 0);
+  assert( CGAL::approximate_dihedral_angle(d, b, c, a) < 0);
+  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.},
@@ -26,11 +67,27 @@ int main() {
     { {  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);
 }

Surface_mesh_segmentation/include/CGAL/Surface_mesh_segmentation/internal/Surface_mesh_segmentation.h

@@ -324,8 +324,8 @@ class Surface_mesh_segmentation
 
     CGAL_precondition( (! (face(edge,mesh)==boost::graph_traits<Polyhedron>::null_face()))
                        && (! (face(opposite(edge,mesh),mesh)==boost::graph_traits<Polyhedron>::null_face())) );
-    const Point a = get(vertex_point_pmap,target(edge,mesh));
-    const Point b = get(vertex_point_pmap,target(prev(edge,mesh),mesh));
+    const Point a = get(vertex_point_pmap,source(edge,mesh));
+    const Point b = get(vertex_point_pmap,target(edge,mesh));
     const Point c = get(vertex_point_pmap,target(next(edge,mesh),mesh));
     const Point d = get(vertex_point_pmap,target(next(opposite(edge,mesh),mesh),mesh));
     // As far as I check: if, say, dihedral angle is 5, this returns 175,