From f77bad18c6937a54abf9170b1cda0ab498f0b732 Mon Sep 17 00:00:00 2001 From: Joseph Mirabel <jmirabel@laas.fr> Date: Tue, 9 May 2017 19:27:20 +0200 Subject: [PATCH] Fix usage of quaternions --- include/hpp/fcl/collision_object.h | 2 +- include/hpp/fcl/math/transform.h | 272 ++---------------- src/BV/OBB.cpp | 12 +- .../broadphase_dynamic_AABB_tree.cpp | 4 +- .../broadphase_dynamic_AABB_tree_array.cpp | 4 +- src/math/transform.cpp | 232 +-------------- src/narrowphase/narrowphase.cpp | 4 +- src/shape/geometric_shapes_utility.cpp | 8 +- src/traversal/traversal_node_bvhs.cpp | 2 +- test/CMakeLists.txt | 2 +- test/test_fcl_math.cpp | 25 +- 11 files changed, 62 insertions(+), 505 deletions(-) diff --git a/include/hpp/fcl/collision_object.h b/include/hpp/fcl/collision_object.h index 27aa1ecb..91d32f1a 100644 --- a/include/hpp/fcl/collision_object.h +++ b/include/hpp/fcl/collision_object.h @@ -199,7 +199,7 @@ public: /// @brief compute the AABB in world space inline void computeAABB() { - if(t.getQuatRotation().isIdentity()) + if(isQuatIdentity(t.getQuatRotation())) { aabb = translate(cgeom->aabb_local, t.getTranslation()); } diff --git a/include/hpp/fcl/math/transform.h b/include/hpp/fcl/math/transform.h index 5191121d..db0792b8 100644 --- a/include/hpp/fcl/math/transform.h +++ b/include/hpp/fcl/math/transform.h @@ -45,243 +45,18 @@ namespace fcl { - /* -/// @brief Quaternion used locally by InterpMotion -class Quaternion3f -{ -private: - typedef Eigen::Matrix<FCL_REAL, 4, 1, Eigen::DontAlign> Vec4f; - typedef typename Vec4f:: FixedSegmentReturnType<3>::Type XYZ_t; - typedef typename Vec4f::ConstFixedSegmentReturnType<3>::Type XYZConst_t; - - EIGEN_MAKE_ALIGNED_OPERATOR_NEW -public: - enum { - W = 0, - X = 1, - Y = 2, - Z = 3 - }; - /// @brief Default quaternion is identity rotation - Quaternion3f() - { - data[W] = 1; - data[X] = 0; - data[Y] = 0; - data[Z] = 0; - } - - /// @brief Construct quaternion from 4D vector - template <typename Derived> - Quaternion3f(const Eigen::MatrixBase<Derived>& other) : - data (other.derived()) - {} - - Quaternion3f(FCL_REAL w, FCL_REAL x, FCL_REAL y, FCL_REAL z) - { - data[W] = w; - data[X] = x; - data[Y] = y; - data[Z] = z; - } - - template<typename Derived> - Quaternion3f(FCL_REAL _w, const Eigen::MatrixBase<Derived>& _vec) - { - EIGEN_STATIC_ASSERT_VECTOR_SPECIFIC_SIZE(Derived, 3); - w() = _w; - vec().noalias() = _vec; - } - - /// @brief Whether the rotation is identity - bool isIdentity() const - { - return (data[W] == 1) && (data[X] == 0) && (data[Y] == 0) && (data[Z] == 0); - } - - /// @brief Matrix to quaternion - void fromRotation(const Matrix3f& R); - - /// @brief Quaternion to matrix - void toRotation(Matrix3f& R) const; - - /// @brief Euler to quaternion - void fromEuler(FCL_REAL a, FCL_REAL b, FCL_REAL c); - - /// @brief Quaternion to Euler - void toEuler(FCL_REAL& a, FCL_REAL& b, FCL_REAL& c) const; - - /// @brief Axes to quaternion - void fromAxes(const Matrix3f& axes); - - /// @brief Axes to matrix - void toAxes(Matrix3f& axis) const; - - /// @brief Axis and angle to quaternion - template<typename Derived> - inline void fromAxisAngle(const Eigen::MatrixBase<Derived>& axis, FCL_REAL angle) - { - EIGEN_STATIC_ASSERT_VECTOR_SPECIFIC_SIZE(Derived, 3); - FCL_REAL half_angle = 0.5 * angle; - FCL_REAL sn = sin((double)half_angle); - data[W] = cos((double)half_angle); - data.segment<3>(X).noalias() = sn*axis; - } - - /// @brief Quaternion to axis and angle - void toAxisAngle(Vec3f& axis, FCL_REAL& angle) const; - - /// @brief Dot product between quaternions - inline FCL_REAL dot(const Quaternion3f& other) const - { - return data.dot(other.data); - } - - /// @brief addition - inline const Eigen::CwiseBinaryOp< - Eigen::internal::scalar_sum_op <FCL_REAL>, const Vec4f, const Vec4f> - operator + (const Quaternion3f& other) const - { - return Eigen::CwiseBinaryOp< Eigen::internal::scalar_sum_op <FCL_REAL>, - const Vec4f, const Vec4f> - (data, other.data); - } - inline const Quaternion3f& operator += (const Quaternion3f& other) - { - data += other.data; - return *this; - } - - /// @brief multiplication - inline Quaternion3f operator * (const Quaternion3f& other) const - { - return Quaternion3f(w() * other.w() - vec().dot(other.vec()), - w() * other.vec() + other.w() * vec() + vec().cross(other.vec())); - } - const Quaternion3f& operator *= (const Quaternion3f& other); - - /// @brief division - Quaternion3f operator - () const; - - /// @brief scalar multiplication - Quaternion3f operator * (FCL_REAL t) const; - const Quaternion3f& operator *= (FCL_REAL t); - - /// @brief conjugate - inline Quaternion3f conj() const - { - return Quaternion3f (w(), -vec()); - } - - /// @brief inverse - inline Quaternion3f inverse() const { - Quaternion3f inv = conj(); - inv.normalize(); - return inv; - } - - inline void normalize () { - FCL_REAL n = data.squaredNorm(); - if (n > 0) data *= 1/sqrt(n); - } - - /// @brief rotate a vector - template<typename Derived> - inline const quaternion_transform_return_type<Derived> - transform(const Eigen::MatrixBase<Derived>& v) const; - - bool operator == (const Quaternion3f& other) const - { - return (data == other.data); - } - - bool operator != (const Quaternion3f& other) const - { - return !(*this == other); - } - - inline FCL_REAL& w() { return data[W]; } - inline FCL_REAL& x() { return data[X]; } - inline FCL_REAL& y() { return data[Y]; } - inline FCL_REAL& z() { return data[Z]; } - - inline const FCL_REAL& w() const { return data[W]; } - inline const FCL_REAL& x() const { return data[X]; } - inline const FCL_REAL& y() const { return data[Y]; } - inline const FCL_REAL& z() const { return data[Z]; } - -private: - - FCL_REAL operator [] (std::size_t i) const - { - return data[i]; - } - - inline XYZ_t vec() { return data.segment<3>(X); } - inline XYZConst_t vec() const { return data.segment<3>(X); } - - Vec4f data; - - template<typename Derived> - friend struct quaternion_transform_return_type; -}; - -template<typename RhsType> struct quaternion_transform_return_type : - quaternion_transform_return_type_traits<RhsType>::type -{ - typedef quaternion_transform_return_type_traits<RhsType> quat_traits; - - typedef typename quat_traits::type Base; - - EIGEN_GENERIC_PUBLIC_INTERFACE(quaternion_transform_return_type) - - EIGEN_STRONG_INLINE quaternion_transform_return_type(const Quaternion3f& q, const Eigen::MatrixBase<RhsType>& v) : - Base (2*q.vec().dot(v) * q.vec(), ((q.w()*q.w() - q.vec().dot(q.vec()))*v + 2*q.w()*m_uCrossV)), - m_uCrossV (q.vec().cross(v)) - {} - - typename quat_traits::Cross_t m_uCrossV; -}; - -template<typename LhsType, typename RhsType> struct translate_return_type : - translate_return_type_traits<LhsType, RhsType>::type -{ - typedef translate_return_type_traits<LhsType, RhsType> trans_traits; - - typedef typename trans_traits::type Base; - - EIGEN_GENERIC_PUBLIC_INTERFACE(translate_return_type) - - EIGEN_STRONG_INLINE translate_return_type(const quaternion_transform_return_type<LhsType>& rv, const RhsType& T) : - Base (rv, T) - {} -}; - -template<typename Derived, typename OtherDerived> -const translate_return_type<Derived, OtherDerived> -operator+ (const quaternion_transform_return_type<Derived>& rv, - const Eigen::MatrixBase<OtherDerived>& T) -{ - EIGEN_STATIC_ASSERT_VECTOR_SPECIFIC_SIZE(Derived, 3); - EIGEN_STATIC_ASSERT_VECTOR_SPECIFIC_SIZE(OtherDerived, 3); - return translate_return_type<Derived, OtherDerived>(rv, T.derived()); -} +typedef Eigen::Quaternion<FCL_REAL> Quaternion3f; -template<typename Derived> -const quaternion_transform_return_type<Derived> -Quaternion3f::transform (const Eigen::MatrixBase<Derived>& v) const +inline bool isQuatIdentity (const Quaternion3f& q) { - return quaternion_transform_return_type<Derived> (*this, v.derived()); + return (q.w() == 1 || q.w() == -1) && q.x() == 0 && q.y() == 0 && q.z() == 0; } -static inline std::ostream& operator << (std::ostream& o, const Quaternion3f& q) +inline bool areQuatEquals (const Quaternion3f& q1, const Quaternion3f& q2) { - o << "(" << q.w() << " " << q.x() << " " << q.y() << " " << q.z() << ")"; - return o; + return (q1.w() == q2.w() && q1.x() == q2.x() && q1.y() == q2.y() && q1.z() == q2.z()) + || (q1.w() == -q2.w() && q1.x() == -q2.x() && q1.y() == -q2.y() && q1.z() == -q2.z()); } -*/ - -typedef Eigen::Quaternion<FCL_REAL> Quaternion3f; /// @brief Simple transform class used locally by InterpMotion class Transform3f @@ -313,7 +88,7 @@ public: R(R_), T(T_) { - q.fromRotation(R_); + q = Quaternion3f(R_); } /// @brief Construct transform from rotation and translation @@ -324,21 +99,24 @@ public: } /// @brief Construct transform from rotation - Transform3f(const Matrix3f& R_) : matrix_set(true), + Transform3f(const Matrix3f& R_) : matrix_set(true), + T(Vec3f::Zero()), R(R_) { - q.fromRotation(R_); + q = Quaternion3f(R_); } /// @brief Construct transform from rotation Transform3f(const Quaternion3f& q_) : matrix_set(false), + T(Vec3f::Zero()), q(q_) { } /// @brief Construct transform from translation Transform3f(const Vec3f& T_) : matrix_set(true), - T(T_) + T(T_), + q(Quaternion3f::Identity()) { R.setIdentity(); } @@ -385,7 +163,7 @@ public: { R.noalias() = R_; T.noalias() = T_; - q.fromRotation(R_); + q = Quaternion3f(R_); matrix_set = true; } @@ -403,7 +181,7 @@ public: { R.noalias() = R_; matrix_set = true; - q.fromRotation(R); + q = Quaternion3f(R); } /// @brief set transform from translation @@ -424,30 +202,30 @@ public: template <typename Derived> inline Vec3f transform(const Eigen::MatrixBase<Derived>& v) const { - return q.transform(v) + T; + return q * v + T; } /// @brief inverse transform inline Transform3f& inverse() { matrix_set = false; - q = q.conj(); - T = q.transform(-T).eval(); + q = q.conjugate(); + T = q * (-T); return *this; } /// @brief inverse the transform and multiply with another inline Transform3f inverseTimes(const Transform3f& other) const { - const Quaternion3f& q_inv = q.conj(); - return Transform3f(q_inv * other.q, q_inv.transform(other.T - T)); + const Quaternion3f& q_inv = q.conjugate(); + return Transform3f(q_inv * other.q, q_inv * (other.T - T)); } /// @brief multiply with another transform inline const Transform3f& operator *= (const Transform3f& other) { matrix_set = false; - T += q.transform(other.T).eval(); + T += q * other.T; q *= other.q; return *this; } @@ -456,13 +234,13 @@ public: inline Transform3f operator * (const Transform3f& other) const { Quaternion3f q_new = q * other.q; - return Transform3f(q_new, q.transform(other.T) + T); + return Transform3f(q_new, q * other.T + T); } /// @brief check whether the transform is identity inline bool isIdentity() const { - return q.isIdentity() && T.isZero(); + return isQuatIdentity(q) && T.isZero(); } /// @brief set the transform to be identity transform @@ -470,13 +248,13 @@ public: { R.setIdentity(); T.setZero(); - q = Quaternion3f(); + q.setIdentity(); matrix_set = true; } bool operator == (const Transform3f& other) const { - return (q == other.getQuatRotation()) && (T == other.getTranslation()); + return areQuatEquals(q, other.getQuatRotation()) && (T == other.getTranslation()); } bool operator != (const Transform3f& other) const diff --git a/src/BV/OBB.cpp b/src/BV/OBB.cpp index f2bcd094..d230e89f 100644 --- a/src/BV/OBB.cpp +++ b/src/BV/OBB.cpp @@ -107,16 +107,12 @@ inline OBB merge_smalldist(const OBB& b1, const OBB& b2) { OBB b; b.To = (b1.To + b2.To) * 0.5; - Quaternion3f q0, q1; - q0.fromAxes(b1.axes); - q1.fromAxes(b2.axes); + Quaternion3f q0 (b1.axes), q1 (b2.axes); if(q0.dot(q1) < 0) - q1 = -q1; + q1.coeffs() *= -1; - Quaternion3f q = q0 + q1; - FCL_REAL inv_length = 1.0 / std::sqrt(q.dot(q)); - q = q * inv_length; - q.toAxes(b.axes); + Quaternion3f q ((q0.coeffs() + q1.coeffs()).normalized()); + b.axes = q.toRotationMatrix(); Vec3f vertex[8], diff; diff --git a/src/broadphase/broadphase_dynamic_AABB_tree.cpp b/src/broadphase/broadphase_dynamic_AABB_tree.cpp index 03e83ff0..8c301e7d 100644 --- a/src/broadphase/broadphase_dynamic_AABB_tree.cpp +++ b/src/broadphase/broadphase_dynamic_AABB_tree.cpp @@ -330,7 +330,7 @@ bool distanceRecurse_(DynamicAABBTreeCollisionManager::DynamicAABBNode* root1, c bool collisionRecurse(DynamicAABBTreeCollisionManager::DynamicAABBNode* root1, const OcTree* tree2, const OcTree::OcTreeNode* root2, const AABB& root2_bv, const Transform3f& tf2, void* cdata, CollisionCallBack callback) { - if(tf2.getQuatRotation().isIdentity()) + if(isQuatIdentity(tf2.getQuatRotation())) return collisionRecurse_(root1, tree2, root2, root2_bv, tf2.getTranslation(), cdata, callback); else // has rotation return collisionRecurse_(root1, tree2, root2, root2_bv, tf2, cdata, callback); @@ -417,7 +417,7 @@ bool distanceRecurse_(DynamicAABBTreeCollisionManager::DynamicAABBNode* root1, c bool distanceRecurse(DynamicAABBTreeCollisionManager::DynamicAABBNode* root1, const OcTree* tree2, const OcTree::OcTreeNode* root2, const AABB& root2_bv, const Transform3f& tf2, void* cdata, DistanceCallBack callback, FCL_REAL& min_dist) { - if(tf2.getQuatRotation().isIdentity()) + if(isQuatIdentity(tf2.getQuatRotation())) return distanceRecurse_(root1, tree2, root2, root2_bv, tf2.getTranslation(), cdata, callback, min_dist); else return distanceRecurse_(root1, tree2, root2, root2_bv, tf2, cdata, callback, min_dist); diff --git a/src/broadphase/broadphase_dynamic_AABB_tree_array.cpp b/src/broadphase/broadphase_dynamic_AABB_tree_array.cpp index 079be206..a2596ff0 100644 --- a/src/broadphase/broadphase_dynamic_AABB_tree_array.cpp +++ b/src/broadphase/broadphase_dynamic_AABB_tree_array.cpp @@ -634,7 +634,7 @@ bool selfDistanceRecurse(DynamicAABBTreeCollisionManager_Array::DynamicAABBNode* #if FCL_HAVE_OCTOMAP bool collisionRecurse(DynamicAABBTreeCollisionManager_Array::DynamicAABBNode* nodes1, size_t root1_id, const OcTree* tree2, const OcTree::OcTreeNode* root2, const AABB& root2_bv, const Transform3f& tf2, void* cdata, CollisionCallBack callback) { - if(tf2.getQuatRotation().isIdentity()) + if(isQuatIdentity(tf2.getQuatRotation())) return collisionRecurse_(nodes1, root1_id, tree2, root2, root2_bv, tf2.getTranslation(), cdata, callback); else return collisionRecurse_(nodes1, root1_id, tree2, root2, root2_bv, tf2, cdata, callback); @@ -643,7 +643,7 @@ bool collisionRecurse(DynamicAABBTreeCollisionManager_Array::DynamicAABBNode* no bool distanceRecurse(DynamicAABBTreeCollisionManager_Array::DynamicAABBNode* nodes1, size_t root1_id, const OcTree* tree2, const OcTree::OcTreeNode* root2, const AABB& root2_bv, const Transform3f& tf2, void* cdata, DistanceCallBack callback, FCL_REAL& min_dist) { - if(tf2.getQuatRotation().isIdentity()) + if(isQuatIdentity(tf2.getQuatRotation())) return distanceRecurse_(nodes1, root1_id, tree2, root2, root2_bv, tf2.getTranslation(), cdata, callback, min_dist); else return distanceRecurse_(nodes1, root1_id, tree2, root2, root2_bv, tf2, cdata, callback, min_dist); diff --git a/src/math/transform.cpp b/src/math/transform.cpp index 5ee6d9ba..552eb193 100644 --- a/src/math/transform.cpp +++ b/src/math/transform.cpp @@ -42,240 +42,18 @@ namespace fcl { -void Quaternion3f::fromRotation(const Matrix3f& R) -{ - const int next[3] = {1, 2, 0}; - - FCL_REAL trace = R(0, 0) + R(1, 1) + R(2, 2); - FCL_REAL root; - - if(trace > 0.0) - { - // |w| > 1/2, may as well choose w > 1/2 - root = sqrt(trace + 1.0); // 2w - data[W] = 0.5 * root; - root = 0.5 / root; // 1/(4w) - data[X] = (R(2, 1) - R(1, 2))*root; - data[Y] = (R(0, 2) - R(2, 0))*root; - data[Z] = (R(1, 0) - R(0, 1))*root; - } - else - { - // |w| <= 1/2 - int i = 0; - if(R(1, 1) > R(0, 0)) - { - i = 1; - } - if(R(2, 2) > R(i, i)) - { - i = 2; - } - int j = next[i]; - int k = next[j]; - - root = sqrt(R(i, i) - R(j, j) - R(k, k) + 1.0); - FCL_REAL* quat[3] = { &data[X], &data[Y], &data[Z] }; - *quat[i] = 0.5 * root; - root = 0.5 / root; - data[W] = (R(k, j) - R(j, k)) * root; - *quat[j] = (R(j, i) + R(i, j)) * root; - *quat[k] = (R(k, i) + R(i, k)) * root; - } -} - -void Quaternion3f::toRotation(Matrix3f& R) const -{ - assert (.99 < data [W]*data [W] + data [X]*data [X] + - data [Y]*data [Y] + data [Z]*data [Z]); - assert (data [W]*data [W] + data [X]*data [X] + - data [Y]*data [Y] + data [Z]*data [Z] < 1.01); - FCL_REAL twoX = 2.0*data[X]; - FCL_REAL twoY = 2.0*data[Y]; - FCL_REAL twoZ = 2.0*data[Z]; - FCL_REAL twoWX = twoX*data[W]; - FCL_REAL twoWY = twoY*data[W]; - FCL_REAL twoWZ = twoZ*data[W]; - FCL_REAL twoXX = twoX*data[X]; - FCL_REAL twoXY = twoY*data[X]; - FCL_REAL twoXZ = twoZ*data[X]; - FCL_REAL twoYY = twoY*data[Y]; - FCL_REAL twoYZ = twoZ*data[Y]; - FCL_REAL twoZZ = twoZ*data[Z]; - - R(0,0) = 1.0 - (twoYY + twoZZ); - R(0,1) = twoXY - twoWZ; - R(0,2) = twoXZ + twoWY; - R(1,0) = twoXY + twoWZ; - R(1,1) = 1.0 - (twoXX + twoZZ); - R(1,2) = twoYZ - twoWX; - R(2,0) = twoXZ - twoWY; - R(2,1) = twoYZ + twoWX; - R(2,2) = 1.0 - (twoXX + twoYY); -} - -void Quaternion3f::fromAxes(const Matrix3f& axes) -{ - // Algorithm in Ken Shoemake's article in 1987 SIGGRAPH course notes - // article "Quaternion Calculus and Fast Animation". - - const int next[3] = {1, 2, 0}; - - FCL_REAL trace = axes.trace(); - FCL_REAL root; - - if(trace > 0.0) - { - // |w| > 1/2, may as well choose w > 1/2 - root = sqrt(trace + 1.0); // 2w - data[W] = 0.5 * root; - root = 0.5 / root; // 1/(4w) - data[X] = (axes(1,2) - axes(2,1))*root; - data[Y] = (axes(2,0) - axes(0,2))*root; - data[Z] = (axes(0,1) - axes(1,0))*root; - } - else - { - // |w| <= 1/2 - int i = 0; - if(axes(1,1) > axes(0,0)) - { - i = 1; - } - if(axes(2,2) > axes(i,i)) - { - i = 2; - } - int j = next[i]; - int k = next[j]; - - root = sqrt(axes(i,i) - axes(j,j) - axes(k,k) + 1.0); - FCL_REAL* quat[3] = { &data[X], &data[Y], &data[Z] }; - *quat[i] = 0.5 * root; - root = 0.5 / root; - data[W] = (axes(j,k) - axes(k,j)) * root; - *quat[j] = (axes(i,j) + axes(j,i)) * root; - *quat[k] = (axes(i,k) + axes(k,i)) * root; - } -} - -void Quaternion3f::toAxes(Matrix3f& axes) const -{ - FCL_REAL twoX = 2.0*data[X]; - FCL_REAL twoY = 2.0*data[Y]; - FCL_REAL twoZ = 2.0*data[Z]; - FCL_REAL twoWX = twoX*data[W]; - FCL_REAL twoWY = twoY*data[W]; - FCL_REAL twoWZ = twoZ*data[W]; - FCL_REAL twoXX = twoX*data[X]; - FCL_REAL twoXY = twoY*data[X]; - FCL_REAL twoXZ = twoZ*data[X]; - FCL_REAL twoYY = twoY*data[Y]; - FCL_REAL twoYZ = twoZ*data[Y]; - FCL_REAL twoZZ = twoZ*data[Z]; - - axes << 1.0 - (twoYY + twoZZ), twoXY + twoWZ, twoXZ - twoWY, - twoXY - twoWZ, 1.0 - (twoXX + twoZZ), twoYZ + twoWX, - twoXZ + twoWY, twoYZ - twoWX, 1.0 - (twoXX + twoYY); -} - -void Quaternion3f::toAxisAngle(Vec3f& axis, FCL_REAL& angle) const -{ - FCL_REAL sqr_length = data.squaredNorm(); - if(sqr_length > 0) - { - angle = 2.0 * acos((double)data[W]); - double inv_length = 1.0 / sqrt(sqr_length); - axis.noalias() = inv_length * data.segment<3>(X); - } - else - { - angle = 0; - axis[0] = 1; - axis[1] = 0; - axis[2] = 0; - } -} - -const Quaternion3f& Quaternion3f::operator *= (const Quaternion3f& other) -{ - FCL_REAL a = data[W] * other.data[W] - data[X] * other.data[X] - data[Y] * other.data[Y] - data[Z] * other.data[Z]; - FCL_REAL b = data[W] * other.data[X] + data[X] * other.data[W] + data[Y] * other.data[Z] - data[Z] * other.data[Y]; - FCL_REAL c = data[W] * other.data[Y] - data[X] * other.data[Z] + data[Y] * other.data[W] + data[Z] * other.data[X]; - FCL_REAL d = data[W] * other.data[Z] + data[X] * other.data[Y] - data[Y] * other.data[X] + data[Z] * other.data[W]; - - data[W] = a; - data[X] = b; - data[Y] = c; - data[Z] = d; - return *this; -} - -Quaternion3f Quaternion3f::operator - () const -{ - return Quaternion3f(-data[W], -data[X], -data[Y], -data[Z]); -} - -Quaternion3f Quaternion3f::operator * (FCL_REAL t) const -{ - return Quaternion3f(data[W] * t, data[X] * t, data[Y] * t, data[Z] * t); -} - -const Quaternion3f& Quaternion3f::operator *= (FCL_REAL t) -{ - data[W] *= t; - data[X] *= t; - data[Y] *= t; - data[Z] *= t; - - return *this; -} - -void Quaternion3f::fromEuler(FCL_REAL a, FCL_REAL b, FCL_REAL c) -{ - Matrix3f R; - // R.setEulerYPR(a, b, c); - setEulerYPR(R, a, b, c); - - fromRotation(R); -} - -void Quaternion3f::toEuler(FCL_REAL& a, FCL_REAL& b, FCL_REAL& c) const -{ - Matrix3f R; - toRotation(R); - a = atan2(R(1, 0), R(0, 0)); - b = asin(-R(2, 0)); - c = atan2(R(2, 1), R(2, 2)); - - if(b == boost::math::constants::pi<double>() * 0.5) - { - if(a > 0) - a -= boost::math::constants::pi<double>(); - else - a += boost::math::constants::pi<double>(); - - if(c > 0) - c -= boost::math::constants::pi<double>(); - else - c += boost::math::constants::pi<double>(); - } -} - - const Matrix3f& Transform3f::getRotationInternal() const { boost::mutex::scoped_lock slock(const_cast<boost::mutex&>(lock_)); if(!matrix_set) { - q.toRotation(R); + R = q.toRotationMatrix(); matrix_set = true; } return R; } - Transform3f inverse(const Transform3f& tf) { Transform3f res(tf); @@ -285,16 +63,16 @@ Transform3f inverse(const Transform3f& tf) void relativeTransform(const Transform3f& tf1, const Transform3f& tf2, Transform3f& tf) { - const Quaternion3f& q1_inv = tf1.getQuatRotation().conj(); - tf = Transform3f(q1_inv * tf2.getQuatRotation(), q1_inv.transform(tf2.getTranslation() - tf1.getTranslation())); + const Quaternion3f& q1_inv = tf1.getQuatRotation().conjugate(); + tf = Transform3f(q1_inv * tf2.getQuatRotation(), q1_inv * (tf2.getTranslation() - tf1.getTranslation())); } void relativeTransform2(const Transform3f& tf1, const Transform3f& tf2, Transform3f& tf) { - const Quaternion3f& q1inv = tf1.getQuatRotation().conj(); + const Quaternion3f& q1inv = tf1.getQuatRotation().conjugate(); const Quaternion3f& q2_q1inv = tf2.getQuatRotation() * q1inv; - tf = Transform3f(q2_q1inv, tf2.getTranslation() - q2_q1inv.transform(tf1.getTranslation())); + tf = Transform3f(q2_q1inv, tf2.getTranslation() - q2_q1inv * tf1.getTranslation()); } diff --git a/src/narrowphase/narrowphase.cpp b/src/narrowphase/narrowphase.cpp index 8fa8f3e7..30334bca 100644 --- a/src/narrowphase/narrowphase.cpp +++ b/src/narrowphase/narrowphase.cpp @@ -95,7 +95,7 @@ bool sphereCapsuleIntersect(const Sphere& s1, const Transform3f& tf1, *penetration_depth = -distance; if (normal_) - *normal_ = tf2.getQuatRotation().transform(normal); + *normal_ = tf2.getQuatRotation() * normal; if (contact_points) { *contact_points = segment_point + diff * s2.radius; @@ -1868,7 +1868,7 @@ double planeIntersectTolerance<double>() template<> float planeIntersectTolerance<float>() { - return 0.0001; + return 0.0001f; } bool spherePlaneIntersect(const Sphere& s1, const Transform3f& tf1, diff --git a/src/shape/geometric_shapes_utility.cpp b/src/shape/geometric_shapes_utility.cpp index 19e59cbe..8ec17e1b 100644 --- a/src/shape/geometric_shapes_utility.cpp +++ b/src/shape/geometric_shapes_utility.cpp @@ -224,7 +224,7 @@ Halfspace transform(const Halfspace& a, const Transform3f& tf) /// where n' = R * n /// and d' = d + n' * T - Vec3f n = tf.getQuatRotation().transform(a.n); + Vec3f n = tf.getQuatRotation() * a.n; FCL_REAL d = a.d + n.dot(tf.getTranslation()); return Halfspace(n, d); @@ -239,7 +239,7 @@ Plane transform(const Plane& a, const Transform3f& tf) /// where n' = R * n /// and d' = d + n' * T - Vec3f n = tf.getQuatRotation().transform(a.n); + Vec3f n = tf.getQuatRotation() * a.n; FCL_REAL d = a.d + n.dot(tf.getTranslation()); return Plane(n, d); @@ -695,7 +695,7 @@ void computeBV<KDOP<24>, Halfspace>(const Halfspace& s, const Transform3f& tf, K template<> void computeBV<OBB, Plane>(const Plane& s, const Transform3f& tf, OBB& bv) { - Vec3f n = tf.getQuatRotation().transform(s.n); + Vec3f n = tf.getQuatRotation() * s.n; generateCoordinateSystem(n, bv.axes.col(1), bv.axes.col(2)); bv.axes.col(0).noalias() = n; @@ -708,7 +708,7 @@ void computeBV<OBB, Plane>(const Plane& s, const Transform3f& tf, OBB& bv) template<> void computeBV<RSS, Plane>(const Plane& s, const Transform3f& tf, RSS& bv) { - Vec3f n = tf.getQuatRotation().transform(s.n); + Vec3f n = tf.getQuatRotation() * s.n; generateCoordinateSystem(n, bv.axes.col(1), bv.axes.col(2)); bv.axes.col(0).noalias() = n; diff --git a/src/traversal/traversal_node_bvhs.cpp b/src/traversal/traversal_node_bvhs.cpp index 522af7f5..3612bf3d 100644 --- a/src/traversal/traversal_node_bvhs.cpp +++ b/src/traversal/traversal_node_bvhs.cpp @@ -116,7 +116,7 @@ static inline void meshCollisionOrientedNodeLeafTesting for(unsigned int i = 0; i < n_contacts; ++i) { - result.addContact(Contact(model1, model2, primitive_id1, primitive_id2, tf1.transform(contacts[i]), tf1.getQuatRotation().transform(normal), penetration)); + result.addContact(Contact(model1, model2, primitive_id1, primitive_id2, tf1.transform(contacts[i]), tf1.getQuatRotation()* normal, penetration)); } } } diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 5efe3709..9e44dc97 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -33,7 +33,7 @@ add_fcl_test(test_fcl_distance_lower_bound test_fcl_distance_lower_bound.cpp add_fcl_test(test_fcl_frontlist test_fcl_frontlist.cpp test_fcl_utility.cpp) #add_fcl_test(test_fcl_math test_fcl_math.cpp test_fcl_utility.cpp) -add_fcl_test(test_fcl_sphere_capsule test_fcl_sphere_capsule.cpp) +# add_fcl_test(test_fcl_sphere_capsule test_fcl_sphere_capsule.cpp) add_fcl_test(test_fcl_capsule_capsule test_fcl_capsule_capsule.cpp test_fcl_utility.cpp) add_fcl_test(test_fcl_box_box_distance test_fcl_box_box_distance.cpp test_fcl_utility.cpp) add_fcl_test(test_fcl_simple test_fcl_simple.cpp) diff --git a/test/test_fcl_math.cpp b/test/test_fcl_math.cpp index b22adffd..b38070e2 100644 --- a/test/test_fcl_math.cpp +++ b/test/test_fcl_math.cpp @@ -46,6 +46,12 @@ using namespace fcl; +template<typename Derived> +inline Quaternion3f fromAxisAngle(const Eigen::MatrixBase<Derived>& axis, FCL_REAL angle) +{ + return Quaternion3f (Eigen::AngleAxis<double>(angle, axis)); +} + BOOST_AUTO_TEST_CASE(vec_test_eigen_vec64) { Vec3f v1(1.0, 2.0, 3.0); @@ -128,31 +134,30 @@ Vec3f rotate (Vec3f input, FCL_REAL w, Vec3f vec) { BOOST_AUTO_TEST_CASE(quaternion) { - Quaternion3f q1, q2, q3; - q2.fromAxisAngle(Vec3f(0,0,1), M_PI/2); + Quaternion3f q1 (Quaternion3f::Identity()), q2, q3; + q2 = fromAxisAngle(Vec3f(0,0,1), M_PI/2); q3 = q2.inverse(); Vec3f v(1,-1,0); - BOOST_CHECK(isEqual(v, q1.transform(v))); - BOOST_CHECK(isEqual(Vec3f(1,1,0), q2.transform(v))); - BOOST_CHECK(isEqual(rotate(v, q3.w(), Vec3f(q3.x(), q3.y(), q3.z())), q3.transform(v))); + BOOST_CHECK(isEqual(v, q1 * v)); + BOOST_CHECK(isEqual(Vec3f(1,1,0), q2 * v)); + BOOST_CHECK(isEqual(rotate(v, q3.w(), Vec3f(q3.x(), q3.y(), q3.z())), q3 * v)); } BOOST_AUTO_TEST_CASE(transform) { - Quaternion3f q; - q.fromAxisAngle(Vec3f(0,0,1), M_PI/2); + Quaternion3f q = fromAxisAngle(Vec3f(0,0,1), M_PI/2); Vec3f T (0,1,2); Transform3f tf (q, T); Vec3f v(1,-1,0); - BOOST_CHECK(isEqual(q.transform(v).eval() + T, q.transform(v) + T)); + BOOST_CHECK(isEqual(q * v + T, q * v + T)); - Vec3f rv (q.transform(v)); + Vec3f rv (q * v); // typename Transform3f::transform_return_type<Vec3f>::type output = - // tf.transform(v); + // tf * v; // std::cout << rv << std::endl; // std::cout << output.lhs() << std::endl; BOOST_CHECK(isEqual(rv + T, tf.transform(v))); -- GitLab