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

奇异值分解（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 Axes3Ddef 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 = 1b = 2c = -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()

代码解释

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

复杂度分析

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