///|
/// An undirected mesh edge with usage information.
pub(all) struct MeshEdge {
  left : Int
  right : Int
  length : Double
  face_count : Int
  boundary : Bool
} derive(Eq, Debug, ToJson)

///|
/// Quality statistics for a triangle mesh.
pub(all) struct MeshQuality {
  vertex_count : Int
  triangle_count : Int
  edge_count : Int
  total_area : Double
  minimum_area : Double
  maximum_area : Double
  mean_area : Double
  minimum_angle : Double
  degenerate_triangles : Int
  boundary_edges : Int
} derive(Eq, Debug, ToJson)

///|
/// Parameters for deterministic Laplacian mesh smoothing.
pub(all) struct MeshRelaxOptions {
  iterations : Int
  amount : Double
  preserve_boundary : Bool
} derive(Eq, Debug, ToJson)

///|
pub fn MeshRelaxOptions::new() -> MeshRelaxOptions {
  { iterations: 1, amount: 0.5, preserve_boundary: true }
}

///|
pub fn MeshRelaxOptions::with_iterations(
  self : MeshRelaxOptions,
  iterations : Int,
) -> MeshRelaxOptions {
  { ..self, iterations: iterations.max(0) }
}

///|
pub fn MeshRelaxOptions::with_amount(
  self : MeshRelaxOptions,
  amount : Double,
) -> MeshRelaxOptions {
  { ..self, amount: amount.max(0.0).min(1.0) }
}

///|
pub fn MeshRelaxOptions::allow_boundary(
  self : MeshRelaxOptions,
) -> MeshRelaxOptions {
  { ..self, preserve_boundary: false }
}

///|
/// A barycentric coordinate triple for interpolation inside one triangle.
pub(all) struct Barycentric {
  u : Double
  v : Double
  w : Double
  inside : Bool
} derive(Eq, Debug, ToJson)

///|
/// Returns a vertex when the index is valid.
pub fn TriangleMesh::vertex(self : TriangleMesh, index : Int) -> Vec2? {
  if index < 0 || index >= self.vertices.length() {
    None
  } else {
    Some(self.vertices[index])
  }
}

///|
pub fn TriangleMesh::vertex_count(self : TriangleMesh) -> Int {
  self.vertices.length()
}

///|
pub fn TriangleMesh::triangle_count(self : TriangleMesh) -> Int {
  self.indices.length()
}

///|
pub fn TriangleMesh::triangle_vertices(
  self : TriangleMesh,
  index : Int,
) -> Array[Vec2] {
  if index < 0 || index >= self.indices.length() {
    []
  } else {
    let triangle = self.indices[index]
    [
      self.vertices[triangle.0],
      self.vertices[triangle.1],
      self.vertices[triangle.2],
    ]
  }
}

///|
pub fn TriangleMesh::triangle_centroid(
  self : TriangleMesh,
  index : Int,
) -> Vec2 {
  let points = self.triangle_vertices(index)
  if points.length() != 3 {
    Vec2::zero()
  } else {
    points[0].add(points[1]).add(points[2]).scale(1.0 / 3.0)
  }
}

///|
pub fn TriangleMesh::triangle_perimeter(
  self : TriangleMesh,
  index : Int,
) -> Double {
  let points = self.triangle_vertices(index)
  if points.length() != 3 {
    0.0
  } else {
    polygon_perimeter(points)
  }
}

///|
pub fn TriangleMesh::triangle_incenter(
  self : TriangleMesh,
  index : Int,
) -> Vec2 {
  let points = self.triangle_vertices(index)
  if points.length() != 3 {
    Vec2::zero()
  } else {
    let a = points[1].distance(points[2])
    let b = points[0].distance(points[2])
    let c = points[0].distance(points[1])
    let total = a + b + c
    if total == 0.0 {
      points[0]
    } else {
      points[0]
      .scale(a / total)
      .add(points[1].scale(b / total))
      .add(points[2].scale(c / total))
    }
  }
}

