SVD_0">奇异值分解（SVD）拟合平面

在三维空间中，使用奇异值分解（SVD）来拟合平面是一种常见且有效的方法。下面将详细介绍其原理、实现步骤，并给出Python代码示例。

原理

给定一组三维空间中的点 P = { p 1 , p 2 , ⋯   , p n } \mathbf{P} = \{\mathbf{p}_1, \mathbf{p}_2, \cdots, \mathbf{p}_n\} P={p1​,p2​,⋯,pn​}，其中 ${\mathbf{p}}_i=[x_i,y_i,z_i{]}^T$。我们的目标是找到一个平面方程 $z=ax+by+c$ 来拟合这些点。

可以通过最小化点到平面的距离平方和来实现平面拟合。具体步骤如下：

1. 计算点集的质心：$\mathbf{c}=\frac{1}{n}{\sum }_{i=1}^n{\mathbf{p}}_i$。
2. 将点集进行中心化：${\mathbf{q}}_i={\mathbf{p}}_i-\mathbf{c}$，$i=1,2,\cdots ,n$。
3. 构建数据矩阵：$A=[{\mathbf{q}}_1,{\mathbf{q}}_2,\cdots ,{\mathbf{q}}_n{]}^T$。
4. 对数据矩阵进行奇异值分解：$A=U\Sigma V^T$，其中 $U$ 是 $n\times n$ 的正交矩阵，$\Sigma$ 是 $n\times 3$ 的对角矩阵，$V$ 是 $3\times 3$ 的正交矩阵。
5. 确定平面的法向量：$V$ 的最后一列对应的向量即为平面的法向量 $\mathbf{n}=[a,b,c{]}^T$。

Python代码实现

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

def plane_fit(x, y, z):
    """
    拟合平面 z = Ax + By + C 到给定的 (x, y, z) 数据
    :param x: 一维数组，x 坐标
    :param y: 一维数组，y 坐标
    :param z: 一维数组，z 坐标
    :return: A, B, C 平面方程的系数
    """
    # 计算数据点的质心
    P = np.array([np.mean(x), np.mean(y), np.mean(z)])
    # 数据中心化
    centered_data = np.column_stack((x - P[0], y - P[1], z - P[2]))
    # 进行奇异值分解
    U, S, Vt = np.linalg.svd(centered_data)
    # 获取平面的法向量
    print(Vt)
    N =-1/Vt[-1,-1]* Vt[-1, :] 
    print(N)
    # 提取平面方程的系数
    A = N[0]
    B = N[1]
    C = -np.dot(P, N)
    print(A,B,C)
    return A, B, C

# 示例代码
if __name__ == "__main__":
    # 生成网格数据
    x_grid, y_grid = np.meshgrid(np.linspace(0, 10, 20), np.linspace(0, 10, 20))
    # 定义平面方程的真实参数
    a = 1
    b = 2
    c = -2
    # 计算理论 z 值
    z_grid = (a * x_grid) + (b * y_grid) + c
    # 将网格数据转换为一维数组
    x = x_grid.flatten()
    y = y_grid.flatten()
    z = z_grid.flatten()
    # 添加高斯噪声
    z = z + np.random.randn(len(z))
    # 拟合平面
    A, B, C = plane_fit(x, y, z)
    # 生成用于绘制拟合平面的网格
    X, Y = np.meshgrid(np.linspace(np.min(x), np.max(x), 20), np.linspace(np.min(y), np.max(y), 20))
    # 计算拟合平面的 z 值
    Z = (A * X) + (B * Y) + C
    # 创建 3D 图形
    fig = plt.figure()
    ax = fig.add_subplot(111, projection='3d')
    # 绘制原始数据点
    ax.scatter(x, y, z, c='r', marker='.')
    # 绘制拟合平面
    ax.plot_surface(X, Y, Z, color='g', alpha=0.5)
    # 设置坐标轴网格
    ax.grid(True)
    # 设置图形标题
    title = f' a={a:.4f},b={b:.4f},c={c:.4f},\nA={A:.4f},B={B:.4f},C={C:.4f}'
    ax.set_title(title)
    # 显示图形
    plt.show()

SVD拟合点云平面" />

代码解释

1. 计算质心：使用 np.mean 函数计算点集的质心。
2. 中心化点集：将每个点减去质心，得到中心化后的点集。
3. 奇异值分解：使用 np.linalg.svd 函数对数据矩阵进行奇异值分解。
4. 确定法向量：取 $V$ 的最后一列作为平面的法向量。

复杂度分析

• 时间复杂度：主要由奇异值分解的复杂度决定，为 $O(n{m}^2)$，其中 $n$ 是点的数量，$m=3$ 是点的维度。
• 空间复杂度：主要用于存储数据矩阵和奇异值分解的结果，为 $O(nm)$。