Ticket #944: geometry.h

File geometry.h, 2.3 KB (added by Tristan Croll, 8 years ago)
Line 
1
2#ifndef isolde_geometry
3#define isolde_geometry
4
5#include <math.h>
6#include <cstdint>
7
8namespace geometry
9{
10
11template <typename T> void cross_product_3D(T *a, T *b, T* result)
12{
13 result[0] = a[1]*b[2] - a[2]*b[1];
14 result[1] = a[2]*b[0] - a[0]*b[2];
15 result[2] = a[0]*b[1] - a[1]*b[0];
16}
17
18template <typename T> T dot_product_3D(T a[3], T b[3])
19{
20 T accum = 0;
21 for (int i=0; i < 3; ++i) {
22 accum += a[i]*b[i];
23 }
24 return accum;
25}
26
27template <typename T> T l2_norm_3d(T a[3])
28{
29 T accum = 0;
30 for (int i = 0; i < 3; i++) {
31 accum += a[i]*a[i];
32 }
33 return sqrt(accum);
34
35
36}
37
38template <typename T> void normalize_vector_3d(T vector[3])
39{
40 T norm = l2_norm_3d<T>(vector);
41 for (int i = 0; i < 3; ++i) {
42 vector[i] /= norm;
43 }
44}
45
46template <typename T> T dihedral_angle(T p0[3], T p1[3], T p2[3], T p3[3])
47{
48 T b0[3];
49 T b1[3];
50 T b2[3];
51
52 for (int i = 0; i < 3; i++) {
53 b0[i] = p0[i] - p1[i];
54 b1[i] = p2[i] - p1[i];
55 b2[i] = p3[i] - p2[i];
56 }
57
58 T nb1 = l2_norm_3d(b1);
59
60 for (int i =0; i < 3; i++) {
61 b1[i] /= nb1;
62 }
63
64 T dp01 = dot_product_3D<T>(b0, b1);
65 T dp21 = dot_product_3D<T>(b2, b1);
66
67 T v[3];
68 T w[3];
69
70 for (int i = 0; i< 3; i++) {
71 v[i] = b0[i] - dp01*b1[i];
72 w[i] = b2[i] - dp21*b1[i];
73 }
74
75 T x = dot_product_3D<T>(v, w);
76 T xp[3];
77 cross_product_3D<T>(b1, v, xp);
78 T y = dot_product_3D<T>(xp, w);
79
80 return atan2(y, x);
81}
82
83// For use with ChimeraX Point objects
84template <typename T, typename R> const R dihedral_angle(const T& p0, const T& p1, const T& p2, const T& p3)
85{
86 R b0[3];
87 R b1[3];
88 R b2[3];
89
90 for (int i = 0; i < 3; i++) {
91 b0[i] = p0[i] - p1[i];
92 b1[i] = p2[i] - p1[i];
93 b2[i] = p3[i] - p2[i];
94 }
95
96 R nb1 = l2_norm_3d(b1);
97
98 for (int i =0; i < 3; i++) {
99 b1[i] /= nb1;
100 }
101
102 R dp01 = dot_product_3D<R>(b0, b1);
103 R dp21 = dot_product_3D<R>(b2, b1);
104
105 R v[3];
106 R w[3];
107
108 for (int i = 0; i< 3; i++) {
109 v[i] = b0[i] - dp01*b1[i];
110 w[i] = b2[i] - dp21*b1[i];
111 }
112
113 R x = dot_product_3D<R>(v, w);
114 R xp[3];
115 cross_product_3D<R>(b1, v, xp);
116 R y = dot_product_3D<R>(xp, w);
117
118 return atan2(y, x);
119}
120
121} // namespace geometry
122
123#endif //isolde_geometry