///|
pub(all) struct BoundingBox {
  minimum_x : Double
  maximum_x : Double
  minimum_y : Double
  maximum_y : Double
} derive(Debug, Eq, ToJson)

///|
pub fn bounding_box(
  points : ArrayView[Point],
) -> BoundingBox raise @domain.ExtractionError {
  if points.length() == 0 {
    raise @domain.ExtractionError::InsufficientData(required=1, actual=0)
  }
  let mut minimum_x = points[0].x
  let mut maximum_x = points[0].x
  let mut minimum_y = points[0].y
  let mut maximum_y = points[0].y
  for point in points {
    minimum_x = minimum_x.min(point.x)
    maximum_x = maximum_x.max(point.x)
    minimum_y = minimum_y.min(point.y)
    maximum_y = maximum_y.max(point.y)
  }
  { minimum_x, maximum_x, minimum_y, maximum_y }
}

///|
pub fn BoundingBox::width(self : BoundingBox) -> Double {
  self.maximum_x - self.minimum_x
}

///|
pub fn BoundingBox::height(self : BoundingBox) -> Double {
  self.maximum_y - self.minimum_y
}

///|
pub fn BoundingBox::area(self : BoundingBox) -> Double {
  self.width().max(0.0) * self.height().max(0.0)
}

///|
pub fn BoundingBox::contains(self : BoundingBox, point : Point) -> Bool {
  point.x >= self.minimum_x &&
  point.x <= self.maximum_x &&
  point.y >= self.minimum_y &&
  point.y <= self.maximum_y
}

///|
pub fn translate_polygon(
  points : ArrayView[Point],
  delta_x~ : Double,
  delta_y~ : Double,
) -> Array[Point] raise @domain.ExtractionError {
  @domain.require_finite("translation x", delta_x)
  @domain.require_finite("translation y", delta_y)
  [
    for point in points => Point::new(x=point.x + delta_x, y=point.y + delta_y)
  ]
}

///|
pub fn scale_polygon(
  points : ArrayView[Point],
  scale_x~ : Double,
  scale_y~ : Double,
  origin? : Point = Point::{ x: 0.0, y: 0.0 },
) -> Array[Point] raise @domain.ExtractionError {
  @domain.require_finite("polygon x scale", scale_x)
  @domain.require_finite("polygon y scale", scale_y)
  [
    for point in points => {
      Point::new(
        x=origin.x + scale_x * (point.x - origin.x),
        y=origin.y + scale_y * (point.y - origin.y),
      )
    }
  ]
}

///|
pub fn rotate_point(
  point : Point,
  angle_sine~ : Double,
  angle_cosine~ : Double,
  origin? : Point = Point::{ x: 0.0, y: 0.0 },
) -> Point raise @domain.ExtractionError {
  @domain.require_finite("rotation sine", angle_sine)
  @domain.require_finite("rotation cosine", angle_cosine)
  let x = point.x - origin.x
  let y = point.y - origin.y
  Point::new(
    x=origin.x + angle_cosine * x - angle_sine * y,
    y=origin.y + angle_sine * x + angle_cosine * y,
  )
}

///|
pub fn rotate_polygon(
  points : ArrayView[Point],
  angle_sine~ : Double,
  angle_cosine~ : Double,
  origin? : Point = Point::{ x: 0.0, y: 0.0 },
) -> Array[Point] raise @domain.ExtractionError {
  [
    for point in points => {
      rotate_point(point, angle_sine~, angle_cosine~, origin~)
    }
  ]
}

///|
pub fn polygon_perimeter(
  points : ArrayView[Point],
) -> Double raise @domain.ExtractionError {
  if points.length() < 2 {
    raise @domain.ExtractionError::InsufficientData(
      required=2,
      actual=points.length(),
    )
  }
  let mut total = 0.0
  for index, point in points {
    total += point.distance(points[(index + 1) % points.length()])
  }
  total
}

///|
pub fn is_convex(points : ArrayView[Point]) -> Bool {
  if points.length() < 3 {
    false
  } else {
    let mut sign = 0
    for index in 0.. 1.0e-12 {
        sign += 1
      }
      if cross < -1.0e-12 {
        sign -= 1
      }
    }
    sign == points.length() || sign == -points.length()
  }
}

///|
pub fn signed_turns(points : ArrayView[Point]) -> Array[Double] {
  [
    for index in 0.. {
      let a = points[index]
      let b = points[(index + 1) % points.length()]
      let c = points[(index + 2) % points.length()]
      (b.x - a.x) * (c.y - b.y) - (b.y - a.y) * (c.x - b.x)
    }
  ]
}

///|
fn distance_to_segment(
  point : Point,
  first : Point,
  second : Point,
) -> (Double, Point) {
  let dx = second.x - first.x
  let dy = second.y - first.y
  let length_squared = dx * dx + dy * dy
  let fraction = if length_squared <= 1.0e-20 {
    0.0
  } else {
    ((point.x - first.x) * dx + (point.y - first.y) * dy) / length_squared
  }
  let clamped = fraction.max(0.0).min(1.0)
  let projected = Point::{
    x: first.x + clamped * dx,
    y: first.y + clamped * dy,
  }
  (point.distance(projected), projected)
}

///|
pub(all) struct NearestBoundaryPoint {
  point : Point
  distance : Double
  segment : Int
} derive(Debug, Eq, ToJson)

