///|
/// Numerical diagnostics for sampled paths used by plotter and print exports.
pub(all) struct GeometryProfile {
  point_count : Int
  segment_count : Int
  length : Double
  area : Double
  bounds : Rect
  centroid : Vec2
  min_segment : Double
  max_segment : Double
  mean_segment : Double
  curvature_mean : Double
  curvature_max : Double
  closed : Bool
} derive(Eq, Debug, ToJson)

///|
pub(all) struct GeometryCompare {
  hausdorff : Double
  mean_error : Double
  length_error : Double
  area_error : Double
  centroid_error : Double
  score : Double
} derive(Eq, Debug, ToJson)

///|
fn ga_points(path : Path, subdivisions : Int) -> Array[Vec2] {
  path.flatten(subdivisions=subdivisions.max(1))
}

///|
pub fn geometry_centroid(points : Array[Vec2]) -> Vec2 {
  if points.length() == 0 {
    Vec2::zero()
  } else {
    let mut x = 0.0
    let mut y = 0.0
    for point in points {
      x += point.x
      y += point.y
    }
    Vec2::new(x / points.length().to_double(), y / points.length().to_double())
  }
}

///|
pub fn geometry_segment_lengths(
  points : Array[Vec2],
  closed : Bool,
) -> Array[Double] {
  let result = []
  if points.length() > 1 {
    for index in 0..<(points.length() - 1) {
      result.push(points[index].distance(points[index + 1]))
    }
    if closed {
      result.push(points[points.length() - 1].distance(points[0]))
    }
  }
  result
}

