6 #ifndef __pinocchio_spatial_explog_hpp__ 
    7 #define __pinocchio_spatial_explog_hpp__ 
    9 #include "pinocchio/fwd.hpp" 
   10 #include "pinocchio/utils/static-if.hpp" 
   11 #include "pinocchio/math/fwd.hpp" 
   12 #include "pinocchio/math/sincos.hpp" 
   13 #include "pinocchio/math/taylor-expansion.hpp" 
   14 #include "pinocchio/spatial/motion.hpp" 
   15 #include "pinocchio/spatial/skew.hpp" 
   16 #include "pinocchio/spatial/se3.hpp" 
   18 #include <Eigen/Geometry> 
   20 #include "pinocchio/spatial/log.hpp" 
   32   template<
typename Vector3Like>
 
   33   typename Eigen::Matrix<
typename Vector3Like::Scalar,3,3,PINOCCHIO_EIGEN_PLAIN_TYPE(Vector3Like)::Options>
 
   34   exp3(const Eigen::MatrixBase<Vector3Like> & v)
 
   36     PINOCCHIO_ASSERT_MATRIX_SPECIFIC_SIZE (Vector3Like, v, 3, 1);
 
   38     typedef typename Vector3Like::Scalar Scalar;
 
   39     typedef typename PINOCCHIO_EIGEN_PLAIN_TYPE(Vector3Like) Vector3LikePlain;
 
   40     typedef Eigen::Matrix<Scalar,3,3,Vector3LikePlain::Options> Matrix3;
 
   42     const Scalar t2 = v.squaredNorm();
 
   44     const Scalar t = math::sqrt(t2);
 
   45     Scalar ct,st; 
SINCOS(t,&st,&ct);
 
   49                                                      Scalar(1)/Scalar(2) - t2/24);
 
   53     Matrix3 res(alpha_vxvx * v * v.transpose());
 
   54     res.coeffRef(0,1) -= alpha_vx * v[2]; res.coeffRef(1,0) += alpha_vx * v[2];
 
   55     res.coeffRef(0,2) += alpha_vx * v[1]; res.coeffRef(2,0) -= alpha_vx * v[1];
 
   56     res.coeffRef(1,2) -= alpha_vx * v[0]; res.coeffRef(2,1) += alpha_vx * v[0];
 
   61     res.diagonal().array() += ct;   
 
   73   template<
typename Matrix3Like>
 
   74   Eigen::Matrix<
typename Matrix3Like::Scalar,3,1,PINOCCHIO_EIGEN_PLAIN_TYPE(Matrix3Like)::Options>
 
   75   log3(const Eigen::MatrixBase<Matrix3Like> & R,
 
   76        typename Matrix3Like::Scalar & theta)
 
   78     typedef typename Matrix3Like::Scalar Scalar;
 
   79     typedef Eigen::Matrix<Scalar,3,1,
 
   80                           PINOCCHIO_EIGEN_PLAIN_TYPE(Matrix3Like)::Options> Vector3;
 
   94   template<
typename Matrix3Like>
 
   95   Eigen::Matrix<
typename Matrix3Like::Scalar,3,1,PINOCCHIO_EIGEN_PLAIN_TYPE(Matrix3Like)::Options>
 
   96   log3(const Eigen::MatrixBase<Matrix3Like> & R)
 
   98     typename Matrix3Like::Scalar theta;
 
   99     return log3(R.derived(),theta);
 
  110   template<AssignmentOperatorType op, 
typename Vector3Like, 
typename Matrix3Like>
 
  111   void Jexp3(
const Eigen::MatrixBase<Vector3Like> & r,
 
  112              const Eigen::MatrixBase<Matrix3Like> & Jexp)
 
  114     PINOCCHIO_ASSERT_MATRIX_SPECIFIC_SIZE (Vector3Like, r   , 3, 1);
 
  115     PINOCCHIO_ASSERT_MATRIX_SPECIFIC_SIZE (Matrix3Like, Jexp, 3, 3);
 
  117     Matrix3Like & Jout = PINOCCHIO_EIGEN_CONST_CAST(Matrix3Like,Jexp);
 
  118     typedef typename Matrix3Like::Scalar Scalar;
 
  120     const Scalar n2 = r.squaredNorm();
 
  121     const Scalar n = math::sqrt(n2);
 
  122     const Scalar n_inv = Scalar(1)/n;
 
  123     const Scalar n2_inv = n_inv * n_inv;
 
  124     Scalar cn,sn; 
