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
bounded_acos(T v)16 template<typename T> T bounded_acos(T v)
17 {
18 using std::acos;
19 using std::min;
20 using std::max;
21 return acos((max)(T(-1),(min)(v,T(1))));
22 }
23
check_slerp(const QuatType & q0,const QuatType & q1)24 template<typename QuatType> void check_slerp(const QuatType& q0, const QuatType& q1)
25 {
26 typedef typename QuatType::Scalar Scalar;
27 typedef Matrix<Scalar,3,1> VectorType;
28 typedef AngleAxis<Scalar> AA;
29
30 Scalar largeEps = test_precision<Scalar>();
31
32 Scalar theta_tot = AA(q1*q0.inverse()).angle();
33 if(theta_tot>M_PI)
34 theta_tot = 2.*M_PI-theta_tot;
35 for(Scalar t=0; t<=1.001; t+=0.1)
36 {
37 QuatType q = q0.slerp(t,q1);
38 Scalar theta = AA(q*q0.inverse()).angle();
39 VERIFY(internal::abs(q.norm() - 1) < largeEps);
40 if(theta_tot==0) VERIFY(theta_tot==0);
41 else VERIFY(internal::abs(theta/theta_tot - t) < largeEps);
42 }
43 }
44
quaternion(void)45 template<typename Scalar, int Options> void quaternion(void)
46 {
47 /* this test covers the following files:
48 Quaternion.h
49 */
50
51 typedef Matrix<Scalar,3,3> Matrix3;
52 typedef Matrix<Scalar,3,1> Vector3;
53 typedef Matrix<Scalar,4,1> Vector4;
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 = 1e-3f;
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(M_PI), Scalar(M_PI)),
69 b = internal::random<Scalar>(-Scalar(M_PI), Scalar(M_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 // concatenation
79 q1 *= q2;
80
81 q1 = AngleAxisx(a, v0.normalized());
82 q2 = AngleAxisx(a, v1.normalized());
83
84 // angular distance
85 Scalar refangle = internal::abs(AngleAxisx(q1.inverse()*q2).angle());
86 if (refangle>Scalar(M_PI))
87 refangle = Scalar(2)*Scalar(M_PI) - refangle;
88
89 if((q1.coeffs()-q2.coeffs()).norm() > 10*largeEps)
90 {
91 VERIFY_IS_MUCH_SMALLER_THAN(internal::abs(q1.angularDistance(q2) - refangle), Scalar(1));
92 }
93
94 // rotation matrix conversion
95 VERIFY_IS_APPROX(q1 * v2, q1.toRotationMatrix() * v2);
96 VERIFY_IS_APPROX(q1 * q2 * v2,
97 q1.toRotationMatrix() * q2.toRotationMatrix() * v2);
98
99 VERIFY( (q2*q1).isApprox(q1*q2, largeEps)
100 || !(q2 * q1 * v2).isApprox(q1.toRotationMatrix() * q2.toRotationMatrix() * v2));
101
102 q2 = q1.toRotationMatrix();
103 VERIFY_IS_APPROX(q1*v1,q2*v1);
104
105
106 // angle-axis conversion
107 AngleAxisx aa = AngleAxisx(q1);
108 VERIFY_IS_APPROX(q1 * v1, Quaternionx(aa) * v1);
109
110 // Do not execute the test if the rotation angle is almost zero, or
111 // the rotation axis and v1 are almost parallel.
112 if (internal::abs(aa.angle()) > 5*test_precision<Scalar>()
113 && (aa.axis() - v1.normalized()).norm() < 1.99
114 && (aa.axis() + v1.normalized()).norm() < 1.99)
115 {
116 VERIFY_IS_NOT_APPROX(q1 * v1, Quaternionx(AngleAxisx(aa.angle()*2,aa.axis())) * v1);
117 }
118
119 // from two vector creation
120 VERIFY_IS_APPROX( v2.normalized(),(q2.setFromTwoVectors(v1, v2)*v1).normalized());
121 VERIFY_IS_APPROX( v1.normalized(),(q2.setFromTwoVectors(v1, v1)*v1).normalized());
122 VERIFY_IS_APPROX(-v1.normalized(),(q2.setFromTwoVectors(v1,-v1)*v1).normalized());
123 if (internal::is_same<Scalar,double>::value)
124 {
125 v3 = (v1.array()+eps).matrix();
126 VERIFY_IS_APPROX( v3.normalized(),(q2.setFromTwoVectors(v1, v3)*v1).normalized());
127 VERIFY_IS_APPROX(-v3.normalized(),(q2.setFromTwoVectors(v1,-v3)*v1).normalized());
128 }
129
130 // from two vector creation static function
131 VERIFY_IS_APPROX( v2.normalized(),(Quaternionx::FromTwoVectors(v1, v2)*v1).normalized());
132 VERIFY_IS_APPROX( v1.normalized(),(Quaternionx::FromTwoVectors(v1, v1)*v1).normalized());
133 VERIFY_IS_APPROX(-v1.normalized(),(Quaternionx::FromTwoVectors(v1,-v1)*v1).normalized());
134 if (internal::is_same<Scalar,double>::value)
135 {
136 v3 = (v1.array()+eps).matrix();
137 VERIFY_IS_APPROX( v3.normalized(),(Quaternionx::FromTwoVectors(v1, v3)*v1).normalized());
138 VERIFY_IS_APPROX(-v3.normalized(),(Quaternionx::FromTwoVectors(v1,-v3)*v1).normalized());
139 }
140
141 // inverse and conjugate
142 VERIFY_IS_APPROX(q1 * (q1.inverse() * v1), v1);
143 VERIFY_IS_APPROX(q1 * (q1.conjugate() * v1), v1);
144
145 // test casting
146 Quaternion<float> q1f = q1.template cast<float>();
147 VERIFY_IS_APPROX(q1f.template cast<Scalar>(),q1);
148 Quaternion<double> q1d = q1.template cast<double>();
149 VERIFY_IS_APPROX(q1d.template cast<Scalar>(),q1);
150
151 // test bug 369 - improper alignment.
152 Quaternionx *q = new Quaternionx;
153 delete q;
154
155 q1 = AngleAxisx(a, v0.normalized());
156 q2 = AngleAxisx(b, v1.normalized());
157 check_slerp(q1,q2);
158
159 q1 = AngleAxisx(b, v1.normalized());
160 q2 = AngleAxisx(b+M_PI, v1.normalized());
161 check_slerp(q1,q2);
162
163 q1 = AngleAxisx(b, v1.normalized());
164 q2 = AngleAxisx(-b, -v1.normalized());
165 check_slerp(q1,q2);
166
167 q1.coeffs() = Vector4::Random().normalized();
168 q2.coeffs() = -q1.coeffs();
169 check_slerp(q1,q2);
170 }
171
mapQuaternion(void)172 template<typename Scalar> void mapQuaternion(void){
173 typedef Map<Quaternion<Scalar>, Aligned> MQuaternionA;
174 typedef Map<Quaternion<Scalar> > MQuaternionUA;
175 typedef Map<const Quaternion<Scalar> > MCQuaternionUA;
176 typedef Quaternion<Scalar> Quaternionx;
177
178 EIGEN_ALIGN16 Scalar array1[4];
179 EIGEN_ALIGN16 Scalar array2[4];
180 EIGEN_ALIGN16 Scalar array3[4+1];
181 Scalar* array3unaligned = array3+1;
182
183 // std::cerr << array1 << " " << array2 << " " << array3 << "\n";
184 MQuaternionA(array1).coeffs().setRandom();
185 (MQuaternionA(array2)) = MQuaternionA(array1);
186 (MQuaternionUA(array3unaligned)) = MQuaternionA(array1);
187
188 Quaternionx q1 = MQuaternionA(array1);
189 Quaternionx q2 = MQuaternionA(array2);
190 Quaternionx q3 = MQuaternionUA(array3unaligned);
191 Quaternionx q4 = MCQuaternionUA(array3unaligned);
192
193 VERIFY_IS_APPROX(q1.coeffs(), q2.coeffs());
194 VERIFY_IS_APPROX(q1.coeffs(), q3.coeffs());
195 VERIFY_IS_APPROX(q4.coeffs(), q3.coeffs());
196 #ifdef EIGEN_VECTORIZE
197 if(internal::packet_traits<Scalar>::Vectorizable)
198 VERIFY_RAISES_ASSERT((MQuaternionA(array3unaligned)));
199 #endif
200 }
201
quaternionAlignment(void)202 template<typename Scalar> void quaternionAlignment(void){
203 typedef Quaternion<Scalar,AutoAlign> QuaternionA;
204 typedef Quaternion<Scalar,DontAlign> QuaternionUA;
205
206 EIGEN_ALIGN16 Scalar array1[4];
207 EIGEN_ALIGN16 Scalar array2[4];
208 EIGEN_ALIGN16 Scalar array3[4+1];
209 Scalar* arrayunaligned = array3+1;
210
211 QuaternionA *q1 = ::new(reinterpret_cast<void*>(array1)) QuaternionA;
212 QuaternionUA *q2 = ::new(reinterpret_cast<void*>(array2)) QuaternionUA;
213 QuaternionUA *q3 = ::new(reinterpret_cast<void*>(arrayunaligned)) QuaternionUA;
214
215 q1->coeffs().setRandom();
216 *q2 = *q1;
217 *q3 = *q1;
218
219 VERIFY_IS_APPROX(q1->coeffs(), q2->coeffs());
220 VERIFY_IS_APPROX(q1->coeffs(), q3->coeffs());
221 #if defined(EIGEN_VECTORIZE) && EIGEN_ALIGN_STATICALLY
222 if(internal::packet_traits<Scalar>::Vectorizable)
223 VERIFY_RAISES_ASSERT((::new(reinterpret_cast<void*>(arrayunaligned)) QuaternionA));
224 #endif
225 }
226
check_const_correctness(const PlainObjectType &)227 template<typename PlainObjectType> void check_const_correctness(const PlainObjectType&)
228 {
229 // there's a lot that we can't test here while still having this test compile!
230 // the only possible approach would be to run a script trying to compile stuff and checking that it fails.
231 // CMake can help with that.
232
233 // verify that map-to-const don't have LvalueBit
234 typedef typename internal::add_const<PlainObjectType>::type ConstPlainObjectType;
235 VERIFY( !(internal::traits<Map<ConstPlainObjectType> >::Flags & LvalueBit) );
236 VERIFY( !(internal::traits<Map<ConstPlainObjectType, Aligned> >::Flags & LvalueBit) );
237 VERIFY( !(Map<ConstPlainObjectType>::Flags & LvalueBit) );
238 VERIFY( !(Map<ConstPlainObjectType, Aligned>::Flags & LvalueBit) );
239 }
240
test_geo_quaternion()241 void test_geo_quaternion()
242 {
243 for(int i = 0; i < g_repeat; i++) {
244 CALL_SUBTEST_1(( quaternion<float,AutoAlign>() ));
245 CALL_SUBTEST_1( check_const_correctness(Quaternionf()) );
246 CALL_SUBTEST_2(( quaternion<double,AutoAlign>() ));
247 CALL_SUBTEST_2( check_const_correctness(Quaterniond()) );
248 CALL_SUBTEST_3(( quaternion<float,DontAlign>() ));
249 CALL_SUBTEST_4(( quaternion<double,DontAlign>() ));
250 CALL_SUBTEST_5(( quaternionAlignment<float>() ));
251 CALL_SUBTEST_6(( quaternionAlignment<double>() ));
252 CALL_SUBTEST_1( mapQuaternion<float>() );
253 CALL_SUBTEST_2( mapQuaternion<double>() );
254 }
255 }
256