diff --git a/lib/graphmath/Mat33.ex b/lib/graphmath/Mat33.ex index eaf894d..2a24896 100644 --- a/lib/graphmath/Mat33.ex +++ b/lib/graphmath/Mat33.ex @@ -1,8 +1,17 @@ defmodule Graphmath.Mat33 do @moduledoc """ - This is the 3D mathematics. + 3x3 matrices for 3D linear algebra and 2D affine transformations. - This submodule handles 3x3 matrices using tuples of floats. + Matrices use tuples of floats. Reflection and shear constructors explicitly + name the spatial dimension: + + * `_2d` constructors act on `{x, y, w}`, where `w` is a homogeneous coordinate. + Use `transform_point/2` for points (`w = 1`) and `transform_vector/2` for + directions (`w = 0`). + * `_3d` constructors act on full spatial vectors `{x, y, z}`. Use + `apply_left/2`; all three coordinates participate in the linear transform. + + For 3D affine transformations, including translation, use `Graphmath.Mat44`. Tuples store matrix rows in order: `{a11,a12,a13,a21,a22,a23,a31,a32,a33}`. `apply(a, v)` computes the column-vector product **A****v**. Graphics @@ -186,6 +195,154 @@ defmodule Graphmath.Mat33 do {ct, st, 0.0, -st, ct, 0.0, 0.0, 0.0, 1.0} end + @doc """ + Creates a 2D affine reflection across `nx*x + ny*y = offset`. + + The normal may have any nonzero length. `offset` is the line 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_2d(vec2, float) :: mat33 + def make_reflect_2d({nx, ny}, offset) + when is_float(nx) and is_float(ny) and is_float(offset) do + scale = max(abs(nx), abs(ny)) + ux = nx / scale + uy = ny / scale + factor = 2.0 / (ux * ux + uy * uy) + scaled_offset = offset / scale + + {1.0 - factor * ux * ux, -factor * ux * uy, 0.0, -factor * uy * ux, 1.0 - factor * uy * uy, + 0.0, factor * scaled_offset * ux, factor * scaled_offset * uy, 1.0} + end + + def make_reflect_2d({nx, ny}, offset) do + scale = max(abs(nx), abs(ny)) + ux = nx / scale + uy = ny / scale + factor = 2.0 / (ux * ux + uy * uy) + scaled_offset = offset / scale + + {1.0 - factor * ux * ux, -factor * ux * uy, 0.0, -factor * uy * ux, 1.0 - factor * uy * uy, + 0.0, factor * scaled_offset * ux, factor * scaled_offset * uy, 1.0} + end + + @doc """ + Creates a 2D affine X shear: `x' = x + k*y`. + + The other spatial coordinates and homogeneous coordinate are unchanged. + Use `transform_point/2` or `transform_vector/2` with 2 spatial components. + """ + @spec make_shear_x_2d(float) :: mat33 + def make_shear_x_2d(k) + when is_float(k) do + {1.0, 0.0, 0.0, k, 1.0, 0.0, 0.0, 0.0, 1.0} + end + + def make_shear_x_2d(k) do + {1.0, 0.0, 0.0, 1.0 * k, 1.0, 0.0, 0.0, 0.0, 1.0} + end + + @doc """ + Creates a 2D affine Y shear: `y' = y + k*x`. + + The other spatial coordinates and homogeneous coordinate are unchanged. + Use `transform_point/2` or `transform_vector/2` with 2 spatial components. + """ + @spec make_shear_y_2d(float) :: mat33 + def make_shear_y_2d(k) + when is_float(k) do + {1.0, k, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0} + end + + def make_shear_y_2d(k) do + {1.0, 1.0 * k, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0} + end + + @doc """ + Creates a 3D linear reflection across the plane through the origin + with normal `{nx, ny, nz}`. + + The normal may have any nonzero length; a zero normal raises `ArithmeticError`. + Apply using `apply_left/2` with a full 3-component vector. + """ + @spec make_reflect_3d(vec3) :: mat33 + def make_reflect_3d({nx, ny, nz}) + when is_float(nx) and is_float(ny) and is_float(nz) 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) + + {1.0 - factor * ux * ux, -factor * ux * uy, -factor * ux * uz, -factor * uy * ux, + 1.0 - factor * uy * uy, -factor * uy * uz, -factor * uz * ux, -factor * uz * uy, + 1.0 - factor * uz * uz} + end + + def make_reflect_3d({nx, ny, nz}) 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) + + {1.0 - factor * ux * ux, -factor * ux * uy, -factor * ux * uz, -factor * uy * ux, + 1.0 - factor * uy * uy, -factor * uy * uz, -factor * uz * ux, -factor * uz * uy, + 1.0 - factor * uz * uz} + end + + @doc """ + Creates a 3D linear X shear: `x' = x + ky*y + kz*z`. + + The other spatial coordinates are unchanged. + Apply using `apply_left/2` with a full 3-component vector. + """ + @spec make_shear_x_3d(float, float) :: mat33 + def make_shear_x_3d(ky, kz) + when is_float(ky) and is_float(kz) do + {1.0, 0.0, 0.0, ky, 1.0, 0.0, kz, 0.0, 1.0} + end + + def make_shear_x_3d(ky, kz) do + {1.0, 0.0, 0.0, 1.0 * ky, 1.0, 0.0, 1.0 * kz, 0.0, 1.0} + end + + @doc """ + Creates a 3D linear Y shear: `y' = y + kx*x + kz*z`. + + The other spatial coordinates are unchanged. + Apply using `apply_left/2` with a full 3-component vector. + """ + @spec make_shear_y_3d(float, float) :: mat33 + def make_shear_y_3d(kx, kz) + when is_float(kx) and is_float(kz) do + {1.0, kx, 0.0, 0.0, 1.0, 0.0, 0.0, kz, 1.0} + end + + def make_shear_y_3d(kx, kz) do + {1.0, 1.0 * kx, 0.0, 0.0, 1.0, 0.0, 0.0, 1.0 * kz, 1.0} + end + + @doc """ + Creates a 3D linear Z shear: `z' = z + kx*x + ky*y`. + + The other spatial coordinates are unchanged. + Apply using `apply_left/2` with a full 3-component vector. + """ + @spec make_shear_z_3d(float, float) :: mat33 + def make_shear_z_3d(kx, ky) + when is_float(kx) and is_float(ky) do + {1.0, 0.0, kx, 0.0, 1.0, ky, 0.0, 0.0, 1.0} + end + + def make_shear_z_3d(kx, ky) do + {1.0, 0.0, 1.0 * kx, 0.0, 1.0, 1.0 * ky, 0.0, 0.0, 1.0} + end + @doc """ `round( a, sigfigs )` rounds every element of a `mat33` to some number of decimal places. @@ -558,17 +715,11 @@ defmodule Graphmath.Mat33 do } @doc """ - `transform_point( a, v )` transforms a `vec2` point by a `mat33`. + Transforms a 2D point `{x, y}` using the row-vector product `{x, y, 1.0} a`. - `a` is a `mat33` used to transform the point. - - `v` is a `vec2` to be transformed. - - This returns a `vec2` representing the application of `a` to `v`. - - The point `a` is internally treated as having a third coordinate equal to 1.0. - - Note that transforming a point will work for all transforms. + Returns the first two coordinates, including the effect of translation. + Use with 2D affine matrices. No perspective division is performed. + For a full 3D linear transformation, use `apply_left/2` with `{x, y, z}`. """ @spec transform_point(mat33, vec2) :: vec2 def transform_point({a11, a21, _, a12, a22, _, a13, a23, _}, {x, y}) @@ -586,17 +737,11 @@ defmodule Graphmath.Mat33 do } @doc """ - `transform_vector( a, v )` transforms a `vec2` vector by a `mat33`. - - `a` is a `mat33` used to transform the point. - - `v` is a `vec2` to be transformed. - - This returns a `vec2` representing the application of `a` to `v`. - - The point `a` is internally treated as having a third coordinate equal to 0.0. + Transforms a 2D direction `{x, y}` using the row-vector product `{x, y, 0.0} a`. - Note that transforming a vector will work for only rotations, scales, and shears. + Returns the first two coordinates. The zero homogeneous coordinate excludes + translation; rotations, scales, reflections and shears still affect the vector. + For a full 3D linear transformation, use `apply_left/2` with `{x, y, z}`. """ @spec transform_vector(mat33, vec2) :: vec2 def transform_vector({a11, a21, _, a12, a22, _, _, _, _}, {x, y}) diff --git a/test/mat33/mat33_reflection_shear_test.exs b/test/mat33/mat33_reflection_shear_test.exs new file mode 100644 index 0000000..ef31793 --- /dev/null +++ b/test/mat33/mat33_reflection_shear_test.exs @@ -0,0 +1,185 @@ +defmodule GraphmathTest.Mat33.ReflectionShear do + use ExUnit.Case, async: true + alias Graphmath.{Mat33, Vec3} + + @moduletag :mat33 + + @tag :make_reflect_2d + test "an offset 2D mirror translates points but not directions" do + for {normal, offset} <- [ + {{1.0, -1.0}, 2.0}, + {{1, -1}, 2}, + {{1.0, -1}, 2.0}, + {{1.0, -1.0}, 2} + ] do + matrix = Mat33.make_reflect_2d(normal, offset) + assert_close(matrix, {0, 1, 0, 1, 0, 0, 2, -2, 1}) + assert Enum.all?(Tuple.to_list(matrix), &is_float/1) + assert_close(Mat33.transform_point(matrix, {3, 4}), {6, 1}) + assert_close(Mat33.transform_vector(matrix, {3, 4}), {4, 3}) + assert_close(Mat33.transform_point(matrix, {4, 2}), {4, 2}) + assert_close(Mat33.apply_left({3, 4, 1}, matrix), {6, 1, 1}) + assert_close(Mat33.apply_left({3, 4, 0}, matrix), {4, 3, 0}) + end + end + + @tag :make_reflect_2d + test "2D reflection handles coordinate lines and a line through the origin" do + x = Mat33.make_reflect_2d({1.0, 0.0}, 2.0) + y = Mat33.make_reflect_2d({0.0, 1.0}, -2.0) + origin = Mat33.make_reflect_2d({1.0, -1.0}, 0.0) + assert_close(Mat33.transform_point(x, {5, 3}), {-1, 3}) + assert_close(Mat33.transform_point(y, {5, 3}), {5, -7}) + assert_close(Mat33.transform_point(origin, {5, 3}), {3, 5}) + assert_close(Mat33.transform_vector(origin, {5, 3}), {3, 5}) + end + + @tag :make_reflect_2d + test "the 2D offset is the line equation constant for a non-unit normal" do + matrix = Mat33.make_reflect_2d({3.0, 4.0}, 10.0) + assert_close(matrix, {7 / 25, -24 / 25, 0, -24 / 25, -7 / 25, 0, 12 / 5, 16 / 5, 1}) + assert_close(Mat33.transform_point(matrix, {2, 1}), {2, 1}) + assert_close(Mat33.transform_point(matrix, {5, 0}), {19 / 5, -8 / 5}) + assert_close(Mat33.transform_vector(matrix, {3, 4}), {-3, -4}) + assert_close(Mat33.multiply(matrix, matrix), Mat33.identity()) + assert_in_delta Mat33.determinant(matrix), -1.0, 1.0e-12 + + for scale <- [2.0, -2.0, 1.0e-200, 1.0e200] do + assert_close(Mat33.make_reflect_2d({3 * scale, 4 * scale}, 10 * scale), matrix) + end + end + + @tag :make_reflect_3d + test "3D reflection uses all three spatial coordinates" do + for normal <- [{1.0, 2.0, 2.0}, {1, 2, 2}, {1.0, 2, 2.0}] do + matrix = Mat33.make_reflect_3d(normal) + assert_close(matrix, {7 / 9, -4 / 9, -4 / 9, -4 / 9, 1 / 9, -8 / 9, -4 / 9, -8 / 9, 1 / 9}) + assert Enum.all?(Tuple.to_list(matrix), &is_float/1) + assert_close(Mat33.apply_left({3, 4, 5}, matrix), {-5 / 3, -16 / 3, -13 / 3}) + assert_close(Mat33.apply_left(normal, matrix), {-1, -2, -2}) + assert_close(Mat33.apply_left({2, -1, 0}, matrix), {2, -1, 0}) + end + end + + @tag :make_reflect_3d + test "coordinate-plane mirrors change only the normal component" do + for {normal, expected} <- [ + {{1.0, 0.0, 0.0}, {-3, 4, 5}}, + {{0.0, 1.0, 0.0}, {3, -4, 5}}, + {{0.0, 0.0, 1.0}, {3, 4, -5}} + ] do + assert_close(Mat33.apply_left({3, 4, 5}, Mat33.make_reflect_3d(normal)), expected) + end + end + + @tag :make_reflect_3d + test "3D reflection preserves length, reverses orientation and is its own inverse" do + matrix = Mat33.make_reflect_3d({1.0, 2.0, 2.0}) + assert_in_delta Vec3.length_squared(Mat33.apply_left({3, 4, 5}, matrix)), 50.0, 1.0e-12 + assert_in_delta Mat33.determinant(matrix), -1.0, 1.0e-12 + assert_close(Mat33.multiply(matrix, matrix), Mat33.identity()) + assert_close(Mat33.inverse(matrix), matrix) + + for scale <- [2.0, -2.0, 1.0e-200, 1.0e200] do + assert_close(Mat33.make_reflect_3d({scale, 2 * scale, 2 * scale}), matrix) + end + end + + @tag :make_reflect_2d + @tag :make_reflect_3d + test "both reflection constructors reject zero normals" do + for normal <- [{0.0, 0.0}, {0, 0}, {0, 0.0}] do + assert_raise ArithmeticError, fn -> Mat33.make_reflect_2d(normal, 1.0) end + end + + for normal <- [{0.0, 0.0, 0.0}, {0, 0, 0}, {0, 0.0, 0}] do + assert_raise ArithmeticError, fn -> Mat33.make_reflect_3d(normal) end + end + end + + for {constructor, expected} <- [make_shear_x_2d: {-5, -4}, make_shear_y_2d: {3, 2}] do + @constructor constructor + @expected expected + @tag constructor + test "#{constructor} shears two spatial coordinates and preserves the homogeneous coordinate" do + for k <- [2.0, 2] do + matrix = apply(Mat33, @constructor, [k]) + assert Enum.all?(Tuple.to_list(matrix), &is_float/1) + assert_close(Mat33.transform_point(matrix, {3, -4}), @expected) + assert_close(Mat33.transform_vector(matrix, {3, -4}), @expected) + {x, y} = @expected + + for w <- [0.0, 1.0, 2.5] do + assert_close(Mat33.apply_left({3, -4, w}, matrix), {x, y, w}) + end + + assert Mat33.determinant(matrix) === 1.0 + assert_close(Mat33.multiply(matrix, apply(Mat33, @constructor, [-k])), Mat33.identity()) + end + + for k <- [0, 0.0], do: assert(apply(Mat33, @constructor, [k]) === Mat33.identity()) + end + end + + @tag :make_shear_x_2d + @tag :make_shear_y_2d + test "2D fractional shears compose with translation in row-vector order" do + x = Mat33.make_shear_x_2d(0.5) + y = Mat33.make_shear_y_2d(-0.5) + assert_close(Mat33.transform_point(x, {3, -4}), {1, -4}) + assert_close(Mat33.transform_point(y, {3, -4}), {3, -5.5}) + translation = Mat33.make_translate(10.0, -2.0) + assert_close(Mat33.transform_point(Mat33.multiply(x, translation), {3, -4}), {11, -6}) + assert_close(Mat33.transform_point(Mat33.multiply(translation, x), {3, -4}), {10, -6}) + assert_close(Mat33.transform_vector(Mat33.multiply(x, translation), {3, -4}), {1, -4}) + end + + for {constructor, expected} <- [ + make_shear_x_3d: {-15, -3, 4}, + make_shear_y_3d: {2, -5, 4}, + make_shear_z_3d: {2, -3, 16} + ] do + @constructor constructor + @expected expected + @tag constructor + test "#{constructor} uses both coefficients and has a known inverse" do + for {a, b} <- [{3.0, -2.0}, {3, -2}, {3, -2.0}, {3.0, -2}] do + matrix = apply(Mat33, @constructor, [a, b]) + assert Enum.all?(Tuple.to_list(matrix), &is_float/1) + assert_close(Mat33.apply_left({2, -3, 4}, matrix), @expected) + assert Mat33.determinant(matrix) === 1.0 + inverse = apply(Mat33, @constructor, [-a, -b]) + assert_close(Mat33.multiply(matrix, inverse), Mat33.identity()) + assert_close(Mat33.inverse(matrix), inverse) + end + + for k <- [0, 0.0], do: assert(apply(Mat33, @constructor, [k, k]) === Mat33.identity()) + end + end + + @tag :make_shear_x_3d + @tag :make_shear_y_3d + @tag :make_shear_z_3d + test "3D fractional coefficients retain their magnitude" do + assert_close(Mat33.apply_left({2, -3, 4}, Mat33.make_shear_x_3d(0.5, -0.25)), {-0.5, -3, 4}) + assert_close(Mat33.apply_left({2, -3, 4}, Mat33.make_shear_y_3d(0.5, -0.25)), {2, -3, 4}) + assert_close(Mat33.apply_left({2, -3, 4}, Mat33.make_shear_z_3d(0.5, -0.25)), {2, -3, 5.75}) + end + + @tag :make_shear_x_2d + @tag :make_shear_x_3d + test "2D and 3D shears explicitly distinguish homogeneous and spatial third coordinates" do + affine = Mat33.make_shear_x_2d(2.0) + linear = Mat33.make_shear_x_3d(2.0, 3.0) + assert_close(Mat33.apply_left({1, 2, 4}, affine), {5, 2, 4}) + assert_close(Mat33.apply_left({1, 2, 4}, linear), {17, 2, 4}) + 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