path: root/src/linear_algebra.rs

#![forbid(unsafe_code)]
use std::mem::{ size_of, MaybeUninit };
use std::ops::{ Add, Sub, Mul, Div, Neg };
use vulkanalia::vk::{ self, HasBuilder };


pub trait Real: Add<Output=Self> + Sub<Output=Self>
                + Mul<Output=Self> + Div<Output=Self>
                + Neg<Output=Self>
                + Sized
                + From<f32>
{
  fn sin_cos(self) -> (Self, Self);
  fn tan(self) -> Self;
  fn recip(self) -> Self;
}


impl Real for f32 {
  fn sin_cos(self) -> (Self, Self) { f32::sin_cos(self) }
  fn tan(self) -> Self { f32::tan(self) }
  fn recip(self) -> Self { f32::recip(self) }
}


pub static VERTICES: [Vertex<f32>; 4] = [
  Vertex::new(Vec2::new(-0.5, -0.5), Vec3::new(1.0, 0.0, 0.0)),
  Vertex::new(Vec2::new( 0.5, -0.5), Vec3::new(0.0, 1.0, 0.0)),
  Vertex::new(Vec2::new( 0.5,  0.5), Vec3::new(0.0, 0.0, 1.0)),
  Vertex::new(Vec2::new(-0.5,  0.5), Vec3::new(1.0, 1.0, 1.0)),
];

pub const INDICES: &[u16] = &[0, 1, 2, 2, 3, 0];


#[repr(C)]
#[derive(Clone, Debug)]
pub struct Vec2<T: Copy>([T; 2]);

//   The GLSL spec version 4.60.8, in section 4.4, defers to the OpenGL spec
// for the meaning of the "std140" qualifier. The OpenGL spec version 4.6
// core profile, in section 7.6.2.2, states that the base alignment of a
// 3-vector is the same as that for a 4-vector. It further states that, when
// a structure contains another structure, the offset of the item
// immediately following the sub-structure will be coerced to a multiple of
// the alignment of the sub-structure. Nothing actually states that a vector
// counts as a sub-structure for this purpose, and other nearby wording
// makes a distinction between vectors and structures, but as-implemented,
// it appears to.
//
//   So, we need padding.
//
//   That padding needs to be of the same size as T, and we want to do this
// with the bare minimum restrictions on what T can be. Furthermore, our
// constructors are const, and that's important to us. This may seem
// impossible, but MaybeUninit comes to the rescue! It does require T to be
// Copy; this is the reason that Copy is threaded everywhere in all these
// types.
#[repr(C)]
#[derive(Clone, Debug)]
pub struct Vec3<T: Copy>([T; 3], MaybeUninit<T>);

#[repr(C)]
#[derive(Clone, Debug)]
#[allow(unused)]
pub struct Vec4<T: Copy>([T; 4]);

// All our matrices are column-major, which is what GLSL expects.
#[repr(C)]
#[derive(Clone, Debug)]
#[allow(unused)]
pub struct Mat2<T: Copy>([Vec2<T>; 2]);

//   Using Vec3 here means we get its padding behavior, which is what we want,
// per that same section of the OpenGL spec.
#[repr(C)]
#[derive(Clone, Debug)]
#[allow(unused)]
pub struct Mat3<T: Copy>([Vec3<T>; 3]);

#[repr(C)]
#[derive(Clone, Debug)]
pub struct Mat4<T: Copy>([Vec4<T>; 4]);


impl<T: Copy> Vec2<T> {
  //   We use letters from the beginning of the alphabet for all of these, so
  // as not to take a stance on any of the TODO finish this sentence
  pub const fn new(a: T, b: T) -> Self {
    Vec2([a, b])
  }
}

impl<T: Copy> Vec3<T> {
  pub const fn new(a: T, b: T, c: T) -> Self {
    Vec3([a, b, c], MaybeUninit::zeroed())
  }
}

impl<T: Copy> Vec4<T> {
  pub const fn new(a: T, b: T, c: T, d: T) -> Self {
    Vec4([a, b, c, d])
  }

  // This is kept carefully the same as quaternion_conjugate() in shader.vert.
  #[allow(unused)]
  pub fn quaternion_conjugate(&self) -> Self
    where T: Real
  {
    Vec4::new(self.0[0], -self.0[1], -self.0[2], -self.0[3])
  }