///|
pub fn TriangleMesh::triangle_circumcenter(
  self : TriangleMesh,
  index : Int,
) -> Vec2? {
  let points = self.triangle_vertices(index)
  if points.length() != 3 {
    None
  } else {
    let a = points[0]
    let b = points[1]
    let c = points[2]
    let denominator = 2.0 *
      (a.x * (b.y - c.y) + b.x * (c.y - a.y) + c.x * (a.y - b.y))
    if denominator.abs() < 0.0000001 {
      None
    } else {
      let aa = a.dot(a)
      let bb = b.dot(b)
      let cc = c.dot(c)
      Some(
        Vec2::new(
          (aa * (b.y - c.y) + bb * (c.y - a.y) + cc * (a.y - b.y)) / denominator,
          (aa * (c.x - b.x) + bb * (a.x - c.x) + cc * (b.x - a.x)) / denominator,
        ),
      )
    }
  }
}

///|
fn canonical_edge(left : Int, right : Int) -> (Int, Int) {
  if left < right {
    (left, right)
  } else {
    (right, left)
  }
}

///|
fn edge_face_count(mesh : TriangleMesh, edge : (Int, Int)) -> Int {
  let mut count = 0
  for triangle in mesh.indices {
    let sides = [
      (triangle.0, triangle.1),
      (triangle.1, triangle.2),
      (triangle.2, triangle.0),
    ]
    for side in sides {
      if canonical_edge(side.0, side.1) == edge {
        count += 1
      }
    }
  }
  count
}

///|
/// Enumerates unique edges and marks boundary edges.
pub fn TriangleMesh::edges(self : TriangleMesh) -> Array[MeshEdge] {
  let unique : Array[(Int, Int)] = []
  for triangle in self.indices {
    for
      side in [
        (triangle.0, triangle.1),
        (triangle.1, triangle.2),
        (triangle.2, triangle.0),
      ] {
      let edge = canonical_edge(side.0, side.1)
      if !unique.contains(edge) {
        unique.push(edge)
      }
    }
  }
  let result = []
  for edge in unique {
    let count = edge_face_count(self, edge)
    let length = if edge.0 < self.vertices.length() &&
      edge.1 < self.vertices.length() {
      self.vertices[edge.0].distance(self.vertices[edge.1])
    } else {
      0.0
    }
    result.push({
      left: edge.0,
      right: edge.1,
      length,
      face_count: count,
      boundary: count == 1,
    })
  }
  result
}

///|
pub fn TriangleMesh::boundary_edges(self : TriangleMesh) -> Array[MeshEdge] {
  let result = []
  for edge in self.edges() {
    if edge.boundary {
      result.push(edge)
    }
  }
  result
}

///|
pub fn TriangleMesh::vertex_neighbors(
  self : TriangleMesh,
  vertex : Int,
) -> Array[Int] {
  let result = []
  for edge in self.edges() {
    if edge.left == vertex {
      result.push(edge.right)
    }
    if edge.right == vertex {
      result.push(edge.left)
    }
  }
  result
}

///|
pub fn TriangleMesh::vertex_faces(
  self : TriangleMesh,
  vertex : Int,
) -> Array[Int] {
  let result = []
  for index, triangle in self.indices {
    if triangle.0 == vertex || triangle.1 == vertex || triangle.2 == vertex {
      result.push(index)
    }
  }
  result
}

