The Higher Education and Research forge

Home My Page Projects Code Snippets Project Openings Garamon
Summary Activity SCM

SCM Repository

1 // Copyright (c) 2018 by University Paris-Est Marne-la-Vallee
2 // qc3ga2Tools.hpp
3 // Authors: Vincent Nozick and Stephane Breuils 
4 // Contact: vincent.nozick@u-pem.fr
7 /// \file qc3ga2Tools.hpp
8 /// \author Vincent Nozick, and Stephane Breuils
9 /// \brief some useful functions when using Quadric Conformal Geometric Algebra of R^3. Use this file if you generated the lib using the "qc3ga.conf" configuration file. Quadric conformal geometric algebra of R^3 is described in the following paper: Stéphane Breuils, Vincent Nozick, Akihiro Sugimoto, Eckhard Hitzer, Quadric Conformal Geometric Algebra of R9,6, Advances in Applied Clifford Algebras, Springer Verlag, 2018, 28 (2), pp.35 
12 // Anti-doublon
13 #ifndef QC3GA_TOOLS_HPP__
14 #define QC3GA_TOOLS_HPP__
15 #pragma once
17 // External includes
18 #include <string>
19 #include <vector>
21 // Internal Includes
22 #include <qc3ga/Mvec.hpp>
23 #include <Eigen/Dense>        // rank in display quadric
24 #include <Eigen/Eigenvalues>  // eigen solver
26 /// \namespace grouping the multivectors object
27 namespace qc3ga{
29     /// \brief build a point from a vector
30     /// \param x vector component related to e1
31     /// \param y vector component related to e2
32     /// \param z vector component related to e3
33     /// \return a multivector corresponding to a point of qc3ga
34     template<typename T>
35     qc3ga::Mvec<T> point(const T &x, const T &y, const T &z){
36         qc3ga::Mvec<T> mv;
37         mv[qc3ga::E1] = x;
38         mv[qc3ga::E2] = y;
39         mv[qc3ga::E3] = z;
40         mv[qc3ga::Ei1] = 0.5*(x*x);
41         mv[qc3ga::Ei2] = 0.5*(y*y);
42         mv[qc3ga::Ei3] = 0.5*(z*z);
43         mv[qc3ga::Ei4] = (x*y);
44         mv[qc3ga::Ei5] = (x*z);
45         mv[qc3ga::Ei6] = (y*z);
46                 mv[qc3ga::E01] = mv[qc3ga::E02] = mv[qc3ga::E03] = 1.0;
47         return mv;
48     }
50     template<typename T>
51     qc3ga::Mvec<T> point(const Eigen::Matrix<T,3,1> &pt){
52         qc3ga::Mvec<T> mv;
53         mv[qc3ga::E1] = pt(0);
54         mv[qc3ga::E2] = pt(1);
55         mv[qc3ga::E3] = pt(2);
56         mv[qc3ga::Ei1] = 0.5*(pt(0)*pt(0));
57         mv[qc3ga::Ei2] = 0.5*(pt(1)*pt(1));
58         mv[qc3ga::Ei3] = 0.5*(pt(2)*pt(2));
59         mv[qc3ga::Ei4] = pt(0)*pt(1);
60         mv[qc3ga::Ei5] = pt(0)*pt(2);
61         mv[qc3ga::Ei6] = pt(1)*pt(2);
62         mv[qc3ga::E01] = mv[qc3ga::E02] = mv[qc3ga::E03] = 1.0;
63         return mv;
64     }
66     template<typename T>
67     qc3ga::Mvec<T> sphere(const T &cx, const T &cy, const T &cz, const T &radius){
68         qc3ga::Mvec<T> mv = qc3ga::point<T>(cx, cy, cz) - radius * radius * (qc3ga::ei1<T>() + qc3ga::ei2<T>() + qc3ga::ei3<T>()) / 6.0;
69         return mv;
70     }
72     /// a quadric has the form (a b c d e f g h i j)
73     /// with ax^2 + by^2 + cz^2 + dxy + exz + fyz + gxw + hyw + izw + jw^2 = 0 
74     template<typename T>
75     std::vector<T> ga2DualQuadric(const qc3ga::Mvec<T> &mv){
77         std::vector<T> quadric(10);
78         quadric[0] = - mv[qc3ga::E01] / 2.0;
79         quadric[1] = - mv[qc3ga::E02] / 2.0;
80         quadric[2] = - mv[qc3ga::E03] / 2.0;
81         quadric[3] = - mv[qc3ga::E04] ;
82         quadric[4] = - mv[qc3ga::E05] ;
83         quadric[5] = - mv[qc3ga::E06] ;
85         quadric[6] =  mv[qc3ga::E1];
86         quadric[7] =  mv[qc3ga::E2];
87         quadric[8] =  mv[qc3ga::E3];
89         quadric[9] = - 3 * mv[qc3ga::Ei1];
91         return quadric;
92     }
95     /// a quadric has the form (a b c d e f g h i j)
96     /// with ax^2 + by^2 + cz^2 + dxy + exz + fyz + gx + hy + iz + j = 0 
97     template<typename T>
98     qc3ga::Mvec<T> dualQuadric2ga(const T a, const T b, const T c, const T d, const T e, const T f, const T g, const T h, const T i, const T j){
99         qc3ga::Mvec<T> mv;
100         mv[qc3ga::E01] = - 2.0 * a; // x^2
101         mv[qc3ga::E02] = - 2.0 * b; // y^2
102         mv[qc3ga::E03] = - 2.0 * c; // z^2
103         mv[qc3ga::E04] = - d;       // xy
104         mv[qc3ga::E05] = - e;       // xz
105         mv[qc3ga::E06] = - f;       // yz
106         mv[qc3ga::E1] = g;          // x
107         mv[qc3ga::E2] = h;          // y
108         mv[qc3ga::E3] = i;          // z
110         mv[qc3ga::Ei1] = mv[qc3ga::Ei2] = mv[qc3ga::Ei3] = - j / 3.0;
111         
112         return mv;
113     }
115     /// a quadric has the form (a b c d e f g h i j)
116     /// with ax^2 + by^2 + cz^2 + dxy + exz + fyz + gx + hy + iz + j = 0 
117     template<typename T>
118     qc3ga::Mvec<T> planeToDualGA(const T a, const T b, const T c, const T d){
119         qc3ga::Mvec<T> mv;
120         mv[qc3ga::E1] = a;          // x
121         mv[qc3ga::E2] = b;          // y
122         mv[qc3ga::E3] = c;          // z
124         mv[qc3ga::Ei1] = mv[qc3ga::Ei2] = mv[qc3ga::Ei3] = - d / 3.0; // 1
125         
126         return mv;
127     }
131     /// a quadric has the form (a b c d e f g h i j)
132     /// with ax^2 + by^2 + cz^2 + dxy + exz + fyz + gx + hy + iz + j = 0 
133     template<typename T>
134     qc3ga::Mvec<T> dualQuadric2ga(const std::vector<T> &quadric){
135         return dualQuadric2ga<T>(quadric[0], quadric[1], quadric[2], quadric[3], quadric[4], quadric[5], quadric[6], quadric[7], quadric[8], quadric[9]);
136     }
141         
142         /// a quadric has the form (a b c d e f g h i j)
143     /// with ax^2 + by^2 + cz^2 + dxy + exz + fyz + gxw + hyw + izw + jw^2 = 0 
144     template<typename T>
145     std::vector<T> ga2Quadric(const qc3ga::Mvec<T> &mv){
147         std::vector<T> quadric(10);
148         quadric[0] = - mv[qc3ga::E1230102i203i304i405i506i6] / 2.0; // dual(e01)
149         quadric[1] = - mv[qc3ga::E12301i10203i304i405i506i6] / 2.0; // dual(e02)
150         quadric[2] = - mv[qc3ga::E12301i102i20304i405i506i6] / 2.0; // dual(e03)
151         quadric[3] = - mv[qc3ga::E12301i102i203i30405i506i6]; // dual(e04)
152         quadric[4] = - mv[qc3ga::E12301i102i203i304i40506i6]; // dual(e05)
153         quadric[5] = - mv[qc3ga::E12301i102i203i304i405i506]; // dual(e06)
154                 
155         quadric[6] = -mv[qc3ga::E2301i102i203i304i405i506i6]; // dual(e1)
156         quadric[7] =  mv[qc3ga::E1301i102i203i304i405i506i6]; // dual(e2)
157         quadric[8] = -mv[qc3ga::E1201i102i203i304i405i506i6]; // dual(e3)
158                 
159         quadric[9] = 3 * mv[qc3ga::E123i102i203i304i405i506i6]; // dual(ei1)
161         return quadric;
162     }
164     /// a quadric has the form (a b c d e f g h i j)
165     /// with ax^2 + by^2 + cz^2 + dxy + exz + fyz + gxw + hyw + izw + jw^2 = 0
166     /// defines the primal form of a conic
167     template<typename T>
168     qc3ga::Mvec<T> quadric2ga(const T a, const T b, const T c, const T d, const T e, const T f, const T g, const T h, const T i, const T j){
169         qc3ga::Mvec<T> mv;
171         mv[qc3ga::E1230102i203i304i405i506i6] = - 2.0 * a; // dual(e01) component
172         mv[qc3ga::E12301i10203i304i405i506i6] = - 2.0 * b; // dual(e02) component
173         mv[qc3ga::E12301i102i20304i405i506i6] = - 2.0 * c; // dual(e03) component
174         mv[qc3ga::E12301i102i203i30405i506i6] = - d; // dual(e04) component
175         mv[qc3ga::E12301i102i203i304i40506i6] = - e; // dual(e05) component
176         mv[qc3ga::E12301i102i203i304i405i506] = - f; // dual(e06) component
179         mv[qc3ga::E2301i102i203i304i405i506i6] = g; // dual(e1) component 
180         mv[qc3ga::E1301i102i203i304i405i506i6] = -h; // dual(e2) component
181         mv[qc3ga::E1201i102i203i304i405i506i6] = i; // dual(e3) component
183         mv[qc3ga::E123i102i203i304i405i506i6] = mv[qc3ga::E12301i1i203i304i405i506i6] = mv[qc3ga::E12301i102i2i304i405i506i6] = j / 3.0; // dual(eij) component 
184         
185         return mv;
186     }
188     /// a quadric has the form (a b c d e f g h i j)
189     /// with ax^2 + by^2 + cz^2 + dxy + exz + fyz + gxw + hyw + izw + jw^2 = 0
190         ///     defines the primal form of a conic
191     template<typename T>
192     qc3ga::Mvec<T> quadric2ga(const std::vector<T> &quadric){
193         return qc3ga::quadric2ga<T>(quadric[0], quadric[1], quadric[2], quadric[3], quadric[4], quadric[5], quadric[6], quadric[7], quadric[8], quadric[9]);
194     }
197     /// \brief compute the dual of a 14-vector into a 1-vector without computing the inner product
198     template<typename T>
199     qc3ga::Mvec<T> primalSurfaceDualization( qc3ga::Mvec<T> &primalSurface){ // const ???
200         qc3ga::Mvec<T> mv;
202         mv[qc3ga::E01] = primalSurface[qc3ga::E1230102i203i304i405i506i6]; // dual(e01) component
203         mv[qc3ga::E02] = primalSurface[qc3ga::E12301i10203i304i405i506i6]; // dual(e02) component
204         mv[qc3ga::E03] = primalSurface[qc3ga::E12301i102i20304i405i506i6]; // dual(e03) component
205         mv[qc3ga::E04] = primalSurface[qc3ga::E12301i102i203i30405i506i6]; // dual(e04) component
206         mv[qc3ga::E05] = primalSurface[qc3ga::E12301i102i203i304i40506i6]; // dual(e05) component
207         mv[qc3ga::E06] = primalSurface[qc3ga::E12301i102i203i304i405i506]; // dual(e06) component
210         mv[qc3ga::E1]  = primalSurface[qc3ga::E2301i102i203i304i405i506i6];  // dual(e1) component 
211         mv[qc3ga::E2]  = -primalSurface[qc3ga::E1301i102i203i304i405i506i6]; // dual(e2) component
212         mv[qc3ga::E3]  = primalSurface[qc3ga::E1201i102i203i304i405i506i6];  // dual(e3) component
214         mv[qc3ga::Ei1] = -primalSurface[qc3ga::E123i102i203i304i405i506i6];
215         mv[qc3ga::Ei2] = -primalSurface[qc3ga::E12301i1i203i304i405i506i6];
216         mv[qc3ga::Ei3] = -primalSurface[qc3ga::E12301i102i2i304i405i506i6];
218         return mv;   
219     }
222     /// \brief compute the dual of a 14-vector into a 1-vector without computing the inner product
223     template<typename T>
224     qc3ga::Mvec<T> primalSurfaceDualization(const qc3ga::Mvec<T> &primalSurface){
225         qc3ga::Mvec<T> mv;
227         mv[qc3ga::E01] = primalSurface[qc3ga::E1230102i203i304i405i506i6]; // dual(e01) component
228         mv[qc3ga::E02] = primalSurface[qc3ga::E12301i10203i304i405i506i6]; // dual(e02) component
229         mv[qc3ga::E03] = primalSurface[qc3ga::E12301i102i20304i405i506i6]; // dual(e03) component
230         mv[qc3ga::E04] = primalSurface[qc3ga::E12301i102i203i30405i506i6]; // dual(e04) component
231         mv[qc3ga::E05] = primalSurface[qc3ga::E12301i102i203i304i40506i6]; // dual(e05) component
232         mv[qc3ga::E06] = primalSurface[qc3ga::E12301i102i203i304i405i506]; // dual(e06) component
235         mv[qc3ga::E1]  = -primalSurface[qc3ga::E2301i102i203i304i405i506i6]; // dual(e1) component 
236         mv[qc3ga::E2]  = primalSurface[qc3ga::E1301i102i203i304i405i506i6]; // dual(e2) component
237         mv[qc3ga::E3]  = -primalSurface[qc3ga::E1201i102i203i304i405i506i6]; // dual(e3) component
239         mv[qc3ga::Ei1] = -primalSurface[qc3ga::E123i102i203i304i405i506i6];
240         mv[qc3ga::Ei2] = -primalSurface[qc3ga::E12301i1i203i304i405i506i6];
241         mv[qc3ga::Ei3] = -primalSurface[qc3ga::E12301i102i2i304i405i506i6];
243         return mv;
244     }
247     /// \brief compute the normal of a surface on a point.
248     /// \param surface is a 14-vector (primal form).
249     /// \param point is a normalized point lying on the surface where the normal is estimated.
250     /// \return a normal vector (e1,e2,e3) with L2 norm = 1 
251     template<typename T>
252     qc3ga::Mvec<T> surfaceNormal(qc3ga::Mvec<T> &surface, qc3ga::Mvec<T> &point){
253     
254       // extract the normal bivector
255       qc3ga::Mvec<T> normalBivector;
256       normalBivector = point ^ (primalSurfaceDualization(surface)) ;
257            
258       // convert the normal bivector in normal vector
259       qc3ga::Mvec<T> normal;
260       normal[qc3ga::E1] = -normalBivector[qc3ga::E101] - normalBivector[qc3ga::E204] - normalBivector[qc3ga::E305];  
261       normal[qc3ga::E2] = -normalBivector[qc3ga::E202] - normalBivector[qc3ga::E104] - normalBivector[qc3ga::E306]; 
262       normal[qc3ga::E3] = -normalBivector[qc3ga::E303] - normalBivector[qc3ga::E105] - normalBivector[qc3ga::E206]; 
263  
264       normal[qc3ga::Ei1] = 0.5*normal[qc3ga::E1]*normal[qc3ga::E1];  
265       normal[qc3ga::Ei2] = 0.5*normal[qc3ga::E2]*normal[qc3ga::E2];  
266       normal[qc3ga::Ei3] = 0.5*normal[qc3ga::E3]*normal[qc3ga::E3];
268       normal[qc3ga::Ei4] = normal[qc3ga::E1]*normal[qc3ga::E2];  
269       normal[qc3ga::Ei5] = normal[qc3ga::E1]*normal[qc3ga::E3];  
270       normal[qc3ga::Ei6] = normal[qc3ga::E2]*normal[qc3ga::E3];  
272       normal[qc3ga::E01] = normal[qc3ga::E02] = normal[qc3ga::E03] = 1.0;  
273         
274       return -normal;
275     }
277     /// \brief compute the normal of a surface on a point.
278     /// \param surface is a 14-vector (primal form).
279     /// \param point is a normalized point lying on the surface where the normal is estimated.
280     /// \return a normal vector (e1,e2,e3), becareful the resulting normal vector is not normalized!! 
281     template<typename T>
282     qc3ga::Mvec<T> surfaceNormalPrimalWithoutWedge(qc3ga::Mvec<T> &surface, qc3ga::Mvec<T> &point){
283     
284       // extract the normal vector of the surface at the point given by point 
285       // qc3ga::Mvec<T> normalBivector;
286       qc3ga::Mvec<T> normal;
287       normal[qc3ga::E1] = -0.5*point[qc3ga::E1]*surface[qc3ga::E01] -point[qc3ga::E2]*surface[qc3ga::E04] -point[qc3ga::E3]*surface[qc3ga::E05] + surface[qc3ga::E1];
288       normal[qc3ga::E2] = -0.5*point[qc3ga::E2]*surface[qc3ga::E02] -point[qc3ga::E1]*surface[qc3ga::E04] -point[qc3ga::E3]*surface[qc3ga::E06] + surface[qc3ga::E2];
289       normal[qc3ga::E3] = -0.5*point[qc3ga::E3]*surface[qc3ga::E03] -point[qc3ga::E1]*surface[qc3ga::E05] -point[qc3ga::E2]*surface[qc3ga::E06] + surface[qc3ga::E3];
290         
291       return -normal;
292   }
295  /// \brief compute the normal of a surface on a point.
296     /// \param surface is a 2-vector (dual form).
297     /// \param point is a normalized point lying on the surface where the normal is estimated.
298     /// \return a normal vector (e1,e2,e3), becareful the resulting normal vector is not normalized!! 
299     template<typename T>
300     qc3ga::Mvec<T> dualSurfaceNormalWithoutWedge(qc3ga::Mvec<T> &surface, qc3ga::Mvec<T> &point){
301     
302       // extract the normal vector of the surface at the point given by point 
303       // qc3ga::Mvec<T> normalBivector;
304       qc3ga::Mvec<T> normal;
305       normal[qc3ga::E1] = -point[qc3ga::E1]*surface[qc3ga::E01] -point[qc3ga::E2]*surface[qc3ga::E04] -point[qc3ga::E3]*surface[qc3ga::E05] + surface[qc3ga::E1];
306       normal[qc3ga::E2] = -point[qc3ga::E2]*surface[qc3ga::E02] -point[qc3ga::E1]*surface[qc3ga::E04] -point[qc3ga::E3]*surface[qc3ga::E06] + surface[qc3ga::E2];
307       normal[qc3ga::E3] = -point[qc3ga::E3]*surface[qc3ga::E03] -point[qc3ga::E1]*surface[qc3ga::E05] -point[qc3ga::E2]*surface[qc3ga::E06] + surface[qc3ga::E3];
308         
309       return -normal;
310     }
316      template<typename T>
317      int quadricLineIntersection(qc3ga::Mvec<T> ptOnLine1, qc3ga::Mvec<T> ptOnLine2, qc3ga::Mvec<T> q, qc3ga::Mvec<T>& pt1, qc3ga::Mvec<T>& pt2, const T epsilon){
319          // extract from the points
320          qc3ga::Mvec<T> m = ptOnLine1-ptOnLine2;
321          m /= std::sqrt(m[qc3ga::E1]*m[qc3ga::E1]+m[qc3ga::E2]*m[qc3ga::E2] + m[qc3ga::E3]*m[qc3ga::E3]); // Normalized support vector of the line
322          ptOnLine2 = ptOnLine1 - m;
324          // constants
325          const T m1   = m[qc3ga::E1];
326          const T m2   = m[qc3ga::E2];
327          const T m3   = m[qc3ga::E3];
328          const T ptl1 = ptOnLine2[qc3ga::E1];
329          const T ptl2 = ptOnLine2[qc3ga::E2];
330          const T ptl3 = ptOnLine2[qc3ga::E3];
331          const T q1   = q[qc3ga::E1];
332          const T q2   = q[qc3ga::E2];
333          const T q3   = q[qc3ga::E3];
334          const T q01  = q[qc3ga::E01];
335          const T q02  = q[qc3ga::E02];
336          const T q03  = q[qc3ga::E03];
337          const T q04  = q[qc3ga::E04];
338          const T q05  = q[qc3ga::E05];
339          const T q06  = q[qc3ga::E06];
340          const T qi1  = q[qc3ga::Ei1];
341          const T qi2  = q[qc3ga::Ei2];
342          const T qi3  = q[qc3ga::Ei3];
346         T a = (( (-0.5*m1*m1*q01)+(-0.5*m2*m2*q02)+(-0.5*m3*m3*q03)
347                 + (-m1*m2*q04) + (-m1*m3*q05) +(-m2*m3*q06)));
349         T c = ((-0.5*ptl1*ptl1*q01)+(-0.5*ptl2*ptl2*q02)+(-0.5*ptl3*ptl3*q03)
350                 +(-ptl1*ptl2*q04)+ (-ptl1*ptl3*q05) +(-ptl2*ptl3*q06) + (ptl1*q1)+ (ptl2*q2)
351                 + (ptl3*q3) + (-(qi1+qi2+qi3)));
353         T b = ((-m1*ptl1*q01)+(-m2*ptl2*q02)+(-m3*ptl3*q03)
354                 + (-(m1*ptl2+m2*ptl1)*q04) + (-(m1*ptl3+m3*ptl1)*q05) 
355                 + (-(m2*ptl3+m3*ptl2)*q06) + (m1*q1 + m2*q2 + m3*q3));
358         // With a epsilon
359         if(std::fabs(a) < epsilon){
360             // Flat point
361             pt1 = m*(-c/b) + ptOnLine2;
362             //pt2 = pt1;
363             return 1;
364         }
366         T squaredDist = (b*b) -(4.0*a*c);
368         if(squaredDist< 0.0)
369             return 0; // complex solutions
371         // Should do it with an epsilon
372         if(std::fabs(squaredDist) < epsilon){
373             // Flat point
374             pt1 =(m*(-b/(2.0*(a))))+ptOnLine2;
375             //pt2 = pt1;
376             return 1;
377         }
379         pt1 =(m*((-b+std::sqrt(squaredDist))/(2.0*a)))+ptOnLine2;
380         pt2 =(m*((-b-std::sqrt(squaredDist))/(2.0*a)))+ptOnLine2;
382         return 2;
383      }
387     template<typename T>
388     void quadricValueDisplay(qc3ga::Mvec<T> mv, const double epsilon){
390         // remove zeros
391         mv.roundZero(epsilon);
393         // check unhomogeneous multivectors
394         if(mv.grades().size() != 1)
395             std::cout << "warning: quadricValueDisplay: not homogeneous multivector" << std::endl;
397         // convert to dual quadric
398         qc3ga::Mvec<T> dualQaudric;
399         switch(mv.grade()) {
400             case 1: // dual quadric
401                 dualQaudric = mv;
402                 break;
403             case 14: // primal quadric
404                 dualQaudric = mv.dual();
405                 break;
406             default: // not a quadric
407                 std::cout << "this multivector is not a quadric:\n" << std::endl;
408                 break;
409         }
411         // data to display
412         std::string quadricComponents;    // like "ax^2 + ..."
413         std::string quadricCoefficients;  // like "a=42 ..."
416         // extract the quadric components
417         T component;
419         // a : x^2
420         component = dualQaudric[qc3ga::E01] / -2.0;
421         if(fabs(component) > epsilon){
422             // handle the sign
423             std::string sign;
424             if(component < 0) sign = "- ";
425             else if(quadricComponents.size() > 0) sign = "+ ";
426             quadricComponents   += " + a.x^2";
427             quadricCoefficients += "a = " + sign + std::to_string(fabs(component)) + "\n";
428         }
430         // b : y^2
431         component = dualQaudric[qc3ga::E02] / -2.0;
432         if(fabs(component) > epsilon){
433             // handle the sign
434             std::string sign;
435             if(component < 0) sign = "- ";
436             else if(quadricComponents.size() > 0) sign = "+ ";
437             quadricComponents   += " + b.y^2";
438             quadricCoefficients += "b = " + sign + std::to_string(fabs(component)) + "\n";
439         }
441         // c : z^2
442         component = dualQaudric[qc3ga::E03] / -2.0;
443         if(fabs(component) > epsilon){
444             // handle the sign
445             std::string sign;
446             if(component < 0) sign = "- ";
447             else if(quadricComponents.size() > 0) sign = "+ ";
448             quadricComponents   += " + c.z^2";
449             quadricCoefficients += "c = " + sign + std::to_string(fabs(component)) + "\n";
450         }
452         // d : xy
453         component = - dualQaudric[qc3ga::E04];
454         if(fabs(component) > epsilon){
455             // handle the sign
456             std::string sign;
457             if(component < 0) sign = "- ";
458             else if(quadricComponents.size() > 0) sign = "+ ";
459             quadricComponents   += " + d.xy";
460             quadricCoefficients += "d = " + sign + std::to_string(fabs(component)) + "\n";
461         }
463         // e : xz
464         component = - dualQaudric[qc3ga::E05];
465         if(fabs(component) > epsilon){
466             // handle the sign
467             std::string sign;
468             if(component < 0) sign = "- ";
469             else if(quadricComponents.size() > 0) sign = "+ ";
470             quadricComponents   += " + e.xz";
471             quadricCoefficients += "e = " + sign + std::to_string(fabs(component)) + "\n";
472         }
474         // f : yz
475         component = - dualQaudric[qc3ga::E06];
476         if(fabs(component) > epsilon){
477             // handle the sign
478             std::string sign;
479             if(component < 0) sign = "- ";
480             else if(quadricComponents.size() > 0) sign = "+ ";
481             quadricComponents   += " + f.yz";
482             quadricCoefficients += "f = " + sign + std::to_string(fabs(component)) + "\n";
483         }
485         // g : x
486         component = dualQaudric[qc3ga::E1];
487         if(fabs(component) > epsilon){
488             // handle the sign
489             std::string sign;
490             if(component < 0) sign = "- ";
491             else if(quadricComponents.size() > 0) sign = "+ ";
492             quadricComponents   += " + g.x";
493             quadricCoefficients += "g = " + sign + std::to_string(fabs(component)) + "\n";
494         }
496         // h : y
497         component = dualQaudric[qc3ga::E2];
498         if(fabs(component) > epsilon){
499             // handle the sign
500             std::string sign;
501             if(component < 0) sign = "- ";
502             else if(quadricComponents.size() > 0) sign = "+ ";
503             quadricComponents   += " + h.y";
504             quadricCoefficients += "h = " + sign + std::to_string(fabs(component)) + "\n";
505         }
507         // i : z
508         component = dualQaudric[qc3ga::E3];
509         if(fabs(component) > epsilon){
510             // handle the sign
511             std::string sign;
512             if(component < 0) sign = "- ";
513             else if(quadricComponents.size() > 0) sign = "+ ";
514             quadricComponents   += " + i.z";
515             quadricCoefficients += "i = " + sign + std::to_string(fabs(component)) + "\n";
516         }
518         // j : constant
519         component = - (dualQaudric[qc3ga::Ei1] + dualQaudric[qc3ga::Ei2] + dualQaudric[qc3ga::Ei3]);
520         if(fabs(component) > epsilon){
521             // handle the sign
522             std::string sign;
523             if(component < 0) sign = "- ";
524             else if(quadricComponents.size() > 0) sign = "+ ";
525             quadricComponents   += " + j";
526             quadricCoefficients += "j = " + sign + std::to_string(fabs(component)) + "\n";
527         }
529         // matrices delta : http://mathworld.wolfram.com/QuadraticSurface.html
530         // delta 1 (rank)
531         Eigen::Matrix3d delta1;
532         delta1 << dualQaudric[qc3ga::E01] / -2.0, dualQaudric[qc3ga::E04] / -2.0, dualQaudric[qc3ga::E05] / -2.0,
533                   dualQaudric[qc3ga::E04] / -2.0, dualQaudric[qc3ga::E02] / -2.0, dualQaudric[qc3ga::E06] / -2.0,
534                   dualQaudric[qc3ga::E05] / -2.0, dualQaudric[qc3ga::E06] / -2.0, dualQaudric[qc3ga::E03] / -2.0;
535         Eigen::FullPivLU<Eigen::Matrix3d> lu_decomp3(delta1);
536         quadricCoefficients += "rho3    = " + std::to_string(lu_decomp3.rank()) + "\n";
538         // delta 2 (rank and sign of det)
539         Eigen::Matrix4d delta2;
540         delta2 << dualQaudric[qc3ga::E01] / -2.0, dualQaudric[qc3ga::E04] / -2.0, dualQaudric[qc3ga::E05] / -2.0, dualQaudric[qc3ga::E1] / 2.0,
541                   dualQaudric[qc3ga::E04] / -2.0, dualQaudric[qc3ga::E02] / -2.0, dualQaudric[qc3ga::E06] / -2.0, dualQaudric[qc3ga::E2] / 2.0,
542                   dualQaudric[qc3ga::E05] / -2.0, dualQaudric[qc3ga::E06] / -2.0, dualQaudric[qc3ga::E03] / -2.0, dualQaudric[qc3ga::E3] / 2.0,
543                   dualQaudric[qc3ga::E1]  /  2.0, dualQaudric[qc3ga::E2]  /  2.0, dualQaudric[qc3ga::E3]  /  2.0, -(dualQaudric[qc3ga::Ei1] + dualQaudric[qc3ga::Ei2] + dualQaudric[qc3ga::Ei3]) * 3.0;
544         Eigen::FullPivLU<Eigen::Matrix4d> lu_decomp4(delta2);
545         quadricCoefficients += "rho4    = " + std::to_string(lu_decomp4.rank()) + "\n";
546         quadricCoefficients += "sign(d) = ";
547         if(delta2.determinant()>0) quadricCoefficients += "+";
548         if(delta2.determinant()<0) quadricCoefficients += "-";
549         quadricCoefficients += "\n";
551         // delta1 eigen value sign
552         Eigen::EigenSolver<Eigen::Matrix3d> es(delta1,false);
553         std::vector<bool> eigenValuesSignPositive;
554         if(es.eigenvalues().real()(0) != 0){
555             if(es.eigenvalues().real()(0) > 0) eigenValuesSignPositive.push_back(true);
556             else eigenValuesSignPositive.push_back(false);
557         }
558         if(es.eigenvalues().real()(1) != 0){
559             if(es.eigenvalues().real()(1) > 0) eigenValuesSignPositive.push_back(true);
560             else eigenValuesSignPositive.push_back(false);
561         }
562         if(es.eigenvalues().real()(2) != 0){
563             if(es.eigenvalues().real()(2) > 0) eigenValuesSignPositive.push_back(true);
564             else eigenValuesSignPositive.push_back(false);
565         }
566         int k = 1;
567         if(eigenValuesSignPositive.size() > 1)
568             for(unsigned int i=1; i<eigenValuesSignPositive.size(); ++i)
569                 if(eigenValuesSignPositive[i] != eigenValuesSignPositive[0])
570                     k = 0;
572         quadricCoefficients += "k       = " + std::to_string(k) + "\n";
574         // print the result
575         std::cout << "\ncurrent quadric : " << quadricComponents << " = 0" << std::endl;
576         std::cout << quadricCoefficients << std::endl;
577     }
580 } // namespace
582 #endif // projection_inclusion_guard