geo_quaternion.cpp
Go to the documentation of this file.
1 // This file is part of Eigen, a lightweight C++ template library
2 // for linear algebra.
3 //
4 // Copyright (C) 2008-2009 Gael Guennebaud <gael.guennebaud@inria.fr>
5 // Copyright (C) 2009 Mathieu Gautier <mathieu.gautier@cea.fr>
6 //
7 // This Source Code Form is subject to the terms of the Mozilla
8 // Public License v. 2.0. If a copy of the MPL was not distributed
9 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
10 
11 #include "main.h"
12 #include <Eigen/Geometry>
13 #include <Eigen/LU>
14 #include <Eigen/SVD>
15 #include "AnnoyingScalar.h"
16 
17 template<typename T> T bounded_acos(T v)
18 {
19  using std::acos;
20  using std::min;
21  using std::max;
22  return acos((max)(T(-1),(min)(v,T(1))));
23 }
24 
25 template<typename QuatType> void check_slerp(const QuatType& q0, const QuatType& q1)
26 {
27  using std::abs;
28  typedef typename QuatType::Scalar Scalar;
29  typedef AngleAxis<Scalar> AA;
30 
31  Scalar largeEps = test_precision<Scalar>();
32 
33  Scalar theta_tot = AA(q1*q0.inverse()).angle();
34  if(theta_tot>Scalar(EIGEN_PI))
35  theta_tot = Scalar(2.)*Scalar(EIGEN_PI)-theta_tot;
36  for(Scalar t=0; t<=Scalar(1.001); t+=Scalar(0.1))
37  {
38  QuatType q = q0.slerp(t,q1);
39  Scalar theta = AA(q*q0.inverse()).angle();
40  VERIFY(abs(q.norm() - 1) < largeEps);
41  if(theta_tot==0) VERIFY(theta_tot==0);
42  else VERIFY(abs(theta - t * theta_tot) < largeEps);
43  }
44 }
45 
46 template<typename Scalar, int Options> void quaternion(void)
47 {
48  /* this test covers the following files:
49  Quaternion.h
50  */
51  using std::abs;
52  typedef Matrix<Scalar,3,1> Vector3;
53  typedef Matrix<Scalar,3,3> Matrix3;
54  typedef Quaternion<Scalar,Options> Quaternionx;
55  typedef AngleAxis<Scalar> AngleAxisx;
56 
57  Scalar largeEps = test_precision<Scalar>();
58  if (internal::is_same<Scalar,float>::value)
59  largeEps = Scalar(1e-3);
60 
61  Scalar eps = internal::random<Scalar>() * Scalar(1e-2);
62 
63  Vector3 v0 = Vector3::Random(),
64  v1 = Vector3::Random(),
65  v2 = Vector3::Random(),
66  v3 = Vector3::Random();
67 
68  Scalar a = internal::random<Scalar>(-Scalar(EIGEN_PI), Scalar(EIGEN_PI)),
69  b = internal::random<Scalar>(-Scalar(EIGEN_PI), Scalar(EIGEN_PI));
70 
71  // Quaternion: Identity(), setIdentity();
72  Quaternionx q1, q2;
73  q2.setIdentity();
74  VERIFY_IS_APPROX(Quaternionx(Quaternionx::Identity()).coeffs(), q2.coeffs());
75  q1.coeffs().setRandom();
76  VERIFY_IS_APPROX(q1.coeffs(), (q1*q2).coeffs());
77 
78 #ifndef EIGEN_NO_IO
79  // Printing
80  std::ostringstream ss;
81  ss << q2;
82  VERIFY(ss.str() == "0i + 0j + 0k + 1");
83 #endif
84 
85  // concatenation
86  q1 *= q2;
87 
88  q1 = AngleAxisx(a, v0.normalized());
89  q2 = AngleAxisx(a, v1.normalized());
90 
91  // angular distance
92  Scalar refangle = abs(AngleAxisx(q1.inverse()*q2).angle());
93  if (refangle>Scalar(EIGEN_PI))
94  refangle = Scalar(2)*Scalar(EIGEN_PI) - refangle;
95 
96  if((q1.coeffs()-q2.coeffs()).norm() > Scalar(10)*largeEps)
97  {
98  VERIFY_IS_MUCH_SMALLER_THAN(abs(q1.angularDistance(q2) - refangle), Scalar(1));
99  }
100 
101  // rotation matrix conversion
102  VERIFY_IS_APPROX(q1 * v2, q1.toRotationMatrix() * v2);
103  VERIFY_IS_APPROX(q1 * q2 * v2,
104  q1.toRotationMatrix() * q2.toRotationMatrix() * v2);
105 
106  VERIFY( (q2*q1).isApprox(q1*q2, largeEps)
107  || !(q2 * q1 * v2).isApprox(q1.toRotationMatrix() * q2.toRotationMatrix() * v2));
108 
109  q2 = q1.toRotationMatrix();
110  VERIFY_IS_APPROX(q1*v1,q2*v1);
111 
112  Matrix3 rot1(q1);
113  VERIFY_IS_APPROX(q1*v1,rot1*v1);
114  Quaternionx q3(rot1.transpose()*rot1);
115  VERIFY_IS_APPROX(q3*v1,v1);
116 
117 
118  // angle-axis conversion
119  AngleAxisx aa = AngleAxisx(q1);
120  VERIFY_IS_APPROX(q1 * v1, Quaternionx(aa) * v1);
121 
122  // Do not execute the test if the rotation angle is almost zero, or
123  // the rotation axis and v1 are almost parallel.
124  if (abs(aa.angle()) > Scalar(5)*test_precision<Scalar>()
125  && (aa.axis() - v1.normalized()).norm() < Scalar(1.99)
126  && (aa.axis() + v1.normalized()).norm() < Scalar(1.99))
127  {
128  VERIFY_IS_NOT_APPROX(q1 * v1, Quaternionx(AngleAxisx(aa.angle()*2,aa.axis())) * v1);
129  }
130 
131  // from two vector creation
132  VERIFY_IS_APPROX( v2.normalized(),(q2.setFromTwoVectors(v1, v2)*v1).normalized());
133  VERIFY_IS_APPROX( v1.normalized(),(q2.setFromTwoVectors(v1, v1)*v1).normalized());
134  VERIFY_IS_APPROX(-v1.normalized(),(q2.setFromTwoVectors(v1,-v1)*v1).normalized());
135  if (internal::is_same<Scalar,double>::value)
136  {
137  v3 = (v1.array()+eps).matrix();
138  VERIFY_IS_APPROX( v3.normalized(),(q2.setFromTwoVectors(v1, v3)*v1).normalized());
139  VERIFY_IS_APPROX(-v3.normalized(),(q2.setFromTwoVectors(v1,-v3)*v1).normalized());
140  }
141 
142  // from two vector creation static function
143  VERIFY_IS_APPROX( v2.normalized(),(Quaternionx::FromTwoVectors(v1, v2)*v1).normalized());
144  VERIFY_IS_APPROX( v1.normalized(),(Quaternionx::FromTwoVectors(v1, v1)*v1).normalized());
145  VERIFY_IS_APPROX(-v1.normalized(),(Quaternionx::FromTwoVectors(v1,-v1)*v1).normalized());
146  if (internal::is_same<Scalar,double>::value)
147  {
148  v3 = (v1.array()+eps).matrix();
149  VERIFY_IS_APPROX( v3.normalized(),(Quaternionx::FromTwoVectors(v1, v3)*v1).normalized());
150  VERIFY_IS_APPROX(-v3.normalized(),(Quaternionx::FromTwoVectors(v1,-v3)*v1).normalized());
151  }
152 
153  // inverse and conjugate
154  VERIFY_IS_APPROX(q1 * (q1.inverse() * v1), v1);
155  VERIFY_IS_APPROX(q1 * (q1.conjugate() * v1), v1);
156 
157  // test casting
158  Quaternion<float> q1f = q1.template cast<float>();
159  VERIFY_IS_APPROX(q1f.template cast<Scalar>(),q1);
160  Quaternion<double> q1d = q1.template cast<double>();
161  VERIFY_IS_APPROX(q1d.template cast<Scalar>(),q1);
162 
163  // test bug 369 - improper alignment.
164  Quaternionx *q = new Quaternionx;
165  delete q;
166 
167  q1 = Quaternionx::UnitRandom();
168  q2 = Quaternionx::UnitRandom();
169  check_slerp(q1,q2);
170 
171  q1 = AngleAxisx(b, v1.normalized());
172  q2 = AngleAxisx(b+Scalar(EIGEN_PI), v1.normalized());
173  check_slerp(q1,q2);
174 
175  q1 = AngleAxisx(b, v1.normalized());
176  q2 = AngleAxisx(-b, -v1.normalized());
177  check_slerp(q1,q2);
178 
179  q1 = Quaternionx::UnitRandom();
180  q2.coeffs() = -q1.coeffs();
181  check_slerp(q1,q2);
182 }
183 
184 template<typename Scalar> void mapQuaternion(void){
185  typedef Map<Quaternion<Scalar>, Aligned> MQuaternionA;
186  typedef Map<const Quaternion<Scalar>, Aligned> MCQuaternionA;
187  typedef Map<Quaternion<Scalar> > MQuaternionUA;
188  typedef Map<const Quaternion<Scalar> > MCQuaternionUA;
189  typedef Quaternion<Scalar> Quaternionx;
190  typedef Matrix<Scalar,3,1> Vector3;
191  typedef AngleAxis<Scalar> AngleAxisx;
192 
193  Vector3 v0 = Vector3::Random(),
194  v1 = Vector3::Random();
195  Scalar a = internal::random<Scalar>(-Scalar(EIGEN_PI), Scalar(EIGEN_PI));
196 
197  EIGEN_ALIGN_MAX Scalar array1[4];
198  EIGEN_ALIGN_MAX Scalar array2[4];
199  EIGEN_ALIGN_MAX Scalar array3[4+1];
200  Scalar* array3unaligned = array3+1;
201 
202  MQuaternionA mq1(array1);
203  MCQuaternionA mcq1(array1);
204  MQuaternionA mq2(array2);
205  MQuaternionUA mq3(array3unaligned);
206  MCQuaternionUA mcq3(array3unaligned);
207 
208 // std::cerr << array1 << " " << array2 << " " << array3 << "\n";
209  mq1 = AngleAxisx(a, v0.normalized());
210  mq2 = mq1;
211  mq3 = mq1;
212 
213  Quaternionx q1 = mq1;
214  Quaternionx q2 = mq2;
215  Quaternionx q3 = mq3;
216  Quaternionx q4 = MCQuaternionUA(array3unaligned);
217 
218  VERIFY_IS_APPROX(q1.coeffs(), q2.coeffs());
219  VERIFY_IS_APPROX(q1.coeffs(), q3.coeffs());
220  VERIFY_IS_APPROX(q4.coeffs(), q3.coeffs());
221 
222  VERIFY_IS_APPROX(mq1 * (mq1.inverse() * v1), v1);
223  VERIFY_IS_APPROX(mq1 * (mq1.conjugate() * v1), v1);
224 
225  VERIFY_IS_APPROX(mcq1 * (mcq1.inverse() * v1), v1);
226  VERIFY_IS_APPROX(mcq1 * (mcq1.conjugate() * v1), v1);
227 
228  VERIFY_IS_APPROX(mq3 * (mq3.inverse() * v1), v1);
229  VERIFY_IS_APPROX(mq3 * (mq3.conjugate() * v1), v1);
230 
231  VERIFY_IS_APPROX(mcq3 * (mcq3.inverse() * v1), v1);
232  VERIFY_IS_APPROX(mcq3 * (mcq3.conjugate() * v1), v1);
233 
234  VERIFY_IS_APPROX(mq1*mq2, q1*q2);
235  VERIFY_IS_APPROX(mq3*mq2, q3*q2);
236  VERIFY_IS_APPROX(mcq1*mq2, q1*q2);
237  VERIFY_IS_APPROX(mcq3*mq2, q3*q2);
238 
239  // Bug 1461, compilation issue with Map<const Quat>::w(), and other reference/constness checks:
240  VERIFY_IS_APPROX(mcq3.coeffs().x() + mcq3.coeffs().y() + mcq3.coeffs().z() + mcq3.coeffs().w(), mcq3.coeffs().sum());
241  VERIFY_IS_APPROX(mcq3.x() + mcq3.y() + mcq3.z() + mcq3.w(), mcq3.coeffs().sum());
242  mq3.w() = 1;
243  const Quaternionx& cq3(q3);
244  VERIFY( &cq3.x() == &q3.x() );
245  const MQuaternionUA& cmq3(mq3);
246  VERIFY( &cmq3.x() == &mq3.x() );
247  // FIXME the following should be ok. The problem is that currently the LValueBit flag
248  // is used to determine whether we can return a coeff by reference or not, which is not enough for Map<const ...>.
249  //const MCQuaternionUA& cmcq3(mcq3);
250  //VERIFY( &cmcq3.x() == &mcq3.x() );
251 
252  // test cast
253  {
254  Quaternion<float> q1f = mq1.template cast<float>();
255  VERIFY_IS_APPROX(q1f.template cast<Scalar>(),mq1);
256  Quaternion<double> q1d = mq1.template cast<double>();
257  VERIFY_IS_APPROX(q1d.template cast<Scalar>(),mq1);
258  }
259 }
260 
261 template<typename Scalar> void quaternionAlignment(void){
262  typedef Quaternion<Scalar,AutoAlign> QuaternionA;
263  typedef Quaternion<Scalar,DontAlign> QuaternionUA;
264 
265  EIGEN_ALIGN_MAX Scalar array1[4];
266  EIGEN_ALIGN_MAX Scalar array2[4];
267  EIGEN_ALIGN_MAX Scalar array3[4+1];
268  Scalar* arrayunaligned = array3+1;
269 
270  QuaternionA *q1 = ::new(reinterpret_cast<void*>(array1)) QuaternionA;
271  QuaternionUA *q2 = ::new(reinterpret_cast<void*>(array2)) QuaternionUA;
272  QuaternionUA *q3 = ::new(reinterpret_cast<void*>(arrayunaligned)) QuaternionUA;
273 
274  q1->coeffs().setRandom();
275  *q2 = *q1;
276  *q3 = *q1;
277 
278  VERIFY_IS_APPROX(q1->coeffs(), q2->coeffs());
279  VERIFY_IS_APPROX(q1->coeffs(), q3->coeffs());
280 }
281 
282 template<typename PlainObjectType> void check_const_correctness(const PlainObjectType&)
283 {
284  // there's a lot that we can't test here while still having this test compile!
285  // the only possible approach would be to run a script trying to compile stuff and checking that it fails.
286  // CMake can help with that.
287 
288  // verify that map-to-const don't have LvalueBit
289  typedef typename internal::add_const<PlainObjectType>::type ConstPlainObjectType;
290  VERIFY( !(internal::traits<Map<ConstPlainObjectType> >::Flags & LvalueBit) );
291  VERIFY( !(internal::traits<Map<ConstPlainObjectType, Aligned> >::Flags & LvalueBit) );
292  VERIFY( !(Map<ConstPlainObjectType>::Flags & LvalueBit) );
293  VERIFY( !(Map<ConstPlainObjectType, Aligned>::Flags & LvalueBit) );
294 }
295 
296 #if EIGEN_HAS_RVALUE_REFERENCES
297 
298 // Regression for bug 1573
299 struct MovableClass {
300  // The following line is a workaround for gcc 4.7 and 4.8 (see bug 1573 comments).
301  static_assert(std::is_nothrow_move_constructible<Quaternionf>::value,"");
302  MovableClass() = default;
303  MovableClass(const MovableClass&) = default;
304  MovableClass(MovableClass&&) noexcept = default;
305  MovableClass& operator=(const MovableClass&) = default;
306  MovableClass& operator=(MovableClass&&) = default;
307  Quaternionf m_quat;
308 };
309 
310 #endif
311 
312 EIGEN_DECLARE_TEST(geo_quaternion)
313 {
314  for(int i = 0; i < g_repeat; i++) {
315  CALL_SUBTEST_1(( quaternion<float,AutoAlign>() ));
316  CALL_SUBTEST_1( check_const_correctness(Quaternionf()) );
317  CALL_SUBTEST_1(( quaternion<float,DontAlign>() ));
318  CALL_SUBTEST_1(( quaternionAlignment<float>() ));
319  CALL_SUBTEST_1( mapQuaternion<float>() );
320 
321  CALL_SUBTEST_2(( quaternion<double,AutoAlign>() ));
322  CALL_SUBTEST_2( check_const_correctness(Quaterniond()) );
323  CALL_SUBTEST_2(( quaternion<double,DontAlign>() ));
324  CALL_SUBTEST_2(( quaternionAlignment<double>() ));
325  CALL_SUBTEST_2( mapQuaternion<double>() );
326 
327 #ifndef EIGEN_TEST_ANNOYING_SCALAR_DONT_THROW
328  AnnoyingScalar::dont_throw = true;
329 #endif
330  CALL_SUBTEST_3(( quaternion<AnnoyingScalar,AutoAlign>() ));
331  }
332 }
Scalar
SCALAR Scalar
Definition: bench_gemm.cpp:46
max
#define max(a, b)
Definition: datatypes.h:20
v2
Vector v2
Definition: testSerializationBase.cpp:39
b
Scalar * b
Definition: benchVecAdd.cpp:17
gtsam::Vector3
Eigen::Vector3d Vector3
Definition: Vector.h:43
v1
Vector v1
Definition: testSerializationBase.cpp:38
EIGEN_PI
#define EIGEN_PI
Definition: Eigen/src/Core/MathFunctions.h:16
min
#define min(a, b)
Definition: datatypes.h:19
Eigen::Map
A matrix or vector expression mapping an existing array of data.
Definition: Map.h:94
CALL_SUBTEST_3
#define CALL_SUBTEST_3(FUNC)
Definition: split_test_helper.h:16
Eigen::Transpose
Expression of the transpose of a matrix.
Definition: Transpose.h:52
Eigen::Quaterniond
Quaternion< double > Quaterniond
Definition: 3rdparty/Eigen/Eigen/src/Geometry/Quaternion.h:366
Eigen::LvalueBit
const unsigned int LvalueBit
Definition: Constants.h:144
ceres::acos
Jet< T, N > acos(const Jet< T, N > &f)
Definition: jet.h:432
EIGEN_ALIGN_MAX
#define EIGEN_ALIGN_MAX
Definition: ConfigureVectorization.h:157
quaternion
void quaternion(void)
Definition: geo_quaternion.cpp:46
AnnoyingScalar::dont_throw
static bool dont_throw
Definition: AnnoyingScalar.h:98
bounded_acos
T bounded_acos(T v)
Definition: geo_quaternion.cpp:17
VERIFY_IS_APPROX
#define VERIFY_IS_APPROX(a, b)
Definition: integer_types.cpp:15
v3
Vector v3
Definition: testSerializationBase.cpp:40
CALL_SUBTEST_1
#define CALL_SUBTEST_1(FUNC)
Definition: split_test_helper.h:4
test_callbacks.value
value
Definition: test_callbacks.py:148
Eigen::g_repeat
static int g_repeat
Definition: main.h:169
v
Array< int, Dynamic, 1 > v
Definition: Array_initializer_list_vector_cxx11.cpp:1
T
Eigen::Triplet< double > T
Definition: Tutorial_sparse_example.cpp:6
Eigen::Aligned
Definition: Constants.h:240
e
Array< double, 1, 3 > e(1./3., 0.5, 2.)
Eigen::numext::q
EIGEN_DEVICE_FUNC const Scalar & q
Definition: SpecialFunctionsImpl.h:1984
ss
static std::stringstream ss
Definition: testBTree.cpp:31
VERIFY_IS_MUCH_SMALLER_THAN
#define VERIFY_IS_MUCH_SMALLER_THAN(a, b)
Definition: main.h:390
mapQuaternion
void mapQuaternion(void)
Definition: geo_quaternion.cpp:184
VERIFY
#define VERIFY(a)
Definition: main.h:380
Eigen::Quaternion
The quaternion class used to represent 3D orientations and rotations.
Definition: ForwardDeclarations.h:293
v0
static const double v0
Definition: testCal3DFisheye.cpp:31
VERIFY_IS_NOT_APPROX
#define VERIFY_IS_NOT_APPROX(a, b)
Definition: integer_types.cpp:17
check_const_correctness
void check_const_correctness(const PlainObjectType &)
Definition: geo_quaternion.cpp:282
a
ArrayXXi a
Definition: Array_initializer_list_23_cxx11.cpp:1
CALL_SUBTEST_2
#define CALL_SUBTEST_2(FUNC)
Definition: split_test_helper.h:10
quaternionAlignment
void quaternionAlignment(void)
Definition: geo_quaternion.cpp:261
matrix
Map< Matrix< T, Dynamic, Dynamic, ColMajor >, 0, OuterStride<> > matrix(T *data, int rows, int cols, int stride)
Definition: gtsam/3rdparty/Eigen/blas/common.h:110
Eigen::internal::isApprox
EIGEN_DEVICE_FUNC bool isApprox(const Scalar &x, const Scalar &y, const typename NumTraits< Scalar >::Real &precision=NumTraits< Scalar >::dummy_precision())
Definition: Eigen/src/Core/MathFunctions.h:1947
abs
#define abs(x)
Definition: datatypes.h:17
i
int i
Definition: BiCGSTAB_step_by_step.cpp:9
check_slerp
void check_slerp(const QuatType &q0, const QuatType &q1)
Definition: geo_quaternion.cpp:25
Eigen::AngleAxis
Represents a 3D rotation as a rotation angle around an arbitrary 3D axis.
Definition: ForwardDeclarations.h:290
align_3::t
Point2 t(10, 10)
Eigen::Quaternionf
Quaternion< float > Quaternionf
Definition: 3rdparty/Eigen/Eigen/src/Geometry/Quaternion.h:363
EIGEN_DECLARE_TEST
EIGEN_DECLARE_TEST(geo_quaternion)
Definition: geo_quaternion.cpp:312

gtsam
Author(s):
autogenerated on Tue Jul 4 2023 02:34:17