///|
pub fn TriangleMesh::quality(self : TriangleMesh) -> MeshQuality {
  let mut total = 0.0
  let mut minimum = if self.indices.length() == 0 {
    0.0
  } else {
    self.triangle_area(0)
  }
  let mut maximum = minimum
  let mut minimum_angle = @math.PI
  let mut degenerate = 0
  for index in 0.. Barycentric {
  let v0 = b.sub(a)
  let v1 = c.sub(a)
  let v2 = point.sub(a)
  let denominator = v0.cross(v1)
  if denominator.abs() < 0.0000001 {
    { u: 0.0, v: 0.0, w: 0.0, inside: false }
  } else {
    let v = v2.cross(v1) / denominator
    let w = v0.cross(v2) / denominator
    let u = 1.0 - v - w
    { u, v, w, inside: u >= -0.000001 && v >= -0.000001 && w >= -0.000001 }
  }
}

///|
pub fn TriangleMesh::locate_triangle(
  self : TriangleMesh,
  point : Vec2,
) -> (Int, Barycentric)? {
  for index in 0.. Double? {
  if index < 0 || index >= mesh.indices.length() {
    None
  } else {
    let triangle = mesh.indices[index]
    if triangle.0 >= values.length() ||
      triangle.1 >= values.length() ||
      triangle.2 >= values.length() {
      None
    } else {
      let barycentric = triangle_barycentric(
        point,
        mesh.vertices[triangle.0],
        mesh.vertices[triangle.1],
        mesh.vertices[triangle.2],
      )
      if !barycentric.inside {
        None
      } else {
        Some(
          values[triangle.0] * barycentric.u +
          values[triangle.1] * barycentric.v +
          values[triangle.2] * barycentric.w,
        )
      }
    }
  }
}

///|
fn mesh_boundary_vertex(mesh : TriangleMesh, index : Int) -> Bool {
  for edge in mesh.boundary_edges() {
    if edge.left == index || edge.right == index {
      return true
    }
  }
  false
}

///|
/// Performs Laplacian smoothing while optionally preserving the boundary.
pub fn TriangleMesh::relax(
  self : TriangleMesh,
  options : MeshRelaxOptions,
) -> TriangleMesh {
  let mut current = self
  for _ in 0.. 0 {
          let mut average = Vec2::zero()
          for neighbor in neighbors {
            average = average.add(current.vertices[neighbor])
          }
          average = average.scale(1.0 / neighbors.length().to_double())
          next[index] = current.vertices[index].lerp(average, options.amount)
        }
      }
    }
    current = { ..current, vertices: next }
  }
  current
}

///|
fn midpoint_index(
  mesh : TriangleMesh,
  vertices : Array[Vec2],
  cache : Array[((Int, Int), Int)],
  left : Int,
  right : Int,
) -> Int {
  let edge = canonical_edge(left, right)
  for entry in cache {
    if entry.0 == edge {
      return entry.1
    }
  }
  let point = mesh.vertices[left].lerp(mesh.vertices[right], 0.5)
  vertices.push(point)
  let index = vertices.length() - 1
  cache.push((edge, index))
  index
}

///|
/// Splits every triangle into four midpoint triangles.
pub fn TriangleMesh::subdivide_once(self : TriangleMesh) -> TriangleMesh {
  let vertices = self.vertices.copy()
  let indices : Array[(Int, Int, Int)] = []
  let cache : Array[((Int, Int), Int)] = []
  for triangle in self.indices {
    let ab = midpoint_index(self, vertices, cache, triangle.0, triangle.1)
    let bc = midpoint_index(self, vertices, cache, triangle.1, triangle.2)
    let ca = midpoint_index(self, vertices, cache, triangle.2, triangle.0)
    indices.push((triangle.0, ab, ca))
    indices.push((ab, triangle.1, bc))
    indices.push((ca, bc, triangle.2))
    indices.push((ab, bc, ca))
  }
  { vertices, indices }
}

