0

I have 3 points on a plane, 3 points on another plane. I want the angle between the planes so I use the points to make vectors, cross product to find the normal, dot product and divide by magnitude and inverse cos to find me the angle.

I was expecting 45 or 0.785398 as the angle between the normal of the planes from the points. I get 0.6154797086703875

There must be something wrong with my logic, can anyone help please? Thank you.

 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")
1
  • You are normalizing by the wrong vectors. You should use the lengths of the normals (which you call plane1/2) instead. Dec 20, 2018 at 12:57

2 Answers 2

2

Use np.linalg.norm:

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

As mentioned in the comments the error is that you are using the wrong vectors for normalization.

I'd like to add that Binet-Cauchy works here as well and is way faster than using cross product.

import math

def angle_BC(plane1, plane2):
    planes = np.concatenate((plane1[:2]-plane1[2], plane2[:2]-plane2[2]), axis=0)
    a = ([email protected]).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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.