BLI_math_matrix_test.cc

Go to the documentation of this file.

/* SPDX-FileCopyrightText: 2023 Blender Authors
 *
 * SPDX-License-Identifier: Apache-2.0 */
 
#include "testing/testing.h"
 
#include "BLI_array.hh"
#include "BLI_math_base.h"
#include "BLI_math_matrix.h"
#include "BLI_math_matrix.hh"
#include "BLI_math_rotation.h"
#include "BLI_math_rotation.hh"
#include "BLI_rand.hh"
 
TEST(math_matrix, interp_m4_m4m4_regular)
{
  /* Test 4x4 matrix interpolation without singularity, i.e. without axis flip. */
 
  /* Transposed matrix, so that the code here is written in the same way as print_m4() outputs. */
  /* This matrix represents T=(0.1, 0.2, 0.3), R=(40, 50, 60) degrees, S=(0.7, 0.8, 0.9) */
  float matrix_a[4][4] = {
      {0.224976f, -0.333770f, 0.765074f, 0.100000f},
      {0.389669f, 0.647565f, 0.168130f, 0.200000f},
      {-0.536231f, 0.330541f, 0.443163f, 0.300000f},
      {0.000000f, 0.000000f, 0.000000f, 1.000000f},
  };
  transpose_m4(matrix_a);
 
  float matrix_i[4][4];
  unit_m4(matrix_i);
 
  float result[4][4];
  const float epsilon = 1e-6;
  interp_m4_m4m4(result, matrix_i, matrix_a, 0.0f);
  EXPECT_M4_NEAR(result, matrix_i, epsilon);
 
  interp_m4_m4m4(result, matrix_i, matrix_a, 1.0f);
  EXPECT_M4_NEAR(result, matrix_a, epsilon);
 
  /* This matrix is based on the current implementation of the code, and isn't guaranteed to be
   * correct. It's just consistent with the current implementation. */
  float matrix_halfway[4][4] = {
      {0.690643f, -0.253244f, 0.484996f, 0.050000f},
      {0.271924f, 0.852623f, 0.012348f, 0.100000f},
      {-0.414209f, 0.137484f, 0.816778f, 0.150000f},
      {0.000000f, 0.000000f, 0.000000f, 1.000000f},
  };
 
  transpose_m4(matrix_halfway);
  interp_m4_m4m4(result, matrix_i, matrix_a, 0.5f);
  EXPECT_M4_NEAR(result, matrix_halfway, epsilon);
}
 
TEST(math_matrix, interp_m3_m3m3_singularity)
{
  /* A singularity means that there is an axis mirror in the rotation component of the matrix.
   * This is reflected in its negative determinant.
   *
   * The interpolation of 4x4 matrices performs linear interpolation on the translation component,
   * and then uses the 3x3 interpolation function to handle rotation and scale. As a result, this
   * test for a singularity in the rotation matrix only needs to test the 3x3 case. */
 
  /* Transposed matrix, so that the code here is written in the same way as print_m4() outputs. */
  /* This matrix represents R=(4, 5, 6) degrees, S=(-1, 1, 1) */
  float matrix_a[3][3] = {
      {-0.990737f, -0.098227f, 0.093759f},
      {-0.104131f, 0.992735f, -0.060286f},
      {0.087156f, 0.069491f, 0.993768f},
  };
  transpose_m3(matrix_a);
  EXPECT_NEAR(-1.0f, determinant_m3_array(matrix_a), 1e-6);
 
  /* This matrix represents R=(0, 0, 0), S=(-1, 1, 1) */
  float matrix_b[3][3] = {
      {-1.0f, 0.0f, 0.0f},
      {0.0f, 1.0f, 0.0f},
      {0.0f, 0.0f, 1.0f},
  };
  transpose_m3(matrix_b);
 
  float result[3][3];
  interp_m3_m3m3(result, matrix_a, matrix_b, 0.0f);
  EXPECT_M3_NEAR(result, matrix_a, 1e-5);
 
  interp_m3_m3m3(result, matrix_a, matrix_b, 1.0f);
  EXPECT_M3_NEAR(result, matrix_b, 1e-5);
 
  interp_m3_m3m3(result, matrix_a, matrix_b, 0.5f);
  float expect[3][3] = {
      {-0.997681f, -0.049995f, 0.046186f},
      {-0.051473f, 0.998181f, -0.031385f},
      {0.044533f, 0.033689f, 0.998440f},
  };
  transpose_m3(expect);
  EXPECT_M3_NEAR(result, expect, 1e-5);
 
  /* Interpolating between a matrix with and without axis flip can cause it to go through a zero
   * point. The determinant det(A) of a matrix represents the change in volume; interpolating
   * between matrices with det(A)=-1 and det(B)=1 will have to go through a point where
   * det(result)=0, so where the volume becomes zero. */
  float matrix_i[3][3];
  unit_m3(matrix_i);
  zero_m3(expect);
  interp_m3_m3m3(result, matrix_a, matrix_i, 0.5f);
  EXPECT_NEAR(0.0f, determinant_m3_array(result), 1e-5);
  EXPECT_M3_NEAR(result, expect, 1e-5);
}
 
