geo_quaternion.cpp revision c981c48f5bc9aefeffc0bcb0cc3934c2fae179dd
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 16template<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 24template<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 45template<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 172template<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 202template<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 227template<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 241void 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