Add Mat44 perspective and orthographic projections

lib/graphmath/Mat44.ex

@@ -1,13 +1,18 @@
 defmodule Graphmath.Mat44 do
   @moduledoc """
-  4x4 matrices and 3D affine transformations, using tuples of floats.
+  4x4 matrices and 3D transformations, 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.
 
+  `perspective/4` and `ortho/6` project right-handed camera coordinates with -Z
+  forward into the OpenGL normalized device depth range [-1, +1]. Perspective
+  projection requires keeping the output of `apply_left/2` and dividing its
+  first three coordinates by `w`; `transform_point/2` does not perform that step.
+
   Tuples store matrix rows in order. `apply(a, v)` computes the column-vector
   product **A****v**. Graphics constructors use row vectors: `transform_point/2`
   and `transform_vector/2` multiply `{x,y,z,1}` and `{x,y,z,0}` by **A** from the
@@ -255,6 +260,87 @@ 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 right-handed perspective projection for row vectors.
+
+  `fov_y` is the vertical field of view in radians, `aspect` is width / height,
+  and `near` and `far` are positive distances along the camera's -Z direction.
+  Requires `0 < fov_y < pi`, `aspect > 0`, and `0 < near < far`.
+  Invalid parameters raise `ArithmeticError`.
+
+  Maps the near and far planes to normalized device depths -1 and +1,
+  respectively, following the OpenGL depth convention. Apply with `apply_left/2`
+  to a homogeneous point `{x, y, z, 1.0}`, then divide the first three output
+  coordinates by the fourth (`w = -z`). No clipping or division is performed
+  by this constructor. `transform_point/2` discards `w` and is unsuitable for
+  completing a perspective projection.
+
+  ## Examples
+
+      iex> projection = Graphmath.Mat44.perspective(:math.pi() / 2, 1.0, 1.0, 3.0)
+      iex> {x, y, z, w} = Graphmath.Mat44.apply_left({0.0, 0.0, -1.0, 1.0}, projection)
+      iex> {x / w, y / w, z / w}
+      {0.0, 0.0, -1.0}
+  """
+  @spec perspective(float, float, float, float) :: mat44
+  def perspective(fov_y, aspect, near, far)
+      when is_float(fov_y) and is_float(aspect) and is_float(near) and is_float(far) do
+    if fov_y <= 0.0 or fov_y >= :math.pi() or aspect <= 0.0 or near <= 0.0 or far <= near do
+      raise ArithmeticError,
+            "perspective requires 0 < fov_y < pi, aspect > 0, and 0 < near < far"
+    end
+
+    cotangent = 1.0 / :math.tan(fov_y / 2.0)
+    depth_ratio = near / (far - near)
+
+    # Equivalent to -(far + near)/(far - near) and -2*far*near/(far - near),
+    # without multiplying the near and far distances together.
+    {cotangent / aspect, 0.0, 0.0, 0.0, 0.0, cotangent, 0.0, 0.0, 0.0, 0.0,
+     -(1.0 + 2.0 * depth_ratio), -1.0, 0.0, 0.0, -2.0 * near * (1.0 + depth_ratio), 0.0}
+  end
+
+  def perspective(fov_y, aspect, near, far) do
+    perspective(1.0 * fov_y, 1.0 * aspect, 1.0 * near, 1.0 * far)
+  end
+
+  @doc """
+  Creates a right-handed orthographic projection for row vectors.
+
+  Maps X bounds `x_min` and `x_max` to -1 and +1, and Y bounds `y_min` and
+  `y_max` to -1 and +1. The camera looks down -Z: eye-space `z = -near`
+  maps to depth -1, and `z = -far` maps to depth +1 (the OpenGL convention).
+
+  Each pair of bounds must have distinct endpoints; equal endpoints raise
+  `ArithmeticError`. Reversed bounds flip the corresponding axis. Unlike
+  `perspective/4`, `near` and `far` may be zero or negative, allowing volumes
+  that extend behind the eye.
+
+  Preserves the homogeneous coordinate, so `transform_point/2` can apply this
+  projection directly. With `apply_left/2`, a point's output `w` remains 1.
+
+  ## Examples
+
+      iex> projection = Graphmath.Mat44.ortho(-2.0, 6.0, -4.0, 4.0, 1.0, 3.0)
+      iex> Graphmath.Mat44.transform_point(projection, {2.0, 0.0, -2.0})
+      {0.0, 0.0, 0.0}
+  """
+  @spec ortho(float, float, float, float, float, float) :: mat44
+  def ortho(x_min, x_max, y_min, y_max, near, far)
+      when is_float(x_min) and is_float(x_max) and is_float(y_min) and is_float(y_max) and
+             is_float(near) and is_float(far) do
+    if x_min == x_max or y_min == y_max or near == far do
+      raise ArithmeticError, "orthographic bounds must have nonzero extent on every axis"
+    end
+
+    {2.0 / (x_max - x_min), 0.0, 0.0, 0.0, 0.0, 2.0 / (y_max - y_min), 0.0, 0.0, 0.0, 0.0,
+     -2.0 / (far - near), 0.0, -(x_max + x_min) / (x_max - x_min),
+     -(y_max + y_min) / (y_max - y_min), -(far + near) / (far - near), 1.0}
+  end
+
+  def ortho(x_min, x_max, y_min, y_max, near, far) do
+    ortho(1.0 * x_min, 1.0 * x_max, 1.0 * y_min, 1.0 * y_max, 1.0 * near, 1.0 * far)
+  end
+
   @doc """
   Creates a 3D affine reflection across `nx*x + ny*y + nz*z = offset`.
 

test/mat44/mat44_projection_test.exs

@@ -0,0 +1,189 @@
+defmodule GraphmathTest.Mat44.Projection do
+  use ExUnit.Case, async: true
+  alias Graphmath.{Mat44, Vec4}
+
+  @moduletag :mat44
+  doctest Graphmath.Mat44, only: [perspective: 4, ortho: 6]
+
+  for kind <- [:integer, :float, :mixed] do
+    @kind kind
+
+    @tag :perspective
+    test "perspective has known entries with #{@kind} parameters" do
+      matrix = apply(Mat44, :perspective, numeric([1, 2, 1, 5], @kind))
+
+      assert_close(
+        matrix,
+        {0.915243860856226, 0, 0, 0, 0, 1.830487721712452, 0, 0, 0, 0, -1.5, -1, 0, 0, -2.5, 0}
+      )
+
+      assert Enum.all?(Tuple.to_list(matrix), &is_float/1)
+      assert_close(ndc(matrix, {0, 0, -1}), {0, 0, -1})
+      assert_close(ndc(matrix, {0, 0, -5}), {0, 0, 1})
+    end
+
+    @tag :ortho
+    test "an asymmetric orthographic volume has known entries with #{@kind} parameters" do
+      matrix = apply(Mat44, :ortho, numeric([-2, 6, -4, 2, 1, 9], @kind))
+      assert_close(matrix, {0.25, 0, 0, 0, 0, 1 / 3, 0, 0, 0, 0, -0.25, 0, -0.5, 1 / 3, -1.25, 1})
+      assert Enum.all?(Tuple.to_list(matrix), &is_float/1)
+      assert_close(Mat44.transform_point(matrix, {-2, -4, -1}), {-1, -1, -1})
+      assert_close(Mat44.transform_point(matrix, {6, 2, -9}), {1, 1, 1})
+      assert_close(Mat44.transform_point(matrix, {2, -1, -5}), {0, 0, 0})
+    end
+  end
+
+  @tag :perspective
+  test "all near and far frustum corners map to the expected clip boundaries" do
+    matrix = Mat44.perspective(:math.pi() / 2, 2.0, 2.0, 10.0)
+
+    for {distance, expected_z} <- [{2.0, -1}, {10.0, 1}], x <- [-1, 1], y <- [-1, 1] do
+      assert_close(ndc(matrix, {x * 2 * distance, y * distance, -distance}), {x, y, expected_z})
+    end
+  end
+
+  @tag :ortho
+  test "all eight orthographic corners map to the expected boundaries" do
+    matrix = Mat44.ortho(-2.0, 6.0, -4.0, 2.0, 1.0, 9.0)
+
+    for {x, nx} <- [{-2, -1}, {6, 1}],
+        {y, ny} <- [{-4, -1}, {2, 1}],
+        {z, nz} <- [{-1, -1}, {-9, 1}] do
+      assert_close(Mat44.apply_left({x, y, z, 1}, matrix), {nx, ny, nz, 1})
+    end
+  end
+
+  @tag :perspective
+  @tag :ortho
+  test "fractional aspect ratios and clipping distances retain their values" do
+    perspective = Mat44.perspective(:math.pi() / 2, 0.5, 0.5, 4.5)
+    assert_close(ndc(perspective, {0.25, 0.5, -0.5}), {1, 1, -1})
+    assert_close(ndc(perspective, {-2.25, -4.5, -4.5}), {-1, -1, 1})
+    ortho = Mat44.ortho(-1.5, 2.5, -0.25, 1.75, 0.5, 4.5)
+    assert_close(Mat44.transform_point(ortho, {0.5, 0.75, -2.5}), {0, 0, 0})
+  end
+
+  @tag :perspective
+  test "perspective foreshortening and nonlinear depth use division by w" do
+    matrix = Mat44.perspective(:math.pi() / 2, 2.0, 1.0, 9.0)
+    assert_close(ndc(matrix, {2, 1, -2}), {0.5, 0.5, 0.125})
+    assert_close(ndc(matrix, {2, 1, -4}), {0.25, 0.25, 0.6875})
+    assert_close(Mat44.apply_left({2, 1, -4, 1}, matrix), {1, 1, 2.75, 4})
+  end
+
+  @tag :perspective
+  test "perspective produces w from negative camera-space z without clipping" do
+    matrix = Mat44.perspective(:math.pi() / 2, 1.0, 1.0, 3.0)
+    assert_close(Mat44.apply_left({0, 0, -2, 1}, matrix), {0, 0, 1, 2})
+    assert_close(Mat44.apply_left({0, 0, 2, 1}, matrix), {0, 0, -7, -2})
+    assert_close(Mat44.apply_left({0, 0, 0, 1}, matrix), {0, 0, -3, 0})
+  end
+
+  @tag :perspective
+  test "depth coefficients avoid overflow or underflow from multiplying near and far" do
+    for scale <- [1.0e-200, 1.0e200] do
+      matrix = Mat44.perspective(:math.pi() / 2, 1.0, scale, 3 * scale)
+      assert_close(ndc(matrix, {scale, scale, -2 * scale}), {0.5, 0.5, 0.5})
+      assert_close(ndc(matrix, {0, 0, -scale}), {0, 0, -1})
+      assert_close(ndc(matrix, {0, 0, -3 * scale}), {0, 0, 1})
+    end
+  end
+
+  @tag :ortho
+  test "reversing a pair of bounds flips that projected axis" do
+    for {bounds, expected} <- [
+          {[6, -2, -4, 2, 1, 9], {1, -1, -1}},
+          {[-2, 6, 2, -4, 1, 9], {-1, 1, -1}},
+          {[-2, 6, -4, 2, 9, 1], {-1, -1, 1}},
+          {[6, -2, 2, -4, 9, 1], {1, 1, 1}}
+        ] do
+      matrix = apply(Mat44, :ortho, bounds)
+      assert_close(Mat44.transform_point(matrix, {-2, -4, -1}), expected)
+    end
+  end
+
+  @tag :ortho
+  test "orthographic near and far bounds may be zero or negative" do
+    for {near, far} <- [{0.0, 10.0}, {-2.0, 2.0}, {-5.0, -1.0}] do
+      matrix = Mat44.ortho(-1.0, 1.0, -1.0, 1.0, near, far)
+      assert_close(Mat44.transform_point(matrix, {0, 0, -near}), {0, 0, -1})
+      assert_close(Mat44.transform_point(matrix, {0, 0, -far}), {0, 0, 1})
+    end
+  end
+
+  @tag :perspective
+  test "invalid perspective volumes raise ArithmeticError through numeric clauses" do
+    for kind <- [:integer, :float, :mixed],
+        args <- [
+          [0, 1, 1, 10],
+          [-1, 1, 1, 10],
+          [4, 1, 1, 10],
+          [1, 0, 1, 10],
+          [1, -1, 1, 10],
+          [1, 1, 0, 10],
+          [1, 1, -1, 10],
+          [1, 1, 2, 2],
+          [1, 1, 2, 1],
+          [1, 1, 1, 0]
+        ] do
+      assert_raise ArithmeticError, fn -> apply(Mat44, :perspective, numeric(args, kind)) end
+    end
+
+    assert_raise ArithmeticError, fn -> Mat44.perspective(:math.pi(), 1.0, 1.0, 10.0) end
+  end
+
+  @tag :ortho
+  test "zero extent on any orthographic axis raises through numeric clauses" do
+    for kind <- [:integer, :float, :mixed],
+        args <- [[1, 1, -2, 2, 1, 9], [-2, 2, 1, 1, 1, 9], [-2, 2, -2, 2, 1, 1]] do
+      assert_raise ArithmeticError, fn -> apply(Mat44, :ortho, numeric(args, kind)) end
+    end
+  end
+
+  @tag :perspective
+  @tag :ortho
+  @tag :multiply
+  test "row-vector composition applies the view before projection" do
+    # A camera at {0, 0, 5}, looking down -Z, has this view translation.
+    view = Mat44.make_translate(0.0, 0.0, -5.0)
+    perspective = Mat44.perspective(:math.pi() / 2, 2.0, 1.0, 9.0)
+    ortho = Mat44.ortho(-2.0, 6.0, -4.0, 2.0, 1.0, 9.0)
+    assert_close(ndc(Mat44.multiply(view, perspective), {2, 1, 1}), {0.25, 0.25, 0.6875})
+    assert_close(Mat44.transform_point(Mat44.multiply(view, ortho), {2, 1, 1}), {0, 2 / 3, -0.25})
+  end
+
+  @tag :perspective
+  @tag :ortho
+  @tag :inverse
+  test "inverse projection recovers the full homogeneous point" do
+    for matrix <- [
+          Mat44.perspective(1.0, 2.0, 1.0, 5.0),
+          Mat44.ortho(-2.0, 6.0, -4.0, 2.0, 1.0, 9.0)
+        ] do
+      point = {0.5, -0.75, -3.0, 1.0}
+      clip = Mat44.apply_left(point, matrix)
+      assert_close(Mat44.apply_left(clip, Mat44.inverse(matrix)), point)
+    end
+  end
+
+  defp ndc(matrix, point) do
+    {x, y, z, w} = Mat44.apply_left(Vec4.from_point3(point), matrix)
+    {x / w, y / w, z / w}
+  end
+
+  defp numeric(values, kind) do
+    values
+    |> Enum.with_index()
+    |> Enum.map(fn {value, i} ->
+      if kind == :float or (kind == :mixed and rem(i, 2) == 0), do: value * 1.0, else: value
+    end)
+  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