///|
pub fn AffineTransform::AffineTransform(
  a~ : Double,
  b~ : Double,
  c~ : Double,
  d~ : Double,
  tx~ : Double,
  ty~ : Double,
) -> AffineTransform {
  { a, b, c, d, tx, ty }
}

///|
pub fn AffineTransform::identity() -> AffineTransform {
  { a: 1.0, b: 0.0, c: 0.0, d: 1.0, tx: 0.0, ty: 0.0 }
}

///|
pub fn AffineTransform::translation(
  dx~ : Double,
  dy~ : Double,
) -> AffineTransform {
  { a: 1.0, b: 0.0, c: 0.0, d: 1.0, tx: dx, ty: dy }
}

///|
pub fn AffineTransform::scaling(sx~ : Double, sy~ : Double) -> AffineTransform {
  { a: sx, b: 0.0, c: 0.0, d: sy, tx: 0.0, ty: 0.0 }
}

///|
pub fn AffineTransform::rotation(angle : Double) -> AffineTransform {
  let cosine = @math.cos(angle)
  let sine = @math.sin(angle)
  { a: cosine, b: -sine, c: sine, d: cosine, tx: 0.0, ty: 0.0 }
}

///|
pub fn AffineTransform::apply(
  transform : AffineTransform,
  point : Point,
) -> Point {
  Point::Point(
    x=transform.a * point.x + transform.b * point.y + transform.tx,
    y=transform.c * point.x + transform.d * point.y + transform.ty,
  )
}

///|
pub fn AffineTransform::determinant(transform : AffineTransform) -> Double {
  transform.a * transform.d - transform.b * transform.c
}

///|
pub fn AffineTransform::inverse(transform : AffineTransform) -> AffineTransform {
  let determinant = transform.determinant()
  if determinant.abs() <= 0.0000000001 {
    AffineTransform::identity()
  } else {
    let a = transform.d / determinant
    let b = -transform.b / determinant
    let c = -transform.c / determinant
    let d = transform.a / determinant
    {
      a,
      b,
      c,
      d,
      tx: -(a * transform.tx + b * transform.ty),
      ty: -(c * transform.tx + d * transform.ty),
    }
  }
}

///|
pub fn AffineTransform::compose(
  first : AffineTransform,
  second : AffineTransform,
) -> AffineTransform {
  {
    a: second.a * first.a + second.b * first.c,
    b: second.a * first.b + second.b * first.d,
    c: second.c * first.a + second.d * first.c,
    d: second.c * first.b + second.d * first.d,
    tx: second.a * first.tx + second.b * first.ty + second.tx,
    ty: second.c * first.tx + second.d * first.ty + second.ty,
  }
}

///|
fn determinant3(matrix : Array[Double]) -> Double {
  matrix[0] * (matrix[4] * matrix[8] - matrix[5] * matrix[7]) -
  matrix[1] * (matrix[3] * matrix[8] - matrix[5] * matrix[6]) +
  matrix[2] * (matrix[3] * matrix[7] - matrix[4] * matrix[6])
}

///|
fn solve3(matrix : Array[Double], rhs : Array[Double]) -> Array[Double] {
  let denominator = determinant3(matrix)
  if denominator.abs() <= 0.0000000001 {
    [0.0, 0.0, 0.0]
  } else {
    let x = determinant3([
        rhs[0],
        matrix[1],
        matrix[2],
        rhs[1],
        matrix[4],
        matrix[5],
        rhs[2],
        matrix[7],
        matrix[8],
      ]) /
      denominator
    let y = determinant3([
        matrix[0],
        rhs[0],
        matrix[2],
        matrix[3],
        rhs[1],
        matrix[5],
        matrix[6],
        rhs[2],
        matrix[8],
      ]) /
      denominator
    let z = determinant3([
        matrix[0],
        matrix[1],
        rhs[0],
        matrix[3],
        matrix[4],
        rhs[1],
        matrix[6],
        matrix[7],
        rhs[2],
      ]) /
      denominator
    [x, y, z]
  }
}

