///|
pub fn point_moments(points : Array[Point]) -> Array[Double] {
  let center = polygon_centroid(points)
  let mut m20 = 0.0
  let mut m02 = 0.0
  let mut m11 = 0.0
  let mut m30 = 0.0
  let mut m03 = 0.0
  let mut m21 = 0.0
  let mut m12 = 0.0
  for point in points {
    let x = point.x - center.x
    let y = point.y - center.y
    m20 += x * x
    m02 += y * y
    m11 += x * y
    m30 += x * x * x
    m03 += y * y * y
    m21 += x * x * y
    m12 += x * y * y
  }
  [m20, m02, m11, m30, m03, m21, m12]
}

///|
pub fn hu_like_descriptors(points : Array[Point]) -> Array[Double] {
  let moments = point_moments(points)
  if moments.length() < 7 {
    []
  } else {
    let first = moments[0] + moments[1]
    let second = (moments[0] - moments[1]) * (moments[0] - moments[1]) +
      4.0 * moments[2] * moments[2]
    let third = (moments[3] - 3.0 * moments[5]) *
      (moments[3] - 3.0 * moments[5]) +
      (3.0 * moments[6] - moments[4]) * (3.0 * moments[6] - moments[4])
    [first, second, third]
  }
}

///|
pub fn radial_profile(
  points : Array[Point],
  center : Point,
  bins : Int,
) -> Array[Double] {
  let result = Array::make(bins.max(1), 0.0)
  let counts = Array::make(bins.max(1), 0)
  let maximum = maximum_value(points.map(fn(point) { distance(point, center) })).max(
    0.0000001,
  )
  for point in points {
    let index = (distance(point, center) /
      maximum *
      Double::from_int(result.length() - 1))
      .to_int()
      .clamp(min=0, max=result.length() - 1)
    result[index] += distance(point, center)
    counts[index] += 1
  }
  for i in 0.. 0 {
      result[i] /= Double::from_int(counts[i])
    }
  }
  result
}

///|
pub fn angular_profile(
  points : Array[Point],
  center : Point,
  bins : Int,
) -> Array[Double] {
  let result = Array::make(bins.max(1), 0.0)
  let counts = Array::make(bins.max(1), 0)
  for point in points {
    let angle = normalize_angle(measure_angle(center, point))
    let index = ((angle + 3.141592653589793) /
      6.283185307179586 *
      Double::from_int(result.length()))
      .to_int()
      .clamp(min=0, max=result.length() - 1)
    result[index] += distance(point, center)
    counts[index] += 1
  }
  for i in 0.. 0 {
      result[i] /= Double::from_int(counts[i])
    }
  }
  result
}

///|
pub fn descriptor_distance(
  first : Array[Double],
  second : Array[Double],
) -> Double {
  if first.length() != second.length() {
    0.0
  } else {
    let mut total = 0.0
    for i in 0.. Double {
  let first_descriptor = hu_like_descriptors(first)
  let second_descriptor = hu_like_descriptors(second)
  1.0 / (1.0 + descriptor_distance(first_descriptor, second_descriptor))
}

///|
pub fn oriented_shape_similarity(
  first : Array[Point],
  second : Array[Point],
) -> Double {
  let first_angle = principal_orientation(first)
  let second_angle = principal_orientation(second)
  shape_similarity(
    rotate_points(first, first_angle),
    rotate_points(second, second_angle),
  )
}

///|
pub fn rotate_points(points : Array[Point], angle : Double) -> Array[Point] {
  let center = polygon_centroid(points)
  points.map(fn(point) { rotate_point(point, center, angle) })
}

///|
pub fn shape_extent(points : Array[Point]) -> Double {
  let bounds = polygon_bounds(points)
  bounding_width(bounds).max(bounding_height(bounds))
}

///|
pub fn shape_compactness(points : Array[Point]) -> Double {
  let metrics = polygon_metrics(points)
  if metrics.perimeter <= 0.0 {
    0.0
  } else {
    4.0 *
    3.141592653589793 *
    metrics.area /
    (metrics.perimeter * metrics.perimeter)
  }
}