与odrpack平面配合

Question

我正在尝试使用scipy.odr为某些x，y，z点获得最佳拟合平面。

我将平面方程隐式地定义为ax + by + cz + d = 0并且我执行最小二乘（使用scipy.linalg.lstsq）以向odr提供初始估计。

由odr返回的β矢量的组成部分（其中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.

Answer 1

据我了解odr它不是为3D数据制作的，但我可能在这里错了。由于这是一个简单的平面拟合，我建议使用简单的leastsq。此外，请注意，你没有真正有4个自由参数，因为你可以划分a * x + b * y + c * z + d = 0，例如由d提供a' * x + b' * y + c' * z + 1 = 0（如果d不为零）。

相反，如果我们以下面的形式写平面：所有点P为which(P - p0) * n = 0我们已经有免费的odr函数。通过假设平面偏移矢量p0 = s * n是缩放的法向量，可以简化。像这样有3个自由参数，比例s和法向量(theta, phi)的方向角。

相应的拟合如下

import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import numpy as np
from scipy.optimize import leastsq 
from random import random

# for rotating test data
def y_mx( theta ):
    out = np.array( [ np.cos( theta ),0, np.sin( theta ),  0, 1, 0, -np.sin( theta ),0, np.cos( theta ) ] )
    return out.reshape( 3, 3)


# for rotating test data
def z_mx( theta ):
    out = np.array( [ np.cos( theta ), np.sin( theta ), 0, -np.sin( theta ), np.cos( theta ), 0, 0, 0, 1 ] )
    return out.reshape( 3, 3)


# for test data
def make_plane( theta, phi, px, py, pz, n=100 ):
    points=[]
    for i in range( n ):
        x = 1 - 2 * random( )
        y = 1 - 2 * random( )
        z = 0.15 * ( 1 - 2 * random() )
        points += [ np.array( [ x, y, z] ) ]
    points = np.array( points)

    points = [ np.array( [px, py, pz ] ) + np.dot( z_mx( phi ),  np.dot( y_mx( theta ) , p ) ) for p in points  ]
    return np.array( points )


# residual function for leastsq
# note the plane equation is (P - p0) n = 0 if P is member of plane
# and n is normal vector of plane directly provides the normal distance function
# moreover p0 can be chosen to be s * n
def residuals( params, points ):
    scale, theta, phi = params
    nVector = np.array( [ np.sin( theta ) * np.cos( phi ), np.sin( theta ) * np.sin( phi ), np.cos( theta ) ] )
    p0 = scale * nVector
    diff = [ np.dot( p - p0, nVector ) for p in points]
    return diff


# some test data
pnts = make_plane( 1.5, 1.49, .15, .2, .33)

#and the fit
guess=[ 0, 0, 0 ]
bestfit, err = leastsq( residuals, guess, pnts )

#the resulting normal vectot and offset
nVectorFit = np.array( [ np.sin( bestfit[1] ) * np.cos( bestfit[2] ), np.sin( bestfit[1] ) * np.sin( bestfit[2] ), np.cos( bestfit[1] ) ] )
p0Fit = bestfit[0] * nVectorFit

# converting to standard plane equation
a = nVectorFit[0] / nVectorFit[1]
c = nVectorFit[2] / nVectorFit[1]
d = bestfit[0] / nVectorFit[1]

# plane equation data
X = np.linspace( -.6, .6, 20 )
Z = np.linspace( -.6, .6, 20 )
XX, ZZ = np.meshgrid( X, Z )
YY = -a * XX - c * ZZ + d

#plotting
fig = plt.figure()
ax = fig.add_subplot( 1, 1, 1, projection='3d')
# original data
ax.scatter( pnts[:,0], pnts[:,1] , pnts[:,2])
# offset vector
ax.plot( [0, p0Fit[0] ], [0, p0Fit[1] ], [0, p0Fit[2] ], color = 'r')
# fitted plane
ax.plot_wireframe(XX, YY, ZZ , color = '#9900bb')

ax.set_xlim( [-1,1] )
ax.set_ylim( [-1,1] )
ax.set_zlim( [-1,1] )
ax.set_xlabel("x")
ax.set_ylabel("y")
ax.set_zlabel("z")
plt.show()

提供

蓝点是噪声数据，紫色是拟合平面，红色是偏移矢量。

很容易看出，这里的情况y = a * x + c * z + d，a, c, d是从拟合结果直接计算出来的。

Answer 2

对于多维数据，ODR完全没问题，你正朝着正确的方向前进。

您只是选择了不好用来使用隐式版本的ODR和f_3定义。问题是你有一个函数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拟合（预期）略有不同。