•  05-07-2019
  •  | 
  •  

题

给定一组 3D 数据点，计算 (x, y, z) 空间中的最小二乘平面的算法是什么？换句话说，如果我有一堆点，例如 (1, 2, 3), (4, 5, 6), (7, 8, 9) 等，那么如何计算最佳拟合平面 f (x, y) = ax + by + c？从一组 3D 点中获取 a、b 和 c 的算法是什么？

有帮助吗？

解决方案

如果你有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空间中，您可能想要了解“几何”。最小二乘。

另请注意，如果您的点在一行中，则会失败，因为您的示例点是。

其他提示

除非有人告诉我如何在这里输入方程式，否则让我写下你必须做的最终计算：

首先，给定点r_i \ n \ R，i = 1..N，计算所有点的质心：

r_G = \frac{\sum_{i=1}^N r_i}{N}

然后，计算法向量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定理。

参见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位程序）

与任何最小二乘法一样，您可以这样进行：

开始编码之前

  1. 写下一些参数化中的平面方程，例如 0 = ax + by + z + d 在你的参数中 (a, b, d).

  2. 找一个表达 D(\vec{v};a, b, d) 到任意点的距离 \vec{v}.

  3. 写下总和 S = \sigma_i=0,n D^2(\vec{x}_i), ，并简化，直到它用以下各分量的简单和来表示 v 喜欢 \sigma v_x, \sigma v_y^2, \sigma v_x*v_z ...

  4. 写下每个参数的最小化表达式 dS/dx_0 = 0, dS/dy_0 = 0 ...它为您提供了一组三个参数的三个方程以及上一步的总和。

  5. 求解这组方程的参数。

（或者对于简单的情况，只需查找表格即可）。使用符号代数包（如 Mathematica）可以让你的生活变得更轻松。

编码

  • 编写代码以形成所需的总和并从上面最后一组中找到参数。

备择方案

请注意，如果您实际上只有三个点，那么您最好只找到穿过它们的平面。

另外，如果解析解不可行（飞机不是这种情况，但一般来说是可能的），您可以执行步骤 1 和 2，并使用 蒙特卡洛最小化器 关于步骤 3 中的总和。

CGAL::linear_least_squares_fitting_3
  

函数linear_least_squares_fitting_3计算最佳拟合3D   一组3D对象的线或平面（在最小二乘意义上）   作为点，段，三角形，球体，球，长方体或四面体。

http://www.cgal.org/Manual/最新/ DOC_HTML / cgal_manual / Principal_component_analysis_ref / Function_linear_least_squares_fitting_3.html

平面的等式是：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()

这减少了最小二乘法问题，可以使用 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;
}

你所要做的就是解决方程组。

如果这些是你的观点： （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

所以你的问题实际应该是：我如何解决方程组？

因此，我建议您阅读这个问题。

如果我误解了你的问题，请告诉我们。

修改：

忽略我的回答，因为你可能意味着别的东西。

许可以下： CC-BY-SA 和 归因
不隶属于 StackOverflow
scroll top