我正在尝试使用scipy.odr来获得一些x，y，z点的最佳拟合平面。
我将平面方程式隐式定义为ax + by + cz + d = 0，然后执行最小二乘（使用scipy.linalg.lstsq）为odr提供初始估计。
odr返回的beta向量的分量（其中beta = [a，b，c，d]）的大小在1e167和1e172之间...这样的结果值得信赖吗？我觉得数字很荒谬。
请注意，这些点来自3D扫描相对平坦的平面，该平面几乎平行于xz平面（几乎垂直）。
这是odr结果的pprint（）：
'
Beta: [ 3.14570111e-170 3.21821458e-169 4.49232028e-172 4.49374557e-167]
Beta Std Error: [ 0. 0. 0. 0.]
Beta Covariance: [[ 6.37459471e-10 -8.57690019e-09 -2.18092934e-11 -1.13009384e-06]
[ -8.57690019e-09 5.11732570e-07 1.30123070e-09 6.74263262e-05]
[ -2.18092934e-11 1.30123070e-09 5.22674068e-12 1.70799469e-07]
[ -1.13009384e-06 6.74263262e-05 1.70799469e-07 8.88444676e-03]]
Residual Variance: 0.0
Inverse Condition #: 0.0010484041422201213
Reason(s) for Halting:
Sum of squares convergence
None
'
我正在使用的代码：
import numpy as np
import scipy.linalg
from scipy import odr
import pickle
def planar_fit(points):
# best-fit linear plane
a = np.c_[points[:, 0], points[:, 1], np.ones(points.shape[0])]
c, _, _, _ = scipy.linalg.lstsq(a, points[:, 2]) # coefficients
# The coefficients are returned as an array beta=[a, b, c, d] from the implicit form 'a*x + b*y + c*z + d = 0'.
beta = np.r_[c[0], c[1], -1, c[2]] / c[2]
return beta
def odr_planar_fit(points):
def f_3(beta, xyz):
""" implicit definition of the plane"""
return beta[0] * xyz[0] + beta[1] * xyz[1] + beta[2] * xyz[2] + beta[3]
# # Coordinates of the 2D points
x = points[:, 0]
y = points[:, 1]
z = points[:, 2]
# Use least squares for initial estimate.
beta0 = planar_fit(points)
# Create the data object for the odr. The equation is given in the implicit form 'a*x + b*y + c*z + d = 0' and
# beta=[a, b, c, d] (beta is the vector to be fitted). The positional argument y=1 means that the dimensionality
# of the fitting is 1.
lsc_data = odr.Data(np.row_stack([x, y, z]), y=1)
# Create the odr model
lsc_model = odr.Model(f_3, implicit=True)
# Create the odr object based on the data, the model and the first estimation vector.
lsc_odr = odr.ODR(lsc_data, lsc_model, beta0)
# run the regression.
lsc_out = lsc_odr.run()
return lsc_out, beta0
def main():
#import from pickle.
with open('./points.pkl', 'rb') as f:
points = np.array(pickle.load(f))
# Perform the ODR
odr_out, lstsq = odr_planar_fit(points)
print(lstsq)
print(odr_out.pprint())
main()
The pickle containing my points.
最佳答案
ODR对多维数据完全适用，您的方向正确。
您只是选择了不好，无法在f_3定义中使用ODR的隐式版本。问题是您有一个函数A*X=0，尝试将其最小化而没有任何其他限制。当然，优化器可以做的最好的事情就是将A的幅度最小化为零-从而最大程度地减少误差！为了使隐式优化起作用，您需要以某种方式引入对A大小的约束，例如除以最后一个数字：
def f_3(beta, xyz):
""" implicit definition of the plane"""
return beta[0]/beta[3] * xyz[0] + beta[1]/beta[3] * xyz[1] + beta[2]/beta[3] * xyz[2] + 1.0
这样，优化器除了执行您想要的操作之外别无选择：)
另外，您也可以将模型转换为显式形式
y = ax + cz + d，该形式不会受到幅度问题的困扰（一直以来都是b == 1）。当然，通过将点移至原点并将其缩放为距离的单位差异，可以获得更高的精度。
由于我也将要使用ODR，因此我对它的属性感到好奇，因此我四处探索以发现它的精确性和敏感性，结果如下：https://gist.github.com/peci1/fb1cea77c41fe8ace6c0db8ef82539a3。
我测试了隐式和显式ODR，是否进行了归一化，并使用LSQ或较差的值进行了初步猜测（以查看对猜测的敏感程度）。您的数据看起来像这样：
基本上，黄色和灰色平面是未经归一化的隐式拟合，结果非常糟糕，而其他ODR拟合大致相同。您会看到ODR拟合与（淡蓝色）LSQ拟合（预期的）有些不同。
关于python - Plannar适合odrpack，我们在Stack Overflow上找到一个类似的问题：https://stackoverflow.com/questions/48115464/