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