Blender V5.0
BLI_math_matrix.hh
Go to the documentation of this file.
1/* SPDX-FileCopyrightText: 2022 Blender Authors
2 *
3 * SPDX-License-Identifier: GPL-2.0-or-later */
4
5#pragma once
6
10
11#include "BLI_math_base.hh"
14#include "BLI_math_vector.hh"
15#include "BLI_unroll.hh"
16
17namespace blender::math {
18
19/* -------------------------------------------------------------------- */
22
27template<typename T, int Size>
28[[nodiscard]] MatBase<T, Size, Size> invert(const MatBase<T, Size, Size> &mat, bool &r_success);
29
33template<typename T, int NumCol, int NumRow>
34[[nodiscard]] MatBase<T, NumCol, NumRow> transpose(const MatBase<T, NumRow, NumCol> &mat);
35
39template<typename T, int NumCol, int NumRow>
40[[nodiscard]] MatBase<T, NumCol, NumRow> normalize(const MatBase<T, NumCol, NumRow> &a);
41
46template<typename T, int NumCol, int NumRow, typename VectorT>
47[[nodiscard]] MatBase<T, NumCol, NumRow> normalize_and_get_size(
48 const MatBase<T, NumCol, NumRow> &a, VectorT &r_size);
49
54template<typename T, int Size> [[nodiscard]] T determinant(const MatBase<T, Size, Size> &mat);
55
59template<typename T, int Size>
60[[nodiscard]] MatBase<T, Size, Size> adjoint(const MatBase<T, Size, Size> &mat);
61
65template<typename T, int NumCol, int NumRow, typename VectorT>
66[[nodiscard]] MatBase<T, NumCol, NumRow> translate(const MatBase<T, NumCol, NumRow> &mat,
67 const VectorT &translation);
68
73template<typename T, int NumCol, int NumRow, typename RotationT>
75 const RotationT &rotation);
76
80template<typename T, int NumCol, int NumRow, typename VectorT>
82 const VectorT &scale);
83
87template<typename T, int NumCol, int NumRow>
90 T t);
91
107template<typename T>
108[[nodiscard]] MatBase<T, 3, 3> interpolate(const MatBase<T, 3, 3> &a,
109 const MatBase<T, 3, 3> &b,
110 T t);
111
120template<typename T>
121[[nodiscard]] MatBase<T, 4, 4> interpolate(const MatBase<T, 4, 4> &a,
122 const MatBase<T, 4, 4> &b,
123 T t);
124
136template<typename T>
138 const MatBase<T, 3, 3> &b,
139 T t);
140
153template<typename T>
155 const MatBase<T, 4, 4> &b,
156 T t);
157
168template<typename T, int Size>
170 T epsilon = 1e-8);
171
173
174/* -------------------------------------------------------------------- */
177
181template<typename MatT> [[nodiscard]] MatT from_location(const typename MatT::loc_type &location);
182
187template<typename MatT, int ScaleDim>
189
193template<typename MatT, typename RotationT>
194[[nodiscard]] MatT from_rotation(const RotationT &rotation);
195
199template<typename MatT, typename RotationT, typename VectorT>
200[[nodiscard]] MatT from_rot_scale(const RotationT &rotation, const VectorT &scale);
201
205template<typename MatT, typename RotationT>
206[[nodiscard]] MatT from_loc_rot(const typename MatT::loc_type &location,
207 const RotationT &rotation);
208
212template<typename MatT, int ScaleDim>
213[[nodiscard]] MatT from_loc_scale(const typename MatT::loc_type &location,
215
219template<typename MatT, typename RotationT, int ScaleDim>
220[[nodiscard]] MatT from_loc_rot_scale(const typename MatT::loc_type &location,
221 const RotationT &rotation,
223
228template<typename T> [[nodiscard]] MatBase<T, 2, 2> from_direction(const VecBase<T, 2> &direction);
229
235template<typename MatT, typename VectorT>
236[[nodiscard]] MatT from_orthonormal_axes(const VectorT forward, const VectorT up);
237
242template<typename MatT, typename VectorT>
243[[nodiscard]] MatT from_orthonormal_axes(const VectorT location,
244 const VectorT forward,
245 const VectorT up);
246
255template<typename MatT, typename VectorT> [[nodiscard]] MatT from_up_axis(const VectorT up);
256
266template<typename MatT> [[nodiscard]] MatT orthogonalize(const MatT &mat, const Axis axis);
267
273template<typename MatT, typename VectorT>
274[[nodiscard]] MatT from_origin_transform(const MatT &transform, const VectorT origin);
275
277
278/* -------------------------------------------------------------------- */
281
286template<typename T> [[nodiscard]] inline AngleRadianBase<T> to_angle(const MatBase<T, 2, 2> &mat);
287template<typename T> [[nodiscard]] inline EulerXYZBase<T> to_euler(const MatBase<T, 3, 3> &mat);
288template<typename T> [[nodiscard]] inline EulerXYZBase<T> to_euler(const MatBase<T, 4, 4> &mat);
289template<typename T>
290[[nodiscard]] inline Euler3Base<T> to_euler(const MatBase<T, 3, 3> &mat, EulerOrder order);
291template<typename T>
292[[nodiscard]] inline Euler3Base<T> to_euler(const MatBase<T, 4, 4> &mat, EulerOrder order);
293
301template<typename T>
302[[nodiscard]] inline EulerXYZBase<T> to_nearest_euler(const MatBase<T, 3, 3> &mat,
303 const EulerXYZBase<T> &reference);
304template<typename T>
305[[nodiscard]] inline EulerXYZBase<T> to_nearest_euler(const MatBase<T, 4, 4> &mat,
306 const EulerXYZBase<T> &reference);
307template<typename T>
308[[nodiscard]] inline Euler3Base<T> to_nearest_euler(const MatBase<T, 3, 3> &mat,
309 const Euler3Base<T> &reference);
310template<typename T>
311[[nodiscard]] inline Euler3Base<T> to_nearest_euler(const MatBase<T, 4, 4> &mat,
312 const Euler3Base<T> &reference);
313
317template<typename T>
318[[nodiscard]] inline QuaternionBase<T> to_quaternion(const MatBase<T, 3, 3> &mat);
319template<typename T>
320[[nodiscard]] inline QuaternionBase<T> to_quaternion(const MatBase<T, 4, 4> &mat);
321
328[[nodiscard]] Quaternion to_quaternion_legacy(const float3x3 &mat);
329
335template<bool AllowNegativeScale = false, typename T, int NumCol, int NumRow>
336[[nodiscard]] inline VecBase<T, 3> to_scale(const MatBase<T, NumCol, NumRow> &mat);
337template<bool AllowNegativeScale = false, typename T>
338[[nodiscard]] inline VecBase<T, 2> to_scale(const MatBase<T, 2, 2> &mat);
339
346template<bool AllowNegativeScale = false, typename T>
347inline void to_rot_scale(const MatBase<T, 2, 2> &mat,
348 AngleRadianBase<T> &r_rotation,
349 VecBase<T, 2> &r_scale);
350template<bool AllowNegativeScale = false, typename T>
351inline void to_loc_rot_scale(const MatBase<T, 3, 3> &mat,
352 VecBase<T, 2> &r_location,
353 AngleRadianBase<T> &r_rotation,
354 VecBase<T, 2> &r_scale);
355template<bool AllowNegativeScale = false, typename T, typename RotationT>
356inline void to_rot_scale(const MatBase<T, 3, 3> &mat,
357 RotationT &r_rotation,
358 VecBase<T, 3> &r_scale);
359template<bool AllowNegativeScale = false, typename T, typename RotationT>
360inline void to_loc_rot_scale(const MatBase<T, 4, 4> &mat,
361 VecBase<T, 3> &r_location,
362 RotationT &r_rotation,
363 VecBase<T, 3> &r_scale);
364
366
367/* -------------------------------------------------------------------- */
370
374template<typename T>
375[[nodiscard]] VecBase<T, 2> transform_point(const MatBase<T, 2, 2> &mat,
376 const VecBase<T, 2> &point);
377
381template<typename T>
382[[nodiscard]] VecBase<T, 2> transform_point(const MatBase<T, 3, 3> &mat,
383 const VecBase<T, 2> &point);
384
388template<typename T>
389[[nodiscard]] VecBase<T, 3> transform_point(const MatBase<T, 3, 3> &mat,
390 const VecBase<T, 3> &point);
391
395template<typename T>
396[[nodiscard]] VecBase<T, 3> transform_point(const MatBase<T, 4, 4> &mat,
397 const VecBase<T, 3> &point);
398
402template<typename T>
403[[nodiscard]] VecBase<T, 3> transform_direction(const MatBase<T, 3, 3> &mat,
404 const VecBase<T, 3> &direction);
405
409template<typename T>
410[[nodiscard]] VecBase<T, 3> transform_direction(const MatBase<T, 4, 4> &mat,
411 const VecBase<T, 3> &direction);
412
416template<typename MatT, typename VectorT>
417[[nodiscard]] VectorT project_point(const MatT &mat, const VectorT &point);
418
420
421/* -------------------------------------------------------------------- */
424
425namespace projection {
426
432template<typename T>
433[[nodiscard]] MatBase<T, 4, 4> orthographic(
434 T left, T right, T bottom, T top, T near_clip, T far_clip);
435
443template<typename T> MatBase<T, 4, 4> orthographic_infinite(T left, T right, T bottom, T top);
444
451template<typename T>
452[[nodiscard]] MatBase<T, 4, 4> perspective(
453 T left, T right, T bottom, T top, T near_clip, T far_clip);
454
461template<typename T>
463 T angle_left, T angle_right, T angle_bottom, T angle_top, T near_clip, T far_clip);
464
472template<typename T>
473[[nodiscard]] MatBase<T, 4, 4> perspective_infinite(T left, T right, T bottom, T top, T near_clip);
474
480template<typename T>
481[[nodiscard]] MatBase<T, 4, 4> translate(const MatBase<T, 4, 4> &mat, const VecBase<T, 2> &offset);
482
483} // namespace projection
484
486
487/* -------------------------------------------------------------------- */
490
497template<typename T> [[nodiscard]] bool is_negative(const MatBase<T, 3, 3> &mat);
498template<typename T> [[nodiscard]] bool is_negative(const MatBase<T, 4, 4> &mat);
499
503template<typename T, int NumCol, int NumRow>
504[[nodiscard]] inline bool is_equal(const MatBase<T, NumCol, NumRow> &a,
506 const T epsilon = T(0))
507{
508 for (int i = 0; i < NumCol; i++) {
509 for (int j = 0; j < NumRow; j++) {
510 if (math::abs(a[i][j] - b[i][j]) > epsilon) {
511 return false;
512 }
513 }
514 }
515 return true;
516}
517
521template<typename T, int NumCol, int NumRow>
522[[nodiscard]] inline bool is_identity(const MatBase<T, NumCol, NumRow> &mat)
523{
524 for (int i = 0; i < NumCol; i++) {
525 for (int j = 0; j < NumRow; j++) {
526 if (mat[i][j] != (i != j ? 0.0f : 1.0f)) {
527 return false;
528 }
529 }
530 }
531 return true;
532}
533
537template<typename MatT> [[nodiscard]] inline bool is_orthogonal(const MatT &mat)
538{
539 if (math::abs(math::dot(mat.x_axis(), mat.y_axis())) > 1e-5f) {
540 return false;
541 }
542 if (math::abs(math::dot(mat.y_axis(), mat.z_axis())) > 1e-5f) {
543 return false;
544 }
545 if (math::abs(math::dot(mat.z_axis(), mat.x_axis())) > 1e-5f) {
546 return false;
547 }
548 return true;
549}
550
554template<typename MatT> [[nodiscard]] inline bool is_orthonormal(const MatT &mat)
555{
556 if (!is_orthogonal(mat)) {
557 return false;
558 }
559 if (math::abs(math::length_squared(mat.x_axis()) - 1) > 1e-5f) {
560 return false;
561 }
562 if (math::abs(math::length_squared(mat.y_axis()) - 1) > 1e-5f) {
563 return false;
564 }
565 if (math::abs(math::length_squared(mat.z_axis()) - 1) > 1e-5f) {
566 return false;
567 }
568 return true;
569}
570
574template<typename MatT> [[nodiscard]] inline bool is_uniformly_scaled(const MatT &mat)
575{
576 if (!is_orthogonal(mat)) {
577 return false;
578 }
579 using T = typename MatT::base_type;
580 const T eps = 1e-7;
581 const T x = math::length_squared(mat.x_axis());
582 const T y = math::length_squared(mat.y_axis());
583 const T z = math::length_squared(mat.z_axis());
584 return (math::abs(x - y) < eps) && math::abs(x - z) < eps;
585}
586
587template<typename T, int NumCol, int NumRow>
588inline bool is_zero(const MatBase<T, NumCol, NumRow> &mat)
589{
590 for (int i = 0; i < NumCol; i++) {
591 if (!is_zero(mat[i])) {
592 return false;
593 }
594 }
595 return true;
596}
597
599
600/* -------------------------------------------------------------------- */
603
604/* Implementation details. */
605namespace detail {
606
607template<typename T, int NumCol, int NumRow>
608[[nodiscard]] MatBase<T, NumCol, NumRow> from_rotation(const AngleRadianBase<T> &rotation);
609
610template<typename T, int NumCol, int NumRow>
611[[nodiscard]] MatBase<T, NumCol, NumRow> from_rotation(const EulerXYZBase<T> &rotation);
612
613template<typename T, int NumCol, int NumRow>
614[[nodiscard]] MatBase<T, NumCol, NumRow> from_rotation(const Euler3Base<T> &rotation);
615
616template<typename T, int NumCol, int NumRow>
617[[nodiscard]] MatBase<T, NumCol, NumRow> from_rotation(const QuaternionBase<T> &rotation);
618
619template<typename T, int NumCol, int NumRow>
621
622template<typename T, int NumCol, int NumRow>
623[[nodiscard]] MatBase<T, NumCol, NumRow> from_rotation(const CartesianBasis &rotation);
624
625template<typename T, int NumCol, int NumRow, typename AngleT>
627
628} // namespace detail
629
630/* Returns true if each individual columns are unit scaled. Mainly for assert usage. */
631template<typename T, int NumCol, int NumRow>
632[[nodiscard]] inline bool is_unit_scale(const MatBase<T, NumCol, NumRow> &m)
633{
634 for (int i = 0; i < NumCol; i++) {
635 if (!is_unit_scale(m[i])) {
636 return false;
637 }
638 }
639 return true;
640}
641
642template<typename T, int Size>
644{
645 bool success;
646 /* Explicit template parameter to please MSVC. */
647 return invert<T, Size>(mat, success);
648}
649
650template<typename T, int NumCol, int NumRow>
652{
654 unroll<NumCol>([&](auto i) { unroll<NumRow>([&](auto j) { result[i][j] = mat[j][i]; }); });
655 return result;
656}
657
658template<typename T, int NumCol, int NumRow, typename VectorT>
660 const VectorT &translation)
661{
662 using MatT = MatBase<T, NumCol, NumRow>;
663 BLI_STATIC_ASSERT(VectorT::type_length <= MatT::col_len - 1,
664 "Translation should be at least 1 column less than the matrix.");
665 constexpr int location_col = MatT::col_len - 1;
666 /* Avoid multiplying the last row if it exists.
667 * Allows using non square matrices like float3x2 and saves computation. */
668 using IntermediateVecT =
669 VecBase<typename MatT::base_type,
670 (MatT::row_len > MatT::col_len - 1) ? (MatT::col_len - 1) : MatT::row_len>;
671
672 MatT result = mat;
673 unroll<VectorT::type_length>([&](auto c) {
674 *reinterpret_cast<IntermediateVecT *>(
675 &result[location_col]) += translation[c] *
676 *reinterpret_cast<const IntermediateVecT *>(&mat[c]);
677 });
678 return result;
679}
680
681template<typename T, int NumCol, int NumRow, typename AngleT>
683 const AxisAngleBase<T, AngleT> &rotation)
684{
685 using MatT = MatBase<T, NumCol, NumRow>;
686 using Vec3T = typename MatT::vec3_type;
687 const T angle_sin = sin(rotation.angle());
688 const T angle_cos = cos(rotation.angle());
689 const Vec3T &axis_vec = rotation.axis();
690
691 MatT result = mat;
692 /* axis_vec is given to be normalized. */
693 if (axis_vec.x == T(1)) {
694 unroll<MatT::row_len>([&](auto c) {
695 result[2][c] = -angle_sin * mat[1][c] + angle_cos * mat[2][c];
696 result[1][c] = angle_cos * mat[1][c] + angle_sin * mat[2][c];
697 });
698 }
699 else if (axis_vec.y == T(1)) {
700 unroll<MatT::row_len>([&](auto c) {
701 result[0][c] = angle_cos * mat[0][c] - angle_sin * mat[2][c];
702 result[2][c] = angle_sin * mat[0][c] + angle_cos * mat[2][c];
703 });
704 }
705 else if (axis_vec.z == T(1)) {
706 unroll<MatT::row_len>([&](auto c) {
707 result[0][c] = angle_cos * mat[0][c] + angle_sin * mat[1][c];
708 result[1][c] = -angle_sin * mat[0][c] + angle_cos * mat[1][c];
709 });
710 }
711 else {
712 /* Un-optimized case. Arbitrary */
713 result *= from_rotation<MatT>(rotation);
714 }
715 return result;
716}
717
718template<typename T, int NumCol, int NumRow, typename VectorT>
720 const VectorT &scale)
721{
722 BLI_STATIC_ASSERT(VectorT::type_length <= NumCol,
723 "Scale should be less or equal to the matrix in column count.");
725 unroll<VectorT::type_length>([&](auto c) { result[c] *= scale[c]; });
726 return result;
727}
728
729template<typename T, int NumCol, int NumRow>
732 T t)
733{
735 unroll<NumCol>([&](auto c) { result[c] = interpolate(a[c], b[c], t); });
736 return result;
737}
738
739template<typename T,
740 int NumCol,
741 int NumRow,
742 int SrcNumCol,
743 int SrcNumRow,
744 int SrcStartCol,
745 int SrcStartRow,
746 int SrcAlignment>
755
756template<typename T,
757 int NumCol,
758 int NumRow,
759 int SrcNumCol,
760 int SrcNumRow,
761 int SrcStartCol,
762 int SrcStartRow,
763 int SrcAlignment,
764 typename VectorT>
767 &a,
768 VectorT &r_size)
769{
770 BLI_STATIC_ASSERT(VectorT::type_length == NumCol,
771 "r_size dimension should be equal to matrix column count.");
773 unroll<NumCol>([&](auto i) { result[i] = math::normalize_and_get_length(a[i], r_size[i]); });
774 return result;
775}
776
777template<typename T, int NumCol, int NumRow>
784
785template<typename T, int NumCol, int NumRow, typename VectorT>
787 const MatBase<T, NumCol, NumRow> &a, VectorT &r_size)
788{
789 BLI_STATIC_ASSERT(VectorT::type_length == NumCol,
790 "r_size dimension should be equal to matrix column count.");
792 unroll<NumCol>([&](auto i) { result[i] = math::normalize_and_get_length(a[i], r_size[i]); });
793 return result;
794}
795
796namespace detail {
797
799{
801 return AngleRadianBase(mat[0][0], mat[0][1]);
802}
803
804template<typename T>
806{
808
809 const T cy = math::hypot(mat[0][0], mat[0][1]);
810 if (cy > T(16) * std::numeric_limits<T>::epsilon()) {
811 eul1.x() = math::atan2(mat[1][2], mat[2][2]);
812 eul1.y() = math::atan2(-mat[0][2], cy);
813 eul1.z() = math::atan2(mat[0][1], mat[0][0]);
814
815 eul2.x() = math::atan2(-mat[1][2], -mat[2][2]);
816 eul2.y() = math::atan2(-mat[0][2], -cy);
817 eul2.z() = math::atan2(-mat[0][1], -mat[0][0]);
818 }
819 else {
820 eul1.x() = math::atan2(-mat[2][1], mat[1][1]);
821 eul1.y() = math::atan2(-mat[0][2], cy);
822 eul1.z() = 0.0f;
823
824 eul2 = eul1;
825 }
826}
827template<typename T>
829{
831 const int i_index = eul1.i_index();
832 const int j_index = eul1.j_index();
833 const int k_index = eul1.k_index();
834
835 const T cy = math::hypot(mat[i_index][i_index], mat[i_index][j_index]);
836 if (cy > T(16) * std::numeric_limits<T>::epsilon()) {
837 eul1.i() = math::atan2(mat[j_index][k_index], mat[k_index][k_index]);
838 eul1.j() = math::atan2(-mat[i_index][k_index], cy);
839 eul1.k() = math::atan2(mat[i_index][j_index], mat[i_index][i_index]);
840
841 eul2.i() = math::atan2(-mat[j_index][k_index], -mat[k_index][k_index]);
842 eul2.j() = math::atan2(-mat[i_index][k_index], -cy);
843 eul2.k() = math::atan2(-mat[i_index][j_index], -mat[i_index][i_index]);
844 }
845 else {
846 eul1.i() = math::atan2(-mat[k_index][j_index], mat[j_index][j_index]);
847 eul1.j() = math::atan2(-mat[i_index][k_index], cy);
848 eul1.k() = 0.0f;
849
850 eul2 = eul1;
851 }
852
853 if (eul1.parity()) {
854 eul1 = -eul1;
855 eul2 = -eul2;
856 }
857}
858
859/* Using explicit template instantiations in order to reduce compilation time. */
860extern template void normalized_to_eul2(const float3x3 &mat,
861 Euler3Base<float> &eul1,
862 Euler3Base<float> &eul2);
863extern template void normalized_to_eul2(const float3x3 &mat,
865 EulerXYZBase<float> &eul2);
866extern template void normalized_to_eul2(const double3x3 &mat,
869
871{
873 /* Caller must ensure matrices aren't negative for valid results, see: #24291, #94231. */
875
877
878 /* Method outlined by Mike Day, ref: https://math.stackexchange.com/a/3183435/220949
879 * with an additional `sqrtf(..)` for higher precision result.
880 * Removing the `sqrt` causes tests to fail unless the precision is set to 1e-6 or larger. */
881
882 if (mat[2][2] < 0.0f) {
883 if (mat[0][0] > mat[1][1]) {
884 const T trace = 1.0f + mat[0][0] - mat[1][1] - mat[2][2];
885 T s = 2.0f * math::sqrt(trace);
886 if (mat[1][2] < mat[2][1]) {
887 /* Ensure W is non-negative for a canonical result. */
888 s = -s;
889 }
890 q.x = 0.25f * s;
891 s = 1.0f / s;
892 q.w = (mat[1][2] - mat[2][1]) * s;
893 q.y = (mat[0][1] + mat[1][0]) * s;
894 q.z = (mat[2][0] + mat[0][2]) * s;
895 if (UNLIKELY((trace == 1.0f) && (q.w == 0.0f && q.y == 0.0f && q.z == 0.0f))) {
896 /* Avoids the need to normalize the degenerate case. */
897 q.x = 1.0f;
898 }
899 }
900 else {
901 const T trace = 1.0f - mat[0][0] + mat[1][1] - mat[2][2];
902 T s = 2.0f * math::sqrt(trace);
903 if (mat[2][0] < mat[0][2]) {
904 /* Ensure W is non-negative for a canonical result. */
905 s = -s;
906 }
907 q.y = 0.25f * s;
908 s = 1.0f / s;
909 q.w = (mat[2][0] - mat[0][2]) * s;
910 q.x = (mat[0][1] + mat[1][0]) * s;
911 q.z = (mat[1][2] + mat[2][1]) * s;
912 if (UNLIKELY((trace == 1.0f) && (q.w == 0.0f && q.x == 0.0f && q.z == 0.0f))) {
913 /* Avoids the need to normalize the degenerate case. */
914 q.y = 1.0f;
915 }
916 }
917 }
918 else {
919 if (mat[0][0] < -mat[1][1]) {
920 const T trace = 1.0f - mat[0][0] - mat[1][1] + mat[2][2];
921 T s = 2.0f * math::sqrt(trace);
922 if (mat[0][1] < mat[1][0]) {
923 /* Ensure W is non-negative for a canonical result. */
924 s = -s;
925 }
926 q.z = 0.25f * s;
927 s = 1.0f / s;
928 q.w = (mat[0][1] - mat[1][0]) * s;
929 q.x = (mat[2][0] + mat[0][2]) * s;
930 q.y = (mat[1][2] + mat[2][1]) * s;
931 if (UNLIKELY((trace == 1.0f) && (q.w == 0.0f && q.x == 0.0f && q.y == 0.0f))) {
932 /* Avoids the need to normalize the degenerate case. */
933 q.z = 1.0f;
934 }
935 }
936 else {
937 /* NOTE(@ideasman42): A zero matrix will fall through to this block,
938 * needed so a zero scaled matrices to return a quaternion without rotation, see: #101848.
939 */
940 const T trace = 1.0f + mat[0][0] + mat[1][1] + mat[2][2];
941 T s = 2.0f * math::sqrt(trace);
942 q.w = 0.25f * s;
943 s = 1.0f / s;
944 q.x = (mat[1][2] - mat[2][1]) * s;
945 q.y = (mat[2][0] - mat[0][2]) * s;
946 q.z = (mat[0][1] - mat[1][0]) * s;
947 if (UNLIKELY((trace == 1.0f) && (q.x == 0.0f && q.y == 0.0f && q.z == 0.0f))) {
948 /* Avoids the need to normalize the degenerate case. */
949 q.w = 1.0f;
950 }
951 }
952 }
953
954 BLI_assert(!(q.w < 0.0f));
955
956 /* Sometimes normalization is necessary due to round-off errors in the above
957 * calculations. The comparison here uses tighter tolerances than
958 * BLI_ASSERT_UNIT_QUAT(), so it's likely that even after a few more
959 * transformations the quaternion will still be considered unit-ish. */
960 const T q_len_squared = q.x * q.x + q.y * q.y + q.z * q.z + q.w * q.w;
961 const T threshold = 0.0002f /* #BLI_ASSERT_UNIT_EPSILON */ * 3;
962 if (math::abs(q_len_squared - 1.0f) >= threshold) {
963 const T q_len_inv = 1.0 / math::sqrt(q_len_squared);
964 q.x *= q_len_inv;
965 q.y *= q_len_inv;
966 q.z *= q_len_inv;
967 q.w *= q_len_inv;
968 }
969
971 return q;
972}
973
975{
976 const T det = math::determinant(mat);
977 if (UNLIKELY(!std::isfinite(det))) {
979 }
980 if (UNLIKELY(det < T(0))) {
981 return normalized_to_quat_fast(-mat);
982 }
983 return normalized_to_quat_fast(mat);
984}
985
986/* Using explicit template instantiations in order to reduce compilation time. */
989
990template<typename T, int NumCol, int NumRow>
992{
993 using MatT = MatBase<T, NumCol, NumRow>;
994 using DoublePrecision = typename TypeTraits<T>::DoublePrecision;
995 const DoublePrecision cos_i = math::cos(DoublePrecision(rotation.x().radian()));
996 const DoublePrecision cos_j = math::cos(DoublePrecision(rotation.y().radian()));
997 const DoublePrecision cos_k = math::cos(DoublePrecision(rotation.z().radian()));
998 const DoublePrecision sin_i = math::sin(DoublePrecision(rotation.x().radian()));
999 const DoublePrecision sin_j = math::sin(DoublePrecision(rotation.y().radian()));
1000 const DoublePrecision sin_k = math::sin(DoublePrecision(rotation.z().radian()));
1001 const DoublePrecision cos_i_cos_k = cos_i * cos_k;
1002 const DoublePrecision cos_i_sin_k = cos_i * sin_k;
1003 const DoublePrecision sin_i_cos_k = sin_i * cos_k;
1004 const DoublePrecision sin_i_sin_k = sin_i * sin_k;
1005
1006 MatT mat = MatT::identity();
1007 mat[0][0] = T(cos_j * cos_k);
1008 mat[1][0] = T(sin_j * sin_i_cos_k - cos_i_sin_k);
1009 mat[2][0] = T(sin_j * cos_i_cos_k + sin_i_sin_k);
1010
1011 mat[0][1] = T(cos_j * sin_k);
1012 mat[1][1] = T(sin_j * sin_i_sin_k + cos_i_cos_k);
1013 mat[2][1] = T(sin_j * cos_i_sin_k - sin_i_cos_k);
1014
1015 mat[0][2] = T(-sin_j);
1016 mat[1][2] = T(cos_j * sin_i);
1017 mat[2][2] = T(cos_j * cos_i);
1018 return mat;
1019}
1020
1021template<typename T, int NumCol, int NumRow>
1023{
1024 using MatT = MatBase<T, NumCol, NumRow>;
1025 const int i_index = rotation.i_index();
1026 const int j_index = rotation.j_index();
1027 const int k_index = rotation.k_index();
1028#if 1 /* Reference. */
1029 EulerXYZBase<T> euler_xyz(rotation.ijk());
1030 const MatT mat = from_rotation<T, NumCol, NumRow>(rotation.parity() ? -euler_xyz : euler_xyz);
1031 MatT result = MatT::identity();
1032 result[i_index][i_index] = mat[0][0];
1033 result[j_index][i_index] = mat[1][0];
1034 result[k_index][i_index] = mat[2][0];
1035 result[i_index][j_index] = mat[0][1];
1036 result[j_index][j_index] = mat[1][1];
1037 result[k_index][j_index] = mat[2][1];
1038 result[i_index][k_index] = mat[0][2];
1039 result[j_index][k_index] = mat[1][2];
1040 result[k_index][k_index] = mat[2][2];
1041#else
1042 /* TODO(fclem): Manually inline and check performance difference. */
1043#endif
1044 return result;
1045}
1046
1047template<typename T, int NumCol, int NumRow>
1049{
1050 using MatT = MatBase<T, NumCol, NumRow>;
1051 using DoublePrecision = typename TypeTraits<T>::DoublePrecision;
1052 const DoublePrecision q0 = numbers::sqrt2 * DoublePrecision(rotation.w);
1053 const DoublePrecision q1 = numbers::sqrt2 * DoublePrecision(rotation.x);
1054 const DoublePrecision q2 = numbers::sqrt2 * DoublePrecision(rotation.y);
1055 const DoublePrecision q3 = numbers::sqrt2 * DoublePrecision(rotation.z);
1056
1057 const DoublePrecision qda = q0 * q1;
1058 const DoublePrecision qdb = q0 * q2;
1059 const DoublePrecision qdc = q0 * q3;
1060 const DoublePrecision qaa = q1 * q1;
1061 const DoublePrecision qab = q1 * q2;
1062 const DoublePrecision qac = q1 * q3;
1063 const DoublePrecision qbb = q2 * q2;
1064 const DoublePrecision qbc = q2 * q3;
1065 const DoublePrecision qcc = q3 * q3;
1066
1067 MatT mat = MatT::identity();
1068 mat[0][0] = T(1.0 - qbb - qcc);
1069 mat[0][1] = T(qdc + qab);
1070 mat[0][2] = T(-qdb + qac);
1071
1072 mat[1][0] = T(-qdc + qab);
1073 mat[1][1] = T(1.0 - qaa - qcc);
1074 mat[1][2] = T(qda + qbc);
1075
1076 mat[2][0] = T(qdb + qac);
1077 mat[2][1] = T(-qda + qbc);
1078 mat[2][2] = T(1.0 - qaa - qbb);
1079 return mat;
1080}
1081
1082/* Not technically speaking a simple rotation, but a whole transform. */
1083template<typename T, int NumCol, int NumRow>
1085{
1086 using MatT = MatBase<T, NumCol, NumRow>;
1087 BLI_assert(is_normalized(rotation));
1094 /* Follow the paper notation. */
1095 const QuaternionBase<T> &c0 = rotation.quat;
1096 const QuaternionBase<T> &ce = rotation.trans;
1097 const T &w0 = c0.w, &x0 = c0.x, &y0 = c0.y, &z0 = c0.z;
1098 const T &we = ce.w, &xe = ce.x, &ye = ce.y, &ze = ce.z;
1099 /* Rotation. */
1100 MatT mat = from_rotation<T, NumCol, NumRow>(c0);
1101 /* Translation. */
1102 mat[3][0] = T(2) * (-we * x0 + xe * w0 - ye * z0 + ze * y0);
1103 mat[3][1] = T(2) * (-we * y0 + xe * z0 + ye * w0 - ze * x0);
1104 mat[3][2] = T(2) * (-we * z0 - xe * y0 + ye * x0 + ze * w0);
1105 /* Scale. */
1106 if (rotation.scale_weight != T(0)) {
1107 mat.template view<4, 4>() = mat * rotation.scale;
1108 }
1109 return mat;
1110}
1111
1112template<typename T, int NumCol, int NumRow>
1114{
1115 using MatT = MatBase<T, NumCol, NumRow>;
1116 MatT mat = MatT::identity();
1117 mat.x_axis() = to_vector<VecBase<T, 3>>(rotation.axes.x);
1118 mat.y_axis() = to_vector<VecBase<T, 3>>(rotation.axes.y);
1119 mat.z_axis() = to_vector<VecBase<T, 3>>(rotation.axes.z);
1120 return mat;
1121}
1122
1123template<typename T, int NumCol, int NumRow, typename AngleT>
1125{
1126 using MatT = MatBase<T, NumCol, NumRow>;
1127 using Vec3T = typename MatT::vec3_type;
1128 const T angle_sin = sin(rotation.angle());
1129 const T angle_cos = cos(rotation.angle());
1130 const Vec3T &axis = rotation.axis();
1131
1133
1134 T ico = (T(1) - angle_cos);
1135 Vec3T nsi = axis * angle_sin;
1136
1137 Vec3T n012 = (axis * axis) * ico;
1138 T n_01 = (axis[0] * axis[1]) * ico;
1139 T n_02 = (axis[0] * axis[2]) * ico;
1140 T n_12 = (axis[1] * axis[2]) * ico;
1141
1142 MatT mat = from_scale<MatT>(n012 + angle_cos);
1143 mat[0][1] = n_01 + nsi[2];
1144 mat[0][2] = n_02 - nsi[1];
1145 mat[1][0] = n_01 - nsi[2];
1146 mat[1][2] = n_12 + nsi[0];
1147 mat[2][0] = n_02 + nsi[1];
1148 mat[2][1] = n_12 - nsi[0];
1149 return mat;
1150}
1151
1152template<typename T, int NumCol, int NumRow>
1154{
1155 using MatT = MatBase<T, NumCol, NumRow>;
1156 const T cos_i = cos(rotation);
1157 const T sin_i = sin(rotation);
1158
1159 MatT mat = MatT::identity();
1160 mat[0][0] = cos_i;
1161 mat[1][0] = -sin_i;
1162
1163 mat[0][1] = sin_i;
1164 mat[1][1] = cos_i;
1165 return mat;
1166}
1167
1168/* Using explicit template instantiations in order to reduce compilation time. */
1169extern template MatBase<float, 2, 2> from_rotation(const AngleRadian &rotation);
1170extern template MatBase<float, 3, 3> from_rotation(const AngleRadian &rotation);
1171extern template MatBase<float, 3, 3> from_rotation(const EulerXYZ &rotation);
1172extern template MatBase<float, 4, 4> from_rotation(const EulerXYZ &rotation);
1173extern template MatBase<float, 3, 3> from_rotation(const Euler3 &rotation);
1174extern template MatBase<float, 4, 4> from_rotation(const Euler3 &rotation);
1175extern template MatBase<float, 3, 3> from_rotation(const Quaternion &rotation);
1176extern template MatBase<float, 4, 4> from_rotation(const Quaternion &rotation);
1177extern template MatBase<float, 3, 3> from_rotation(const AxisAngle &rotation);
1178extern template MatBase<float, 4, 4> from_rotation(const AxisAngle &rotation);
1180extern template MatBase<float, 4, 4> from_rotation(const AxisAngleCartesian &rotation);
1181
1182} // namespace detail
1183
1184template<typename T> [[nodiscard]] inline AngleRadianBase<T> to_angle(const MatBase<T, 2, 2> &mat)
1185{
1186 return detail::normalized_to_angle(mat);
1187}
1188
1189template<typename T>
1190[[nodiscard]] inline Euler3Base<T> to_euler(const MatBase<T, 3, 3> &mat, EulerOrder order)
1191{
1192 Euler3Base<T> eul1(order), eul2(order);
1193 detail::normalized_to_eul2(mat, eul1, eul2);
1194 /* Return best, which is just the one with lowest values in it. */
1195 return (length_manhattan(VecBase<T, 3>(eul1)) > length_manhattan(VecBase<T, 3>(eul2))) ? eul2 :
1196 eul1;
1197}
1198
1199template<typename T> [[nodiscard]] inline EulerXYZBase<T> to_euler(const MatBase<T, 3, 3> &mat)
1200{
1201 EulerXYZBase<T> eul1, eul2;
1202 detail::normalized_to_eul2(mat, eul1, eul2);
1203 /* Return best, which is just the one with lowest values in it. */
1204 return (length_manhattan(VecBase<T, 3>(eul1)) > length_manhattan(VecBase<T, 3>(eul2))) ? eul2 :
1205 eul1;
1206}
1207
1208template<typename T>
1209[[nodiscard]] inline Euler3Base<T> to_euler(const MatBase<T, 4, 4> &mat, EulerOrder order)
1210{
1211 /* TODO(fclem): Avoid the copy with 3x3 ref. */
1212 return to_euler<T>(MatBase<T, 3, 3>(mat), order);
1213}
1214
1215template<typename T> [[nodiscard]] inline EulerXYZBase<T> to_euler(const MatBase<T, 4, 4> &mat)
1216{
1217 /* TODO(fclem): Avoid the copy with 3x3 ref. */
1218 return to_euler<T>(MatBase<T, 3, 3>(mat));
1219}
1220
1221template<typename T>
1222[[nodiscard]] inline Euler3Base<T> to_nearest_euler(const MatBase<T, 3, 3> &mat,
1223 const Euler3Base<T> &reference)
1224{
1225 Euler3Base<T> eul1(reference.order()), eul2(reference.order());
1226 detail::normalized_to_eul2(mat, eul1, eul2);
1227 eul1 = eul1.wrapped_around(reference);
1228 eul2 = eul2.wrapped_around(reference);
1229 /* Return best, which is just the one with lowest values it in. */
1230 return (length_manhattan(VecBase<T, 3>(eul1) - VecBase<T, 3>(reference)) >
1231 length_manhattan(VecBase<T, 3>(eul2) - VecBase<T, 3>(reference))) ?
1232 eul2 :
1233 eul1;
1234}
1235
1236template<typename T>
1237[[nodiscard]] inline EulerXYZBase<T> to_nearest_euler(const MatBase<T, 3, 3> &mat,
1238 const EulerXYZBase<T> &reference)
1239{
1240 EulerXYZBase<T> eul1, eul2;
1241 detail::normalized_to_eul2(mat, eul1, eul2);
1242 eul1 = eul1.wrapped_around(reference);
1243 eul2 = eul2.wrapped_around(reference);
1244 /* Return best, which is just the one with lowest values it in. */
1245 return (length_manhattan(VecBase<T, 3>(eul1) - VecBase<T, 3>(reference)) >
1246 length_manhattan(VecBase<T, 3>(eul2) - VecBase<T, 3>(reference))) ?
1247 eul2 :
1248 eul1;
1249}
1250
1251template<typename T>
1252[[nodiscard]] inline Euler3Base<T> to_nearest_euler(const MatBase<T, 4, 4> &mat,
1253 const Euler3Base<T> &reference)
1254{
1255 /* TODO(fclem): Avoid the copy with 3x3 ref. */
1256 return to_euler<T>(MatBase<T, 3, 3>(mat), reference);
1257}
1258
1259template<typename T>
1260[[nodiscard]] inline EulerXYZBase<T> to_nearest_euler(const MatBase<T, 4, 4> &mat,
1261 const EulerXYZBase<T> &reference)
1262{
1263 /* TODO(fclem): Avoid the copy with 3x3 ref. */
1264 return to_euler<T>(MatBase<T, 3, 3>(mat), reference);
1265}
1266
1267template<typename T>
1268[[nodiscard]] inline QuaternionBase<T> to_quaternion(const MatBase<T, 3, 3> &mat)
1269{
1271}
1272
1273template<typename T>
1274[[nodiscard]] inline QuaternionBase<T> to_quaternion(const MatBase<T, 4, 4> &mat)
1275{
1276 /* TODO(fclem): Avoid the copy with 3x3 ref. */
1278}
1279
1284template<typename T>
1286{
1287 /* Conversion to quaternion asserts when the matrix contains some kinds of shearing, conversion
1288 * to euler does not. */
1289 /* TODO: Find a better algorithm that can convert untrusted matrices to quaternions directly. */
1290 return to_quaternion(to_euler(mat));
1291}
1292
1293template<bool AllowNegativeScale, typename T, int NumCol, int NumRow>
1294[[nodiscard]] inline VecBase<T, 3> to_scale(const MatBase<T, NumCol, NumRow> &mat)
1295{
1296 VecBase<T, 3> result = {length(mat.x_axis()), length(mat.y_axis()), length(mat.z_axis())};
1297 if constexpr (AllowNegativeScale) {
1298 if (UNLIKELY(is_negative(mat))) {
1299 result = -result;
1300 }
1301 }
1302 return result;
1303}
1304
1305template<bool AllowNegativeScale, typename T>
1306[[nodiscard]] inline VecBase<T, 2> to_scale(const MatBase<T, 2, 2> &mat)
1307{
1308 VecBase<T, 2> result = {length(mat.x), length(mat.y)};
1309 if constexpr (AllowNegativeScale) {
1310 if (UNLIKELY(is_negative(mat))) {
1311 result = -result;
1312 }
1313 }
1314 return result;
1315}
1316
1317/* Implementation details. Use `to_euler` and `to_quaternion` instead. */
1318namespace detail {
1319
1320template<typename T>
1321inline void to_rotation(const MatBase<T, 2, 2> &mat, AngleRadianBase<T> &r_rotation)
1322{
1323 r_rotation = to_angle<T>(mat);
1324}
1325
1326template<typename T>
1327inline void to_rotation(const MatBase<T, 3, 3> &mat, QuaternionBase<T> &r_rotation)
1328{
1329 r_rotation = to_quaternion<T>(mat);
1330}
1331
1332template<typename T>
1333inline void to_rotation(const MatBase<T, 3, 3> &mat, EulerXYZBase<T> &r_rotation)
1334{
1335 r_rotation = to_euler<T>(mat);
1336}
1337
1338template<typename T>
1339inline void to_rotation(const MatBase<T, 3, 3> &mat, Euler3Base<T> &r_rotation)
1340{
1341 r_rotation = to_euler<T>(mat, r_rotation.order());
1342}
1343
1344} // namespace detail
1345
1346template<bool AllowNegativeScale, typename T>
1347inline void to_rot_scale(const MatBase<T, 2, 2> &mat,
1348 AngleRadianBase<T> &r_rotation,
1349 VecBase<T, 2> &r_scale)
1350{
1351 MatBase<T, 2, 2> normalized_mat = normalize_and_get_size(mat, r_scale);
1352 if constexpr (AllowNegativeScale) {
1353 if (UNLIKELY(is_negative(normalized_mat))) {
1354 normalized_mat = -normalized_mat;
1355 r_scale = -r_scale;
1356 }
1357 }
1358 detail::to_rotation<T>(normalized_mat, r_rotation);
1359}
1360
1361template<bool AllowNegativeScale, typename T>
1362inline void to_loc_rot_scale(const MatBase<T, 3, 3> &mat,
1363 VecBase<T, 2> &r_location,
1364 AngleRadianBase<T> &r_rotation,
1365 VecBase<T, 2> &r_scale)
1366{
1367 r_location = mat.location();
1368 to_rot_scale<AllowNegativeScale>(MatBase<T, 2, 2>(mat), r_rotation, r_scale);
1369}
1370
1371template<bool AllowNegativeScale, typename T, typename RotationT>
1372inline void to_rot_scale(const MatBase<T, 3, 3> &mat,
1373 RotationT &r_rotation,
1374 VecBase<T, 3> &r_scale)
1375{
1376 MatBase<T, 3, 3> normalized_mat = normalize_and_get_size(mat, r_scale);
1377 if constexpr (AllowNegativeScale) {
1378 if (UNLIKELY(is_negative(normalized_mat))) {
1379 normalized_mat = -normalized_mat;
1380 r_scale = -r_scale;
1381 }
1382 }
1383 detail::to_rotation<T>(normalized_mat, r_rotation);
1384}
1385
1386template<bool AllowNegativeScale, typename T, typename RotationT>
1387inline void to_loc_rot_scale(const MatBase<T, 4, 4> &mat,
1388 VecBase<T, 3> &r_location,
1389 RotationT &r_rotation,
1390 VecBase<T, 3> &r_scale)
1391{
1392 r_location = mat.location();
1393 to_rot_scale<AllowNegativeScale>(MatBase<T, 3, 3>(mat), r_rotation, r_scale);
1394}
1395
1400template<bool AllowNegativeScale, typename T, typename RotationT>
1402 VecBase<T, 3> &r_location,
1403 RotationT &r_rotation,
1404 VecBase<T, 3> &r_scale)
1405{
1406 EulerXYZBase<T> euler_rotation;
1407 to_loc_rot_scale<AllowNegativeScale>(mat, r_location, euler_rotation, r_scale);
1408 if constexpr (std::is_same_v<std::decay_t<RotationT>, QuaternionBase<T>>) {
1409 r_rotation = to_quaternion(euler_rotation);
1410 }
1411 else {
1412 r_rotation = RotationT(euler_rotation);
1413 }
1414}
1415
1416template<typename MatT> [[nodiscard]] MatT from_location(const typename MatT::loc_type &location)
1417{
1418 MatT mat = MatT::identity();
1419 mat.location() = location;
1420 return mat;
1421}
1422
1423template<typename MatT, int ScaleDim>
1425{
1426 BLI_STATIC_ASSERT(ScaleDim <= MatT::min_dim,
1427 "Scale dimension should fit the matrix diagonal length.");
1428 MatT result{};
1430 [&](auto i) { result[i][i] = (i < ScaleDim) ? scale[i] : typename MatT::base_type(1); });
1431 return result;
1432}
1433
1434template<typename MatT, typename RotationT>
1435[[nodiscard]] MatT from_rotation(const RotationT &rotation)
1436{
1438}
1439
1440template<typename MatT, typename RotationT, typename VectorT>
1441[[nodiscard]] MatT from_rot_scale(const RotationT &rotation, const VectorT &scale)
1442{
1443 return from_rotation<MatT>(rotation) * from_scale<MatT>(scale);
1444}
1445
1446template<typename MatT, typename RotationT, int ScaleDim>
1447[[nodiscard]] MatT from_loc_rot_scale(const typename MatT::loc_type &location,
1448 const RotationT &rotation,
1450{
1451 using MatRotT =
1453 MatT mat = MatT(from_rot_scale<MatRotT>(rotation, scale));
1454 mat.location() = location;
1455 return mat;
1456}
1457
1458template<typename MatT, typename RotationT>
1459[[nodiscard]] MatT from_loc_rot(const typename MatT::loc_type &location, const RotationT &rotation)
1460{
1461 using MatRotT =
1463 MatT mat = MatT(from_rotation<MatRotT>(rotation));
1464 mat.location() = location;
1465 return mat;
1466}
1467
1468template<typename MatT, int ScaleDim>
1469[[nodiscard]] MatT from_loc_scale(const typename MatT::loc_type &location,
1471{
1472 MatT mat = MatT(from_scale<MatT>(scale));
1473 mat.location() = location;
1474 return mat;
1475}
1476
1477template<typename T> MatBase<T, 2, 2> from_direction(const VecBase<T, 2> &direction)
1478{
1479 BLI_assert(is_unit_scale(direction));
1480 return MatBase<T, 2, 2>(direction,
1481 VecBase<T, 2>(direction.y, direction.x) * VecBase<T, 2>(-1, 1));
1482}
1483
1484template<typename MatT, typename VectorT>
1485[[nodiscard]] MatT from_orthonormal_axes(const VectorT forward, const VectorT up)
1486{
1487 BLI_assert(is_unit_scale(forward));
1489
1490 /* TODO(fclem): This is wrong. Forward is Y. */
1491 MatT matrix;
1492 matrix.x_axis() = forward;
1493 /* Beware of handedness! Blender uses right-handedness.
1494 * Resulting matrix should have determinant of 1. */
1495 matrix.y_axis() = math::cross(up, forward);
1496 matrix.z_axis() = up;
1497 return matrix;
1498}
1499
1500template<typename MatT, typename VectorT>
1501[[nodiscard]] MatT from_orthonormal_axes(const VectorT location,
1502 const VectorT forward,
1503 const VectorT up)
1504{
1506 MatT matrix = MatT(from_orthonormal_axes<Mat3x3>(forward, up));
1507 matrix.location() = location;
1508 return matrix;
1509}
1510
1511template<typename MatT, typename VectorT> [[nodiscard]] MatT from_up_axis(const VectorT up)
1512{
1514 using T = typename MatT::base_type;
1515 using Vec3T = VecBase<T, 3>;
1516 /* Duff, Tom, et al. "Building an orthonormal basis, revisited." JCGT 6.1 (2017). */
1517 T sign = up.z >= T(0) ? T(1) : T(-1);
1518 T a = T(-1) / (sign + up.z);
1519 T b = up.x * up.y * a;
1520
1521 MatBase<T, 3, 3> basis;
1522 basis.x_axis() = Vec3T(1.0f + sign * square(up.x) * a, sign * b, -sign * up.x);
1523 basis.y_axis() = Vec3T(b, sign + square(up.y) * a, -up.y);
1524 basis.z_axis() = up;
1525 return MatT(basis);
1526}
1527
1528template<typename MatT> [[nodiscard]] MatT orthogonalize(const MatT &mat, const Axis axis)
1529{
1530 using T = typename MatT::base_type;
1531 using Vec3T = VecBase<T, 3>;
1532 Vec3T scale;
1534 R.x = normalize_and_get_length(mat.x_axis(), scale.x);
1535 R.y = normalize_and_get_length(mat.y_axis(), scale.y);
1536 R.z = normalize_and_get_length(mat.z_axis(), scale.z);
1537 /* NOTE(fclem) This is a direct port from `orthogonalize_m4()`.
1538 * To select the secondary axis, it checks if the candidate axis is not colinear.
1539 * The issue is that the candidate axes are not normalized so this dot product
1540 * check is kind of pointless.
1541 * Because of this, the target axis could still be colinear but pass the check. */
1542#if 1 /* Reproduce C API behavior. Do not normalize other axes. */
1543 switch (axis) {
1544 case Axis::X:
1545 R.y = mat.y_axis();
1546 R.z = mat.z_axis();
1547 break;
1548 case Axis::Y:
1549 R.x = mat.x_axis();
1550 R.z = mat.z_axis();
1551 break;
1552 case Axis::Z:
1553 R.x = mat.x_axis();
1554 R.y = mat.y_axis();
1555 break;
1556 }
1557#endif
1558
1564 switch (axis) {
1565 case Axis::X:
1566 if (dot(R.x, R.y) < T(1)) {
1567 R.z = normalize(cross(R.x, R.y));
1568 R.y = cross(R.z, R.x);
1569 }
1570 else if (dot(R.x, R.z) < T(1)) {
1571 R.y = normalize(cross(R.z, R.x));
1572 R.z = cross(R.x, R.y);
1573 }
1574 else {
1575 R.z = normalize(cross(R.x, Vec3T(R.x.y, R.x.z, R.x.x)));
1576 R.y = cross(R.z, R.x);
1577 }
1578 break;
1579 case Axis::Y:
1580 if (dot(R.y, R.x) < T(1)) {
1581 R.z = normalize(cross(R.x, R.y));
1582 R.x = cross(R.y, R.z);
1583 }
1584 /* FIXME(fclem): THIS IS WRONG. Should be dot(R.y, R.z). Following C code for now... */
1585 else if (dot(R.x, R.z) < T(1)) {
1586 R.x = normalize(cross(R.y, R.z));
1587 R.z = cross(R.x, R.y);
1588 }
1589 else {
1590 R.x = normalize(cross(R.y, Vec3T(R.y.y, R.y.z, R.y.x)));
1591 R.z = cross(R.x, R.y);
1592 }
1593 break;
1594 case Axis::Z:
1595 if (dot(R.z, R.x) < T(1)) {
1596 R.y = normalize(cross(R.z, R.x));
1597 R.x = cross(R.y, R.z);
1598 }
1599 else if (dot(R.z, R.y) < T(1)) {
1600 R.x = normalize(cross(R.y, R.z));
1601 R.y = cross(R.z, R.x);
1602 }
1603 else {
1604 R.x = normalize(cross(Vec3T(R.z.y, R.z.z, R.z.x), R.z));
1605 R.y = cross(R.z, R.x);
1606 }
1607 break;
1608 default:
1610 break;
1611 }
1612 /* Reapply the lost scale. */
1613 R.x *= scale.x;
1614 R.y *= scale.y;
1615 R.z *= scale.z;
1616
1617 MatT result(R);
1618 result.location() = mat.location();
1619 return result;
1620}
1621
1622template<typename MatT, typename VectorT>
1623[[nodiscard]] MatT from_origin_transform(const MatT &transform, const VectorT origin)
1624{
1625 return from_location<MatT>(origin) * transform * from_location<MatT>(-origin);
1626}
1627
1628template<typename T>
1630{
1631 return mat * point;
1632}
1633
1634template<typename T>
1636{
1637 return mat.template view<2, 2>() * point + mat.location();
1638}
1639
1640template<typename T>
1642{
1643 return mat * point;
1644}
1645
1646template<typename T>
1648{
1649 return mat.template view<3, 3>() * point + mat.location();
1650}
1651
1652template<typename T>
1654{
1655 return mat * direction;
1656}
1657
1658template<typename T>
1660{
1661 return mat.template view<3, 3>() * direction;
1662}
1663
1664template<typename T, int N, int NumRow>
1666{
1667 VecBase<T, N + 1> tmp(point, T(1));
1668 tmp = mat * tmp;
1669 /* Absolute value to not flip the frustum upside down behind the camera. */
1670 return VecBase<T, N>(tmp) / math::abs(tmp[N]);
1671}
1672
1673extern template float3 transform_point(const float3x3 &mat, const float3 &point);
1674extern template float3 transform_point(const float4x4 &mat, const float3 &point);
1675extern template float3 transform_direction(const float3x3 &mat, const float3 &direction);
1676extern template float3 transform_direction(const float4x4 &mat, const float3 &direction);
1677extern template float3 project_point(const float4x4 &mat, const float3 &point);
1678extern template float2 project_point(const float3x3 &mat, const float2 &point);
1679
1680namespace projection {
1681
1682template<typename T>
1683MatBase<T, 4, 4> orthographic(T left, T right, T bottom, T top, T near_clip, T far_clip)
1684{
1685 const T x_delta = right - left;
1686 const T y_delta = top - bottom;
1687 const T z_delta = far_clip - near_clip;
1688
1690 if (x_delta != 0 && y_delta != 0 && z_delta != 0) {
1691 mat[0][0] = T(2.0) / x_delta;
1692 mat[3][0] = -(right + left) / x_delta;
1693 mat[1][1] = T(2.0) / y_delta;
1694 mat[3][1] = -(top + bottom) / y_delta;
1695 mat[2][2] = -T(2.0) / z_delta; /* NOTE: negate Z. */
1696 mat[3][2] = -(far_clip + near_clip) / z_delta;
1697 }
1698 return mat;
1699}
1700
1701template<typename T>
1702MatBase<T, 4, 4> orthographic_infinite(T left, T right, T bottom, T top, T near_clip)
1703{
1704 const T x_delta = right - left;
1705 const T y_delta = top - bottom;
1706
1708 if (x_delta != 0 && y_delta != 0) {
1709 mat[0][0] = T(2.0) / x_delta;
1710 mat[3][0] = -(right + left) / x_delta;
1711 mat[1][1] = T(2.0) / y_delta;
1712 mat[3][1] = -(top + bottom) / y_delta;
1713 /* Page 17. Choosing an epsilon for 32 bit floating-point precision. */
1714 constexpr float eps = 2.4e-7f;
1715 /* From "Projection Matrix Tricks" by Eric Lengyel GDC 2007.
1716 * Following same procedure as the reference but for orthographic matrix.
1717 * This avoids degenerate matrix (0 determinant). */
1718 mat[2][2] = -eps;
1719 mat[3][2] = -1.0f - eps * near_clip;
1720 }
1721 return mat;
1722}
1723
1724template<typename T>
1725MatBase<T, 4, 4> perspective(T left, T right, T bottom, T top, T near_clip, T far_clip)
1726{
1727 const T x_delta = right - left;
1728 const T y_delta = top - bottom;
1729 const T z_delta = far_clip - near_clip;
1730
1732 if (x_delta != 0 && y_delta != 0 && z_delta != 0) {
1733 mat[0][0] = near_clip * T(2.0) / x_delta;
1734 mat[1][1] = near_clip * T(2.0) / y_delta;
1735 mat[2][0] = (right + left) / x_delta; /* NOTE: negate Z. */
1736 mat[2][1] = (top + bottom) / y_delta;
1737 mat[2][2] = -(far_clip + near_clip) / z_delta;
1738 mat[2][3] = -1.0f;
1739 mat[3][2] = (-2.0f * near_clip * far_clip) / z_delta;
1740 mat[3][3] = 0.0f;
1741 }
1742 return mat;
1743}
1744
1745template<typename T>
1746MatBase<T, 4, 4> perspective_infinite(T left, T right, T bottom, T top, T near_clip)
1747{
1748 const T x_delta = right - left;
1749 const T y_delta = top - bottom;
1750
1751 /* From "Projection Matrix Tricks" by Eric Lengyel GDC 2007. */
1753 if (x_delta != 0 && y_delta != 0) {
1754 mat[0][0] = near_clip * T(2.0) / x_delta;
1755 mat[1][1] = near_clip * T(2.0) / y_delta;
1756 mat[2][0] = (right + left) / x_delta; /* NOTE: negate Z. */
1757 mat[2][1] = (top + bottom) / y_delta;
1758 /* Page 17. Choosing an epsilon for 32 bit floating-point precision. */
1759 constexpr float eps = 2.4e-7f;
1760 mat[2][2] = -1.0f;
1761 mat[2][3] = (eps - 1.0f);
1762 mat[3][2] = (eps - 2.0f) * near_clip;
1763 mat[3][3] = 0.0f;
1764 }
1765 return mat;
1766}
1767
1768template<typename T>
1770 T angle_left, T angle_right, T angle_bottom, T angle_top, T near_clip, T far_clip)
1771{
1772 MatBase<T, 4, 4> mat = perspective(math::tan(angle_left),
1773 math::tan(angle_right),
1774 math::tan(angle_bottom),
1775 math::tan(angle_top),
1776 near_clip,
1777 far_clip);
1778 mat[0][0] /= near_clip;
1779 mat[1][1] /= near_clip;
1780 return mat;
1781}
1782
1783template<typename T>
1784[[nodiscard]] MatBase<T, 4, 4> translate(const MatBase<T, 4, 4> &mat, const VecBase<T, 2> &offset)
1785{
1787 const bool is_perspective = mat[2][3] == -1.0f;
1788 const bool is_perspective_infinite = mat[2][2] == -1.0f;
1789 if (is_perspective || is_perspective_infinite) {
1790 result[2][0] -= mat[0][0] * offset.x / math::length(float3(mat[0][0], mat[1][0], mat[2][0]));
1791 result[2][1] -= mat[1][1] * offset.y / math::length(float3(mat[0][1], mat[1][1], mat[2][1]));
1792 }
1793 else {
1794 result[3][0] += offset.x;
1795 result[3][1] += offset.y;
1796 }
1797 return result;
1798}
1799
1800extern template float4x4 orthographic(
1801 float left, float right, float bottom, float top, float near_clip, float far_clip);
1802extern template float4x4 perspective(
1803 float left, float right, float bottom, float top, float near_clip, float far_clip);
1804
1805} // namespace projection
1806
1808
1815
1818 MutableSpan<float3> points,
1819 bool use_threading = true);
1821 const float4x4 &transform,
1823 bool use_threading = true);
1824
1825} // namespace blender::math
#define BLI_assert_unreachable()
Definition BLI_assert.h:93
#define BLI_STATIC_ASSERT(a, msg)
Definition BLI_assert.h:83
#define BLI_assert(a)
Definition BLI_assert.h:46
#define UNLIKELY(x)
static AppView * view
ATTR_WARN_UNUSED_RESULT const BMVert const BMEdge * e
SIMD_FORCE_INLINE btVector3 transform(const btVector3 &point) const
btMatrix3x3 adjoint() const
Return the adjoint of the matrix.
btScalar determinant() const
Return the determinant of the matrix.
SIMD_FORCE_INLINE const btScalar & z() const
Return the z value.
Definition btQuadWord.h:117
SIMD_FORCE_INLINE btScalar length() const
Return the length of the vector.
Definition btVector3.h:257
static constexpr Value Z
static constexpr Value X
static constexpr Value Y
static float normals[][3]
uint top
#define sin
#define cos
VecBase< float, D > normalize(VecOp< float, D >) RET
MatBase< R, C > transpose(MatBase< C, R >) RET
VecBase< float, 3 > cross(VecOp< float, 3 >, VecOp< float, 3 >) RET
static int left
#define N
#define T
#define R
QuaternionBase< T > normalized_to_quat_with_checks(const MatBase< T, 3, 3 > &mat)
void normalized_to_eul2(const MatBase< T, 3, 3 > &mat, EulerXYZBase< T > &eul1, EulerXYZBase< T > &eul2)
void to_rotation(const MatBase< T, 2, 2 > &mat, AngleRadianBase< T > &r_rotation)
QuaternionBase< T > normalized_to_quat_fast(const MatBase< T, 3, 3 > &mat)
AngleRadianBase< T > normalized_to_angle(const MatBase< T, 2, 2 > &mat)
MatBase< T, NumCol, NumRow > from_rotation(const AngleRadianBase< T > &rotation)
MatBase< T, 4, 4 > orthographic(T left, T right, T bottom, T top, T near_clip, T far_clip)
Create an orthographic projection matrix using OpenGL coordinate convention: Maps each axis range to ...
MatBase< T, 4, 4 > orthographic_infinite(T left, T right, T bottom, T top)
Create an orthographic projection matrix using OpenGL coordinate convention: Maps each axis range to ...
MatBase< T, 4, 4 > perspective_fov(T angle_left, T angle_right, T angle_bottom, T angle_top, T near_clip, T far_clip)
Create a perspective projection matrix using OpenGL coordinate convention: Maps each axis range to [-...
MatBase< T, 4, 4 > translate(const MatBase< T, 4, 4 > &mat, const VecBase< T, 2 > &offset)
Translate a projection matrix after creation in the screen plane. Usually used for anti-aliasing jitt...
MatBase< T, 4, 4 > perspective_infinite(T left, T right, T bottom, T top, T near_clip)
Create a perspective projection matrix using OpenGL coordinate convention: Maps each axis range to [-...
MatBase< T, 4, 4 > perspective(T left, T right, T bottom, T top, T near_clip, T far_clip)
Create a perspective projection matrix using OpenGL coordinate convention: Maps each axis range to [-...
AngleRadianBase< float > AngleRadian
AngleRadianBase< T > to_angle(const MatBase< T, 2, 2 > &mat)
QuaternionBase< float > Quaternion
T length_squared(const VecBase< T, Size > &a)
T cos(const AngleRadianBase< T > &a)
bool is_normalized(const DualQuaternionBase< T > &dq)
T sqrt(const T &a)
QuaternionBase< T > to_quaternion(const AxisAngleBase< T, AngleT > &axis_angle)
bool is_negative(const MatBase< T, 3, 3 > &mat)
MatBase< T, NumCol, NumRow > scale(const MatBase< T, NumCol, NumRow > &mat, const VectorT &scale)
bool is_orthonormal(const MatT &mat)
T sign(const T &a)
Quaternion to_quaternion_legacy(const float3x3 &mat)
MatT from_up_axis(const VectorT up)
T to_vector(const Axis axis)
T length(const VecBase< T, Size > &a)
MatT from_loc_rot(const typename MatT::loc_type &location, const RotationT &rotation)
MatT from_loc_scale(const typename MatT::loc_type &location, const VecBase< typename MatT::base_type, ScaleDim > &scale)
AxisAngleBase< float, AngleRadianBase< float > > AxisAngle
bool is_uniformly_scaled(const MatT &mat)
QuaternionBase< T > normalize_and_get_length(const QuaternionBase< T > &q, T &out_length)
bool is_unit_scale(const MatBase< T, NumCol, NumRow > &m)
MatBase< T, Size, Size > pseudo_invert(const MatBase< T, Size, Size > &mat, T epsilon=1e-8)
T length_manhattan(const VecBase< T, Size > &a)
EulerXYZBase< float > EulerXYZ
T dot(const QuaternionBase< T > &a, const QuaternionBase< T > &b)
bool is_zero(const T &a)
CartesianBasis invert(const CartesianBasis &basis)
AxisSigned cross(const AxisSigned a, const AxisSigned b)
VectorT project_point(const MatT &mat, const VectorT &point)
MatT from_scale(const VecBase< typename MatT::base_type, ScaleDim > &scale)
T interpolate(const T &a, const T &b, const FactorT &t)
void transform_normals(const float3x3 &transform, MutableSpan< float3 > normals)
EulerXYZBase< T > to_nearest_euler(const MatBase< T, 3, 3 > &mat, const EulerXYZBase< T > &reference)
MatBase< T, 2, 2 > from_direction(const VecBase< T, 2 > &direction)
T atan2(const T &y, const T &x)
MatBase< T, 3, 3 > interpolate_fast(const MatBase< T, 3, 3 > &a, const MatBase< T, 3, 3 > &b, T t)
MatBase< T, NumCol, NumRow > translate(const MatBase< T, NumCol, NumRow > &mat, const VectorT &translation)
QuaternionBase< T > normalized_to_quaternion_safe(const MatBase< T, 3, 3 > &mat)
MatBase< T, NumCol, NumRow > normalize_and_get_size(const MatBase< T, NumCol, NumRow > &a, VectorT &r_size)
MatT from_origin_transform(const MatT &transform, const VectorT origin)
MatBase< T, NumCol, NumRow > normalize(const MatBase< T, NumCol, NumRow > &a)
VecBase< T, 3 > to_scale(const MatBase< T, NumCol, NumRow > &mat)
AxisAngleBase< float, AngleCartesianBase< float > > AxisAngleCartesian
T square(const T &a)
bool is_equal(const MatBase< T, NumCol, NumRow > &a, const MatBase< T, NumCol, NumRow > &b, const T epsilon=T(0))
T sin(const AngleRadianBase< T > &a)
void to_rot_scale(const MatBase< T, 2, 2 > &mat, AngleRadianBase< T > &r_rotation, VecBase< T, 2 > &r_scale)
MatT orthogonalize(const MatT &mat, const Axis axis)
Euler3Base< float > Euler3
bool is_identity(const MatBase< T, NumCol, NumRow > &mat)
T abs(const T &a)
MatBase< T, NumCol, NumRow > interpolate_linear(const MatBase< T, NumCol, NumRow > &a, const MatBase< T, NumCol, NumRow > &b, T t)
T tan(const AngleRadianBase< T > &a)
void to_loc_rot_scale_safe(const MatBase< T, 4, 4 > &mat, VecBase< T, 3 > &r_location, RotationT &r_rotation, VecBase< T, 3 > &r_scale)
MatT from_location(const typename MatT::loc_type &location)
void to_loc_rot_scale(const MatBase< T, 3, 3 > &mat, VecBase< T, 2 > &r_location, AngleRadianBase< T > &r_rotation, VecBase< T, 2 > &r_scale)
VecBase< T, 3 > transform_direction(const MatBase< T, 3, 3 > &mat, const VecBase< T, 3 > &direction)
MatT from_rot_scale(const RotationT &rotation, const VectorT &scale)
MatBase< T, NumCol, NumRow > rotate(const MatBase< T, NumCol, NumRow > &mat, const RotationT &rotation)
CartesianBasis from_orthonormal_axes(const AxisSigned forward, const AxisSigned up)
MatT from_rotation(const RotationT &rotation)
T hypot(const T &y, const T &x)
T determinant(const MatBase< T, Size, Size > &mat)
Euler3Base< T > to_euler(const AxisAngleBase< T, AngleT > &axis_angle, EulerOrder order)
void transform_points(const float4x4 &transform, MutableSpan< float3 > points, bool use_threading=true)
MatT from_loc_rot_scale(const typename MatT::loc_type &location, const RotationT &rotation, const VecBase< typename MatT::base_type, ScaleDim > &scale)
VecBase< T, 3 > transform_point(const CartesianBasis &basis, const VecBase< T, 3 > &v)
bool is_orthogonal(const MatT &mat)
MatBase< float, 4, 4 > float4x4
MatBase< double, 3, 3 > double3x3
VecBase< float, 2 > float2
MatBase< float, 3, 3 > float3x3
void unroll(Fn fn)
Definition BLI_unroll.hh:25
VecBase< float, 3 > float3
const btScalar eps
Definition poly34.cpp:11
static MatBase identity()
Euler3Base wrapped_around(const Euler3Base &reference) const
const EulerOrder & order() const
EulerXYZBase wrapped_around(const EulerXYZBase &reference) const
i
Definition text_draw.cc:230