///|
pub fn geometry_angles(points : Array[Vec2], closed : Bool) -> Array[Double] {
  let result = []
  let count = points.length()
  if count >= 3 {
    let start = if closed { 0 } else { 1 }
    let end = if closed { count } else { count - 1 }
    for index in start.. Array[Double] {
  let angles = geometry_angles(points, closed)
  let lengths = geometry_segment_lengths(points, closed)
  let result = []
  for index, angle in angles {
    let denominator = if lengths.length() == 0 {
      1.0
    } else {
      lengths[index.min(lengths.length() - 1)].max(0.0001)
    }
    result.push(angle / denominator)
  }
  result
}

///|
pub fn geometry_profile(
  path : Path,
  subdivisions : Int,
  closed : Bool,
) -> GeometryProfile {
  let points = ga_points(path, subdivisions)
  let lengths = geometry_segment_lengths(points, closed)
  let curvatures = geometry_curvatures(points, closed)
  let mut min_segment = if lengths.length() == 0 { 0.0 } else { lengths[0] }
  let mut max_segment = min_segment
  let mut total = 0.0
  for length in lengths {
    min_segment = min_segment.min(length)
    max_segment = max_segment.max(length)
    total += length
  }
  let mut curvature_total = 0.0
  let mut curvature_max = 0.0
  for curvature in curvatures {
    curvature_total += curvature
    curvature_max = curvature_max.max(curvature)
  }
  let count = curvatures.length()
  {
    point_count: points.length(),
    segment_count: lengths.length(),
    length: total,
    area: polygon_area(points).abs(),
    bounds: bounds_of_points(points),
    centroid: geometry_centroid(points),
    min_segment,
    max_segment,
    mean_segment: if lengths.length() == 0 {
      0.0
    } else {
      total / lengths.length().to_double()
    },
    curvature_mean: if count == 0 {
      0.0
    } else {
      curvature_total / count.to_double()
    },
    curvature_max,
    closed,
  }
}

///|
pub fn geometry_resample(
  points : Array[Vec2],
  count : Int,
  closed : Bool,
) -> Array[Vec2] {
  let result = []
  if points.length() == 0 || count <= 0 {
    return result
  }
  let lengths = geometry_segment_lengths(points, closed)
  let total = lengths.fold(init=0.0, fn(sum, value) { sum + value })
  let denominator = if closed { count } else { (count - 1).max(1) }
  for index in 0.. Array[Vec2] {
  if points.length() <= 2 {
    points.copy()
  } else {
    let result = [points[0]]
    let threshold = tolerance.max(0.0)
    for index in 1..<(points.length() - 1) {
      let previous = result[result.length() - 1]
      if previous.distance(points[index]) >= threshold {
        result.push(points[index])
      }
    }
    result.push(points[points.length() - 1])
    result
  }
}

///|
pub fn geometry_smooth(
  points : Array[Vec2],
  radius : Int,
  amount : Double,
  closed : Bool,
) -> Array[Vec2] {
  let result = points.copy()
  let strength = amount.max(0.0).min(1.0)
  if points.length() > 2 {
    for index in 0.. Double {
  fn directed(source : Array[Vec2], target : Array[Vec2]) -> Double {
    let mut maximum = 0.0
    for point in source {
      let mut nearest = 1000000.0
      for candidate in target {
        nearest = nearest.min(point.distance(candidate))
      }
      maximum = maximum.max(nearest)
    }
    maximum
  }
  if left.length() == 0 && right.length() == 0 {
    0.0
  } else if left.length() == 0 || right.length() == 0 {
    1000000.0
  } else {
    directed(left, right).max(directed(right, left))
  }
}

///|
pub fn geometry_compare(
  left : Path,
  right : Path,
  subdivisions : Int,
  closed : Bool,
) -> GeometryCompare {
  let left_points = ga_points(left, subdivisions)
  let right_points = ga_points(right, subdivisions)
  let left_profile = geometry_profile(left, subdivisions, closed)
  let right_profile = geometry_profile(right, subdivisions, closed)
  let left_samples = geometry_resample(left_points, 32, closed)
  let right_samples = geometry_resample(right_points, 32, closed)
  let hausdorff = geometry_hausdorff(left_samples, right_samples)
  let mut mean_error = 0.0
  for index, point in left_samples {
    mean_error += point.distance(
      right_samples[index.min(right_samples.length() - 1)],
    )
  }
  mean_error /= left_samples.length().max(1).to_double()
  let length_error = (left_profile.length - right_profile.length).abs()
  let area_error = (left_profile.area - right_profile.area).abs()
  let centroid_error = left_profile.centroid.distance(right_profile.centroid)
  let scale = left_profile.bounds.width.max(left_profile.bounds.height).max(1.0)
  let score = (1.0 -
    (
      hausdorff +
      mean_error +
      length_error / scale +
      area_error / scale / scale +
      centroid_error / scale
    ) /
    5.0)
    .max(0.0)
    .min(1.0)
  { hausdorff, mean_error, length_error, area_error, centroid_error, score }
}

///|
pub fn geometry_profile_csv(profile : GeometryProfile) -> String {
  "point_count,segment_count,length,area,min_segment,max_segment,mean_segment,curvature_mean,curvature_max,closed\n" +
  profile.point_count.to_string() +
  "," +
  profile.segment_count.to_string() +
  "," +
  profile.length.to_string() +
  "," +
  profile.area.to_string() +
  "," +
  profile.min_segment.to_string() +
  "," +
  profile.max_segment.to_string() +
  "," +
  profile.mean_segment.to_string() +
  "," +
  profile.curvature_mean.to_string() +
  "," +
  profile.curvature_max.to_string() +
  "," +
  profile.closed.to_string() +
  "\n"
}

///|
pub fn geometry_compare_report(compare : GeometryCompare) -> String {
  "hausdorff=" +
  compare.hausdorff.to_string() +
  " mean_error=" +
  compare.mean_error.to_string() +
  " length_error=" +
  compare.length_error.to_string() +
  " area_error=" +
  compare.area_error.to_string() +
  " centroid_error=" +
  compare.centroid_error.to_string() +
  " score=" +
  compare.score.to_string()
}

///|
pub fn geometry_polyline_to_path(points : Array[Vec2], closed : Bool) -> Path {
  if points.length() == 0 {
    Path::new()
  } else {
    let path = Path::new().move_to(points[0].x, points[0].y)
    let mut result = path
    for point in points[1:] {
      result = result.line_to(point.x, point.y)
    }
    if closed {
      result.close()
    } else {
      result
    }
  }
}

///|
pub fn geometry_offset(
  points : Array[Vec2],
  distance : Double,
  closed : Bool,
) -> Array[Vec2] {
  let result = []
  if points.length() == 0 {
    result
  } else {
    for index, point in points {
      let previous = if index == 0 {
        if closed {
          points[points.length() - 1]
        } else {
          point
        }
      } else {
        points[index - 1]
      }
      let next = if index + 1 >= points.length() {
        if closed {
          points[0]
        } else {
          point
        }
      } else {
        points[index + 1]
      }
      let tangent = next.sub(previous).normalized()
      let normal = Vec2::new(-tangent.y, tangent.x)
      result.push(point.add(normal.scale(distance)))
    }
    result
  }
}

///|
pub fn geometry_bounds_area(points : Array[Vec2]) -> Double {
  let bounds = bounds_of_points(points)
  bounds.width * bounds.height
}

///|
pub fn geometry_is_degenerate(
  points : Array[Vec2],
  area_tolerance : Double,
  length_tolerance : Double,
) -> Bool {
  polygon_area(points).abs() <= area_tolerance.max(0.0) &&
  polygon_perimeter(points) <= length_tolerance.max(0.0)
}

///|
pub fn geometry_normalize(points : Array[Vec2], bounds : Rect) -> Array[Vec2] {
  let result = []
  for point in points {
    result.push(
      Vec2::new(
        (point.x - bounds.x) / bounds.width.max(0.0001),
        (point.y - bounds.y) / bounds.height.max(0.0001),
      ),
    )
  }
  result
}

///|
pub fn geometry_denormalize(points : Array[Vec2], bounds : Rect) -> Array[Vec2] {
  let result = []
  for point in points {
    result.push(
      Vec2::new(
        bounds.x + point.x * bounds.width,
        bounds.y + point.y * bounds.height,
      ),
    )
  }
  result
}