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