diff --git a/source/blender/blenlib/BLI_float3x3.hh b/source/blender/blenlib/BLI_float3x3.hh
new file mode 100644
index 00000000000..62478556d9b
--- /dev/null
+++ b/source/blender/blenlib/BLI_float3x3.hh
@@ -0,0 +1,192 @@
+/* SPDX-License-Identifier: GPL-2.0-or-later */
+
+#pragma once
+
+#include <cmath>
+#include <cstdint>
+
+#include "BLI_assert.h"
+#include "BLI_math_base.h"
+#include "BLI_math_matrix.h"
+#include "BLI_math_vec_types.hh"
+#include "BLI_math_vector.h"
+
+namespace blender {
+
+struct float3x3 {
+  /* A 3x3 matrix in column major order. */
+  float values[3][3];
+
+  float3x3() = default;
+
+  float3x3(const float *matrix)
+  {
+    memcpy(values, matrix, sizeof(float) * 3 * 3);
+  }
+
+  float3x3(const float matrix[3][3]) : float3x3(static_cast<const float *>(matrix[0]))
+  {
+  }
+
+  static float3x3 zero()
+  {
+    float3x3 result;
+    zero_m3(result.values);
+    return result;
+  }
+
+  static float3x3 identity()
+  {
+    float3x3 result;
+    unit_m3(result.values);
+    return result;
+  }
+
+  static float3x3 from_translation(const float2 translation)
+  {
+    float3x3 result = identity();
+    result.values[2][0] = translation.x;
+    result.values[2][1] = translation.y;
+    return result;
+  }
+
+  static float3x3 from_rotation(float rotation)
+  {
+    float3x3 result = zero();
+    const float cosine = std::cos(rotation);
+    const float sine = std::sin(rotation);
+    result.values[0][0] = cosine;
+    result.values[0][1] = sine;
+    result.values[1][0] = -sine;
+    result.values[1][1] = cosine;
+    result.values[2][2] = 1.0f;
+    return result;
+  }
+
+  static float3x3 from_translation_rotation_scale(const float2 translation,
+                                                  float rotation,
+                                                  const float2 scale)
+  {
+    float3x3 result;
+    const float cosine = std::cos(rotation);
+    const float sine = std::sin(rotation);
+    result.values[0][0] = scale.x * cosine;
+    result.values[0][1] = scale.x * sine;
+    result.values[0][2] = 0.0f;
+    result.values[1][0] = scale.y * -sine;
+    result.values[1][1] = scale.y * cosine;
+    result.values[1][2] = 0.0f;
+    result.values[2][0] = translation.x;
+    result.values[2][1] = translation.y;
+    result.values[2][2] = 1.0f;
+    return result;
+  }
+
+  static float3x3 from_normalized_axes(const float2 translation,
+                                       const float2 horizontal,
+                                       const float2 vertical)
+  {
+    BLI_ASSERT_UNIT_V2(horizontal);
+    BLI_ASSERT_UNIT_V2(vertical);
+
+    float3x3 result;
+    result.values[0][0] = horizontal.x;
+    result.values[0][1] = horizontal.y;
+    result.values[0][2] = 0.0f;
+    result.values[1][0] = vertical.x;
+    result.values[1][1] = vertical.y;
+    result.values[1][2] = 0.0f;
+    result.values[2][0] = translation.x;
+    result.values[2][1] = translation.y;
+    result.values[2][2] = 1.0f;
+    return result;
+  }
+
+  /* Construct a transformation that is pivoted around the given origin point. So for instance,
+   * from_origin_transformation(from_rotation(M_PI_2), float2(0.0f, 2.0f))
+   * will construct a transformation representing a 90 degree rotation around the point (0, 2). */
+  static float3x3 from_origin_transformation(const float3x3 &transformation, const float2 origin)
+  {
+    return from_translation(origin) * transformation * from_translation(-origin);
+  }
+
+  operator float *()
+  {
+    return &values[0][0];
+  }
+
+  operator const float *() const
+  {
+    return &values[0][0];
+  }
+
+  float *operator[](const int64_t index)
+  {
+    BLI_assert(index >= 0);
+    BLI_assert(index < 3);
+    return &values[index][0];
+  }
+
+  const float *operator[](const int64_t index) const
+  {
+    BLI_assert(index >= 0);
+    BLI_assert(index < 3);
+    return &values[index][0];
+  }
+
+  using c_style_float3x3 = float[3][3];
+  c_style_float3x3 &ptr()
+  {
+    return values;
+  }
+
+  const c_style_float3x3 &ptr() const
+  {
+    return values;
+  }
+
+  friend float3x3 operator*(const float3x3 &a, const float3x3 &b)
+  {
+    float3x3 result;
+    mul_m3_m3m3(result.values, a.values, b.values);
+    return result;
+  }
+
+  void operator*=(const float3x3 &other)
+  {
+    mul_m3_m3_post(values, other.values);
+  }
+
+  friend float2 operator*(const float3x3 &transformation, const float2 &vector)
+  {
+    float2 result;
+    mul_v2_m3v2(result, transformation.values, vector);
+    return result;
+  }
+
+  friend float2 operator*(const float3x3 &transformation, const float (*vector)[2])
+  {
+    return transformation * float2(vector);
+  }
+
+  float3x3 transposed() const
+  {
+    float3x3 result;
+    transpose_m3_m3(result.values, values);
+    return result;
+  }
+
+  float3x3 inverted() const
+  {
+    float3x3 result;
+    invert_m3_m3(result.values, values);
+    return result;
+  }
+
+  friend bool operator==(const float3x3 &a, const float3x3 &b)
+  {
+    return equals_m3m3(a.values, b.values);
+  }
+};
+
+}  // namespace blender
diff --git a/source/blender/blenlib/CMakeLists.txt b/source/blender/blenlib/CMakeLists.txt
index 4e76ae3e855..109230ebfa7 100644
--- a/source/blender/blenlib/CMakeLists.txt
+++ b/source/blender/blenlib/CMakeLists.txt
@@ -199,6 +199,7 @@ set(SRC
   BLI_fileops.hh
   BLI_fileops_types.h
   BLI_filereader.h
+  BLI_float3x3.hh
   BLI_float4x4.hh
   BLI_fnmatch.h
   BLI_function_ref.hh
@@ -431,6 +432,7 @@ if(WITH_GTESTS)
     tests/BLI_edgehash_test.cc
     tests/BLI_expr_pylike_eval_test.cc
     tests/BLI_fileops_test.cc
+    tests/BLI_float3x3_test.cc
     tests/BLI_function_ref_test.cc
     tests/BLI_generic_array_test.cc
     tests/BLI_generic_span_test.cc
diff --git a/source/blender/blenlib/tests/BLI_float3x3_test.cc b/source/blender/blenlib/tests/BLI_float3x3_test.cc
new file mode 100644
index 00000000000..d22993ee69e
--- /dev/null
+++ b/source/blender/blenlib/tests/BLI_float3x3_test.cc
@@ -0,0 +1,119 @@
+/* SPDX-License-Identifier: Apache-2.0 */
+
+#include "testing/testing.h"
+
+#include "BLI_float3x3.hh"
+#include "BLI_math_base.h"
+#include "BLI_math_vec_types.hh"
+
+namespace blender::tests {
+
+TEST(float3x3, Identity)
+{
+  float2 point(1.0f, 2.0f);
+  float3x3 transformation = float3x3::identity();
+  float2 result = transformation * point;
+  EXPECT_EQ(result, point);
+}
+
+TEST(float3x3, Translation)
+{
+  float2 point(1.0f, 2.0f);
+  float3x3 transformation = float3x3::from_translation(float2(5.0f, 3.0f));
+  float2 result = transformation * point;
+  EXPECT_FLOAT_EQ(result[0], 6.0f);
+  EXPECT_FLOAT_EQ(result[1], 5.0f);
+}
+
+TEST(float3x3, Rotation)
+{
+  float2 point(1.0f, 2.0f);
+  float3x3 transformation = float3x3::from_rotation(M_PI_2);
+  float2 result = transformation * point;
+  EXPECT_FLOAT_EQ(result[0], -2.0f);
+  EXPECT_FLOAT_EQ(result[1], 1.0f);
+}
+
+TEST(float3x3, TranslationRotationScale)
+{
+  float2 point(1.0f, 2.0f);
+  float3x3 transformation = float3x3::from_translation_rotation_scale(
+      float2(1.0f, 3.0f), M_PI_2, float2(2.0f, 3.0f));
+  float2 result = transformation * point;
+  EXPECT_FLOAT_EQ(result[0], -5.0f);
+  EXPECT_FLOAT_EQ(result[1], 5.0f);
+}
+
+TEST(float3x3, NormalizedAxes)
+{
+  float2 point(1.0f, 2.0f);
+
+  /* The horizontal is aligned with (1, 1) and vertical is aligned with (-1, 1), in other words, a
+   * Pi / 4 rotation. */
+  float value = std::sqrt(2.0f) / 2.0f;
+  float3x3 transformation = float3x3::from_normalized_axes(
+      float2(1.0f, 3.0f), float2(value), float2(-value, value));
+  float2 result = transformation * point;
+
+  float3x3 expected_transformation = float3x3::from_translation_rotation_scale(
+      float2(1.0f, 3.0f), M_PI_4, float2(1.0f));
+  float2 expected = expected_transformation * point;
+
+  EXPECT_FLOAT_EQ(result[0], expected[0]);
+  EXPECT_FLOAT_EQ(result[1], expected[1]);
+}
+
+TEST(float3x3, PostTransformationMultiplication)
+{
+  float2 point(1.0f, 2.0f);
+  float3x3 translation = float3x3::from_translation(float2(5.0f, 3.0f));
+  float3x3 rotation = float3x3::from_rotation(M_PI_2);
+  float3x3 transformation = translation * rotation;
+  float2 result = transformation * point;
+  EXPECT_FLOAT_EQ(result[0], 3.0f);
+  EXPECT_FLOAT_EQ(result[1], 4.0f);
+}
+
+TEST(float3x3, PreTransformationMultiplication)
+{
+  float2 point(1.0f, 2.0f);
+  float3x3 translation = float3x3::from_translation(float2(5.0f, 3.0f));
+  float3x3 rotation = float3x3::from_rotation(M_PI_2);
+  float3x3 transformation = rotation * translation;
+  float2 result = transformation * point;
+  EXPECT_FLOAT_EQ(result[0], -5.0f);
+  EXPECT_FLOAT_EQ(result[1], 6.0f);
+}
+
+TEST(float3x3, TransformationMultiplicationAssignment)
+{
+  float2 point(1.0f, 2.0f);
+  float3x3 transformation = float3x3::from_translation(float2(5.0f, 3.0f));
+  transformation *= float3x3::from_rotation(M_PI_2);
+  float2 result = transformation * point;
+  EXPECT_FLOAT_EQ(result[0], 3.0f);
+  EXPECT_FLOAT_EQ(result[1], 4.0f);
+}
+
+TEST(float3x3, Inverted)
+{
+  float2 point(1.0f, 2.0f);
+  float3x3 transformation = float3x3::from_translation_rotation_scale(
+      float2(1.0f, 3.0f), M_PI_4, float2(1.0f));
+  transformation *= transformation.inverted();
+  float2 result = transformation * point;
+  EXPECT_FLOAT_EQ(result[0], 1.0f);
+  EXPECT_FLOAT_EQ(result[1], 2.0f);
+}
+
+TEST(float3x3, Origin)
+{
+  float2 point(1.0f, 2.0f);
+  float3x3 rotation = float3x3::from_rotation(M_PI_2);
+  float3x3 transformation = float3x3::from_origin_transformation(rotation, float2(0.0f, 2.0f));
+  float2 result = transformation * point;
+  EXPECT_FLOAT_EQ(result[0], 0.0f);
+  EXPECT_FLOAT_EQ(result[1], 3.0f);
+}
+
+}  // namespace blender::tests