///|
/// Computes perpendicular distance from point p to line segment (a, b).
pub fn point_to_segment_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 = (dy * p.x - dx * p.y + b.x * a.y - b.y * a.x).abs()
  num / len_sq.sqrt()
}

///|
/// Visvalingam-Whyatt area-based polygon simplification algorithm.
/// Progressively eliminates vertices forming the triangle of minimum effective area.
pub fn visvalingam_whyatt_simplify(
  points : Array[Point2D],
  min_area_threshold : Double,
) -> Array[Point2D] {
  let n = points.length()
  if n <= 3 {
    return points.copy()
  }

  let pts = points.copy()
  let active = Array::make(n, true)
  let mut active_count = n

  // Triangle area function
  let tri_area = fn(p0 : Point2D, p1 : Point2D, p2 : Point2D) -> Double {
    let a = p0.x * (p1.y - p2.y) + p1.x * (p2.y - p0.y) + p2.x * (p0.y - p1.y)
    (a * 0.5).abs()
  }

  while active_count > 3 {
    let mut min_area = 1.0e30
    let mut min_idx = -1

    for i in 0.. min_area_threshold || min_idx == -1 {
      break
    }

    active[min_idx] = false
    active_count -= 1
  }

  let res : Array[Point2D] = []
  for i in 0.. Array[Point2D] {
  let mut current = points.copy()
  for _iter in 0.. (Point2D, Double) {
  let n = points.length()
  if n == 0 {
    return ({ x: 0.0, y: 0.0 }, 0.0)
  }
  if n == 1 {
    return (points[0], 0.0)
  }
  if n == 2 {
    let p0 = points[0]
    let p1 = points[1]
    let center = Point2D::new((p0.x + p1.x) * 0.5, (p0.y + p1.y) * 0.5)
    let radius = p0.distance(p1) * 0.5
    return (center, radius)
  }

  // Exact circumcircle of 3 points
  let circle_from_3 = fn(
    a : Point2D,
    b : Point2D,
    c : Point2D,
  ) -> (Point2D, Double)? {
    let d = 2.0 * (a.x * (b.y - c.y) + b.x * (c.y - a.y) + c.x * (a.y - b.y))
    if d.abs() < 1.0e-12 {
      return None // Collinear
    }
    let a_sq = a.x * a.x + a.y * a.y
    let b_sq = b.x * b.x + b.y * b.y
    let c_sq = c.x * c.x + c.y * c.y

    let ux = (a_sq * (b.y - c.y) + b_sq * (c.y - a.y) + c_sq * (a.y - b.y)) / d
    let uy = (a_sq * (c.x - b.x) + b_sq * (a.x - c.x) + c_sq * (b.x - a.x)) / d
    let center = Point2D::new(ux, uy)
    let radius = center.distance(a)
    Some((center, radius))
  }

  // Welzl iterative algorithm with bounding point sets
  let mut center = Point2D::new(
    (points[0].x + points[1].x) * 0.5,
    (points[0].y + points[1].y) * 0.5,
  )
  let mut radius = points[0].distance(points[1]) * 0.5

  for i in 2.. radius + 1.0e-9 {
      // Point outside, recompute circle with p on boundary
      center = Point2D::new(
        (points[0].x + p.x) * 0.5,
        (points[0].y + p.y) * 0.5,
      )
      radius = points[0].distance(p) * 0.5
      for j in 1.. radius + 1.0e-9 {
          center = Point2D::new((p.x + q.x) * 0.5, (p.y + q.y) * 0.5)
          radius = p.distance(q) * 0.5
          for k in 0.. radius + 1.0e-9 {
              match circle_from_3(p, q, r) {
                Some((c, rad)) => {
                  center = c
                  radius = rad
                }
                None => ()
              }
            }
          }
        }
      }
    }
  }

  (center, radius)
}

///|
/// Computes shape circularity (4 * pi * Area / Perimeter^2). Circle = 1.0, Square ~ 0.785.
pub fn shape_circularity(area : Double, perimeter : Double) -> Double {
  if perimeter <= 1.0e-12 {
    0.0
  } else {
    4.0 * 3.141592653589793 * area / (perimeter * perimeter)
  }
}