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 atan2(self, other: Self) -> Self;
  fn recip(self) -> Self;
  fn powi(self, power: i32) -> Self;
  fn sqrt(self) -> Self;
}


impl Real for f32 {
  fn sin_cos(self) -> (Self, Self) { f32::sin_cos(self) }
  fn tan(self) -> Self { f32::tan(self) }
  fn atan2(self, other: Self) -> Self { f32::atan2(self, other) }
  fn recip(self) -> Self { f32::recip(self) }
  fn powi(self, power: i32) -> Self { f32::powi(self, power) }
  fn sqrt(self) -> Self { f32::sqrt(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 preference any specific use of these types to represent points,
  // colors, rotations etc over others. There's a lot of variations in which
  // coordinate is which, and picking "proper" names for them would inevitably
  // wind up disagreeing with how other code does it.
  pub const fn new(a: T, b: T) -> Self {
    Vec2([a, b])
  }

  #[allow(unused)]
  pub fn negate(&self) -> Self
      where T: Real
  {
    Vec2::new(-self.0[0], -self.0[1])
  }

  #[allow(unused)]
  pub fn add(&self, other: &Self) -> Self
      where T: Real
  {
    Vec2::new(self.0[0] + other.0[0],
              self.0[1] + other.0[1])
  }

  #[allow(unused)]
  pub fn subtract(&self, other: &Self) -> Self
      where T: Real
  {
    Vec2::new(self.0[0] - other.0[0],
              self.0[1] - other.0[1])
  }

  #[allow(unused)]
  pub fn scalar_product(&self, scalar: T) -> Self
      where T: Real
  {
    Vec2::new(self.0[0] * scalar,
              self.0[1] * scalar)
  }

  #[allow(unused)]
  pub fn scalar_quotient(&self, scalar: T) -> Self
      where T: Real
  {
    Vec2::new(self.0[0] / scalar,
              self.0[1] / scalar)
  }

  #[allow(unused)]
  pub fn magnitude(&self) -> T
      where T: Real
  {
    (self.0[0].powi(2) + self.0[1].powi(2)).sqrt()
  }

  #[allow(unused)]
  pub fn normalize(&self) -> Self
      where T: Real
  {
    self.scalar_quotient(self.magnitude())
  }

  #[allow(unused)]
  pub fn dot_product(&self, other: &Self) -> T
      where T: Real
  {
    let a = self.0;
    let b = other.0;

    (a[0]*b[0]) + (a[1]*b[1])
  }

  #[allow(unused)]
  pub fn distance(&self, other: &Self) -> T
      where T: Real
  {
    self.subtract(other).magnitude()
  }
}

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

  pub fn negate(&self) -> Self
      where T: Real
  {
    Vec3::new(-self.0[0], -self.0[1], -self.0[2])
  }

  #[allow(unused)]
  pub fn add(&self, other: &Self) -> Self
      where T: Real
  {
    Vec3::new(self.0[0] + other.0[0],
              self.0[1] + other.0[1],
              self.0[2] + other.0[2])
  }

  pub fn subtract(&self, other: &Self) -> Self
      where T: Real
  {
    Vec3::new(self.0[0] - other.0[0],
              self.0[1] - other.0[1],
              self.0[2] - other.0[2])
  }

  pub fn scalar_product(&self, scalar: T) -> Self
      where T: Real
  {
    Vec3::new(self.0[0] * scalar,
              self.0[1] * scalar,
              self.0[2] * scalar)
  }

  pub fn scalar_quotient(&self, scalar: T) -> Self
      where T: Real
  {
    Vec3::new(self.0[0] / scalar,
              self.0[1] / scalar,
              self.0[2] / scalar)
  }

  pub fn magnitude(&self) -> T
      where T: Real
  {
    (self.0[0].powi(2) + self.0[1].powi(2) + self.0[2].powi(2)).sqrt()
  }

  pub fn normalize(&self) -> Self
      where T: Real
  {
    self.scalar_quotient(self.magnitude())
  }

  pub fn dot_product(&self, other: &Self) -> T
      where T: Real
  {
    let a = self.0;
    let b = other.0;

    (a[0]*b[0]) + (a[1]*b[1]) + (a[2]*b[2])
  }

  pub fn cross_product(&self, other: &Self) -> Self
      where T: Real
  {
    let a = self.0;
    let b = other.0;

    Vec3::new(a[1]*b[2] - a[2]*b[1],
              a[2]*b[0] - a[0]*b[2],
              a[0]*b[1] - a[1]*b[0])
  }

  pub fn distance(&self, other: &Self) -> T
      where T: Real
  {
    self.subtract(other).magnitude()
  }

  pub fn rotate_position_by_quaternion(&self, rotation: &Vec4<T>) -> Self
      where T: Real
  {
    let position = Vec4::new(0.0.into(), self.0[0], self.0[1], self.0[2]);
    let result = rotation.quaternion_product(&position)
                         .quaternion_product(&rotation.quaternion_conjugate());
    Vec3::new(result.0[1], result.0[2], result.0[3])
  }
}

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

