From 6b570f767e9631bc02be240fd61deb6b15e2d693 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?S=C3=A9bastien=20Loriot?= Date: Tue, 28 Nov 2023 17:45:54 +0100 Subject: [PATCH 1/9] fix sign --- .../include/CGAL/Kernel/function_objects.h | 2 +- .../test_approximate_dihedral_angle_3.cpp | 35 +++++++++++++++---- 2 files changed, 29 insertions(+), 8 deletions(-) diff --git a/Kernel_23/include/CGAL/Kernel/function_objects.h b/Kernel_23/include/CGAL/Kernel/function_objects.h index 9b6af428e25..bb476ded1f2 100644 --- a/Kernel_23/include/CGAL/Kernel/function_objects.h +++ b/Kernel_23/include/CGAL/Kernel/function_objects.h @@ -934,7 +934,7 @@ namespace CommonKernelFunctors { 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)); - return FT(std::atan2(y, x) * 180 / CGAL_PI ); + return -FT(std::atan2(y, x) * 180 / CGAL_PI ); } }; diff --git a/Kernel_23/test/Kernel_23/test_approximate_dihedral_angle_3.cpp b/Kernel_23/test/Kernel_23/test_approximate_dihedral_angle_3.cpp index b39ce5c01e7..dfa703b86f7 100644 --- a/Kernel_23/test/Kernel_23/test_approximate_dihedral_angle_3.cpp +++ b/Kernel_23/test/Kernel_23/test_approximate_dihedral_angle_3.cpp @@ -10,20 +10,41 @@ 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); +} + int main() { + sign_test(); + Point_3 a = {0, 0, 0}; Point_3 b = {0, 1, 0}; Point_3 c = {1, 0, 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) { From 0145bafbc5079ea98e91f82d528220a7ce91123a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?S=C3=A9bastien=20Loriot?= Date: Wed, 29 Nov 2023 12:22:46 +0100 Subject: [PATCH 2/9] update formula will add comments in a upcoming commit --- Kernel_23/include/CGAL/Kernel/function_objects.h | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/Kernel_23/include/CGAL/Kernel/function_objects.h b/Kernel_23/include/CGAL/Kernel/function_objects.h index bb476ded1f2..ab23552b682 100644 --- a/Kernel_23/include/CGAL/Kernel/function_objects.h +++ b/Kernel_23/include/CGAL/Kernel/function_objects.h @@ -930,11 +930,12 @@ 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); 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 ); } }; From 0e3f3a33d1165fa1462f4cabd6973c9a87cb6e4c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?S=C3=A9bastien=20Loriot?= Date: Wed, 29 Nov 2023 12:27:34 +0100 Subject: [PATCH 3/9] test to be fixed --> tetra orientation is 0 --- Kernel_23/test/Kernel_23/test_approximate_dihedral_angle_3.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Kernel_23/test/Kernel_23/test_approximate_dihedral_angle_3.cpp b/Kernel_23/test/Kernel_23/test_approximate_dihedral_angle_3.cpp index dfa703b86f7..ec352fb21dd 100644 --- a/Kernel_23/test/Kernel_23/test_approximate_dihedral_angle_3.cpp +++ b/Kernel_23/test/Kernel_23/test_approximate_dihedral_angle_3.cpp @@ -41,7 +41,7 @@ int main() { { { 1, 0, 1}, -45.}, { { 0, 0, 1}, -90.}, { { -1, 0, 1}, -135.}, - { { -1, 0, 0}, -180.}, + //{ { -1, 0, 0}, -180.}, { { -1, 0, -1}, 135.}, { { 0, 0, -1}, 90.}, { { 1, 0, -1}, 45.}, From 5879bb72c63a5a8f0a15a75da806778c6ad20dbb Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?S=C3=A9bastien=20Loriot?= Date: Wed, 29 Nov 2023 13:11:21 +0100 Subject: [PATCH 4/9] add comments about the formula --- .../include/CGAL/Kernel/function_objects.h | 17 +++++++++++++++++ 1 file changed, 17 insertions(+) diff --git a/Kernel_23/include/CGAL/Kernel/function_objects.h b/Kernel_23/include/CGAL/Kernel/function_objects.h index ab23552b682..73f6c672d15 100644 --- a/Kernel_23/include/CGAL/Kernel/function_objects.h +++ b/Kernel_23/include/CGAL/Kernel/function_objects.h @@ -931,6 +931,23 @@ 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 + // (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 vector are orthogonal. + // We have ab * abad = 0, so in the plane (zy) of the new basis + // the dihedral angle is the angle between the z axis and abad + // which is the arctan of y/z (up to normalization) + // (Note that ab^abac is in the plane abc, pointing inside the tetra if its orientation is positive and outside 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. 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)); From 4bed66e82f142e02066ade584e27e93f5f136559 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?S=C3=A9bastien=20Loriot?= Date: Wed, 29 Nov 2023 13:17:10 +0100 Subject: [PATCH 5/9] fix description --- Kernel_23/include/CGAL/Kernel/function_objects.h | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/Kernel_23/include/CGAL/Kernel/function_objects.h b/Kernel_23/include/CGAL/Kernel/function_objects.h index 73f6c672d15..55620fd2ed3 100644 --- a/Kernel_23/include/CGAL/Kernel/function_objects.h +++ b/Kernel_23/include/CGAL/Kernel/function_objects.h @@ -937,11 +937,11 @@ namespace CommonKernelFunctors { // (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 vector are orthogonal. - // We have ab * abad = 0, so in the plane (zy) of the new basis - // the dihedral angle is the angle between the z axis and abad + // 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 inside the tetra if its orientation is positive and outside otherwise). + // (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. From e3e5bcd34448a62c7e8d5e129ea8041dff359560 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?S=C3=A9bastien=20Loriot?= Date: Wed, 29 Nov 2023 15:25:44 +0100 Subject: [PATCH 6/9] accomodate handle sign face orientation was inversed and somehow matched the bug in the dihedral angle computation function --- .../internal/Surface_mesh_segmentation.h | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Surface_mesh_segmentation/include/CGAL/Surface_mesh_segmentation/internal/Surface_mesh_segmentation.h b/Surface_mesh_segmentation/include/CGAL/Surface_mesh_segmentation/internal/Surface_mesh_segmentation.h index d0a77fa1a78..a2a2228f6f5 100644 --- a/Surface_mesh_segmentation/include/CGAL/Surface_mesh_segmentation/internal/Surface_mesh_segmentation.h +++ b/Surface_mesh_segmentation/include/CGAL/Surface_mesh_segmentation/internal/Surface_mesh_segmentation.h @@ -324,8 +324,8 @@ private: CGAL_precondition( (! (face(edge,mesh)==boost::graph_traits::null_face())) && (! (face(opposite(edge,mesh),mesh)==boost::graph_traits::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, From df4eed9302f873bf92f30941c9d7de374b69be14 Mon Sep 17 00:00:00 2001 From: Laurent Rineau Date: Wed, 29 Nov 2023 15:27:57 +0100 Subject: [PATCH 7/9] add tests for the dihedral angle --- Data/data/meshes/regular_tetrahedron.off | 11 ++++ .../include/CGAL/Kernel/function_objects.h | 37 +++++++---- .../test_approximate_dihedral_angle_3.cpp | 66 ++++++++++++++----- 3 files changed, 87 insertions(+), 27 deletions(-) create mode 100644 Data/data/meshes/regular_tetrahedron.off diff --git a/Data/data/meshes/regular_tetrahedron.off b/Data/data/meshes/regular_tetrahedron.off new file mode 100644 index 00000000000..08a25ad43b3 --- /dev/null +++ b/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 + diff --git a/Kernel_23/include/CGAL/Kernel/function_objects.h b/Kernel_23/include/CGAL/Kernel/function_objects.h index 55620fd2ed3..6bc98833dab 100644 --- a/Kernel_23/include/CGAL/Kernel/function_objects.h +++ b/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 + // (, , ) + // where denote the normalized vector u/|u|. + // + // In this orthonormal basis, the vector adab has the coordinates + // x = * abad + // y = * abad + // z = * 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)); diff --git a/Kernel_23/test/Kernel_23/test_approximate_dihedral_angle_3.cpp b/Kernel_23/test/Kernel_23/test_approximate_dihedral_angle_3.cpp index ec352fb21dd..9ad255befac 100644 --- a/Kernel_23/test/Kernel_23/test_approximate_dihedral_angle_3.cpp +++ b/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); } From 46877d41341d7426ec67aeded87530988a882299 Mon Sep 17 00:00:00 2001 From: Laurent Rineau Date: Wed, 29 Nov 2023 15:34:36 +0100 Subject: [PATCH 8/9] fix two typos --- Kernel_23/include/CGAL/Kernel/function_objects.h | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Kernel_23/include/CGAL/Kernel/function_objects.h b/Kernel_23/include/CGAL/Kernel/function_objects.h index 6bc98833dab..7848f55e94c 100644 --- a/Kernel_23/include/CGAL/Kernel/function_objects.h +++ b/Kernel_23/include/CGAL/Kernel/function_objects.h @@ -945,10 +945,10 @@ namespace CommonKernelFunctors { // x = * abad // y = * abad // z = * abad - // We have x == 0, because abad and ab are orthogonal, and then abad is in + // We have x == 0, because abad and ab are orthogonal, and thus 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 + // 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 From 9822b6b0855d0a577af2b9b225c8925680c03628 Mon Sep 17 00:00:00 2001 From: Sebastien Loriot Date: Thu, 30 Nov 2023 09:20:53 +0100 Subject: [PATCH 9/9] protect from macro substitution Co-authored-by: Andreas Fabri --- Kernel_23/test/Kernel_23/test_approximate_dihedral_angle_3.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Kernel_23/test/Kernel_23/test_approximate_dihedral_angle_3.cpp b/Kernel_23/test/Kernel_23/test_approximate_dihedral_angle_3.cpp index 9ad255befac..dc3ccdaf989 100644 --- a/Kernel_23/test/Kernel_23/test_approximate_dihedral_angle_3.cpp +++ b/Kernel_23/test/Kernel_23/test_approximate_dihedral_angle_3.cpp @@ -30,7 +30,7 @@ void sign_test() } auto almost_equal_angle(double a, double b) { - return std::min(std::abs(a - b), std::abs(a + 360 - b)) < 0.1; + return (std::min)(std::abs(a - b), std::abs(a + 360 - b)) < 0.1; } void test_regular_tetrahedron()