geo_quaternion.cpp 11 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332
  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. #include "main.h"
  11. #include <Eigen/Geometry>
  12. #include <Eigen/LU>
  13. #include <Eigen/SVD>
  14. #include "AnnoyingScalar.h"
  15. template<typename T> T bounded_acos(T v)
  16. {
  17. using std::acos;
  18. using std::min;
  19. using std::max;
  20. return acos((max)(T(-1),(min)(v,T(1))));
  21. }
  22. template<typename QuatType> void check_slerp(const QuatType& q0, const QuatType& q1)
  23. {
  24. using std::abs;
  25. typedef typename QuatType::Scalar Scalar;
  26. typedef AngleAxis<Scalar> AA;
  27. Scalar largeEps = test_precision<Scalar>();
  28. Scalar theta_tot = AA(q1*q0.inverse()).angle();
  29. if(theta_tot>Scalar(EIGEN_PI))
  30. theta_tot = Scalar(2.)*Scalar(EIGEN_PI)-theta_tot;
  31. for(Scalar t=0; t<=Scalar(1.001); t+=Scalar(0.1))
  32. {
  33. QuatType q = q0.slerp(t,q1);
  34. Scalar theta = AA(q*q0.inverse()).angle();
  35. VERIFY(abs(q.norm() - 1) < largeEps);
  36. if(theta_tot==0) VERIFY(theta_tot==0);
  37. else VERIFY(abs(theta - t * theta_tot) < largeEps);
  38. }
  39. }
  40. template<typename Scalar, int Options> void quaternion(void)
  41. {
  42. /* this test covers the following files:
  43. Quaternion.h
  44. */
  45. using std::abs;
  46. typedef Matrix<Scalar,3,1> Vector3;
  47. typedef Matrix<Scalar,3,3> Matrix3;
  48. typedef Quaternion<Scalar,Options> Quaternionx;
  49. typedef AngleAxis<Scalar> AngleAxisx;
  50. Scalar largeEps = test_precision<Scalar>();
  51. if (internal::is_same<Scalar,float>::value)
  52. largeEps = Scalar(1e-3);
  53. Scalar eps = internal::random<Scalar>() * Scalar(1e-2);
  54. Vector3 v0 = Vector3::Random(),
  55. v1 = Vector3::Random(),
  56. v2 = Vector3::Random(),
  57. v3 = Vector3::Random();
  58. Scalar a = internal::random<Scalar>(-Scalar(EIGEN_PI), Scalar(EIGEN_PI)),
  59. b = internal::random<Scalar>(-Scalar(EIGEN_PI), Scalar(EIGEN_PI));
  60. // Quaternion: Identity(), setIdentity();
  61. Quaternionx q1, q2;
  62. q2.setIdentity();
  63. VERIFY_IS_APPROX(Quaternionx(Quaternionx::Identity()).coeffs(), q2.coeffs());
  64. q1.coeffs().setRandom();
  65. VERIFY_IS_APPROX(q1.coeffs(), (q1*q2).coeffs());
  66. #ifndef EIGEN_NO_IO
  67. // Printing
  68. std::ostringstream ss;
  69. ss << q2;
  70. VERIFY(ss.str() == "0i + 0j + 0k + 1");
  71. #endif
  72. // concatenation
  73. q1 *= q2;
  74. q1 = AngleAxisx(a, v0.normalized());
  75. q2 = AngleAxisx(a, v1.normalized());
  76. // angular distance
  77. Scalar refangle = abs(AngleAxisx(q1.inverse()*q2).angle());
  78. if (refangle>Scalar(EIGEN_PI))
  79. refangle = Scalar(2)*Scalar(EIGEN_PI) - refangle;
  80. if((q1.coeffs()-q2.coeffs()).norm() > Scalar(10)*largeEps)
  81. {
  82. VERIFY_IS_MUCH_SMALLER_THAN(abs(q1.angularDistance(q2) - refangle), Scalar(1));
  83. }
  84. // rotation matrix conversion
  85. VERIFY_IS_APPROX(q1 * v2, q1.toRotationMatrix() * v2);
  86. VERIFY_IS_APPROX(q1 * q2 * v2,
  87. q1.toRotationMatrix() * q2.toRotationMatrix() * v2);
  88. VERIFY( (q2*q1).isApprox(q1*q2, largeEps)
  89. || !(q2 * q1 * v2).isApprox(q1.toRotationMatrix() * q2.toRotationMatrix() * v2));
  90. q2 = q1.toRotationMatrix();
  91. VERIFY_IS_APPROX(q1*v1,q2*v1);
  92. Matrix3 rot1(q1);
  93. VERIFY_IS_APPROX(q1*v1,rot1*v1);
  94. Quaternionx q3(rot1.transpose()*rot1);
  95. VERIFY_IS_APPROX(q3*v1,v1);
  96. // angle-axis conversion
  97. AngleAxisx aa = AngleAxisx(q1);
  98. VERIFY_IS_APPROX(q1 * v1, Quaternionx(aa) * v1);
  99. // Do not execute the test if the rotation angle is almost zero, or
  100. // the rotation axis and v1 are almost parallel.
  101. if (abs(aa.angle()) > Scalar(5)*test_precision<Scalar>()
  102. && (aa.axis() - v1.normalized()).norm() < Scalar(1.99)
  103. && (aa.axis() + v1.normalized()).norm() < Scalar(1.99))
  104. {
  105. VERIFY_IS_NOT_APPROX(q1 * v1, Quaternionx(AngleAxisx(aa.angle()*2,aa.axis())) * v1);
  106. }
  107. // from two vector creation
  108. VERIFY_IS_APPROX( v2.normalized(),(q2.setFromTwoVectors(v1, v2)*v1).normalized());
  109. VERIFY_IS_APPROX( v1.normalized(),(q2.setFromTwoVectors(v1, v1)*v1).normalized());
  110. VERIFY_IS_APPROX(-v1.normalized(),(q2.setFromTwoVectors(v1,-v1)*v1).normalized());
  111. if (internal::is_same<Scalar,double>::value)
  112. {
  113. v3 = (v1.array()+eps).matrix();
  114. VERIFY_IS_APPROX( v3.normalized(),(q2.setFromTwoVectors(v1, v3)*v1).normalized());
  115. VERIFY_IS_APPROX(-v3.normalized(),(q2.setFromTwoVectors(v1,-v3)*v1).normalized());
  116. }
  117. // from two vector creation static function
  118. VERIFY_IS_APPROX( v2.normalized(),(Quaternionx::FromTwoVectors(v1, v2)*v1).normalized());
  119. VERIFY_IS_APPROX( v1.normalized(),(Quaternionx::FromTwoVectors(v1, v1)*v1).normalized());
  120. VERIFY_IS_APPROX(-v1.normalized(),(Quaternionx::FromTwoVectors(v1,-v1)*v1).normalized());
  121. if (internal::is_same<Scalar,double>::value)
  122. {
  123. v3 = (v1.array()+eps).matrix();
  124. VERIFY_IS_APPROX( v3.normalized(),(Quaternionx::FromTwoVectors(v1, v3)*v1).normalized());
  125. VERIFY_IS_APPROX(-v3.normalized(),(Quaternionx::FromTwoVectors(v1,-v3)*v1).normalized());
  126. }
  127. // inverse and conjugate
  128. VERIFY_IS_APPROX(q1 * (q1.inverse() * v1), v1);
  129. VERIFY_IS_APPROX(q1 * (q1.conjugate() * v1), v1);
  130. // test casting
  131. Quaternion<float> q1f = q1.template cast<float>();
  132. VERIFY_IS_APPROX(q1f.template cast<Scalar>(),q1);
  133. Quaternion<double> q1d = q1.template cast<double>();
  134. VERIFY_IS_APPROX(q1d.template cast<Scalar>(),q1);
  135. // test bug 369 - improper alignment.
  136. Quaternionx *q = new Quaternionx;
  137. delete q;
  138. q1 = Quaternionx::UnitRandom();
  139. q2 = Quaternionx::UnitRandom();
  140. check_slerp(q1,q2);
  141. q1 = AngleAxisx(b, v1.normalized());
  142. q2 = AngleAxisx(b+Scalar(EIGEN_PI), v1.normalized());
  143. check_slerp(q1,q2);
  144. q1 = AngleAxisx(b, v1.normalized());
  145. q2 = AngleAxisx(-b, -v1.normalized());
  146. check_slerp(q1,q2);
  147. q1 = Quaternionx::UnitRandom();
  148. q2.coeffs() = -q1.coeffs();
  149. check_slerp(q1,q2);
  150. }
  151. template<typename Scalar> void mapQuaternion(void){
  152. typedef Map<Quaternion<Scalar>, Aligned> MQuaternionA;
  153. typedef Map<const Quaternion<Scalar>, Aligned> MCQuaternionA;
  154. typedef Map<Quaternion<Scalar> > MQuaternionUA;
  155. typedef Map<const Quaternion<Scalar> > MCQuaternionUA;
  156. typedef Quaternion<Scalar> Quaternionx;
  157. typedef Matrix<Scalar,3,1> Vector3;
  158. typedef AngleAxis<Scalar> AngleAxisx;
  159. Vector3 v0 = Vector3::Random(),
  160. v1 = Vector3::Random();
  161. Scalar a = internal::random<Scalar>(-Scalar(EIGEN_PI), Scalar(EIGEN_PI));
  162. EIGEN_ALIGN_MAX Scalar array1[4];
  163. EIGEN_ALIGN_MAX Scalar array2[4];
  164. EIGEN_ALIGN_MAX Scalar array3[4+1];
  165. Scalar* array3unaligned = array3+1;
  166. MQuaternionA mq1(array1);
  167. MCQuaternionA mcq1(array1);
  168. MQuaternionA mq2(array2);
  169. MQuaternionUA mq3(array3unaligned);
  170. MCQuaternionUA mcq3(array3unaligned);
  171. // std::cerr << array1 << " " << array2 << " " << array3 << "\n";
  172. mq1 = AngleAxisx(a, v0.normalized());
  173. mq2 = mq1;
  174. mq3 = mq1;
  175. Quaternionx q1 = mq1;
  176. Quaternionx q2 = mq2;
  177. Quaternionx q3 = mq3;
  178. Quaternionx q4 = MCQuaternionUA(array3unaligned);
  179. VERIFY_IS_APPROX(q1.coeffs(), q2.coeffs());
  180. VERIFY_IS_APPROX(q1.coeffs(), q3.coeffs());
  181. VERIFY_IS_APPROX(q4.coeffs(), q3.coeffs());
  182. VERIFY_IS_APPROX(mq1 * (mq1.inverse() * v1), v1);
  183. VERIFY_IS_APPROX(mq1 * (mq1.conjugate() * v1), v1);
  184. VERIFY_IS_APPROX(mcq1 * (mcq1.inverse() * v1), v1);
  185. VERIFY_IS_APPROX(mcq1 * (mcq1.conjugate() * v1), v1);
  186. VERIFY_IS_APPROX(mq3 * (mq3.inverse() * v1), v1);
  187. VERIFY_IS_APPROX(mq3 * (mq3.conjugate() * v1), v1);
  188. VERIFY_IS_APPROX(mcq3 * (mcq3.inverse() * v1), v1);
  189. VERIFY_IS_APPROX(mcq3 * (mcq3.conjugate() * v1), v1);
  190. VERIFY_IS_APPROX(mq1*mq2, q1*q2);
  191. VERIFY_IS_APPROX(mq3*mq2, q3*q2);
  192. VERIFY_IS_APPROX(mcq1*mq2, q1*q2);
  193. VERIFY_IS_APPROX(mcq3*mq2, q3*q2);
  194. // Bug 1461, compilation issue with Map<const Quat>::w(), and other reference/constness checks:
  195. VERIFY_IS_APPROX(mcq3.coeffs().x() + mcq3.coeffs().y() + mcq3.coeffs().z() + mcq3.coeffs().w(), mcq3.coeffs().sum());
  196. VERIFY_IS_APPROX(mcq3.x() + mcq3.y() + mcq3.z() + mcq3.w(), mcq3.coeffs().sum());
  197. mq3.w() = 1;
  198. const Quaternionx& cq3(q3);
  199. VERIFY( &cq3.x() == &q3.x() );
  200. const MQuaternionUA& cmq3(mq3);
  201. VERIFY( &cmq3.x() == &mq3.x() );
  202. // FIXME the following should be ok. The problem is that currently the LValueBit flag
  203. // is used to determine whether we can return a coeff by reference or not, which is not enough for Map<const ...>.
  204. //const MCQuaternionUA& cmcq3(mcq3);
  205. //VERIFY( &cmcq3.x() == &mcq3.x() );
  206. // test cast
  207. {
  208. Quaternion<float> q1f = mq1.template cast<float>();
  209. VERIFY_IS_APPROX(q1f.template cast<Scalar>(),mq1);
  210. Quaternion<double> q1d = mq1.template cast<double>();
  211. VERIFY_IS_APPROX(q1d.template cast<Scalar>(),mq1);
  212. }
  213. }
  214. template<typename Scalar> void quaternionAlignment(void){
  215. typedef Quaternion<Scalar,AutoAlign> QuaternionA;
  216. typedef Quaternion<Scalar,DontAlign> QuaternionUA;
  217. EIGEN_ALIGN_MAX Scalar array1[4];
  218. EIGEN_ALIGN_MAX Scalar array2[4];
  219. EIGEN_ALIGN_MAX Scalar array3[4+1];
  220. Scalar* arrayunaligned = array3+1;
  221. QuaternionA *q1 = ::new(reinterpret_cast<void*>(array1)) QuaternionA;
  222. QuaternionUA *q2 = ::new(reinterpret_cast<void*>(array2)) QuaternionUA;
  223. QuaternionUA *q3 = ::new(reinterpret_cast<void*>(arrayunaligned)) QuaternionUA;
  224. q1->coeffs().setRandom();
  225. *q2 = *q1;
  226. *q3 = *q1;
  227. VERIFY_IS_APPROX(q1->coeffs(), q2->coeffs());
  228. VERIFY_IS_APPROX(q1->coeffs(), q3->coeffs());
  229. }
  230. template<typename PlainObjectType> void check_const_correctness(const PlainObjectType&)
  231. {
  232. // there's a lot that we can't test here while still having this test compile!
  233. // the only possible approach would be to run a script trying to compile stuff and checking that it fails.
  234. // CMake can help with that.
  235. // verify that map-to-const don't have LvalueBit
  236. typedef typename internal::add_const<PlainObjectType>::type ConstPlainObjectType;
  237. VERIFY( !(internal::traits<Map<ConstPlainObjectType> >::Flags & LvalueBit) );
  238. VERIFY( !(internal::traits<Map<ConstPlainObjectType, Aligned> >::Flags & LvalueBit) );
  239. VERIFY( !(Map<ConstPlainObjectType>::Flags & LvalueBit) );
  240. VERIFY( !(Map<ConstPlainObjectType, Aligned>::Flags & LvalueBit) );
  241. }
  242. #if EIGEN_HAS_RVALUE_REFERENCES
  243. // Regression for bug 1573
  244. struct MovableClass {
  245. // The following line is a workaround for gcc 4.7 and 4.8 (see bug 1573 comments).
  246. static_assert(std::is_nothrow_move_constructible<Quaternionf>::value,"");
  247. MovableClass() = default;
  248. MovableClass(const MovableClass&) = default;
  249. MovableClass(MovableClass&&) noexcept = default;
  250. MovableClass& operator=(const MovableClass&) = default;
  251. MovableClass& operator=(MovableClass&&) = default;
  252. Quaternionf m_quat;
  253. };
  254. #endif
  255. EIGEN_DECLARE_TEST(geo_quaternion)
  256. {
  257. for(int i = 0; i < g_repeat; i++) {
  258. CALL_SUBTEST_1(( quaternion<float,AutoAlign>() ));
  259. CALL_SUBTEST_1( check_const_correctness(Quaternionf()) );
  260. CALL_SUBTEST_1(( quaternion<float,DontAlign>() ));
  261. CALL_SUBTEST_1(( quaternionAlignment<float>() ));
  262. CALL_SUBTEST_1( mapQuaternion<float>() );
  263. CALL_SUBTEST_2(( quaternion<double,AutoAlign>() ));
  264. CALL_SUBTEST_2( check_const_correctness(Quaterniond()) );
  265. CALL_SUBTEST_2(( quaternion<double,DontAlign>() ));
  266. CALL_SUBTEST_2(( quaternionAlignment<double>() ));
  267. CALL_SUBTEST_2( mapQuaternion<double>() );
  268. #ifndef EIGEN_TEST_ANNOYING_SCALAR_DONT_THROW
  269. AnnoyingScalar::dont_throw = true;
  270. #endif
  271. CALL_SUBTEST_3(( quaternion<AnnoyingScalar,AutoAlign>() ));
  272. }
  273. }