///|
pub fn nearest_boundary_point(
  polygon : ArrayView[Point],
  point : Point,
) -> NearestBoundaryPoint raise @domain.ExtractionError {
  if polygon.length() < 2 {
    raise @domain.ExtractionError::InsufficientData(
      required=2,
      actual=polygon.length(),
    )
  }
  let mut best_distance = 1.0e300
  let mut best_point = polygon[0]
  let mut best_segment = 0
  for index, start in polygon {
    let (distance, projection) = distance_to_segment(
      point,
      start,
      polygon[(index + 1) % polygon.length()],
    )
    if distance < best_distance {
      best_distance = distance
      best_point = projection
      best_segment = index
    }
  }
  { point: best_point, distance: best_distance, segment: best_segment }
}

///|
pub(all) struct LineIntersection {
  intersects : Bool
  point : Point
  first_fraction : Double
  second_fraction : Double
} derive(Debug, Eq, ToJson)

///|
pub fn line_intersection(
  first_start : Point,
  first_end : Point,
  second_start : Point,
  second_end : Point,
  tolerance? : Double = 1.0e-12,
) -> LineIntersection raise @domain.ExtractionError {
  @domain.require_positive("line intersection tolerance", tolerance)
  let dx1 = first_end.x - first_start.x
  let dy1 = first_end.y - first_start.y
  let dx2 = second_end.x - second_start.x
  let dy2 = second_end.y - second_start.y
  let denominator = dx1 * dy2 - dy1 * dx2
  if denominator.abs() <= tolerance {
    {
      intersects: false,
      point: first_start,
      first_fraction: 0.0,
      second_fraction: 0.0,
    }
  } else {
    let dx = second_start.x - first_start.x
    let dy = second_start.y - first_start.y
    let first_fraction = (dx * dy2 - dy * dx2) / denominator
    let second_fraction = (dx * dy1 - dy * dx1) / denominator
    let point = Point::{
      x: first_start.x + first_fraction * dx1,
      y: first_start.y + first_fraction * dy1,
    }
    {
      intersects: first_fraction >= 0.0 &&
      first_fraction <= 1.0 &&
      second_fraction >= 0.0 &&
      second_fraction <= 1.0,
      point,
      first_fraction,
      second_fraction,
    }
  }
}

///|
pub fn clip_polygon_half_plane(
  polygon : ArrayView[Point],
  line_start : Point,
  line_end : Point,
  keep_left? : Bool = true,
) -> Array[Point] raise @domain.ExtractionError {
  if polygon.length() < 3 {
    raise @domain.ExtractionError::InsufficientData(
      required=3,
      actual=polygon.length(),
    )
  }
  let output : Array[Point] = []
  let side = (point : Point) => {
    (line_end.x - line_start.x) * (point.y - line_start.y) -
    (line_end.y - line_start.y) * (point.x - line_start.x)
  }
  let inside = (point : Point) => {
    if keep_left {
      side(point) >= 0.0
    } else {
      side(point) <= 0.0
    }
  }
  for index, current in polygon {
    let previous = polygon[(index + polygon.length() - 1) % polygon.length()]
    let current_inside = inside(current)
    let previous_inside = inside(previous)
    if current_inside && previous_inside {
      output.push(current)
    } else if current_inside && !previous_inside {
      let crossing = line_intersection(previous, current, line_start, line_end)
      if crossing.intersects ||
        (crossing.first_fraction >= 0.0 && crossing.first_fraction <= 1.0) {
        output.push(crossing.point)
      }
      output.push(current)
    } else if !current_inside && previous_inside {
      let crossing = line_intersection(previous, current, line_start, line_end)
      if crossing.intersects ||
        (crossing.first_fraction >= 0.0 && crossing.first_fraction <= 1.0) {
        output.push(crossing.point)
      }
    }
  }
  output
}

///|
pub fn polygon_compact(
  points : ArrayView[Point],
  tolerance~ : Double,
) -> Array[Point] raise @domain.ExtractionError {
  @domain.require_positive("polygon compact tolerance", tolerance)
  if points.length() < 3 {
    points.to_owned()
  } else {
    let result : Array[Point] = [points[0]]
    for index in 1.. tolerance {
        result.push(points[index])
      }
    }
    if result.length() > 1 &&
      result[0].distance(result[result.length() - 1]) <= tolerance {
      ignore(result.pop())
    }
    result
  }
}

///|
pub fn polygon_overlap_ratio(
  polygon : ArrayView[Point],
  box : BoundingBox,
) -> Double raise @domain.ExtractionError {
  let clipped_left = clip_polygon_half_plane(
    polygon,
    Point::{ x: box.minimum_x, y: box.minimum_y },
    Point::{ x: box.minimum_x, y: box.maximum_y },
  )
  let clipped_right = clip_polygon_half_plane(
    clipped_left,
    Point::{ x: box.maximum_x, y: box.maximum_y },
    Point::{ x: box.maximum_x, y: box.minimum_y },
  )
  let clipped_bottom = clip_polygon_half_plane(
    clipped_right,
    Point::{ x: box.minimum_x, y: box.minimum_y },
    Point::{ x: box.maximum_x, y: box.minimum_y },
    keep_left=false,
  )
  let clipped = clip_polygon_half_plane(
    clipped_bottom,
    Point::{ x: box.maximum_x, y: box.maximum_y },
    Point::{ x: box.minimum_x, y: box.maximum_y },
    keep_left=false,
  )
  if clipped.length() < 3 {
    0.0
  } else {
    polygon_area(clipped) / box.area()
  }
}