///|
/// Computes perpendicular distance from point p to line segment (a, b).
pub fn perpendicular_distance(p : Point2D, a : Point2D, b : Point2D) -> Double {
  let dx = b.x - a.x
  let dy = b.y - a.y
  let len_sq = dx * dx + dy * dy
  if len_sq < 1.0e-12 {
    return p.distance(a)
  }

  let num = ((b.y - a.y) * p.x - (b.x - a.x) * p.y + b.x * a.y - b.y * a.x).abs()
  num / len_sq.sqrt()
}

///|
/// Recursive helper for Ramer-Douglas-Peucker (RDP) algorithm.
fn rdp_helper(
  pts : Array[Point2D],
  start_idx : Int,
  end_idx : Int,
  epsilon : Double,
  out : Array[Point2D],
) -> Unit {
  if end_idx <= start_idx + 1 {
    return
  }

  let mut max_dist = 0.0
  let mut max_idx = start_idx

  let a = pts[start_idx]
  let b = pts[end_idx]

  for i in (start_idx + 1).. max_dist {
      max_dist = d
      max_idx = i
    }
  }

  if max_dist > epsilon {
    rdp_helper(pts, start_idx, max_idx, epsilon, out)
    out.push(pts[max_idx])
    rdp_helper(pts, max_idx, end_idx, epsilon, out)
  }
}

///|
/// Simplifies a polyline using Ramer-Douglas-Peucker (RDP) algorithm.
pub fn rdp_simplify(pts : Array[Point2D], epsilon : Double) -> Array[Point2D] {
  let n = pts.length()
  if n <= 2 {
    return pts.copy()
  }

  let out = Array::new()
  out.push(pts[0])
  rdp_helper(pts, 0, n - 1, epsilon, out)
  out.push(pts[n - 1])
  out
}

///|
/// Checks if a 4-point polygon is strictly convex.
pub fn is_quad_convex(quad : Quad) -> Bool {
  if quad.corners.length() != 4 {
    return false
  }

  let mut sign_pos = false
  let mut sign_neg = false

  for i in 0..<4 {
    let p0 = quad.corners[i]
    let p1 = quad.corners[(i + 1) % 4]
    let p2 = quad.corners[(i + 2) % 4]

    let dx1 = p1.x - p0.x
    let dy1 = p1.y - p0.y
    let dx2 = p2.x - p1.x
    let dy2 = p2.y - p1.y

    let cross = dx1 * dy2 - dy1 * dx2
    if cross > 1.0e-6 {
      sign_pos = true
    }
    if cross < -1.0e-6 {
      sign_neg = true
    }
  }

  !(sign_pos && sign_neg)
}

///|
/// Reorders 4 corners of a quad into canonical order: TL (0), TR (1), BR (2), BL (3).
pub fn order_quad_corners(pts : Array[Point2D]) -> Quad {
  if pts.length() != 4 {
    return { corners: pts.copy() }
  }

  // Sort by sum/diff of coordinates to find corners
  let mut tl_idx = 0
  let mut tr_idx = 0
  let mut br_idx = 0
  let mut bl_idx = 0

  let mut min_sum = 1.0e18
  let mut max_sum = -1.0e18
  let mut min_diff = 1.0e18
  let mut max_diff = -1.0e18

  for i in 0..<4 {
    let sum = pts[i].x + pts[i].y
    let diff = pts[i].y - pts[i].x

    if sum < min_sum {
      min_sum = sum
      tl_idx = i
    }
    if sum > max_sum {
      max_sum = sum
      br_idx = i
    }
    if diff < min_diff {
      min_diff = diff
      tr_idx = i
    }
    if diff > max_diff {
      max_diff = diff
      bl_idx = i
    }
  }

  let ordered = [pts[tl_idx], pts[tr_idx], pts[br_idx], pts[bl_idx]]
  let quad = { corners: ordered }

  // Ensure top-left corner is top-most / left-most
  if is_quad_convex(quad) {
    quad
  } else {
    { corners: pts.copy() }
  }
}

///|
/// Filters contours and extracts valid convex quadrilateral candidates.
pub fn extract_quad_candidates(
  contours : Array[Contour],
  min_area : Double,
  min_perimeter : Double,
) -> Array[Quad] {
  let quads = Array::new()

  for i in 0..= min_area && perim >= min_perimeter {
      let epsilon = perim * 0.03
      let simplified = rdp_simplify(c.points, epsilon)

      if simplified.length() == 4 || simplified.length() == 5 {
        let four_pts = if simplified.length() == 5 {
          [simplified[0], simplified[1], simplified[2], simplified[3]]
        } else {
          simplified
        }

        let quad = order_quad_corners(four_pts)
        if is_quad_convex(quad) {
          quads.push(quad)
        }
      }
    }
  }

  quads
}