///|
/// Polygon triangulation, barycentric coordinates, and interpolation helpers.
pub(all) struct Triangle {
  first : Point
  second : Point
  third : Point
  area : Double
} derive(Debug, Eq, ToJson)

///|
pub(all) struct Barycentric {
  first : Double
  second : Double
  third : Double
  inside : Bool
} derive(Debug, Eq, ToJson)

///|
pub(all) struct MeshStatistics {
  vertices : Int
  triangles : Int
  area : Double
  minimum_area : Double
  maximum_area : Double
  boundary_length : Double
} derive(Debug, Eq, ToJson)

///|
pub fn triangle_area(first : Point, second : Point, third : Point) -> Double {
  ((second.x - first.x) * (third.y - first.y) -
  (third.x - first.x) * (second.y - first.y)).abs() /
  2.0
}

///|
pub fn Triangle::new(
  first : Point,
  second : Point,
  third : Point,
) -> Triangle raise @domain.ExtractionError {
  let area = triangle_area(first, second, third)
  if area <= 1.0e-15 {
    raise @domain.ExtractionError::InvalidArgument(
      message="triangle is degenerate",
    )
  }
  { first, second, third, area }
}

///|
pub fn Triangle::centroid(self : Triangle) -> Point {
  {
    x: (self.first.x + self.second.x + self.third.x) / 3.0,
    y: (self.first.y + self.second.y + self.third.y) / 3.0,
  }
}

///|
pub fn Triangle::contains(
  self : Triangle,
  point : Point,
  tolerance : Double,
) -> Bool {
  let barycentric = barycentric_coordinates(self, point, tolerance)
  barycentric.inside
}

///|
pub fn Triangle::perimeter(self : Triangle) -> Double {
  self.first.distance(self.second) +
  self.second.distance(self.third) +
  self.third.distance(self.first)
}

///|
pub fn barycentric_coordinates(
  triangle : Triangle,
  point : Point,
  tolerance : Double,
) -> Barycentric {
  let denominator = (triangle.second.y - triangle.third.y) *
    (triangle.first.x - triangle.third.x) +
    (triangle.third.x - triangle.second.x) *
    (triangle.first.y - triangle.third.y)
  if denominator.abs() <= 1.0e-15 {
    return { first: 0.0, second: 0.0, third: 0.0, inside: false }
  }
  let first = (
      (triangle.second.y - triangle.third.y) * (point.x - triangle.third.x) +
      (triangle.third.x - triangle.second.x) * (point.y - triangle.third.y)
    ) /
    denominator
  let second = (
      (triangle.third.y - triangle.first.y) * (point.x - triangle.third.x) +
      (triangle.first.x - triangle.third.x) * (point.y - triangle.third.y)
    ) /
    denominator
  let third = 1.0 - first - second
  {
    first,
    second,
    third,
    inside: first >= -tolerance && second >= -tolerance && third >= -tolerance,
  }
}

///|
pub fn interpolate_triangle(
  triangle : Triangle,
  first_value : Double,
  second_value : Double,
  third_value : Double,
  point : Point,
  tolerance : Double,
) -> Double raise @domain.ExtractionError {
  let coordinates = barycentric_coordinates(triangle, point, tolerance)
  if !coordinates.inside {
    raise @domain.ExtractionError::OutOfRange(
      name="point",
      value=0.0,
      lower=0.0,
      upper=1.0,
    )
  }
  coordinates.first * first_value +
  coordinates.second * second_value +
  coordinates.third * third_value
}

///|
pub fn fan_triangulate(
  polygon : ArrayView[Point],
) -> Array[Triangle] raise @domain.ExtractionError {
  if polygon.length() < 3 {
    raise @domain.ExtractionError::InsufficientData(
      required=3,
      actual=polygon.length(),
    )
  }
  let triangles = []
  let anchor = polygon[0]
  for index in 1..<(polygon.length() - 1) {
    triangles.push(Triangle::new(anchor, polygon[index], polygon[index + 1]))
  }
  triangles
}

///|
pub fn triangle_mesh_area(triangles : ArrayView[Triangle]) -> Double {
  triangles.fold(init=0.0, (sum, triangle) => sum + triangle.area)
}

///|
pub fn mesh_statistics(
  polygon : ArrayView[Point],
  triangles : ArrayView[Triangle],
) -> MeshStatistics {
  let mut minimum = if triangles.length() == 0 {
    0.0
  } else {
    triangles[0].area
  }
  let mut maximum = minimum
  for triangle in triangles {
    minimum = minimum.min(triangle.area)
    maximum = maximum.max(triangle.area)
  }
  let mut boundary = 0.0
  if polygon.length() > 1 {
    for index in 0.. Int {
  for index in 0.. Double raise @domain.ExtractionError {
  if triangles.length() != values.length() {
    raise @domain.ExtractionError::InvalidArgument(
      message="triangle values do not match mesh",
    )
  }
  let index = locate_triangle(triangles, point, tolerance)
  if index < 0 {
    raise @domain.ExtractionError::OutOfRange(
      name="point",
      value=0.0,
      lower=0.0,
      upper=1.0,
    )
  }
  let value = values[index]
  interpolate_triangle(
    triangles[index],
    value.0,
    value.1,
    value.2,
    point,
    tolerance,
  )
}

///|
pub fn triangle_aspect_ratio(triangle : Triangle) -> Double {
  let first = triangle.first.distance(triangle.second)
  let second = triangle.second.distance(triangle.third)
  let third = triangle.third.distance(triangle.first)
  let longest = first.max(second).max(third)
  let shortest = first.min(second).min(third)
  if shortest == 0.0 {
    0.0
  } else {
    longest / shortest
  }
}

///|
pub fn mesh_maximum_aspect_ratio(triangles : ArrayView[Triangle]) -> Double {
  let mut result = 0.0
  for triangle in triangles {
    result = result.max(triangle_aspect_ratio(triangle))
  }
  result
}

///|
pub fn mesh_quality_label(triangles : ArrayView[Triangle]) -> String {
  let aspect = mesh_maximum_aspect_ratio(triangles)
  if aspect <= 2.0 {
    "good"
  } else if aspect <= 5.0 {
    "usable"
  } else if aspect <= 20.0 {
    "thin"
  } else {
    "poor"
  }
}

///|
pub fn point_in_polygon(point : Point, polygon : ArrayView[Point]) -> Bool {
  if polygon.length() < 3 {
    return false
  }
  let mut inside = false
  let mut previous = polygon.length() - 1
  for index in 0.. point.y) != (previous_point.y > point.y)
    let intersection = if previous_point.y == current_point.y {
      false
    } else {
      point.x <
      (previous_point.x - current_point.x) *
      (point.y - current_point.y) /
      (previous_point.y - current_point.y) +
      current_point.x
    }
    if crosses && intersection {
      inside = !inside
    }
    previous = index
  }
  inside
}

///|
pub fn polygon_grid(
  bounds : BoundingBox,
  columns : Int,
  rows : Int,
) -> Array[Point] raise @domain.ExtractionError {
  if columns < 1 || rows < 1 {
    raise @domain.ExtractionError::InvalidArgument(
      message="grid dimensions must be positive",
    )
  }
  let points = []
  for row in 0.. String {
  let lines = ["index,x1,y1,x2,y2,x3,y3,area,aspect_ratio"]
  for index in 0.. Double raise @domain.ExtractionError {
  (triangle_mesh_area(triangles) - polygon_area(polygon)).abs()
}