SINCOS(n,&sn,&cn);
 
  127                                             Scalar(1) - n2/Scalar(6),
 
  130                                             - Scalar(1)/Scalar(2) - n2/Scalar(24),
 
  133                                             Scalar(1)/Scalar(6) - n2/Scalar(120),
 
  139         Jout.diagonal().setConstant(a);
 
  140         Jout(0,1) = -b*r[2]; Jout(1,0) = -Jout(0,1);
 
  141         Jout(0,2) =  b*r[1]; Jout(2,0) = -Jout(0,2);
 
  142         Jout(1,2) = -b*r[0]; Jout(2,1) = -Jout(1,2); 
 
  143         Jout.noalias() += c * r * r.transpose();
 
  146         Jout.diagonal().array() += a;
 
  147         Jout(0,1) += -b*r[2]; Jout(1,0) += b*r[2];
 
  148         Jout(0,2) +=  b*r[1]; Jout(2,0) += -b*r[1];
 
  149         Jout(1,2) += -b*r[0]; Jout(2,1) += b*r[0]; 
 
  150         Jout.noalias() += c * r * r.transpose();
 
  153         Jout.diagonal().array() -= a;
 
  154         Jout(0,1) -= -b*r[2]; Jout(1,0) -= b*r[2];
 
  155         Jout(0,2) -=  b*r[1]; Jout(2,0) -= -b*r[1];
 
  156         Jout(1,2) -= -b*r[0]; Jout(2,1) -= b*r[0]; 
 
  157         Jout.noalias() -= c * r * r.transpose();
 
  160         assert(
false && 
"Wrong Op requesed value");
 
  173   template<
typename Vector3Like, 
typename Matrix3Like>
 
  174   void Jexp3(
const Eigen::MatrixBase<Vector3Like> & r,
 
  175              const Eigen::MatrixBase<Matrix3Like> & Jexp)
 
  177     Jexp3<SETTO>(r, Jexp);
 
  222   template<
typename Scalar, 
typename Vector3Like, 
typename Matrix3Like>
 
  224              const Eigen::MatrixBase<Vector3Like> & log,
 
  225              const Eigen::MatrixBase<Matrix3Like> & Jlog)
 
  228                             PINOCCHIO_EIGEN_CONST_CAST(Matrix3Like,Jlog));
 
  243   template<
typename Matrix3Like1, 
typename Matrix3Like2>
 
  244   void Jlog3(
const Eigen::MatrixBase<Matrix3Like1> & R,
 
  245              const Eigen::MatrixBase<Matrix3Like2> & Jlog)
 
  247     typedef typename Matrix3Like1::Scalar Scalar;
 
  248     typedef Eigen::Matrix<Scalar,3,1,PINOCCHIO_EIGEN_PLAIN_TYPE(Matrix3Like1)::Options> Vector3;
 
  251     Vector3 w(
log3(R,t));
 
  252     Jlog3(t,w,PINOCCHIO_EIGEN_CONST_CAST(Matrix3Like2,Jlog));
 
  255   template<
typename Scalar, 
typename Vector3Like1, 
typename Vector3Like2, 
typename Matrix3Like>
 
  256   void Hlog3(
const Scalar & theta,
 
  257              const Eigen::MatrixBase<Vector3Like1> & log,
 
  258              const Eigen::MatrixBase<Vector3Like2> & v,
 
  259              const Eigen::MatrixBase<Matrix3Like> & vt_Hlog)
 
  261     typedef Eigen::Matrix<Scalar,3,1,PINOCCHIO_EIGEN_PLAIN_TYPE(Matrix3Like)::Options> Vector3;
 
  262     Matrix3Like & vt_Hlog_ = PINOCCHIO_EIGEN_CONST_CAST(Matrix3Like,vt_Hlog);
 
  268     Scalar ctheta,stheta; 
SINCOS(theta,&stheta,&ctheta);
 
  270     Scalar denom = .5 / (1-ctheta),
 
  271            a = theta * stheta * denom,
 
  272            da_dt = (stheta - theta) * denom, 
 
  273            b = ( 1 - a ) / (theta*theta),
 
  275            db_dt = - (2 / theta - (theta + stheta) * denom) / (theta*theta); 
 
  280     Vector3 dl_dv_v (a*v + .5*log.cross(v) + b*log*log.transpose()*v);
 
  282     Scalar dt_dv_v = log.dot(dl_dv_v) / theta;
 
  285     vt_Hlog_.noalias() = db_dt * dt_dv_v * log * log.transpose();
 
  286     vt_Hlog_.noalias() += b * dl_dv_v * log.transpose();
 
  287     vt_Hlog_.noalias() += b * log * dl_dv_v.transpose();
 
  289     addSkew(.5 * dl_dv_v, vt_Hlog_);
 
  291     vt_Hlog_.diagonal().array() += da_dt * dt_dv_v;
 
  302   template<
