title: Getting relative angular velocity from quaternions
body: <p>I want to simulate gyroscope readings (angular velocity in local frame) by supplied vectors.
  I have this vectors (x,y,z) with sampling 1s:</p>
  <pre><code>[0, 0, 1000]
  [469.47156279, 0.0, 882.94759286]
  [981.62718345, 0.0, 190.80899538]
  </code></pre>
  <p>Rotating Z by 28 degrees by Y axis and then by 51 degrees again around Y.
  Vectors looks like this:
  <a href="https://i.stack.imgur.com/opNSA.jpg" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/opNSA.jpg" alt="enter image description here" /></a></p>
  <p>I am getting quaternions describing relative rotations using</p>
  <pre><code>scipy.spatial.transform.Rotation.align_vectors
  </code></pre>
  <p>In a form that</p>
  <pre><code>1. vector_orientation = quaternion_1
  2. vector_orientation = quaternion_2 * quaternion_1
  3. vector_orientation = quaternion_3 * quaternion_2 * quaternion_1
  </code></pre>
  <p>I found this <a href="https://stackoverflow.com/questions/56790186/how-to-compute-angular-velocity-using-numpy-quaternion">How to compute angular velocity using numpy-quaternion</a> question, so I tried to implement suggested code for calculating angular velocity between orientation 2 and 3:</p>
  <pre><code># quaternion_3 is relative to quaternion_2, so I guess I don't need quaternion difference
  delta_q = quaternion_3
  delta_q_len = np.linalg.norm(delta_q[1:])
  delta_q_angle = 2*np.arctan2(delta_q_len, delta_q[0])
  w = delta_q[1:] * delta_q_angle * 1
  </code></pre>
  <p>I got some velocities - <code>[0.0, 0.38320564, 0.0]</code>.</p>
  <p>However if I want to check the result by applying it to the quaternion_2 by obtaining derivative using this equation:
  <a href="https://i.stack.imgur.com/H8Qqr.jpg" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/H8Qqr.jpg" alt="enter image description here" /></a></p>
  <p>So in Python:</p>
  <pre><code>omega = np.array([[0, -w[0], -w[1], -w[2]], 
                    [w[0], 0, w[2], -w[1]], 
                    [w[1], -w[2], 0, w[0]], 
                    [w[2], w[1], -w[0], 0]])
  quaternion_derivative = 0.5 * omega @ quaternion_2
  </code></pre>
  <p>And integrate it with</p>
  <pre><code>quaternion_result = (2 * quaternion_derivative / quaternion_2).normalized()
  </code></pre>
  <p>Problem is <code>quaternion_result != quaternion_3</code>. Does anybody has some idea what I am doing wrong?</p>
  <p>Edit:
  Adding (not)working example. Also changed angles for better demonstration. I am using <a href="https://github.com/moble/quaternion" rel="nofollow noreferrer">https://github.com/moble/quaternion</a> for work with quaternions:</p>
  <pre><code>import numpy as np
  from scipy.spatial.transform import Rotation as R
  import quaternion as qt
  
  def get_quaternion_rotation(v1, v2):
      rotation,_ = R.align_vectors([v1], [v2])
      qt = rotation.as_quat()
      return np.array([qt[3], qt[0], qt[1], qt[2]])
  
  def rotate_vector(vector, axis, angle):
      # normalize
      unit_axis = axis/np.linalg.norm(axis)
      angle_radians = np.radians(angle)
      rotation_vector = angle_radians * unit_axis
      rotation = R.from_rotvec(rotation_vector)
      return rotation.apply(vector)
  
  # 1. vector is just 1000 in Z
  vector_1 = np.array([0, 0, 1000])
  # 2. vector rotates vector_1 28 degrees by Y
  vector_2 = rotate_vector(vector_1, [0,1,0], 28)
  # 3. vector rotates vector_2 51 degrees by Y
  vector_3 = rotate_vector(vector_2, [0,1,0], 51)
  
  # rotations between vectors
  quaternion_1 = get_quaternion_rotation([0,0,0], vector_1)
  quaternion_2 = get_quaternion_rotation(vector_2, vector_1)
  quaternion_3 = get_quaternion_rotation(vector_3, vector_2)
  
  delta_q = quaternion_3
  delta_q_len = np.linalg.norm(delta_q[1:])
  delta_q_angle = 2 * np.arctan2(delta_q_len, delta_q[0])
  w = delta_q[1:] * delta_q_angle * 1
  
  omega = np.array([[0, -w[0], -w[1], -w[2]], 
                    [w[0], 0, w[2], -w[1]], 
                    [w[1], -w[2], 0, w[0]], 
                    [w[2], w[1], -w[0], 0]])
  
  
  qt_derivative = 0.5 * omega @ quaternion_2
  # converting to quaternion type for normalizing
  qt_derivative = qt.quaternion(qt_derivative[0], qt_derivative[1], qt_derivative[2], qt_derivative[3]).normalized()
  # apply derivative to quaternion_2
  qt_result = (2 * qt_derivative / qt.quaternion(quaternion_2[0], quaternion_2[1], quaternion_2[2], quaternion_2[3])).normalized()
  
  print(f&quot;Quaternion 3 = {quaternion_3}&quot;)
  print(f&quot;Quaternion result = {qt.as_float_array(qt_result)}&quot;)
  
  
  rotation_quaternion_3 = R.from_quat([quaternion_3[1], quaternion_3[2], quaternion_3[3], quaternion_3[0]])
  print(f&quot;Applied rotation from quaternion_3 to vector_2: {rotation_quaternion_3.apply(vector_2)}&quot;)
  rotation_result = R.from_quat([qt_result.x, qt_result.y, qt_result.z, qt_result.w])
  print(f&quot;Applied rotation from quaternion result to vector_2: {rotation_result.apply(vector_2)}&quot;)
  </code></pre>
  <p>Result is:</p>
  <pre><code>Quaternion 3 = [0.90258528 0.         0.4305111  0.        ]
  Quaternion result = [-2.77555756e-17 -2.76809689e-17 -1.00000000e+00 -1.11022302e-16]
  Applied rotation from quaternion_3 to vector_2: [981.62718345   0.         190.80899538]
  Applied rotation from quaternion result to vector_2: [-4.69471563e+02  2.22044605e-13 -8.82947593e+02]
  </code></pre>