///|
pub struct CalibrationObservation {
  object : Point
  image : Point
  weight : Int
} derive(Debug, Eq)

///|
pub fn CalibrationObservation::new(
  object~ : Point,
  image~ : Point,
  weight? : Int = 1,
) -> CalibrationObservation {
  { object, image, weight: max_i(1, weight) }
}

///|
pub struct CalibrationSet {
  name : String
  observations : Array[CalibrationObservation]
} derive(Debug, Eq)

///|
pub fn CalibrationSet::new(name~ : String) -> CalibrationSet {
  { name, observations: [] }
}

///|
pub fn CalibrationSet::add(
  self : CalibrationSet,
  observation : CalibrationObservation,
) -> Unit {
  self.observations.push(observation)
}

///|
pub fn CalibrationSet::count(self : CalibrationSet) -> Int {
  self.observations.length()
}

///|
pub fn CalibrationSet::object_bounds(self : CalibrationSet) -> Rect {
  let points : Array[Point] = []
  for observation in self.observations {
    points.push(observation.object)
  }
  Polygon::new(points).bounds()
}

///|
pub fn CalibrationSet::image_bounds(self : CalibrationSet) -> Rect {
  let points : Array[Point] = []
  for observation in self.observations {
    points.push(observation.image)
  }
  Polygon::new(points).bounds()
}

///|
pub struct CalibrationEstimate {
  transform : Affine2D
  observations : Int
  error_x1000 : Int
  valid : Bool
} derive(Debug, Eq)

///|
pub fn CalibrationEstimate::new(
  transform~ : Affine2D,
  observations~ : Int,
  error_x1000~ : Int,
  valid~ : Bool,
) -> CalibrationEstimate {
  { transform, observations, error_x1000, valid }
}

///|
pub fn estimate_affine_translation(set : CalibrationSet) -> CalibrationEstimate {
  if set.observations.length() == 0 {
    CalibrationEstimate::new(
      transform=Affine2D::identity(),
      observations=0,
      error_x1000=0,
      valid=false,
    )
  } else {
    let mut weighted_dx = 0
    let mut weighted_dy = 0
    let mut weight = 0
    for observation in set.observations {
      weighted_dx += (observation.image.x - observation.object.x) *
        observation.weight
      weighted_dy += (observation.image.y - observation.object.y) *
        observation.weight
      weight += observation.weight
    }
    let transform = Affine2D::new(
      a=1000,
      b=0,
      c=0,
      d=1000,
      tx=weighted_dx / max_i(1, weight),
      ty=weighted_dy / max_i(1, weight),
    )
    let error = mean_reprojection_error(set, transform)
    CalibrationEstimate::new(
      transform~,
      observations=set.observations.length(),
      error_x1000=error,
      valid=set.observations.length() >= 3,
    )
  }
}

///|
pub fn mean_reprojection_error(
  set : CalibrationSet,
  transform : Affine2D,
) -> Int {
  if set.observations.length() == 0 {
    0
  } else {
    let mut total = 0
    let mut weight = 0
    for observation in set.observations {
      let projected = transform.apply(observation.object)
      total += projected.distance_l1(observation.image) * observation.weight
      weight += observation.weight
    }
    total * 1000 / max_i(1, weight)
  }
}

///|
pub fn calibration_residuals(
  set : CalibrationSet,
  transform : Affine2D,
) -> Array[Point] {
  let out : Array[Point] = []
  for observation in set.observations {
    let predicted = transform.apply(observation.object)
    out.push(
      Point::new(
        x=observation.image.x - predicted.x,
        y=observation.image.y - predicted.y,
      ),
    )
  }
  out
}

///|
pub fn calibration_quality(
  set : CalibrationSet,
  transform : Affine2D,
) -> String {
  let error = mean_reprojection_error(set, transform)
  let bounds = set.image_bounds()
  let sb = StringBuilder()
  sb.write_string("{\"observations\":\{set.count()},")
  sb.write_string("\"mean_l1_error_x1000\":\{error},\"image_bounds\":{")
  sb.write_string("\"x\":\{bounds.x},\"y\":\{bounds.y},")
  sb.write_string("\"width\":\{bounds.width},\"height\":\{bounds.height}}}")
  sb.to_string()
}

///|
pub struct Projective2D {
  h11 : Int
  h12 : Int
  h13 : Int
  h21 : Int
  h22 : Int
  h23 : Int
  h31 : Int
  h32 : Int
  denominator : Int
} derive(Debug, Eq)

///|
pub fn Projective2D::identity() -> Projective2D {
  {
    h11: 1000,
    h12: 0,
    h13: 0,
    h21: 0,
    h22: 1000,
    h23: 0,
    h31: 0,
    h32: 0,
    denominator: 1000,
  }
}

///|
pub fn Projective2D::new(
  h11~ : Int,
  h12~ : Int,
  h13~ : Int,
  h21~ : Int,
  h22~ : Int,
  h23~ : Int,
  h31~ : Int,
  h32~ : Int,
  denominator? : Int = 1000,
) -> Projective2D {
  { h11, h12, h13, h21, h22, h23, h31, h32, denominator: max_i(1, denominator) }
}

///|
pub fn Projective2D::apply(self : Projective2D, point : Point) -> Point {
  let denominator = self.h31 * point.x + self.h32 * point.y + self.denominator
  if denominator == 0 {
    point
  } else {
    Point::new(
      x=(self.h11 * point.x + self.h12 * point.y + self.h13) / denominator,
      y=(self.h21 * point.x + self.h22 * point.y + self.h23) / denominator,
    )
  }
}

///|
pub fn Projective2D::inverse_approx(
  self : Projective2D,
  point : Point,
) -> Point {
  let mut estimate = point
  for _ in 0..<5 {
    let projected = self.apply(estimate)
    estimate = estimate.translate(point.x - projected.x, point.y - projected.y)
  }
  estimate
}

///|
pub fn transform_box_projective(
  transform : Projective2D,
  box : Box2D,
  size : Size,
) -> Box2D {
  let points = [
    transform.apply(Point::new(x=box.rect.x, y=box.rect.y)),
    transform.apply(Point::new(x=box.rect.right() - 1, y=box.rect.y)),
    transform.apply(Point::new(x=box.rect.x, y=box.rect.bottom() - 1)),
    transform.apply(Point::new(x=box.rect.right() - 1, y=box.rect.bottom() - 1)),
  ]
  Box2D::new(
    label=box.label,
    rect=clamp_rect(Polygon::new(points).bounds(), size),
    occluded=box.occluded,
  )
}

///|
pub fn GrayImage::warp_projective(
  self : GrayImage,
  transform : Projective2D,
  fill? : Gray = Gray::black(),
) -> GrayImage {
  let out = GrayImage::new(width=self.width, height=self.height, fill~)
  for y in 0.. Int {
  if points.length() == 0 || size.area() <= 0 {
    0
  } else {
    Polygon::new(points).bounds().intersection(mask_bounds(size)).area() *
    1000 /
    size.area()
  }
}

///|
pub fn validate_calibration_grid(
  points : Array[Point],
  cols : Int,
  rows : Int,
  minimum_spacing : Int,
) -> ValidationReport {
  let report = ValidationReport::new()
  let expected = max_i(0, cols) * max_i(0, rows)
  if points.length() != expected {
    report.error("calibration point count does not match grid dimensions")
  }
  if minimum_spacing < 1 {
    report.error("minimum calibration spacing must be positive")
  }
  for i in 1..