///|
pub fn line_residuals(line : Line, points : Array[Point]) -> Array[Double] {
  points.map(fn(point) { line.distance_to(point) })
}

///|
pub fn line_fit_error(points : Array[Point]) -> Double {
  if points.length() < 2 {
    0.0
  } else {
    describe_samples(line_residuals(fit_line(points), points)).standard_deviation
  }
}

///|
pub fn robust_line_fit(
  points : Array[Point],
  iterations : Int,
  cutoff : Double,
) -> Line {
  if points.length() < 2 {
    Line::Line(a=0.0, b=1.0, c=0.0)
  } else {
    let active = points
    let mut line = fit_line(active)
    for _ in 0..= 2 {
        line = fit_line(selected)
      }
    }
    line
  }
}

///|
pub fn try_robust_line_fit(
  points : Array[Point],
  iterations : Int,
  cutoff : Double,
) -> Line? {
  if points.length() < 2 {
    None
  } else {
    Some(robust_line_fit(points, iterations, cutoff))
  }
}

///|
pub fn circle_residuals(
  circle : Circle,
  points : Array[Point],
) -> Array[Double] {
  points.map(fn(point) {
    (distance(circle.center, point) - circle.radius).abs()
  })
}

///|
pub fn circle_fit_error(points : Array[Point]) -> Double {
  if points.length() < 3 {
    0.0
  } else {
    describe_samples(
      circle_residuals(detect_circle_from_points(points), points),
    ).mean
  }
}

///|
pub fn robust_circle_fit(
  points : Array[Point],
  iterations : Int,
  cutoff : Double,
) -> Circle {
  let initial = detect_circle_from_points(points)
  let mut current = initial
  for _ in 0..= 3 {
      current = detect_circle_from_points(selected)
    }
  }
  current
}

///|
pub fn try_robust_circle_fit(
  points : Array[Point],
  iterations : Int,
  cutoff : Double,
) -> Circle? {
  if points.length() < 3 {
    None
  } else {
    Some(robust_circle_fit(points, iterations, cutoff))
  }
}

///|
pub fn ellipse_from_points(points : Array[Point]) -> Ellipse {
  if points.is_empty() {
    {
      center: Point::Point(x=0.0, y=0.0),
      major_radius: 0.0,
      minor_radius: 0.0,
      angle: 0.0,
      rms_error: 0.0,
    }
  } else {
    let center = polygon_centroid(points)
    let mut xx = 0.0
    let mut yy = 0.0
    let mut xy = 0.0
    for point in points {
      let dx = point.x - center.x
      let dy = point.y - center.y
      xx += dx * dx
      yy += dy * dy
      xy += dx * dy
    }
    let n = Double::from_int(points.length())
    let covariance_x = xx / n
    let covariance_y = yy / n
    let covariance_xy = xy / n
    let trace = covariance_x + covariance_y
    let determinant = covariance_x * covariance_y -
      covariance_xy * covariance_xy
    let root = (trace * trace / 4.0 - determinant).max(0.0).sqrt()
    let major_variance = (trace / 2.0 + root).max(0.0)
    let minor_variance = (trace / 2.0 - root).max(0.0)
    let angle = if covariance_xy.abs() <= 0.0000001 {
      0.0
    } else {
      @math.atan2(major_variance - covariance_x, covariance_xy)
    }
    let major = major_variance.sqrt()
    let minor = minor_variance.sqrt()
    let errors : Array[Double] = []
    for point in points {
      let rotated = rotate_point(point, center, -angle)
      let value = if major <= 0.0000001 || minor <= 0.0000001 {
        0.0
      } else {
        (rotated.x - center.x) * (rotated.x - center.x) / (major * major) +
        (rotated.y - center.y) * (rotated.y - center.y) / (minor * minor)
      }
      errors.push((value - 1.0).abs())
    }
    {
      center,
      major_radius: major,
      minor_radius: minor,
      angle,
      rms_error: describe_samples(errors).standard_deviation,
    }
  }
}

///|
pub fn ellipse_area(ellipse : Ellipse) -> Double {
  3.141592653589793 * ellipse.major_radius * ellipse.minor_radius
}

///|
pub fn ellipse_perimeter_approximation(ellipse : Ellipse) -> Double {
  let a = ellipse.major_radius.max(ellipse.minor_radius)
  let b = ellipse.major_radius.min(ellipse.minor_radius)
  let h = if a + b <= 0.0 {
    0.0
  } else {
    (a - b) * (a - b) / ((a + b) * (a + b))
  }
  3.141592653589793 *
  (a + b) *
  (1.0 + 3.0 * h / (10.0 + (4.0 - 3.0 * h).sqrt()))
}

///|
pub fn oriented_bounds(points : Array[Point], angle : Double) -> DoubleBounds {
  if points.is_empty() {
    { min: 0.0, max: 0.0 }
  } else {
    let center = polygon_centroid(points)
    let mut minimum = 0.0
    let mut maximum = 0.0
    let axis = Vector2::Vector2(x=@math.cos(angle), y=@math.sin(angle))
    for i in 0.. Double {
  if points.length() < 2 {
    0.0
  } else {
    let center = polygon_centroid(points)
    let mut xx = 0.0
    let mut yy = 0.0
    let mut xy = 0.0
    for point in points {
      let dx = point.x - center.x
      let dy = point.y - center.y
      xx += dx * dx
      yy += dy * dy
      xy += dx * dy
    }
    0.5 * @math.atan2(2.0 * xy, xx - yy)
  }
}

///|
pub fn geometric_median(points : Array[Point], iterations : Int) -> Point {
  if points.is_empty() {
    Point::Point(x=0.0, y=0.0)
  } else {
    let mut result = polygon_centroid(points)
    for _ in 0.. Double {
  let directed = fn(a : Array[Point], b : Array[Point]) {
    let mut result = 0.0
    for point in a {
      let mut best = 1.0e300
      for candidate in b {
        best = best.min(distance(point, candidate))
      }
      result = result.max(best)
    }
    result
  }
  directed(first, second).max(directed(second, first))
}

///|
pub fn chamfer_distance(first : Array[Point], second : Array[Point]) -> Double {
  if first.is_empty() || second.is_empty() {
    0.0
  } else {
    let mut total = 0.0
    for point in first {
      let mut best = 1.0e300
      for candidate in second {
        best = best.min(distance(point, candidate))
      }
      total += best
    }
    total / Double::from_int(first.length())
  }
}