From 5247f0fbdbdb9b754050bd03e19e646267f23740 Mon Sep 17 00:00:00 2001 From: Andreas Fabri Date: Tue, 19 Sep 2023 15:43:02 +0100 Subject: [PATCH] Bounding_volumes: Use Eigen library in demo --- .../demo/Bounding_volumes/CMakeLists.txt | 10 +- GraphicsView/demo/Bounding_volumes/Ellipse.h | 92 ++++++++++--------- .../demo/Bounding_volumes/data.pts.cgal | 3 + .../demo/Bounding_volumes/rect.pts.cgal | 4 + .../demo/Bounding_volumes/slanted.pts.cgal | 4 + 5 files changed, 68 insertions(+), 45 deletions(-) create mode 100644 GraphicsView/demo/Bounding_volumes/data.pts.cgal create mode 100644 GraphicsView/demo/Bounding_volumes/rect.pts.cgal create mode 100644 GraphicsView/demo/Bounding_volumes/slanted.pts.cgal diff --git a/GraphicsView/demo/Bounding_volumes/CMakeLists.txt b/GraphicsView/demo/Bounding_volumes/CMakeLists.txt index 78f97078c03e..ac02fc6554c9 100644 --- a/GraphicsView/demo/Bounding_volumes/CMakeLists.txt +++ b/GraphicsView/demo/Bounding_volumes/CMakeLists.txt @@ -18,6 +18,10 @@ find_package(CGAL REQUIRED OPTIONAL_COMPONENTS Qt5) find_package(Qt5 QUIET COMPONENTS Widgets) include_directories(BEFORE ./include) +find_package(Eigen3 3.1.91 QUIET) #(requires 3.1.91 or greater) +include(CGAL_Eigen3_support) + +if(TARGET CGAL::Eigen3_support) if(CGAL_Qt5_FOUND AND Qt5_FOUND) set(CMAKE_INCLUDE_CURRENT_DIR ON) @@ -40,7 +44,8 @@ if(CGAL_Qt5_FOUND AND Qt5_FOUND) add_to_cached_list(CGAL_EXECUTABLE_TARGETS Bounding_volumes) target_link_libraries(Bounding_volumes PRIVATE CGAL::CGAL CGAL::CGAL_Qt5 - Qt5::Widgets) + Qt5::Widgets + PUBLIC CGAL::Eigen3_support) include(${CGAL_MODULES_DIR}/CGAL_add_test.cmake) cgal_add_compilation_test(Bounding_volumes) @@ -53,3 +58,6 @@ else() message("NOTICE: This demo requires CGAL and Qt5, and will not be compiled.") endif() +else() + message("NOTICE: This demo requires CGAL and Eigen, and will not be compiled.") +endif() diff --git a/GraphicsView/demo/Bounding_volumes/Ellipse.h b/GraphicsView/demo/Bounding_volumes/Ellipse.h index ce5e09a556c3..d950b8fd8061 100644 --- a/GraphicsView/demo/Bounding_volumes/Ellipse.h +++ b/GraphicsView/demo/Bounding_volumes/Ellipse.h @@ -1,6 +1,9 @@ #ifndef ELLIPSE_H #define ELLIPSE_H +#include +#include + template class Ellipse_2 { @@ -20,50 +23,51 @@ template CGAL::Simple_cartesian::Conic_2 co; me.ellipse().double_conic(co); - double a11 = co.r(); - double a12 = co.t()/2.0; - double a22 = co.s(); - double b1 = co.u(); - double b2 = co.v(); - double c = co.w(); - - assert(a11*a22-a12*a12 >= 0); - - double delta = (a11>a22)? a11-a22 : a22-a11; - if( (delta < 0.00000000001) && (a12 == 0)){ - ce = Point_2(-b1/(2*a11), -b2/(2*a11)); - a_ = b_ = std::sqrt(b1*b1+b2*b2-4*a11*c)/(2*a11); - va_ = Vector_2(a_, 0); - vb_ = Vector_2(0, b_); - return; - } - - double kden = 2.0*(a12*a12 - a11*a22); - double k1 = (a22*b1 - a12*b2)/kden; - double k2 = (a11*b2 - a12*b1)/kden; - ce = Point_2(k1,k2); - double mu = 1/(a11*k1*k1 + 2.0*a12*k1*k2 + a22*k2*k2 - c); - double m11 = mu*a11; - double m12 = mu*a12; - double m22 = mu*a22; - - double r = std::sqrt((m11-m22)*(m11-m22) + 4.0*m12*m12); - double lambda1 = ((m11+m22) + r)/2.0; - double lambda2 = ((m11+m22) - r)/2.0; - - b_ = 1.0/sqrt(lambda1); - a_ = 1.0/sqrt(lambda2); - double omega = 1.0/ - std::sqrt((lambda1-m22)*(lambda1-m22) + m12*m12); - - double u1x = - a_ * omega*m12; - double u1y = a_ * omega * (lambda1-m22); - - double u2x = b_* omega * (lambda1-m22); - double u2y = b_* omega*m12; - - va_ = Vector_2(u1x,u1y); - vb_ = Vector_2(u2x,u2y); + // https://en.wikipedia.org/wiki/Matrix_representation_of_conic_sections + double A = co.r(); + double B = co.t(); + double C = co.s(); + double D = co.u(); + double E = co.v(); + double F = co.w(); + + assert(A * C - B*B/4 >= 0); + + Eigen::Matrix2d A33, A33i; + A33 << A, B/2, B/2, C; + A33i = A33.inverse(); + + Eigen::Vector2d v; + v << -D/2 , -E/2; + + v = A33i * v; + + ce = Point_2(v[0],v[1]); + + Eigen::EigenSolver es(A33); + + Eigen::Matrix2cd ev; + ev = es.eigenvectors(); + double x0 = ev.col(0)[0].real(); + double y0 = ev.col(0)[1].real(); + double x1 = ev.col(1)[0].real(); + double y1 = ev.col(1)[1].real(); + + double lambda1 = es.eigenvalues()[0].real(); + double lambda2 = es.eigenvalues()[1].real(); + + Eigen::Matrix3d AQ; + AQ << A, B/2, D/2, + B/2, C, E/2, + D/2, E/2, F; + + double K = AQ.determinant() / A33.determinant(); + + a_ = sqrt(-K/lambda1); + b_ = sqrt(-K/lambda2); + + va_ = a_ * Vector_2(x0,y0); + vb_ = b_ * Vector_2(x1,y1); } diff --git a/GraphicsView/demo/Bounding_volumes/data.pts.cgal b/GraphicsView/demo/Bounding_volumes/data.pts.cgal new file mode 100644 index 000000000000..49b66d2f77a8 --- /dev/null +++ b/GraphicsView/demo/Bounding_volumes/data.pts.cgal @@ -0,0 +1,3 @@ +0 2 +2 1 +0 0 diff --git a/GraphicsView/demo/Bounding_volumes/rect.pts.cgal b/GraphicsView/demo/Bounding_volumes/rect.pts.cgal new file mode 100644 index 000000000000..ab69a7198547 --- /dev/null +++ b/GraphicsView/demo/Bounding_volumes/rect.pts.cgal @@ -0,0 +1,4 @@ +-2 1 +-2 -1 +2 1 +2 -1 diff --git a/GraphicsView/demo/Bounding_volumes/slanted.pts.cgal b/GraphicsView/demo/Bounding_volumes/slanted.pts.cgal new file mode 100644 index 000000000000..0af8ad54f6b2 --- /dev/null +++ b/GraphicsView/demo/Bounding_volumes/slanted.pts.cgal @@ -0,0 +1,4 @@ +-2 -1 +-1 -2 +1 2 +2 1