Eigen3d.cpp
Go to the documentation of this file.
1 /*
2  * Copyright (c) 2008, AIST, the University of Tokyo and General Robotix Inc.
3  * All rights reserved. This program is made available under the terms of the
4  * Eclipse Public License v1.0 which accompanies this distribution, and is
5  * available at http://www.eclipse.org/legal/epl-v10.html
6  * Contributors:
7  * National Institute of Advanced Industrial Science and Technology (AIST)
8  * General Robotix Inc.
9  */
10 
11 #include <iostream>
12 #include <iomanip>
13 #include "Eigen3d.h"
14 
15 using namespace std;
16 using namespace hrp;
17 
18 
19 static const double PI = 3.14159265358979323846;
20 
21 
22 void hrp::calcRodrigues(Matrix33& out_R, const Vector3& axis, double q)
23 {
24  // E + a_hat*sin(q) + a_hat*a_hat*(1-cos(q))
25  //
26  // | 0 -az ay|
27  // =E+| az 0 -ax|*s + a_hat*a_hat*v
28  // |-ay ax 0|
29  //
30  // | 0 -az ay| |-az*az-ay*ay ax*ay az*ax|
31  // =E+| az 0 -ax|*s + | ax*ay -az*az-ax*ax ay*az|*v
32  // |-ay ax 0| | az*ax ay*az -ax*ax-ay*ay|
33  //
34  // |1-az*az*v-ay*ay*v -az*s+ax*ay*v ay*s+az*ax*v|
35  // =| az*s+ax*ay*v 1-az*az*v-ax*ax*v -ax*s+ay+az*v|
36  // | -ay*s+az*ax*v ax*s+ay*az*v 1-ax*ax*v-ay*ay*v|
37  //
38 
39  const double sth = sin(q);
40  const double vth = 1.0 - cos(q);
41 
42  double ax = axis(0);
43  double ay = axis(1);
44  double az = axis(2);
45 
46  const double axx = ax*ax*vth;
47  const double ayy = ay*ay*vth;
48  const double azz = az*az*vth;
49  const double axy = ax*ay*vth;
50  const double ayz = ay*az*vth;
51  const double azx = az*ax*vth;
52 
53  ax *= sth;
54  ay *= sth;
55  az *= sth;
56 
57  out_R << 1.0 - azz - ayy, -az + axy, ay + azx,
58  az + axy, 1.0 - azz - axx, -ax + ayz,
59  -ay + azx, ax + ayz, 1.0 - ayy - axx;
60 }
61 
62 
64 {
65  using ::std::numeric_limits;
66 
67  double alpha = (r(0,0) + r(1,1) + r(2,2) - 1.0) / 2.0;
68 
69  if (alpha > 1.0) {
70  if (alpha > 1.0 + 1.0e-6) {
71  cerr << scientific << setprecision(16);
72  cerr << "alpha exceeded the upper limit in omegaFromRot" << endl;
73  cerr << "alpha=" << alpha << endl;
74  cerr << "r=" << r << endl;
75  }
76  alpha = 1.0;
77  }
78 
79  if(fabs(alpha - 1.0) < 1.0e-12) { //th=0,2PI;
80  return Vector3::Zero();
81 
82  } else {
83  double th = acos(alpha);
84  double s = sin(th);
85 
86  if (s < numeric_limits<double>::epsilon()) { //th=PI
87  return Vector3( sqrt((r(0,0)+1)*0.5)*th, sqrt((r(1,1)+1)*0.5)*th, sqrt((r(2,2)+1)*0.5)*th );
88  }
89 
90  double k = -0.5 * th / s;
91 
92  return Vector3( (r(1,2) - r(2,1)) * k,
93  (r(2,0) - r(0,2)) * k,
94  (r(0,1) - r(1,0)) * k );
95  }
96 }
97 
98 
100 {
101  double roll, pitch, yaw;
102 
103  if ((fabs(m(0,0))<fabs(m(2,0))) && (fabs(m(1,0))<fabs(m(2,0)))) {
104  // cos(p) is nearly = 0
105  double sp = -m(2,0);
106  if (sp < -1.0) {
107  sp = -1;
108  } else if (sp > 1.0) {
109  sp = 1;
110  }
111  pitch = asin(sp); // -pi/2< p < pi/2
112 
113  roll = atan2(sp*m(0,1)+m(1,2), // -cp*cp*sr*cy
114  sp*m(0,2)-m(1,1)); // -cp*cp*cr*cy
115 
116  if (m(0,0)>0.0) { // cy > 0
117  (roll < 0.0) ? (roll += PI) : (roll -= PI);
118  }
119  double sr=sin(roll), cr=cos(roll);
120  if (sp > 0.0) {
121  yaw = atan2(sr*m(1,1)+cr*m(1,2), //sy*sp
122  sr*m(0,1)+cr*m(0,2));//cy*sp
123  } else {
124  yaw = atan2(-sr*m(1,1)-cr*m(1,2),
125  -sr*m(0,1)-cr*m(0,2));
126  }
127  } else {
128  yaw = atan2(m(1,0), m(0,0));
129  const double sa = sin(yaw);
130  const double ca = cos(yaw);
131  pitch = atan2(-m(2,0), ca*m(0,0)+sa*m(1,0));
132  roll = atan2(sa*m(0,2)-ca*m(1,2), -sa*m(0,1)+ca*m(1,1));
133  }
134  return Vector3(roll, pitch, yaw);
135 }
136 
137 
138 void hrp::calcRotFromRpy(Matrix33& out_R, double r, double p, double y)
139 {
140  const double cr = cos(r), sr = sin(r), cp = cos(p), sp = sin(p), cy = cos(y), sy = sin(y);
141  out_R(0,0)= cp*cy;
142  out_R(0,1)= sr*sp*cy - cr*sy;
143  out_R(0,2)= cr*sp*cy + sr*sy;
144  out_R(1,0)= cp*sy;
145  out_R(1,1)= sr*sp*sy + cr*cy;
146  out_R(1,2)= cr*sp*sy - sr*cy;
147  out_R(2,0)= -sp;
148  out_R(2,1)= sr*cp;
149  out_R(2,2)= cr*cp;
150 }
151 
153  Matrix33 w(m * m.transpose() - Matrix33::Identity());
154  return (abs(w.array())<1.0e-12).all();
155  //return all_elements( m * m.transpose() == Matrix33::Identity() );
156 }
157 
HRP_UTIL_EXPORT bool isOrthogonalMatrix(Matrix33 &m)
Definition: Eigen3d.cpp:152
HRP_UTIL_EXPORT void calcRotFromRpy(Matrix33 &out_R, double r, double p, double y)
Definition: Eigen3d.cpp:138
w
Eigen::Vector3d Vector3
Definition: EigenTypes.h:11
Eigen::Matrix3d Matrix33
Definition: EigenTypes.h:12
OpenHRP::vector3 Vector3
static const double PI
Definition: Eigen3d.cpp:19
HRP_UTIL_EXPORT Vector3 rpyFromRot(const Matrix33 &m)
Definition: Eigen3d.cpp:99
HRP_UTIL_EXPORT void calcRodrigues(Matrix33 &out_R, const Vector3 &axis, double q)
Definition: Eigen3d.cpp:22
HRP_UTIL_EXPORT Vector3 omegaFromRot(const Matrix33 &r)
Definition: Eigen3d.cpp:63
* y
Definition: IceUtils.h:97


openhrp3
Author(s): AIST, General Robotix Inc., Nakamura Lab of Dept. of Mechano Informatics at University of Tokyo
autogenerated on Thu Sep 8 2022 02:24:03