typename Matrix3Like1, 
typename Vector3Like, 
typename Matrix3Like2>
 
  303   void Hlog3(
const Eigen::MatrixBase<Matrix3Like1> & R,
 
  304              const Eigen::MatrixBase<Vector3Like> & v,
 
  305              const Eigen::MatrixBase<Matrix3Like2> & vt_Hlog)
 
  307     typedef typename Matrix3Like1::Scalar Scalar;
 
  308     typedef Eigen::Matrix<Scalar,3,1,PINOCCHIO_EIGEN_PLAIN_TYPE(Matrix3Like1)::Options> Vector3;
 
  311     Vector3 w(
log3(R,t));
 
  312     Hlog3(t,w,v,PINOCCHIO_EIGEN_CONST_CAST(Matrix3Like2,vt_Hlog));
 
  324   template<
typename MotionDerived>
 
  325   SE3Tpl<
typename MotionDerived::Scalar,PINOCCHIO_EIGEN_PLAIN_TYPE(
typename MotionDerived::Vector3)::Options>
 
  328     typedef typename MotionDerived::Scalar Scalar;
 
  329     enum { Options = PINOCCHIO_EIGEN_PLAIN_TYPE(
typename MotionDerived::Vector3)::Options };
 
  334     typename SE3::LinearType & trans = res.translation();
 
  335     typename SE3::AngularType & rot = res.rotation();
 
  337     const typename MotionDerived::ConstAngularType & w = nu.angular();
 
  338     const typename MotionDerived::ConstLinearType & v = nu.linear();
 
  340     Scalar alpha_wxv, alpha_v, alpha_w, diagonal_term;
 
  341     const Scalar t2 = w.squaredNorm();
 
  342     const Scalar t = math::sqrt(t2);
 
  343     Scalar ct,st; 
SINCOS(t,&st,&ct);
 
  344     const Scalar inv_t2 = Scalar(1)/t2;
 
  347                                        Scalar(1)/Scalar(2) - t2/24,
 
  348                                        (Scalar(1) - ct)*inv_t2);
 
  355                                      (Scalar(1)/Scalar(6) - t2/120),
 
  356                                      (Scalar(1) - alpha_v)*inv_t2);
 
  363     trans.noalias() = (alpha_v*v + (alpha_w*w.dot(v))*w + alpha_wxv*w.cross(v));
 
  366     rot.noalias() = alpha_wxv * w * w.transpose();
 
  367     rot.coeffRef(0,1) -= alpha_v * w[2]; rot.coeffRef(1,0) += alpha_v * w[2];
 
  368     rot.coeffRef(0,2) += alpha_v * w[1]; rot.coeffRef(2,0) -= alpha_v * w[1];
 
  369     rot.coeffRef(1,2) -= alpha_v * w[0]; rot.coeffRef(2,1) += alpha_v * w[0];
 
  370     rot.diagonal().array() += diagonal_term;
 
  383   template<
typename Vector6Like>
 
  384   SE3Tpl<
typename Vector6Like::Scalar,PINOCCHIO_EIGEN_PLAIN_TYPE(Vector6Like)::Options>
 
  385   exp6(const Eigen::MatrixBase<Vector6Like> & v)
 
  387     PINOCCHIO_ASSERT_MATRIX_SPECIFIC_SIZE (Vector6Like, v, 6, 1);
 
  401   template<
typename Scalar, 
int Options>
 
  402   MotionTpl<Scalar,Options>
 
  419   template<
typename Matrix4Like>
 
  420   MotionTpl<typename Matrix4Like::Scalar,Eigen::internal::traits<Matrix4Like>::Options>
 
  421   log6(
const Eigen::MatrixBase<Matrix4Like> & M)
 
  423     PINOCCHIO_ASSERT_MATRIX_SPECIFIC_SIZE(Matrix4Like, M, 4, 4);
 
  425     typedef typename Matrix4Like::Scalar Scalar;
 
  426     enum {Options = Eigen::internal::traits<Matrix4Like>::Options};
 
  438   template<AssignmentOperatorType op, 
