6 #ifndef __pinocchio_inertia_hpp__
7 #define __pinocchio_inertia_hpp__
11 #include "pinocchio/math/fwd.hpp"
12 #include "pinocchio/spatial/symmetric3.hpp"
13 #include "pinocchio/spatial/force.hpp"
14 #include "pinocchio/spatial/motion.hpp"
15 #include "pinocchio/spatial/skew.hpp"
20 template<
class Derived>
25 typedef Derived Derived_t;
26 SPATIAL_TYPEDEF_TEMPLATE(Derived_t);
29 Derived_t & derived() {
return *
static_cast<Derived_t*
>(
this); }
30 const Derived_t & derived()
const {
return *
static_cast<const Derived_t*
>(
this); }
32 Scalar mass()
const {
return static_cast<const Derived_t*
>(
this)->mass(); }
33 Scalar & mass() {
return static_cast<const Derived_t*
>(
this)->mass(); }
34 const Vector3 & lever()
const {
return static_cast<const Derived_t*
>(
this)->lever(); }
35 Vector3 & lever() {
return static_cast<const Derived_t*
>(
this)->lever(); }
36 const Symmetric3 & inertia()
const {
return static_cast<const Derived_t*
>(
this)->inertia(); }
37 Symmetric3 & inertia() {
return static_cast<const Derived_t*
>(
this)->inertia(); }
39 Matrix6 matrix()
const {
return derived().matrix_impl(); }
40 operator Matrix6 ()
const {
return matrix(); }
42 Derived_t& operator= (
const Derived_t& clone){
return derived().__equl__(clone);}
43 bool operator==(
const Derived_t & other)
const {
return derived().isEqual(other);}
44 bool operator!=(
const Derived_t & other)
const {
return !(*
this == other); }
46 Derived_t& operator+= (
const Derived_t & Yb) {
return derived().__pequ__(Yb); }
47 Derived_t operator+(
const Derived_t & Yb)
const {
return derived().__plus__(Yb); }
49 template<
typename MotionDerived>
52 {
return derived().__mult__(v); }
54 Scalar vtiv(
const Motion & v)
const {
return derived().vtiv_impl(v); }
55 Matrix6 variation(
const Motion & v)
const {
return derived().variation_impl(v); }
65 static void vxi(
const Motion & v,
const Derived & I,
const Eigen::MatrixBase<M6> & Iout)
67 EIGEN_STATIC_ASSERT_SAME_MATRIX_SIZE(M6, Matrix6);
68 Derived::vxi_impl(v,I,Iout);
74 vxi(v,derived(),Iout);
86 static void ivx(
const Motion & v,
const Derived & I,
const Eigen::MatrixBase<M6> & Iout)
88 EIGEN_STATIC_ASSERT_SAME_MATRIX_SIZE(M6, Matrix6);
89 Derived::ivx_impl(v,I,Iout);
95 ivx(v,derived(),Iout);
99 void setZero() { derived().setZero(); }
100 void setIdentity() { derived().setIdentity(); }
101 void setRandom() { derived().setRandom(); }
103 bool isApprox(
const Derived & other,
const Scalar & prec = Eigen::NumTraits<Scalar>::dummy_precision())
const
104 {
return derived().isApprox_impl(other, prec); }
106 bool isZero(
const Scalar & prec = Eigen::NumTraits<Scalar>::dummy_precision())
const
107 {
return derived().isZero_impl(prec); }
110 Derived_t
se3Action(
const SE3 & M)
const {
return derived().se3Action_impl(M); }
115 void disp(std::ostream & os)
const {
static_cast<const Derived_t*
>(
this)->disp_impl(os); }
116 friend std::ostream & operator << (std::ostream & os,
const InertiaBase<Derived_t> & X)
125 template<
typename T,
int U>
129 typedef Eigen::Matrix<T,3,1,U> Vector3;
130 typedef Eigen::Matrix<T,4,1,U> Vector4;
131 typedef Eigen::Matrix<T,6,1,U> Vector6;
132 typedef Eigen::Matrix<T,3,3,U> Matrix3;
133 typedef Eigen::Matrix<T,4,4,U> Matrix4;
134 typedef Eigen::Matrix<T,6,6,U> Matrix6;
135 typedef Matrix6 ActionMatrix_t;
136 typedef Vector3 Angular_t;
137 typedef Vector3 Linear_t;
138 typedef const Vector3 ConstAngular_t;
139 typedef const Vector3 ConstLinear_t;
140 typedef Eigen::Quaternion<T,U> Quaternion_t;
151 template<
typename _Scalar,
int _Options>
157 enum { Options = _Options };
158 EIGEN_MAKE_ALIGNED_OPERATOR_NEW
160 typedef typename Symmetric3::AlphaSkewSquare AlphaSkewSquare;
162 typedef typename Eigen::Matrix<_Scalar, 10, 1, _Options> Vector10;
169 InertiaTpl(
const Scalar mass,
const Vector3 & com,
const Matrix3 & rotational_inertia)
170 : m_mass(mass), m_com(com), m_inertia(rotational_inertia)
173 InertiaTpl(
const Matrix6 & I6)
175 assert((I6 - I6.transpose()).isMuchSmallerThan(I6));
176 mass() = I6(LINEAR, LINEAR);
177 const Matrix3 & mc_cross = I6.template block <3,3>(ANGULAR,LINEAR);
178 lever() =
unSkew(mc_cross);
181 Matrix3 I3 (mc_cross * mc_cross);
183 I3 += I6.template block<3,3>(ANGULAR,ANGULAR);
184 const Symmetric3 S3(I3);
188 InertiaTpl(Scalar mass,
const Vector3 & com,
const Symmetric3 & rotational_inertia)
189 : m_mass(mass), m_com(com), m_inertia(rotational_inertia)
192 InertiaTpl(
const InertiaTpl & clone)
193 : m_mass(clone.mass()), m_com(clone.lever()), m_inertia(clone.inertia())
196 InertiaTpl& operator=(
const InertiaTpl & clone)
198 m_mass = clone.mass();
199 m_com = clone.lever();
200 m_inertia = clone.inertia();
205 InertiaTpl(
const InertiaTpl<Scalar,O2> & clone)
206 : m_mass(clone.mass())
207 , m_com(clone.lever())
208 , m_inertia(clone.inertia().matrix())
212 static InertiaTpl Zero()
214 return InertiaTpl(Scalar(0),
219 void setZero() { mass() = Scalar(0); lever().setZero(); inertia().setZero(); }
221 static InertiaTpl Identity()
223 return InertiaTpl(Scalar(1),
225 Symmetric3::Identity());
230 mass() = Scalar(1); lever().setZero(); inertia().setIdentity();
233 static InertiaTpl Random()
236 return InertiaTpl(Eigen::internal::random<Scalar>()+1,
238 Symmetric3::RandomPositive());
265 const Scalar a = mass * (y*y + z*z) / Scalar(5);
266 const Scalar b = mass * (x*x + z*z) / Scalar(5);
267 const Scalar c = mass * (y*y + x*x) / Scalar(5);
269 Scalar(0), Scalar(0), c));
283 const Scalar radius_square = radius * radius;
284 const Scalar a = mass * (radius_square / Scalar(4) + length*length / Scalar(12));
285 const Scalar c = mass * (radius_square / Scalar(2));
287 Scalar(0), Scalar(0), c));
303 const Scalar a = mass * (y*y + z*z) / Scalar(12);
304 const Scalar b = mass * (x*x + z*z) / Scalar(12);
305 const Scalar c = mass * (y*y + x*x) / Scalar(12);
307 Scalar(0), Scalar(0), c));
312 mass() =
static_cast<Scalar
>(std::rand())/
static_cast<Scalar
>(RAND_MAX);
313 lever().setRandom(); inertia().setRandom();
316 Matrix6 matrix_impl()
const
320 M.template block<3,3>(LINEAR, LINEAR ).setZero();
321 M.template block<3,3>(LINEAR, LINEAR ).diagonal().fill (mass());
322 M.template block<3,3>(ANGULAR,LINEAR ) =
alphaSkew(mass(),lever());
323 M.template block<3,3>(LINEAR, ANGULAR) = -M.template block<3,3>(ANGULAR, LINEAR);
324 M.template block<3,3>(ANGULAR,ANGULAR) = (inertia() - AlphaSkewSquare(mass(),lever())).matrix();
340 v.template segment<3>(1).noalias() = mass() * lever();
341 v.template segment<6>(4) = (inertia() - AlphaSkewSquare(mass(),lever())).data();
354 template<
typename Vector10Like>
357 PINOCCHIO_ASSERT_MATRIX_SPECIFIC_SIZE(Vector10Like, params, 10, 1);
359 const Scalar & mass = params[0];
360 Vector3 lever = params.template segment<3>(1);
363 return InertiaTpl(mass, lever,
Symmetric3(params.template segment<6>(4)) + AlphaSkewSquare(mass,lever));
369 mass() = clone.mass(); lever() = clone.lever(); inertia() = clone.inertia();
374 bool isEqual(
const InertiaTpl& Y2 )
const
376 return (mass()==Y2.mass()) && (lever()==Y2.lever()) && (inertia()==Y2.inertia());
379 bool isApprox_impl(
const InertiaTpl & other,
380 const Scalar & prec = Eigen::NumTraits<Scalar>::dummy_precision())
const
383 return fabs(
static_cast<Scalar
>(mass() - other.mass())) <= prec
384 && lever().isApprox(other.lever(),prec)
385 && inertia().isApprox(other.inertia(),prec);
388 bool isZero_impl(
const Scalar & prec = Eigen::NumTraits<Scalar>::dummy_precision())
const
391 return fabs(mass()) <= prec
392 && lever().isZero(prec)
393 && inertia().isZero(prec);
396 InertiaTpl __plus__(
const InertiaTpl & Yb)
const
403 const Scalar eps = ::Eigen::NumTraits<Scalar>::epsilon();
405 const Scalar & mab = mass()+Yb.mass();
406 const Scalar mab_inv = Scalar(1)/math::max((Scalar)(mass()+Yb.mass()),eps);
407 const Vector3 & AB = (lever()-Yb.lever()).eval();
408 return InertiaTpl(mab,
409 (mass()*lever()+Yb.mass()*Yb.lever())*mab_inv,
410 inertia()+Yb.inertia() - (mass()*Yb.mass()*mab_inv)*
typename Symmetric3::SkewSquare(AB));
413 InertiaTpl& __pequ__(
const InertiaTpl & Yb)
415 const Scalar eps = ::Eigen::NumTraits<Scalar>::epsilon();
417 const InertiaTpl& Ya = *
this;
418 const Scalar & mab = mass()+Yb.mass();
419 const Scalar mab_inv = Scalar(1)/math::max((Scalar)(mass()+Yb.mass()),eps);
420 const Vector3 & AB = (Ya.lever()-Yb.lever()).eval();
421 lever() *= (mass()*mab_inv); lever() += (Yb.mass()*mab_inv)*Yb.lever();
422 inertia() += Yb.inertia(); inertia() -= (Ya.mass()*Yb.mass()*mab_inv)*
typename Symmetric3::SkewSquare(AB);
427 template<
typename MotionDerived>
428 ForceTpl<typename traits<MotionDerived>::Scalar,traits<MotionDerived>::Options>
429 __mult__(
const MotionDense<MotionDerived> & v)
const
431 typedef ForceTpl<typename traits<MotionDerived>::Scalar,traits<MotionDerived>::Options> ReturnType;
437 template<
typename MotionDerived,
typename ForceDerived>
438 void __mult__(
const MotionDense<MotionDerived> & v, ForceDense<ForceDerived> & f)
const
440 f.linear().noalias() = mass()*(v.linear() - lever().cross(v.angular()));
441 Symmetric3::rhsMult(inertia(),v.angular(),f.angular());
442 f.angular() += lever().cross(f.linear());
446 Scalar vtiv_impl(
const Motion & v)
const
448 const Vector3 cxw (lever().
cross(v.angular()));
449 Scalar res = mass() * (v.linear().squaredNorm() - Scalar(2)*v.linear().dot(cxw));
450 const Vector3 mcxcxw (-mass()*lever().
cross(cxw));
451 res += v.angular().dot(mcxcxw);
452 res += inertia().vtiv(v.angular());
457 Matrix6 variation(
const Motion & v)
const
460 const Motion mv(v*mass());
462 res.template block<3,3>(LINEAR,ANGULAR) = -
skew(mv.linear()) -
skewSquare(mv.angular(),lever()) +
skewSquare(lever(),mv.angular());
463 res.template block<3,3>(ANGULAR,LINEAR) = res.template block<3,3>(LINEAR,ANGULAR).transpose();
467 res.template block<3,3>(ANGULAR,ANGULAR) = -
skewSquare(mv.linear(),lever()) -
skewSquare(lever(),mv.linear());
469 res.template block<3,3>(LINEAR,LINEAR) = (inertia() - AlphaSkewSquare(mass(),lever())).matrix();
471 res.template block<3,3>(ANGULAR,ANGULAR) -= res.template block<3,3>(LINEAR,LINEAR) *
skew(v.angular());
472 res.template block<3,3>(ANGULAR,ANGULAR) +=
cross(v.angular(),res.template block<3,3>(LINEAR,LINEAR));
474 res.template block<3,3>(LINEAR,LINEAR).setZero();
478 template<
typename M6>
479 static void vxi_impl(
const Motion & v,
480 const InertiaTpl & I,
481 const Eigen::MatrixBase<M6> & Iout)
483 EIGEN_STATIC_ASSERT_MATRIX_SPECIFIC_SIZE(M6,6,6);
484 M6 & Iout_ = PINOCCHIO_EIGEN_CONST_CAST(M6,Iout);
487 alphaSkew(I.mass(),v.angular(),Iout_.template block<3,3>(LINEAR,LINEAR));
489 const Vector3 mc(I.mass()*I.lever());
492 skewSquare(-v.angular(),mc,Iout_.template block<3,3>(LINEAR,ANGULAR));
496 alphaSkew(I.mass(),v.linear(),Iout_.template block<3,3>(ANGULAR,LINEAR));
497 Iout_.template block<3,3>(ANGULAR,LINEAR) -= Iout_.template block<3,3>(LINEAR,ANGULAR);
500 skewSquare(-v.linear(),mc,Iout_.template block<3,3>(ANGULAR,ANGULAR));
507 Symmetric3 mcxcx(
typename Symmetric3::AlphaSkewSquare(I.mass(),I.lever()));
508 Iout_.template block<3,3>(ANGULAR,ANGULAR) += I.inertia().vxs(v.angular());
509 Iout_.template block<3,3>(ANGULAR,ANGULAR) -= mcxcx.vxs(v.angular());
513 template<
typename M6>
514 static void ivx_impl(
const Motion & v,
515 const InertiaTpl & I,
516 const Eigen::MatrixBase<M6> & Iout)
518 EIGEN_STATIC_ASSERT_MATRIX_SPECIFIC_SIZE(M6,6,6);
519 M6 & Iout_ = PINOCCHIO_EIGEN_CONST_CAST(M6,Iout);
522 alphaSkew(I.mass(),v.angular(),Iout_.template block<3,3>(LINEAR,LINEAR));
525 const Vector3 mc(I.mass()*I.lever());
526 skewSquare(mc,v.angular(),Iout_.template block<3,3>(ANGULAR,LINEAR));
529 alphaSkew(I.mass(),v.linear(),Iout_.template block<3,3>(LINEAR,ANGULAR));
532 cross(-I.lever(),Iout_.template block<3,3>(ANGULAR,LINEAR),Iout_.template block<3,3>(ANGULAR,ANGULAR));
533 Iout_.template block<3,3>(ANGULAR,ANGULAR) += I.inertia().svx(v.angular());
534 for(
int k = 0; k < 3; ++k)
535 Iout_.template block<3,3>(ANGULAR,ANGULAR).col(k) += I.lever().cross(Iout_.template block<3,3>(LINEAR,ANGULAR).col(k));
538 Iout_.template block<3,3>(LINEAR,ANGULAR) -= Iout_.template block<3,3>(ANGULAR,LINEAR);
543 Scalar mass()
const {
return m_mass; }
544 const Vector3 & lever()
const {
return m_com; }
545 const Symmetric3 & inertia()
const {
return m_inertia; }
547 Scalar & mass() {
return m_mass; }
548 Vector3 & lever() {
return m_com; }
549 Symmetric3 & inertia() {
return m_inertia; }
558 M.translation()+M.rotation()*lever(),
559 inertia().rotate(M.rotation()));
566 M.rotation().transpose()*(lever()-M.translation()),
567 inertia().rotate(M.rotation().transpose()) );
572 const Vector3 & mcxw = mass()*lever().cross(v.angular());
573 const Vector3 & mv_mcxw = mass()*v.linear()-mcxw;
574 return Force( v.angular().cross(mv_mcxw),
575 v.angular().cross(lever().
cross(mv_mcxw)+inertia()*v.angular())-v.linear().cross(mcxw) );
578 void disp_impl(std::ostream & os)
const
581 <<
" m = " << mass() <<
"\n"
582 <<
" c = " << lever().transpose() <<
"\n"
583 <<
" I = \n" << inertia().matrix() <<
"";
587 template<
typename NewScalar>
591 lever().
template cast<NewScalar>(),
592 inertia().
template cast<NewScalar>());
613 #endif // ifndef __pinocchio_inertia_hpp__