平面拟合到 4 个（或更多）XYZ 点

Question

我有 4 个点，它们非常接近一个平面 - 这是 1,4-二氢吡啶循环。

我需要计算从 C3 和 N1 到由 C1-C2-C4-C5 组成的平面的距离。计算距离还可以，但是拟合平面对我来说很难。

1,4-DHP循环：

1,4-DHP循环，另一种观点：

 from array import *
from numpy import *
from scipy import *

# coordinates (XYZ) of C1, C2, C4 and C5
x = [0.274791784, -1.001679346, -1.851320839, 0.365840754]
y = [-1.155674199, -1.215133985, 0.053119249, 1.162878076]
z = [1.216239624, 0.764265677, 0.956099579, 1.198231236]

# plane equation Ax + By + Cz = D
# non-fitted plane
abcd = [0.506645455682, -0.185724560275, -1.43998120646, 1.37626378129]

# creating distance variable
distance =  zeros(4, float)

# calculating distance from point to plane
for i in range(4):
    distance[i] = (x[i]*abcd[0]+y[i]*abcd[1]+z[i]*abcd[2]+abcd[3])/sqrt(abcd[0]**2 + abcd[1]**2 + abcd[2]**2)

print distance

# calculating squares
squares = distance**2

print squares

如何使 sum(squares) 最小化？我试过最小二乘法，但这对我来说太难了。

原文由 XuMuK 发布，翻译遵循 CC BY-SA 4.0 许可协议

Answer 1

你适合飞机的事实在这里只是稍微相关。您要做的是最小化从猜测开始的 特定 函数。为此使用 scipy.optimize 。请注意，不能保证这是 全局最优 解，只能是 局部最优 解。不同的初始条件可能会收敛到不同的结果，如果您从接近您正在寻找的局部最小值开始，这会很有效。

我冒昧地利用 numpy 的广播来清理你的代码：

 import numpy as np

# coordinates (XYZ) of C1, C2, C4 and C5
XYZ = np.array([
        [0.274791784, -1.001679346, -1.851320839, 0.365840754],
        [-1.155674199, -1.215133985, 0.053119249, 1.162878076],
        [1.216239624, 0.764265677, 0.956099579, 1.198231236]])

# Inital guess of the plane
p0 = [0.506645455682, -0.185724560275, -1.43998120646, 1.37626378129]

def f_min(X,p):
    plane_xyz = p[0:3]
    distance = (plane_xyz*X.T).sum(axis=1) + p[3]
    return distance / np.linalg.norm(plane_xyz)

def residuals(params, signal, X):
    return f_min(X, params)

from scipy.optimize import leastsq
sol = leastsq(residuals, p0, args=(None, XYZ))[0]

print("Solution: ", sol)
print("Old Error: ", (f_min(XYZ, p0)**2).sum())
print("New Error: ", (f_min(XYZ, sol)**2).sum())

这给出：

 Solution:  [  14.74286241    5.84070802 -101.4155017   114.6745077 ]
Old Error:  0.441513295404
New Error:  0.0453564286112

原文由 Hooked 发布，翻译遵循 CC BY-SA 4.0 许可协议

Answer 2

这听起来不错，但您应该用 SVD 代替非线性优化。下面创建惯性张量 M，然后使用 SVD 得到平面的法线。这应该是最小二乘拟合的近似值，并且速度更快，更可预测。它返回点云中心和法线。

 def planeFit(points):
    """
    p, n = planeFit(points)

    Given an array, points, of shape (d,...)
    representing points in d-dimensional space,
    fit an d-dimensional plane to the points.
    Return a point, p, on the plane (the point-cloud centroid),
    and the normal, n.
    """
    import numpy as np
    from numpy.linalg import svd
    points = np.reshape(points, (np.shape(points)[0], -1)) # Collapse trialing dimensions
    assert points.shape[0] <= points.shape[1], "There are only {} points in {} dimensions.".format(points.shape[1], points.shape[0])
    ctr = points.mean(axis=1)
    x = points - ctr[:,np.newaxis]
    M = np.dot(x, x.T) # Could also use np.cov(x) here.
    return ctr, svd(M)[0][:,-1]

例如：在 (10, 100) 处构造一个二维云，它在 x 方向上很薄，在 y 方向上大 100 倍：

 >>> pts = np.diag((.1, 10)).dot(randn(2,1000)) + np.reshape((10, 100),(2,-1))

拟合平面非常接近于 (10, 100)，法线非常接近于 x 轴。

 >>> planeFit(pts)

    (array([ 10.00382471,  99.48404676]),
     array([  9.99999881e-01,   4.88824145e-04]))

原文由 Ben 发布，翻译遵循 CC BY-SA 3.0 许可协议