三维空间中离散点直线拟合方法详解

三维空间中离散点直线拟合方法详解

引用

CSDN

1.

https://blog.csdn.net/qq_35635374/article/details/145533470

在计算机图形学中，根据离散点拟合直线是一个常见的问题。本文将介绍两种方法：高中教材中的最小二乘代数方法和优化最小二乘模型目标函数的方法，并给出详细的数学推导和代码实现。

问题描述

给定一系列的三维空间点，需要拟合得到迫近的直线方程。

方法一：高中教材中的最小二乘的代数方法

这种方法基于最小二乘原理，通过求解线性方程组来得到直线方程。具体来说，就是通过求解矩阵方程Ax=b，其中A是由点的坐标构成的矩阵，b是由点的坐标构成的向量，x是直线方程的系数向量。

但是这种方法在高维空间中计算量较大，且容易受到噪声点的影响，因此在实际应用中较少使用。

方法二：优化最小二乘模型目标函数的方法

这种方法通过优化最小二乘模型目标函数来求解直线方程。具体来说，就是通过求解矩阵S的最小特征值对应的特征向量来得到直线的方向向量v，然后通过计算所有点的坐标平均值来得到直线经过的一点L0。

这种方法在高维空间中计算量较小，且对噪声点的鲁棒性较好，因此在实际应用中较为常用。

（1）推导

设直线的点向式方程为：

由式(1)，得到直线的参数方程为：

式(2)写成向量形式为：

其中，

如下图，红色点Li(xi,yi,zi)为给定的一系列三维空间点，根据给定三维空间点，拟合直线方程(3)，也就是计算L0和v，使得在某种“距离”的度量下，达到最佳的直线拟合效果。

点Li到直线距离的平方为：

L0Li在直线的投影的平方为：

在最小二乘准则下，可以建立直线拟合的优化模型目标函数：

可以得到结论：待拟合的直线经过一个点L0，该点的坐标为所有给定点的坐标平均值。如下图所示，一旦确定直线的单位方向向量v，则直线的方程便确定。

f的最小值为矩阵S最小特征值对应的特征向量。直线方向向量v的求解问题转化为矩阵最小特征值对应的特征向量的求解问题！

（2）代码实现

%{
Function: line_fitting
Description: 直线拟合
Input: 任意维直线点数据points，行数为点个数，列数为点的维数
Output: 拟合得到的直线经过的一点L0,直线的单位方向向量v
Author: Marc Pony(marc_pony@163.com)
%}
function [L0, v] = line_fitting(points)
n = size(points, 1);
x = points(:, 1);
y = points(:, 2);
z = points(:, 3);
L0 = [mean(x); mean(y); mean(z)];
S = zeros(3,3);
for i = 1 : n
    Yi = [x(i) - L0(1); y(i) - L0(2); z(i) - L0(3)];
    S = S + (Yi' * Yi * eye(3, 3) - Yi * Yi');
end
[V, ~] = eig(S);
v = V(:, 1); %矩阵S最小特征值对应的特征向量
end

%{
Function: generate_line_points
Description: 直线路径点生成
Input: 直线经过的一点L0,直线的单位方向向量v,点个数n，路径标量最小值minS，路径标量最大值maxS
Output: 任意维直线点数据points，行数为点个数，列数为点的维数
Author: Marc Pony(marc_pony@163.com)
%}
function points = generate_line_points(L0, v, n, minS, maxS)
points = zeros(n, length(v));
s = linspace(minS, maxS, n);
for i = 1 : n
    points(i, :) = (L0 + v * s(i))';
end
end

clear
clc
close all
%% 验证恒等式: v'*Yi*v = v*v'*Yi
syms v1 v2 v3 y1 y2 y3 real
v = [v1; v2; v3];
Yi = [y1; y2; y3];
res1 = simplify(v'*Yi*v - v*v'*Yi)
%% 验证恒等式: Yi'*Yi = v'*(Yi'*Yi)*v, 其中v'*v=1
res2 = [Yi'*Yi; simplify(v'*(Yi'*Yi)*v)]
%% 验证恒等式: (v'*Yi)^2 = v'*(Yi*Yi')*v
res3 = simplify((v'*Yi)^2 - v'*(Yi*Yi')*v)
% points = [1 0 0
%     1 10 0
%     1 20 0
%     ];
% points = [0 1 0
%     10 1 0
%     200 1 0
%     ];
% points = [1 1 1
%     2 1 2
%     ];
figure
axis([-10, 10, -10, 10])
hold on
pointCount = 6;
points = zeros(pointCount, 3);
for i = 1 : pointCount
    [points(i, 1), points(i, 2)] = ginput(1);
    plot(points(i, 1), points(i, 2), '+')
end
[L0, v] = line_fitting(points)
n = 100;
len = sqrt((max(points(:,1)) - min(points(:,1)))^2 + (max(points(:,2)) - min(points(:,2)))^2 + (max(points(:,3)) - min(points(:,3)))^2);
minS = -0.6 * len;
maxS = 0.6 * len;
p = generate_line_points(L0, v, n, minS, maxS);
plot3(p(:,1), p(:,2), p(:,3), '-')

（3）效果

（此处应有拟合效果的示意图，但由于原文中未提供具体的效果图，因此无法展示。）