Pteros  2.0
Molecular modeling library for human beings!
primitives-3d.h
1 // Copyright Maarten L. Hekkelman, Radboud University 2008-2011.
2 // Distributed under the Boost Software License, Version 1.0.
3 // (See accompanying file LICENSE_1_0.txt or copy at
4 // http://www.boost.org/LICENSE_1_0.txt)
5 //
6 // some data types and routines for working with 3d data
7 
8 #pragma once
9 
10 #include <vector>
11 
12 #include <boost/tr1/tuple.hpp>
13 #include <boost/math/quaternion.hpp>
14 
15 typedef boost::math::quaternion<double> MQuaternion;
16 
17 extern const double kPI;
18 
19 // --------------------------------------------------------------------
20 // The basic point type, can be used to store vectors in 3d space as well of course
21 
22 struct MPoint
23 {
24  MPoint();
25  MPoint(double x, double y, double z);
26  MPoint(const MPoint& rhs);
27 
28  MPoint& operator=(const MPoint& rhs);
29 
30  MPoint& operator+=(const MPoint& rhs);
31  MPoint& operator-=(const MPoint& rhs);
32 
33  MPoint& operator+=(double f);
34  MPoint& operator-=(double f);
35 
36  MPoint& operator*=(double f);
37  MPoint& operator/=(double f);
38 
39  double Normalize();
40  void Rotate(const MQuaternion& q);
41 
42  double mX, mY, mZ;
43 };
44 
45 std::ostream& operator<<(std::ostream& os, const MPoint& pt);
46 MPoint operator+(const MPoint& lhs, const MPoint& rhs);
47 MPoint operator-(const MPoint& lhs, const MPoint& rhs);
48 MPoint operator-(const MPoint& pt);
49 MPoint operator*(const MPoint& pt, double f);
50 MPoint operator/(const MPoint& pt, double f);
51 
52 // --------------------------------------------------------------------
53 // several standard 3d operations
54 
55 double Distance(const MPoint& a, const MPoint& b);
56 double DistanceSquared(const MPoint& a, const MPoint& b);
57 double DotProduct(const MPoint& p1, const MPoint& p2);
58 MPoint CrossProduct(const MPoint& p1, const MPoint& p2);
59 double DihedralAngle(const MPoint& p1, const MPoint& p2, const MPoint& p3, const MPoint& p4);
60 double CosinusAngle(const MPoint& p1, const MPoint& p2, const MPoint& p3, const MPoint& p4);
61 
62 // --------------------------------------------------------------------
63 // We use quaternions to do rotations in 3d space
64 
65 MQuaternion Normalize(MQuaternion q);
66 
67 std::tr1::tuple<double,MPoint> QuaternionToAngleAxis(MQuaternion q);
68 MPoint Centroid(std::vector<MPoint>& points);
69 MPoint CenterPoints(std::vector<MPoint>& points);
70 MQuaternion AlignPoints(const std::vector<MPoint>& a, const std::vector<MPoint>& b);
71 double RMSd(const std::vector<MPoint>& a, const std::vector<MPoint>& b);
72 
73 // --------------------------------------------------------------------
74 // inlines
75 
76 inline
77 MPoint::MPoint()
78  : mX(0)
79  , mY(0)
80  , mZ(0)
81 {
82 }
83 
84 inline
85 MPoint::MPoint(double x, double y, double z)
86  : mX(x)
87  , mY(y)
88  , mZ(z)
89 {
90 }
91 
92 inline
93 MPoint::MPoint(const MPoint& rhs)
94  : mX(rhs.mX)
95  , mY(rhs.mY)
96  , mZ(rhs.mZ)
97 {
98 }
99 
100 inline
101 MPoint& MPoint::operator=(const MPoint& rhs)
102 {
103  mX = rhs.mX;
104  mY = rhs.mY;
105  mZ = rhs.mZ;
106 
107  return *this;
108 }
109 
110 inline
111 MPoint& MPoint::operator+=(const MPoint& rhs)
112 {
113  mX += rhs.mX;
114  mY += rhs.mY;
115  mZ += rhs.mZ;
116 
117  return *this;
118 }
119 
120 inline
121 MPoint& MPoint::operator-=(const MPoint& rhs)
122 {
123  mX -= rhs.mX;
124  mY -= rhs.mY;
125  mZ -= rhs.mZ;
126 
127  return *this;
128 }
129 
130 inline
131 MPoint& MPoint::operator+=(double f)
132 {
133  mX += f;
134  mY += f;
135  mZ += f;
136 
137  return *this;
138 }
139 
140 inline
141 MPoint& MPoint::operator-=(double f)
142 {
143  mX -= f;
144  mY -= f;
145  mZ -= f;
146 
147  return *this;
148 }
149 
150 inline
151 MPoint& MPoint::operator*=(double f)
152 {
153  mX *= f;
154  mY *= f;
155  mZ *= f;
156 
157  return *this;
158 }
159 
160 inline
161 MPoint& MPoint::operator/=(double f)
162 {
163  mX /= f;
164  mY /= f;
165  mZ /= f;
166 
167  return *this;
168 }
169 
170 inline
171 void MPoint::Rotate(const MQuaternion& q)
172 {
173  MQuaternion p(0, mX, mY, mZ);
174 
175  p = q * p * conj(q);
176 
177  mX = p.R_component_2();
178  mY = p.R_component_3();
179  mZ = p.R_component_4();
180 }
181 
182 inline double DotProduct(const MPoint& a, const MPoint& b)
183 {
184  return a.mX * b.mX + a.mY * b.mY + a.mZ * b.mZ;
185 }
186 
187 inline MPoint CrossProduct(const MPoint& a, const MPoint& b)
188 {
189  return MPoint(a.mY * b.mZ - b.mY * a.mZ,
190  a.mZ * b.mX - b.mZ * a.mX,
191  a.mX * b.mY - b.mX * a.mY);
192 }
193 
194 inline double DistanceSquared(const MPoint& a, const MPoint& b)
195 {
196  return
197  (a.mX - b.mX) * (a.mX - b.mX) +
198  (a.mY - b.mY) * (a.mY - b.mY) +
199  (a.mZ - b.mZ) * (a.mZ - b.mZ);
200 }
201 
202 inline double Distance(const MPoint& a, const MPoint& b)
203 {
204  return sqrt(
205  (a.mX - b.mX) * (a.mX - b.mX) +
206  (a.mY - b.mY) * (a.mY - b.mY) +
207  (a.mZ - b.mZ) * (a.mZ - b.mZ));
208 }