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

Question

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

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

1,4-DHP循环http://i.stack.imgur.com/dhNDo.png

1,4-DHP循环，另一个视图http://i.stack.imgur.com/6Xs0z.png

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

如何使和（平方）最小化？我试过最小二乘法，但我也试过。

Answer 1

听起来不错，但是你应该用奇异值分解来代替非线性优化。下面创建转动惯量张量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]

例如：在（10100）处构建一个二维云，该云在x方向上很薄，在y方向上大100倍：

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

拟合平面非常接近于（10100），法线非常接近于x轴。

>>> planeFit(pts)

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

Answer 2

你适合一架飞机这一事实在这里只是有点关系。你要做的是从猜测开始最小化一个特殊的函数。为此，scipy.optimize。注意，不能保证这是全局最优的解决方案，只能保证局部最优的解决方案。一个不同的初始条件可能会收敛到不同的结果，如果你开始接近你正在寻找的局部极小值，这个效果很好。

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

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

Answer 3

最小二乘法应该很容易适应一个平面。一个平面的方程是：ax+by+c=z。所以用所有数据建立这样的矩阵：

    x_0   y_0   1  
A = x_1   y_1   1  
          ... 
    x_n   y_n   1  

以及

    a  
x = b  
    c

以及

    z_0   
B = z_1   
    ...   
    z_n

换句话说：Ax=B，现在求x的系数。但由于你有3个以上的点，系统被过度确定，所以你需要使用左伪逆。所以答案是：

a 
b = (A^T A)^-1 A^T B
c

下面是一些简单的Python代码和一个例子：

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

N_POINTS = 10
TARGET_X_SLOPE = 2
TARGET_y_SLOPE = 3
TARGET_OFFSET  = 5
EXTENTS = 5
NOISE = 5

# create random data
xs = [np.random.uniform(2*EXTENTS)-EXTENTS for i in range(N_POINTS)]
ys = [np.random.uniform(2*EXTENTS)-EXTENTS for i in range(N_POINTS)]
zs = []
for i in range(N_POINTS):
    zs.append(xs[i]*TARGET_X_SLOPE + \
              ys[i]*TARGET_y_SLOPE + \
              TARGET_OFFSET + np.random.normal(scale=NOISE))

# plot raw data
plt.figure()
ax = plt.subplot(111, projection='3d')
ax.scatter(xs, ys, zs, color='b')

# do fit
tmp_A = []
tmp_b = []
for i in range(len(xs)):
    tmp_A.append([xs[i], ys[i], 1])
    tmp_b.append(zs[i])
b = np.matrix(tmp_b).T
A = np.matrix(tmp_A)
fit = (A.T * A).I * A.T * b
errors = b - A * fit
residual = np.linalg.norm(errors)

print "solution:"
print "%f x + %f y + %f = z" % (fit[0], fit[1], fit[2])
print "errors:"
print errors
print "residual:"
print residual

# plot plane
xlim = ax.get_xlim()
ylim = ax.get_ylim()
X,Y = np.meshgrid(np.arange(xlim[0], xlim[1]),
                  np.arange(ylim[0], ylim[1]))
Z = np.zeros(X.shape)
for r in range(X.shape[0]):
    for c in range(X.shape[1]):
        Z[r,c] = fit[0] * X[r,c] + fit[1] * Y[r,c] + fit[2]
ax.plot_wireframe(X,Y,Z, color='k')

ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('z')
plt.show()

你的观点的解决方案是：

0.143509 x + 0.057196 y + 1.129595 = z