geo_quaternion.cpp (11568B)
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 }