说明
包围体是一个简单的几何空间，里面包含着复杂形状的物体。为物体添加包围体的目的是快速的进行碰撞检测或者进行精确的碰撞检测之前进行过滤（即当包围体碰撞，才进行精确碰撞检测和处理）。包围体类型包括球体、轴对齐包围盒（AABB）、有向包围盒（OBB）等。
大致区别可参考下面的网络图片，OBB 盒子的包裹更紧密。
在前面的文章《S01E22：利用矩阵 SVD 分解，拟合直线与平面》中，我们讲过矩阵 SVD 的几何意义和代码实现。但是 SVD 分解不仅可以用来拟合直线和平面，还可以用来计算 OBB 包围盒。
几何
首先，我们知道，SVD 分解与最小二乘法在数学上是等价的。那也就是说，我们分解后得到的特征向量，实际上就是最小二乘法拟合出的直线。在 3D 中也是这样，我们可以认为，特征值最大的特征向量（一般默认第一个）就是点云分散最广的方向；第二个则是垂直于第一个的情况下，点云分散最广的方向；第三个在 3D 中，直接垂直于前两个方向。
下面这个图中的例子比较夸张，但很能说明几何含义：全部的点分布在一个平面上，红色的点 M 是它们的几何平均点（相加再除以总数量）。显然分散最广的方向是红色箭头（ x 轴方向），在垂直于 x 轴的情况下分散最广的是绿色箭头方向（y 轴方向）。
这三个向量代表了点的分布情况，用来确定 OBB 盒子的朝向没有问题，但问题是原有的几何中心点并不是我们想要的 OBB 盒子中心点。我们真正需要的是如下图中的红色 Center 点，它是各个轴上最大值和最小值的平均值。（由于原点附近的点比较多，几何中心点 M 显然也会靠近点密集的地方，但这并不是我们想要的中心点）
所以当我们用 SVD 得到三个方向后，需要将所有点转换到新的坐标系中，求得各个轴的最大值和最小值，然后得到三边长及中心点，最后再将其转换回原始的世界坐标系中。这样就得到了 OBB 包围盒的真正中心点和三边长。
这里我们还可以顺便估算一下包围球的大小，也就是用 OBB 的中心点当做球心，用 OBB 对角线当成球的直径。估算的意思，也就是说：我们可以拿 OBB 包围盒的外接球，来顶替一下真正的包围球。一般来说，AABB 包围盒的外接球计算最简单，结果会稍大一些；而真正的最小包围球计算复杂，结果是最小的。
代码
extension Collection where Self.Element == simd_float3 {
///SVD 拟合直线、平面、外接球、包围盒
func estimateSVD() -> (line:Line?, plane:Plane?, sphere:Sphere?, boundingBox:simd_float4x4?) {
if self.count < 3 {
return (nil, nil, nil, nil)
}
var source:[Float] = Array(repeating: 0, count: self.count*3)
var position = simd_float3.zero
for point in self {
position += (point / Float(self.count))
}
for (row,point) in self.enumerated() {
source[row] = point.x - position.x
source[row+self.count] = point.y - position.y
source[row+2*self.count] = point.z - position.z
}
// source 中数据顺序为[x0,x1,x2.....y0,y1,y2.....z0,z1,z3....]，竖向按列依次填充到 rowCount行*3列 的矩阵中
let ss = Matrix(source: source, rowCount: self.count, columnCount:3)
let svdResult = Matrix.svd(a: ss)
let vt = svdResult.vt
// let sigma = svdResult.sigma
var offset = 0
let direction = simd_float3(vt[offset], vt[offset+vt.columnCount], vt[offset+2*vt.columnCount])
let line = Line(position: position, direction: direction)
offset = (vt.rowCount - 1)
let normal = simd_float3(vt[offset], vt[offset+vt.columnCount], vt[offset+2*vt.columnCount])
let plane = Plane(position: position, normal: normal)
offset = 1
let xAxis = direction
let yAxis = simd_float3(vt[offset], vt[offset+vt.columnCount], vt[offset+2*vt.columnCount])
let zAxis = normal
let matrix = simd_float4x4([simd_float4(xAxis, 0),
simd_float4(yAxis, 0),
simd_float4(zAxis, 0),
simd_float4(position, 1)])
var maxEdge = simd_float4(-Float.infinity, -Float.infinity, -Float.infinity, 1)
var minEdge = simd_float4(Float.infinity, Float.infinity, Float.infinity, 1)
for point in self {
let p = matrix.inverse * simd_float4(point, 1)
maxEdge = simd.max(maxEdge, p)
minEdge = simd.min(minEdge, p)
}
let centerPoint = matrix * (maxEdge + minEdge)/2
let lengthBox = maxEdge - minEdge
let centerP = simd_float3(centerPoint.x, centerPoint.y, centerPoint.z)
let radius = simd_length(simd_float3(lengthBox.x, lengthBox.y, lengthBox.z)) / 2
let sphere = Sphere(position: centerP, radius: radius)
let boundingBox = simd_float4x4([simd_float4(xAxis*lengthBox.x, 0),
simd_float4(yAxis*lengthBox.y, 0),
simd_float4(zAxis*lengthBox.z, 0),
simd_float4(centerP, 1)])
return (line:line, plane:plane, sphere:sphere, boundingBox:boundingBox)
}
}
