У меня есть набор трехмерных точек, в которые я хотел бы подогнать плоскость, используя уравнение плоскости: ax + by + c = z. Поскольку система перегружена, я использовал матрицы как Ax = b, где A = [xi, yi, 1], x = [a, b, c], B = [zi]. Это код, который я использовал.
Проблема, с которой я столкнулся, заключается в том, что плоскость не строится с использованием plot_wireframe. Кто-нибудь может помочь?
Спасибо!
import numpy as np
import math as m
import random
import scipy.optimize as optimization
import matplotlib.pyplot as plt
#the cartesian cordinates of targets in the CCD frame is known
CCD_target_1 = np.matrix([[0.034125],[0.039125],[0],[1]])
CCD_target_2 = np.matrix([[-0.034125],[0.039125],[0],[1]])
CCD_target_3 = np.matrix([[-0.034125],[-0.039125],[0],[1]])
CCD_target_4 = np.matrix([[0.034125],[-0.039125],[0],[1]])
#transformation matrix for converting the CCD frame to the LT frame
CCD_to_LT = np.matrix([[1,0,0,1.2],[0,1,0,0],[0,0,1,11],[0,0,0,1]])
#converting the target points from CCD to LT frame
LT_target_1 = CCD_to_LT*CCD_target_1
LT_target_2 = CCD_to_LT*CCD_target_2
LT_target_3 = CCD_to_LT*CCD_target_3
LT_target_4 = CCD_to_LT*CCD_target_4
#fuction for converting cartesian to spherical co ordinate system
def cart_to_sphere(A):
rho = m.sqrt(A[0,0]**2 + A[1,0]**2 + A[2,0]**2)
psi = m.acos(A[1,0]/A[0,0])
theta = m.acos(A[2,0]/rho)
return np.matrix([[rho],[theta],[psi],[1]])
#function for converting spherical to cartesian co ordinate system
def sphere_to_cart(A):
y = A[0,0]*m.sin(A[1,0])*m.cos(A[2,0])
x = A[0,0]*m.sin(A[1,0])*m.sin(A[2,0])
z = A[0,0]*m.cos(A[1,0])
return np.matrix([[x],[y],[z],[1]])
#function to find the centre through mid point of diagonal
def centre(A,B):
x = ((A[0,0]+B[0,0])/2)
y = ((A[1,0]+B[1,0])/2)
return (x,y)
#The equation for a plane is: ax+by+c=z
xc = []
yc = []
zc = []
for i in range(100):
#converting the LT targets from cartesian coordinates to spherical coordinates and introducing errors
spherical_target_1 = cart_to_sphere(LT_target_1) + np.matrix([[random.uniform(-81e-6,81e-6)],[random.uniform(-7.2722E-6,0.00000727)],[random.uniform(-7.2722E-6,0.00000727)],[0]])
spherical_target_2 = cart_to_sphere(LT_target_2) + np.matrix([[random.uniform(-81e-6,81e-6)],[random.uniform(-7.2722E-6,0.00000727)],[random.uniform(-7.2722E-6,0.00000727)],[0]])
spherical_target_3 = cart_to_sphere(LT_target_3) + np.matrix([[random.uniform(-81e-6,81e-6)],[random.uniform(-7.2722E-6,0.00000727)],[random.uniform(-7.2722E-6,0.00000727)],[0]])
spherical_target_4 = cart_to_sphere(LT_target_4) + np.matrix([[random.uniform(-81e-6,81e-6)],[random.uniform(-7.2722E-6,0.00000727)],[random.uniform(-7.2722E-6,0.00000727)],[0]])
#converting the LT targets back to cartesian co ordinate system
cartesian_target_1 = sphere_to_cart(spherical_target_1)
cartesian_target_2 = sphere_to_cart(spherical_target_2)
cartesian_target_3 = sphere_to_cart(spherical_target_3)
cartesian_target_4 = sphere_to_cart(spherical_target_4)
#making matrices xc, yc, zc to represent combined x,y,z coordinates of the targets in the cartesian LT frame
xc.append(cartesian_target_1[0,0])
xc.append(cartesian_target_2[0,0])
xc.append(cartesian_target_3[0,0])
xc.append(cartesian_target_4[0,0])
yc.append(cartesian_target_1[1,0])
yc.append(cartesian_target_2[1,0])
yc.append(cartesian_target_3[1,0])
yc.append(cartesian_target_4[1,0])
zc.append(cartesian_target_1[2,0])
zc.append(cartesian_target_2[2,0])
zc.append(cartesian_target_3[2,0])
zc.append(cartesian_target_4[2,0])
# plot raw data
plt.figure()
ax = plt.subplot(111, projection='3d')
ax.scatter(xc, yc, zc, color='b')
# do a plane fit for the set of 3D points
tmp_A = []
tmp_b = []
for i in range(len(xc)):
tmp_A.append([xc[i], yc[i], 1])
tmp_b.append(zc[i])
b = np.matrix(tmp_b).T
A = np.matrix(tmp_A)
fit = (A.T * A).I * A.T * b
# plot plane for the 3D points
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()