  #[allow(unused)]
  pub fn negate(&self) -> Self
      where T: Real
  {
    Vec4::new(-self.0[0], -self.0[1], -self.0[2], -self.0[3])
  }

  #[allow(unused)]
  pub fn add(&self, other: &Self) -> Self
      where T: Real
  {
    Vec4::new(self.0[0] + other.0[0],
              self.0[1] + other.0[1],
              self.0[2] + other.0[2],
              self.0[3] + other.0[3])
  }

  #[allow(unused)]
  pub fn subtract(&self, other: &Self) -> Self
      where T: Real
  {
    Vec4::new(self.0[0] - other.0[0],
              self.0[1] - other.0[1],
              self.0[2] - other.0[2],
              self.0[3] - other.0[3])
  }

  #[allow(unused)]
  pub fn scalar_product(&self, scalar: T) -> Self
      where T: Real
  {
    Vec4::new(self.0[0] * scalar,
              self.0[1] * scalar,
              self.0[2] * scalar,
              self.0[3] * scalar)
  }

  #[allow(unused)]
  pub fn scalar_quotient(&self, scalar: T) -> Self
      where T: Real
  {
    Vec4::new(self.0[0] / scalar,
              self.0[1] / scalar,
              self.0[2] / scalar,
              self.0[3] / scalar)
  }

  #[allow(unused)]
  pub fn magnitude(&self) -> T
      where T: Real
  {
    (self.0[0].powi(2) + self.0[1].powi(2) + self.0[2].powi(2)
     + self.0[3].powi(2)).sqrt()
  }

  #[allow(unused)]
  pub fn normalize(&self) -> Self
      where T: Real
  {
    self.scalar_quotient(self.magnitude())
  }

  #[allow(unused)]
  pub fn dot_product(&self, other: &Self) -> T
      where T: Real
  {
    let a = self.0;
    let b = other.0;

    (a[0]*b[0]) + (a[1]*b[1]) + (a[2]*b[2]) + (a[3]*b[3])
  }

  #[allow(unused)]
  pub fn distance(&self, other: &Self) -> T
      where T: Real
  {
    self.subtract(other).magnitude()
  }

  // 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 down, and positive z is deeper
  // into the screen. Such a system is right-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 Transformation<T: Copy> {
  pub rotation: Vec4<T>,
  pub translation: Vec3<T>,
}

