Quat.cpp 5.4 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277
  1. #include "Quat.hpp"
  2. #include <cmath>
  3. #include "Defs.hpp"
  4. #include "Vector3.hpp"
  5. #include "Basis.hpp"
  6. namespace godot {
  7. real_t Quat::length() const
  8. {
  9. return ::sqrt(length_squared());
  10. }
  11. void Quat::normalize()
  12. {
  13. *this /= length();
  14. }
  15. Quat Quat::normalized() const
  16. {
  17. return *this / length();
  18. }
  19. Quat Quat::inverse() const
  20. {
  21. return Quat( -x, -y, -z, w );
  22. }
  23. void Quat::set_euler(const Vector3& p_euler)
  24. {
  25. real_t half_a1 = p_euler.x * 0.5;
  26. real_t half_a2 = p_euler.y * 0.5;
  27. real_t half_a3 = p_euler.z * 0.5;
  28. // R = X(a1).Y(a2).Z(a3) convention for Euler angles.
  29. // Conversion to quaternion as listed in https://ntrs.nasa.gov/archive/nasa/casi.ntrs.nasa.gov/19770024290.pdf (page A-2)
  30. // a3 is the angle of the first rotation, following the notation in this reference.
  31. real_t cos_a1 = ::cos(half_a1);
  32. real_t sin_a1 = ::sin(half_a1);
  33. real_t cos_a2 = ::cos(half_a2);
  34. real_t sin_a2 = ::sin(half_a2);
  35. real_t cos_a3 = ::cos(half_a3);
  36. real_t sin_a3 = ::sin(half_a3);
  37. set(sin_a1*cos_a2*cos_a3 + sin_a2*sin_a3*cos_a1,
  38. -sin_a1*sin_a3*cos_a2 + sin_a2*cos_a1*cos_a3,
  39. sin_a1*sin_a2*cos_a3 + sin_a3*cos_a1*cos_a2,
  40. -sin_a1*sin_a2*sin_a3 + cos_a1*cos_a2*cos_a3);
  41. }
  42. Quat Quat::slerp(const Quat& q, const real_t& t) const {
  43. Quat to1;
  44. real_t omega, cosom, sinom, scale0, scale1;
  45. // calc cosine
  46. cosom = dot(q);
  47. // adjust signs (if necessary)
  48. if ( cosom <0.0 ) {
  49. cosom = -cosom;
  50. to1.x = - q.x;
  51. to1.y = - q.y;
  52. to1.z = - q.z;
  53. to1.w = - q.w;
  54. } else {
  55. to1.x = q.x;
  56. to1.y = q.y;
  57. to1.z = q.z;
  58. to1.w = q.w;
  59. }
  60. // calculate coefficients
  61. if ( (1.0 - cosom) > CMP_EPSILON ) {
  62. // standard case (slerp)
  63. omega = ::acos(cosom);
  64. sinom = ::sin(omega);
  65. scale0 = ::sin((1.0 - t) * omega) / sinom;
  66. scale1 = ::sin(t * omega) / sinom;
  67. } else {
  68. // "from" and "to" quaternions are very close
  69. // ... so we can do a linear interpolation
  70. scale0 = 1.0 - t;
  71. scale1 = t;
  72. }
  73. // calculate final values
  74. return Quat(
  75. scale0 * x + scale1 * to1.x,
  76. scale0 * y + scale1 * to1.y,
  77. scale0 * z + scale1 * to1.z,
  78. scale0 * w + scale1 * to1.w
  79. );
  80. }
  81. Quat Quat::slerpni(const Quat& q, const real_t& t) const {
  82. const Quat &from = *this;
  83. real_t dot = from.dot(q);
  84. if (::fabs(dot) > 0.9999) return from;
  85. real_t theta = ::acos(dot),
  86. sinT = 1.0 / ::sin(theta),
  87. newFactor = ::sin(t * theta) * sinT,
  88. invFactor = ::sin((1.0 - t) * theta) * sinT;
  89. return Quat(invFactor * from.x + newFactor * q.x,
  90. invFactor * from.y + newFactor * q.y,
  91. invFactor * from.z + newFactor * q.z,
  92. invFactor * from.w + newFactor * q.w);
  93. }
  94. Quat Quat::cubic_slerp(const Quat& q, const Quat& prep, const Quat& postq,const real_t& t) const
  95. {
  96. //the only way to do slerp :|
  97. real_t t2 = (1.0-t)*t*2;
  98. Quat sp = this->slerp(q,t);
  99. Quat sq = prep.slerpni(postq,t);
  100. return sp.slerpni(sq,t2);
  101. }
  102. void Quat::get_axis_and_angle(Vector3& r_axis, real_t &r_angle) const {
  103. r_angle = 2 * ::acos(w);
  104. r_axis.x = x / ::sqrt(1-w*w);
  105. r_axis.y = y / ::sqrt(1-w*w);
  106. r_axis.z = z / ::sqrt(1-w*w);
  107. }
  108. Quat Quat::operator*(const Vector3& v) const
  109. {
  110. return Quat( w * v.x + y * v.z - z * v.y,
  111. w * v.y + z * v.x - x * v.z,
  112. w * v.z + x * v.y - y * v.x,
  113. -x * v.x - y * v.y - z * v.z);
  114. }
  115. Vector3 Quat::xform(const Vector3& v) const {
  116. Quat q = *this * v;
  117. q *= this->inverse();
  118. return Vector3(q.x,q.y,q.z);
  119. }
  120. Quat::operator String() const
  121. {
  122. return String(); // @Todo
  123. }
  124. Quat::Quat(const Vector3& axis, const real_t& angle)
  125. {
  126. real_t d = axis.length();
  127. if (d==0)
  128. set(0,0,0,0);
  129. else {
  130. real_t sin_angle = ::sin(angle * 0.5);
  131. real_t cos_angle = ::cos(angle * 0.5);
  132. real_t s = sin_angle / d;
  133. set(axis.x * s, axis.y * s, axis.z * s,
  134. cos_angle);
  135. }
  136. }
  137. Quat::Quat(const Vector3& v0, const Vector3& v1) // shortest arc
  138. {
  139. Vector3 c = v0.cross(v1);
  140. real_t d = v0.dot(v1);
  141. if (d < -1.0 + CMP_EPSILON) {
  142. x=0;
  143. y=1;
  144. z=0;
  145. w=0;
  146. } else {
  147. real_t s = ::sqrt((1.0 + d) * 2.0);
  148. real_t rs = 1.0 / s;
  149. x=c.x*rs;
  150. y=c.y*rs;
  151. z=c.z*rs;
  152. w=s * 0.5;
  153. }
  154. }
  155. real_t Quat::dot(const Quat& q) const {
  156. return x * q.x+y * q.y+z * q.z+w * q.w;
  157. }
  158. real_t Quat::length_squared() const {
  159. return dot(*this);
  160. }
  161. void Quat::operator+=(const Quat& q) {
  162. x += q.x; y += q.y; z += q.z; w += q.w;
  163. }
  164. void Quat::operator-=(const Quat& q) {
  165. x -= q.x; y -= q.y; z -= q.z; w -= q.w;
  166. }
  167. void Quat::operator*=(const Quat& q) {
  168. x *= q.x; y *= q.y; z *= q.z; w *= q.w;
  169. }
  170. void Quat::operator*=(const real_t& s) {
  171. x *= s; y *= s; z *= s; w *= s;
  172. }
  173. void Quat::operator/=(const real_t& s) {
  174. *this *= 1.0 / s;
  175. }
  176. Quat Quat::operator+(const Quat& q2) const {
  177. const Quat& q1 = *this;
  178. return Quat( q1.x+q2.x, q1.y+q2.y, q1.z+q2.z, q1.w+q2.w );
  179. }
  180. Quat Quat::operator-(const Quat& q2) const {
  181. const Quat& q1 = *this;
  182. return Quat( q1.x-q2.x, q1.y-q2.y, q1.z-q2.z, q1.w-q2.w);
  183. }
  184. Quat Quat::operator*(const Quat& q2) const {
  185. Quat q1 = *this;
  186. q1 *= q2;
  187. return q1;
  188. }
  189. Quat Quat::operator-() const {
  190. const Quat& q2 = *this;
  191. return Quat( -q2.x, -q2.y, -q2.z, -q2.w);
  192. }
  193. Quat Quat::operator*(const real_t& s) const {
  194. return Quat(x * s, y * s, z * s, w * s);
  195. }
  196. Quat Quat::operator/(const real_t& s) const {
  197. return *this * (1.0 / s);
  198. }
  199. bool Quat::operator==(const Quat& p_quat) const {
  200. return x==p_quat.x && y==p_quat.y && z==p_quat.z && w==p_quat.w;
  201. }
  202. bool Quat::operator!=(const Quat& p_quat) const {
  203. return x!=p_quat.x || y!=p_quat.y || z!=p_quat.z || w!=p_quat.w;
  204. }
  205. Vector3 Quat::get_euler() const
  206. {
  207. Basis m(*this);
  208. return m.get_euler();
  209. }
  210. }