///|
fn contour_crosses(a : Double, b : Double, level : Double) -> Bool {
  is_finite(a) &&
  is_finite(b) &&
  ((a < level && b >= level) || (a >= level && b < level))
}

///|
fn contour_intersection(
  a : GeoPoint,
  b : GeoPoint,
  av : Double,
  bv : Double,
  level : Double,
) -> GeoPoint {
  let denominator = bv - av
  let t = if denominator == 0.0 { 0.5 } else { (level - av) / denominator }
  interpolate_point(a, b, t.clamp(min=0.0, max=1.0))
}

///|
fn contour_edges(
  points : Array[GeoPoint],
  values : Array[Double],
  level : Double,
) -> Array[GeoPoint] {
  let intersections : Array[GeoPoint] = []
  for edge in [(0, 1), (1, 2), (2, 3), (3, 0)] {
    let (a, b) = edge
    if contour_crosses(values[a], values[b], level) {
      intersections.push(
        contour_intersection(points[a], points[b], values[a], values[b], level),
      )
    }
  }
  intersections
}

///|
/// Extract marching-squares segments for one contour level.
pub fn Grid::contour(self : Grid, level : Double) -> Array[ContourSegment] {
  let output : Array[ContourSegment] = []
  if self.rows < 2 || self.cols < 2 {
    output
  } else {
    for row in 0..<(self.rows - 1) {
      for col in 0..<(self.cols - 1) {
        match
          (
            self.cell(row, col),
            self.cell(row, col + 1),
            self.cell(row + 1, col + 1),
            self.cell(row + 1, col),
          ) {
          (Some(nw), Some(ne), Some(se), Some(sw)) => {
            let points = [nw.point, ne.point, se.point, sw.point]
            let values = [nw.value, ne.value, se.value, sw.value]
            let intersections = contour_edges(points, values, level)
            if intersections.length() == 2 {
              output.push({
                start: intersections[0],
                finish: intersections[1],
                level,
              })
            } else if intersections.length() == 4 {
              output.push({
                start: intersections[0],
                finish: intersections[1],
                level,
              })
              output.push({
                start: intersections[2],
                finish: intersections[3],
                level,
              })
            }
          }
          _ => ()
        }
      }
    }
    output
  }
}

///|
pub fn Grid::contours(
  self : Grid,
  levels : ArrayView[Double],
) -> Array[ContourSegment] {
  let output : Array[ContourSegment] = []
  for level in levels {
    for segment in self.contour(level) {
      output.push(segment)
    }
  }
  output
}

///|
pub fn contours_to_geojson(segments : ArrayView[ContourSegment]) -> String {
  let features : Array[String] = []
  for segment in segments {
    let coordinates = "[[\{segment.start.lon},\{segment.start.lat}],[\{segment.finish.lon},\{segment.finish.lat}]]"
    features.push(
      "{\"type\":\"Feature\",\"geometry\":{\"type\":\"LineString\",\"coordinates\":\{coordinates}},\"properties\":{\"level\":\{segment.level}}}",
    )
  }
  let encoded = features.join(",")
  "{\"type\":\"FeatureCollection\",\"features\":[\{encoded}]}"
}

///|
pub fn contour_level_range(grid : Grid, count : Int) -> Array[Double] {
  let output : Array[Double] = []
  if count > 0 {
    let range = grid.value_range()
    for i in 0.. Double {
  haversine_distance_m(segment.start, segment.finish)
}

///|
pub fn contours_total_length(segments : ArrayView[ContourSegment]) -> Double {
  let mut total = 0.0
  for segment in segments {
    total += contour_segment_length(segment)
  }
  total
}