  // This is kept carefully the same as quaternion_product() in shader.vert.
  pub fn quaternion_product(&self, other: &Self) -> Self
      where T: Real
  {
    let a = &self.0;
    let b = &other.0;

    Vec4::new((a[0] * b[0]) - (a[1] * b[1]) - (a[2] * b[2]) - (a[3] * b[3]),
              (a[0] * b[1]) + (a[1] * b[0]) + (a[2] * b[3]) - (a[3] * b[2]),
              (a[0] * b[2]) - (a[1] * b[3]) + (a[2] * b[0]) + (a[3] * b[1]),
              (a[0] * b[3]) + (a[1] * b[2]) - (a[2] * b[1]) + (a[3] * b[0]))
  }

  //   The input angle is in radians. The resulting rotation is clockwise when
  // looking in the positive direction of the axis, which is sometimes
  // intuitive and sometimes not.
  pub fn rotation_quaternion(axis: &Vec3<T>, angle: T) -> Self
      where T: Real
  {
    let (sine, cosine) = (angle / 2.0.into()).sin_cos();

    Self::new(cosine, sine * axis.0[0], sine * axis.0[1], sine * axis.0[2])
  }

  // A convenience method.
  pub fn rotate(&self, axis: &Vec3<T>, angle: T) -> Self
      where T: Real
  {
    //   Rotation quaternions are composed by multiplying them. The one on the
    // left side of the product is applied second, and the one on the right
    // side is applied first. So, we do this in the "reverse" order because
    // we're notionally composing the new rotation after the one we already
    // have.
    Self::rotation_quaternion(axis, angle).quaternion_product(self)
  }
}

impl<T: Copy> Mat2<T> {
  #[allow(unused)]
  pub const fn new(a: Vec2<T>, b: Vec2<T>) -> Self {
    Mat2([a, b])
  }
}

//   You may be tempted to reimplement Mat3<T> as a synonym for Vec3<Vec3<T>>.
// This would be incorrect, because the padding behavior only applies at the
// inner layer.
impl<T: Copy> Mat3<T> {
  #[allow(unused)]
  pub const fn new(a: Vec3<T>, b: Vec3<T>, c: Vec3<T>) -> Self {
    Mat3([a, b, c])
  }
}

impl<T: Copy> Mat4<T> {
  pub const fn new(a: Vec4<T>, b: Vec4<T>, c: Vec4<T>, d: Vec4<T>) -> Self {
    Mat4([a, b, c, d])
  }