///|
pub fn AffineTransform::from_correspondences(
  source : Array[Point],
  target : Array[Point],
) -> AffineTransform {
  match AffineTransform::try_from_correspondences(source, target) {
    Some(transform) => transform
    None => AffineTransform::identity()
  }
}

///|
pub fn AffineTransform::try_from_correspondences(
  source : Array[Point],
  target : Array[Point],
) -> AffineTransform? {
  if source.length() < 3 ||
    target.length() < 3 ||
    source.length() != target.length() {
    None
  } else {
    let matrix = [
      source[0].x,
      source[0].y,
      1.0,
      source[1].x,
      source[1].y,
      1.0,
      source[2].x,
      source[2].y,
      1.0,
    ]
    if determinant3(matrix).abs() <= 0.0000000001 {
      None
    } else {
      let x = solve3(matrix, [target[0].x, target[1].x, target[2].x])
      let y = solve3(matrix, [target[0].y, target[1].y, target[2].y])
      Some({ a: x[0], b: x[1], c: y[0], d: y[1], tx: x[2], ty: y[2] })
    }
  }
}

///|
pub fn fit_affine(
  source : Array[Point],
  target : Array[Point],
) -> CalibrationFit {
  let transform = AffineTransform::from_correspondences(source, target)
  let residuals : Array[Double] = []
  if source.length() == target.length() {
    for i in 0.. Array[Point] {
  points.map(fn(point) { transform.apply(point) })
}

///|
pub fn rotate_point(point : Point, center : Point, angle : Double) -> Point {
  AffineTransform::translation(dx=center.x, dy=center.y)
  .compose(
    AffineTransform::rotation(angle).compose(
      AffineTransform::translation(dx=-center.x, dy=-center.y),
    ),
  )
  .apply(point)
}

///|
pub fn reflect_point(point : Point, line : Line) -> Point {
  let signed = line.a * point.x + line.b * point.y + line.c
  Point::Point(
    x=point.x - 2.0 * signed * line.a,
    y=point.y - 2.0 * signed * line.b,
  )
}

///|
pub fn project_point(point : Point, line : Line) -> Point {
  let signed = line.a * point.x + line.b * point.y + line.c
  Point::Point(x=point.x - signed * line.a, y=point.y - signed * line.b)
}

///|
pub fn line_intersection(first : Line, second : Line) -> Point {
  let denominator = first.a * second.b - second.a * first.b
  if denominator.abs() <= 0.0000000001 {
    Point::Point(x=0.0, y=0.0)
  } else {
    Point::Point(
      x=(first.b * second.c - second.b * first.c) / denominator,
      y=(second.a * first.c - first.a * second.c) / denominator,
    )
  }
}

///|
pub fn line_from_points(first : Point, second : Point) -> Line {
  let dx = second.x - first.x
  let dy = second.y - first.y
  if dx.abs() >= dy.abs() {
    Line::Line(a=-dy, b=dx, c=dy * first.x - dx * first.y)
  } else {
    Line::Line(a=-dy, b=dx, c=dy * first.x - dx * first.y)
  }
}

///|
pub fn angle_between_lines(first : Line, second : Line) -> Double {
  let dot = first.a * second.a + first.b * second.b
  @math.acos(dot.clamp(min=-1.0, max=1.0))
}

///|
pub fn normalize_angle(angle : Double) -> Double {
  let two_pi = 6.283185307179586
  let mut result = angle % two_pi
  if result <= -3.141592653589793 {
    result += two_pi
  }
  if result > 3.141592653589793 {
    result -= two_pi
  }
  result
}

///|
pub fn align_centroids(
  source : Array[Point],
  target : Array[Point],
) -> AffineTransform {
  if source.is_empty() || target.is_empty() {
    AffineTransform::identity()
  } else {
    let source_center = polygon_centroid(source)
    let target_center = polygon_centroid(target)
    AffineTransform::translation(
      dx=target_center.x - source_center.x,
      dy=target_center.y - source_center.y,
    )
  }
}

///|
pub fn transform_residuals(
  transform : AffineTransform,
  source : Array[Point],
  target : Array[Point],
) -> Array[Double] {
  if source.length() != target.length() {
    []
  } else {
    let result : Array[Double] = []
    for i in 0..