diff --git a/lib/graphmath/Mat44.ex b/lib/graphmath/Mat44.ex index 8c6a3f7..6157014 100644 --- a/lib/graphmath/Mat44.ex +++ b/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( diff --git a/test/mat44/mat44_reflection_shear_test.exs b/test/mat44/mat44_reflection_shear_test.exs new file mode 100644 index 0000000..7a1047c --- /dev/null +++ b/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