GteMatrix4x4.h
Go to the documentation of this file.
1 // David Eberly, Geometric Tools, Redmond WA 98052
2 // Copyright (c) 1998-2017
3 // Distributed under the Boost Software License, Version 1.0.
4 // http://www.boost.org/LICENSE_1_0.txt
5 // http://www.geometrictools.com/License/Boost/LICENSE_1_0.txt
6 // File Version: 3.0.0 (2016/06/19)
7 
8 #pragma once
9 
10 #include <Mathematics/GteMatrix.h>
11 #include <Mathematics/GteVector4.h>
12 
13 namespace gte
14 {
15 
16 // Template alias for convenience.
17 template <typename Real>
18 using Matrix4x4 = Matrix<4, 4, Real>;
19 
20 // Geometric operations.
21 template <typename Real>
22 Matrix4x4<Real> Inverse(Matrix4x4<Real> const& M,
23  bool* reportInvertibility = nullptr);
24 
25 template <typename Real>
26 Matrix4x4<Real> Adjoint(Matrix4x4<Real> const& M);
27 
28 template <typename Real>
29 Real Determinant(Matrix4x4<Real> const& M);
30 
31 template <typename Real>
32 Real Trace(Matrix4x4<Real> const& M);
33 
34 // Special matrices. In the comments, the matrices are shown using the
35 // GTE_USE_MAT_VEC multiplication convention.
36 
37 // The projection plane is Dot(N,X-P) = 0 where N is a 3-by-1 unit-length
38 // normal vector and P is a 3-by-1 point on the plane. The projection is
39 // oblique to the plane, in the direction of the 3-by-1 vector D. Necessarily
40 // Dot(N,D) is not zero for this projection to make sense. Given a 3-by-1
41 // point U, compute the intersection of the line U+t*D with the plane to
42 // obtain t = -Dot(N,U-P)/Dot(N,D); then
43 //
44 // projection(U) = P + [I - D*N^T/Dot(N,D)]*(U-P)
45 //
46 // A 4-by-4 homogeneous transformation representing the projection is
47 //
48 // +- -+
49 // M = | D*N^T - Dot(N,D)*I -Dot(N,P)D |
50 // | 0^T -Dot(N,D) |
51 // +- -+
52 //
53 // where M applies to [U^T 1]^T by M*[U^T 1]^T. The matrix is chosen so
54 // that M[3][3] > 0 whenever Dot(N,D) < 0; the projection is onto the
55 // "positive side" of the plane.
56 template <typename Real>
57 Matrix4x4<Real> MakeObliqueProjection(Vector4<Real> const& origin,
58  Vector4<Real> const& normal, Vector4<Real> const& direction);
59 
60 // The perspective projection of a point onto a plane is
61 //
62 // +- -+
63 // M = | Dot(N,E-P)*I - E*N^T -(Dot(N,E-P)*I - E*N^T)*E |
64 // | -N^t Dot(N,E) |
65 // +- -+
66 //
67 // where E is the eye point, P is a point on the plane, and N is a
68 // unit-length plane normal.
69 template <typename Real>
70 Matrix4x4<Real> MakePerspectiveProjection(Vector4<Real> const& origin,
71  Vector4<Real> const& normal, Vector4<Real> const& eye);
72 
73 // The reflection of a point through a plane is
74 // +- -+
75 // M = | I-2*N*N^T 2*Dot(N,P)*N |
76 // | 0^T 1 |
77 // +- -+
78 //
79 // where P is a point on the plane and N is a unit-length plane normal.
80 
81 template <typename Real>
82 Matrix4x4<Real> MakeReflection(Vector4<Real> const& origin,
83  Vector4<Real> const& normal);
84 
85 
86 template <typename Real>
87 Matrix4x4<Real> Inverse(Matrix4x4<Real> const& M, bool* reportInvertibility)
88 {
89  Matrix4x4<Real> inverse;
90  bool invertible;
91  Real a0 = M(0, 0) * M(1, 1) - M(0, 1) * M(1, 0);
92  Real a1 = M(0, 0) * M(1, 2) - M(0, 2) * M(1, 0);
93  Real a2 = M(0, 0) * M(1, 3) - M(0, 3) * M(1, 0);
94  Real a3 = M(0, 1) * M(1, 2) - M(0, 2) * M(1, 1);
95  Real a4 = M(0, 1) * M(1, 3) - M(0, 3) * M(1, 1);
96  Real a5 = M(0, 2) * M(1, 3) - M(0, 3) * M(1, 2);
97  Real b0 = M(2, 0) * M(3, 1) - M(2, 1) * M(3, 0);
98  Real b1 = M(2, 0) * M(3, 2) - M(2, 2) * M(3, 0);
99  Real b2 = M(2, 0) * M(3, 3) - M(2, 3) * M(3, 0);
100  Real b3 = M(2, 1) * M(3, 2) - M(2, 2) * M(3, 1);
101  Real b4 = M(2, 1) * M(3, 3) - M(2, 3) * M(3, 1);
102  Real b5 = M(2, 2) * M(3, 3) - M(2, 3) * M(3, 2);
103  Real det = a0 * b5 - a1 * b4 + a2 * b3 + a3 * b2 - a4 * b1 + a5 * b0;
104  if (det != (Real)0)
105  {
106  Real invDet = ((Real)1) / det;
107  inverse = Matrix4x4<Real>
108  {
109  (+M(1, 1) * b5 - M(1, 2) * b4 + M(1, 3) * b3) * invDet,
110  (-M(0, 1) * b5 + M(0, 2) * b4 - M(0, 3) * b3) * invDet,
111  (+M(3, 1) * a5 - M(3, 2) * a4 + M(3, 3) * a3) * invDet,
112  (-M(2, 1) * a5 + M(2, 2) * a4 - M(2, 3) * a3) * invDet,
113  (-M(1, 0) * b5 + M(1, 2) * b2 - M(1, 3) * b1) * invDet,
114  (+M(0, 0) * b5 - M(0, 2) * b2 + M(0, 3) * b1) * invDet,
115  (-M(3, 0) * a5 + M(3, 2) * a2 - M(3, 3) * a1) * invDet,
116  (+M(2, 0) * a5 - M(2, 2) * a2 + M(2, 3) * a1) * invDet,
117  (+M(1, 0) * b4 - M(1, 1) * b2 + M(1, 3) * b0) * invDet,
118  (-M(0, 0) * b4 + M(0, 1) * b2 - M(0, 3) * b0) * invDet,
119  (+M(3, 0) * a4 - M(3, 1) * a2 + M(3, 3) * a0) * invDet,
120  (-M(2, 0) * a4 + M(2, 1) * a2 - M(2, 3) * a0) * invDet,
121  (-M(1, 0) * b3 + M(1, 1) * b1 - M(1, 2) * b0) * invDet,
122  (+M(0, 0) * b3 - M(0, 1) * b1 + M(0, 2) * b0) * invDet,
123  (-M(3, 0) * a3 + M(3, 1) * a1 - M(3, 2) * a0) * invDet,
124  (+M(2, 0) * a3 - M(2, 1) * a1 + M(2, 2) * a0) * invDet
125  };
126  invertible = true;
127  }
128  else
129  {
130  inverse.MakeZero();
131  invertible = false;
132  }
133 
134  if (reportInvertibility)
135  {
136  *reportInvertibility = invertible;
137  }
138  return inverse;
139 }
140 
141 template <typename Real>
142 Matrix4x4<Real> Adjoint(Matrix4x4<Real> const& M)
143 {
144  Real a0 = M(0, 0) * M(1, 1) - M(0, 1) * M(1, 0);
145  Real a1 = M(0, 0) * M(1, 2) - M(0, 2) * M(1, 0);
146  Real a2 = M(0, 0) * M(1, 3) - M(0, 3) * M(1, 0);
147  Real a3 = M(0, 1) * M(1, 2) - M(0, 2) * M(1, 1);
148  Real a4 = M(0, 1) * M(1, 3) - M(0, 3) * M(1, 1);
149  Real a5 = M(0, 2) * M(1, 3) - M(0, 3) * M(1, 2);
150  Real b0 = M(2, 0) * M(3, 1) - M(2, 1) * M(3, 0);
151  Real b1 = M(2, 0) * M(3, 2) - M(2, 2) * M(3, 0);
152  Real b2 = M(2, 0) * M(3, 3) - M(2, 3) * M(3, 0);
153  Real b3 = M(2, 1) * M(3, 2) - M(2, 2) * M(3, 1);
154  Real b4 = M(2, 1) * M(3, 3) - M(2, 3) * M(3, 1);
155  Real b5 = M(2, 2) * M(3, 3) - M(2, 3) * M(3, 2);
156 
157  return Matrix4x4<Real>
158  {
159  +M(1, 1) * b5 - M(1, 2) * b4 + M(1, 3) * b3,
160  -M(0, 1) * b5 + M(0, 2) * b4 - M(0, 3) * b3,
161  +M(3, 1) * a5 - M(3, 2) * a4 + M(3, 3) * a3,
162  -M(2, 1) * a5 + M(2, 2) * a4 - M(2, 3) * a3,
163  -M(1, 0) * b5 + M(1, 2) * b2 - M(1, 3) * b1,
164  +M(0, 0) * b5 - M(0, 2) * b2 + M(0, 3) * b1,
165  -M(3, 0) * a5 + M(3, 2) * a2 - M(3, 3) * a1,
166  +M(2, 0) * a5 - M(2, 2) * a2 + M(2, 3) * a1,
167  +M(1, 0) * b4 - M(1, 1) * b2 + M(1, 3) * b0,
168  -M(0, 0) * b4 + M(0, 1) * b2 - M(0, 3) * b0,
169  +M(3, 0) * a4 - M(3, 1) * a2 + M(3, 3) * a0,
170  -M(2, 0) * a4 + M(2, 1) * a2 - M(2, 3) * a0,
171  -M(1, 0) * b3 + M(1, 1) * b1 - M(1, 2) * b0,
172  +M(0, 0) * b3 - M(0, 1) * b1 + M(0, 2) * b0,
173  -M(3, 0) * a3 + M(3, 1) * a1 - M(3, 2) * a0,
174  +M(2, 0) * a3 - M(2, 1) * a1 + M(2, 2) * a0
175  };
176 }
177 
178 template <typename Real>
179 Real Determinant(Matrix4x4<Real> const& M)
180 {
181  Real a0 = M(0, 0) * M(1, 1) - M(0, 1) * M(1, 0);
182  Real a1 = M(0, 0) * M(1, 2) - M(0, 2) * M(1, 0);
183  Real a2 = M(0, 0) * M(1, 3) - M(0, 3) * M(1, 0);
184  Real a3 = M(0, 1) * M(1, 2) - M(0, 2) * M(1, 1);
185  Real a4 = M(0, 1) * M(1, 3) - M(0, 3) * M(1, 1);
186  Real a5 = M(0, 2) * M(1, 3) - M(0, 3) * M(1, 2);
187  Real b0 = M(2, 0) * M(3, 1) - M(2, 1) * M(3, 0);
188  Real b1 = M(2, 0) * M(3, 2) - M(2, 2) * M(3, 0);
189  Real b2 = M(2, 0) * M(3, 3) - M(2, 3) * M(3, 0);
190  Real b3 = M(2, 1) * M(3, 2) - M(2, 2) * M(3, 1);
191  Real b4 = M(2, 1) * M(3, 3) - M(2, 3) * M(3, 1);
192  Real b5 = M(2, 2) * M(3, 3) - M(2, 3) * M(3, 2);
193  Real det = a0 * b5 - a1 * b4 + a2 * b3 + a3 * b2 - a4 * b1 + a5 * b0;
194  return det;
195 }
196 
197 template <typename Real>
198 Real Trace(Matrix4x4<Real> const& M)
199 {
200  Real trace = M(0, 0) + M(1, 1) + M(2, 2) + M(3, 3);
201  return trace;
202 }
203 
204 template <typename Real>
205 Matrix4x4<Real> MakeObliqueProjection(Vector4<Real> const& origin,
206  Vector4<Real> const& normal, Vector4<Real> const& direction)
207 {
208  Matrix4x4<Real> M;
209 
210  Real const zero = (Real)0;
211  Real dotND = Dot(normal, direction);
212  Real dotNO = Dot(origin, normal);
213 
214 #if defined(GTE_USE_MAT_VEC)
215  M(0, 0) = direction[0] * normal[0] - dotND;
216  M(0, 1) = direction[0] * normal[1];
217  M(0, 2) = direction[0] * normal[2];
218  M(0, 3) = -dotNO * direction[0];
219  M(1, 0) = direction[1] * normal[0];
220  M(1, 1) = direction[1] * normal[1] - dotND;
221  M(1, 2) = direction[1] * normal[2];
222  M(1, 3) = -dotNO * direction[1];
223  M(2, 0) = direction[2] * normal[0];
224  M(2, 1) = direction[2] * normal[1];
225  M(2, 2) = direction[2] * normal[2] - dotND;
226  M(2, 3) = -dotNO * direction[2];
227  M(3, 0) = zero;
228  M(3, 1) = zero;
229  M(3, 2) = zero;
230  M(3, 3) = -dotND;
231 #else
232  M(0, 0) = direction[0] * normal[0] - dotND;
233  M(1, 0) = direction[0] * normal[1];
234  M(2, 0) = direction[0] * normal[2];
235  M(3, 0) = -dotNO * direction[0];
236  M(0, 1) = direction[1] * normal[0];
237  M(1, 1) = direction[1] * normal[1] - dotND;
238  M(2, 1) = direction[1] * normal[2];
239  M(3, 1) = -dotNO * direction[1];
240  M(0, 2) = direction[2] * normal[0];
241  M(1, 2) = direction[2] * normal[1];
242  M(2, 2) = direction[2] * normal[2] - dotND;
243  M(3, 2) = -dotNO * direction[2];
244  M(0, 2) = zero;
245  M(1, 3) = zero;
246  M(2, 3) = zero;
247  M(3, 3) = -dotND;
248 #endif
249 
250  return M;
251 }
252 
253 template <typename Real>
254 Matrix4x4<Real> MakePerspectiveProjection(Vector4<Real> const& origin,
255  Vector4<Real> const& normal, Vector4<Real> const& eye)
256 {
257  Matrix4x4<Real> M;
258 
259  Real dotND = Dot(normal, eye - origin);
260 
261 #if defined(GTE_USE_MAT_VEC)
262  M(0, 0) = dotND - eye[0] * normal[0];
263  M(0, 1) = -eye[0] * normal[1];
264  M(0, 2) = -eye[0] * normal[2];
265  M(0, 3) = -(M(0, 0) * eye[0] + M(0, 1) * eye[1] + M(0, 2) * eye[2]);
266  M(1, 0) = -eye[1] * normal[0];
267  M(1, 1) = dotND - eye[1] * normal[1];
268  M(1, 2) = -eye[1] * normal[2];
269  M(1, 3) = -(M(1, 0) * eye[0] + M(1, 1) * eye[1] + M(1, 2) * eye[2]);
270  M(2, 0) = -eye[2] * normal[0];
271  M(2, 1) = -eye[2] * normal[1];
272  M(2, 2) = dotND - eye[2] * normal[2];
273  M(2, 3) = -(M(2, 0) * eye[0] + M(2, 1) * eye[1] + M(2, 2) * eye[2]);
274  M(3, 0) = -normal[0];
275  M(3, 1) = -normal[1];
276  M(3, 2) = -normal[2];
277  M(3, 3) = Dot(eye, normal);
278 #else
279  M(0, 0) = dotND - eye[0] * normal[0];
280  M(1, 0) = -eye[0] * normal[1];
281  M(2, 0) = -eye[0] * normal[2];
282  M(3, 0) = -(M(0, 0) * eye[0] + M(0, 1) * eye[1] + M(0, 2) * eye[2]);
283  M(0, 1) = -eye[1] * normal[0];
284  M(1, 1) = dotND - eye[1] * normal[1];
285  M(2, 1) = -eye[1] * normal[2];
286  M(3, 1) = -(M(1, 0) * eye[0] + M(1, 1) * eye[1] + M(1, 2) * eye[2]);
287  M(0, 2) = -eye[2] * normal[0];
288  M(1, 2) = -eye[2] * normal[1];
289  M(2, 2) = dotND - eye[2] * normal[2];
290  M(3, 2) = -(M(2, 0) * eye[0] + M(2, 1) * eye[1] + M(2, 2) * eye[2]);
291  M(0, 3) = -normal[0];
292  M(1, 3) = -normal[1];
293  M(2, 3) = -normal[2];
294  M(3, 3) = Dot(eye, normal);
295 #endif
296 
297  return M;
298 }
299 
300 template <typename Real>
301 Matrix4x4<Real> MakeReflection(Vector4<Real> const& origin,
302  Vector4<Real> const& normal)
303 {
304  Matrix4x4<Real> M;
305 
306  Real const zero = (Real)0, one = (Real)1, two = (Real)2;
307  Real twoDotNO = two * Dot(origin, normal);
308 
309 #if defined(GTE_USE_MAT_VEC)
310  M(0, 0) = one - two * normal[0] * normal[0];
311  M(0, 1) = -two * normal[0] * normal[1];
312  M(0, 2) = -two * normal[0] * normal[2];
313  M(0, 3) = twoDotNO * normal[0];
314  M(1, 0) = M(0, 1);
315  M(1, 1) = one - two * normal[1] * normal[1];
316  M(1, 2) = -two * normal[1] * normal[2];
317  M(1, 3) = twoDotNO * normal[1];
318  M(2, 0) = M(0, 2);
319  M(2, 1) = M(1, 2);
320  M(2, 2) = one - two * normal[2] * normal[2];
321  M(2, 3) = twoDotNO * normal[2];
322  M(3, 0) = zero;
323  M(3, 1) = zero;
324  M(3, 2) = zero;
325  M(3, 3) = one;
326 #else
327  M(0, 0) = one - two * normal[0] * normal[0];
328  M(1, 0) = -two * normal[0] * normal[1];
329  M(2, 0) = -two * normal[0] * normal[2];
330  M(3, 0) = twoDotNO * normal[0];
331  M(0, 1) = M(1, 0);
332  M(1, 1) = one - two * normal[1] * normal[1];
333  M(2, 1) = -two * normal[1] * normal[2];
334  M(3, 1) = twoDotNO * normal[1];
335  M(0, 2) = M(2, 0);
336  M(1, 2) = M(2, 1);
337  M(2, 2) = one - two * normal[2] * normal[2];
338  M(3, 2) = twoDotNO * normal[2];
339  M(0, 3) = zero;
340  M(1, 3) = zero;
341  M(2, 3) = zero;
342  M(3, 3) = one;
343 #endif
344 
345  return M;
346 }
347 
348 
349 }
gte::Adjoint
Matrix2x2< Real > Adjoint(Matrix2x2< Real > const &M)
Definition: GteMatrix2x2.h:108
gte::Trace
Real Trace(Matrix2x2< Real > const &M)
Definition: GteMatrix2x2.h:125
gte::Dot
DualQuaternion< Real > Dot(DualQuaternion< Real > const &d0, DualQuaternion< Real > const &d1)
gte::MakeReflection
Matrix4x4< Real > MakeReflection(Vector4< Real > const &origin, Vector4< Real > const &normal)
Definition: GteMatrix4x4.h:301
gte::MakeObliqueProjection
Matrix4x4< Real > MakeObliqueProjection(Vector4< Real > const &origin, Vector4< Real > const &normal, Vector4< Real > const &direction)
Definition: GteMatrix4x4.h:205
gte::Matrix::MakeZero
void MakeZero()
Definition: GteMatrix.h:442
gte::Inverse
Quaternion< Real > Inverse(Quaternion< Real > const &d)
Definition: GteQuaternion.h:227
gte::Determinant
Real Determinant(GMatrix< Real > const &M)
Definition: GteGMatrix.h:618
gte::MakePerspectiveProjection
Matrix4x4< Real > MakePerspectiveProjection(Vector4< Real > const &origin, Vector4< Real > const &normal, Vector4< Real > const &eye)
Definition: GteMatrix4x4.h:254

geometric_tools_engine
Author(s): Yijiang Huang
autogenerated on Thu Jul 18 2019 04:00:01