  //   This produces a 4x4 matrix suitable for use in a shader to transform
  // coordinates. The matrix expects input in a Cartesian coordinate system
  // where positive x is right, positive y is up, and positive z is deeper
  // into the screen. Such a system is left-handed. The matrix's output
  // coordinate system ranges from -1 to 1 in the x and y axes, and from 0 to
  // 1 in the z axis, with the x and y axes being Euclidean and the z axis
  // being homogenous - see below.
  //
  //   To apply the matrix, right-multiply with it. That is, compute a dot
  // product where the right input is the matrix, and the left input is a
  // position vector. The position must be a 4-vector with coordinates x, y,
  // and z followed by a fourth value that is always 1.0. This fixed final
  // coordinate of 1.0 allows the matrix to perform translation in addition to
  // scaling and rotation.
  //
  //   The result of the multiplication is a 4-vector with coordinates x, y,
  // z, and w, with the expectation that the other three coordinates will be
  // divided by w as part of the final viewport transformation performed
  // internally to the graphics API. This is necessary in order to produce the
  // homogenous behavior of the z-coordinate; it also acts as a convenient way
  // to scale the x and y coordinates by depth, implementing what artists call
  // foreshortening.
  //
  //   The field of view is in radians, and is used as the angle measured
  // vertically from top to bottom of the screen. The computed y scale will be
  // directly based on this field of view, while the x scale will be based on
  // the field of view and the aspect ratio. Thus, in setups where the screen
  // is wider than it is tall, the horizontal field of view will be larger
  // than the one given.
  //
  //   The final z value, after the matrix multiplication and w-division have
  // both been applied, is homogenous rather than Cartesian: coordinate values
  // designate points further apart in space the farther they get from the
  // origin. This is important because floating-point numbers have limited
  // precision, and we want to spend most of that precision on the most
  // visually noticeable things, which are the closest ones. Anyway, graphics
  // APIs don't actually give us a choice about that. :)
  //
  //   The near and far planes are used as the clipping boundaries for z
  // values. An input position with z equal to the near value given here will
  // have a final (multiplied and w-divided) z-value of 0.0; an input with z
  // equal to the far value will have a final z of 1.0.
  pub fn perspective(vertical_field_of_view: T, aspect_ratio: T,
                     near: T, far: T)
      -> Self
      where T: Real
  {
    //   We compute some coefficients at the top here, which are best
    // understood in the context of the full matrix, below.
    let y_scale = (vertical_field_of_view / 2.0.into()).tan().recip();
    let x_scale = y_scale / aspect_ratio;
    let z_scale = far / (far - near);
    let z_offset = -near * z_scale;

    //   This matrix will be treated as being in column-major order. That is,
    // each Vec4 becomes a column. In the shader, we will then right-multiply
    // with it, meaning we'll compute the dot product of the point with the
    // matrix in that order. The result of this is that the items in the
    // first vector are summed to give the final "x", the second vector gives
    // the final "y", and so on. It expects the input point to have a value of
    // 1.0 for w (the final item), which allows using the final item of each
    // Vec4 to represent a constant term.
    Mat4::new(Vec4::new(x_scale, 0.0.into(), 0.0.into(), 0.0.into()),
              Vec4::new(0.0.into(), y_scale, 0.0.into(), 0.0.into()),
              Vec4::new(0.0.into(), 0.0.into(), z_scale, z_offset),
              Vec4::new(0.0.into(), 0.0.into(), 1.0.into(), 0.0.into()))

    //   These coefficients require significant explanation.
    // Counterintuitively, z_offset is the one that ends up contributing the
    // part that varies with the input z, though it is not multiplied by z
    // during the matrix multiplication. Notice the 1.0 in the bottom row,
    // which sets the output w equal to the input z. When the w-division later
    // happens, the z factor that was multiplied with z_scale cancels out,
    // while z_offset winds up being inversely proportionate to the input z.
    //
    //   Writing the full formula for the final z after matrix multiplication
    // and w-division in terms of the input z, we get:
    //
    //   z_final = [far / (far-near)] - [ near*far / z*(far-near)]
    //
    //   Notice how both the constant part (from z_scale) and the part that
    // varies with z (from z_offset) have (far - near) in the denominator.
    // This is the part responsible for scaling the coordinates to the range
    // 0.0 to 1.0. To see why z = near maps to z_final = 0.0 and z = far maps
    // to z_final = 1.0, try substituting and simplifying.
    //
    //   Also consider the subexpression near*far / z. This is the way the
    // non-linear behavior needed for a homogenous coordinate is produced.
    // When z is at its minimum, this subexpression reduces to just "far";
    // when z is at its maximum, the subexpression reduces to just "near".
    // In between, the values are clustered densely around the "near" end. In
    // the context of the full expression, this becomes the zero-to-one range.
  }
}


#[repr(C)]
#[derive(Clone, Debug)]
pub struct Vertex<T: Copy> {
  pub position: Vec2<T>,
  pub color: Vec3<T>,
}

#[repr(C)]
#[derive(Clone, Debug)]
pub struct Transformation<T: Copy> {
  pub rotation: Vec4<T>,
  pub translation: Vec3<T>,
}

#[repr(C)]
#[derive(Clone, Debug)]
pub struct UniformBlock<T: Copy> {
  pub scale: Vec3<T>,
  pub model: Transformation<T>,
  pub view: Transformation<T>,
  pub projection: Mat4<T>,
}

impl<T: Copy> Vertex<T> {
  const fn new(position: Vec2<T>, color: Vec3<T>) -> Self {
    Self { position, color }
  }

  pub fn binding_description() -> vk::VertexInputBindingDescription {
    vk::VertexInputBindingDescription::builder()
        .binding(0)
        .stride(size_of::<Vertex<f32>>() as u32)
        .input_rate(vk::VertexInputRate::VERTEX)
        .build()
  }

  pub fn attribute_descriptions()
      -> [vk::VertexInputAttributeDescription; 2]
  {
    let position = vk::VertexInputAttributeDescription::builder()
                       .binding(0)
                       .location(0)
                       .format(vk::Format::R32G32_SFLOAT)
                       .offset(0)
                       .build();

    let color = vk::VertexInputAttributeDescription::builder()
                       .binding(0)
                       .location(1)
                       .format(vk::Format::R32G32B32_SFLOAT)
                       .offset(size_of::<Vec2<f32>>() as u32)
                       .build();

    [position, color]
  }
}