Merge pull request #7725 from afabri/Bounding_volumes-fix_Ellipse-GF

Bounding_volumes: Use Eigen library in demo

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()

GraphicsView/demo/Bounding_volumes/Ellipse.h

@@ -1,6 +1,9 @@
 #ifndef ELLIPSE_H
 #define ELLIPSE_H
 
+#include <Eigen/Dense>
+#include <Eigen/Eigenvalues>
+
 template <class R>
  class Ellipse_2 {
 
@@ -20,50 +23,51 @@ template <class R>
       CGAL::Simple_cartesian<double>::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<Eigen::Matrix2d> 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);
     }
 
 

GraphicsView/demo/Bounding_volumes/data.pts.cgal

@@ -0,0 +1,3 @@
+0 2
+2 1
+0 0

GraphicsView/demo/Bounding_volumes/rect.pts.cgal

@@ -0,0 +1,4 @@
+-2 1
+-2 -1
+2 1
+2 -1

GraphicsView/demo/Bounding_volumes/slanted.pts.cgal

@@ -0,0 +1,4 @@
+-2 -1
+-1 -2
+1 2
+2 1