S01E22：利用矩阵 SVD 分解，拟合直线与平面在平面相关问题中，讲到将多个点拟合成一个平面的算法，缺点是需要注意点

说明

在平面相关问题中，讲到将多个点拟合成一个平面的算法，缺点是需要注意点的顺序。那么有没有办法拟合一堆无顺序的点呢？或者拟合成一条直线可以吗？

当然是可以的，其实用最小二乘法就可以解决些问题。但是，在 3D 中直接使用最小二乘法非常不方便。因此常见解法是，用投影方式，投影到两个平面上，再进行 2D 的求解，但这样也不方便，而且有时会出现拟合结果与平面垂直的情况，需要换个平面反复计算。

SVD 分解

奇异值分解(Singular Value Decomposition，以下简称SVD)就是解决最小二乘法的利器，它不仅可以拟合直线、平面，还可以得到点云的最小包围盒。关于 SVD 与最小二乘的数学原理和关联，可以直接网上搜索查找，资料大把。本文主要讲解其几何意义和代码实现。

首先，我们可以将一大堆需要拟合点的坐标(x,y,z)，组成一个 n 行三列的矩阵 A（需要提前减去平均值 $(x_0,y_0,z_0)$ ，即点云的中心与(0,0,0)对齐），然后进行 SVD 分解，得到三个矩阵 U𝚺 $V^T$ ，矩阵 V 的第一列（Column）即是直线的方向向量，而第三列就是平面的法线向量。

A = \left[ \begin{matrix} x_1-x_0 & y_1-y_0 & z_1-z_0 \\ x_2-x_0 & y_2-y_0 & z_2-z_0 \\ \vdots & \vdots & \vdots \\ x_n-x_0 & y_n-y_0 & z_n-z_0 \\ \end{matrix} \right]

其中矩阵 V 中的向量被称为特征向量，它是个 3x3 的正交矩阵；而矩阵 𝚺 中的值则是特征值，矩阵 𝚺 只有对角线上有值，因此也可以省略其它的 0，直接写在一个向量里。

注意： $V^T$ 是 V 的转置矩阵，因此矩阵V的第一列实际上是矩阵 $V^T$ 的第一行。另外，对于正交矩阵 V 来说，它的转置等于它的逆，也就是 $V^{-1} = V^T$ ，下文出现的 $V^*$ 也是同样意义。

几何

那么这个 SVD 分解几何意义是什么，我们又该怎么理解记忆呢？让我们先从二维上一个圆说起：假设一堆点，组成了一个椭圆（图右上角），如果我们要拟合直线，就需要找到这个椭圆的长轴。而矩阵 A 的含义就是，将一个标准圆，经过变换得到我们要拟合的椭圆。

那么根据矩阵乘法定义： $V^T$ 矩阵的含义就是，先将这个标准圆旋转，矩阵 𝚺 含义是将其缩放，矩阵 U 表示再旋转一下。

所以， $V^T$ 是一个正交矩阵，代表旋转；𝚺 是一个对角矩阵，代表拉伸缩放，因为只有对角线有值，所以可以直接写到一个向量里，以节省空间；U 也是一个正交矩阵，代表旋转。

在三维中，也是同样道理，只不过标准圆变成了一个球体，而椭圆则变成了椭球体。而平面的拟合，则是需要找到椭球体的最短轴（即最扁的方向）。

代码

那么SVD 代码如何写呢？最常见的是使用LAPACK，这是一个专门做线性代数的库。苹果的Accelerate框架直接内置它，只需要稍做封装就能直接使用。

苹果官方的示例程序Denoising an Image Using Linear Algebra则直接提供了 SVD 相关函数的封装，包括对任意维度矩阵的封装Matrix，以及 SVD 函数sgesdd_()，我们可以直接拿来用。

sgesdd_()是利用分治算法实现的 SVD 算法，sgesvd_()是普通版。详见www.netlib.org/lapack/lug/…

Type of problem	Function and storage scheme	Single precision Real	Single precision Complex	Double precision Real	Double precision Complex
SVD	simple driver	SGESVD	CGESVD	DGESVD	ZGESVD
	divide and conquer driver	SGESDD	CGESDD	DGESDD	ZGESDD

具体代码如下：

static func estimateLineSVD(from points:[simd_float3]) -> Line? {
    if points.count < 2 {
        return nil
    }
    var source:[Float] = Array(repeating: 0, count: points.count*3)
    var position = simd_float3.zero
    for point in points {
        position += (point / Float(points.count))
    }
    for (row,point) in points.enumerated() {
        source[row] = point.x - position.x
        source[row+points.count] = point.y - position.y
        source[row+2*points.count] = point.z - position.z
    }
    // source 中数据顺序为[x0,x1,x2.....y0,y1,y2.....z0,z1,z3....]，竖向按列依次填充到 rowCount行*3列 的矩阵中
    let ss = Matrix(source: source, rowCount: points.count, columnCount:3)
    let svdResult = Matrix.svd(a: ss)
    let vt = svdResult.vt
    let offset = 0
    // vt 中第一行，即 v 中第一列
    let direction = simd_float3(vt[offset], vt[offset+vt.columnCount], vt[offset+2*vt.columnCount])
    return Line(position: position, direction: direction)
}


static func estimatePlaneSVD(from points:[simd_float3]) -> Plane? {
    if points.count < 3 {
        return nil
    }
    var source:[Float] = Array(repeating: 0, count: points.count*3)
    var position = simd_float3.zero
    for point in points {
        position += (point / Float(points.count))
    }
    for (row,point) in points.enumerated() {
        source[row] = point.x - position.x
        source[row+points.count] = point.y - position.y
        source[row+2*points.count] = point.z - position.z
    }
    // source 中数据顺序为[x0,x1,x2.....y0,y1,y2.....z0,z1,z3....]，竖向按列依次填充到 rowCount行*3列 的矩阵中
    let ss = Matrix(source: source, rowCount: points.count, columnCount:3)
    let svdResult = Matrix.svd(a: ss)
    let vt = svdResult.vt
    let offset = (vt.rowCount - 1)
    // vt 中第三行
    let normal = simd_float3(vt[offset], vt[offset+vt.columnCount], vt[offset+2*vt.columnCount])
    return Plane(position: position, normal: normal)
}

其他

云算网-SVD 在线分解

Denoising an Image Using Linear Algebra