1 // This file is part of Eigen, a lightweight C++ template library
2 // for linear algebra.
3 //
4 // Copyright (C) 2008-2009 Gael Guennebaud <[email protected]>
5 // Copyright (C) 2009 Mathieu Gautier <[email protected]>
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
bounded_acos(T v)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
check_slerp(const QuatType & q0,const QuatType & q1)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
quaternion(void)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
mapQuaternion(void)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
quaternionAlignment(void)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
check_const_correctness(const PlainObjectType &)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
EIGEN_DECLARE_TEST(geo_quaternion)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 }
333