通过Python中的正交回归拟合平面

问题描述

我想在Python中使平面适合一组点(x,y,z)。我发现了各种答案,如果相对于z轴测量了误差,该如何执行拟合,但是我想考虑正交方向的误差。我发现以下问题(

enter image description here

解决了相同的问题-但我不清楚如何在Python中实现此问题(可能使用NumPy / SciPy)。有关数学推导的更多详细信息,也可以在这里找到:Best fit plane by minimizing orthogonal distances(第2节)。

解决方法

您给出的第一个链接确实描述了正交距离拟合的算法,但是非常简洁。如果有帮助,这里是一个更复杂的描述:

我想您有个点(在您的情况下为3d,但是维数对算法没有影响)P [i],i = 1..N 您想要找到一个距点最小正交距离的(超)平面。

p可以用单位矢量n和标量d描述超平面。平面上的点集是
{ P | n.P + d = 0 }

,并且点P与平面的(正交)距离为

n.P + d

所以我们想找到n和d来最小化

Q(n,d) = Sum{ i | (n.P[i]+d)*(n.P[i]+d) } /N

(除以N并不是必须的,并且对找到的n和d的值没有影响,但在我看来使代数更整洁)

首先要注意的是,如果我们知道n,则使Q最小的d将是

d = -n.Pbar where
Pbar = Sum{ i | P[i]}/N,the mean of the P[]

我们也可以使用d的这个值,这样,在经过一点代数运算后,问题就会减少到最小化Q ^:

Q^(n) = Sum{ i | (n.P[i]-n.Pbar)*(n.P[i]-n.Pbar) } /N
      = n' * C * n
where
C = Sum{ i | (P[i]-Pbar)*(P[i]-Pbar) } /N

Q ^的形式告诉我们,使Q ^最小的n的值将是C的特征向量,对应于最小特征值。

所以(对不起,我不能提供代码,但是我的python是可鄙的):

a /计算

Pbar = Sum{ i | P[i]}/N,the mean of the points

b /计算

C = Sum{ i | (P[i]-Pbar)*(P[i]-Pbar) } /N,the covariance matrix of the points

c /对角线C,并选取一个最小特征值和相应的特征向量n

d /计算

d = -Pbar.n

然后n,d定义所需的超平面。

,

基于您的second refernce [enter image description here]

假设您有n个样本(x,y,z)

我将这3个术语称为M*A=V,并定义数组

X=[ x_0,x_1 .. x_n ]'
Y=[ y_0,y_1 .. y_n ]'
Z=[ z_0,z_1 .. z_n ]'

定义(n乘3)矩阵XY1=[X,Y,1n]

      [[x_0,y_0,1],XY1=   [x_1,y_1,...
       [x_n,y_n,1]]

矩阵M可以通过

获得
M = XY1' * XY1

撇号(')是换位运算符,(*)是矩阵乘积。

数组V

V = XY1'*Z

最小二乘解可以通过摩尔-penrose伪逆数获得:[(M'*M)^-1 * M']

~A  = [(M'*M)^-1 * M'] * V

示例代码:

import numpy as np
from mpl_toolkits import mplot3d
import matplotlib.pyplot as plt

#Input your values
A=3
B=2
C=1

#reserve memory
xy1=np.ones([n,3])

#Make random data,n  ( x,y ) tuples.
n=30 #samples
xy1[:,:2]=np.random.rand(n,2)

#plane: A*x+B*y+C = z,the z coord is calculated from random x,y
z=xy1.dot (np.array([[A,B,C],]).transpose() )
    
#addnoise
xy1[:,:2]+=np.random.normal(scale=0.05,size=[n,2])
z+=np.random.normal(scale=0.05,1])

#calculate M and V    
M=xy1.transpose().dot(xy1)
V=xy1.transpose().dot(z)

#pseudoinverse:
Mp=np.linalg.inv(M.transpose().dot(M)).dot(M.transpose())

#Least-squares Solution
ABC= Mp.dot(V) 

输出

In [24]: ABC
Out[24]: 
array([[3.11395111],[2.02909874],[1.01340411]])
,

我也不得不处理这种情况,起初数学符号可能会让人不知所措,但最终解决方案非常简单。

一旦您直觉到定义最佳拟合平面 Ax + By + Cz + D = 0 的向量(A,B,C)是一个解释了一组坐标的最小方差,那么解决方案很简单。

首先要做的是使坐标居中(这样,在平面方程中 D 将为0)

coords -= coords.mean(axis=0)

然后,您有2个选项来获取您感兴趣的向量:(1)使用sklearnscipy中的PCA实现来获取解释最小方差的向量

pca = PCA(n_components=3)
pca.fit(coords)
# The last component/vector is the one with minimal variance,see PCA documentation
normal_vector = pca.components_[-1]

(2)重新实现已链接的“几何工具”参考中描述的过程。

@njit
def get_best_fitting_plane_vector(coords):

    # Calculate the covariance matrix of the coordinates
    covariance_matrix = np.cov(coords,rowvar=False) # Variables = columns

    # Calculate the eigenvalues & eigenvectors of the covariance matrix
    e_val,e_vect = np.linalg.eig(covariance_matrix)

    # The normal vector to the plane is the eigenvector associated to the minimum eigenvalue
    min_eval = np.argmin(e_val)
    normal_vector = e_vect[:,min_eval]

    return normal_vector

在速度方面,重新实现的过程比使用PCA更快,并且如果您使用numba(可以用@njit装饰函数)则可以更快。