QGIS API Documentation 4.1.0-Master (64029b4150b)
Loading...
Searching...
No Matches
qgsmatrix4x4.cpp
Go to the documentation of this file.
1/***************************************************************************
2 qgsmatrix4x4.cpp
3 --------------------------------------
4 Date : July 2023
5 Copyright : (C) 2023 by Martin Dobias
6 Email : wonder dot sk at gmail dot com
7 ***************************************************************************
8 * *
9 * This program is free software; you can redistribute it and/or modify *
10 * it under the terms of the GNU General Public License as published by *
11 * the Free Software Foundation; either version 2 of the License, or *
12 * (at your option) any later version. *
13 * *
14 ***************************************************************************/
15
16#include "qgsmatrix4x4.h"
17
18// the implementation is partially based on Qt's QMatrix4x4 (simplified)
19
20
21// clang-format off
23 double m11, double m12, double m13, double m14,
24 double m21, double m22, double m23, double m24,
25 double m31, double m32, double m33, double m34,
26 double m41, double m42, double m43, double m44
27)
28// clang-format on
29{
30 m[0][0] = m11;
31 m[0][1] = m21;
32 m[0][2] = m31;
33 m[0][3] = m41;
34 m[1][0] = m12;
35 m[1][1] = m22;
36 m[1][2] = m32;
37 m[1][3] = m42;
38 m[2][0] = m13;
39 m[2][1] = m23;
40 m[2][2] = m33;
41 m[2][3] = m43;
42 m[3][0] = m14;
43 m[3][1] = m24;
44 m[3][2] = m34;
45 m[3][3] = m44;
46}
47
48bool QgsMatrix4x4::fuzzyEqual( const QgsMatrix4x4 &other, double epsilon ) const
49{
50 const double *data = *m;
51 const double *otherData = *( other.m );
52 for ( int i = 0; i < 16; ++i, data++, otherData++ )
53 {
54 if ( !qgsDoubleNear( *data, *otherData, epsilon ) )
55 {
56 return false;
57 }
58 }
59
60 return true;
61}
62
64{
65 m[3][0] += m[0][0] * vector.x() + m[1][0] * vector.y() + m[2][0] * vector.z();
66 m[3][1] += m[0][1] * vector.x() + m[1][1] * vector.y() + m[2][1] * vector.z();
67 m[3][2] += m[0][2] * vector.x() + m[1][2] * vector.y() + m[2][2] * vector.z();
68 m[3][3] += m[0][3] * vector.x() + m[1][3] * vector.y() + m[2][3] * vector.z();
69}
70
71QList< double > QgsMatrix4x4::dataList() const
72{
73 QList< double > res;
74 res.reserve( 16 );
75 for ( int i = 0; i < 16; ++i )
76 {
77 res.append( m[i / 4][i % 4] );
78 }
79 return res;
80}
81
82QgsVector3D operator*( const QgsMatrix4x4 &matrix, const QgsVector3D &vector )
83{
84 double x, y, z, w;
85
86 x = vector.x() * matrix.m[0][0] + vector.y() * matrix.m[1][0] + vector.z() * matrix.m[2][0] + matrix.m[3][0];
87 y = vector.x() * matrix.m[0][1] + vector.y() * matrix.m[1][1] + vector.z() * matrix.m[2][1] + matrix.m[3][1];
88 z = vector.x() * matrix.m[0][2] + vector.y() * matrix.m[1][2] + vector.z() * matrix.m[2][2] + matrix.m[3][2];
89 w = vector.x() * matrix.m[0][3] + vector.y() * matrix.m[1][3] + vector.z() * matrix.m[2][3] + matrix.m[3][3];
90 if ( w == 1.0 )
91 return QgsVector3D( x, y, z );
92 else
93 return QgsVector3D( x / w, y / w, z / w );
94}
95
96// Simplified from Qt's QMatrix4x4::mapVector implementation.
97// Copyright (C) The Qt Company Ltd.
99{
100 return QgsVector3D(
101 vector.x() * m[0][0] + vector.y() * m[1][0] + vector.z() * m[2][0],
102
103 vector.x() * m[0][1] + vector.y() * m[1][1] + vector.z() * m[2][1],
104
105 vector.x() * m[0][2] + vector.y() * m[1][2] + vector.z() * m[2][2]
106 );
107}
108
110{
111 if ( m[0][0] != 1.0 || m[0][1] != 0.0 || m[0][2] != 0.0 )
112 return false;
113 if ( m[0][3] != 0.0 || m[1][0] != 0.0 || m[1][1] != 1.0 )
114 return false;
115 if ( m[1][2] != 0.0 || m[1][3] != 0.0 || m[2][0] != 0.0 )
116 return false;
117 if ( m[2][1] != 0.0 || m[2][2] != 1.0 || m[2][3] != 0.0 )
118 return false;
119 if ( m[3][0] != 0.0 || m[3][1] != 0.0 || m[3][2] != 0.0 )
120 return false;
121 return ( m[3][3] == 1.0 );
122}
123
125{
126 m[0][0] = 1.0;
127 m[0][1] = 0.0;
128 m[0][2] = 0.0;
129 m[0][3] = 0.0;
130 m[1][0] = 0.0;
131 m[1][1] = 1.0;
132 m[1][2] = 0.0;
133 m[1][3] = 0.0;
134 m[2][0] = 0.0;
135 m[2][1] = 0.0;
136 m[2][2] = 1.0;
137 m[2][3] = 0.0;
138 m[3][0] = 0.0;
139 m[3][1] = 0.0;
140 m[3][2] = 0.0;
141 m[3][3] = 1.0;
142}
143
144
146{
147 QgsMatrix4x4 m( 1 );
148 m.m[0][0] = m1.m[0][0] * m2.m[0][0] + m1.m[1][0] * m2.m[0][1] + m1.m[2][0] * m2.m[0][2] + m1.m[3][0] * m2.m[0][3];
149 m.m[0][1] = m1.m[0][1] * m2.m[0][0] + m1.m[1][1] * m2.m[0][1] + m1.m[2][1] * m2.m[0][2] + m1.m[3][1] * m2.m[0][3];
150 m.m[0][2] = m1.m[0][2] * m2.m[0][0] + m1.m[1][2] * m2.m[0][1] + m1.m[2][2] * m2.m[0][2] + m1.m[3][2] * m2.m[0][3];
151 m.m[0][3] = m1.m[0][3] * m2.m[0][0] + m1.m[1][3] * m2.m[0][1] + m1.m[2][3] * m2.m[0][2] + m1.m[3][3] * m2.m[0][3];
152
153 m.m[1][0] = m1.m[0][0] * m2.m[1][0] + m1.m[1][0] * m2.m[1][1] + m1.m[2][0] * m2.m[1][2] + m1.m[3][0] * m2.m[1][3];
154 m.m[1][1] = m1.m[0][1] * m2.m[1][0] + m1.m[1][1] * m2.m[1][1] + m1.m[2][1] * m2.m[1][2] + m1.m[3][1] * m2.m[1][3];
155 m.m[1][2] = m1.m[0][2] * m2.m[1][0] + m1.m[1][2] * m2.m[1][1] + m1.m[2][2] * m2.m[1][2] + m1.m[3][2] * m2.m[1][3];
156 m.m[1][3] = m1.m[0][3] * m2.m[1][0] + m1.m[1][3] * m2.m[1][1] + m1.m[2][3] * m2.m[1][2] + m1.m[3][3] * m2.m[1][3];
157
158 m.m[2][0] = m1.m[0][0] * m2.m[2][0] + m1.m[1][0] * m2.m[2][1] + m1.m[2][0] * m2.m[2][2] + m1.m[3][0] * m2.m[2][3];
159 m.m[2][1] = m1.m[0][1] * m2.m[2][0] + m1.m[1][1] * m2.m[2][1] + m1.m[2][1] * m2.m[2][2] + m1.m[3][1] * m2.m[2][3];
160 m.m[2][2] = m1.m[0][2] * m2.m[2][0] + m1.m[1][2] * m2.m[2][1] + m1.m[2][2] * m2.m[2][2] + m1.m[3][2] * m2.m[2][3];
161 m.m[2][3] = m1.m[0][3] * m2.m[2][0] + m1.m[1][3] * m2.m[2][1] + m1.m[2][3] * m2.m[2][2] + m1.m[3][3] * m2.m[2][3];
162
163 m.m[3][0] = m1.m[0][0] * m2.m[3][0] + m1.m[1][0] * m2.m[3][1] + m1.m[2][0] * m2.m[3][2] + m1.m[3][0] * m2.m[3][3];
164 m.m[3][1] = m1.m[0][1] * m2.m[3][0] + m1.m[1][1] * m2.m[3][1] + m1.m[2][1] * m2.m[3][2] + m1.m[3][1] * m2.m[3][3];
165 m.m[3][2] = m1.m[0][2] * m2.m[3][0] + m1.m[1][2] * m2.m[3][1] + m1.m[2][2] * m2.m[3][2] + m1.m[3][2] * m2.m[3][3];
166 m.m[3][3] = m1.m[0][3] * m2.m[3][0] + m1.m[1][3] * m2.m[3][1] + m1.m[2][3] * m2.m[3][2] + m1.m[3][3] * m2.m[3][3];
167 return m;
168}
169
170void QgsMatrix4x4::scale( const QgsVector3D &vector )
171{
172 m[0][0] *= vector.x();
173 m[0][1] *= vector.x();
174 m[0][2] *= vector.x();
175 m[0][3] *= vector.x();
176
177 m[1][0] *= vector.y();
178 m[1][1] *= vector.y();
179 m[1][2] *= vector.y();
180 m[1][3] *= vector.y();
181
182 m[2][0] *= vector.z();
183 m[2][1] *= vector.z();
184 m[2][2] *= vector.z();
185 m[2][3] *= vector.z();
186}
187
188void QgsMatrix4x4::rotate( double angle, double x, double y, double z )
189{
190 const double length = std::sqrt( x * x + y * y + z * z );
191 if ( qgsDoubleNear( length, 0.0 ) )
192 {
193 return;
194 }
195
196 const double nx = x / length;
197 const double ny = y / length;
198 const double nz = z / length;
199
200 const double angleRad = angle * M_PI / 180.0;
201 const double angleCos = std::cos( angleRad );
202 const double angleSin = std::sin( angleRad );
203 const double angleCosInv = 1.0 - angleCos;
204
205 const double nxny = nx * ny;
206 const double nxnz = nx * nz;
207 const double nynz = ny * nz;
208
209 // Rotation matrix
210 QgsMatrix4x4 rot;
211 rot.m[0][0] = angleCosInv * nx * nx + angleCos;
212 rot.m[0][1] = angleCosInv * nxny + nz * angleSin;
213 rot.m[0][2] = angleCosInv * nxnz - ny * angleSin;
214 rot.m[0][3] = 0.0;
215
216 rot.m[1][0] = angleCosInv * nxny - nz * angleSin;
217 rot.m[1][1] = angleCosInv * ny * ny + angleCos;
218 rot.m[1][2] = angleCosInv * nynz + nx * angleSin;
219 rot.m[1][3] = 0.0;
220
221 rot.m[2][0] = angleCosInv * nxnz + ny * angleSin;
222 rot.m[2][1] = angleCosInv * nynz - nx * angleSin;
223 rot.m[2][2] = angleCosInv * nz * nz + angleCos;
224 rot.m[2][3] = 0.0;
225
226 rot.m[3][0] = 0.0;
227 rot.m[3][1] = 0.0;
228 rot.m[3][2] = 0.0;
229 rot.m[3][3] = 1.0;
230
231 *this = *this * rot;
232}
233
234void QgsMatrix4x4::rotate( double angle, const QgsVector3D &vector )
235{
236 rotate( angle, vector.x(), vector.y(), vector.z() );
237}
238
239static inline double matrixDet2( const double m[4][4], int col0, int col1, int row0, int row1 )
240{
241 return m[col0][row0] * m[col1][row1] - m[col0][row1] * m[col1][row0];
242}
243
244// The 4x4 matrix inverse algorithm is based on that described at:
245// http://www.j3d.org/matrix_faq/matrfaq_latest.html#Q24
246// Some optimization has been done to avoid making copies of 3x3
247// sub-matrices and to unroll the loops.
248// Calculate the determinant of a 3x3 sub-matrix.
249// | A B C |
250// M = | D E F | det(M) = A * (EI - HF) - B * (DI - GF) + C * (DH - GE)
251// | G H I |
252static inline double matrixDet3( const double m[4][4], int col0, int col1, int col2, int row0, int row1, int row2 )
253{
254 return m[col0][row0] * matrixDet2( m, col1, col2, row1, row2 ) - m[col1][row0] * matrixDet2( m, col0, col2, row1, row2 ) + m[col2][row0] * matrixDet2( m, col0, col1, row1, row2 );
255}
256
257// Calculate the determinant of a 4x4 matrix.
258static inline double matrixDet4( const double m[4][4] )
259{
260 double det;
261 det = m[0][0] * matrixDet3( m, 1, 2, 3, 1, 2, 3 );
262 det -= m[1][0] * matrixDet3( m, 0, 2, 3, 1, 2, 3 );
263 det += m[2][0] * matrixDet3( m, 0, 1, 3, 1, 2, 3 );
264 det -= m[3][0] * matrixDet3( m, 0, 1, 2, 1, 2, 3 );
265 return det;
266}
267
268// Simplified from Qt's QMatrix4x4::determinant implementation.
269// Copyright (C) The Qt Company Ltd.
271{
272 return matrixDet4( m );
273}
QgsVector3D mapVector(const QgsVector3D &vector) const
Maps vector by multiplying the top 3x3 portion of this matrix by vector.
void scale(const QgsVector3D &vector)
Multiplies this matrix by another that scales coordinates by the components of a vector.
QgsMatrix4x4()
Initializes identity matrix.
void setToIdentity()
Sets matrix to be identity matrix.
bool isIdentity() const
Returns whether this matrix is an identity matrix.
void translate(const QgsVector3D &vector)
Multiplies this matrix by another that translates coordinates by the components of a vector.
bool fuzzyEqual(const QgsMatrix4x4 &other, double epsilon=1e-8) const
Performs fuzzy comparison between this matrix and other using an epsilon.
void rotate(double angle, double x, double y, double z)
Multiples this matrix by another that rotates coordinates through angle degrees about the vector (x,...
double determinant() const
Returns the determinant of this matrix.
double * data()
Returns pointer to the matrix data (stored in column-major order).
QList< double > dataList() const
Returns matrix data (in column-major order).
A 3D vector (similar to QVector3D) with the difference that it uses double precision instead of singl...
Definition qgsvector3d.h:33
double y() const
Returns Y coordinate.
Definition qgsvector3d.h:60
double z() const
Returns Z coordinate.
Definition qgsvector3d.h:62
double x() const
Returns X coordinate.
Definition qgsvector3d.h:58
bool qgsDoubleNear(double a, double b, double epsilon=4 *std::numeric_limits< double >::epsilon())
Compare two doubles (but allow some difference).
Definition qgis.h:7209
QgsVector3D operator*(const QgsMatrix4x4 &matrix, const QgsVector3D &vector)
Matrix-vector multiplication (vector is converted to homogeneous coordinates [X,Y,...