3D 最小二乘平面
Posted
技术标签:
【中文标题】3D 最小二乘平面【英文标题】：3D Least Squares Plane 【发布时间】：2010-11-26 21:04:54 【问题描述】：给定一组 3D 数据点，在 (x, y, z) 空间中计算最小二乘平面的算法是什么？换句话说，如果我有一堆点，如 (1, 2, 3), (4, 5, 6), (7, 8, 9) 等，如何计算最佳拟合平面 f (x, y) = ax + by + c?从一组 3D 点中获取 a、b 和 c 的算法是什么？
【问题讨论】：
您应该定义“最小二乘”的确切含义。参见例如en.wikipedia.org/wiki/Least_squares 用于定义它的各种方式。 这是一条评论。如果有人愿意将其转移到斯蒂芬佳能的答案中，那就太好了。我希望这能澄清他所说的“解向量的三个分量是最小二乘拟合平面 a,b,c 的系数”到底是什么意思。首先，给定 Ax = b 是初等矩阵代数，其中 A 是一个矩阵，b 和 x 是向量，只有当 A 有一个非零行列式。这意味着您需要超过 3 个点，并且至少有一个点不应该在平面上（不准确）。第二个澄清，以及原因 【参考方案1】：如果您有 n 个数据点 (x[i], y[i], z[i])，请计算 3x3 对称矩阵 A，其条目为：
sum_i x[i]*x[i], sum_i x[i]*y[i], sum_i x[i]
sum_i x[i]*y[i], sum_i y[i]*y[i], sum_i y[i]
sum_i x[i], sum_i y[i], n
同时计算 3 元素向量 b：
sum_i x[i]*z[i], sum_i y[i]*z[i], sum_i z[i]
然后对给定的 A 和 b 求解 Ax = b。解向量的三个分量是最小二乘拟合平面 a,b,c 的系数。
请注意，这是“普通最小二乘法”拟合，仅当预期 z 是 x 和 y 的线性函数时才适用。如果您更一般地在 3 空间中寻找“最佳拟合平面”，您可能需要了解“几何”最小二乘法。
另请注意，如果您的点在一条线上，就像您的示例点一样，这将失败。
【讨论】：
常用的方法有3种。您提到的基于 SVD 的方法对于某些问题是首选，但比我使用的相当基本的“正规方程”更难解释（和理解）。基于 SVD 的方法实际上比求解正规方程慢，但由于它在数值上更稳定，因此更受欢迎。 如果您想了解所有不同的 LLS 算法，en.wikipedia.org/wiki/… 有更多详细信息。 @AKE：不，当你使用 QR 或 SVD 时，你不使用正规方程（意味着你不形成我描述的 3x3 矩阵，而是直接在 nx3 矩阵上操作测量）。这导致更好的稳定性 (a) 因为正交分解具有良好的稳定性属性，更重要的是 (b) 因为 XtX 的条件数是 X 的条件数的平方。 我在哪里可以找到这种方法的解释，特别是矩阵的制作方式以及为什么解向量等于平面的系数？我看到了 Wikipedia 文章，但它涵盖了对我没有帮助的非常笼统的方程式。 我想你可以在这篇文章中找到一些解释ilikebigbits.com/blog/2015/3/2/plane-from-points【参考方案2】：平面的方程是：ax + by + c = z。所以用你的所有数据设置这样的矩阵：
x_0 y_0 1
A = x_1 y_1 1
...
x_n y_n 1
和
a
x = b
c
和
z_0
B = z_1
...
z_n
换句话说：Ax = B。现在求解x，它们是你的系数。但是因为（我假设）你有超过 3 个点，系统是超定的，所以你需要使用左伪逆。所以答案是：
a
b = (A^T A)^-1 A^T B
c
下面是一些简单的 Python 代码和示例：
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import numpy as np
N_POINTS = 10
TARGET_X_SLOPE = 2
TARGET_y_SLOPE = 3
TARGET_OFFSET = 5
EXTENTS = 5
NOISE = 5
# create random data
xs = [np.random.uniform(2*EXTENTS)-EXTENTS for i in range(N_POINTS)]
ys = [np.random.uniform(2*EXTENTS)-EXTENTS for i in range(N_POINTS)]
zs = []
for i in range(N_POINTS):
zs.append(xs[i]*TARGET_X_SLOPE + \
ys[i]*TARGET_y_SLOPE + \
TARGET_OFFSET + np.random.normal(scale=NOISE))
# plot raw data
plt.figure()
ax = plt.subplot(111, projection='3d')
ax.scatter(xs, ys, zs, color='b')
# do fit
tmp_A = []
tmp_b = []
for i in range(len(xs)):
tmp_A.append([xs[i], ys[i], 1])
tmp_b.append(zs[i])
b = np.matrix(tmp_b).T
A = np.matrix(tmp_A)
fit = (A.T * A).I * A.T * b
errors = b - A * fit
residual = np.linalg.norm(errors)
print "solution:"
print "%f x + %f y + %f = z" % (fit[0], fit[1], fit[2])
print "errors:"
print errors
print "residual:"
print residual
# plot plane
xlim = ax.get_xlim()
ylim = ax.get_ylim()
X,Y = np.meshgrid(np.arange(xlim[0], xlim[1]),
np.arange(ylim[0], ylim[1]))
Z = np.zeros(X.shape)
for r in range(X.shape[0]):
for c in range(X.shape[1]):
Z[r,c] = fit[0] * X[r,c] + fit[1] * Y[r,c] + fit[2]
ax.plot_wireframe(X,Y,Z, color='k')
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('z')
plt.show()
【讨论】：
如果出现“奇异矩阵”错误：xs = [] ys = [] zs = [] for i in range(N_POINTS): xs.append(i * i) ys.append(0) zs.append(1) @Alexey，您的输入点在一条直线上。所以有无限的平面同样适合数据。这是一个退化的情况，所以最小二乘解决方案不起作用。实际上，我不确定任何通用解决方案都可以使用它。如果您有具体问题，请打开一个新问题，可能会参考这个问题。【参考方案3】：除非有人告诉我如何在这里输入方程式，否则让我写下你必须做的最终计算：
首先，给定点r_i \n \R，i=1..N，计算所有点的质心：
r_G = \frac\sum_i=1^N r_iN
然后，计算法向量n，它与基向量r_G一起通过计算3x3矩阵A来定义平面
A = \sum_i=1^N (r_i - r_G)(r_i - r_G)^T
有了这个矩阵，法线向量 n 现在由 A 的特征向量给出，对应于 A 的最小特征值。
要了解特征向量/特征值对，请使用您选择的任何线性代数库。
此解决方案基于 Hermitian 矩阵 A 的 Rayleight-Ritz 定理。
【讨论】：
你能告诉我你从哪里得到的吗？我看不出与 Rayleight-Ritz 定理的联系？？ @josch 这里没有 LaTeX。见meta.stackexchange.com/questions/30559/latex-on-stack-overflow【参考方案4】：请参阅 David Eberly 的“数据的最小二乘拟合”，了解我是如何提出这个以最小化几何拟合（从点到平面的正交距离）的。
bool Geom_utils::Fit_plane_direct(const arma::mat& pts_in, Plane& plane_out)
bool success(false);
int K(pts_in.n_cols);
if(pts_in.n_rows == 3 && K > 2) // check for bad sizing and indeterminate case
plane_out._p_3 = (1.0/static_cast<double>(K))*arma::sum(pts_in,1);
arma::mat A(pts_in);
A.each_col() -= plane_out._p_3; //[x1-p, x2-p, ..., xk-p]
arma::mat33 M(A*A.t());
arma::vec3 D;
arma::mat33 V;
if(arma::eig_sym(D,V,M))
// diagonalization succeeded
plane_out._n_3 = V.col(0); // in ascending order by default
if(plane_out._n_3(2) < 0)
plane_out._n_3 = -plane_out._n_3; // upward pointing
success = true;
return success;
时间为 37 微秒将平面拟合到 1000 个点（Windows 7、i7、32 位程序）
【讨论】：
【参考方案5】：与任何最小二乘法一样，您可以这样进行：
开始编码之前
在参数化(a, b, d) 中写下一个平面方程，例如0 = ax + by + z + d。
找到一个表达式D(\vecv;a, b, d)，表示到任意点\vecv的距离。
写下 S = \sigma_i=0,n D^2(\vecx_i) 的总和，并进行化简，直到用 v 的各分量的简单总和表示，例如 \sigma v_x、\sigma v_y^2、\sigma v_x*v_z ...
写下每个参数的最小化表达式dS/dx_0 = 0、dS/dy_0 = 0 ... 这会为您提供一组三个参数中的三个方程以及上一步的总和。
求解这组方程的参数。
（或者对于简单的情况，只需查找表单）。使用符号代数包（如 Mathematica）可以让您的生活更轻松。
编码
编写代码以形成所需的总和并从上面最后一组中找到参数。替代方案
请注意，如果您实际上只有三个点，您最好只找到穿过它们的平面。
此外，如果解析解不可行（不是飞机的情况，但一般情况下可能），您可以执行步骤 1 和 2，并在步骤 3 的总和上使用 Monte Carlo minimizer。
【讨论】：
【参考方案6】：CGAL::linear_least_squares_fitting_3
函数 linear_least_squares_fitting_3 计算最佳拟合 3D 一组 3D 对象的线或平面（在最小二乘意义上），例如 作为点、线段、三角形、球体、球、长方体或四面体。
http://www.cgal.org/Manual/latest/doc_html/cgal_manual/Principal_component_analysis_ref/Function_linear_least_squares_fitting_3.html
【讨论】：
【参考方案7】：这简化为Total Least Squares 问题，可以使用SVD 分解来解决。
使用 OpenCV 的 C++ 代码：
float fitPlaneToSetOfPoints(const std::vector<cv::Point3f> &pts, cv::Point3f &p0, cv::Vec3f &nml)
const int SCALAR_TYPE = CV_32F;
typedef float ScalarType;
// Calculate centroid
p0 = cv::Point3f(0,0,0);
for (int i = 0; i < pts.size(); ++i)
p0 = p0 + conv<cv::Vec3f>(pts[i]);
p0 *= 1.0/pts.size();
// Compose data matrix subtracting the centroid from each point
cv::Mat Q(pts.size(), 3, SCALAR_TYPE);
for (int i = 0; i < pts.size(); ++i)
Q.at<ScalarType>(i,0) = pts[i].x - p0.x;
Q.at<ScalarType>(i,1) = pts[i].y - p0.y;
Q.at<ScalarType>(i,2) = pts[i].z - p0.z;
// Compute SVD decomposition and the Total Least Squares solution, which is the eigenvector corresponding to the least eigenvalue
cv::SVD svd(Q, cv::SVD::MODIFY_A|cv::SVD::FULL_UV);
nml = svd.vt.row(2);
// Calculate the actual RMS error
float err = 0;
for (int i = 0; i < pts.size(); ++i)
err += powf(nml.dot(pts[i] - p0), 2);
err = sqrtf(err / pts.size());
return err;
【讨论】：
conv 函数是什么？这不是 OpenCV 之一，对吧？
对，我不认为是opencv。只是Point3f到Vec3f的转换，和减去cv::Point3f(0,0,0)一样。【参考方案8】：
听起来您想要做的就是使用 2 个回归器进行线性回归。关于这个主题的wikipedia page 应该告诉你所有你需要知道的，然后是一些。
【讨论】：
【参考方案9】：你所要做的就是解方程组。
如果这些是你的观点： (1, 2, 3), (4, 5, 6), (7, 8, 9)
这给了你方程式：
3=a*1 + b*2 + c
6=a*4 + b*5 + c
9=a*7 + b*8 + c
所以你的问题实际上应该是：我如何求解方程组？
因此我建议阅读this SO question。
如果我误解了您的问题，请告诉我们。
编辑：
忽略我的回答，因为你可能是别的意思。
【讨论】：
OPs 的例子很糟糕，因为他给出了三分，但他想解决给定 n 分的一般情况，因此是一个过约束系统。如果不是所有的点都在一个平面上，他想找到最佳拟合，即最小二乘法中所有点到平面的距离最小化的平面。以上是关于3D 最小二乘平面的主要内容，如果未能解决你的问题，请参考以下文章