阐明

前面咱们讲过《S03E11： 用 SVD 分解创建点云的 OBB 围住盒》，不仅可以求出 OBB 围住盒，还可以用围住盒的外接球来近似替代最小围住球。但事实上，可以用 Welzl算法 来直接求出准确的最小围住球。

S03E14： 求多个点的最小围住球（Welzl算法）

几何

Welzl算法 很容易用递归算法来实现：

  • 先依据 np-1 个点，生成一个球体（递归）；
  • 判别第 np 个点是否在球体内，是，则坚持球体；
  • 否，需要依据旧的球体和第 np 个点，生成新的球体（递归）。其间第 np 个点必在新生成的球面上；
  • 递归结束条件：当一个点、两个点时，可直接生成球体；三个点都在球面上，则生成三角形的外接球；

三个点都在球面上，其实有个条件：第三个点不在前两个点生成的球体内，第二个点不在一三点生成的球体内，第一个点不在二三点生成的球体内。在递归算法中，这是必然可以保证的，故此时外接球便是最小围住球。从几何角度来看，便是说三个点组成的三角形不能为钝角三角形，由于关于钝角三角形，直接拿最长边当直径，可以得到更小的球体。

S03E14： 求多个点的最小围住球（Welzl算法）
关于三角形外接圆（外接球）的圆心和半径，可以参看《S02E16：三角形外心（外接圆圆心）》中的办法。 伪代码如下：

// pt：一切的点集；np：要处理的前 np 个点；bnd：已确定在球面上的点，初始为空
func minSphere(pt:[simd_float3], np:Int, bnd:[simd_float3] = []) -> Sphere {
    if 满足递归中止条件 {
         return 球面上的点得到的球体（依据一个点、两个点或三个点，生成一个球体）
    }
   // 先求 np-1 个点的最小围住球
    let D = minSphere(pt: pt, np: np-1, bnd: bnd)
    if 第 np 个点在球内  {
        return D
    }
    // 将第 np 个点参加球面集合 bnd 中
    var newbnd = bnd
    newbnd.append(pt[np-1])
    // 依据前 np-1 个点的最小围住球，与第 np 个点，生成新的球体
    return minSphere(pt: pt, np: np-1, bnd: newbnd)
}

代码

///多个点的最小围住球，Welzl算法
///http://www.sunshine2k.de/coding/java/Welzl/Welzl.html
static func minSphere(points:[simd_float3]) -> Sphere {
    return minSphere(pt: points, np: points.count)
}
private static func minSphere(pt:[simd_float3], np:Int, bnd:[simd_float3] = []) -> Sphere {
    if np == 1 {
        if bnd.isEmpty {
            return sphere1pt(pt[0])
        } else if bnd.count == 1 {
            return sphere2pts(p1: pt[0], p2: bnd[0])
        }
    } else if np == 0 {
        if bnd.isEmpty {
            return sphere1pt(.zero)
        } else if bnd.count == 1 {
            return sphere1pt(bnd[0])
        } else if bnd.count == 2 {
            return sphere2pts(p1: bnd[0], p2: bnd[1])
        }
    }
    if bnd.count == 3 {
        return sphere3pts(p1: bnd[0], p2: bnd[1], p3: bnd[2])
    }
    let D = minSphere(pt: pt, np: np-1, bnd: bnd)
    if distanceBetween(point: pt[np-1], sphere: D) <= 0  {
        return D
    }
    var newbnd = bnd
    newbnd.append(pt[np-1])
    return minSphere(pt: pt, np: np-1, bnd: newbnd)
}
private static func sphere1pt(_ p:simd_float3) ->Sphere {
    return Sphere(position: p, radius: 0)
}
private static func sphere2pts(p1:simd_float3, p2:simd_float3) ->Sphere {
    return Sphere(position: (p1+p2)/2, radius: simd_distance(p1, p2)/2)
}
///三个点的外接球
private static func sphere3pts(p1:simd_float3, p2:simd_float3,p3:simd_float3) ->Sphere {
    let t = Triangle(points: float3x3(p1, p2, p3))
    let r = Triangle.circumcenterRadius(triangle: t)
    let c = Triangle.circumcenter(triangle: t)
    return Sphere(position: c, radius: r)
}