Plan de montage à 4 (ou plus) points XYZ

Le fait que vous êtes à côté d'un avion est peu pertinent ici. Ce que vous essayez de faire est de minimiser une particulier fonction à partir d'une supposition. Pour que l'utilisation scipy.optimize. Notez qu'il n'y a aucune garantie que cela est le globalement optimale solution, seulement localement optimale. Une autre condition initiale peut converger vers un résultat différent, cela fonctionne bien si vous commencez à fermer les minima locaux, vous êtes à la recherche d'.

J'ai pris la liberté de nettoyer votre code en profitant de numpy de radiodiffusion:

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())

Cela donne:

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

Des moindres carrés doivent s'adapter à un avion facilement. L'équation d'un plan est: ax + by + c = z. Donc, mettre en place des matrices comme ça avec toutes vos données:

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

Et

    a  
x = b  
    c

Et

    z_0   
B = z_1   
    ...   
    z_n

En d'autres termes: Ax = B. Maintenant résoudre pour x qui sont vos coefficients. Mais puisque vous avez plus de 3 points, le système est sur-déterminé de sorte que vous devez utiliser la gauche pseudo-inverse. Donc la réponse est:

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

Et ici est un simple code Python avec un exemple:

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()

La solution pour vos points:

0.143509 x + 0.057196 y + 1.129595 = z

Un autre moyen, hormis la svd d'atteindre rapidement une solution tout en traitant avec des valeurs aberrantes ( lorsque vous avez un grand ensemble de données ) est ransac :

     def fit_plane(voxels, iterations=50, inlier_thresh=10):  # voxels : x,y,z
inliers, planes = [], []
xy1 = np.concatenate([voxels[:, :-1], np.ones((voxels.shape[0], 1))], axis=1)
z = voxels[:, -1].reshape(-1, 1)
for _ in range(iterations):
random_pts = voxels[np.random.choice(voxels.shape[0], voxels.shape[1] * 10, replace=False), :]
plane_transformation, residual = fit_pts_to_plane(random_pts)
inliers.append(((z - np.matmul(xy1, plane_transformation)) <= inlier_thresh).sum())
planes.append(plane_transformation)
return planes[np.array(inliers).argmax()]
def fit_pts_to_plane(voxels):  # x y z  (m x 3)
# https: //math.stackexchange.com /questions /99299 /best - fitting - plane - given - a - set - of - points
xy1 = np.concatenate([voxels[:, :-1], np.ones((voxels.shape[0], 1))], axis=1)
z = voxels[:, -1].reshape(-1, 1)
fit = np.matmul(np.matmul(np.linalg.inv(np.matmul(xy1.T, xy1)), xy1.T), z)
errors = z - np.matmul(xy1, fit)
residual = np.linalg.norm(errors)
return fit, residual