1 // Copyright (c) 2006  GeometryFactory (France). All rights reserved.
2 //
3 // This file is part of CGAL (www.cgal.org).
4 //
5 // $URL: https://github.com/CGAL/cgal/blob/v5.3/Surface_mesh_simplification/include/CGAL/Cartesian/MatrixC33.h $
6 // $Id: MatrixC33.h ff09c5d 2019-10-25T16:35:53+02:00 Mael Rouxel-Labbé
7 // SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial
8 //
9 // Author(s)     : Fernando Cacciola <fernando.cacciola@geometryfactory.com>
10 //
11 #ifndef CGAL_CARTESIAN_MATRIXC33_H
12 #define CGAL_CARTESIAN_MATRIXC33_H
13 
14 #include <CGAL/license/Surface_mesh_simplification.h>
15 
16 #include <CGAL/determinant.h>
17 #include <CGAL/Null_matrix.h>
18 #include <CGAL/number_utils.h>
19 #include <CGAL/Vector_3.h>
20 
21 #include <boost/optional/optional.hpp>
22 
23 namespace CGAL {
24 
25 template <class R_>
26 class MatrixC33
27 {
28 public:
29   typedef R_                                                R;
30 
31   typedef typename R::FT                                    RT;
32   typedef typename R::Vector_3                              Vector_3;
33 
MatrixC33(Null_matrix)34   MatrixC33(Null_matrix)
35    : mR0(NULL_VECTOR),
36      mR1(NULL_VECTOR),
37      mR2(NULL_VECTOR)
38   {}
39 
MatrixC33(const RT & r0x,const RT & r0y,const RT & r0z,const RT & r1x,const RT & r1y,const RT & r1z,const RT & r2x,const RT & r2y,const RT & r2z)40   MatrixC33(const RT& r0x, const RT& r0y, const RT& r0z,
41             const RT& r1x, const RT& r1y, const RT& r1z,
42             const RT& r2x, const RT& r2y, const RT& r2z)
43     : mR0(r0x,r0y,r0z),
44       mR1(r1x,r1y,r1z),
45       mR2(r2x,r2y,r2z)
46   {}
47 
MatrixC33(const Vector_3 & r0,const Vector_3 & r1,const Vector_3 & r2)48   MatrixC33(const Vector_3& r0, const Vector_3& r1, const Vector_3& r2)
49     : mR0(r0),
50       mR1(r1),
51       mR2(r2)
52   {}
53 
r0()54   const Vector_3& r0() const { return mR0; }
r1()55   const Vector_3& r1() const { return mR1; }
r2()56   const Vector_3& r2() const { return mR2; }
57 
r0()58   Vector_3& r0() { return mR0; }
r1()59   Vector_3& r1() { return mR1; }
r2()60   Vector_3& r2() { return mR2; }
61 
62   const Vector_3& operator[](int row) const { return row == 0 ? mR0 : (row == 1 ? mR1 : mR2); }
63   Vector_3& operator[](int row) { return row == 0 ? mR0 : (row == 1 ? mR1 : mR2); }
64 
65   MatrixC33& operator+=(const MatrixC33& m)
66   {
67     mR0 = mR0 + m.r0();
68     mR1 = mR1 + m.r1();
69     mR2 = mR2 + m.r2();
70     return *this;
71   }
72 
73   MatrixC33& operator-=(const MatrixC33& m)
74   {
75     mR0 = mR0 - m.r0();
76     mR1 = mR1 - m.r1();
77     mR2 = mR2 - m.r2();
78     return *this;
79   }
80 
81   MatrixC33& operator*=(const RT& c)
82   {
83     mR0 = mR0 * c;
84     mR1 = mR1 * c;
85     mR2 = mR2 * c;
86     return *this;
87   }
88 
89   MatrixC33& operator/=(const RT& c)
90   {
91     mR0 = mR0 / c;
92     mR1 = mR1 / c;
93     mR2 = mR2 / c;
94     return *this;
95   }
96 
97   friend MatrixC33 operator+(const MatrixC33& a, const MatrixC33& b)
98   {
99     return MatrixC33(a.r0() + b.r0(),
100                      a.r1() + b.r1(),
101                      a.r2() + b.r2());
102   }
103 
104   friend MatrixC33 operator-(const MatrixC33& a, const MatrixC33& b)
105   {
106     return MatrixC33(a.r0() - b.r0(),
107                      a.r1() - b.r1(),
108                      a.r2() - b.r2());
109   }
110 
111   friend MatrixC33 operator*(const MatrixC33& m, const RT& c)
112   {
113     return MatrixC33(m.r0()*c, m.r1()*c, m.r2()*c);
114   }
115   friend MatrixC33 operator*(const RT& c, const MatrixC33& m)
116   {
117     return MatrixC33(m.r0()*c, m.r1()*c, m.r2()*c);
118   }
119 
120   friend MatrixC33 operator/(const MatrixC33& m, const RT& c)
121   {
122     return MatrixC33(m.r0()/c, m.r1()/c, m.r2()/c);
123   }
124 
125   friend Vector_3 operator*(const MatrixC33& m, const Vector_3& v)
126   {
127     return Vector_3(m.r0()*v, m.r1()*v, m.r2()*v);
128   }
129   friend Vector_3 operator*(const Vector_3& v, const MatrixC33& m)
130   {
131     return Vector_3(v*m.r0(), v*m.r1(), v*m.r2());
132   }
133 
determinant()134   RT determinant() const
135   {
136     return CGAL::determinant(r0().x(), r0().y(), r0().z(),
137                              r1().x(), r1().y(), r1().z(),
138                              r2().x(), r2().y(), r2().z());
139   }
140 
transpose()141   MatrixC33& transpose()
142   {
143     mR0 = Vector_3(r0().x(),r1().x(),r2().x());
144     mR1 = Vector_3(r0().y(),r1().y(),r2().y());
145     mR2 = Vector_3(r0().z(),r1().z(),r2().z());
146     return *this;
147   }
148 
149 private:
150 
151   Vector_3 mR0;
152   Vector_3 mR1;
153   Vector_3 mR2;
154 };
155 
156 template<class R>
157 inline
direct_product(const Vector_3<R> & u,const Vector_3<R> & v)158 MatrixC33<R> direct_product(const Vector_3<R>& u,
159                             const Vector_3<R>& v)
160 {
161   return MatrixC33<R>(v * u.x(),
162                       v * u.y(),
163                       v * u.z());
164 }
165 
166 template<class R>
transposed_matrix(const MatrixC33<R> & m)167 MatrixC33<R> transposed_matrix(const MatrixC33<R>& m)
168 {
169   MatrixC33<R> copy = m;
170   copy.Transpose();
171   return copy;
172 }
173 
174 template<class R>
cofactors_matrix(const MatrixC33<R> & m)175 MatrixC33<R> cofactors_matrix(const MatrixC33<R>& m)
176 {
177   typedef typename R::RT RT;
178 
179   RT c00 =  determinant(m.r1().y(),m.r1().z(),m.r2().y(),m.r2().z());
180   RT c01 = -determinant(m.r1().x(),m.r1().z(),m.r2().x(),m.r2().z());
181   RT c02 =  determinant(m.r1().x(),m.r1().y(),m.r2().x(),m.r2().y());
182 
183   RT c10 = -determinant(m.r0().y(),m.r0().z(),m.r2().y(),m.r2().z());
184   RT c11 =  determinant(m.r0().x(),m.r0().z(),m.r2().x(),m.r2().z());
185   RT c12 = -determinant(m.r0().x(),m.r0().y(),m.r2().x(),m.r2().y());
186 
187   RT c20 =  determinant(m.r0().y(),m.r0().z(),m.r1().y(),m.r1().z());
188   RT c21 = -determinant(m.r0().x(),m.r0().z(),m.r1().x(),m.r1().z());
189   RT c22 =  determinant(m.r0().x(),m.r0().y(),m.r1().x(),m.r1().y());
190 
191   return MatrixC33<R>(c00,c01,c02,
192                       c10,c11,c12,
193                       c20,c21,c22);
194 }
195 
196 template<class R>
adjoint_matrix(const MatrixC33<R> & m)197 MatrixC33<R> adjoint_matrix(const MatrixC33<R>& m)
198 {
199   return cofactors_matrix(m).transpose();
200 }
201 
202 template<class R>
inverse_matrix(const MatrixC33<R> & m)203 boost::optional< MatrixC33<R> > inverse_matrix(const MatrixC33<R>& m)
204 {
205   typedef typename R::RT                                        RT;
206   typedef MatrixC33<R>                                          Matrix;
207   typedef boost::optional<Matrix>                               result_type;
208 
209   result_type rInverse;
210 
211   RT det = m.determinant();
212 
213   if(! CGAL_NTS is_zero(det))
214   {
215     RT c00 = (m.r1().y()*m.r2().z() - m.r1().z()*m.r2().y()) / det;
216     RT c01 = (m.r2().y()*m.r0().z() - m.r0().y()*m.r2().z()) / det;
217     RT c02 = (m.r0().y()*m.r1().z() - m.r1().y()*m.r0().z()) / det;
218 
219     RT c10 = (m.r1().z()*m.r2().x() - m.r1().x()*m.r2().z()) / det;
220     RT c11 = (m.r0().x()*m.r2().z() - m.r2().x()*m.r0().z()) / det;
221     RT c12 = (m.r1().x()*m.r0().z() - m.r0().x()*m.r1().z()) / det;
222 
223     RT c20 = (m.r1().x()*m.r2().y() - m.r2().x()*m.r1().y()) / det;
224     RT c21 = (m.r2().x()*m.r0().y() - m.r0().x()*m.r2().y()) / det;
225     RT c22 = (m.r0().x()*m.r1().y() - m.r0().y()*m.r1().x()) / det;
226 
227     rInverse = result_type(Matrix(c00,c01,c02,
228                                   c10,c11,c12,
229                                   c20,c21,c22));
230   }
231 
232   return rInverse;
233 }
234 
235 } // namespace CGAL
236 
237 #endif // CGAL_CARTESIAN_MATRIXC33_H //
238 // EOF //
239 
240 
241