Python Numpy并没有给出平面之间角度的预期结果

Question

我在一个平面上有3个点，在另一个平面中有3个。我想要平面之间的角度，所以我用点来做向量，叉积来求法线，点积，除以幅度和反余弦来求角度。

我预计45或0.785398是这些点的平面法线之间的角度。我得到0.6154797086703875

我的逻辑一定有问题，有人能帮忙吗？非常感谢。

import numpy as np
"""
y axis up page
x axis to the right on page
z axis towards user out of page
point1 -1,1,0
point2 1,1,0
point3 1,1,1
so horizontal plane at y=1
point4 0,0,0
point5 1,1,0
point6 1,1,1
plane at 45 degrees meeting other plane at 1,1,n
"""
#use points to make vectors
#point2 - point1
vx1 = 2#1--1
vy1 = 0#1-1
vz1 = 0#0-0
#point3 - point1
vx2 = 2#1--1
vy2 = 0#1-1
vz2 = 1#1-0
#cross product to find normal
plane1 = np.cross([vx1,vy1,vz1],[vx2,vy2,vz2])
#use points to make vectors
#point5 - point4
vx1 = 1#1-0
vy1 = 1#1-0
vz1 = 0#0-0
#point6 - point4
vx2 = 1#1-0
vy2 = 1#1-0
vz2 = 1#1-0
#cross product to find normal
plane2 = np.cross([vx1,vy1,vz1],[vx2,vy2,vz2])
#angle between the two normals
#dot product
ang = np.dot(plane1,plane2)
#divide by magnitude of vectors
ang = ang / (np.sqrt((vx1*vx1)+(vy1*vy1)+(vz1*vz1))*np.sqrt((vx2*vx2)+(vy2*vy2)+(vz2*vz2)))
#inverse cos to find angle
ang = np.arccos(ang)
#should be 45 or 0.785398
print("Angle calculated",ang)
if (ang < 95) and (ang > 85):
print("RIGHT ANGLE")
else:
print("OTHER ANGLE")

Answer 1

使用np.linalg.norm:

ang = np.arccos(np.dot(plane1, plane2) / (np.linalg.norm(plane1) * np.linalg.norm(plane2)))
ang
>> 0.78539816

Answer 2

如注释中所述，错误在于使用了错误的矢量进行规范化。

我想补充一点，Binet Cauchy在这里也能工作，而且比使用叉积快得多。

import math
def angle_BC(plane1, plane2):
planes = np.concatenate((plane1[:2]-plane1[2], plane2[:2]-plane2[2]), axis=0)
a = (planes@planes.T).ravel().tolist()
return math.acos((a[2]*a[7]-a[3]*a[6]) / math.sqrt((a[0]*a[5]-a[1]*a[1])*(a[10]*a[15]-a[11]*a[11])))... 
def angle_CP(plane1, plane2):
n1 = np.cross(*plane1[:2]-plane1[2])
n2 = np.cross(*plane2[:2]-plane2[2])
return math.acos(n1@n2 / math.sqrt(n1@n1 * n2@n2))
plane1, plane2 = np.array(((-1,1,0),(1,1,0),(1,1,1))), np.array(((0,0,0),(1,1,0),(1,1,1)))
angle_BC(plane1, plane2), angle_CP(plane1, plane2), np.pi/4
# (0.7853981633974484, 0.7853981633974484, 0.7853981633974483)
from timeit import repeat
min(repeat(lambda: angle_BC(plane1, plane2), number=10000))
# 0.08035301300697029
min(repeat(lambda: angle_CP(plane1, plane2), number=10000))
# 0.9213669009623118