typename MotionDerived, 
typename Matrix6Like>
 
  440              const Eigen::MatrixBase<Matrix6Like> & Jexp)
 
  442     PINOCCHIO_ASSERT_MATRIX_SPECIFIC_SIZE (Matrix6Like, Jexp, 6, 6);
 
  444     typedef typename MotionDerived::Scalar Scalar;
 
  445     typedef typename MotionDerived::Vector3 Vector3;
 
  446     typedef Eigen::Matrix<Scalar, 3, 3, Vector3::Options> Matrix3;
 
  447     Matrix6Like & Jout = PINOCCHIO_EIGEN_CONST_CAST(Matrix6Like,Jexp);
 
  449     const typename MotionDerived::ConstLinearType  & v = nu.linear();
 
  450     const typename MotionDerived::ConstAngularType & w = nu.angular();
 
  451     const Scalar t2 = w.squaredNorm();
 
  452     const Scalar t = math::sqrt(t2);
 
  454     const Scalar tinv = Scalar(1)/t,
 
  456     Scalar st,ct; 
SINCOS (t, &st, &ct);
 
  457     const Scalar inv_2_2ct = Scalar(1)/(Scalar(2)*(Scalar(1)-ct));
 
  461                                                Scalar(1)/Scalar(12) + t2/Scalar(720),
 
  462                                                t2inv - st*tinv*inv_2_2ct);
 
  465                                                               Scalar(1)/Scalar(360),
 
  466                                                               -Scalar(2)*t2inv*t2inv + (Scalar(1) + st*tinv) * t2inv * inv_2_2ct);
 
  472         Jexp3<SETTO>(w, Jout.template bottomRightCorner<3,3>());
 
  473         Jout.template topLeftCorner<3,3>() = Jout.template bottomRightCorner<3,3>();
 
  474         const Vector3 p = Jout.template topLeftCorner<3,3>().transpose() * v;
 
  475         const Scalar wTp (w.dot (p));
 
  477                          (beta_dot_over_theta*wTp)                *w*w.transpose()
 
  478                          - (t2*beta_dot_over_theta+Scalar(2)*beta)*p*w.transpose()
 
  479                          + wTp * beta                             * Matrix3::Identity()
 
  480                          + beta                                   *w*p.transpose());
 
  481         Jout.template topRightCorner<3,3>().noalias() =
 
  482           - Jout.template topLeftCorner<3,3>() * J;
 
  483         Jout.template bottomLeftCorner<3,3>().setZero();
 
  489         Jexp3<SETTO>(w, Jtmp3);
 
  490         Jout.template bottomRightCorner<3,3>() += Jtmp3;
 
  491         Jout.template topLeftCorner<3,3>() += Jtmp3;
 
  492         const Vector3 p = Jtmp3.transpose() * v;
 
  493         const Scalar wTp (w.dot (p));
 
  495                          (beta_dot_over_theta*wTp)                *w*w.transpose()
 
  496                          - (t2*beta_dot_over_theta+Scalar(2)*beta)*p*w.transpose()
 
  497                          + wTp * beta                             * Matrix3::Identity()
 
  498                          + beta                                   *w*p.transpose());
 
  499         Jout.template topRightCorner<3,3>().noalias() +=
 
  506         Jexp3<SETTO>(w, Jtmp3);
 
  507         Jout.template bottomRightCorner<3,3>() -= Jtmp3;
 
  508         Jout.template topLeftCorner<3,3>() -= Jtmp3;
 
  509         const Vector3 p = Jtmp3.transpose() * v;
 
  510         const Scalar wTp (w.dot (p));
 
  512                          (beta_dot_over_theta*wTp)                *w*w.transpose()
 
  513                          - (t2*beta_dot_over_theta+Scalar(2)*beta)*p*w.transpose()
 
  514                          + wTp * beta                             * Matrix3::Identity()
 
  515                          + beta                                   *w*p.transpose());
 
  516         Jout.template topRightCorner<3,3>().noalias() -=
 
  521         assert(
false && 
"Wrong Op requesed value");
 
  528   template<
typename MotionDerived, 
typename Matrix6Like>
 
  530              const Eigen::MatrixBase<Matrix6Like> & Jexp)
 
  532     Jexp6<SETTO>(nu, Jexp);
 
  593   template<
typename Scalar, 
int Options, 
typename Matrix6Like>
 
  595              const Eigen::MatrixBase<Matrix6Like> & Jlog)
 
  600   template<
typename Scalar, 
int Options>
 
  601   template<
typename OtherScalar>
 
  604                                                              const OtherScalar & alpha)
 
  606     typedef SE3Tpl<Scalar,Options> ReturnType;
 
  607     typedef MotionTpl<Scalar,Options> Motion;
 
  609     Motion dv = log6(A.actInv(B));
 
  610     ReturnType res = A * exp6(alpha*dv);
 
  616 #include "pinocchio/spatial/explog-quaternion.hpp" 
  617 #include "pinocchio/spatial/log.hxx" 
  619 #endif //#ifndef __pinocchio_spatial_explog_hpp__