///|
/// Applies midpoint subdivision for a bounded number of levels.
pub fn TriangleMesh::subdivide(
  self : TriangleMesh,
  levels : Int,
) -> TriangleMesh {
  let mut result = self
  for _ in 0.. Array[Int] {
  let result = []
  for edge in self.boundary_edges() {
    if !result.contains(edge.left) {
      result.push(edge.left)
    }
    if !result.contains(edge.right) {
      result.push(edge.right)
    }
  }
  result
}

///|
/// Converts every edge into a line path.
pub fn TriangleMesh::mesh_wireframe(self : TriangleMesh) -> Array[Path] {
  let result = []
  for edge in self.edges() {
    if edge.left < self.vertices.length() && edge.right < self.vertices.length() {
      let left = self.vertices[edge.left]
      let right = self.vertices[edge.right]
      result.push(Path::new().move_to(left.x, left.y).line_to(right.x, right.y))
    }
  }
  result
}

///|
/// Converts boundary edges into a line path collection.
pub fn TriangleMesh::mesh_boundary_paths(self : TriangleMesh) -> Array[Path] {
  let result = []
  for edge in self.boundary_edges() {
    let left = self.vertices[edge.left]
    let right = self.vertices[edge.right]
    result.push(Path::new().move_to(left.x, left.y).line_to(right.x, right.y))
  }
  result
}

///|
/// Renders a mesh wireframe with a consistent style.
pub fn mesh_wireframe_document(
  mesh : TriangleMesh,
  width : Double,
  height : Double,
  background : Color,
  style : SvgStyle,
) -> SvgDocument {
  let document = SvgDocument::new(width, height)
  document.add(
    SvgElement::SvgPath(
      Path::rect(width, height),
      SvgStyle::new(fill=background.to_hex(), stroke="none"),
    ),
  )
  for path in mesh.mesh_wireframe() {
    document.add(SvgElement::SvgPath(path, style))
  }
  document
}

///|
/// Computes a scalar value at every mesh centroid.
pub fn mesh_centroid_values(
  mesh : TriangleMesh,
  values : Array[Double],
) -> Array[Double] {
  let result = []
  for index in 0.. Array[Color] {
  let result = []
  for value in mesh_centroid_values(mesh, values) {
    result.push(gradient.at(value))
  }
  result
}

///|
/// Fills each triangle with a separate gradient-derived color.
pub fn mesh_color_document(
  mesh : TriangleMesh,
  values : Array[Double],
  gradient : Gradient,
  width : Double,
  height : Double,
) -> SvgDocument {
  let document = SvgDocument::new(width, height)
  for index, triangle in mesh.indices {
    let polygon = [
      mesh.vertices[triangle.0],
      mesh.vertices[triangle.1],
      mesh.vertices[triangle.2],
    ]
    document.add(
      SvgElement::SvgPath(
        polyline(polygon, closed=true),
        SvgStyle::new(
          fill=mesh_triangle_colors(mesh, values, gradient)[index].to_hex(),
          stroke="none",
        ),
      ),
    )
  }
  document
}

///|
/// Reports the maximum vertex displacement between two compatible meshes.
pub fn mesh_displacement(left : TriangleMesh, right : TriangleMesh) -> Double {
  let mut maximum = 0.0
  let count = left.vertices.length().min(right.vertices.length())
  for index in 0.. String {
  "vertices,triangles,edges,area,min_area,max_area,mean_area,min_angle,degenerate,boundary\n" +
  "\{quality.vertex_count},\{quality.triangle_count},\{quality.edge_count},\{quality.total_area},\{quality.minimum_area},\{quality.maximum_area},\{quality.mean_area},\{quality.minimum_angle},\{quality.degenerate_triangles},\{quality.boundary_edges}\n"
}

///|
/// Creates a fan mesh with optional radial refinement.
pub fn refined_fan_mesh(
  center : Vec2,
  boundary : Array[Vec2],
  levels : Int,
) -> TriangleMesh {
  fan_mesh(center, boundary).subdivide(levels)
}

///|
/// Creates a normalized grid mesh and applies a scalar height function.
pub fn displaced_grid_mesh(
  columns : Int,
  rows : Int,
  width : Double,
  height : Double,
  displacement : (Double, Double) -> Double,
) -> TriangleMesh {
  let mesh = grid_mesh(columns, rows, width, height)
  for index, vertex in mesh.vertices {
    mesh.vertices[index] = Vec2::new(
      vertex.x,
      vertex.y + displacement(vertex.x, vertex.y),
    )
  }
  mesh
}

///|
/// Returns a deterministic checksum for mesh topology and coordinates.
pub fn mesh_signature(mesh : TriangleMesh) -> Int {
  let mut value = mesh.vertices.length() * 31 + mesh.indices.length()
  for vertex in mesh.vertices {
    value = (value * 31 + vertex.x.to_int() * 7 + vertex.y.to_int()) %
      1000000007
  }
  for triangle in mesh.indices {
    value = (value * 31 + triangle.0 * 7 + triangle.1 * 11 + triangle.2 * 13) %
      1000000007
  }
  value
}