///|
pub fn simplify_route(
  points : ArrayView[Point],
  tolerance_m? : Double = 10.0,
) -> Result[Array[Point], RouteError] {
  if tolerance_m.is_nan() || tolerance_m.is_inf() || tolerance_m < 0.0 {
    return Err(InvalidTolerance(tolerance_m~))
  }
  match require_points(points) {
    Ok(_) => ()
    Err(err) => return Err(err)
  }
  if points.length() <= 2 || tolerance_m == 0.0 {
    return Ok(points.to_owned())
  }
  let keep = Array::make(points.length(), false)
  keep[0] = true
  keep[points.length() - 1] = true
  mark_simplify(points, 0, points.length() - 1, tolerance_m, keep)
  let out = Array::new()
  for i, point in points {
    if keep[i] {
      out.push(point)
    }
  }
  Ok(out)
}

///|
fn mark_simplify(
  points : ArrayView[Point],
  start : Int,
  end : Int,
  tolerance_m : Double,
  keep : Array[Bool],
) -> Unit {
  if end <= start + 1 {
    return
  }
  let mut best_index = -1
  let mut best_distance = -1.0
  for i in (start + 1).. best_distance {
      best_distance = distance
      best_index = i
    }
  }
  if best_index >= 0 && best_distance > tolerance_m {
    keep[best_index] = true
    mark_simplify(points, start, best_index, tolerance_m, keep)
    mark_simplify(points, best_index, end, tolerance_m, keep)
  }
}

///|
fn point_segment_distance_m(point : Point, a : Point, b : Point) -> Double {
  let ref_lat = (a.lat + b.lat) / 2.0
  let (px, py) = project_for_distance(point, ref_lat)
  let (ax, ay) = project_for_distance(a, ref_lat)
  let (bx, by) = project_for_distance(b, ref_lat)
  let dx = bx - ax
  let dy = by - ay
  if dx == 0.0 && dy == 0.0 {
    return ((px - ax) * (px - ax) + (py - ay) * (py - ay)).sqrt()
  }
  let t = (((px - ax) * dx + (py - ay) * dy) / (dx * dx + dy * dy)).clamp(
    min=0.0,
    max=1.0,
  )
  let nearest_x = ax + t * dx
  let nearest_y = ay + t * dy
  ((px - nearest_x) * (px - nearest_x) + (py - nearest_y) * (py - nearest_y)).sqrt()
}

///|
fn project_for_distance(point : Point, ref_lat : Double) -> (Double, Double) {
  let x = deg_to_rad(point.lon) *
    EARTH_RADIUS_M *
    @math.cos(deg_to_rad(ref_lat))
  let y = deg_to_rad(point.lat) * EARTH_RADIUS_M
  (x, y)
}