impl<T: Copy> Transformation<T> {
  //   This computes a Transformation suitable for use as a "view" transform,
  // describing a camera situated at the position given by "vantage" and
  // pointed at the position given by "target". Because this would be
  // under-constrained on its own, we also require an "up" vector which says
  // which axis, and which direction along it, should be oriented towards the
  // top of the screen.
  //
  //   We make the assumption that the transform we return will be applied to
  // a starting situation in which the camera is looking along the positive z
  // axis, with the top of the screen towards the negative y axis.
  //
  //   TODO we should really add one final rotation which fixes the camera to
  // actually have the up vector as its up-axis. Right now we use it for
  // pitching the camera but not for anything else, and that feels incorrect.
  pub fn look_at(vantage: &Vec3<T>, target: &Vec3<T>, up: &Vec3<T>) -> Self
      where T: Real + std::fmt::Debug
  {
    //   We're given an up vector as absolute orientation. We start by using
    // dot product to separate our target and vantage positions into scalar
    // values along that axis, and vectors representing everything BUT that
    // axis. You can think of these "planar" vectors as representing the
    // shadows of the two positions on a notional ground plane.
    let target_up = target.dot_product(&up);
    let vantage_up = vantage.dot_product(&up);
    let target_planar = target.subtract(&up.scalar_product(target_up));
    let vantage_planar = vantage.subtract(&up.scalar_product(vantage_up));

    //   Now we can use the "planar" positions to compute a vector measured
    // along that plane. We'll use that vector as an offset in the angle
    // computations below, and we'll use its normalized version right here as
    // an axis for a cross-product that will give us a "right" vector that
    // is used in one of our rotations. More on its meaning below.
    let planar_offset = target_planar.subtract(&vantage_planar);
    let forward = planar_offset.normalize();
    let right = forward.cross_product(&up).normalize();

    //   We do two rotations. Here, we figure out the first one, which is the
    // "yaw" to orient ourselves left-right. Our caller has given us a camera
    // pointing in the positive-z direction, and we need to rotate from that
    // global reference frame into one pointing directly at the object. We
    // compute the legs of a right triangle using the x and z ades in global
    // coordinates, and use the arctangent to get the angle at its center.
    //
    //   Notice that this doesn't take the "up" vector into account. It's not
    // clear that doing so is a sensible thing for this step.
    //
    //   By the way, atan2() has less precision loss than doing the division
    // ourselves and using plain atan(). If you didn't know that, now you
    // know. It's a very thoughtful function in that way.
    let right_distance = planar_offset.0[0];
    let forward_distance = planar_offset.0[2];
    let horizontal_angle = T::atan2(right_distance, forward_distance);

    //   Now the second rotation, the "pitch" to orient up-down. Again we
    // compute two legs of a triangle and use the arctangent to get an angle
    // from it.
    let vertical_distance = vantage_up - target_up;
    let planar_distance = vantage_planar.distance(&target_planar);
    let vertical_angle = T::atan2(vertical_distance, planar_distance);

    //   Here, we compute the two quaternions for our rotations and compose
    // them into a single one. Doing this here on the CPU on a per-frame basis
    // will mean the GPU has less per-vertex work to do.
    //
    //   The axis for the "yaw" rotation is the up vector. The axis for the
    // "pitch" rotation is the "right" vector we computed. We promised to
    // explain that "right" vector: It's pointing directly out to the right of
    // the camera, in the position it NOW has, which is already yawed. Thus,
    // rotating about it now will only affect the up-down direction.
    let rotation = Vec4::rotation_quaternion(&up, horizontal_angle)
                       .rotate(&right, vertical_angle);

    //   The semantics of a Transformation are that the rotation is always
    // applied first, followed by the translation. That's natural and
    // intuitive for model transformations, since the transformation is
    // applied to the coordinates of the thing being drawn, and the model's
    // local origin will presumably be at some logical point that corresponds
    // to people's expectations.
    //
    //   It is not intuitive for view transformations, such as the one we're
    // returning, since the transformation is STILL applied to the thing being
    // drawn, though notionally what we care about is instead the camera. We
    // resolve this by applying the rotation we compute to the translation,
    // so that the combined Transformation will move the model in such a way
    // that it's as if the camera has moved.
    //
    //   Notice that, since this translation describes the camera's position
    // in world space, it's the only calculation we do that uses the vantage
    // by itself, rather than the offset between the vantage and the target.
    let translation = vantage.negate()
                             .rotate_position_by_quaternion(&rotation);

    Transformation { rotation, translation }
  }
}



#[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 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]
  }
}