Add Mat44 affine reflection and shear constructors

lib/graphmath/Mat44.ex

@@ -1,8 +1,12 @@
 defmodule Graphmath.Mat44 do
   @moduledoc """
-  This is the 3D mathematics.
+  4x4 matrices and 3D affine transformations, using tuples of floats.
 
-  This submodule handles 4x4 matrices using tuples of floats.
+  Reflection and shear constructors act on three spatial coordinates and
+  preserve the fourth, homogeneous coordinate. Use `transform_point/2` for
+  points (`w = 1`) and `transform_vector/2` for directions (`w = 0`), so that
+  translation affects only points. A shear can change X, Y or Z; `w` is not
+  another spatial axis.
 
   Tuples store matrix rows in order. `apply(a, v)` computes the column-vector
   product **A****v**. Graphics constructors use row vectors: `transform_point/2`
@@ -251,6 +255,95 @@ defmodule Graphmath.Mat44 do
     {ct, st, 0.0, 0.0, -st, ct, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 1.0}
   end
 
+  @doc """
+  Creates a 3D affine reflection across `nx*x + ny*y + nz*z = offset`.
+
+  The normal may have any nonzero length. `offset` is the plane equation
+  constant; it is a signed distance only when the normal has unit length.
+  Scaling the normal and offset by the same nonzero factor leaves the mirror
+  unchanged. A zero normal raises `ArithmeticError`.
+
+  Use `transform_point/2` for points and `transform_vector/2` for directions.
+  The homogeneous coordinate is preserved; translation affects only points.
+  """
+  @spec make_reflect(vec3, float) :: mat44
+  def make_reflect({nx, ny, nz}, offset)
+      when is_float(nx) and is_float(ny) and is_float(nz) and is_float(offset) do
+    scale = max(abs(nx), max(abs(ny), abs(nz)))
+    ux = nx / scale
+    uy = ny / scale
+    uz = nz / scale
+    factor = 2.0 / (ux * ux + uy * uy + uz * uz)
+    scaled_offset = offset / scale
+
+    {1.0 - factor * ux * ux, -factor * ux * uy, -factor * ux * uz, 0.0, -factor * uy * ux,
+     1.0 - factor * uy * uy, -factor * uy * uz, 0.0, -factor * uz * ux, -factor * uz * uy,
+     1.0 - factor * uz * uz, 0.0, factor * scaled_offset * ux, factor * scaled_offset * uy,
+     factor * scaled_offset * uz, 1.0}
+  end
+
+  def make_reflect({nx, ny, nz}, offset) do
+    scale = max(abs(nx), max(abs(ny), abs(nz)))
+    ux = nx / scale
+    uy = ny / scale
+    uz = nz / scale
+    factor = 2.0 / (ux * ux + uy * uy + uz * uz)
+    scaled_offset = offset / scale
+
+    {1.0 - factor * ux * ux, -factor * ux * uy, -factor * ux * uz, 0.0, -factor * uy * ux,
+     1.0 - factor * uy * uy, -factor * uy * uz, 0.0, -factor * uz * ux, -factor * uz * uy,
+     1.0 - factor * uz * uz, 0.0, factor * scaled_offset * ux, factor * scaled_offset * uy,
+     factor * scaled_offset * uz, 1.0}
+  end
+
+  @doc """
+  Creates a 3D affine X shear: `x' = x + ky*y + kz*z`.
+
+  The other spatial coordinates and homogeneous coordinate are unchanged.
+  Use `transform_point/2` or `transform_vector/2` with 3 spatial components.
+  """
+  @spec make_shear_x(float, float) :: mat44
+  def make_shear_x(ky, kz)
+      when is_float(ky) and is_float(kz) do
+    {1.0, 0.0, 0.0, 0.0, ky, 1.0, 0.0, 0.0, kz, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 1.0}
+  end
+
+  def make_shear_x(ky, kz) do
+    {1.0, 0.0, 0.0, 0.0, 1.0 * ky, 1.0, 0.0, 0.0, 1.0 * kz, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 1.0}
+  end
+
+  @doc """
+  Creates a 3D affine Y shear: `y' = y + kx*x + kz*z`.
+
+  The other spatial coordinates and homogeneous coordinate are unchanged.
+  Use `transform_point/2` or `transform_vector/2` with 3 spatial components.
+  """
+  @spec make_shear_y(float, float) :: mat44
+  def make_shear_y(kx, kz)
+      when is_float(kx) and is_float(kz) do
+    {1.0, kx, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, kz, 1.0, 0.0, 0.0, 0.0, 0.0, 1.0}
+  end
+
+  def make_shear_y(kx, kz) do
+    {1.0, 1.0 * kx, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0 * kz, 1.0, 0.0, 0.0, 0.0, 0.0, 1.0}
+  end
+
+  @doc """
+  Creates a 3D affine Z shear: `z' = z + kx*x + ky*y`.
+
+  The other spatial coordinates and homogeneous coordinate are unchanged.
+  Use `transform_point/2` or `transform_vector/2` with 3 spatial components.
+  """
+  @spec make_shear_z(float, float) :: mat44
+  def make_shear_z(kx, ky)
+      when is_float(kx) and is_float(ky) do
+    {1.0, 0.0, kx, 0.0, 0.0, 1.0, ky, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 1.0}
+  end
+
+  def make_shear_z(kx, ky) do
+    {1.0, 0.0, 1.0 * kx, 0.0, 0.0, 1.0, 1.0 * ky, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 1.0}
+  end
+
   @doc """
   `round( a, sigfigs )` rounds every element of a `mat44` to some number of decimal places.
 
