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

wblong_cs

1487人浏览 · 2025-02-24 21:49:54

wblong_cs · 2025-02-24 21:49:54 发布

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

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

原理

给定一组三维空间中的点 $,pn}\mathbf{P} = \{\mathbf{p}_1, \mathbf{p}_2, \cdots, \mathbf{p}_n\}$ ，其中 $pi=[xi,yi,zi]T\mathbf{p}_i = [x_i, y_i, z_i]^T$ 。我们的目标是找到一个平面方程 $z = a x + b y + c$ 来拟合这些点。

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

计算点集的质心： $c=1n∑i=1npi\mathbf{c} = \frac{1}{n}\sum_{i = 1}^{n}\mathbf{p}_i$ 。
将点集进行中心化： $qi=pi−c\mathbf{q}_i = \mathbf{p}_i - \mathbf{c}$ ， $\cdots, n$ 。
构建数据矩阵： $[\mathbf{q}_1, \mathbf{q}_2, \cdots, \mathbf{q}_n]^T$ 。
对数据矩阵进行奇异值分解： $U\Sigma V^T$ ，其中 $U$ 是 $n×nn\times n$ 的正交矩阵， $Σ\Sigma$ 是 $n×3n\times 3$ 的对角矩阵， $V$ 是 $3×33\times 3$ 的正交矩阵。
确定平面的法向量： $V$ 的最后一列对应的向量即为平面的法向量 $n=[a,b,c]T\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()