TEST(math_matrix, mul_m3_series)
{
  float matrix[3][3] = {
      {2.0f, 0.0f, 0.0f},
      {0.0f, 3.0f, 0.0f},
      {0.0f, 0.0f, 5.0f},
  };
  mul_m3_series(matrix, matrix, matrix, matrix);
  float const expect[3][3] = {
      {8.0f, 0.0f, 0.0f},
      {0.0f, 27.0f, 0.0f},
      {0.0f, 0.0f, 125.0f},
  };
  EXPECT_M3_NEAR(matrix, expect, 1e-5);
}
 
TEST(math_matrix, mul_m4_series)
{
  float matrix[4][4] = {
      {2.0f, 0.0f, 0.0f, 0.0f},
      {0.0f, 3.0f, 0.0f, 0.0f},
      {0.0f, 0.0f, 5.0f, 0.0f},
      {0.0f, 0.0f, 0.0f, 7.0f},
  };
  mul_m4_series(matrix, matrix, matrix, matrix);
  float const expect[4][4] = {
      {8.0f, 0.0f, 0.0f, 0.0f},
      {0.0f, 27.0f, 0.0f, 0.0f},
      {0.0f, 0.0f, 125.0f, 0.0f},
      {0.0f, 0.0f, 0.0f, 343.0f},
  };
  EXPECT_M4_NEAR(matrix, expect, 1e-5);
}
 
