MeshLib C++ Docs
Loading...
Searching...
No Matches
MRMatrix4.h
Go to the documentation of this file.
1#pragma once
2
3#include "MRVector4.h"
4#include <limits>
5#include <cassert>
6
7namespace MR
8{
9
10#ifdef _MSC_VER
11#pragma warning(push)
12#pragma warning(disable: 4804) // unsafe use of type 'bool' in operation
13#pragma warning(disable: 4146) // unary minus operator applied to unsigned type, result still unsigned
14#endif
15
18template <typename T>
19struct Matrix4
20{
21 using ValueType = T;
23
25 Vector4<T> x{ 1, 0, 0, 0 };
26 Vector4<T> y{ 0, 1, 0, 0 };
27 Vector4<T> z{ 0, 0, 1, 0 };
28 Vector4<T> w{ 0, 0, 0, 1 };
29
30 constexpr Matrix4() noexcept = default;
32 constexpr Matrix4( const Vector4<T>& x, const Vector4<T>& y, const Vector4<T>& z, const Vector4<T>& w ) : x( x ), y( y ), z( z ), w( w ) { }
33
35 constexpr Matrix4( const Matrix3<T>& r, const Vector3<T>& t )
36 {
37 x = Vector4<T>( r.x.x, r.x.y, r.x.z, t.x );
38 y = Vector4<T>( r.y.x, r.y.y, r.y.z, t.y );
39 z = Vector4<T>( r.z.x, r.z.y, r.z.z, t.z );
40 w = Vector4<T>( 0, 0, 0, 1 );
41 }
42
43 // Currently `AffineXf3<long long>` doesn't seem to compile, so we disable this constructor for `Matrix4<long long>`, because otherwise
44 // mrbind instantiates the entire `AffineXf3<long long>` and chokes on it.
45 template <MR_SAME_TYPE_TEMPLATE_PARAM(T, TT)>
46 constexpr Matrix4( const AffineXf3<TT>& xf ) MR_REQUIRES_IF_SUPPORTED( std::floating_point<T> ) : Matrix4( xf.A, xf.b ) {}
47
48 template <typename U>
49 constexpr explicit Matrix4( const Matrix4<U> & m ) : x( m.x ), y( m.y ), z( m.z ), w( m.w ) { }
50 static constexpr Matrix4 zero() noexcept { return Matrix4( Vector4<T>(), Vector4<T>(), Vector4<T>(), Vector4<T>() ); }
51 static constexpr Matrix4 identity() noexcept { return Matrix4(); }
53 static constexpr Matrix4 scale( T s ) noexcept { return Matrix4( { s, T(0), T(0), T(0) }, { T(0), s, T(0), T(0) }, { T(0), T(0), s, T(0) }, { T(0), T(0), T(0), s } ); }
54
56 constexpr const T& operator ()( int row, int col ) const noexcept { return operator[]( row )[col]; }
57 constexpr T& operator ()( int row, int col ) noexcept { return operator[]( row )[col]; }
58
60 constexpr const Vector4<T> & operator []( int row ) const noexcept { return *( &x + row ); }
61 constexpr Vector4<T> & operator []( int row ) noexcept { return *( &x + row ); }
62
64 constexpr Vector4<T> col( int i ) const noexcept { return { x[i], y[i], z[i], w[i] }; }
65
67 constexpr T trace() const noexcept { return x.x + y.y + z.z + w.w; }
69 constexpr T normSq() const noexcept { return x.lengthSq() + y.lengthSq() + z.lengthSq() + w.lengthSq(); }
70 constexpr auto norm() const noexcept
71 {
72 // Calling `sqrt` this way to hopefully support boost.multiprecision numbers.
73 // Returning `auto` to not break on integral types.
74 using std::sqrt;
75 return sqrt( normSq() );
76 }
78 Matrix3<T> submatrix3( int i, int j ) const noexcept;
80 T det() const noexcept;
82 constexpr Matrix4<T> inverse() const noexcept MR_REQUIRES_IF_SUPPORTED( !std::is_integral_v<T> );
84 constexpr Matrix4<T> transposed() const noexcept;
85
86 constexpr Matrix3<T> getRotation() const noexcept;
87 void setRotation( const Matrix3<T>& rot) noexcept;
88 constexpr Vector3<T> getTranslation() const noexcept;
89 void setTranslation( const Vector3<T>& t ) noexcept;
90
91 constexpr T* data() { return (T*) (&x); };
92 constexpr const T* data() const { return (T*) (&x); };
93
94 template <MR_SAME_TYPE_TEMPLATE_PARAM(T, TT)>
95 operator AffineXf3<TT>() const MR_REQUIRES_IF_SUPPORTED( std::floating_point<T> )
96 {
97 assert( std::abs( w.x ) < std::numeric_limits<T>::epsilon() * 1000 );
98 assert( std::abs( w.y ) < std::numeric_limits<T>::epsilon() * 1000 );
99 assert( std::abs( w.z ) < std::numeric_limits<T>::epsilon() * 1000 );
100 assert( std::abs( 1 - w.w ) < std::numeric_limits<T>::epsilon() * 1000 );
101 AffineXf3<T> res;
102 res.A.x.x = x.x; res.A.x.y = x.y; res.A.x.z = x.z; res.b.x = x.w;
103 res.A.y.x = y.x; res.A.y.y = y.y; res.A.y.z = y.z; res.b.y = y.w;
104 res.A.z.x = z.x; res.A.z.y = z.y; res.A.z.z = z.z; res.b.z = z.w;
105 return res;
106 }
107
110 Vector3<T> operator ()( const Vector3<T> & b ) const MR_REQUIRES_IF_SUPPORTED( !std::is_integral_v<T> );
111
112 [[nodiscard]] friend constexpr bool operator ==( const Matrix4<T> & a, const Matrix4<T> & b ) { return a.x == b.x && a.y == b.y && a.z == b.z && a.w == b.w; }
113 [[nodiscard]] friend constexpr bool operator !=( const Matrix4<T> & a, const Matrix4<T> & b ) { return !( a == b ); }
114
115 // NOTE: We use `std::declval()` in the operators below because libclang 18 in our binding generator is bugged and chokes on decltyping `a.x` and such. TODO fix this when we update libclang.
116
117 [[nodiscard]] friend constexpr auto operator +( const Matrix4<T> & a, const Matrix4<T> & b ) -> Matrix4<decltype( std::declval<T>() + std::declval<T>() )> { return { a.x + b.x, a.y + b.y, a.z + b.z, a.w + b.w }; }
118 [[nodiscard]] friend constexpr auto operator -( const Matrix4<T> & a, const Matrix4<T> & b ) -> Matrix4<decltype( std::declval<T>() - std::declval<T>() )> { return { a.x - b.x, a.y - b.y, a.z - b.z, a.w - b.w }; }
119 [[nodiscard]] friend constexpr auto operator *( T a, const Matrix4<T> & b ) -> Matrix4<decltype( std::declval<T>() * std::declval<T>() )> { return { a * b.x, a * b.y, a * b.z, a * b.w }; }
120 [[nodiscard]] friend constexpr auto operator *( const Matrix4<T> & b, T a ) -> Matrix4<decltype( std::declval<T>() * std::declval<T>() )> { return { a * b.x, a * b.y, a * b.z, a * b.z }; }
121 [[nodiscard]] friend constexpr auto operator /( Matrix4<T> b, T a ) -> Matrix4<decltype( std::declval<T>() / std::declval<T>() )>
122 {
123 if constexpr ( std::is_integral_v<T> )
124 return { b.x / a, b.y / a, b.z / a, b.w / a };
125 else
126 return b * ( 1 / a );
127 }
128
129 friend constexpr Matrix4<T> & operator +=( Matrix4<T> & a, const Matrix4<T> & b ) MR_REQUIRES_IF_SUPPORTED( requires{ a + b; } ) { a.x += b.x; a.y += b.y; a.z += b.z; a.w += b.w; return a; }
130 friend constexpr Matrix4<T> & operator -=( Matrix4<T> & a, const Matrix4<T> & b ) MR_REQUIRES_IF_SUPPORTED( requires{ a - b; } ) { a.x -= b.x; a.y -= b.y; a.z -= b.z; a.w -= b.w; return a; }
131 friend constexpr Matrix4<T> & operator *=( Matrix4<T> & a, T b ) MR_REQUIRES_IF_SUPPORTED( requires{ a * b; } ) { a.x *= b; a.y *= b; a.z *= b; a.w *= b; return a; }
132 friend constexpr Matrix4<T> & operator /=( Matrix4<T> & a, T b ) MR_REQUIRES_IF_SUPPORTED( requires{ a / b; } )
133 {
134 if constexpr ( std::is_integral_v<T> )
135 { a.x /= b; a.y /= b; a.z /= b; a.w /= b; return a; }
136 else
137 return a *= ( 1 / b );
138 }
139
141 [[nodiscard]] friend constexpr auto operator *( const Matrix4<T> & a, const Vector4<T> & b ) -> Vector4<decltype( dot( std::declval<Vector4<T>>(), std::declval<Vector4<T>>() ) )>
142 {
143 return { dot( a.x, b ), dot( a.y, b ), dot( a.z, b ), dot( a.w, b ) };
144 }
145
147 [[nodiscard]] friend constexpr auto operator *( const Matrix4<T> & a, const Matrix4<T> & b ) -> Matrix4<decltype( dot( std::declval<Vector4<T>>(), std::declval<Vector4<T>>() ) )>
148 {
149 Matrix4<decltype( dot( std::declval<Vector4<T>>(), std::declval<Vector4<T>>() ) )> res;
150 for ( int i = 0; i < 4; ++i )
151 for ( int j = 0; j < 4; ++j )
152 res[i][j] = dot( a[i], b.col(j) );
153 return res;
154 }
155};
156
159
161template <typename T>
162inline auto dot( const Matrix4<T> & a, const Matrix4<T> & b ) -> decltype( dot( a.x, b.x ) )
163{
164 return dot( a.x, b.x ) + dot( a.y, b.y ) + dot( a.z, b.z ) + dot( a.w, b.w );
165}
166
167template <typename T>
168inline Vector3<T> Matrix4<T>::operator ()( const Vector3<T> & b ) const MR_REQUIRES_IF_SUPPORTED( !std::is_integral_v<T> )
169{
170 return ( *this * Vector4<T>{ b.x, b.y, b.z, T(1) } ).proj3d();
171}
172
174template <typename T>
175inline Matrix4<T> outer( const Vector4<T> & a, const Vector4<T> & b )
176{
177 return { a.x * b, a.y * b, a.z * b, a.w * b };
178}
179
180template <typename T>
181Matrix3<T> Matrix4<T>::submatrix3( int i, int j ) const noexcept
182{
183 Matrix3<T> res;
184 auto* resM = (T*) &res.x;
185 int cur = 0;
186 for ( int m = 0; m < 4; m++ )
187 {
188 if ( m == i )
189 continue;
190 for ( int n = 0; n < 4; n++ )
191 {
192 if ( n == j )
193 continue;
194 resM[cur++] = (*this)[m][n];
195 }
196 }
197 assert( cur == 9 );
198 return res;
199}
200
201template <typename T>
202T Matrix4<T>::det() const noexcept
203{
204 return
205 x.x * submatrix3( 0, 0 ).det()
206 - x.y * submatrix3( 0, 1 ).det()
207 + x.z * submatrix3( 0, 2 ).det()
208 - x.w * submatrix3( 0, 3 ).det();
209}
210
211template <typename T>
212constexpr Matrix4<T> Matrix4<T>::transposed() const noexcept
213{
214 return Matrix4<T>
215 {
216 { x.x, y.x, z.x, w.x },
217 { x.y, y.y, z.y, w.y },
218 { x.z, y.z, z.z, w.z },
219 { x.w, y.w, z.w, w.w },
220 };
221}
222
223template <typename T>
224constexpr Matrix4<T> Matrix4<T>::inverse() const noexcept MR_REQUIRES_IF_SUPPORTED( !std::is_integral_v<T> )
225{
226 Matrix4<T> inv;
227 T* m = (T*) (&x);
228 T det;
229
230 inv[0][0] = m[5] * m[10] * m[15] -
231 m[5] * m[11] * m[14] -
232 m[9] * m[6] * m[15] +
233 m[9] * m[7] * m[14] +
234 m[13] * m[6] * m[11] -
235 m[13] * m[7] * m[10];
236
237 inv[1][0] = -m[4] * m[10] * m[15] +
238 m[4] * m[11] * m[14] +
239 m[8] * m[6] * m[15] -
240 m[8] * m[7] * m[14] -
241 m[12] * m[6] * m[11] +
242 m[12] * m[7] * m[10];
243
244 inv[2][0] = m[4] * m[9] * m[15] -
245 m[4] * m[11] * m[13] -
246 m[8] * m[5] * m[15] +
247 m[8] * m[7] * m[13] +
248 m[12] * m[5] * m[11] -
249 m[12] * m[7] * m[9];
250
251 inv[3][0] = -m[4] * m[9] * m[14] +
252 m[4] * m[10] * m[13] +
253 m[8] * m[5] * m[14] -
254 m[8] * m[6] * m[13] -
255 m[12] * m[5] * m[10] +
256 m[12] * m[6] * m[9];
257
258 inv[0][1] = -m[1] * m[10] * m[15] +
259 m[1] * m[11] * m[14] +
260 m[9] * m[2] * m[15] -
261 m[9] * m[3] * m[14] -
262 m[13] * m[2] * m[11] +
263 m[13] * m[3] * m[10];
264
265 inv[1][1] = m[0] * m[10] * m[15] -
266 m[0] * m[11] * m[14] -
267 m[8] * m[2] * m[15] +
268 m[8] * m[3] * m[14] +
269 m[12] * m[2] * m[11] -
270 m[12] * m[3] * m[10];
271
272 inv[2][1] = -m[0] * m[9] * m[15] +
273 m[0] * m[11] * m[13] +
274 m[8] * m[1] * m[15] -
275 m[8] * m[3] * m[13] -
276 m[12] * m[1] * m[11] +
277 m[12] * m[3] * m[9];
278
279 inv[3][1] = m[0] * m[9] * m[14] -
280 m[0] * m[10] * m[13] -
281 m[8] * m[1] * m[14] +
282 m[8] * m[2] * m[13] +
283 m[12] * m[1] * m[10] -
284 m[12] * m[2] * m[9];
285
286 inv[0][2] = m[1] * m[6] * m[15] -
287 m[1] * m[7] * m[14] -
288 m[5] * m[2] * m[15] +
289 m[5] * m[3] * m[14] +
290 m[13] * m[2] * m[7] -
291 m[13] * m[3] * m[6];
292
293 inv[1][2] = -m[0] * m[6] * m[15] +
294 m[0] * m[7] * m[14] +
295 m[4] * m[2] * m[15] -
296 m[4] * m[3] * m[14] -
297 m[12] * m[2] * m[7] +
298 m[12] * m[3] * m[6];
299
300 inv[2][2] = m[0] * m[5] * m[15] -
301 m[0] * m[7] * m[13] -
302 m[4] * m[1] * m[15] +
303 m[4] * m[3] * m[13] +
304 m[12] * m[1] * m[7] -
305 m[12] * m[3] * m[5];
306
307 inv[3][2] = -m[0] * m[5] * m[14] +
308 m[0] * m[6] * m[13] +
309 m[4] * m[1] * m[14] -
310 m[4] * m[2] * m[13] -
311 m[12] * m[1] * m[6] +
312 m[12] * m[2] * m[5];
313
314 inv[0][3] = -m[1] * m[6] * m[11] +
315 m[1] * m[7] * m[10] +
316 m[5] * m[2] * m[11] -
317 m[5] * m[3] * m[10] -
318 m[9] * m[2] * m[7] +
319 m[9] * m[3] * m[6];
320
321 inv[1][3] = m[0] * m[6] * m[11] -
322 m[0] * m[7] * m[10] -
323 m[4] * m[2] * m[11] +
324 m[4] * m[3] * m[10] +
325 m[8] * m[2] * m[7] -
326 m[8] * m[3] * m[6];
327
328 inv[2][3] = -m[0] * m[5] * m[11] +
329 m[0] * m[7] * m[9] +
330 m[4] * m[1] * m[11] -
331 m[4] * m[3] * m[9] -
332 m[8] * m[1] * m[7] +
333 m[8] * m[3] * m[5];
334
335 inv[3][3] = m[0] * m[5] * m[10] -
336 m[0] * m[6] * m[9] -
337 m[4] * m[1] * m[10] +
338 m[4] * m[2] * m[9] +
339 m[8] * m[1] * m[6] -
340 m[8] * m[2] * m[5];
341
342 det = m[0] * inv[0][0] + m[1] * inv[1][0] + m[2] * inv[2][0] + m[3] * inv[3][0];
343
344 if( det == 0 )
345 {
346 // impossible to invert singular matrix
347 assert( false );
348 return Matrix4<T>();
349 }
350
351 inv /= det;
352
353 return inv;
354}
355
356template <typename T>
357constexpr Matrix3<T> Matrix4<T>::getRotation() const noexcept
358{
359 return Matrix3<T>
360 {
361 { x.x, x.y, x.z },
362 { y.x, y.y, y.z },
363 { z.x, z.y, z.z }
364 };
365}
366
367template <typename T>
368void Matrix4<T>::setRotation( const Matrix3<T>& rot ) noexcept
369{
370 x.x = rot.x.x; x.y = rot.x.y; x.z = rot.x.z;
371 y.x = rot.y.x; y.y = rot.y.y; y.z = rot.y.z;
372 z.x = rot.z.x; z.y = rot.z.y; z.z = rot.z.z;
373}
374
375template <typename T>
376constexpr Vector3<T> Matrix4<T>::getTranslation() const noexcept
377{
378 return Vector3<T>{ x.w, y.w, z.w };
379}
380
381template <typename T>
382void Matrix4<T>::setTranslation( const Vector3<T>& t ) noexcept
383{
384 x.w = t.x; y.w = t.y; z.w = t.z;
385}
386
388
389#ifdef _MSC_VER
390#pragma warning(pop)
391#endif
392
393} // namespace MR
#define MR_REQUIRES_IF_SUPPORTED(...)
Definition MRMacros.h:31
MRMESH_CLASS Vector3< double > Matrix2< double > Matrix4
Definition MRMesh/MRMeshFwd.h:202
MRMESH_CLASS Vector3
Definition MRMesh/MRMeshFwd.h:170
AffineXf< Vector3< T > > AffineXf3
Definition MRMesh/MRMeshFwd.h:241
Definition MRMesh/MRAffineXf.h:14
V b
Definition MRMesh/MRAffineXf.h:19
M A
Definition MRMesh/MRAffineXf.h:18
Definition MRMesh/MRMatrix3.h:19
Vector3< T > x
rows, identity matrix by default
Definition MRMesh/MRMatrix3.h:24
Vector3< T > y
Definition MRMesh/MRMatrix3.h:25
Vector3< T > z
Definition MRMesh/MRMatrix3.h:26
Definition MRMatrix4.h:20
constexpr Matrix4(const AffineXf3< TT > &xf)
Definition MRMatrix4.h:46
friend constexpr bool operator!=(const Matrix4< T > &a, const Matrix4< T > &b)
Definition MRMatrix4.h:113
Vector4< T > z
Definition MRMatrix4.h:27
constexpr T * data()
Definition MRMatrix4.h:91
T ValueType
Definition MRMatrix4.h:21
constexpr T trace() const noexcept
computes trace of the matrix
Definition MRMatrix4.h:67
void setTranslation(const Vector3< T > &t) noexcept
friend constexpr auto operator*(T a, const Matrix4< T > &b) -> Matrix4< decltype(std::declval< T >() *std::declval< T >())>
Definition MRMatrix4.h:119
Matrix3< T > submatrix3(int i, int j) const noexcept
computes submatrix of the matrix with excluded i-th row and j-th column
static constexpr Matrix4 zero() noexcept
Definition MRMatrix4.h:50
constexpr T normSq() const noexcept
compute sum of squared matrix elements
Definition MRMatrix4.h:69
constexpr const T & operator()(int row, int col) const noexcept
element access
Definition MRMatrix4.h:56
Vector4< T > y
Definition MRMatrix4.h:26
static constexpr Matrix4 scale(T s) noexcept
returns a matrix that scales uniformly
Definition MRMatrix4.h:53
void setRotation(const Matrix3< T > &rot) noexcept
constexpr Matrix3< T > getRotation() const noexcept
constexpr Matrix4(const Matrix3< T > &r, const Vector3< T > &t)
construct from rotation matrix and translation vector
Definition MRMatrix4.h:35
friend constexpr auto operator/(Matrix4< T > b, T a) -> Matrix4< decltype(std::declval< T >()/std::declval< T >())>
Definition MRMatrix4.h:121
T det() const noexcept
computes determinant of the matrix
Vector4< T > w
Definition MRMatrix4.h:28
constexpr Matrix4< T > transposed() const noexcept
computes transposed matrix
constexpr auto norm() const noexcept
Definition MRMatrix4.h:70
auto dot(const Matrix4< T > &a, const Matrix4< T > &b) -> decltype(dot(a.x, b.x))
double-dot product: x = a : b
Definition MRMatrix4.h:162
friend constexpr auto operator+(const Matrix4< T > &a, const Matrix4< T > &b) -> Matrix4< decltype(std::declval< T >()+std::declval< T >())>
Definition MRMatrix4.h:117
constexpr Matrix4(const Matrix4< U > &m)
Definition MRMatrix4.h:49
constexpr const Vector4< T > & operator[](int row) const noexcept
row access
Definition MRMatrix4.h:60
friend constexpr Matrix4< T > & operator+=(Matrix4< T > &a, const Matrix4< T > &b)
Definition MRMatrix4.h:129
constexpr Vector4< T > col(int i) const noexcept
column access
Definition MRMatrix4.h:64
friend constexpr Matrix4< T > & operator*=(Matrix4< T > &a, T b)
Definition MRMatrix4.h:131
friend constexpr Matrix4< T > & operator-=(Matrix4< T > &a, const Matrix4< T > &b)
Definition MRMatrix4.h:130
friend constexpr auto operator-(const Matrix4< T > &a, const Matrix4< T > &b) -> Matrix4< decltype(std::declval< T >() - std::declval< T >())>
Definition MRMatrix4.h:118
constexpr Matrix4() noexcept=default
constexpr const T * data() const
Definition MRMatrix4.h:92
friend constexpr bool operator==(const Matrix4< T > &a, const Matrix4< T > &b)
Definition MRMatrix4.h:112
static constexpr Matrix4 identity() noexcept
Definition MRMatrix4.h:51
constexpr Matrix4< T > inverse() const noexcept
computes inverse matrix
friend constexpr Matrix4< T > & operator/=(Matrix4< T > &a, T b)
Definition MRMatrix4.h:132
constexpr Vector3< T > getTranslation() const noexcept
Vector4< T > x
rows, identity matrix by default
Definition MRMatrix4.h:25
Matrix4< T > outer(const Vector4< T > &a, const Vector4< T > &b)
x = a * b^T
Definition MRMatrix4.h:175
Definition MRMesh/MRVector3.h:26
T x
Definition MRMesh/MRVector3.h:32
T y
Definition MRMesh/MRVector3.h:32
T z
Definition MRMesh/MRVector3.h:32
Definition MRVector4.h:20
T y
Definition MRVector4.h:26
T z
Definition MRVector4.h:26
T x
Definition MRVector4.h:26
T w
Definition MRVector4.h:26