IMP logo
IMP Reference Guide  develop.d247202c1c,2019/10/23
The Integrative Modeling Platform
Rotation3D.h
Go to the documentation of this file.
1 /**
2  * \file IMP/algebra/Rotation3D.h \brief Simple 3D rotation class.
3  *
4  * Copyright 2007-2019 IMP Inventors. All rights reserved.
5  *
6  */
7 
8 #ifndef IMPALGEBRA_ROTATION_3D_H
9 #define IMPALGEBRA_ROTATION_3D_H
10 
11 #include <IMP/algebra/algebra_config.h>
12 #include "Vector3D.h"
13 #include "utility.h"
14 #include "constants.h"
15 #include "GeometricPrimitiveD.h"
16 #include <Eigen/Dense>
17 
18 #include <IMP/log.h>
19 #include <cmath>
20 #include <iostream>
21 #include <algorithm>
22 
23 IMPALGEBRA_BEGIN_NAMESPACE
24 
25 typedef Vector4D Rotation3DAdjoint;
26 typedef std::pair<Vector3D, Rotation3DAdjoint> RotatedVector3DAdjoint;
27 typedef std::pair<Rotation3DAdjoint, Rotation3DAdjoint>
28  ComposeRotation3DAdjoint;
29 
30 #if !defined(IMP_DOXYGEN) && !defined(SWIG)
31 class Rotation3D;
32 Rotation3D compose(const Rotation3D &a, const Rotation3D &b);
33 #endif
34 
35 //! 3D rotation class.
36 /** Rotations are currently represented using quaternions and a cached
37  copy of the rotation matrix. The quaternion allows for fast and
38  stable composition and the cached rotation matrix means that
39  rotations are performed quickly. See
40  http://en.wikipedia.org/wiki/Quaternions_and_spatial_rotation for
41  a comparison of different implementations of rotations.
42 
43  Currently the rotation can be initialized from either:
44  - XYZ Euler angles
45  - Rotation Matrix
46  - Quaternion
47  - angle/axis representation
48 
49  \geometry
50 */
51 class IMPALGEBRAEXPORT Rotation3D : public GeometricPrimitiveD<3> {
52  VectorD<4> v_;
53  mutable bool has_cache_;
54  mutable Vector3D matrix_[3];
55  IMP_NO_SWIG(friend Rotation3D compose(const Rotation3D &a,
56  const Rotation3D &b));
57  void fill_cache() const {
58  IMP_USAGE_CHECK(get_is_valid(),
59  "Attempting to apply uninitialized rotation");
60  has_cache_ = true;
61  double v0s = get_squared(v_[0]);
62  double v1s = get_squared(v_[1]);
63  double v2s = get_squared(v_[2]);
64  double v3s = get_squared(v_[3]);
65  double v12 = v_[1] * v_[2];
66  double v01 = v_[0] * v_[1];
67  double v02 = v_[0] * v_[2];
68  double v23 = v_[2] * v_[3];
69  double v03 = v_[0] * v_[3];
70  double v13 = v_[1] * v_[3];
71  matrix_[0] =
72  Vector3D(v0s + v1s - v2s - v3s, 2 * (v12 - v03), 2 * (v13 + v02));
73  matrix_[1] =
74  Vector3D(2 * (v12 + v03), v0s - v1s + v2s - v3s, 2 * (v23 - v01));
75  matrix_[2] =
76  Vector3D(2 * (v13 - v02), 2 * (v23 + v01), v0s - v1s - v2s + v3s);
77  }
78 
79  // When copying preserve the cached rotation matrix, if present
80  void copy_cache(const Rotation3D &rot) {
81  has_cache_ = rot.has_cache_;
82  if (has_cache_) {
83  matrix_[0] = rot.matrix_[0];
84  matrix_[1] = rot.matrix_[1];
85  matrix_[2] = rot.matrix_[2];
86  }
87  }
88 
89  public:
90  //! Rotation3D copy constructor
91  Rotation3D(const Rotation3D &rot) : v_(rot.v_) {
92  copy_cache(rot);
93  }
94 
95  Rotation3D &operator=(const Rotation3D &rot) {
96  v_ = rot.v_;
97  copy_cache(rot);
98  return *this;
99  }
100 
101  //! Create a rotation from a vector of 4 quaternion coefficients.
102  //! @note: use assume_normalized with care - inputting an unnormalized
103  //! vector would result in unexpected results if it is true
104  explicit Rotation3D(const VectorD<4> &v,
105  bool assume_normalized=false)
106  : v_(assume_normalized ? v : v.get_unit_vector()),
107  has_cache_(false) {}
108 
109  //! Create an invalid rotation
110  Rotation3D() : v_(0, 0, 0, 0), has_cache_(false) {}
111  //! Create a rotation from a quaternion
112  /** \throw ValueException if the rotation is not a unit vector.
113  */
114  Rotation3D(double a, double b, double c, double d)
115  : v_(a, b, c, d), has_cache_(false) {
117  v_.get_squared_magnitude(), 1.0,
118  "Attempting to construct a rotation from a "
119  << " non-quaternion value. The coefficient vector"
120  << " must have a length of 1. Got: " << a << " " << b << " " << c
121  << " " << d << " gives " << v_.get_squared_magnitude());
122  if (a < 0) {
123  // make them canonical
124  v_ = -v_;
125  }
126  }
127  ~Rotation3D(){}
128 
129 #ifndef IMP_DOXYGEN
130  Vector3D get_rotated_no_cache(const Vector3D &o) const {
131  IMP_USAGE_CHECK(get_is_valid(),
132  "Attempting to access uninitialized rotation");
133  return Vector3D(
134  (v_[0] * v_[0] + v_[1] * v_[1] - v_[2] * v_[2] - v_[3] * v_[3]) * o[0] +
135  2 * (v_[1] * v_[2] - v_[0] * v_[3]) * o[1] +
136  2 * (v_[1] * v_[3] + v_[0] * v_[2]) * o[2],
137  2 * (v_[1] * v_[2] + v_[0] * v_[3]) * o[0] +
138  (v_[0] * v_[0] - v_[1] * v_[1] + v_[2] * v_[2] - v_[3] * v_[3]) *
139  o[1] +
140  2 * (v_[2] * v_[3] - v_[0] * v_[1]) * o[2],
141  2 * (v_[1] * v_[3] - v_[0] * v_[2]) * o[0] +
142  2 * (v_[2] * v_[3] + v_[0] * v_[1]) * o[1] +
143  (v_[0] * v_[0] - v_[1] * v_[1] - v_[2] * v_[2] + v_[3] * v_[3]) *
144  o[2]);
145  }
146 
147  //! Get only the requested rotation coordinate of the vector
148  double get_rotated_one_coordinate_no_cache(const Vector3D &o,
149  unsigned int coord) const {
150  IMP_USAGE_CHECK(get_is_valid(),
151  "Attempting to apply uninitialized rotation");
152  switch (coord) {
153  case 0:
154  return (v_[0] * v_[0] + v_[1] * v_[1] - v_[2] * v_[2] - v_[3] * v_[3]) *
155  o[0] +
156  2 * (v_[1] * v_[2] - v_[0] * v_[3]) * o[1] +
157  2 * (v_[1] * v_[3] + v_[0] * v_[2]) * o[2];
158  break;
159  case 1:
160  return 2 * (v_[1] * v_[2] + v_[0] * v_[3]) * o[0] +
161  (v_[0] * v_[0] - v_[1] * v_[1] + v_[2] * v_[2] - v_[3] * v_[3]) *
162  o[1] +
163  2 * (v_[2] * v_[3] - v_[0] * v_[1]) * o[2];
164 
165  break;
166  case 2:
167  return 2 * (v_[1] * v_[3] - v_[0] * v_[2]) * o[0] +
168  2 * (v_[2] * v_[3] + v_[0] * v_[1]) * o[1] +
169  (v_[0] * v_[0] - v_[1] * v_[1] - v_[2] * v_[2] + v_[3] * v_[3]) *
170  o[2];
171  break;
172  default:
173  IMP_THROW("Out of range coordinate " << coord, IndexException);
174  }
175  }
176 #endif
177  //! Rotate a vector around the origin
178  Vector3D get_rotated(const Vector3D &o) const {
179  if (!has_cache_) fill_cache();
180  return Vector3D(o * matrix_[0], o * matrix_[1], o * matrix_[2]);
181  }
182 
183 #ifndef SWIG
184  //! Get adjoint of inputs to `get_rotated` from adjoint of output
185  /** Compute the adjoint (reverse-mode sensitivity) of input vector
186  to `get_rotated` and this rotation from the adjoint of the
187  output vector.
188  */
189  void get_rotated_adjoint(const Vector3D &v, const Vector3D &Dw,
190  Vector3D *Dv, Rotation3DAdjoint *Dr) const;
191 #endif
192 
193  //! Get adjoint of inputs to `get_rotated` from adjoint of output
194  /** Compute the adjoint (reverse-mode sensitivity) of input vector
195  to `get_rotated` and this rotation from the adjoint of the
196  output vector.
197  */
198  RotatedVector3DAdjoint
199  get_rotated_adjoint(const Vector3D &v, const Vector3D &Dw) const;
200 
201  //! Get only the requested rotation coordinate of the vector
203  unsigned int coord) const {
204  if (!has_cache_) fill_cache();
205  return o * matrix_[coord];
206  }
207 
208  //! Rotate a vector around the origin
209  Vector3D operator*(const Vector3D &v) const { return get_rotated(v); }
210  Vector3D get_rotation_matrix_row(int i) const {
211  IMP_USAGE_CHECK((i >= 0) && (i <= 2), "row index out of range");
212  if (!has_cache_) fill_cache();
213  return matrix_[i];
214  }
215  IMP_SHOWABLE_INLINE(Rotation3D,
216  { out << v_[0] << " " << v_[1] << " " << v_[2] << " " << v_[3]; });
217 
218  //! Return the rotation which undoes this rotation.
219  inline Rotation3D get_inverse() const {
220  IMP_USAGE_CHECK(get_is_valid(),
221  "Attempting to invert uninitialized rotation");
222  Rotation3D ret(v_[0], -v_[1], -v_[2], -v_[3]);
223  return ret;
224  }
225 
226  //! Return the quaternion so that it can be stored
227  /** Note that there is no guarantee on which of the two
228  equivalent quaternions is returned.
229  */
230  const Vector4D &get_quaternion() const {
231  IMP_USAGE_CHECK(get_is_valid(),
232  "Attempting to access uninitialized rotation");
233  return v_;
234  }
235 
236  //! Multiply two rotations
237  Rotation3D operator*(const Rotation3D &q) const {
238  IMP_USAGE_CHECK(get_is_valid(),
239  "Attempting to compose uninitialized rotation");
240  return compose(*this, q);
241  }
242 
243  //! Compute the rotation which when composed with r gives this
244  Rotation3D operator/(const Rotation3D &r) const {
245  IMP_USAGE_CHECK(get_is_valid(),
246  "Attempting to compose uninitialized rotation");
247  return compose(*this, r.get_inverse());
248  }
249 
250  const Rotation3D &operator/=(const Rotation3D &r) {
251  *this = *this / r;
252  return *this;
253  }
254 
255  //! Return the gradient of rotated vector wrt the ith quaternion element.
256  /** Given the rotation \f$x = R(q) v\f$, where \f$v\f$ is a vector,
257  \f$q=(q_0,q_1,q_2,q_3)\f$ is the quaternion of the rotation with
258  elements \f$q_i\f$, and \f$R(q)\f$ is the corresponding rotation matrix,
259  the function returns the gradient \f$\nabla_{q_i} x\f$.
260 
261  This function just returns a single column from `get_jacobian()`, so it is
262  more efficient to call that function if all columns are needed.
263 
264  \param[in] v vector to be rotated by rotation \f$R(q)\f$
265  \param[in] wrt_unnorm Gradient is computed wrt unnormalized quaternion.
266  Rotation includes a normalization operation, and
267  the gradient is projected to the tangent space at
268  \f$q\f$.
269  */
270  Vector3D get_gradient_of_rotated(const Vector3D &v, unsigned int i,
271  bool wrt_unnorm = false) const;
272 
273  IMPALGEBRA_DEPRECATED_METHOD_DECL(2.12)
274  Vector3D get_derivative(const Vector3D &v, unsigned int i,
275  bool wrt_unnorm = true) const;
276 
277  //! Return the Jacobian of rotated vector wrt the quaternion.
278  /** Given the rotation \f$x = R(q) v\f$, where \f$v\f$ is a vector,
279  \f$q\f$ is the quaternion of the rotation, and \f$R(q)\f$ is the
280  corresponding rotation matrix, the function returns the 3x4 matrix
281  \f$J\f$ with elements \f$J_{ij}=\frac{\partial x_i}{\partial q_j}\f$.
282 
283  \param[in] v vector to be rotated by rotation \f$R(q)\f$
284  \param[in] wrt_unnorm Jacobian is computed wrt unnormalized quaternion.
285  Rotation includes a normalization operation, and
286  the columns are projected to the tangent space at
287  \f$q\f$.
288  */
289  Eigen::MatrixXd get_jacobian_of_rotated(
290  const Eigen::Vector3d &v, bool wrt_unnorm = false) const;
291 
292  IMPALGEBRA_DEPRECATED_METHOD_DECL(2.12)
293  Eigen::MatrixXd get_gradient(
294  const Eigen::Vector3d &v, bool wrt_unnorm = true) const;
295 
296  /** Return true is the rotation is valid, false if
297  invalid or not initialized (e.g., only initialized by
298  the empty constructor)
299  */
300  bool get_is_valid() const {
301  return v_.get_squared_magnitude() > 0; // TODO: add that magnitude ~ 1?
302  }
303 };
304 
306 
307 
308 //! Get gradient of quaternion product with respect to first quaternion.
309 /** Given the rotation \f$R(p)\f$ followed by \f$R(q)\f$, where \f$p\f$ and
310  \f$q\f$ are quaternions, the quaternion of the composed rotation
311  \f$R(s)=R(q) R(p)\f$ can be expressed through the Hamilton product of the
312  two quaternions \f$s(q,p) = q p\f$. This function returns the matrix
313  \f$J\f$ with elements \f$J_{ij}=\frac{\partial s_i}{\partial q_j}\f$.
314 
315  \param[in] q rotation corresponding to first quaternion
316  \param[in] p rotation corresponding to second quaternion
317  \param[in] wrt_unnorm Jacobian is computed wrt unnormalized quaternion.
318  Rotation includes a normalization operation, and
319  the columns are projected to the tangent space at
320  \f$q\f$.
321  */
322 IMPALGEBRAEXPORT Eigen::MatrixXd get_jacobian_of_composed_wrt_first(
323  const Rotation3D &q, const Rotation3D &p, bool wrt_unnorm = false);
324 
325 IMPALGEBRA_DEPRECATED_FUNCTION_DECL(2.12)
326 IMPALGEBRAEXPORT Eigen::MatrixXd
327  get_gradient_of_composed_with_respect_to_first(
328  const Rotation3D &q, const Rotation3D &p, bool wrt_unnorm = true);
329 
330 
331 //! Get Jacobian of quaternion product with respect to second quaternion.
332 /** Given the rotation \f$R(p)\f$ followed by \f$R(q)\f$, where \f$p\f$ and
333  \f$q\f$ are quaternions, the quaternion of the composed rotation
334  \f$R(s)=R(q) R(p)\f$ can be expressed through the Hamilton product of the
335  two quaternions \f$s(q,p) = q p\f$. This function returns the matrix
336  \f$J\f$ with elements \f$J_{ij}=\frac{\partial s_i}{\partial p_j}\f$.
337 
338  \param[in] q rotation corresponding to first quaternion
339  \param[in] p rotation corresponding to second quaternion
340  \param[in] wrt_unnorm Jacobian is computed wrt unnormalized quaternion.
341  Rotation includes a normalization operation, and
342  the columns are projected to the tangent space at
343  \f$p\f$.
344  */
345 IMPALGEBRAEXPORT Eigen::MatrixXd get_jacobian_of_composed_wrt_second(
346  const Rotation3D &q, const Rotation3D &p, bool wrt_unnorm = false);
347 
348 IMPALGEBRA_DEPRECATED_FUNCTION_DECL(2.12)
349 IMPALGEBRAEXPORT Eigen::MatrixXd
350  get_gradient_of_composed_with_respect_to_second(
351  const Rotation3D &q, const Rotation3D &p, bool wrt_unnorm = true);
352 
353 
354 //! Return a rotation that does not do anything
355 /** \see Rotation3D */
356 inline Rotation3D get_identity_rotation_3d() { return Rotation3D(1, 0, 0, 0); }
357 
358 //! Return a distance between the two rotations
359 /** The distance runs between 0 and 1. More precisely,
360  the distance returned is distance between the two
361  quaternion vectors properly normalized, divided
362  by sqrt(2).
363 
364  A vector with distance d from the unit vector
365  represents a rotation of
366 
367  \f$ \theta= \cos^{-1}\left(1-\sqrt{2}d\right)\f$
368  \see Rotation3D
369 */
370 inline double get_distance(const Rotation3D &r0, const Rotation3D &r1) {
371  double dot =
372  (r0.get_quaternion() - r1.get_quaternion()).get_squared_magnitude();
373  double odot =
374  (r0.get_quaternion() + r1.get_quaternion()).get_squared_magnitude();
375  double ans = std::min(dot, odot);
376  // TODO: barak - added static for efficiency
377  static const double s2 = std::sqrt(2.0);
378  double ret = ans / s2;
379  return std::max(std::min(ret, 1.0), 0.0);
380 }
381 
382 //! Generate a Rotation3D object from a rotation around an axis
383 //! that is assumed to be normalized
384 /**
385  \param[in] axis_norm the normalized rotation axis passing through (0,0,0)
386  \param[in] angle the rotation angle in radians in the
387  clockwise direction
388  \note http://en.wikipedia.org/wiki/Rotation_matrix
389  \note www.euclideanspace.com/maths/geometry/rotations/conversions/
390  angleToQuaternion/index.htm
391  \see Rotation3D
392 */
393 inline Rotation3D
395 (const Vector3D &axis_norm, double angle)
396 {
397  IMP_USAGE_CHECK(axis_norm.get_magnitude() - 1.0 < 1e-6,
398  "expected normalized vector as axis of rotation");
399  double s = std::sin(angle / 2);
400  double a, b, c, d;
401  a = std::cos(angle / 2);
402  b = axis_norm[0] * s;
403  c = axis_norm[1] * s;
404  d = axis_norm[2] * s;
405  return Rotation3D(a, b, c, d);
406 }
407 
408 
409 //! Generate a Rotation3D object from a rotation around an axis
410 /**
411  \param[in] axis the rotation axis passes through (0,0,0)
412  \param[in] angle the rotation angle in radians in the
413  clockwise direction
414  \note http://en.wikipedia.org/wiki/Rotation_matrix
415  \note www.euclideanspace.com/maths/geometry/rotations/conversions/
416  angleToQuaternion/index.htm
417  \see Rotation3D
418 */
419 inline Rotation3D
421 (const Vector3D &axis, double angle)
422 {
423  // normalize the vector
424  Vector3D axis_norm = axis.get_unit_vector();
425  return get_rotation_about_normalized_axis(axis_norm, angle);
426 }
427 
428 //! Create a rotation from the first vector to the second one.
429 /** \see Rotation3D
430  */
431 IMPALGEBRAEXPORT Rotation3D
433 
434 //! Generate a Rotation3D object from a rotation matrix
435 /**
436  \see Rotation3D
437 */
438 IMPALGEBRAEXPORT Rotation3D
439  get_rotation_from_matrix(double m00, double m01, double m02, double m10,
440  double m11, double m12, double m20, double m21,
441  double m22);
442 
443 //! Generate a Rotation3D object from a rotation matrix
444 /**
445  \see Rotation3D
446 */
447 IMPALGEBRAEXPORT Rotation3D get_rotation_from_matrix(Eigen::Matrix3d m);
448 
449 //! Pick a rotation at random from all possible rotations
450 /** \see Rotation3D */
451 IMPALGEBRAEXPORT Rotation3D get_random_rotation_3d();
452 
453 //! Pick a rotation at random near the provided one
454 /** This method generates a rotation that is within the provided
455  distance of center.
456  \param[in] center The center of the rotational volume
457  \param[in] distance See
458  get_distance(const Rotation3D&,const Rotation3D&)
459  for a full definition.
460 
461  \note The cost of this operation increases as distance goes to 0.
462 
463  \see Rotation3D
464 */
465 IMPALGEBRAEXPORT Rotation3D
466  get_random_rotation_3d(const Rotation3D &center, double distance);
467 
468 //! Cover the space of rotations evenly
469 /** If you care about the distance between samples instead of the number
470  of samples, the "surface area" of the set of rotations is pi^2. If
471  you allocate each sample a volume of 4/3 pi d^3 (to space them d apart),
472  Then you want 3/4 pi/d^3 points.
473 
474  Creates at least num_points rotations.
475 */
476 IMPALGEBRAEXPORT Rotation3Ds
477  get_uniform_cover_rotations_3d(unsigned int num_points);
478 
479 //! Generates a nondegenerate set of Euler angles with a delta resolution
480 /**
481 \param[in] delta sample every delta angles in radians.
482  */
484  double delta);
485 
486 //! Compute a rotation from an unnormalized quaternion
487 /** \see Rotation3D */
489  VectorD<4> uv = v.get_unit_vector();
490  return Rotation3D(uv[0], uv[1], uv[2], uv[3]);
491 }
492 
493 /** \see Rotation3D
494  */
495 inline Rotation3D compose(const Rotation3D &a, const Rotation3D &b) {
496  return Rotation3D(a.v_[0] * b.v_[0] - a.v_[1] * b.v_[1] - a.v_[2] * b.v_[2] -
497  a.v_[3] * b.v_[3],
498  a.v_[0] * b.v_[1] + a.v_[1] * b.v_[0] + a.v_[2] * b.v_[3] -
499  a.v_[3] * b.v_[2],
500  a.v_[0] * b.v_[2] - a.v_[1] * b.v_[3] + a.v_[2] * b.v_[0] +
501  a.v_[3] * b.v_[1],
502  a.v_[0] * b.v_[3] + a.v_[1] * b.v_[2] - a.v_[2] * b.v_[1] +
503  a.v_[3] * b.v_[0]);
504 }
505 
506 #ifndef SWIG
507 //! Get adjoint of inputs to `compose` from adjoint of output
508 /** Compute the adjoint (reverse-mode sensitivity) of input rotations
509  to `compose` from the adjoint of the output rotation.
510  */
511 IMPALGEBRAEXPORT void
512 compose_adjoint(const Rotation3D &A, const Rotation3D &B, Vector4D DC,
513  Rotation3DAdjoint *DA, Rotation3DAdjoint *DB);
514 #endif
515 
516 //! Get adjoint of inputs to `compose` from adjoint of output
517 /** Compute the adjoint (reverse-mode sensitivity) of input rotations
518  to `compose` from the adjoint of the output rotation.
519  */
520 IMPALGEBRAEXPORT ComposeRotation3DAdjoint
521 compose_adjoint(const Rotation3D &A, const Rotation3D &B, const Rotation3DAdjoint &DC);
522 
523 /** \name Euler Angles
524  There are many conventions for how to define Euler angles, based on choices
525  of which of the x,y,z axis to use in what order and whether the rotation
526  axis is in the body frame (and hence affected by previous rotations) or in
527  in a fixed frame. See
528  http://en.wikipedia.org/wiki/Conversion_between_quaternions_and_Euler_angles
529  for a general description.
530 
531  - All Euler angles are specified in radians.
532  - Only x-y-z order is currently supported.
533  - To convert Euler angles in a different order, one can compose a Rotation3D
534  from three rotations using get_rotation_about_axis function.
535  @{
536 */
537 
538 //! Initialize a rotation in x-y-z order from three angles
539 /** \param[in] xr Rotation around the X axis in radians
540  \param[in] yr Rotation around the Y axis in radians
541  \param[in] zr Rotation around the Z axis in radians
542  \note The three rotations are represented in the original (fixed)
543  coordinate frame.
544  \see Rotation3D
545  \see FixedXYZ
546 */
547 IMPALGEBRAEXPORT Rotation3D
548  get_rotation_from_fixed_xyz(double xr, double yr, double zr);
549 
550 //! Initialize a rotation from Euler angles
551 /**
552  \param[in] phi Rotation around the Z axis in radians
553  \param[in] theta Rotation around the X axis in radians
554  \param[in] psi Rotation around the Z axis in radians
555  \note The first rotation is by an angle phi about the z-axis.
556  The second rotation is by an angle theta in [0,pi] about the
557  former x-axis , and the third rotation is by an angle psi
558  about the former z-axis.
559  \see Rotation3D
560 */
561 IMPALGEBRAEXPORT Rotation3D
562  get_rotation_from_fixed_zxz(double phi, double theta, double psi);
563 
564 //! Generate a rotation object from Euler Angles
565 /** \note The first rotation is by an angle about the z-axis.
566  The second rotation is by an angle about the new y-axis.
567  The third rotation is by an angle about the new z-axis.
568  \param[in] Rot First Euler angle (radians) defining the rotation (Z axis)
569  \param[in] Tilt Second Euler angle (radians) defining the rotation (Y axis)
570  \param[in] Psi Third Euler angle (radians) defining the rotation (Z axis)
571  \see Rotation3D
572 */
573 IMPALGEBRAEXPORT Rotation3D
574  get_rotation_from_fixed_zyz(double Rot, double Tilt, double Psi);
575 
576 //! A simple class for returning XYZ Euler angles
577 class FixedXYZ : public GeometricPrimitiveD<3> {
578  double v_[3];
579 
580  public:
581  FixedXYZ() {}
582  FixedXYZ(double x, double y, double z) {
583  v_[0] = x;
584  v_[1] = y;
585  v_[2] = z;
586  }
587  double get_x() const { return v_[0]; }
588  double get_y() const { return v_[1]; }
589  double get_z() const { return v_[2]; }
591  { out << v_[0] << " " << v_[1] << " " << v_[2]; });
592 };
593 
595 
596 //! The inverse of rotation_from_fixed_xyz()
597 /**
598  \see rotation_from_fixed_xyz()
599  \see Rotation3D
600  \see FixesXYZ
601 */
602 IMPALGEBRAEXPORT FixedXYZ get_fixed_xyz_from_rotation(const Rotation3D &r);
603 
604 /** @}*/
605 
606 //! Interpolate between two rotations
607 /** It f ==0, return b, if f==1 return a.
608  \see Rotation3D */
610  double f) {
611  VectorD<4> bq = b.get_quaternion(), aq = a.get_quaternion();
612  if (bq * aq < 0) bq = -bq;
613  return Rotation3D(f * aq + (1 - f) * bq);
614 }
615 
616 /** Return the rotation which takes the native x and y axes to the
617  given x and y axes.
618  The two axes must be perpendicular unit vectors.
619 */
620 IMPALGEBRAEXPORT Rotation3D
621  get_rotation_from_x_y_axes(const Vector3D &x, const Vector3D &y);
622 
623 //! Decompose a Rotation3D object into a rotation around an axis
624 /** For all identity rotations, returns the axis [1,0,0] and the angle 0.0.
625 
626  \note http://en.wikipedia.org/wiki/Rotation_matrix
627  \note www.euclideanspace.com/maths/geometry/rotations/conversions/
628  angleToQuaternion/index.htm
629  \see Rotation3D
630 
631  @return axis direction and rotation about the axis in Radians
632 */
633 IMPALGEBRAEXPORT std::pair<Vector3D, double> get_axis_and_angle(
634  const Rotation3D &rot);
635 
636 typedef std::pair<Vector3D, double> AxisAnglePair;
637 #ifndef IMP_DOXYGEN
638 typedef Vector<AxisAnglePair> AxisAnglePairs;
639 #endif
640 
641 IMPALGEBRA_END_NAMESPACE
642 #endif /* IMPALGEBRA_ROTATION_3D_H */
Vector3D get_rotated(const Vector3D &o) const
Rotate a vector around the origin.
Definition: Rotation3D.h:178
Rotation3D get_rotation_about_normalized_axis(const Vector3D &axis_norm, double angle)
Definition: Rotation3D.h:395
Rotation3D(const Rotation3D &rot)
Rotation3D copy constructor.
Definition: Rotation3D.h:91
Base class for geometric types.
#define IMP_SHOWABLE_INLINE(Name, how_to_show)
Declare the methods needed by an object that can be printed.
Rotation3D get_rotation_from_x_y_axes(const Vector3D &x, const Vector3D &y)
Rotation3D get_rotation_from_fixed_xyz(double xr, double yr, double zr)
Initialize a rotation in x-y-z order from three angles.
Rotation3D()
Create an invalid rotation.
Definition: Rotation3D.h:110
#define IMP_USAGE_CHECK_FLOAT_EQUAL(expra, exprb, message)
Definition: check_macros.h:178
Rotation2D compose(const Rotation2D &a, const Rotation2D &b)
Compose two rotations a and b.
Definition: Rotation2D.h:108
algebra::Rotation3Ds get_uniformly_sampled_rotations(double delta)
Generates a nondegenerate set of Euler angles with a delta resolution.
Rotation3D get_random_rotation_3d(const Rotation3D &center, double distance)
Pick a rotation at random near the provided one.
Rotation3D operator/(const Rotation3D &r) const
Compute the rotation which when composed with r gives this.
Definition: Rotation3D.h:244
A simple class for returning XYZ Euler angles.
Definition: Rotation3D.h:577
Rotation3D(const VectorD< 4 > &v, bool assume_normalized=false)
Definition: Rotation3D.h:104
Vector3D operator*(const Vector3D &v) const
Rotate a vector around the origin.
Definition: Rotation3D.h:209
Rotation3D compose(const Rotation3D &a, const Rotation3D &b)
Definition: Rotation3D.h:495
double get_rotated_one_coordinate(const Vector3D &o, unsigned int coord) const
Get only the requested rotation coordinate of the vector.
Definition: Rotation3D.h:202
VT get_unit_vector(VT vt)
Returns a unit vector pointing at the same direction as this vector.
Definition: VectorBaseD.h:229
Eigen::MatrixXd get_jacobian_of_composed_wrt_first(const Rotation3D &q, const Rotation3D &p, bool wrt_unnorm=false)
Get gradient of quaternion product with respect to first quaternion.
const Vector4D & get_quaternion() const
Return the quaternion so that it can be stored.
Definition: Rotation3D.h:230
VectorD< 4 > Vector4D
Definition: VectorD.h:425
Rotation3D get_rotation_about_axis(const Vector3D &axis, double angle)
Generate a Rotation3D object from a rotation around an axis.
Definition: Rotation3D.h:421
#define IMP_VALUES(Name, PluralName)
Define the type for storing sets of values.
Definition: value_macros.h:23
Base class for geometric types.
Functions to deal with very common math operations.
A Cartesian vector in D-dimensions.
Definition: VectorD.h:52
FixedXYZ get_fixed_xyz_from_rotation(const Rotation3D &r)
The inverse of rotation_from_fixed_xyz()
Rotation3D get_inverse() const
Return the rotation which undoes this rotation.
Definition: Rotation3D.h:219
Rotation3D get_rotation_from_fixed_zyz(double Rot, double Tilt, double Psi)
Generate a rotation object from Euler Angles.
std::pair< Vector3D, double > get_axis_and_angle(const Rotation3D &rot)
Decompose a Rotation3D object into a rotation around an axis.
double get_distance(const Rotation3D &r0, const Rotation3D &r1)
Return a distance between the two rotations.
Definition: Rotation3D.h:370
Rotation3D get_interpolated(const Rotation3D &a, const Rotation3D &b, double f)
Interpolate between two rotations.
Definition: Rotation3D.h:609
#define IMP_NO_SWIG(x)
Hide the line when SWIG is compiled or parses it.
Definition: swig_macros.h:18
Rotation3D get_rotation_from_vector4d(const VectorD< 4 > &v)
Compute a rotation from an unnormalized quaternion.
Definition: Rotation3D.h:488
Rotation3D get_rotation_taking_first_to_second(const Vector3D &v1, const Vector3D &v2)
Create a rotation from the first vector to the second one.
Rotation3Ds get_uniform_cover_rotations_3d(unsigned int num_points)
Cover the space of rotations evenly.
Various useful constants.
Rotation3D operator*(const Rotation3D &q) const
Multiply two rotations.
Definition: Rotation3D.h:237
3D rotation class.
Definition: Rotation3D.h:51
Eigen::MatrixXd get_jacobian_of_composed_wrt_second(const Rotation3D &q, const Rotation3D &p, bool wrt_unnorm=false)
Get Jacobian of quaternion product with respect to second quaternion.
#define IMP_THROW(message, exception_name)
Throw an exception with a message.
Definition: check_macros.h:50
Rotation3D(double a, double b, double c, double d)
Create a rotation from a quaternion.
Definition: Rotation3D.h:114
VectorD< 3 > Vector3D
Definition: VectorD.h:421
Rotation3D get_identity_rotation_3d()
Return a rotation that does not do anything.
Definition: Rotation3D.h:356
ComposeRotation3DAdjoint compose_adjoint(const Rotation3D &A, const Rotation3D &B, const Rotation3DAdjoint &DC)
Get adjoint of inputs to compose from adjoint of output.
Simple 3D vector class.
#define IMP_USAGE_CHECK(expr, message)
A runtime test for incorrect usage of a class or method.
Definition: check_macros.h:168
Logging and error reporting support.
Rotation3D get_rotation_from_fixed_zxz(double phi, double theta, double psi)
Initialize a rotation from Euler angles.
Rotation3D get_rotation_from_matrix(Eigen::Matrix3d m)
Generate a Rotation3D object from a rotation matrix.
IMP::Vector< Rotation3D > Rotation3Ds
Definition: Rotation3D.h:305