namespace blender::tests {
 
using namespace blender::math;
 
TEST(math_matrix, MatrixInverse)
{
  float3x3 mat = float3x3::diagonal(2);
  float3x3 inv = invert(mat);
  float3x3 expect = float3x3({0.5f, 0.0f, 0.0f}, {0.0f, 0.5f, 0.0f}, {0.0f, 0.0f, 0.5f});
  EXPECT_M3_NEAR(inv, expect, 1e-5f);
 
  bool success;
  float3x3 mat2 = float3x3::all(1);
  float3x3 inv2 = invert(mat2, success);
  float3x3 expect2 = float3x3::all(0);
  EXPECT_M3_NEAR(inv2, expect2, 1e-5f);
  EXPECT_FALSE(success);
}
 
TEST(math_matrix, MatrixPseudoInverse)
{
  float4x4 mat = transpose(float4x4({0.224976f, -0.333770f, 0.765074f, 0.100000f},
                                    {0.389669f, 0.647565f, 0.168130f, 0.200000f},
                                    {-0.536231f, 0.330541f, 0.443163f, 0.300000f},
                                    {0.000000f, 0.000000f, 0.000000f, 1.000000f}));
  float4x4 expect = transpose(float4x4({0.224976f, -0.333770f, 0.765074f, 0.100000f},
                                       {0.389669f, 0.647565f, 0.168130f, 0.200000f},
                                       {-0.536231f, 0.330541f, 0.443163f, 0.300000f},
                                       {0.000000f, 0.000000f, 0.000000f, 1.000000f}));
  /* MSVC 2019 has issues deducing the right template parameters, use an explicit template
   * instantiating to sidestep the issue. */
  float4x4 inv = pseudo_invert<float, 4>(mat);
  pseudoinverse_m4_m4(expect.ptr(), mat.ptr(), 1e-8f);
  EXPECT_M4_NEAR(inv, expect, 1e-5f);
 
  float4x4 mat2 = transpose(float4x4({0.000000f, -0.333770f, 0.765074f, 0.100000f},
                                     {0.000000f, 0.647565f, 0.168130f, 0.200000f},
                                     {0.000000f, 0.330541f, 0.443163f, 0.300000f},
                                     {0.000000f, 0.000000f, 0.000000f, 1.000000f}));
  float4x4 expect2 = transpose(float4x4({0.000000f, 0.000000f, 0.000000f, 0.000000f},
                                        {-0.51311f, 1.02638f, 0.496437f, -0.302896f},
                                        {0.952803f, 0.221885f, 0.527413f, -0.297881f},
                                        {-0.0275438f, -0.0477073f, 0.0656508f, 0.9926f}));
  /* MSVC 2019 has issues deducing the right template parameters, use an explicit template
   * instantiating to sidestep the issue. */
  float4x4 inv2 = pseudo_invert<float, 4>(mat2);
  EXPECT_M4_NEAR(inv2, expect2, 1e-5f);
}
 
TEST(math_matrix, MatrixDeterminant)
{
  float2x2 m2({1, 2}, {3, 4});
  float3x3 m3({1, 2, 3}, {-3, 4, -5}, {5, -6, 7});
  float4x4 m4({1, 2, -3, 3}, {3, 4, -5, 3}, {5, 6, 7, -3}, {5, 6, 7, 1});
  EXPECT_NEAR(determinant(m2), -2.0f, 1e-8f);
  EXPECT_NEAR(determinant(m3), -16.0f, 1e-8f);
  EXPECT_NEAR(determinant(m4), -112.0f, 1e-8f);
  EXPECT_NEAR(determinant(double2x2(m2)), -2.0f, 1e-8f);
  EXPECT_NEAR(determinant(double3x3(m3)), -16.0f, 1e-8f);
  EXPECT_NEAR(determinant(double4x4(m4)), -112.0f, 1e-8f);
}
 
TEST(math_matrix, MatrixAdjoint)
{
  float2x2 m2({1, 2}, {3, 4});
  float3x3 m3({1, 2, 3}, {-3, 4, -5}, {5, -6, 7});
  float4x4 m4({1, 2, -3, 3}, {3, 4, -5, 3}, {5, 6, 7, -3}, {5, 6, 7, 1});
  float2x2 expect2 = transpose(float2x2({4, -3}, {-2, 1}));
  float3x3 expect3 = transpose(float3x3({-2, -4, -2}, {-32, -8, 16}, {-22, -4, 10}));
  float4x4 expect4 = transpose(
      float4x4({232, -184, -8, -0}, {-128, 88, 16, 0}, {80, -76, 4, 28}, {-72, 60, -12, -28}));
  EXPECT_M2_NEAR(adjoint(m2), expect2, 1e-8f);
  EXPECT_M3_NEAR(adjoint(m3), expect3, 1e-8f);
  EXPECT_M4_NEAR(adjoint(m4), expect4, 1e-8f);
}
 
TEST(math_matrix, MatrixAccess)
{
  float4x4 m({1, 2, 3, 4}, {5, 6, 7, 8}, {9, 1, 2, 3}, {4, 5, 6, 7});
  EXPECT_EQ(m.x_axis(), float3(1, 2, 3));
  EXPECT_EQ(m.y_axis(), float3(5, 6, 7));
  EXPECT_EQ(m.z_axis(), float3(9, 1, 2));
  EXPECT_EQ(m.location(), float3(4, 5, 6));
}
 
TEST(math_matrix, MatrixInit)
{
  float4x4 expect;
 
  float4x4 m = from_location<float4x4>({1, 2, 3});
  expect = float4x4({1, 0, 0, 0}, {0, 1, 0, 0}, {0, 0, 1, 0}, {1, 2, 3, 1});
  EXPECT_TRUE(is_equal(m, expect, 0.00001f));
 
  expect = transpose(float4x4({0.411982, -0.833738, -0.36763, 0},
                              {-0.0587266, -0.426918, 0.902382, 0},
                              {-0.909297, -0.350175, -0.224845, 0},
                              {0, 0, 0, 1}));
  EulerXYZ euler(1, 2, 3);
  Quaternion quat = to_quaternion(euler);
  AxisAngle axis_angle = to_axis_angle(euler);
  m = from_rotation<float4x4>(euler);
  EXPECT_M3_NEAR(m, expect, 1e-5);
  m = from_rotation<float4x4>(quat);
  EXPECT_M3_NEAR(m, expect, 1e-5);
  m = from_rotation<float4x4>(axis_angle);
  EXPECT_M3_NEAR(m, expect, 1e-5);
 
  expect = transpose(float4x4({0.823964, -1.66748, -0.735261, 3.28334},
                              {-0.117453, -0.853835, 1.80476, 5.44925},
                              {-1.81859, -0.700351, -0.44969, -0.330972},
                              {0, 0, 0, 1}));
  DualQuaternion dual_quat(quat, Quaternion(0.5f, 0.5f, 0.5f, 1.5f), float4x4::diagonal(2.0f));
  m = from_rotation<float4x4>(dual_quat);
  EXPECT_M3_NEAR(m, expect, 1e-5);
 
  m = from_scale<float4x4>(float4(1, 2, 3, 4));
  expect = float4x4({1, 0, 0, 0}, {0, 2, 0, 0}, {0, 0, 3, 0}, {0, 0, 0, 4});
  EXPECT_TRUE(is_equal(m, expect, 0.00001f));
 
  m = from_scale<float4x4>(float3(1, 2, 3));
  expect = float4x4({1, 0, 0, 0}, {0, 2, 0, 0}, {0, 0, 3, 0}, {0, 0, 0, 1});
  EXPECT_TRUE(is_equal(m, expect, 0.00001f));
 
  m = from_scale<float4x4>(float2(1, 2));
  expect = float4x4({1, 0, 0, 0}, {0, 2, 0, 0}, {0, 0, 1, 0}, {0, 0, 0, 1});
  EXPECT_TRUE(is_equal(m, expect, 0.00001f));
 
  m = from_loc_rot<float4x4>({1, 2, 3}, EulerXYZ{1, 2, 3});
  expect = float4x4({0.411982, -0.0587266, -0.909297, 0},
                    {-0.833738, -0.426918, -0.350175, 0},
                    {-0.36763, 0.902382, -0.224845, 0},
                    {1, 2, 3, 1});
  EXPECT_TRUE(is_equal(m, expect, 0.00001f));
 
  m = from_loc_rot_scale<float4x4>({1, 2, 3}, EulerXYZ{1, 2, 3}, float3{1, 2, 3});
  expect = float4x4({0.411982, -0.0587266, -0.909297, 0},
                    {-1.66748, -0.853835, -0.700351, 0},
                    {-1.10289, 2.70714, -0.674535, 0},
                    {1, 2, 3, 1});
  EXPECT_TRUE(is_equal(m, expect, 0.00001f));
 
  float3 up = normalize(float3(1, 2, 3));
  m = from_up_axis<float4x4>(up);
  /* Output is not expected to be stable. Just test if they satisfy the expectations. */
  EXPECT_EQ(m.z_axis(), up);
  EXPECT_LE(abs(dot(m.z_axis(), m.x_axis())), 0.000001f);
  EXPECT_LE(abs(dot(m.y_axis(), m.x_axis())), 0.000001f);
  EXPECT_LE(abs(dot(m.z_axis(), m.y_axis())), 0.000001f);
  EXPECT_NEAR(1.0f, determinant(m), 1e-6);
}
 
TEST(math_matrix, MatrixModify)
{
  const float epsilon = 1e-6;
  float4x4 result, expect;
  float4x4 m1 = float4x4({0, 3, 0, 0}, {2, 0, 0, 0}, {0, 0, 2, 0}, {0, 0, 0, 1});
 
  expect = float4x4({0, 3, 0, 0}, {2, 0, 0, 0}, {0, 0, 2, 0}, {4, 9, 2, 1});
  result = translate(m1, float3(3, 2, 1));
  EXPECT_M4_NEAR(result, expect, epsilon);
 
  expect = float4x4({0, 3, 0, 0}, {2, 0, 0, 0}, {0, 0, 2, 0}, {4, 0, 0, 1});
  result = translate(m1, float2(0, 2));
  EXPECT_M4_NEAR(result, expect, epsilon);
 
  expect = float4x4({0, 0, -2, 0}, {2, 0, 0, 0}, {0, 3, 0, 0}, {0, 0, 0, 1});
  result = rotate(m1, AxisAngle({0, 1, 0}, M_PI_2));
  EXPECT_M4_NEAR(result, expect, epsilon);
 
  expect = float4x4({0, 9, 0, 0}, {4, 0, 0, 0}, {0, 0, 8, 0}, {0, 0, 0, 1});
  result = scale(m1, float3(3, 2, 4));
  EXPECT_M4_NEAR(result, expect, epsilon);
 
  expect = float4x4({0, 9, 0, 0}, {4, 0, 0, 0}, {0, 0, 2, 0}, {0, 0, 0, 1});
  result = scale(m1, float2(3, 2));
  EXPECT_M4_NEAR(result, expect, epsilon);
}
 
TEST(math_matrix, MatrixCompareTest)
{
  float4x4 m1 = float4x4({0, 3, 0, 0}, {2, 0, 0, 0}, {0, 0, 2, 0}, {0, 0, 0, 1});
  float4x4 m2 = float4x4({0, 3.001, 0, 0}, {1.999, 0, 0, 0}, {0, 0, 2.001, 0}, {0, 0, 0, 1.001});
  float4x4 m3 = float4x4({0, 3.001, 0, 0}, {1, 1, 0, 0}, {0, 0, 2.001, 0}, {0, 0, 0, 1.001});
  float4x4 m4 = float4x4({0, 1, 0, 0}, {1, 0, 0, 0}, {0, 0, 1, 0}, {0, 0, 0, 1});
  float4x4 m5 = float4x4({0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0});
  float4x4 m6 = float4x4({1, 0, 0, 0}, {0, 1, 0, 0}, {0, 0, 1, 0}, {0, 0, 0, 1});
  EXPECT_TRUE(is_equal(m1, m2, 0.01f));
  EXPECT_FALSE(is_equal(m1, m2, 0.0001f));
  EXPECT_FALSE(is_equal(m1, m3, 0.01f));
  EXPECT_TRUE(is_orthogonal(m1));
  EXPECT_FALSE(is_orthogonal(m3));
  EXPECT_TRUE(is_orthonormal(m4));
  EXPECT_FALSE(is_orthonormal(m1));
  EXPECT_FALSE(is_orthonormal(m3));
  EXPECT_FALSE(is_uniformly_scaled(m1));
  EXPECT_TRUE(is_uniformly_scaled(m4));
  EXPECT_FALSE(is_zero(m4));
  EXPECT_TRUE(is_zero(m5));
  EXPECT_TRUE(is_negative(m4));
  EXPECT_FALSE(is_negative(m5));
  EXPECT_FALSE(is_negative(m6));
}
 
TEST(math_matrix, MatrixMultiply)
{
 
  {
    const float4x4 matrix_a = {
        {1.0, 2.0, 3.0, 4.0},
        {5.0, 6.0, 7.0, 8.0},
        {9.0, 10.0, 11.0, 12.0},
        {13.0, 14.0, 15.0, 16.0},
    };
    const float4x4 matrix_b = {
        {0.1f, 0.2f, 0.3f, 0.4f},
        {0.5f, 0.6f, 0.7f, 0.8f},
        {0.9f, 1.0f, 1.1f, 1.2f},
        {1.3f, 1.4f, 1.5f, 1.6f},
    };
 
    const float4x4 expected = {
        {9.0f, 10.0f, 11.0f, 12.0f},
        {20.2f, 22.8f, 25.4f, 28.0f},
        {31.4f, 35.6f, 39.8f, 44.0f},
        {42.6f, 48.4f, 54.2f, 60.0f},
    };
 
    const float4x4 result = matrix_a * matrix_b;
 
    EXPECT_M4_NEAR(result, expected, 1e-5f);
  }
}
 
TEST(math_matrix, MatrixToNearestEuler)
{
  EulerXYZ eul1 = EulerXYZ(225.08542, -1.12485, -121.23738);
  Euler3 eul2 = Euler3(float3{4.06112, 100.561928, -18.9063}, EulerOrder::ZXY);
 
  float3x3 mat = {{0.808309, -0.578051, -0.111775},
                  {0.47251, 0.750174, -0.462572},
                  {0.351241, 0.321087, 0.879507}};
 
  EXPECT_V3_NEAR(float3(to_nearest_euler(mat, eul1)), float3(225.71, 0.112009, -120.001), 1e-3);
  EXPECT_V3_NEAR(float3(to_nearest_euler(mat, eul2)), float3(5.95631, 100.911, -19.5061), 1e-3);
}
 
TEST(math_matrix, MatrixMethods)
{
  float4x4 m = float4x4({0, 3, 0, 0}, {2, 0, 0, 0}, {0, 0, 2, 0}, {0, 1, 0, 1});
  auto expect_eul = EulerXYZ(0, 0, M_PI_2);
  auto expect_qt = Quaternion(0, -M_SQRT1_2, M_SQRT1_2, 0);
  float3 expect_scale = float3(3, 2, 2);
  float3 expect_location = float3(0, 1, 0);
 
  EXPECT_EQ(to_scale(m), expect_scale);
 
  float4 expect_size = {3, 2, 2, M_SQRT2};
  float4 size;
  float4x4 m1 = normalize_and_get_size(m, size);
  EXPECT_TRUE(is_unit_scale(m1));
  EXPECT_V4_NEAR(size, expect_size, 0.0002f);
 
  float4x4 m2 = normalize(m);
  EXPECT_TRUE(is_unit_scale(m2));
 
  EXPECT_V3_NEAR(float3(to_euler(m1)), float3(expect_eul), 0.0002f);
  EXPECT_V4_NEAR(float4(to_quaternion(m1)), float4(expect_qt), 0.0002f);
 
  EulerXYZ eul;
  Quaternion qt;
  float3 scale;
  to_rot_scale(float3x3(m), eul, scale);
  to_rot_scale(float3x3(m), qt, scale);
  EXPECT_V3_NEAR(scale, expect_scale, 0.00001f);
  EXPECT_V4_NEAR(float4(qt), float4(expect_qt), 0.0002f);
  EXPECT_V3_NEAR(float3(eul), float3(expect_eul), 0.0002f);
 
  float3 loc;
  to_loc_rot_scale(m, loc, eul, scale);
  to_loc_rot_scale(m, loc, qt, scale);
  EXPECT_V3_NEAR(scale, expect_scale, 0.00001f);
  EXPECT_V3_NEAR(loc, expect_location, 0.00001f);
  EXPECT_V4_NEAR(float4(qt), float4(expect_qt), 0.0002f);
  EXPECT_V3_NEAR(float3(eul), float3(expect_eul), 0.0002f);
}
 
TEST(math_matrix, Transformation2DMatrixDecomposition)
{
  const float2 translation = float2(1.0f, 2.0f);
  const AngleRadian rotation = AngleRadian(0.5f);
  const float2 scale = float2(5.0f, 3.0f);
 
  const float3x3 transformation = from_loc_rot_scale<float3x3>(translation, rotation, scale);
 
  AngleRadian decomposed_rotation;
  float2 decomposed_translation, decomposed_scale;
  to_loc_rot_scale(transformation, decomposed_translation, decomposed_rotation, decomposed_scale);
 
  EXPECT_V2_NEAR(decomposed_translation, translation, 0.00001f);
  EXPECT_V2_NEAR(decomposed_scale, scale, 0.00001f);
  EXPECT_NEAR(decomposed_rotation.radian(), rotation.radian(), 0.00001f);
}
 
TEST(math_matrix, MatrixToQuaternionLegacy)
{
  float3x3 mat = {{0.808309, -0.578051, -0.111775},
                  {0.47251, 0.750174, -0.462572},
                  {0.351241, 0.321087, 0.879507}};
 
  EXPECT_V4_NEAR(float4(to_quaternion_legacy(mat)),
                 float4(0.927091f, -0.211322f, 0.124857f, -0.283295f),
                 1e-5f);
}
 
TEST(math_matrix, MatrixTranspose)
{
  float4x4 m({1, 2, 3, 4}, {5, 6, 7, 8}, {9, 1, 2, 3}, {2, 5, 6, 7});
  float4x4 expect({1, 5, 9, 2}, {2, 6, 1, 5}, {3, 7, 2, 6}, {4, 8, 3, 7});
  EXPECT_EQ(transpose(m), expect);
}
 
TEST(math_matrix, MatrixInterpolationRegular)
{
  /* Test 4x4 matrix interpolation without singularity, i.e. without axis flip. */
 
  /* Transposed matrix, so that the code here is written in the same way as print_m4() outputs. */
  /* This matrix represents T=(0.1, 0.2, 0.3), R=(40, 50, 60) degrees, S=(0.7, 0.8, 0.9) */
  float4x4 m2 = transpose(float4x4({0.224976f, -0.333770f, 0.765074f, 0.100000f},
                                   {0.389669f, 0.647565f, 0.168130f, 0.200000f},
                                   {-0.536231f, 0.330541f, 0.443163f, 0.300000f},
                                   {0.000000f, 0.000000f, 0.000000f, 1.000000f}));
  float4x4 m1 = float4x4::identity();
  float4x4 result;
  const float epsilon = 1e-6;
  result = interpolate(m1, m2, 0.0f);
  EXPECT_M4_NEAR(result, m1, epsilon);
  result = interpolate(m1, m2, 1.0f);
  EXPECT_M4_NEAR(result, m2, epsilon);
 
  /* This matrix is based on the current implementation of the code, and isn't guaranteed to be
   * correct. It's just consistent with the current implementation. */
  float4x4 expect = transpose(float4x4({0.690643f, -0.253244f, 0.484996f, 0.050000f},
                                       {0.271924f, 0.852623f, 0.012348f, 0.100000f},
                                       {-0.414209f, 0.137484f, 0.816778f, 0.150000f},
                                       {0.000000f, 0.000000f, 0.000000f, 1.000000f}));
  result = interpolate(m1, m2, 0.5f);
  EXPECT_M4_NEAR(result, expect, epsilon);
 
  result = interpolate_fast(m1, m2, 0.5f);
  EXPECT_M4_NEAR(result, expect, epsilon);
}
 
TEST(math_matrix, MatrixInterpolationSingularity)
{
  /* A singularity means that there is an axis mirror in the rotation component of the matrix.
   * This is reflected in its negative determinant.
   *
   * The interpolation of 4x4 matrices performs linear interpolation on the translation component,
   * and then uses the 3x3 interpolation function to handle rotation and scale. As a result, this
   * test for a singularity in the rotation matrix only needs to test the 3x3 case. */
 
  /* Transposed matrix, so that the code here is written in the same way as print_m4() outputs. */
  /* This matrix represents R=(4, 5, 6) degrees, S=(-1, 1, 1) */
  float3x3 matrix_a = transpose(float3x3({-0.990737f, -0.098227f, 0.093759f},
                                         {-0.104131f, 0.992735f, -0.060286f},
                                         {0.087156f, 0.069491f, 0.993768f}));
  EXPECT_NEAR(-1.0f, determinant(matrix_a), 1e-6);
 
  /* This matrix represents R=(0, 0, 0), S=(-1, 1 1) */
  float3x3 matrix_b = transpose(
      float3x3({-1.0f, 0.0f, 0.0f}, {0.0f, 1.0f, 0.0f}, {0.0f, 0.0f, 1.0f}));
 
  float3x3 result = interpolate(matrix_a, matrix_b, 0.0f);
  EXPECT_M3_NEAR(result, matrix_a, 1e-5);
 
  result = interpolate(matrix_a, matrix_b, 1.0f);
  EXPECT_M3_NEAR(result, matrix_b, 1e-5);
 
  result = interpolate(matrix_a, matrix_b, 0.5f);
 
  float3x3 expect = transpose(float3x3({-0.997681f, -0.049995f, 0.046186f},
                                       {-0.051473f, 0.998181f, -0.031385f},
                                       {0.044533f, 0.033689f, 0.998440f}));
  EXPECT_M3_NEAR(result, expect, 1e-5);
 
  result = interpolate_fast(matrix_a, matrix_b, 0.5f);
  EXPECT_M3_NEAR(result, expect, 1e-5);
 
  /* Interpolating between a matrix with and without axis flip can cause it to go through a zero
   * point. The determinant det(A) of a matrix represents the change in volume; interpolating
   * between matrices with det(A)=-1 and det(B)=1 will have to go through a point where
   * det(result)=0, so where the volume becomes zero. */
  float3x3 matrix_i = float3x3::identity();
  expect = float3x3::zero();
  result = interpolate(matrix_a, matrix_i, 0.5f);
  EXPECT_NEAR(0.0f, determinant(result), 1e-5);
  EXPECT_M3_NEAR(result, expect, 1e-5);
}
 
TEST(math_matrix, MatrixTransform)
{
  float3 expect, result;
  const float3 p(1, 2, 3);
  float4x4 m4 = from_loc_rot<float4x4>({10, 0, 0}, EulerXYZ(M_PI_2, M_PI_2, M_PI_2));
  float3x3 m3 = from_rotation<float3x3>(EulerXYZ(M_PI_2, M_PI_2, M_PI_2));
  float4x4 pers4 = projection::perspective(-0.1f, 0.1f, -0.1f, 0.1f, -0.1f, -1.0f);
  float3x3 pers3 = float3x3({1, 0, 0.1f}, {0, 1, 0.1f}, {0, 0.1f, 1});
 
  expect = {13, 2, -1};
  result = transform_point(m4, p);
  EXPECT_V3_NEAR(result, expect, 1e-2);
 
  expect = {3, 2, -1};
  result = transform_point(m3, p);
  EXPECT_V3_NEAR(result, expect, 1e-5);
 
  result = transform_direction(m4, p);
  EXPECT_V3_NEAR(result, expect, 1e-5);
 
  result = transform_direction(m3, p);
  EXPECT_V3_NEAR(result, expect, 1e-5);
 
  expect = {-0.333333, -0.666666, -1.14814};
  result = project_point(pers4, p);
  EXPECT_V3_NEAR(result, expect, 1e-5);
 
  float2 expect2 = {0.76923, 1.61538};
  float2 result2 = project_point(pers3, float2(p));
  EXPECT_V2_NEAR(result2, expect2, 1e-5);
}
 
TEST(math_matrix, MatrixTransform2D)
{
  const float2 sample_point = float2(2.0f, 3.0f);
  const float3x3 transformation = math::from_loc_rot_scale<float3x3>(
      float2(2.0f, 0.5f), AngleRadian(M_PI_2), float2(1.5f, 1.0f));
 
  EXPECT_V2_NEAR(transform_point(transformation, sample_point), float2(-1.0f, 3.5f), 1e-3f);
 
  const float2x2 transformation_2d = float2x2(transformation);
  EXPECT_V2_NEAR(transform_point(transformation_2d, sample_point), float2(-3.0f, 3.0f), 1e-3f);
}
 
TEST(math_matrix, MatrixProjection)
{
  using namespace math::projection;
  float4x4 expect;
  float4x4 ortho = orthographic(-0.2f, 0.3f, -0.2f, 0.4f, -0.2f, -0.5f);
  float4x4 pers1 = perspective(-0.2f, 0.3f, -0.2f, 0.4f, -0.2f, -0.5f);
  float4x4 pers2 = perspective_fov(
      math::atan(-0.2f), math::atan(0.3f), math::atan(-0.2f), math::atan(0.4f), -0.2f, -0.5f);
 
  expect = transpose(float4x4({4.0f, 0.0f, 0.0f, -0.2f},
                              {0.0f, 3.33333f, 0.0f, -0.333333f},
                              {0.0f, 0.0f, 6.66667f, -2.33333f},
                              {0.0f, 0.0f, 0.0f, 1.0f}));
  EXPECT_M4_NEAR(ortho, expect, 1e-5);
 
  expect = transpose(float4x4({-0.8f, 0.0f, 0.2f, 0.0f},
                              {0.0f, -0.666667f, 0.333333f, 0.0f},
                              {0.0f, 0.0f, -2.33333f, 0.666667f},
                              {0.0f, 0.0f, -1.0f, 0.0f}));
  EXPECT_M4_NEAR(pers1, expect, 1e-5);
 
  expect = transpose(float4x4({4.0f, 0.0f, 0.2f, 0.0f},
                              {0.0f, 3.33333f, 0.333333f, 0.0f},
                              {0.0f, 0.0f, -2.33333f, 0.666667f},
                              {0.0f, 0.0f, -1.0f, 0.0f}));
  EXPECT_M4_NEAR(pers2, expect, 1e-5);
}
 
TEST(math_matrix, ToQuaternionSafe)
{
  float3x3 mat;
  mat[0] = {0.493316412f, -0.0f, 0.869849861f};
  mat[1] = {-0.0f, 1.0f, 0.0f};
  mat[2] = {-0.0176299568f, -0.0f, 0.999844611f};
 
  float3x3 expect;
  expect[0] = {0.493316f, 0.000000f, 0.869850f};
  expect[1] = {-0.000000f, 1.000000f, 0.000000f};
  expect[2] = {-0.869850f, -0.000000f, 0.493316f};
 
  /* This is mainly testing if there are any asserts hit because the matrix has shearing. */
  Quaternion rotation = math::normalized_to_quaternion_safe(normalize(mat));
  EXPECT_M3_NEAR(from_rotation<float3x3>(rotation), expect, 1e-5);
}
 
static void transform_normals_test(const float3x3 &transform, const Span<float3> normals)
{
  Array<float3, 10> transformed_normals(normals.size());
  transform_normals(normals, transform, transformed_normals);
  /* Input unit vectors should still be unit length. */
  for (const int64_t i : transformed_normals.index_range()) {
    EXPECT_NEAR(math::length(transformed_normals[i]), 1.0f, 1e-6f);
  }
  /* Make sure any optimization inside #transform_normals that skips the final vector normalization
   * doesn't change the result. */
  const float3x3 normal_transform = math::transpose(math::invert(transform));
  for (const int64_t i : normals.index_range()) {
    const float3 transformed_then_normalized = math::normalize(normal_transform * normals[i]);
    EXPECT_V3_NEAR(transformed_normals[i], transformed_then_normalized, 1e-6f);
  }
}
 
TEST(math_matrix, TransformNormals)
{
  RandomNumberGenerator rng;
  Array<float3, 10> normals(10);
  for (const int64_t i : normals.index_range()) {
    normals[i] = rng.get_unit_float3();
  }
  for (const int64_t i : normals.index_range()) {
    EXPECT_NEAR(math::length(normals[i]), 1.0f, 1e-6f);
  }
 
  transform_normals_test(float3x3::identity(), normals);
  transform_normals_test(from_rotation<float3x3>(EulerXYZ(4.2f, 0.2f, -3.2f)), normals);
  transform_normals_test(float3x3::diagonal(0.5f), normals);
  transform_normals_test(float3x3::diagonal(10.5f), normals);
  transform_normals_test(from_scale<float3x3>(float3(0.5f, 1.5f, 5.7f)), normals);
  transform_normals_test(from_rot_scale<float3x3>(EulerXYZ(0.0f, 7.1f, 3.14f), float3(2.15f)),
                         normals);
  transform_normals_test(
      from_rot_scale<float3x3>(EulerXYZ(0.0f, 7.1f, 3.14f), float3(0.5f, 1.5f, 5.7f)), normals);
}
 
}  // namespace blender::tests