问题描述
我有一个使用函数 z = f(x,y) 生成的曲面,其中 x,y 数组传递给 numpy meshgrid 以生成曲面点。从任意两个点,我使用线方程和平面方程来找到交点。但是,我的脚本未能找到交集,特别是在我求解 t 的步骤中。
3d 线方程:t = x-xo/l = y-yo/m =z-zo/n 平面方程:ax+by+cz=d
我的测试脚本在下面,仅生成 NaN 或 Inf。我期望的输出是一个 t 值,我可以将它插回到线方程中以找到我的交点 (x,y,z)
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
# generate the surface
x = np.linspace(0,100,10)
y = np.linspace(0,10)
xv,yv = np.meshgrid(x,y)
zv = (xv**2+yv)
pos = np.vstack([xv.ravel(),yv.ravel(),zv.ravel()]).T
fig = plt.figure()
ax = fig.gca(projection='3d')
surf = ax.plot_surface(xv,yv,zv)
plt.show()
zv_n1 = []
zv_n2 = []
p3 = []
i_range = zv.shape[0]
j_range = zv.shape[1]
# find AB,AC vectors,a surface point to later compute the normal vec and then plane equation
for i in range(i_range):
for j in range(j_range):
if i<i_range-1 and j<j_range-1:
zv_a = zv[i,j+1]-zv[i,j]
zv_b = zv[i+1,j]-zv[i,j]
p3_arr = np.array([xv[i+1,j],yv[i+1,zv[i+1,j]])
if i == i_range-1 and j < j_range-1:
zv_a = zv[i,j]
zv_b = zv[i-1,j]
p3_arr = np.array([xv[i-1,yv[i-1,zv[i-1,j]])
if i < i_range-1 and j == j_range-1:
zv_a = zv[i,j-1]-zv[i,j]])
zv_n1.append(zv_a)
zv_n2.append(zv_b)
p3.append(p3_arr)
# sorting in pandas df
temp_df = pd.DataFrame(pd.Series(pos[:,0]),columns = ['x_val'])
temp_df['y_val'] = pd.DataFrame(pd.Series(pos[:,1]))
temp_df['z_val'] = pd.DataFrame(pd.Series(pos[:,2]))
temp_df['x1_del'] = x[1] # delta x is same for panels
temp_df['y1_del'] = y[1] # delta y is same for panels
temp_df['z1_del'] = pd.DataFrame(pd.Series(zv_n1))
temp_df['x2_del'] = x[0] # delta x is same for panels
temp_df['y2_del'] = y[0] # delta y is same for panels
temp_df['z2_del'] = pd.DataFrame(pd.Series(zv_n2))
temp_df['p3'] = pd.DataFrame(pd.Series(p3))
# compute normal vec (and thus find a,b,c coeff,and then find d value)
cross_list = []
d_val_list = []
for i in range(len(temp_df)):
a_vec = np.array([temp_df.iloc[i]['x1_del'],temp_df.iloc[i]['y1_del'],temp_df.iloc[i]['z1_del']])
b_vec = np.array([temp_df.iloc[i]['x2_del'],temp_df.iloc[i]['y2_del'],temp_df.iloc[i]['z2_del']])
nor_vec = np.cross(a_vec,b_vec)
cross_list.append(nor_vec)
d_val = np.dot(nor_vec,temp_df.iloc[i]['p3'])
d_val_list.append(d_val)
temp_df['normal_vec'] = pd.DataFrame(pd.Series(cross_list))
temp_df['plane_d_val'] = pd.DataFrame(pd.Series(d_val_list))
# create two points and find the line vector
p1 = np.array([0,60,0])
p2 = np.array([0,5000])
p_vec = p2-p1
# solve for t - rearrange line equations (t = X-x/l = Y-y/m = Z-z/n) to (X =tl+x) and so on..
# and plug in plane equation ax+by+cz=d
for i in range(len(temp_df)):
t = (temp_df.iloc[i]['plane_d_val'] - \
temp_df.iloc[i]['normal_vec'][0] * temp_df.iloc[i]['p3'][0] - \
temp_df.iloc[i]['normal_vec'][1] * temp_df.iloc[i]['p3'][1] - \
temp_df.iloc[i]['normal_vec'][2] * temp_df.iloc[i]['p3'][2])/(\
temp_df.iloc[i]['normal_vec'][0] * p_vec[0] +\
temp_df.iloc[i]['normal_vec'][1] * p_vec[1] +\
temp_df.iloc[i]['normal_vec'][2] * p_vec[2])
print(t)
# after finding t,plug in line equation to solve for X
我希望得到有关此问题的指导。也欢迎任何关于替代方法/库的建议,因为唯一的目的是确定交叉点。
谢谢,
解决方法
所以我想出了我的代码问题。使用平面上的点而不是直线上的点。
此外,此测试脚本对曲面函数也没有用。只用数学方法求解方程似乎是明智的。
无论如何更新的代码如下:
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
x = np.linspace(0,100,10)
y = np.linspace(0,10)
xv,yv = np.meshgrid(x,y)
zv = (xv**2+yv)
fig = plt.figure()
ax = fig.gca(projection='3d')
surf = ax.plot_surface(xv,yv,zv)
plt.show()
xv_n1 = []
xv_n2 = []
yv_n1 = []
yv_n2 = []
zv_n1 = []
zv_n2 = []
p3 = []
i_range = zv.shape[0]
j_range = zv.shape[1]
count = -1
for i in range(i_range):
for j in range(j_range):
count += 1
if i<i_range-1 and j<j_range-1:
# now x
xv_a = xv[i,j+1]-xv[i,j]
xv_b = xv[i+1,j]-xv[i,j]
# now y
yv_a = yv[i,j+1]-yv[i,j]
yv_b = yv[i+1,j]-yv[i,j]
#now z
zv_a = zv[i,j+1]-zv[i,j]
zv_b = zv[i+1,j]-zv[i,j]
p3_arr = np.array([xv[i+1,j],yv[i+1,zv[i+1,j]])
if i == i_range-1 and j < j_range-1:
# now x
xv_a = xv[i,j]
xv_b = xv[i-1,j]
yv_b = yv[i-1,j]
zv_b = zv[i-1,j]
p3_arr = np.array([xv[i-1,yv[i-1,zv[i-1,j]])
if i < i_range-1 and j == j_range-1:
# now x
xv_a = xv[i,j-1]-xv[i,j-1]-yv[i,j-1]-zv[i,j]])
xv_n1.append(xv_a)
xv_n2.append(xv_b)
yv_n1.append(yv_a)
yv_n2.append(yv_b)
zv_n1.append(zv_a)
zv_n2.append(zv_b)
p3.append(p3_arr)
pos = np.vstack([xv.ravel(),yv.ravel(),zv.ravel()]).T
temp_df = pd.DataFrame(pd.Series(pos[:,0]),columns = ['x_val'])
temp_df['y_val'] = pd.DataFrame(pd.Series(pos[:,1]))
temp_df['z_val'] = pd.DataFrame(pd.Series(pos[:,2]))
temp_df['x1_del'] = x[1]
temp_df['y1_del'] = y[1]
temp_df['z1_del'] = pd.DataFrame(pd.Series(zv_n1))
temp_df['x2_del'] = x[0]
temp_df['y2_del'] = y[0]
temp_df['z2_del'] = pd.DataFrame(pd.Series(zv_n2))
temp_df['p3'] = pd.DataFrame(pd.Series(p3))
cross_list = []
d_val_list = []
for i in range(len(temp_df)):
a_vec = np.array([temp_df.iloc[i]['x1_del'],temp_df.iloc[i]['y1_del'],temp_df.iloc[i]['z1_del']])
b_vec = np.array([temp_df.iloc[i]['x2_del'],temp_df.iloc[i]['y2_del'],temp_df.iloc[i]['z2_del']])
nor_vec = np.cross(a_vec,b_vec)
cross_list.append(nor_vec)
d_val = np.dot(nor_vec,temp_df.iloc[i]['p3'])
d_val_list.append(d_val)
temp_df['normal_vec'] = pd.DataFrame(pd.Series(cross_list))
temp_df['plane_d_val'] = pd.DataFrame(pd.Series(d_val_list))
p1 = np.array([0,60,6000])
p2 = np.array([100,6000])
p_vec = p2-p1
for i in range(len(temp_df)):
t = (temp_df.iloc[i]['plane_d_val'] - \
temp_df.iloc[i]['normal_vec'][0] * p1[0] - \
temp_df.iloc[i]['normal_vec'][1] * p1[1] - \
temp_df.iloc[i]['normal_vec'][2] * p1[2])/(\
temp_df.iloc[i]['normal_vec'][0] * p_vec[0] +\
temp_df.iloc[i]['normal_vec'][1] * p_vec[1] +\
temp_df.iloc[i]['normal_vec'][2] * p_vec[2])
p1 = np.array([0,6000])
p_vec = p2-p1
for i in range(len(temp_df)):
t = (temp_df.iloc[i]['plane_d_val'] - \
temp_df.iloc[i]['normal_vec'][0] * p2[0] - \
temp_df.iloc[i]['normal_vec'][1] * p2[1] - \
temp_df.iloc[i]['normal_vec'][2] * p2[2])/(\
temp_df.iloc[i]['normal_vec'][0] * p_vec[0] +\
temp_df.iloc[i]['normal_vec'][1] * p_vec[1] +\
temp_df.iloc[i]['normal_vec'][2] * p_vec[2])
if np.isnan(t) == True:
q = np.nan
elif t == 0:
q = np.nan
elif np.isinf(t) == True:
q = np.nan
else:
x = t*p_vec[0]+p2[0]
y = t*p_vec[1]+p2[1]
z = t*p_vec[2]+p2[2]
q = np.array([x,y,z])
print(q)