python - 找到两条相交线之间的等分平面
问题描述
给定在 3D 空间中的一点相交的两个向量,我想确定(并显示)平分这些线的平面。
到目前为止,我已经完成了以下工作:
0)绘制向量(红色)
1)计算两个向量之间的叉积(蓝色)和点积角度。
2) 使用叉积向量作为旋转轴,将向量之一旋转它们之间的 1/2 角度。结果垂直于平面(绿色)。
3) 使用法线,在平面的 a x+b y+c*z+d = 0 方程中求解“d”。
4) 绘制合成平面。
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import numpy as np
import math
def unit_vector(vector):
""" Returns the unit vector of the vector. """
return vector / np.linalg.norm(vector)
def angle_between(v1, v2):
v1_u = unit_vector(v1)
v2_u = unit_vector(v2)
return np.arccos(np.clip(np.dot(v1_u, v2_u), -1.0, 1.0))
def rotation_matrix(axis, theta):
axis = np.asarray(axis)
axis = axis / math.sqrt(np.dot(axis, axis))
a = math.cos(theta / 2.0)
b, c, d = -axis * math.sin(theta / 2.0)
aa, bb, cc, dd = a * a, b * b, c * c, d * d
bc, ad, ac, ab, bd, cd = b * c, a * d, a * c, a * b, b * d, c * d
return np.array([[aa + bb - cc - dd, 2 * (bc + ad), 2 * (bd - ac)],
[2 * (bc - ad), aa + cc - bb - dd, 2 * (cd + ab)],
[2 * (bd + ac), 2 * (cd - ab), aa + dd - bb - cc]])
#Path of points to plot
pathPoints = np.array ([[ 3, 0, 30],
[13, 7, 25],
[23, 12, 12]])
#Initialize flight path array
flightPath = []
lineSegment = np.array ([0,0,0,0,0,0])
for cntr, point in enumerate (pathPoints):
#skip the last point in the sequence
if pathPoints.shape[0] == cntr + 1: #shape is 1-based
break
print ('Point #%d' % cntr, point)
nextPoint = pathPoints[cntr+1]
print ('Next Point #%d' % cntr, nextPoint)
lineSegment = np.zeros(6)
lineSegment[:3] = point
diff = np.subtract (nextPoint, point)
print ('Diff #%d' % cntr, diff)
lineSegment[3:] = diff
print ('LineSeg #%d' % cntr, lineSegment)
flightPath.append (lineSegment)
print ('Flight Path array ', flightPath)
#Create the plot
plotSize = 15
plt3d = plt.figure(figsize=(plotSize,plotSize)).gca(projection='3d')
#set the plot limits
axisMin=0
axisMax=30
plt3d.set_xlim([axisMin,axisMax])
plt3d.set_ylim([axisMin,axisMax])
plt3d.set_zlim([axisMin,axisMax])
plt3d.set_xlabel('x')
plt3d.set_ylabel('y')
plt3d.set_zlabel('z')
for vector in flightPath:
v = np.array([vector[3],vector[4],vector[5]])
vlength = np.linalg.norm(v)
plt3d.quiver(vector[0],vector[1],vector[2],vector[3],vector[4],vector[5],
pivot='tail', arrow_length_ratio=2.0/vlength,color='red')
#Compute the cross product at the intersection points, and then a plane which
#bisects the intersection point.
flightPlanes = []
flightPlaneOffsets = []
for cntr, currVec in enumerate (flightPath):
#skip the last point in the sequence
if len (flightPath) == cntr + 1: #shape is 1-based
break
print ('Vector #%d' % cntr, currVec)
#Compute the cross product normal
nextVec = flightPath[cntr+1]
print ('Next Vector #%d' % cntr, nextVec)
#Compute the cross product of the differences
crossProd = np.cross (np.negative (currVec[3:]), nextVec[3:])
vlength = np.linalg.norm(crossProd)
print ('Cross Product #%d' % cntr, crossProd)
#Scale the length of the cross product in order to display on the
#plot. Determining the plane doesn't depend on length.
crossProd = np.multiply (crossProd, 15/vlength)
print ('Scaled Cross Product #%d' % cntr, crossProd)
#Recompute the scaled length
vlength = np.linalg.norm(crossProd)
plt3d.quiver(nextVec[0], nextVec[1], nextVec[2],
crossProd[0], crossProd[1], crossProd[2],
pivot='tail', arrow_length_ratio=2.0/vlength,color='blue')
#Generate the bisecting plane
#Compute the angle between the vectors
bisectAngle = angle_between (np.negative (currVec[3:]), nextVec[3:]) / 2
print ('Bisect angle between #%d %.2f deg' % (cntr, np.degrees(bisectAngle)))
#Determining normal vector
#Compute the normal vector to the plane
#See https://stackoverflow.com/questions/6802577/rotation-of-3d-vector
normalVecA = np.dot(rotation_matrix(crossProd, 2*np.pi-bisectAngle), nextVec[:3])
#Scale the length of the normal vector in order to display on the
#plot. Determining the plane doesn't depend on length.
vlength = np.linalg.norm(normalVecA)
plt3d.quiver(nextVec[0], nextVec[1], nextVec[2],
normalVecA[0], normalVecA[1], normalVecA[2],
pivot='tail', arrow_length_ratio=2.0/vlength,color='green')
print ('Normal vector A #%d' % cntr, normalVecA)
#Create a view of one of the normal vectors
normalVec = normalVecA
#Plane point is the intersection point
planePoint = np.array(nextVec[:3])
print ('Plane point #%d' % cntr, planePoint)
# A plane is a*x+b*y+c*z+d=0
# [a,b,c] is the normal. Thus, we have to calculate
# d and we're set
d = -planePoint.dot(normalVec)
print ('D offset is #%d' % cntr, d)
# create x,y
gridSize = 3
ptsSpacing = 10 / abs (normalVec[2])
xRange = np.array(np.arange(nextVec[0]-gridSize,
nextVec[0]+gridSize+1, ptsSpacing))
yRange = np.array(np.arange(nextVec[1]-gridSize,
nextVec[1]+gridSize+1, ptsSpacing))
print ('Xrange #%d' % cntr, xRange)
xx, yy = np.meshgrid(xRange, yRange)
# calculate corresponding z
z = (-normalVec[0] * xx - normalVec[1] * yy - d) * 1. /normalVec[2]
flightPlanes.append(normalVec)
flightPlaneOffsets.append(d)
# plot the surface
plt3d.plot_surface(xx, yy, z, alpha=0.3)
plt.show()
我希望作为旋转轴的叉积向量(蓝色)也应该在二等分平面上(但事实并非如此)。以下代码行是关键:
normalVecA = np.dot(rotation_matrix(crossProd, 2*np.pi-bisectAngle), nextVec[:3])
我尝试将角度设置为“bisectAngle”和“-bisectAngle”,但似乎无法解决问题。
解决方案
在这种特殊情况下,新平面的法线向量只是两个向量的串联。我会这样:
- 给定三点
p1,p2,p3
计算a1 = unit(p2-p1)
,a2=unit(p3-p2)
- 计算平面的法向量为
n = (a2+a1)
p2
构造通过normal的平面n
。
推荐阅读
- java - 如果 editText 包含文本,则发送文本或使用 ontouchlistener 录制音频
- karate - 无法将两个参数作为参数传递给 Javascript 函数
- java - 如何使用从数据库中获取的路径作为字符串播放原始文件夹中的歌曲?
- aem - AEM 中的 YUI 压缩器从颜色值中删除 % 符号
- html - 如何在html中将列向右对齐
- java - 从 Java 应用程序中而不是从命令行启动时系统命令挂起
- c# - 您必须先加载或分配此属性,然后才能读取它
- puppet - Puppet:我尝试在 mainifest 中获取内存事实并得到错误 Operator '[]' is not applicable to an Undef Value
- python - 无法从“sklearn.preprocessing”导入名称“StandardScalar”
- python - Python条带函数