当前位置：首页 > news >正文

python用最小二乘法实现平面拟合

news 2025/8/27 15:57:26

文章目录

- 数学原理
- 代码实现
- 测试

数学原理

平面方程可写为

$A x + B y + C z + D = 0$

假设 $C$ 不为0，则上式可以改写为

$z = a x + b y + d$

则现有一组点 ${p_i\}$ ，则根据 $x_i,y_i$ 以及平面方程，可以得到其对应的 $\hat z_i$

$\hat z_i=ax_i+by_i+d$

从而平面拟合就转换为了最小二乘问题

$\argmin \sum \vert z_i-\hat z_i\vert$

将其转换为矩阵形式，记

$A=\begin{bmatrix} x_1&y_1&1\\x_2&y_2&1\\\vdots&\vdots\\x_n&y_n&1 \end{bmatrix}, x=\begin{bmatrix}a\\b\\d\end{bmatrix}, b=\begin{bmatrix}z_1\\z_2\\\vdots\\z_n\end{bmatrix}$

则拟合方程变为

$A x = b$

相应地 $x$ 可写为

$x=(A^TA)^{-1}A^Tb$

最小二乘法的原理可见：python实现最小二乘法

代码实现

其代码实现如下

import numpy as np
def planefit(points):xs, ys, zs = list(zip(*points))A = np.vstack([xs,ys,np.ones_like(xs)]).Tb = np.reshape(zs, [-1, 1])abd = np.linalg.inv(A.T @ A) @ A.T @ breturn abd.reshape(-1)

其中输入参数points为一组点。第一步将这些点进行坐标拆分，得到一一对应的xs, ys, zs。然后通过这些点，构造矩阵A和向量b，最后输出 $A^TA)^{-1}A^Tb$ 。

测试

首先做一个初始化平面的函数，其功能是随机生成一组在平面中的点，并且为其添加一些噪声。

# 初始化平面
def initPlane(a, b, d, N, err=0.1):xs,ys = np.random.rand(2, N)zs = a*xs + b*ys + d + np.random.rand(N)*errreturn list(zip(xs, ys, zs))

然后用planefit函数对这些点进行拟合，通过对比二者之间的差异，来证实算法的有效性

import matplotlib.pyplot as pltpts = initPlane(2,3,4,100,1)
a,b,d = planefit(pts)xs, ys = np.indices([100,100])/100
zs = a*xs + b*ys + dax = plt.subplot(projection='3d')
ax.plot_surface(xs, ys, zs, cmap='jet')
ax.scatter(*np.array(pts).T, marker='*')
plt.show()

在这里插入图片描述

随着加入的噪声逐渐变大，拟合误差也越来越大

errs = [0.01, 0.1, 0.2, 0.5, 1, 3, 5]
fits = []
for err in errs:pts = initPlane(2,3,4,100,1)a,b,d = planefit(pts)fits.append([abs(a-2),abs(b-3),abs(d-4)])import pprint
pprint.pprint(fits)

[[0.09377971025135245, 0.023025216275622373, 0.4933931906466551],
[0.044310965250572654, 0.05681830483294226, 0.47952260969370997],
[0.051813469166934745, 0.017914573861143257, 0.47553046120193176],
[0.08578595894551588, 0.0464898508775029, 0.42791269232718054],
[0.011569662177250528, 0.15976404558136714, 0.4886516489062753],
[0.006829071411009302, 0.04832062421804073, 0.5193494695593301],
[0.1651263679674586, 0.0736367910618192, 0.44103226768552073]]

查看全文

http://www.lryc.cn/news/236503.html