Ocean
Loading...
Searching...
No Matches
SquareMatrix3.h
Go to the documentation of this file.
1/*
2 * Copyright (c) Meta Platforms, Inc. and affiliates.
3 *
4 * This source code is licensed under the MIT license found in the
5 * LICENSE file in the root directory of this source tree.
6 */
7
8#ifndef META_OCEAN_MATH_SQUARE_MATRIX_3_H
9#define META_OCEAN_MATH_SQUARE_MATRIX_3_H
10
11#include "ocean/math/Math.h"
12#include "ocean/math/Equation.h"
14#include "ocean/math/Vector2.h"
15#include "ocean/math/Vector3.h"
16
18
19#include <vector>
20
21namespace Ocean
22{
23
24// Forward declaration.
25template <typename T> class EulerT;
26// Forward declaration.
27template <typename T> class HomogenousMatrixT4;
28// Forward declaration.
29template <typename T> class QuaternionT;
30// Forward declaration.
31template <typename T> class RotationT;
32// Forward declaration.
33template <typename T> class SquareMatrixT4;
34
35// Forward declaration.
36template <typename T> class SquareMatrixT3;
37
38/**
39 * Definition of the SquareMatrix3 object, depending on the OCEAN_MATH_USE_SINGLE_PRECISION either with single or double precision float data type.
40 * @see SquareMatrixT3
41 * @ingroup math
42 */
44
45/**
46 * Instantiation of the SquareMatrixT3 template class using a double precision float data type.
47 * @see SquareMatrixT3
48 * @ingroup math
49 */
51
52/**
53 * Instantiation of the SquareMatrixT3 template class using a single precision float data type.
54 * @see SquareMatrixT3
55 * @ingroup math
56 */
58
59/**
60 * Definition of a typename alias for vectors with SquareMatrixT3 objects.
61 * @see SquareMatrixT3
62 * @ingroup math
63 */
64template <typename T>
65using SquareMatricesT3 = std::vector<SquareMatrixT3<T>>;
66
67/**
68 * Definition of a vector holding SquareMatrix3 objects.
69 * @see SquareMatrix3
70 * @ingroup math
71 */
72using SquareMatrices3 = std::vector<SquareMatrix3>;
73
74/**
75 * This class implements a 3x3 square matrix.
76 * The matrix can be applied as e.g., rotation matrix for 3D vectors or can represent a Homography and so on.<br>
77 * The values are stored in a column aligned order with indices:
78 * <pre>
79 * | 0 3 6 |
80 * | 1 4 7 |
81 * | 2 5 8 |
82 * </pre>
83 * @tparam T Data type of matrix elements
84 * @see SquareMatrix3, SquareMatrixF3, SquareMatrixD3, Rotation, Euler, Quaternion, ExponentialMap.
85 * @ingroup math
86 */
87template <typename T>
89{
90 template <typename U> friend class SquareMatrixT3;
91
92 public:
93
94 /**
95 * Definition of the used data type.
96 */
97 using Type = T;
98
99 public:
100
101 /**
102 * Creates a new SquareMatrixT3 object with undefined elements.
103 * Beware: This matrix is neither a zero nor an entity matrix!
104 */
106
107 /**
108 * Copy constructor.
109 * @param matrix The matrix to copy
110 */
111 SquareMatrixT3(const SquareMatrixT3<T>& matrix) = default;
112
113 /**
114 * Copy constructor for a matrix with difference element data type than T.
115 * @param matrix The matrix to copy
116 * @tparam U The element data type of the second matrix
117 */
118 template <typename U>
119 explicit inline SquareMatrixT3(const SquareMatrixT3<U>& matrix);
120
121 /**
122 * Creates a new SquareMatrixT3 object.
123 * @param setToIdentity Determines whether a entity matrix will be created, otherwise the matrix is initialized with zeros
124 */
125 explicit SquareMatrixT3(const bool setToIdentity);
126
127 /**
128 * Creates a new SquareMatrixT3 rotation matrix by a given Euler rotation.
129 * @param euler Euler rotation to create a rotation matrix from, must be valid
130 */
131 explicit SquareMatrixT3(const EulerT<T>& euler);
132
133 /**
134 * Creates a new 3x3 matrix object by a given angle-axis rotation.
135 * @param rotation The angle-axis rotation to create a matrix from, must be valid
136 */
137 explicit SquareMatrixT3(const RotationT<T>& rotation);
138
139 /**
140 * Creates a new 3x3 matrix object by a given quaternion rotation.
141 * @param quaternion The quaternion rotation to create a matrix from, must be valid
142 */
143 explicit SquareMatrixT3(const QuaternionT<T>& quaternion);
144
145 /**
146 * Creates a new SquareMatrixT3 object by three axes.
147 * @param xAxis First column of the matrix
148 * @param yAxis Middle column of the matrix
149 * @param zAxis Last column of the matrix
150 */
152
153 /**
154 * Creates a new SquareMatrixT3 object by a given diagonal vector.
155 * @param diagonal The diagonal vector for the new matrix
156 */
158
159 /**
160 * Creates a new SquareMatrixT3 object by nine elements of float type U.
161 * @param arrayValues The nine matrix elements defining the new matrix, must be valid
162 * @tparam U The floating point type of the given elements
163 */
164 template <typename U>
165 explicit SquareMatrixT3(const U* arrayValues);
166
167 /**
168 * Creates a new SquareMatrixT3 object by nine elements.
169 * @param arrayValues The nine matrix elements defining the new matrix, must be valid
170 */
171 explicit SquareMatrixT3(const T* arrayValues);
172
173 /**
174 * Creates a new SquareMatrixT3 object by nine elements.
175 * @param arrayValues The nine matrix elements defining the new matrix, must be valid
176 * @param valuesRowAligned True, if the given values are stored in a row aligned order; False, if the values are stored in a column aligned order (which is the default case for this matrix)
177 * @tparam U The floating point type of the given elements
178 */
179 template <typename U>
180 SquareMatrixT3(const U* arrayValues, const bool valuesRowAligned);
181
182 /**
183 * Creates a new SquareMatrixT3 object by nine elements.
184 * @param arrayValues The nine matrix elements defining the new matrix, must be valid
185 * @param valuesRowAligned True, if the given values are stored in a row aligned order; False, if the values are stored in a column aligned order (which is the default case for this matrix)
186 */
187 SquareMatrixT3(const T* arrayValues, const bool valuesRowAligned);
188
189 /**
190 * Creates a 3x3 rotation matrix by a given 4x4 homogeneous transformation.
191 * @param transformation The transformation to create a 3x3 rotation matrix from
192 */
193 explicit SquareMatrixT3(const HomogenousMatrixT4<T>& transformation);
194
195 /**
196 * Creates a 3x3 square matrix by a given 4x4 square transformation.
197 * @param transformation The transformation to create a 3x3 square matrix from
198 */
199 explicit SquareMatrixT3(const SquareMatrixT4<T>& transformation);
200
201 /**
202 * Creates a 3x3 rotation matrix by 9 given matrix elements.
203 * @param m00 Element of the first row and first column
204 * @param m10 Element of the second row and first column
205 * @param m20 Element of the third row and first column
206 * @param m01 Element of the first row and second column
207 * @param m11 Element of the second row and second column
208 * @param m21 Element of the third row and second column
209 * @param m02 Element of the first row and third column
210 * @param m12 Element of the second row and third column
211 * @param m22 Element of the third row and third column
212 */
213 explicit SquareMatrixT3(const T& m00, const T& m10, const T& m20, const T& m01, const T& m11, const T& m21, const T& m02, const T& m12, const T& m22);
214
215 /**
216 * Returns the transposed of this matrix.
217 * @return Transposed matrix
218 */
220
221 /**
222 * Transposes the matrix.
223 */
224 void transpose();
225
226 /**
227 * Returns the inverted matrix of this matrix.
228 * This matrix must not be singular.<br>
229 * Beware: This function does not throw an exception if the matrix cannot be inverted.<br>
230 * Thus, ensure that the matrix is invertible before calling this function.<br>
231 * Even better: avoid the usage of this function and call invert() instead.<br>
232 * In case, this matrix is not invertible, this matrix will be returned instead.
233 * @return The inverted matrix
234 * @see invert(), isSingular().
235 */
237
238 /**
239 * Inverts this matrix in place.
240 * @return True, if the matrix is not singular and could be inverted
241 * @see inverted(), solve().
242 */
243 bool invert();
244
245 /**
246 * Inverts the matrix and returns the result as parameter.
247 * @param invertedMatrix The resulting inverted matrix
248 * @return True, if the matrix is not singular and could be inverted
249 * @see inverted(), solve().
250 */
251 bool invert(SquareMatrixT3<T>& invertedMatrix) const;
252
253 /**
254 * Returns the determinant of the matrix.
255 * @return Matrix determinant
256 */
257 T determinant() const;
258
259 /**
260 * Returns the trace of the matrix which is the sum of the diagonal elements.
261 * @return Trace of the matrix
262 */
263 T trace() const;
264
265 /**
266 * Sets the matrix to the identity matrix.
267 * @see isIdentity().
268 */
269 inline void toIdentity();
270
271 /**
272 * Sets the matrix to a zero matrix.
273 * @see isNull();
274 */
275 inline void toNull();
276
277 /**
278 * Returns whether this matrix is an identity matrix.
279 * @return True, if so
280 * @see toIdentity().
281 */
282 bool isIdentity() const;
283
284 /**
285 * Returns whether this matrix is a zero matrix.
286 * @return True, if so
287 * @see toNull().
288 */
289 bool isNull() const;
290
291 /**
292 * Returns whether this matrix is singular (and thus cannot be inverted).
293 * A matrix is singular if the determinant of a matrix is 0.<br>
294 * @return True, if so
295 */
296 inline bool isSingular() const;
297
298 /**
299 * Returns true if this matrix is a similarity transformation.
300 * A similarity transformation has four degrees of freedom and contains a rotation, a scale, and a 2D translation and is not singular.<br>
301 * The 3x3 matrix representing the similarity transformation has the following layout:
302 * <pre>
303 * | a -b tx |
304 * | b a ty |
305 * | 0 0 1 |
306 * </pre>
307 * @return True, if this matrix is a similarity transformation, otherwise false
308 */
309 inline bool isSimilarity() const;
310
311 /**
312 * Returns true if this matrix is a affine transformation.
313 * In order to be considered affine, the matrix mustn't be singular and the last row must be equivalent to [0 0 1].
314 * @return True, if this matrix is an affine transformation, otherwise false
315 */
316 inline bool isAffine() const;
317
318 /**
319 * Returns true if this matrix is perspective transform/homography.
320 * In order to be considered a homography, the matrix mustn't be singular and the bottom-right matrix element must be nonzero.
321 * @return True, if this matrix is a homography, otherwise false
322 */
323 inline bool isHomography() const;
324
325 /**
326 * Returns whether this matrix is an orthonormal matrix.
327 * @param epsilon The epsilon threshold to be used, with range [0, infinity)
328 * @return True, if so
329 */
330 bool isOrthonormal(const T epsilon = NumericT<T>::eps()) const;
331
332 /**
333 * Returns whether this matrix is symmetric.
334 * @param epsilon The epsilon threshold to be used, with range [0, infinity)
335 * @return True, if so
336 */
337 inline bool isSymmetric(const T epsilon = NumericT<T>::eps()) const;
338
339 /**
340 * Returns whether two matrices are almost identical up to a specified epsilon.
341 * @param matrix Second matrix that will be checked
342 * @param eps The epsilon threshold to be used, with range [0, infinity)
343 * @return True, if so
344 */
345 inline bool isEqual(const SquareMatrixT3<T>& matrix, const T eps = NumericT<T>::eps()) const;
346
347 /**
348 * Returns the x axis which is the first column of the matrix.
349 * @return Vector with the first column
350 */
352
353 /**
354 * Returns the y axis which is the middle column of the matrix.
355 * @return Vector with the middle column
356 */
358
359 /**
360 * Returns the z axis which is the last column of the matrix.
361 * @return Vector with the last column
362 */
364
365 /**
366 * Returns the orthonormal matrix of this matrix by scaling the x-axis and adjusting y- and z-axis.
367 * This matrix must not be singular.
368 * @return An orthonormal version of this matrix
369 */
371
372 /**
373 * Determines the eigen values of this matrix.
374 * @param eigenValues The three resulting eigen values, sorted: the highest first
375 * @return True, if succeeded
376 */
377 bool eigenValues(T* eigenValues) const;
378
379 /**
380 * Performs an eigen value analysis.
381 * @param eigenValues The three resulting eigen values, sorted: the highest first
382 * @param eigenVectors The three corresponding eigen vectors
383 * @return True, if succeeded
384 */
385 bool eigenSystem(T* eigenValues,VectorT3<T>* eigenVectors) const;
386
387 /**
388 * Returns a 3d vector with values of the matrix diagonal.
389 * @return Vector with diagonal values
390 */
392
393 /**
394 * Solve a simple 3x3 system of linear equations: Ax = b
395 * Beware: The system of linear equations is assumed to be fully determined.
396 * @param b The right-hand side vector
397 * @param x The resulting solution vector
398 * @return True, if the solution is valid, otherwise false
399 * @see invert(), inverted().
400 */
401 inline bool solve(const VectorT3<T>& b, VectorT3<T>& x) const;
402
403 /**
404 * Calculates the sum of absolute value of elements
405 * @return Sum of absolute value of elements
406 */
407 inline T absSum() const;
408
409 /**
410 * Calculates the sum of elements
411 * @return Sum of elements
412 */
413 inline T sum() const;
414
415 /**
416 * Returns a pointer to the internal values.
417 * @return Pointer to the internal values
418 */
419 inline const T* data() const;
420
421 /**
422 * Returns a pointer to the internal values.
423 * @return Pointer to the internal values
424 */
425 inline T* data();
426
427 /**
428 * Copies the elements of this matrix to an array with floating point values of type U.
429 * @param arrayValues Array with 9 floating point values of type U receiving the elements of this matrix
430 * @param valuesRowAligned True, if the target values are stored in a row aligned order; False, if the target values are stored in a column aligned order
431 * @tparam U Floating point type
432 */
433 template <typename U>
434 void copyElements(U* arrayValues, const bool valuesRowAligned = false) const;
435
436 /**
437 * Copies the elements of this matrix to an array with floating point values of type T.
438 * @param arrayValues Array with 9 floating point values of type T receiving the elements of this matrix
439 * @param valuesRowAligned True, if the target values are stored in a row aligned order; False, if the target values are stored in a column aligned order
440 */
441 void copyElements(T* arrayValues, const bool valuesRowAligned = false) const;
442
443 /**
444 * Creates a skew symmetric 3x3 matrix for a specific vector.
445 * The skew symmetric matrix allows to calculate the cross product of the specified vector with a second vector by a matrix multiplication.<br>
446 * That means: skewSymmetricMatrix(vectorA) * vectorB == vectorA.cross(vectorB)<br>
447 * The final matrix has the following form for a vector (v0, v1, v2):
448 * <pre>
449 * | 0 -v2 v1 |
450 * | v2 0 -v0 |
451 * | -v1 v0 0 |
452 * </pre>
453 * @param vector The vector for which the skew symmetric matrix will be created
454 * @return Resulting matrix
455 */
457
458 /**
459 * Multiplies a 2D vector with this matrix (from the right).
460 * The 2D vector is interpreted as a 3D vector with third component equal to 1.<br>
461 * This function is equivalent with the corresponding multiplication operator but returns False if the dot product between the augmented vector and the last row is zero.<br>
462 * The multiplication result will be de-homogenized to provide a 2D vector result, if possible.<br>
463 * Actually this function does:
464 * @param vector The vector to be multiplied/transformed
465 * @param result The de-homogenized resulting 2D vector, if this function succeeds
466 * @return True, if the dot product between the augmented vector and the last row is non-zero; False, otherwise
467 * @see operator*(const VectorT2<T>& vector).
468 */
469 inline bool multiply(const VectorT2<T>& vector, VectorT2<T>& result) const;
470
471 /**
472 * Multiplies this transposed matrix with a second matrix.
473 * Actually, the following matrix will be returned: (*this).transposed() * right.<br>
474 * @param right Matrix to multiply
475 * @return Matrix product
476 */
478
479 /**
480 * Default copy assignment operator.
481 * @return Reference to this object
482 */
484
485 /**
486 * Returns whether two matrices are identical up to a small epsilon.
487 * @param matrix Right operand
488 * @return True, if so
489 */
490 bool operator==(const SquareMatrixT3<T>& matrix) const;
491
492 /**
493 * Returns whether two matrices are not identical up to a small epsilon.
494 * @param matrix Right operand
495 * @return True, if so
496 */
497 inline bool operator!=(const SquareMatrixT3<T>& matrix) const;
498
499 /**
500 * Adds two matrices.
501 * @param matrix Right operand
502 * @return Sum matrix
503 */
505
506 /**
507 * Adds and assigns two matrices.
508 * @param matrix Right operand
509 * @return Reference to this object
510 */
512
513 /**
514 * Subtracts two matrices.
515 * @param matrix Right operand
516 * @return Difference matrix
517 */
519
520 /**
521 * Subtracts and assigns two matrices.
522 * @param matrix Right operand
523 * @return Reference to this object
524 */
526
527 /**
528 * Returns the negative matrix of this matrix (all matrix elements are multiplied by -1).
529 * @return Resulting negative matrix
530 */
532
533 /**
534 * Multiplies two matrices.
535 * @param matrix Right operand
536 * @return Product matrix
537 */
538 OCEAN_FORCE_INLINE SquareMatrixT3<T> operator*(const SquareMatrixT3<T>& matrix) const;
539
540 /**
541 * Multiplies and assigns two matrices.
542 * @param matrix Right operand
543 * @return Reference to this object
544 */
545 OCEAN_FORCE_INLINE SquareMatrixT3<T>& operator*=(const SquareMatrixT3<T>& matrix);
546
547 /**
548 * Multiply operator for a 2D vector.
549 * The 2D vector is interpreted as a 3D vector with third component equal to 1.<br>
550 * The final result will be de-homogenized to provide a 2D vector result.<br>
551 * Beware the dot product between the last row and the (augmented) vector must not be zero!
552 * @param vector Right operand, vector to be multiplied from the right
553 * @return Resulting 2D vector
554 * @see multiply().
555 */
556 OCEAN_FORCE_INLINE VectorT2<T> operator*(const VectorT2<T>& vector) const;
557
558 /**
559 * Multiply operator for a 3D vector.
560 * @param vector Right operand
561 * @return Resulting 3D vector
562 */
563 OCEAN_FORCE_INLINE VectorT3<T> operator*(const VectorT3<T>& vector) const;
564
565 /**
566 * Multiplies this matrix with a scalar value.
567 * @param value Right operand
568 * @return Resulting matrix
569 */
570 OCEAN_FORCE_INLINE SquareMatrixT3<T> operator*(const T value) const;
571
572 /**
573 * Multiplies and assigns this matrix with a scalar value.
574 * @param value right operand
575 * @return Reference to this object
576 */
578
579 /**
580 * Element operator.
581 * Beware: No range check will be done!
582 * @param index The index of the element to return [0, 8]
583 * @return Specified element
584 */
585 inline T operator[](const unsigned int index) const;
586
587 /**
588 * Element operator.
589 * Beware: No range check will be done!
590 * @param index The index of the element to return [0, 8]
591 * @return Specified element
592 */
593 inline T& operator[](const unsigned int index);
594
595 /**
596 * Element operator.
597 * Beware: No range check will be done!
598 * @param row The row of the element to return [0, 2]
599 * @param column The column of the element to return [0, 2]
600 * @return Specified element
601 */
602 inline T operator()(const unsigned int row, const unsigned int column) const;
603
604 /**
605 * Element operator.
606 * Beware: No range check will be done!
607 * @param row The row of the element to return [0, 2]
608 * @param column The column of the element to return [0, 2]
609 * @return Specified element
610 */
611 inline T& operator()(const unsigned int row, const unsigned int column);
612
613 /**
614 * Element operator.
615 * Beware: No range check will be done!
616 * @param index The index of the element to return [0, 8]
617 * @return Specified element
618 */
619 inline T operator()(const unsigned int index) const;
620
621 /**
622 * Element operator.
623 * Beware: No range check will be done!
624 * @param index The index of the element to return [0, 8]
625 * @return Specified element
626 */
627 inline T& operator()(const unsigned int index);
628
629 /**
630 * Access operator.
631 * @return Pointer to the internal values
632 */
633 inline const T* operator()() const;
634
635 /**
636 * Access operator.
637 * @return Pointer to the internal values
638 */
639 inline T* operator()();
640
641 /**
642 * Hash function.
643 * @param matrix The matrix for which the hash value will be determined
644 * @return The resulting hash value
645 */
646 inline size_t operator()(const SquareMatrixT3<T>& matrix) const;
647
648 /**
649 * Returns the number of elements this matrix has.
650 * @return The number of elements, always 9
651 */
652 static inline size_t elements();
653
654 /**
655 * Multiplies several 2D vectors with a given 3x3 matrix.
656 * Each 2D vector is interpreted as a 3D vector with third component equal to 1.<br>
657 * The final result will be de-homogenized to provide a 2D vector result.
658 * @param matrix The matrix to be used for multiplication
659 * @param vectors The input vectors that will be multiplied, may be nullptr if number is 0
660 * @param results The resulting output (multiplied/transformed) vectors, with same number as the provided input vectors
661 * @param number The number of provided vectors (input and output), with range [0, infinity)
662 */
663 static void multiply(const SquareMatrixT3<T>& matrix, const VectorT2<T>* vectors, VectorT2<T>* results, const size_t number);
664
665 /**
666 * Multiplies several 3D vectors with a given 3x3 matrix.
667 * @param matrix The matrix to be used for multiplication
668 * @param vectors The input vectors that will be multiplied, may be nullptr if number is 0
669 * @param results The resulting output (multiplied/transformed) vectors, with same number as the provided input vectors
670 * @param number The number of provided vectors (input and output), with range [0, infinity)
671 */
672 static void multiply(const SquareMatrixT3<T>& matrix, const VectorT3<T>* vectors, VectorT3<T>* results, const size_t number);
673
674 /**
675 * Converts matrices with specific data type to matrices with different data type.
676 * @param matrices The matrices to convert
677 * @return The converted matrices
678 * @tparam U The element data type of the matrices to convert
679 */
680 template <typename U>
682
683 /**
684 * Converts matrices with specific data type to matrices with different data type.
685 * @param matrices The matrices to convert
686 * @param size The number of matrices to convert
687 * @return The converted matrices
688 * @tparam U The element data type of the matrices to convert
689 */
690 template <typename U>
691 static inline SquareMatricesT3<T> matrices2matrices(const SquareMatrixT3<U>* matrices, const size_t size);
692
693 protected:
694
695 /// The nine values of the matrix.
697};
698
699template <typename T>
701{
702 // nothing to do here
703}
704
705template <typename T>
706template <typename U>
708{
709 values_[0] = T(matrix.values_[0]);
710 values_[1] = T(matrix.values_[1]);
711 values_[2] = T(matrix.values_[2]);
712 values_[3] = T(matrix.values_[3]);
713 values_[4] = T(matrix.values_[4]);
714 values_[5] = T(matrix.values_[5]);
715 values_[6] = T(matrix.values_[6]);
716 values_[7] = T(matrix.values_[7]);
717 values_[8] = T(matrix.values_[8]);
718}
719
720template <typename T>
721SquareMatrixT3<T>::SquareMatrixT3(const bool setToIdentity)
722{
723 if (setToIdentity)
724 {
725 values_[0] = T(1.0);
726 values_[1] = T(0.0);
727 values_[2] = T(0.0);
728 values_[3] = T(0.0);
729 values_[4] = T(1.0);
730 values_[5] = T(0.0);
731 values_[6] = T(0.0);
732 values_[7] = T(0.0);
733 values_[8] = T(1.0);
734 }
735 else
736 {
737 values_[0] = T(0.0);
738 values_[1] = T(0.0);
739 values_[2] = T(0.0);
740 values_[3] = T(0.0);
741 values_[4] = T(0.0);
742 values_[5] = T(0.0);
743 values_[6] = T(0.0);
744 values_[7] = T(0.0);
745 values_[8] = T(0.0);
746 }
747}
748
749template <typename T>
751{
752 /**
753 * Rotation matrix around x-axis R(x):
754 * [ 1 0 0 ]
755 * [ 0 cos -sin ]
756 * [ 0 sin cos ]
757 */
758
759 /**
760 * Rotation matrix around y-axis R(y):
761 * [ cos 0 sin ]
762 * [ 0 1 0 ]
763 * [ -sin 0 cos ]
764 */
765
766 /**
767 * Rotation matrix around z-axis R(z):
768 * [ cos -sin 0 ]
769 * [ sin cos 0 ]
770 * [ 0 0 1 ]
771 */
772
773 /**
774 * Combined rotation matrix for R(y)R(x)R(z)
775 * [ cy cz + sx sy sz cz sx sy - cy sz cx sy ]
776 * [ cx sz cx cz -sx ]
777 * [ -cz sy + cy sx sz cy cz sx + sy sz cx cy ]
778 */
779
780 const T cx = NumericT<T>::cos(euler.pitch());
781 const T sx = NumericT<T>::sin(euler.pitch());
782
783 const T cy = NumericT<T>::cos(euler.yaw());
784 const T sy = NumericT<T>::sin(euler.yaw());
785
786 const T cz = NumericT<T>::cos(euler.roll());
787 const T sz = NumericT<T>::sin(euler.roll());
788
789 values_[0] = cy * cz + sx * sy * sz;
790 values_[1] = cx * sz;
791 values_[2] = -cz * sy + cy * sx * sz;
792 values_[3] = cz * sx * sy - cy * sz;
793 values_[4] = cx * cz;
794 values_[5] = cy * cz * sx + sy * sz;
795 values_[6] = cx * sy;
796 values_[7] = -sx;
797 values_[8] = cx * cy;
798
799 ocean_assert(NumericT<T>::isEqual(determinant(), T(1.0)));
800}
801
802template <typename T>
804{
805 // R(n, angle) = cos(angle) * I + (1 - cos(angle) * nn^T - sin(angle) * X(n)
806
807 ocean_assert(rotation.isValid());
808
809 const T cosValue = NumericT<T>::cos(rotation.angle());
810 const T cosValue1 = T(1.0) - cosValue;
811 const T sinValue = NumericT<T>::sin(rotation.angle());
812
813 const VectorT3<T> axis(rotation.axis());
814
815 const T xx = axis.x() * axis.x() * cosValue1;
816 const T yy = axis.y() * axis.y() * cosValue1;
817 const T zz = axis.z() * axis.z() * cosValue1;
818 const T xy = axis.x() * axis.y() * cosValue1;
819 const T xz = axis.x() * axis.z() * cosValue1;
820 const T yz = axis.y() * axis.z() * cosValue1;
821
822 const T nx = axis.x() * sinValue;
823 const T ny = axis.y() * sinValue;
824 const T nz = axis.z() * sinValue;
825
826 values_[0] = xx + cosValue;
827 values_[1] = xy + nz;
828 values_[2] = xz - ny;
829
830 values_[3] = xy - nz;
831 values_[4] = yy + cosValue;
832 values_[5] = yz + nx;
833
834 values_[6] = xz + ny;
835 values_[7] = yz - nx;
836 values_[8] = zz + cosValue;
837
838 ocean_assert(NumericT<T>::isEqual(determinant(), T(1.0)));
839}
840
841template <typename T>
843{
844 ocean_assert(quaternion.isValid());
845
846 const T xx = quaternion.x() * quaternion.x();
847 const T yy = quaternion.y() * quaternion.y();
848 const T zz = quaternion.z() * quaternion.z();
849
850 const T wx = quaternion.w() * quaternion.x();
851 const T wy = quaternion.w() * quaternion.y();
852 const T wz = quaternion.w() * quaternion.z();
853 const T xy = quaternion.x() * quaternion.y();
854 const T xz = quaternion.x() * quaternion.z();
855 const T yz = quaternion.y() * quaternion.z();
856
857 values_[0] = T(1.0) - T(2.0) * (yy + zz);
858 values_[1] = T(2.0) * (wz + xy);
859 values_[2] = T(2.0) * (xz - wy);
860
861 values_[3] = T(2.0) * (xy - wz);
862 values_[4] = T(1.0) - T(2.0) * (xx + zz);
863 values_[5] = T(2.0) * (wx + yz);
864
865 values_[6] = T(2.0) * (wy + xz);
866 values_[7] = T(2.0) * (yz - wx);
867 values_[8] = T(1.0) - T(2.0) * (xx + yy);
868
869 ocean_assert(NumericT<T>::isWeakEqual(determinant(), T(1.0)) && "the quaternion is not normalized");
870}
871
872template <typename T>
873SquareMatrixT3<T>::SquareMatrixT3(const T& m00, const T& m10, const T& m20, const T& m01, const T& m11, const T& m21, const T& m02, const T& m12, const T& m22)
874{
875 values_[0] = m00;
876 values_[1] = m10;
877 values_[2] = m20;
878
879 values_[3] = m01;
880 values_[4] = m11;
881 values_[5] = m21;
882
883 values_[6] = m02;
884 values_[7] = m12;
885 values_[8] = m22;
886}
887
888template <typename T>
890{
891 memcpy(values_, xAxis(), sizeof(T) * 3);
892 memcpy(values_ + 3, yAxis(), sizeof(T) * 3);
893 memcpy(values_ + 6, zAxis(), sizeof(T) * 3);
894}
895
896template <typename T>
898{
899 values_[0] = diagonal[0];
900 values_[1] = T(0.0);
901 values_[2] = T(0.0);
902 values_[3] = T(0.0);
903 values_[4] = diagonal[1];
904 values_[5] = T(0.0);
905 values_[6] = T(0.0);
906 values_[7] = T(0.0);
907 values_[8] = diagonal[2];
908}
909
910template <typename T>
911template <typename U>
913{
914 ocean_assert(arrayValues);
915
916 for (unsigned int n = 0u; n < 9u; ++n)
917 {
918 values_[n] = T(arrayValues[n]);
919 }
920}
921
922template <typename T>
924{
925 ocean_assert(arrayValues);
926 memcpy(values_, arrayValues, sizeof(T) * 9);
927}
928
929template <typename T>
930template <typename U>
931SquareMatrixT3<T>::SquareMatrixT3(const U* arrayValues, const bool valuesRowAligned)
932{
933 ocean_assert(arrayValues);
934
935 if (valuesRowAligned)
936 {
937 values_[0] = T(arrayValues[0]);
938 values_[1] = T(arrayValues[3]);
939 values_[2] = T(arrayValues[6]);
940 values_[3] = T(arrayValues[1]);
941 values_[4] = T(arrayValues[4]);
942 values_[5] = T(arrayValues[7]);
943 values_[6] = T(arrayValues[2]);
944 values_[7] = T(arrayValues[5]);
945 values_[8] = T(arrayValues[8]);
946
947 }
948 else
949 {
950 for (unsigned int n = 0u; n < 9u; ++n)
951 {
952 values_[n] = T(arrayValues[n]);
953 }
954 }
955}
956
957template <typename T>
958SquareMatrixT3<T>::SquareMatrixT3(const T* arrayValues, const bool valuesRowAligned)
959{
960 ocean_assert(arrayValues);
961
962 if (valuesRowAligned)
963 {
964 values_[0] = arrayValues[0];
965 values_[1] = arrayValues[3];
966 values_[2] = arrayValues[6];
967 values_[3] = arrayValues[1];
968 values_[4] = arrayValues[4];
969 values_[5] = arrayValues[7];
970 values_[6] = arrayValues[2];
971 values_[7] = arrayValues[5];
972 values_[8] = arrayValues[8];
973 }
974 else
975 {
976 memcpy(values_, arrayValues, sizeof(T) * 9);
977 }
978}
979
980template <typename T>
982{
983 memcpy(values_, transformation(), sizeof(T) * 3);
984 memcpy(values_ + 3, transformation() + 4, sizeof(T) * 3);
985 memcpy(values_ + 6, transformation() + 8, sizeof(T) * 3);
986}
987
988template <typename T>
990{
991 memcpy(values_, transformation(), sizeof(T) * 3);
992 memcpy(values_ + 3, transformation() + 4, sizeof(T) * 3);
993 memcpy(values_ + 6, transformation() + 8, sizeof(T) * 3);
994}
995
996template<typename T>
997inline bool SquareMatrixT3<T>::solve(const VectorT3<T>& b, VectorT3<T>& x) const
998{
999 // Solve the system of linear equations, Ax=b, using Cramer's rule
1000 //
1001 // [a0 a3 a6] [b0]
1002 // A = [a1 a4 a7], b = [b1]
1003 // [a2 a5 a8] [b2]
1004 //
1005 // d = det(A)
1006 //
1007 // [b0 a3 a6] [a0 b0 a6] [a0 a3 b0]
1008 // d0 = det( [b1 a4 a7] ), d1 = det( [a1 b1 a7] ), d2 = det( [a1 a4 b1] )
1009 // [b2 a5 a8] [a2 b2 a8] [a2 a5 b2]
1010 //
1011 // [d0 / d]
1012 // x = [d1 / d]
1013 // [d2 / d]
1014 const T d = determinant();
1015
1017 {
1018 const T d0 = b[0] * (values_[4] * values_[8] - values_[5] * values_[7]) + b[1] * (values_[5] * values_[6] - values_[3] * values_[8]) + b[2] * (values_[3] * values_[7] - values_[4] * values_[6]);
1019 const T d1 = values_[0] * ( b[1] * values_[8] - b[2] * values_[7]) + values_[1] * ( b[2] * values_[6] - b[0] * values_[8]) + values_[2] * ( b[0] * values_[7] - b[1] * values_[6]);
1020 const T d2 = values_[0] * (values_[4] * b[2] - values_[5] * b[1]) + values_[1] * (values_[5] * b[0] - values_[3] * b[2]) + values_[2] * (values_[3] * b[1] - values_[4] * b[0]);
1021
1022 const T invD = T(1) / d;
1023
1024 x[0] = d0 * invD;
1025 x[1] = d1 * invD;
1026 x[2] = d2 * invD;
1027
1028 return true;
1029 }
1030
1031 return false;
1032}
1033
1034template<typename T>
1036{
1037 return (NumericT<T>::abs(values_[0]) + NumericT<T>::abs(values_[1]) + NumericT<T>::abs(values_[2]) + NumericT<T>::abs(values_[3]) + NumericT<T>::abs(values_[4]) + NumericT<T>::abs(values_[5]) + NumericT<T>::abs(values_[6]) + NumericT<T>::abs(values_[7]) + NumericT<T>::abs(values_[8]));
1038}
1039
1040template<typename T>
1042{
1043 return (values_[0] + values_[1] + values_[2] + values_[3] + values_[4] + values_[5] + values_[6] + values_[7] + values_[8]);
1044}
1045
1046template <typename T>
1047inline const T* SquareMatrixT3<T>::data() const
1048{
1049 return values_;
1050}
1051
1052template <typename T>
1054{
1055 return values_;
1056}
1057
1058template <typename T>
1059inline bool SquareMatrixT3<T>::operator!=(const SquareMatrixT3<T>& matrix) const
1060{
1061 return !(*this == matrix);
1062}
1063
1064template <typename T>
1066{
1067 *this = *this * matrix;
1068 return *this;
1069}
1070
1071template <typename T>
1072inline T SquareMatrixT3<T>::operator[](const unsigned int index) const
1073{
1074 ocean_assert(index < 9u);
1075 return values_[index];
1076}
1077
1078template <typename T>
1079inline T& SquareMatrixT3<T>::operator[](const unsigned int index)
1080{
1081 ocean_assert(index < 9u);
1082 return values_[index];
1083}
1084
1085template <typename T>
1086inline T SquareMatrixT3<T>::operator()(const unsigned int row, const unsigned int column) const
1087{
1088 ocean_assert(row < 3u && column < 3u);
1089 return values_[column * 3 + row];
1090}
1091
1092template <typename T>
1093inline T& SquareMatrixT3<T>::operator()(const unsigned int row, const unsigned int column)
1094{
1095 ocean_assert(row < 3u && column < 3u);
1096 return values_[column * 3 + row];
1097}
1098
1099template <typename T>
1100inline T SquareMatrixT3<T>::operator()(const unsigned int index) const
1101{
1102 ocean_assert(index < 9u);
1103 return values_[index];
1104}
1105
1106template <typename T>
1107inline T& SquareMatrixT3<T>::operator()(const unsigned int index)
1108{
1109 ocean_assert(index < 9u);
1110 return values_[index];
1111}
1112
1113template <typename T>
1114inline const T* SquareMatrixT3<T>::operator()() const
1115{
1116 return values_;
1117}
1118
1119template <typename T>
1121{
1122 return values_;
1123}
1124
1125template <typename T>
1126inline size_t SquareMatrixT3<T>::operator()(const SquareMatrixT3<T>& matrix) const
1127{
1128 size_t seed = std::hash<T>{}(matrix.values_[0]);
1129
1130 for (unsigned int n = 1u; n < 9u; ++n)
1131 {
1132 seed ^= std::hash<T>{}(matrix.values_[n]) + 0x9e3779b9 + (seed << 6) + (seed >> 2);
1133 }
1134
1135 return seed;
1136}
1137
1138template <typename T>
1140{
1141 return 9;
1142}
1143
1144template <typename T>
1146{
1147 SquareMatrixT3<T> result(*this);
1148
1149 result.values_[1] = values_[3];
1150 result.values_[3] = values_[1];
1151
1152 result.values_[2] = values_[6];
1153 result.values_[6] = values_[2];
1154
1155 result.values_[5] = values_[7];
1156 result.values_[7] = values_[5];
1157
1158 return result;
1159}
1160
1161template <typename T>
1163{
1164 SquareMatrixT3<T> tmp(*this);
1165
1166 values_[3] = tmp.values_[1];
1167 values_[1] = tmp.values_[3];
1168
1169 values_[6] = tmp.values_[2];
1170 values_[2] = tmp.values_[6];
1171
1172 values_[7] = tmp.values_[5];
1173 values_[5] = tmp.values_[7];
1174}
1175
1176template <typename T>
1178{
1179 SquareMatrixT3<T> invertedMatrix;
1180
1181 if (!invert(invertedMatrix))
1182 {
1183 ocean_assert(false && "Could not invert the matrix.");
1184 return *this;
1185 }
1186
1187 return invertedMatrix;
1188}
1189
1190template <typename T>
1192{
1193 SquareMatrixT3<T> invertedMatrix;
1194
1195 if (!invert(invertedMatrix))
1196 {
1197 return false;
1198 }
1199
1200 *this = invertedMatrix;
1201
1202 return true;
1203}
1204
1205template <typename T>
1207{
1208 // calculate determinant
1209 const T v48 = values_[4] * values_[8];
1210 const T v57 = values_[5] * values_[7];
1211 const T v56 = values_[5] * values_[6];
1212 const T v38 = values_[3] * values_[8];
1213 const T v37 = values_[3] * values_[7];
1214 const T v46 = values_[4] * values_[6];
1215
1216 const T v48_57 = v48 - v57;
1217 const T v56_38 = v56 - v38;
1218 const T v37_46 = v37 - v46;
1219
1220 const T det = values_[0] * v48_57 + values_[1] * v56_38 + values_[2] * v37_46;
1221
1222 if (NumericT<T>::isEqualEps(det))
1223 {
1224 return false;
1225 }
1226
1227 const T factor = T(1.0) / det;
1228
1229 invertedMatrix.values_[0] = (v48_57) * factor;
1230 invertedMatrix.values_[1] = (values_[2] * values_[7] - values_[1] * values_[8]) * factor;
1231 invertedMatrix.values_[2] = (values_[1] * values_[5] - values_[2] * values_[4]) * factor;
1232
1233 invertedMatrix.values_[3] = (v56_38) * factor;
1234 invertedMatrix.values_[4] = (values_[0] * values_[8] - values_[2] * values_[6]) * factor;
1235 invertedMatrix.values_[5] = (values_[2] * values_[3] - values_[0] * values_[5]) * factor;
1236
1237 invertedMatrix.values_[6] = (v37_46) * factor;
1238 invertedMatrix.values_[7] = (values_[1] * values_[6] - values_[0] * values_[7]) * factor;
1239 invertedMatrix.values_[8] = (values_[0] * values_[4] - values_[1] * values_[3]) * factor;
1240
1241#ifdef OCEAN_INTENSIVE_DEBUG
1242 if (!std::is_same<T, float>::value)
1243 {
1244 const SquareMatrixT3<T> test(*this * invertedMatrix);
1245 const SquareMatrixT3<T> entity(true);
1246
1247 T sqrDistance = T(0);
1248 for (unsigned int n = 0; n < 9u; ++n)
1249 {
1250 sqrDistance += NumericT<T>::sqr(test[n] - entity[n]);
1251 }
1252
1253 const T distance = NumericT<T>::sqrt(sqrDistance * T(0.111111111111111111)); // 1 / 9
1254
1255 if (NumericT<T>::isWeakEqualEps(distance) == false)
1256 {
1257 T absoluteAverageEnergy = 0;
1258 for (unsigned int n = 0u; n < 9u; ++n)
1259 {
1260 absolusteAverageEnergy += NumericT<T>::abs(values[n]);
1261 }
1262
1263 absoluteAverageEnergy *= T(0.111111111111111111); // 1 / 9
1264
1265 // we expect/accept for each magnitude (larger than 1) a zero-inaccuracy of one magnitude (and we again comare it with the weak eps)
1266
1267 if (absoluteAverageEnergy <= 1)
1268 {
1269 ocean_assert_accuracy(!"This should never happen!");
1270 }
1271 else
1272 {
1273 const T adjustedDistance = distance / absoluteAverageEnergy;
1274 ocean_assert_accuracy(NumericT<T>::isWeakEqualEps(adjustedDistance));
1275 }
1276 }
1277 }
1278#endif // OCEAN_DEBUG
1279
1280 return true;
1281}
1282
1283template <typename T>
1285{
1286 return values_[0] * (values_[4] * values_[8] - values_[5] * values_[7])
1287 + values_[1] * (values_[5] * values_[6] - values_[3] * values_[8])
1288 + values_[2] * (values_[3] * values_[7] - values_[4] * values_[6]);
1289}
1290
1291template <typename T>
1293{
1294 return values_[0] + values_[4] + values_[8];
1295}
1296
1297template <typename T>
1299{
1300 values_[0] = T(1.0);
1301 values_[1] = T(0.0);
1302 values_[2] = T(0.0);
1303 values_[3] = T(0.0);
1304 values_[4] = T(1.0);
1305 values_[5] = T(0.0);
1306 values_[6] = T(0.0);
1307 values_[7] = T(0.0);
1308 values_[8] = T(1.0);
1309}
1310
1311template <typename T>
1313{
1314 values_[0] = T(0.0);
1315 values_[1] = T(0.0);
1316 values_[2] = T(0.0);
1317 values_[3] = T(0.0);
1318 values_[4] = T(0.0);
1319 values_[5] = T(0.0);
1320 values_[6] = T(0.0);
1321 values_[7] = T(0.0);
1322 values_[8] = T(0.0);
1323}
1324
1325template <typename T>
1327{
1328 return NumericT<T>::isEqual(values_[0], 1) && NumericT<T>::isEqualEps(values_[1]) && NumericT<T>::isEqualEps(values_[2])
1329 && NumericT<T>::isEqualEps(values_[3]) && NumericT<T>::isEqual(values_[4], 1) && NumericT<T>::isEqualEps(values_[5])
1330 && NumericT<T>::isEqualEps(values_[6]) && NumericT<T>::isEqualEps(values_[7]) && NumericT<T>::isEqual(values_[8], 1);
1331}
1332
1333template <typename T>
1335{
1336 return NumericT<T>::isEqualEps(values_[0]) && NumericT<T>::isEqualEps(values_[1]) && NumericT<T>::isEqualEps(values_[2])
1337 && NumericT<T>::isEqualEps(values_[3]) && NumericT<T>::isEqualEps(values_[4]) && NumericT<T>::isEqualEps(values_[5])
1338 && NumericT<T>::isEqualEps(values_[6]) && NumericT<T>::isEqualEps(values_[7]) && NumericT<T>::isEqualEps(values_[8]);
1339}
1340
1341template <typename T>
1343{
1344 return NumericT<T>::isEqualEps(determinant());
1345}
1346
1347template <typename T>
1349{
1350 return NumericT<T>::isEqual(values_[0], values_[4]) && NumericT<T>::isEqual(values_[1], -values_[3]) && NumericT<T>::isEqualEps(values_[2]) && NumericT<T>::isEqualEps(values_[5]) && NumericT<T>::isEqual(values_[8], 1) && !isSingular();
1351}
1352
1353template <typename T>
1355{
1356 return NumericT<T>::isEqualEps(values_[2]) && NumericT<T>::isEqualEps(values_[5]) && NumericT<T>::isEqual(values_[8], 1) && !isSingular();
1357}
1358
1359template <typename T>
1361{
1362 return NumericT<T>::isNotEqualEps(values_[8]) && !isSingular();
1363}
1364
1365template <typename T>
1366bool SquareMatrixT3<T>::isOrthonormal(const T epsilon) const
1367{
1368 ocean_assert(epsilon >= 0);
1369
1370 const VectorT3<T> xAxis(values_);
1371 const VectorT3<T> yAxis(values_ + 3);
1372 const VectorT3<T> zAxis(values_ + 6);
1373
1374 return NumericT<T>::isEqual(xAxis * yAxis, 0, epsilon) && NumericT<T>::isEqual(xAxis * zAxis, 0, epsilon) && NumericT<T>::isEqual(yAxis * zAxis, 0, epsilon)
1375 && NumericT<T>::isEqual(xAxis.length(), 1, epsilon) && NumericT<T>::isEqual(yAxis.length(), 1, epsilon) && NumericT<T>::isEqual(zAxis.length(), 1, epsilon);
1376}
1377
1378template <typename T>
1379bool SquareMatrixT3<T>::isSymmetric(const T epsilon) const
1380{
1381 ocean_assert(epsilon >= T(0));
1382
1383 return NumericT<T>::isEqual(values_[1], values_[3], epsilon) && NumericT<T>::isEqual(values_[2], values_[6], epsilon) && NumericT<T>::isEqual(values_[5], values_[7], epsilon);
1384}
1385
1386template <typename T>
1387inline bool SquareMatrixT3<T>::isEqual(const SquareMatrixT3<T>& matrix, const T eps) const
1388{
1389 return NumericT<T>::isEqual(values_[0], matrix.values_[0], eps) && NumericT<T>::isEqual(values_[1], matrix.values_[1], eps)
1390 && NumericT<T>::isEqual(values_[2], matrix.values_[2], eps) && NumericT<T>::isEqual(values_[3], matrix.values_[3], eps)
1391 && NumericT<T>::isEqual(values_[4], matrix.values_[4], eps) && NumericT<T>::isEqual(values_[5], matrix.values_[5], eps)
1392 && NumericT<T>::isEqual(values_[6], matrix.values_[6], eps) && NumericT<T>::isEqual(values_[7], matrix.values_[7], eps)
1393 && NumericT<T>::isEqual(values_[8], matrix.values_[8], eps);
1394}
1395
1396template <typename T>
1398{
1399 return VectorT3<T>(values_);
1400}
1401
1402template <typename T>
1404{
1405 return VectorT3<T>(values_ + 3);
1406}
1407
1408template <typename T>
1410{
1411 return VectorT3<T>(values_ + 6);
1412}
1413
1414template <typename T>
1416{
1417 ocean_assert(!isSingular());
1418
1419 VectorT3<T> xAxis(values_);
1420 VectorT3<T> yAxis(values_ + 3);
1421 VectorT3<T> zAxis(values_ + 6);
1422
1423 // X scale factor
1424 const T xScale = xAxis.length();
1425 // normalize x axis
1426 xAxis /= xScale;
1427
1428 // xy shear factor
1429 const T xyShear = xAxis * yAxis;
1430
1431 // compute orthogonal y axis
1432 yAxis -= xAxis * xyShear;
1433
1434 // y scale factor
1435 const T yScale = yAxis.length();
1436 // normalize y axis
1437 yAxis /= yScale;
1438
1439 // xz shear
1440 const T xzShear = xAxis * zAxis;
1441 // compute orthogonal z axis
1442 zAxis -= xAxis * xzShear;
1443
1444 // yz shear
1445 const T yzShear = yAxis * zAxis;
1446 // compute orthogonal y axis
1447 zAxis -= yAxis * yzShear;
1448
1449 // z scale factor
1450 const T zScale = zAxis.length();
1451 // normalize z axis
1452 zAxis /= zScale;
1453
1454#ifdef OCEAN_INTENSIVE_DEBUG
1455 if (std::is_same<T, double>::value)
1456 {
1457 ocean_assert(NumericT<T>::isEqualEps(xAxis * yAxis));
1458 ocean_assert(NumericT<T>::isEqualEps(xAxis * zAxis));
1459 ocean_assert(NumericT<T>::isEqualEps(yAxis * zAxis));
1460 ocean_assert(SquareMatrixT3<T>(xAxis, yAxis, zAxis).isOrthonormal());
1461 }
1462#endif
1463
1464 return SquareMatrixT3<T>(xAxis, yAxis, zAxis);
1465}
1466
1467template <typename T>
1468bool SquareMatrixT3<T>::eigenValues(T* eigenValues) const
1469{
1470 ocean_assert(eigenValues);
1471
1472 /**
1473 * Computation of the characteristic polynomial
1474 * <pre>
1475 *
1476 * [ a b c ]
1477 * A = [ d e f ]
1478 * [ g h i ]
1479 *
1480 * [ a-x b c ]
1481 * A - x * E = [ d e-x f ]
1482 * [ g h i-x ]
1483 *
1484 * Polynomial = Det|A - x * E| = 0
1485 * = (a - x) * (e - x) * (i - x) + bfg + cdh - g * (e - x) * c - h * f * (a - x) - (i - x) * d * b
1486 * = (ae - a * x - e * x + x^2) * (i - x) + bfg + cdh - g * (ec - c * x) - h * (fa - f * x) - d * (ib - x * b)
1487 * = aei - ae * x - ai * x + a * x^2 - ei * x + e * x^2 + i * x^2 - x^3 + bfg + cdh - gec + gc * x - hfa + hf * x - dib + db * x
1488 * = -x^3 + (a + e + i) * x^2 + (-ae - ai - ei + gc + hf + db) * x + aei + bfg + cdh - gec - hfa - dib
1489 * = x^3 - (a + e + i) * x^2 - (-ae - ai - ei + gc + hf + db) * x - aei - bfg - cdh + gec + hfa + dib
1490 * = a1x^3 + a2x^2 + a3x + a4 = 0
1491 * </pre>
1492 */
1493
1494 const T a = values_[0];
1495 const T d = values_[1];
1496 const T g = values_[2];
1497
1498 const T b = values_[3];
1499 const T e = values_[4];
1500 const T h = values_[5];
1501
1502 const T c = values_[6];
1503 const T f = values_[7];
1504 const T i = values_[8];
1505
1506 const T a1 = T(1.0);
1507 const T a2 = -(a + e + i);
1508 const T a3 = -(-a * e - a * i - e * i + g * c + h * f + d * b);
1509 const T a4 = -a * e * i - b * f * g - c * d * h + g * e * c + h * f * a + d * i * b;
1510
1511 if (EquationT<T>::solveCubic(a1, a2, a3, a4, eigenValues[0], eigenValues[1], eigenValues[2]) != 3u)
1512 {
1513 return false;
1514 }
1515
1516 Utilities::sortHighestToFront3(eigenValues[0], eigenValues[1], eigenValues[2]);
1517 return true;
1518}
1519
1520template <typename T>
1521bool SquareMatrixT3<T>::eigenSystem(T* eigenValues, VectorT3<T>* eigenVectors) const
1522{
1523 ocean_assert(eigenValues != nullptr && eigenVectors != nullptr);
1524
1525 /**
1526 * Computation of the characteristic polynomial
1527 * <pre>
1528 * [ a b c ]
1529 * A = [ d e f ]
1530 * [ g h i ]
1531 *
1532 * [ a-x b c ]
1533 * A - x * E = [ d e-x f ]
1534 * [ g h i-x ]
1535 *
1536 * Polynomial = Det|A - x * E| = 0
1537 * = (a - x) * (e - x) * (i - x) + bfg + cdh - g * (e - x) * c - h * f * (a - x) - (i - x) * d * b
1538 * = (ae - a * x - e * x + x^2) * (i - x) + bfg + cdh - g * (ec - c * x) - h * (fa - f * x) - d * (ib - x * b)
1539 * = aei - ae * x - ai * x + a * x^2 - ei * x + e * x^2 + i * x^2 - x^3 + bfg + cdh - gec + gc * x - hfa + hf * x - dib + db * x
1540 * = -x^3 + (a + e + i) * x^2 + (-ae - ai - ei + gc + hf + db) * x + aei + bfg + cdh - gec - hfa - dib
1541 * = x^3 - (a + e + i) * x^2 - (-ae - ai - ei + gc + hf + db) * x - aei - bfg - cdh + gec + hfa + dib
1542 * = a1x^3 + a2x^2 + a3x + a4 = 0
1543 * </pre>
1544 */
1545
1546 const T a = values_[0];
1547 const T d = values_[1];
1548 const T g = values_[2];
1549
1550 const T b = values_[3];
1551 const T e = values_[4];
1552 const T h = values_[5];
1553
1554 const T c = values_[6];
1555 const T f = values_[7];
1556 const T i = values_[8];
1557
1558 const T a1 = T(1.0);
1559 const T a2 = -(a + e + i);
1560 const T a3 = -(-a * e - a * i - e * i + g * c + h * f + d * b);
1561 const T a4 = -a * e * i - b * f * g - c * d * h + g * e * c + h * f * a + d * i * b;
1562
1563 if (EquationT<T>::solveCubic(a1, a2, a3, a4, eigenValues[0], eigenValues[1], eigenValues[2]) != 3u)
1564 {
1565 return false;
1566 }
1567
1568 Utilities::sortHighestToFront3(eigenValues[0], eigenValues[1], eigenValues[2]);
1569
1570 /**
1571 * <pre>
1572 * Determination of the eigen vectors (vx, vy, vz):
1573 * [ a-x b c ] [ vx ]
1574 * A - x * E = [ d e-x f ] * [ vy ] = 0
1575 * [ g h i-x ] [ vz ]
1576 * </pre>
1577 */
1578
1579 for (unsigned int n = 0u; n < 3u; ++n)
1580 {
1581 const VectorT3<T> row0(a - eigenValues[n], b, c);
1582 const VectorT3<T> row1(d, e - eigenValues[n], f);
1583 const VectorT3<T> row2(g, h, i - eigenValues[n]);
1584
1585 VectorT3<T> candidate0(row0.cross(row1));
1586 VectorT3<T> candidate1(row0.cross(row2));
1587 VectorT3<T> candidate2(row1.cross(row2));
1588
1589 T sqrCandidate0(candidate0.sqr());
1590 T sqrCandidate1(candidate1.sqr());
1591 T sqrCandidate2(candidate2.sqr());
1592
1593 Utilities::sortHighestToFront3(sqrCandidate0, sqrCandidate1, sqrCandidate2, candidate0, candidate1, candidate2);
1594
1595 // if all rows (row0, row1 and row2) are parallel, we can decide any vector that is perpendicular to these rows
1596 if (sqrCandidate0 < NumericT<T>::eps() * NumericT<T>::eps())
1597 {
1598 // find one row that is not a null row
1599
1600 candidate0 = row0;
1601 candidate1 = row1;
1602 candidate2 = row2;
1603
1604 sqrCandidate0 = candidate0.sqr();
1605 sqrCandidate1 = candidate1.sqr();
1606 sqrCandidate2 = candidate2.sqr();
1607
1608 Utilities::sortHighestToFront3(sqrCandidate0, sqrCandidate1, sqrCandidate2, candidate0, candidate1, candidate2);
1609
1610 ocean_assert(NumericT<T>::isNotEqualEps(candidate0.length()));
1611 eigenVectors[n] = candidate0.perpendicular();
1612 }
1613 else
1614 {
1615 eigenVectors[n] = candidate0;
1616 }
1617
1618 eigenVectors[n].normalize();
1619 }
1620
1621 return true;
1622}
1623
1624template <typename T>
1626{
1627 return VectorT3<T>(values_[0], values_[4], values_[8]);
1628}
1629
1630template <typename T>
1631template <typename U>
1632void SquareMatrixT3<T>::copyElements(U* arrayValues, const bool valuesRowAligned) const
1633{
1634 ocean_assert(arrayValues != nullptr);
1635
1636 if (valuesRowAligned)
1637 {
1638 // this matrix and the provided array are both column aligned
1639 // thus, we can simply copy the data
1640
1641 arrayValues[0] = U(values_[0]);
1642 arrayValues[1] = U(values_[3]);
1643 arrayValues[2] = U(values_[6]);
1644
1645 arrayValues[3] = U(values_[1]);
1646 arrayValues[4] = U(values_[4]);
1647 arrayValues[5] = U(values_[7]);
1648
1649 arrayValues[6] = U(values_[2]);
1650 arrayValues[7] = U(values_[5]);
1651 arrayValues[8] = U(values_[8]);
1652 }
1653 else
1654 {
1655 // this matrix and the provided array are both column aligned
1656 // thus, we can simply copy the data
1657
1658 arrayValues[0] = U(values_[0]);
1659 arrayValues[1] = U(values_[1]);
1660 arrayValues[2] = U(values_[2]);
1661 arrayValues[3] = U(values_[3]);
1662 arrayValues[4] = U(values_[4]);
1663 arrayValues[5] = U(values_[5]);
1664 arrayValues[6] = U(values_[6]);
1665 arrayValues[7] = U(values_[7]);
1666 arrayValues[8] = U(values_[8]);
1667 }
1668}
1669
1670template <typename T>
1671void SquareMatrixT3<T>::copyElements(T* arrayValues, const bool valuesRowAligned) const
1672{
1673 ocean_assert(arrayValues != nullptr);
1674
1675 if (valuesRowAligned)
1676 {
1677 // this matrix and the provided array are both column aligned
1678 // thus, we can simply copy the data
1679
1680 arrayValues[0] = values_[0];
1681 arrayValues[1] = values_[3];
1682 arrayValues[2] = values_[6];
1683
1684 arrayValues[3] = values_[1];
1685 arrayValues[4] = values_[4];
1686 arrayValues[5] = values_[7];
1687
1688 arrayValues[6] = values_[2];
1689 arrayValues[7] = values_[5];
1690 arrayValues[8] = values_[8];
1691 }
1692 else
1693 {
1694 // this matrix and the provided array are both column aligned
1695 // thus, we can simply copy the data
1696
1697 memcpy(arrayValues, values_, sizeof(T) * 9);
1698 }
1699}
1700
1701template <typename T>
1703{
1704 return SquareMatrixT3<T>(0, vector(2), -vector(1), -vector(2), 0, vector(0), vector(1), -vector(0), 0);
1705}
1706
1707template <typename T>
1708inline bool SquareMatrixT3<T>::multiply(const VectorT2<T>& vector, VectorT2<T>& result) const
1709{
1710 /**
1711 * | x' | | 0 3 6 | | x |
1712 * | y' | = | 1 4 7 | * | y |
1713 * | 1 | | 2 5 8 | | 1 |
1714 */
1715
1716 const T z = values_[2] * vector[0] + values_[5] * vector[1] + values_[8];
1717
1719 {
1720 const T factor = T(1) / z;
1721
1722 result = VectorT2<T>((values_[0] * vector[0] + values_[3] * vector[1] + values_[6]) * factor,
1723 (values_[1] * vector[0] + values_[4] * vector[1] + values_[7]) * factor);
1724
1725 return true;
1726 }
1727
1728 return false;
1729}
1730
1731template <typename T>
1733{
1734 SquareMatrixT3<T> result;
1735
1736 result.values_[0] = values_[0] * matrix.values_[0] + values_[1] * matrix.values_[1] + values_[2] * matrix.values_[2];
1737 result.values_[1] = values_[3] * matrix.values_[0] + values_[4] * matrix.values_[1] + values_[5] * matrix.values_[2];
1738 result.values_[2] = values_[6] * matrix.values_[0] + values_[7] * matrix.values_[1] + values_[8] * matrix.values_[2];
1739 result.values_[3] = values_[0] * matrix.values_[3] + values_[1] * matrix.values_[4] + values_[2] * matrix.values_[5];
1740 result.values_[4] = values_[3] * matrix.values_[3] + values_[4] * matrix.values_[4] + values_[5] * matrix.values_[5];
1741 result.values_[5] = values_[6] * matrix.values_[3] + values_[7] * matrix.values_[4] + values_[8] * matrix.values_[5];
1742 result.values_[6] = values_[0] * matrix.values_[6] + values_[1] * matrix.values_[7] + values_[2] * matrix.values_[8];
1743 result.values_[7] = values_[3] * matrix.values_[6] + values_[4] * matrix.values_[7] + values_[5] * matrix.values_[8];
1744 result.values_[8] = values_[6] * matrix.values_[6] + values_[7] * matrix.values_[7] + values_[8] * matrix.values_[8];
1745
1746#ifdef OCEAN_INTENSIVE_DEBUG
1747 const SquareMatrixT3<T> debugMatrix(transposed() * matrix);
1748 ocean_assert(debugMatrix == result);
1749#endif
1750
1751 return result;
1752}
1753
1754template <typename T>
1756{
1757 return isEqual(matrix);
1758}
1759
1760template <typename T>
1762{
1763 SquareMatrixT3<T> result(*this);
1764
1765 result.values_[0] += matrix.values_[0];
1766 result.values_[1] += matrix.values_[1];
1767 result.values_[2] += matrix.values_[2];
1768 result.values_[3] += matrix.values_[3];
1769 result.values_[4] += matrix.values_[4];
1770 result.values_[5] += matrix.values_[5];
1771 result.values_[6] += matrix.values_[6];
1772 result.values_[7] += matrix.values_[7];
1773 result.values_[8] += matrix.values_[8];
1774
1775 return result;
1776}
1777
1778template <typename T>
1780{
1781 values_[0] += matrix.values_[0];
1782 values_[1] += matrix.values_[1];
1783 values_[2] += matrix.values_[2];
1784 values_[3] += matrix.values_[3];
1785 values_[4] += matrix.values_[4];
1786 values_[5] += matrix.values_[5];
1787 values_[6] += matrix.values_[6];
1788 values_[7] += matrix.values_[7];
1789 values_[8] += matrix.values_[8];
1790
1791 return *this;
1792}
1793
1794template <typename T>
1796{
1797 SquareMatrixT3<T> result(*this);
1798
1799 result.values_[0] -= matrix.values_[0];
1800 result.values_[1] -= matrix.values_[1];
1801 result.values_[2] -= matrix.values_[2];
1802 result.values_[3] -= matrix.values_[3];
1803 result.values_[4] -= matrix.values_[4];
1804 result.values_[5] -= matrix.values_[5];
1805 result.values_[6] -= matrix.values_[6];
1806 result.values_[7] -= matrix.values_[7];
1807 result.values_[8] -= matrix.values_[8];
1808
1809 return result;
1810}
1811
1812template <typename T>
1814{
1815 values_[0] -= matrix.values_[0];
1816 values_[1] -= matrix.values_[1];
1817 values_[2] -= matrix.values_[2];
1818 values_[3] -= matrix.values_[3];
1819 values_[4] -= matrix.values_[4];
1820 values_[5] -= matrix.values_[5];
1821 values_[6] -= matrix.values_[6];
1822 values_[7] -= matrix.values_[7];
1823 values_[8] -= matrix.values_[8];
1824
1825 return *this;
1826}
1827
1828template <typename T>
1830{
1831 SquareMatrixT3<T> result;
1832
1833 result.values_[0] = -values_[0];
1834 result.values_[1] = -values_[1];
1835 result.values_[2] = -values_[2];
1836 result.values_[3] = -values_[3];
1837 result.values_[4] = -values_[4];
1838 result.values_[5] = -values_[5];
1839 result.values_[6] = -values_[6];
1840 result.values_[7] = -values_[7];
1841 result.values_[8] = -values_[8];
1842
1843 return result;
1844}
1845
1846template <typename T>
1848{
1849 SquareMatrixT3<T> result;
1850
1851 result.values_[0] = values_[0] * matrix.values_[0] + values_[3] * matrix.values_[1] + values_[6] * matrix.values_[2];
1852 result.values_[1] = values_[1] * matrix.values_[0] + values_[4] * matrix.values_[1] + values_[7] * matrix.values_[2];
1853 result.values_[2] = values_[2] * matrix.values_[0] + values_[5] * matrix.values_[1] + values_[8] * matrix.values_[2];
1854 result.values_[3] = values_[0] * matrix.values_[3] + values_[3] * matrix.values_[4] + values_[6] * matrix.values_[5];
1855 result.values_[4] = values_[1] * matrix.values_[3] + values_[4] * matrix.values_[4] + values_[7] * matrix.values_[5];
1856 result.values_[5] = values_[2] * matrix.values_[3] + values_[5] * matrix.values_[4] + values_[8] * matrix.values_[5];
1857 result.values_[6] = values_[0] * matrix.values_[6] + values_[3] * matrix.values_[7] + values_[6] * matrix.values_[8];
1858 result.values_[7] = values_[1] * matrix.values_[6] + values_[4] * matrix.values_[7] + values_[7] * matrix.values_[8];
1859 result.values_[8] = values_[2] * matrix.values_[6] + values_[5] * matrix.values_[7] + values_[8] * matrix.values_[8];
1860
1861 return result;
1862}
1863
1864template <typename T>
1865OCEAN_FORCE_INLINE VectorT2<T> SquareMatrixT3<T>::operator*(const VectorT2<T>& vector) const
1866{
1867 /**
1868 * | x' | | 0 3 6 | | x |
1869 * | y' | = | 1 4 7 | * | y |
1870 * | 1 | | 2 5 8 | | 1 |
1871 */
1872
1873 const T z = values_[2] * vector[0] + values_[5] * vector[1] + values_[8];
1874 ocean_assert(NumericT<T>::isNotEqualEps(z) && "Division by zero!");
1875
1876 const T factor = T(1) / z;
1877
1878 return VectorT2<T>((values_[0] * vector[0] + values_[3] * vector[1] + values_[6]) * factor,
1879 (values_[1] * vector[0] + values_[4] * vector[1] + values_[7]) * factor);
1880}
1881
1882template <typename T>
1883OCEAN_FORCE_INLINE VectorT3<T> SquareMatrixT3<T>::operator*(const VectorT3<T>& vector) const
1884{
1885 return VectorT3<T>(values_[0] * vector[0] + values_[3] * vector[1] + values_[6] * vector[2],
1886 values_[1] * vector[0] + values_[4] * vector[1] + values_[7] * vector[2],
1887 values_[2] * vector[0] + values_[5] * vector[1] + values_[8] * vector[2]);
1888}
1889
1890#if defined(OCEAN_HARDWARE_SSE_VERSION) && OCEAN_HARDWARE_SSE_VERSION >= 20
1891
1892template <>
1894{
1895 // the following code uses the following SSE instructions, and needs SSE2 or higher
1896
1897 // SSE2:
1898 // _mm_load1_pd
1899 // _mm_loadu_pd
1900 // _mm_mul_pd
1901 // _mm_add_pd
1902 // _mm_storeu_pd
1903
1904
1905 // we use the following strategy:
1906 // the values of the matrix are column aligned so that we normally would need to transpose the matrix before we can apply simple SIMD instructions
1907 // however, we do not transpose the matrix (we avoid the shuffle instructions) and instead multiply the matrix column-wise:
1908 // finally we sum the four columns and have the result, compared to the transpose-based approach this approach is approx. two times faster
1909 //
1910 // A D G a Aa + Db + Gc
1911 // B E H b Ba + Eb + Hc
1912 // C F I * c = Ca + Fb + Ic
1913
1914 const double* const vectorValues = vector.data();
1915
1916 // first we load the first vector element in all 64bit elements of the 128 bit register, so that we receive [a, a]
1917 const __m128d v0 = _mm_load1_pd(vectorValues + 0);
1918
1919 const __m128d c0 = _mm_loadu_pd(values_ + 0);
1920
1921 // now we multiply the 128 bit register [A, B] * [a, a] = [Aa, Ba]
1922 const __m128d r0 = _mm_mul_pd(c0, v0);
1923
1924
1925 // now we proceed with the second column
1926 const __m128d v1 = _mm_load1_pd(vectorValues + 1);
1927 const __m128d c1 = _mm_loadu_pd(values_ + 3);
1928 const __m128d r1 = _mm_mul_pd(c1, v1);
1929
1930 // and we sum the result of the first column with the result of the second column
1931 __m128d result_f64x2 = _mm_add_pd(r0, r1);
1932
1933
1934 // now we proceed with the third column
1935 const __m128d v2 = _mm_load1_pd(vectorValues + 2);
1936 const __m128d c2 = _mm_loadu_pd(values_ + 6);
1937 const __m128d r2 = _mm_mul_pd(c2, v2);
1938
1939 // we sum the results
1940 result_f64x2 = _mm_add_pd(result_f64x2, r2);
1941
1942 VectorT3<double> result;
1943 result[2] = values_[2] * vector[0] + values_[5] * vector[1] + values_[8] * vector[2];
1944 _mm_storeu_pd(result.data(), result_f64x2);
1945
1946 return result;
1947}
1948
1949template <>
1951{
1952 // the following code uses the following SSE instructions, and needs SSE2 or higher
1953
1954 // SSE1:
1955 // _mm_load1_ps
1956 // _mm_loadu_ps
1957 // _mm_mul_ps
1958 // _mm_add_ps
1959 // _mm_storeu_ps
1960
1961 // SSE2:
1962 // _mm_castps_si128
1963 // _mm_srli_si128
1964 // _mm_castsi128_ps
1965
1966
1967 // we use the following strategy:
1968 // the values of the matrix are column aligned so that we normally would need to transpose the matrix before we can apply simple SIMD instructions
1969 // however, we do not transpose the matrix (we avoid the shuffle instructions) and instead multiply the matrix column-wise:
1970 // finally we sum the four columns and have the result, compared to the transpose-based approach this approach is approx. two times faster
1971 //
1972 // A D G a Aa + Db + Gc
1973 // B E H b Ba + Eb + Hc
1974 // C F I * c = Ca + Fb + Ic
1975
1976
1977 // first we load the first vector element in all 32bit elements of the 128 bit register, so that we receive [a, a, a, a]
1978 const __m128 v0 = _mm_load1_ps(vector.data() + 0);
1979
1980 // now we load the first column to receive: [A, B, C, D], beware: D will be unused
1981 const __m128 c0 = _mm_loadu_ps(values_ + 0);
1982
1983 // now we multiply the 128 bit register [A, B, C, D] * [a, a, a, a] = [Aa, Ba, Ca, Da]
1984 const __m128 r0 = _mm_mul_ps(c0, v0);
1985
1986
1987 // now we proceed with the second column
1988 const __m128 v1 = _mm_load1_ps(vector.data() + 1);
1989 const __m128 c1 = _mm_loadu_ps(values_ + 3);
1990 const __m128 r1 = _mm_mul_ps(c1, v1);
1991
1992 // and we sum the result of the first column with the result of the second column
1993 __m128 result_f32x4 = _mm_add_ps(r0, r1);
1994
1995
1996 // now we proceed with the third column, this time we have to load [F, G, H, I] and we shift the values
1997 const __m128 v2 = _mm_load1_ps(vector.data() + 2);
1998 __m128 c2 = _mm_loadu_ps(values_ + 5);
1999 c2 = _mm_castsi128_ps(_mm_srli_si128(_mm_castps_si128(c2), 4)); // shift the 128 bit register by 4 bytes to the right to receive [G H I 0]
2000 const __m128 r2 = _mm_mul_ps(c2, v2);
2001
2002 // we sum the results
2003 result_f32x4 = _mm_add_ps(result_f32x4, r2);
2004
2005 float resultValues[4];
2006 _mm_storeu_ps(resultValues, result_f32x4);
2007
2008 // and finally we store the results back to the vector
2009 return VectorT3<float>(resultValues);
2010}
2011
2012#endif // OCEAN_HARDWARE_SSE_VERSION >= 20
2013
2014template <typename T>
2016{
2017 SquareMatrixT3<T> result(*this);
2018 result.values_[0] *= value;
2019 result.values_[1] *= value;
2020 result.values_[2] *= value;
2021 result.values_[3] *= value;
2022 result.values_[4] *= value;
2023 result.values_[5] *= value;
2024 result.values_[6] *= value;
2025 result.values_[7] *= value;
2026 result.values_[8] *= value;
2027
2028 return result;
2029}
2030
2031template <typename T>
2033{
2034 values_[0] *= value;
2035 values_[1] *= value;
2036 values_[2] *= value;
2037 values_[3] *= value;
2038 values_[4] *= value;
2039 values_[5] *= value;
2040 values_[6] *= value;
2041 values_[7] *= value;
2042 values_[8] *= value;
2043
2044 return *this;
2045}
2046
2047template <typename T>
2048void SquareMatrixT3<T>::multiply(const SquareMatrixT3<T>& matrix, const VectorT2<T>* vectors, VectorT2<T>* results, const size_t number)
2049{
2050 ocean_assert((vectors && results) || number == 0);
2051
2052 for (size_t n = 0; n < number; ++n)
2053 {
2054 results[n] = matrix * vectors[n];
2055 }
2056}
2057
2058template <typename T>
2059void SquareMatrixT3<T>::multiply(const SquareMatrixT3<T>& matrix, const VectorT3<T>* vectors, VectorT3<T>* results, const size_t number)
2060{
2061 ocean_assert((vectors && results) || number == 0);
2062
2063 for (size_t n = 0; n < number; ++n)
2064 {
2065 results[n] = matrix * vectors[n];
2066 }
2067}
2068
2069#if defined(OCEAN_HARDWARE_SSE_VERSION) && OCEAN_HARDWARE_SSE_VERSION >= 20
2070
2071template <>
2072inline void SquareMatrixT3<float>::multiply(const SquareMatrixT3<float>& matrix, const VectorT3<float>* vectors, VectorT3<float>* results, const size_t number)
2073{
2074 // the following code uses the following SSE instructions, and needs SSE2 or higher
2075
2076 // SSE1:
2077 // _mm_load1_ps
2078 // _mm_loadu_ps
2079 // _mm_mul_ps
2080 // _mm_add_ps
2081 // _mm_storeu_ps
2082
2083 // SSE2:
2084 // _mm_castps_si128
2085 // _mm_srli_si128
2086 // _mm_castsi128_ps
2087
2088 // now we load the three columns (and ignore the last entry)
2089 const __m128 c0 = _mm_loadu_ps(matrix.values_ + 0);
2090 const __m128 c1 = _mm_loadu_ps(matrix.values_ + 3);
2091 __m128 c2 = _mm_loadu_ps(matrix.values_ + 5);
2092 c2 = _mm_castsi128_ps(_mm_srli_si128(_mm_castps_si128(c2), 4)); // shift the 128 bit register by 4 bytes to the right to receive [G H I 0]
2093
2094 for (size_t n = 0u; n < number - 1; ++n)
2095 {
2096 const __m128 v0 = _mm_load1_ps(vectors[n].data() + 0);
2097 const __m128 v1 = _mm_load1_ps(vectors[n].data() + 1);
2098 const __m128 v2 = _mm_load1_ps(vectors[n].data() + 2);
2099
2100 const __m128 r0 = _mm_mul_ps(c0, v0);
2101 const __m128 r1 = _mm_mul_ps(c1, v1);
2102
2103 __m128 result_f32x4= _mm_add_ps(r0, r1);
2104
2105 const __m128 r2 = _mm_mul_ps(c2, v2);
2106
2107 result_f32x4 = _mm_add_ps(result_f32x4, r2);
2108
2109 // now we write all four values (although the target vector holds 3 elements - we can due that as there is still one 3D vector target vector left)
2110 _mm_storeu_ps(results[n].data(), result_f32x4);
2111 }
2112
2113 // now we handle the last vector explicitly
2114
2115 const __m128 v0 = _mm_load1_ps(vectors[number - 1].data() + 0);
2116 const __m128 v1 = _mm_load1_ps(vectors[number - 1].data() + 1);
2117 const __m128 v2 = _mm_load1_ps(vectors[number - 1].data() + 2);
2118
2119 const __m128 r0 = _mm_mul_ps(c0, v0);
2120 const __m128 r1 = _mm_mul_ps(c1, v1);
2121
2122 __m128 result_f32x4 = _mm_add_ps(r0, r1);
2123
2124 const __m128 r2 = _mm_mul_ps(c2, v2);
2125
2126 result_f32x4 = _mm_add_ps(result_f32x4, r2);
2127
2128 float resultValues[4];
2129 _mm_storeu_ps(resultValues, result_f32x4);
2130
2131 results[number - 1] = VectorT3<float>(resultValues);
2132}
2133
2134#endif // OCEAN_HARDWARE_SSE_VERSION >= 20
2135
2136template <typename T>
2137template <typename U>
2139{
2140 SquareMatricesT3<T> result;
2141 result.reserve(matrices.size());
2142
2143 for (const SquareMatrixT3<U>& matrix : matrices)
2144 {
2145 result.emplace_back(matrix);
2146 }
2147
2148 return result;
2149}
2150
2151template <>
2152template <>
2153inline std::vector< SquareMatrixT3<float> > SquareMatrixT3<float>::matrices2matrices(const std::vector< SquareMatrixT3<float> >& matrices)
2154{
2155 return matrices;
2156}
2157
2158template <>
2159template <>
2160inline std::vector< SquareMatrixT3<double> > SquareMatrixT3<double>::matrices2matrices(const std::vector< SquareMatrixT3<double> >& matrices)
2161{
2162 return matrices;
2163}
2164
2165template <typename T>
2166template <typename U>
2167inline std::vector< SquareMatrixT3<T> > SquareMatrixT3<T>::matrices2matrices(const SquareMatrixT3<U>* matrices, const size_t size)
2168{
2169 std::vector< SquareMatrixT3<T> > result;
2170 result.reserve(size);
2171
2172 for (size_t n = 0; n < size; ++n)
2173 {
2174 result.emplace_back(matrices[n]);
2175 }
2176
2177 return result;
2178}
2179
2180template <typename T>
2181std::ostream& operator<<(std::ostream& stream, const SquareMatrixT3<T>& matrix)
2182{
2183 stream << "|" << matrix(0, 0) << ", " << matrix(0, 1) << ", " << matrix(0, 2) << "|" << std::endl;
2184 stream << "|" << matrix(1, 0) << ", " << matrix(1, 1) << ", " << matrix(1, 2) << "|" << std::endl;
2185 stream << "|" << matrix(2, 0) << ", " << matrix(2, 1) << ", " << matrix(2, 2) << "|";
2186
2187 return stream;
2188}
2189
2190template <bool tActive, typename T>
2191MessageObject<tActive>& operator<<(MessageObject<tActive>& messageObject, const SquareMatrixT3<T>& matrix)
2192{
2193 return messageObject << "|" << matrix(0, 0) << ", " << matrix(0, 1) << ", " << matrix(0, 2) << "|\n|"
2194 << matrix(1, 0) << ", " << matrix(1, 1) << ", " << matrix(1, 2) << "|\n|"
2195 << matrix(2, 0) << ", " << matrix(2, 1) << ", " << matrix(2, 2) << "|";
2196}
2197
2198template <bool tActive, typename T>
2199MessageObject<tActive>& operator<<(MessageObject<tActive>&& messageObject, const SquareMatrixT3<T>& matrix)
2200{
2201 return messageObject << "|" << matrix(0, 0) << ", " << matrix(0, 1) << ", " << matrix(0, 2) << "|\n|"
2202 << matrix(1, 0) << ", " << matrix(1, 1) << ", " << matrix(1, 2) << "|\n|"
2203 << matrix(2, 0) << ", " << matrix(2, 1) << ", " << matrix(2, 2) << "|";
2204}
2205
2206}
2207
2208#endif // META_OCEAN_MATH_SQUARE_MATRIX_3_H
This class provides several functions to solve equations with different degree using floating point v...
Definition Equation.h:53
This class implements an euler rotation with angles: yaw, pitch and roll.
Definition Euler.h:80
const T & yaw() const
Returns the yaw angle.
Definition Euler.h:315
const T & roll() const
Returns the roll angle.
Definition Euler.h:339
const T & pitch() const
Returns the pitch angle.
Definition Euler.h:327
This class implements a 4x4 homogeneous transformation matrix using floating point values with the pr...
Definition HomogenousMatrix4.h:110
This class provides basic numeric functionalities.
Definition Numeric.h:57
static T sin(const T value)
Returns the sine of a given value.
Definition Numeric.h:1568
static T abs(const T value)
Returns the absolute value of a given value.
Definition Numeric.h:1220
static T sqrt(const T value)
Returns the square root of a given value.
Definition Numeric.h:1533
static bool isEqual(const T first, const T second)
Returns whether two values are equal up to a small epsilon.
Definition Numeric.h:2389
static constexpr T sqr(const T value)
Returns the square of a given value.
Definition Numeric.h:1495
static constexpr bool isEqualEps(const T value)
Returns whether a value is smaller than or equal to a small epsilon.
Definition Numeric.h:2090
static T cos(const T value)
Returns the cosine of a given value.
Definition Numeric.h:1584
static constexpr bool isNotEqualEps(const T value)
Returns whether a value is not smaller than or equal to a small epsilon.
Definition Numeric.h:2240
This class implements a unit quaternion rotation.
Definition Quaternion.h:100
const T & x() const
Returns the x value of the quaternion.
Definition Quaternion.h:919
bool isValid() const
Returns whether this quaternion is a valid unit quaternion.
Definition Quaternion.h:901
const T & w() const
Returns the w value of the quaternion.
Definition Quaternion.h:907
const T & y() const
Returns the y value of the quaternion.
Definition Quaternion.h:931
const T & z() const
Returns the z value of the quaternion.
Definition Quaternion.h:943
This class implements a axis-angle rotation using floating point values.
Definition Rotation.h:79
VectorT3< T > axis() const
Returns the axis of the rotation.
Definition Rotation.h:748
T angle() const
Returns the angle of the rotation.
Definition Rotation.h:754
bool isValid() const
Returns whether this rotation has valid parameters.
Definition Rotation.h:673
size_t operator()(const SquareMatrixT3< T > &matrix) const
Hash function.
Definition SquareMatrix3.h:1126
VectorT3< T > xAxis() const
Returns the x axis which is the first column of the matrix.
Definition SquareMatrix3.h:1397
bool eigenValues(T *eigenValues) const
Determines the eigen values of this matrix.
Definition SquareMatrix3.h:1468
T * operator()()
Access operator.
Definition SquareMatrix3.h:1120
SquareMatrixT3< T > orthonormalMatrix() const
Returns the orthonormal matrix of this matrix by scaling the x-axis and adjusting y- and z-axis.
Definition SquareMatrix3.h:1415
SquareMatrixT3(const SquareMatrixT3< U > &matrix)
Copy constructor for a matrix with difference element data type than T.
Definition SquareMatrix3.h:707
OCEAN_FORCE_INLINE VectorT2< T > operator*(const VectorT2< T > &vector) const
Multiply operator for a 2D vector.
Definition SquareMatrix3.h:1865
static SquareMatrixT3< T > skewSymmetricMatrix(const VectorT3< T > &vector)
Creates a skew symmetric 3x3 matrix for a specific vector.
Definition SquareMatrix3.h:1702
T operator[](const unsigned int index) const
Element operator.
Definition SquareMatrix3.h:1072
SquareMatrixT3(const U *arrayValues)
Creates a new SquareMatrixT3 object by nine elements of float type U.
Definition SquareMatrix3.h:912
bool eigenSystem(T *eigenValues, VectorT3< T > *eigenVectors) const
Performs an eigen value analysis.
Definition SquareMatrix3.h:1521
bool isEqual(const SquareMatrixT3< T > &matrix, const T eps=NumericT< T >::eps()) const
Returns whether two matrices are almost identical up to a specified epsilon.
Definition SquareMatrix3.h:1387
bool invert()
Inverts this matrix in place.
Definition SquareMatrix3.h:1191
T & operator()(const unsigned int index)
Element operator.
Definition SquareMatrix3.h:1107
bool isNull() const
Returns whether this matrix is a zero matrix.
Definition SquareMatrix3.h:1334
static size_t elements()
Returns the number of elements this matrix has.
Definition SquareMatrix3.h:1139
SquareMatrixT3()
Creates a new SquareMatrixT3 object with undefined elements.
Definition SquareMatrix3.h:700
bool isAffine() const
Returns true if this matrix is a affine transformation.
Definition SquareMatrix3.h:1354
SquareMatrixT3< T > operator+(const SquareMatrixT3< T > &matrix) const
Adds two matrices.
Definition SquareMatrix3.h:1761
SquareMatrixT3(const VectorT3< T > &xAxis, const VectorT3< T > &yAxis, const VectorT3< T > &zAxis)
Creates a new SquareMatrixT3 object by three axes.
Definition SquareMatrix3.h:889
T values_[9]
The nine values of the matrix.
Definition SquareMatrix3.h:696
SquareMatrixT3(const U *arrayValues, const bool valuesRowAligned)
Creates a new SquareMatrixT3 object by nine elements.
Definition SquareMatrix3.h:931
SquareMatrixT3< T > transposedMultiply(const SquareMatrixT3< T > &right) const
Multiplies this transposed matrix with a second matrix.
Definition SquareMatrix3.h:1732
T Type
Definition of the used data type.
Definition SquareMatrix3.h:97
SquareMatrixT3(const T *arrayValues, const bool valuesRowAligned)
Creates a new SquareMatrixT3 object by nine elements.
Definition SquareMatrix3.h:958
void transpose()
Transposes the matrix.
Definition SquareMatrix3.h:1162
SquareMatrixT3< T > & operator=(const SquareMatrixT3< T > &)=default
Default copy assignment operator.
OCEAN_FORCE_INLINE SquareMatrixT3< T > operator*(const T value) const
Multiplies this matrix with a scalar value.
Definition SquareMatrix3.h:2015
OCEAN_FORCE_INLINE VectorT3< T > operator*(const VectorT3< T > &vector) const
Multiply operator for a 3D vector.
Definition SquareMatrix3.h:1883
SquareMatrixT3< T > transposed() const
Returns the transposed of this matrix.
Definition SquareMatrix3.h:1145
static SquareMatricesT3< T > matrices2matrices(const SquareMatricesT3< U > &matrices)
Converts matrices with specific data type to matrices with different data type.
Definition SquareMatrix3.h:2138
void copyElements(U *arrayValues, const bool valuesRowAligned=false) const
Copies the elements of this matrix to an array with floating point values of type U.
Definition SquareMatrix3.h:1632
bool multiply(const VectorT2< T > &vector, VectorT2< T > &result) const
Multiplies a 2D vector with this matrix (from the right).
Definition SquareMatrix3.h:1708
SquareMatrixT3(const SquareMatrixT3< T > &matrix)=default
Copy constructor.
T absSum() const
Calculates the sum of absolute value of elements.
Definition SquareMatrix3.h:1035
T & operator[](const unsigned int index)
Element operator.
Definition SquareMatrix3.h:1079
SquareMatrixT3< T > operator-() const
Returns the negative matrix of this matrix (all matrix elements are multiplied by -1).
Definition SquareMatrix3.h:1829
friend class SquareMatrixT3
Definition SquareMatrix3.h:90
OCEAN_FORCE_INLINE SquareMatrixT3< T > operator*(const SquareMatrixT3< T > &matrix) const
Multiplies two matrices.
Definition SquareMatrix3.h:1847
SquareMatrixT3< T > & operator-=(const SquareMatrixT3< T > &matrix)
Subtracts and assigns two matrices.
Definition SquareMatrix3.h:1813
static void multiply(const SquareMatrixT3< T > &matrix, const VectorT3< T > *vectors, VectorT3< T > *results, const size_t number)
Multiplies several 3D vectors with a given 3x3 matrix.
Definition SquareMatrix3.h:2059
bool operator!=(const SquareMatrixT3< T > &matrix) const
Returns whether two matrices are not identical up to a small epsilon.
Definition SquareMatrix3.h:1059
static void multiply(const SquareMatrixT3< T > &matrix, const VectorT2< T > *vectors, VectorT2< T > *results, const size_t number)
Multiplies several 2D vectors with a given 3x3 matrix.
Definition SquareMatrix3.h:2048
SquareMatrixT3(const bool setToIdentity)
Creates a new SquareMatrixT3 object.
Definition SquareMatrix3.h:721
SquareMatrixT3(const EulerT< T > &euler)
Creates a new SquareMatrixT3 rotation matrix by a given Euler rotation.
Definition SquareMatrix3.h:750
SquareMatrixT3(const HomogenousMatrixT4< T > &transformation)
Creates a 3x3 rotation matrix by a given 4x4 homogeneous transformation.
Definition SquareMatrix3.h:981
VectorT3< T > zAxis() const
Returns the z axis which is the last column of the matrix.
Definition SquareMatrix3.h:1409
SquareMatrixT3(const RotationT< T > &rotation)
Creates a new 3x3 matrix object by a given angle-axis rotation.
Definition SquareMatrix3.h:803
void toNull()
Sets the matrix to a zero matrix.
Definition SquareMatrix3.h:1312
const T * data() const
Returns a pointer to the internal values.
Definition SquareMatrix3.h:1047
T * data()
Returns a pointer to the internal values.
Definition SquareMatrix3.h:1053
const T * operator()() const
Access operator.
Definition SquareMatrix3.h:1114
void copyElements(T *arrayValues, const bool valuesRowAligned=false) const
Copies the elements of this matrix to an array with floating point values of type T.
Definition SquareMatrix3.h:1671
bool isHomography() const
Returns true if this matrix is perspective transform/homography.
Definition SquareMatrix3.h:1360
bool isSimilarity() const
Returns true if this matrix is a similarity transformation.
Definition SquareMatrix3.h:1348
VectorT3< T > yAxis() const
Returns the y axis which is the middle column of the matrix.
Definition SquareMatrix3.h:1403
void toIdentity()
Sets the matrix to the identity matrix.
Definition SquareMatrix3.h:1298
bool solve(const VectorT3< T > &b, VectorT3< T > &x) const
Solve a simple 3x3 system of linear equations: Ax = b Beware: The system of linear equations is assum...
Definition SquareMatrix3.h:997
SquareMatrixT3(const T &m00, const T &m10, const T &m20, const T &m01, const T &m11, const T &m21, const T &m02, const T &m12, const T &m22)
Creates a 3x3 rotation matrix by 9 given matrix elements.
Definition SquareMatrix3.h:873
T operator()(const unsigned int index) const
Element operator.
Definition SquareMatrix3.h:1100
T determinant() const
Returns the determinant of the matrix.
Definition SquareMatrix3.h:1284
SquareMatrixT3< T > operator-(const SquareMatrixT3< T > &matrix) const
Subtracts two matrices.
Definition SquareMatrix3.h:1795
SquareMatrixT3< T > & operator+=(const SquareMatrixT3< T > &matrix)
Adds and assigns two matrices.
Definition SquareMatrix3.h:1779
SquareMatrixT3(const T *arrayValues)
Creates a new SquareMatrixT3 object by nine elements.
Definition SquareMatrix3.h:923
bool invert(SquareMatrixT3< T > &invertedMatrix) const
Inverts the matrix and returns the result as parameter.
Definition SquareMatrix3.h:1206
SquareMatrixT3(const SquareMatrixT4< T > &transformation)
Creates a 3x3 square matrix by a given 4x4 square transformation.
Definition SquareMatrix3.h:989
T sum() const
Calculates the sum of elements.
Definition SquareMatrix3.h:1041
SquareMatrixT3< T > & operator*=(const T value)
Multiplies and assigns this matrix with a scalar value.
Definition SquareMatrix3.h:2032
bool isIdentity() const
Returns whether this matrix is an identity matrix.
Definition SquareMatrix3.h:1326
VectorT3< T > diagonal() const
Returns a 3d vector with values of the matrix diagonal.
Definition SquareMatrix3.h:1625
SquareMatrixT3(const VectorT3< T > &diagonal)
Creates a new SquareMatrixT3 object by a given diagonal vector.
Definition SquareMatrix3.h:897
OCEAN_FORCE_INLINE SquareMatrixT3< T > & operator*=(const SquareMatrixT3< T > &matrix)
Multiplies and assigns two matrices.
Definition SquareMatrix3.h:1065
static SquareMatricesT3< T > matrices2matrices(const SquareMatrixT3< U > *matrices, const size_t size)
Converts matrices with specific data type to matrices with different data type.
SquareMatrixT3< T > inverted() const
Returns the inverted matrix of this matrix.
Definition SquareMatrix3.h:1177
bool isOrthonormal(const T epsilon=NumericT< T >::eps()) const
Returns whether this matrix is an orthonormal matrix.
Definition SquareMatrix3.h:1366
bool isSymmetric(const T epsilon=NumericT< T >::eps()) const
Returns whether this matrix is symmetric.
Definition SquareMatrix3.h:1379
bool operator==(const SquareMatrixT3< T > &matrix) const
Returns whether two matrices are identical up to a small epsilon.
Definition SquareMatrix3.h:1755
T & operator()(const unsigned int row, const unsigned int column)
Element operator.
Definition SquareMatrix3.h:1093
bool isSingular() const
Returns whether this matrix is singular (and thus cannot be inverted).
Definition SquareMatrix3.h:1342
T trace() const
Returns the trace of the matrix which is the sum of the diagonal elements.
Definition SquareMatrix3.h:1292
SquareMatrixT3(const QuaternionT< T > &quaternion)
Creates a new 3x3 matrix object by a given quaternion rotation.
Definition SquareMatrix3.h:842
T operator()(const unsigned int row, const unsigned int column) const
Element operator.
Definition SquareMatrix3.h:1086
This class implements a 4x4 square matrix.
Definition SquareMatrix4.h:85
static void sortHighestToFront3(T &value0, T &value1, T &value2)
Sorts three values so that the highest value will finally be the first value.
Definition base/Utilities.h:724
This class implements a vector with two elements.
Definition Vector2.h:96
This class implements a vector with three elements.
Definition Vector3.h:97
VectorT3< T > perpendicular() const
Returns a vector that is perpendicular to this vector.
Definition Vector3.h:768
const T & y() const noexcept
Returns the y value.
Definition Vector3.h:824
bool normalize()
Normalizes this vector.
Definition Vector3.h:659
VectorT3< T > cross(const VectorT3< T > &vector) const
Returns the cross product of two vectors.
Definition Vector3.h:609
const T & x() const noexcept
Returns the x value.
Definition Vector3.h:812
const T & z() const noexcept
Returns the z value.
Definition Vector3.h:836
const T * data() const noexcept
Returns an pointer to the vector elements.
Definition Vector3.h:854
T length() const
Returns the length of the vector.
Definition Vector3.h:676
T sqr() const
Returns the square of the vector length.
Definition Vector3.h:682
unsigned int sqrDistance(const char first, const char second)
Returns the square distance between two values.
Definition base/Utilities.h:1113
std::vector< SquareMatrixT3< T > > SquareMatricesT3
Definition of a typename alias for vectors with SquareMatrixT3 objects.
Definition SquareMatrix3.h:65
std::vector< SquareMatrix3 > SquareMatrices3
Definition of a vector holding SquareMatrix3 objects.
Definition SquareMatrix3.h:72
The namespace covering the entire Ocean framework.
Definition Accessor.h:15
std::ostream & operator<<(std::ostream &stream, const HighPerformanceStatistic &highPerformanceStatistic)
Definition HighPerformanceTimer.h:963