@@ -738,17 +831,11 @@ defmodule Graphmath.Mat44 do
       }
 
   @doc """
-  `transform_point( a, v )` transforms a `vec3` point by a `mat44`.
-
-  `a` is a `mat44` used to transform the point.
-
-  `v` is a `vec3` to be transformed.
-
-  This returns a `vec3` representing the application of `a` to `v`.
-
-  The point `a` is internally treated as having a fourth coordinate equal to 1.0.
+  Transforms a 3D point `{x, y, z}` using the row-vector product `{x, y, z, 1.0} a`.
 
-  Note that transforming a point will work for all transforms.
+  Returns the first three coordinates, including the effect of translation.
+  Use with 3D affine matrices. No perspective division is performed.
+  Use `apply_left/2` to supply and retain all four coordinates explicitly.
   """
   @spec transform_point(mat44, vec3) :: vec3
   def transform_point(
@@ -777,17 +864,11 @@ defmodule Graphmath.Mat44 do
       }
 
   @doc """
-  `transform_vector( a, v )` transforms a `vec3` vector by a `mat44`.
-
-  `a` is a `mat44` used to transform the point.
-
-  `v` is a `vec3` to be transformed.
-
-  This returns a `vec3` representing the application of `a` to `v`.
-
-  The point `a` is internally treated as having a fourth coordinate equal to 0.0.
+  Transforms a 3D direction `{x, y, z}` using the row-vector product `{x, y, z, 0.0} a`.
 
-  Note that transforming a vector will work for only rotations, scales, and shears.
+  Returns the first three coordinates. The zero homogeneous coordinate excludes
+  translation; rotations, scales, reflections and shears still affect the vector.
+  Use `apply_left/2` to supply and retain all four coordinates explicitly.
   """
   @spec transform_vector(mat44, vec3) :: vec3
   def transform_vector(

test/mat44/mat44_reflection_shear_test.exs

@@ -0,0 +1,209 @@
+defmodule GraphmathTest.Mat44.ReflectionShear do
+  use ExUnit.Case, async: true
+  alias Graphmath.{Mat44, Vec3}
+
+  @moduletag :mat44
+
+  @tag :make_reflect
+  test "an offset plane reflects points and directions with the same linear part" do
+    for {normal, offset} <- [
+          {{1.0, 2.0, 2.0}, 3.0},
+          {{1, 2, 2}, 3},
+          {{1.0, 2, 2.0}, 3.0},
+          {{1.0, 2.0, 2.0}, 3}
+        ] do
+      matrix = Mat44.make_reflect(normal, offset)
+
+      assert_close(matrix, {
+        7 / 9,
+        -4 / 9,
+        -4 / 9,
+        0,
+        -4 / 9,
+        1 / 9,
+        -8 / 9,
+        0,
+        -4 / 9,
+        -8 / 9,
+        1 / 9,
+        0,
+        2 / 3,
+        4 / 3,
+        4 / 3,
+        1
+      })
+
+      assert Enum.all?(Tuple.to_list(matrix), &is_float/1)
+      assert_close(Mat44.transform_point(matrix, {3, 4, 5}), {-1, -4, -3})
+      assert_close(Mat44.transform_vector(matrix, {3, 4, 5}), {-5 / 3, -16 / 3, -13 / 3})
+      assert_close(Mat44.transform_point(matrix, {3, 0, 0}), {3, 0, 0})
+      assert_close(Mat44.transform_vector(matrix, normal), {-1, -2, -2})
+      assert_close(Mat44.transform_vector(matrix, {2, -1, 0}), {2, -1, 0})
+    end
+  end
+
+  @tag :make_reflect
+  test "reflection preserves the fourth coordinate and weights translation by it" do
+    matrix = Mat44.make_reflect({1.0, 2.0, 2.0}, 3.0)
+    assert_close(Mat44.apply_left({3, 4, 5, 1}, matrix), {-1, -4, -3, 1})
+    assert_close(Mat44.apply_left({3, 4, 5, 0}, matrix), {-5 / 3, -16 / 3, -13 / 3, 0})
+    assert_close(Mat44.apply_left({3, 4, 5, 2.5}, matrix), {0, -2, -1, 2.5})
+  end
+
+  @tag :make_reflect
+  test "coordinate-plane reflections change only the normal component" do
+    for {normal, offset, expected} <- [
+          {{1.0, 0.0, 0.0}, 2.0, {-1, 3, 6}},
+          {{0.0, 1.0, 0.0}, -2.0, {5, -7, 6}},
+          {{0.0, 0.0, 1.0}, 2.0, {5, 3, -2}}
+        ] do
+      matrix = Mat44.make_reflect(normal, offset)
+      assert_close(Mat44.transform_point(matrix, {5, 3, 6}), expected)
+    end
+  end
+
+  @tag :make_reflect
+  test "normal and offset scale together without changing the plane" do
+    expected = Mat44.make_reflect({1.0, 2.0, 2.0}, 3.0)
+
+    for scale <- [2.0, -2.0, 1.0e-200, 1.0e200] do
+      assert_close(Mat44.make_reflect({scale, 2 * scale, 2 * scale}, 3 * scale), expected)
+    end
+  end
+
+  @tag :make_reflect
+  test "reflection reverses signed distance and preserves distances between points" do
+    normal = {1.0, 2.0, 2.0}
+    matrix = Mat44.make_reflect(normal, 3.0)
+    a = {3.0, 4.0, 5.0}
+    b = {-2.0, 1.0, -3.0}
+    reflected_a = Mat44.transform_point(matrix, a)
+    reflected_b = Mat44.transform_point(matrix, b)
+    assert_in_delta Vec3.dot(normal, reflected_a) - 3.0, -(Vec3.dot(normal, a) - 3.0), 1.0e-12
+    assert_in_delta Vec3.dot(normal, reflected_b) - 3.0, -(Vec3.dot(normal, b) - 3.0), 1.0e-12
+
+    assert_in_delta Vec3.length_squared(Vec3.subtract(reflected_a, reflected_b)),
+                    Vec3.length_squared(Vec3.subtract(a, b)),
+                    1.0e-12
+  end
+
+  @tag :make_reflect
+  test "reflection reverses orientation and is its own inverse" do
+    matrix = Mat44.make_reflect({1.0, 2.0, 2.0}, 3.0)
+    assert_in_delta Mat44.determinant(matrix), -1.0, 1.0e-12
+    assert_close(Mat44.multiply(matrix, matrix), Mat44.identity())
+    assert_close(Mat44.inverse(matrix), matrix)
+  end
+
+  for {plane, normals, expected, on_plane} <- [
+        {"YZ plane (X normal)", [{1.0, 0.0, 0.0}, {1, 0, 0}, {1.0, 0, 0.0}], {-3, 4, 5},
+         {0, -2, 7}},
+        {"XZ plane (Y normal)", [{0.0, 1.0, 0.0}, {0, 1, 0}, {0, 1.0, 0.0}], {3, -4, 5},
+         {2, 0, 7}},
+        {"XY plane (Z normal)", [{0.0, 0.0, 1.0}, {0, 0, 1}, {0.0, 0, 1.0}], {3, 4, -5},
+         {2, -3, 0}},
+        {"oblique plane", [{1.0, 2.0, 2.0}, {1, 2, 2}, {1.0, 2, 2.0}], {-5 / 3, -16 / 3, -13 / 3},
+         {2, -1, 0}}
+      ] do
+    @normals normals
+    @expected expected
+    @on_plane on_plane
+
+    @tag :make_reflect
+    test "zero offsets reflect across the #{plane} without translation" do
+      for normal <- @normals, offset <- [0.0, 0] do
+        matrix = Mat44.make_reflect(normal, offset)
+        assert_close(Mat44.row3(matrix), {0, 0, 0, 1})
+        assert_close(Mat44.transform_point(matrix, {0, 0, 0}), {0, 0, 0})
+        assert_close(Mat44.transform_point(matrix, @on_plane), @on_plane)
+        assert_close(Mat44.transform_point(matrix, {3, 4, 5}), @expected)
+        assert_close(Mat44.transform_vector(matrix, {3, 4, 5}), @expected)
+        {x, y, z} = @expected
+
+        for w <- [0.0, 1.0] do
+          assert_close(Mat44.apply_left({3, 4, 5, w}, matrix), {x, y, z, w})
+        end
+      end
+    end
+  end
+
+  @tag :make_reflect
+  test "a zero normal raises ArithmeticError" do
+    for normal <- [{0.0, 0.0, 0.0}, {0, 0, 0}, {0, 0.0, 0}], offset <- [0.0, 1] do
+      assert_raise ArithmeticError, fn -> Mat44.make_reflect(normal, offset) end
+    end
+  end
+
+  for {constructor, expected} <- [
+        make_shear_x: {-15, -3, 4},
+        make_shear_y: {2, -5, 4},
+        make_shear_z: {2, -3, 16}
+      ] do
+    @constructor constructor
+    @expected expected
+    @tag constructor
+    test "#{constructor} uses both spatial coefficients and preserves the fourth coordinate" do
+      for {a, b} <- [{3.0, -2.0}, {3, -2}, {3, -2.0}, {3.0, -2}] do
+        matrix = apply(Mat44, @constructor, [a, b])
+        assert Enum.all?(Tuple.to_list(matrix), &is_float/1)
+        assert_close(Mat44.transform_point(matrix, {2, -3, 4}), @expected)
+        assert_close(Mat44.transform_vector(matrix, {2, -3, 4}), @expected)
+        {x, y, z} = @expected
+
+        for w <- [0.0, 1.0, 2.5] do
+          assert_close(Mat44.apply_left({2, -3, 4, w}, matrix), {x, y, z, w})
+        end
+
+        assert Mat44.determinant(matrix) === 1.0
+        inverse = apply(Mat44, @constructor, [-a, -b])
+        assert_close(Mat44.multiply(matrix, inverse), Mat44.identity())
+        assert_close(Mat44.inverse(matrix), inverse)
+      end
+
+      for k <- [0, 0.0], do: assert(apply(Mat44, @constructor, [k, k]) === Mat44.identity())
+    end
+  end
+
+  @tag :make_shear_x
+  @tag :make_shear_y
+  @tag :make_shear_z
+  test "fractional coefficients retain their magnitude" do
+    assert_close(Mat44.transform_point(Mat44.make_shear_x(0.5, -0.25), {2, -3, 4}), {-0.5, -3, 4})
+    assert_close(Mat44.transform_point(Mat44.make_shear_y(0.5, -0.25), {2, -3, 4}), {2, -3, 4})
+    assert_close(Mat44.transform_point(Mat44.make_shear_z(0.5, -0.25), {2, -3, 4}), {2, -3, 5.75})
+  end
+
+  @tag :make_shear_x
+  @tag :make_shear_y
+  test "shear composition follows row-vector order" do
+    x = Mat44.make_shear_x(2.0, -1.0)
+    y = Mat44.make_shear_y(-3.0, 2.0)
+    assert_close(Mat44.transform_point(Mat44.multiply(x, y), {2, -3, 4}), {-8, 29, 4})
+    assert_close(Mat44.transform_point(Mat44.multiply(y, x), {2, -3, 4}), {-4, -1, 4})
+  end
+
+  @tag :make_shear_x
+  @tag :make_translate
+  test "composing shear and translation distinguishes points from directions" do
+    shear = Mat44.make_shear_x(3.0, -2.0)
+    translation = Mat44.make_translate(10.0, 2.0, -1.0)
+    composed = Mat44.multiply(shear, translation)
+    assert_close(Mat44.transform_point(composed, {2, -3, 4}), {-5, -1, 3})
+    assert_close(Mat44.transform_vector(composed, {2, -3, 4}), {-15, -3, 4})
+    assert_close(Mat44.apply_left({2, -3, 4, 1}, composed), {-5, -1, 3, 1})
+    assert_close(Mat44.apply_left({2, -3, 4, 0}, composed), {-15, -3, 4, 0})
+
+    assert_close(
+      Mat44.transform_point(Mat44.multiply(translation, shear), {2, -3, 4}),
+      {3, -1, 3}
+    )
+  end
+
+  defp assert_close(actual, expected) do
+    assert tuple_size(actual) == tuple_size(expected)
+
+    for {a, b} <- Enum.zip(Tuple.to_list(actual), Tuple.to_list(expected)) do
+      assert_in_delta a, b, 1.0e-12
+    end
+  end
+end