python code to calculate angle between three point using their 3D coordinates

Question

I have write down a code to calculate angle between three points using their 3D coordinates.

import  numpy as np

a = np.array([32.49, -39.96,-3.86])

b = np.array([31.39, -39.28, -4.66])

c = np.array([31.14, -38.09,-4.49])

f = a-b # normalization of vectors
e = b-c # normalization of vectors

angle = dot(f, e) # calculates dot product 
print degrees(cos(angle))  # calculated angle in radians to degree 

output of the code:

degree 33.4118214995

but when i used one of the software to calculate the same it gives output bit different 120 degree. please help

reference i have used to write the program:

(How to calculate bond angle in protein db file?)

Solution 1

Your original code is pretty close. Adomas.m's answer is not very idiomatic numpy:

import numpy as np

a = np.array([32.49, -39.96,-3.86])
b = np.array([31.39, -39.28, -4.66])
c = np.array([31.14, -38.09,-4.49])

ba = a - b
bc = c - b

cosine_angle = np.dot(ba, bc) / (np.linalg.norm(ba) * np.linalg.norm(bc))
angle = np.arccos(cosine_angle)

print np.degrees(angle)

Solution 2

I guess numpy is quite enough:

    from numpy import *
    from numpy.linalg import norm
    a = array([32.49, -39.96,-3.86])
    b = array([31.39, -39.28, -4.66])
    c = array([31.14, -38.09,-4.49])
    f = b-a 
    e = b-c 
    abVec = norm(f)
    bcVec = norm(e)
    abNorm = f / abVec;
    bcNorm = e / bcVec;
    res = abNorm[0] * bcNorm[0] + abNorm[1] * bcNorm[1] + abNorm[2] * bcNorm[2];
    angle = arccos(res)*180.0/ pi
    print angle

also the res can be calculated with dot:

    res = abNorm[0] * bcNorm[0] + abNorm[1] * bcNorm[1] + abNorm[2] * bcNorm[2];
    res = dot(abNorm, bcNorm)

Comments

[MSeifert]

1.) You need to subtract the same point to get the vectors you are looking for (see answer to your other question). 2.) You have to normalize the vectors (that's something else than subtraction!) 3.) what other software did you use? 4.) There are python scripts avaiable at google where you can compare your solution to.

[MB-F]

5.) You need to use the inverse of the cosine function (acos or arccos). 6.) You have no idea what your code does, right? :)

[jax]

@kazemakase yes you are right before writing the code i just follow the steps mentioned in the reference that i have added but now a read a lot and have idea what is going on. Actually i m not from the mathematics background :)

[jax]

thanks for the correction you have made. Kindly explain the lines from 10-13 if possible it would be great thanks a lot .

[Eric]

array([f[0] / abVec, f[1] / abVec, f[2] / abVec]) is better spelt f / abVec

[adomas.m]

Thanks Eric...adding to answer.

[adomas.m]

jax lines 10-13 means that vectors ab and bc are normalized and the distance between those are found which then with the help of arccos translates to angle. Maybe someone has a more detailed and elaborate explanation :)

[jax]

thanks Eric for your great help. kindly explain what is going on in line 10 (cosine_angle = -----), and how this code is more idiomatic numpy than adomas.m. I am very new in python. thanks

[Eric]

More idiomatic because of dot and array division, and because important everything into the global namespace (from x import *) is generally frowned upon. The variables in the other answer are poorly named, Division distributes over dot product, so dot(ba / q, bc / r) == dot(ba, bc) / (q * r), which also has the bonus of being faster