geof_shapes

defmodule GEOF.Shapes do
  @moduledoc """
    Shared utility functions for working with 3D geometry.
  """

  import :math

  @typedoc "A point in three dimensions."
  @type vector3 :: {number, number, number}
  @typedoc "A line in three dimensions represented by two points on the line."
  @type line3 :: {vector3, vector3}
  @typedoc "A triangle in three dimensions represented by its vertices."
  @type triangle3 :: {vector3, vector3, vector3}

  @doc "Given a `line` and a `triangle` in 3D space, determines whether the line intersects with the triangle."

  @spec line_intersects_triangle?(line3, triangle3) :: boolean

  def line_intersects_triangle?({l_a, l_b}, {p_0, p_1, p_2}) do
    l_ab = Vector.subtract(l_a, l_b)
    inverse_l_ab = Vector.reverse(l_ab)

    p_01 = Vector.subtract(p_1, p_0)
    p_02 = Vector.subtract(p_2, p_0)

    p_01x02 = Vector.cross(p_01, p_02)
    l_a_minus_p_0 = Vector.subtract(l_a, p_0)

    determinant = Vector.dot(inverse_l_ab, p_01x02)

    if determinant == 0 do
      false
    else
      t =
        Vector.dot(
          p_01x02,
          l_a_minus_p_0
        ) / determinant

      u =
        Vector.dot(
          Vector.cross(
            p_02,
            inverse_l_ab
          ),
          l_a_minus_p_0
        ) / determinant

      v =
        Vector.dot(
          Vector.cross(
            inverse_l_ab,
            p_01
          ),
          l_a_minus_p_0
        ) / determinant

      u >= 0 and v >= 0 and u + v <= 1 and t <= 0
    end
  end

  @ir3 1 / sqrt(3)

  @tetrahedron_in_unit_sphere {
    [
      # vectors
      {@ir3, @ir3, @ir3},
      {-@ir3, @ir3, -@ir3},
      {-@ir3, -@ir3, @ir3},
      {@ir3, -@ir3, -@ir3}
    ],
    [
      # faces
      [1, 3, 2],
      [0, 2, 3],
      [0, 3, 1],
      [0, 2, 1]
    ]
  }

  @tetrahedron_triangles (fn ->
                            {tetr_vects, tetr_faces} = @tetrahedron_in_unit_sphere

                            Enum.map(tetr_faces, fn face ->
                              {
                                Enum.at(tetr_vects, Enum.at(face, 0)),
                                Enum.at(tetr_vects, Enum.at(face, 1)),
                                Enum.at(tetr_vects, Enum.at(face, 2))
                              }
                            end)
                          end).()

  @doc "Given a `position` on a planetary sphere, determines which face of the unit tetrahedron the position lies within."

  @spec face_of_4_hedron(GEOF.Planet.Geometry.position()) :: pos_integer

  def face_of_4_hedron({:pos, lat, lon}) do
    x = cos(lat) * cos(lon)
    y = cos(lat) * sin(lon)
    z = sin(lat)

    line = {
      {0, 0, 0},
      Vector.multiply({x, y, z}, 1.1)
    }

    Enum.find_index(@tetrahedron_triangles, fn triangle ->
      line_intersects_triangle?(line, triangle)
    end)
  end

  @doc "Given a `position` on a planetary sphere, determines which face of the unit octahedron the position lies within."

  @spec face_of_8_hedron(GEOF.Planet.Geometry.position()) :: pos_integer

  def face_of_8_hedron({:pos, lat, lon}) do
    x_offset = if cos(lat) * cos(lon) > 0, do: 1, else: 0
    z_offset = if cos(lat) * sin(lon) > 0, do: 2, else: 0
    y_offset = if sin(lat) > 0, do: 4, else: 0

    x_offset + z_offset + y_offset
  end
end