// 有向包围盒（OBB）。
//
// 对应 trimesh 的 `mesh.bounding_box_oriented`。
//
// 轴对齐包围盒（AABB）好算，但对斜着摆的模型特别不准：一根斜放的长棍，
// AABB 会是一个巨大的立方体，而 OBB 能贴着棍子包上去。
//
// 做法是主成分分析（PCA）：
//
//   1. 求所有顶点的重心和协方差矩阵
//   2. 协方差矩阵的特征向量就是「顶点分布最开的三个方向」
//   3. 把顶点投影到这三个方向上，取各自的极值
//
// PCA 不保证最小体积；重复特征值下轴向不唯一。
// 通过一次 Jacobi 分解取得完整正交基，再包围输入顶点。

///|
/// 有向包围盒。
pub(all) struct Obb {
  /// 盒子的中心
  center : Vec3
  /// 三个轴，已单位化。按对应边的长度从大到小排
  axes : Array[Vec3]
  /// 沿三个轴的半长，和 `axes` 一一对应
  half : Array[Double]
}

///|
/// 三条边的长度。
pub fn Obb::size(self : Obb) -> Vec3 {
  Vec3::new(self.half[0] * 2.0, self.half[1] * 2.0, self.half[2] * 2.0)
}

///|
/// 盒子的体积。真实体积的**上界**——模型装在里面，一般填不满。
pub fn Obb::volume(self : Obb) -> Double {
  self.half[0] * 8.0 * self.half[1] * self.half[2]
}

///|
/// 求有向包围盒。空网格返回 `None`。
pub fn Mesh::oriented_bounding_box(self : Mesh) -> Obb? {
  let pts = unique_positions_exact(self)
  if pts.length() == 0 {
    return None
  }

  // 局部坐标中补偿求和，并归一化协方差以避免平方溢出/下溢。
  let origin = pts[0]
  let cx = MomentSum::new()
  let cy = MomentSum::new()
  let cz = MomentSum::new()
  let mut scale = 0.0
  for p in pts {
    if !finite_vertex(p) {
      return None
    }
    let d = p.sub(origin)
    if !finite_vertex(d) {
      return None
    }
    scale = scale.max(d.x.abs()).max(d.y.abs()).max(d.z.abs())
  }
  let scale = if scale > 0.0 { scale } else { 1.0 }
  for p in pts {
    let d = p.sub(origin)
    cx.add(d.x / scale)
    cy.add(d.y / scale)
    cz.add(d.z / scale)
  }
  let n = pts.length().to_double()
  let local_center = Vec3::new(cx.value() / n, cy.value() / n, cz.value() / n)
  let centroid = Vec3::new(
    origin.x + local_center.x * scale,
    origin.y + local_center.y * scale,
    origin.z + local_center.z * scale,
  )

  // 协方差矩阵。对称，所以 9 个数里只有 6 个独立
  let mut xx = 0.0
  let mut yy = 0.0
  let mut zz = 0.0
  let mut xy = 0.0
  let mut yz = 0.0
  let mut zx = 0.0
  for p in pts {
    let d = p.sub(origin)
    let dx = d.x / scale - local_center.x
    let dy = d.y / scale - local_center.y
    let dz = d.z / scale - local_center.z
    xx = xx + dx * dx
    yy = yy + dy * dy
    zz = zz + dz * dz
    xy = xy + dx * dy
    yz = yz + dy * dz
    zx = zx + dz * dx
  }
  let cov : FixedArray[Double] = FixedArray::make(9, 0.0)
  cov[0] = xx / n
  cov[1] = xy / n
  cov[2] = zx / n
  cov[3] = xy / n
  cov[4] = yy / n
  cov[5] = yz / n
  cov[6] = zx / n
  cov[7] = yz / n
  cov[8] = zz / n

  // 一次分解返回完整正交基；重复特征值不能分别取零空间向量。
  let axes = symmetric_eigendecomposition(cov).vectors

  // 把顶点投影到三个轴上，各取极值。
  // 投影的基准取重心——选哪里都行，取重心数值上稳一些
  let lo : Array[Double] = [0.0, 0.0, 0.0]
  let hi : Array[Double] = [0.0, 0.0, 0.0]
  for i, p in pts {
    for k = 0; k < 3; k = k + 1 {
      let shifted = p.sub(origin)
      let local_point = Vec3::new(
        shifted.x / scale - local_center.x,
        shifted.y / scale - local_center.y,
        shifted.z / scale - local_center.z,
      )
      let d = local_point.dot(axes[k]) * scale
      if not_finite(d) {
        return None
      }
      if i == 0 || d < lo[k] {
        lo[k] = d
      }
      if i == 0 || d > hi[k] {
        hi[k] = d
      }
    }
  }

  // 盒子中心 = 重心 + 各轴上的中点偏移
  let mut center = centroid
  let half : Array[Double] = []
  for k = 0; k < 3; k = k + 1 {
    let mid = lo[k] / 2.0 + hi[k] / 2.0
    center = Vec3::new(
      center.x + axes[k].x * mid,
      center.y + axes[k].y * mid,
      center.z + axes[k].z * mid,
    )
    let radius = hi[k] / 2.0 - lo[k] / 2.0
    if not_finite(radius) {
      return None
    }
    half.push(radius)
  }

  // 按边长从大到小排。主轴按特征值排是「分布最开的方向优先」，
  // 但用户看的是「盒子有多长多宽多高」，按长度排更直观
  let order = sort_axes_by_length(half)
  let sorted_axes : Array[Vec3] = []
  let sorted_half : Array[Double] = []
  for k in order {
    sorted_axes.push(axes[k])
    sorted_half.push(half[k])
  }
  Some({ center, axes: sorted_axes, half: sorted_half, })
}

///|
/// 三个轴按半长从大到小排序，返回下标顺序。
///
/// 只有三个元素，直接手写比较比调排序省事
fn sort_axes_by_length(half : Array[Double]) -> Array[Int32] {
  let idx : Array[Int32] = [0, 1, 2]
  // 冒泡两轮，从大到小
  for _pass = 0; _pass < 2; _pass = _pass + 1 {
    for i = 0; i + 1 < 3; i = i + 1 {
      if half[idx[i]] < half[idx[i + 1]] {
        let t = idx[i]
        idx[i] = idx[i + 1]
        idx[i + 1] = t
      }
    }
  }
  idx
}

///|
/// 去重后的顶点位置。
///
/// 几何统计用精确坐标去重，不让默认拓扑焊接容差压扁微小模型。
fn unique_positions_exact(mesh : Mesh) -> Array[Vec3] {
  let table : Map[PositionKey, Int32] = Map([])
  let out : Array[Vec3] = []
  for t in mesh.triangles {
    for v in [t.a, t.b, t.c] {
      let key = position_key(v, 0.0)
      match table.get(key) {
        Some(_) => ()
        None => {
          table.set(key, out.length())
          out.push(v)
        }
      }
    }
  }
  out
}