code: def quaternion_from_matrix(matrix):
      
      q = numpy.empty((4, ), dtype=numpy.float64)
      M = numpy.array(matrix, dtype=numpy.float64, copy=False)[:4, :4]
      t = numpy.trace(M)
      if t > M[3, 3]:
          q[3] = t
          q[2] = M[1, 0] - M[0, 1]
          q[1] = M[0, 2] - M[2, 0]
          q[0] = M[2, 1] - M[1, 2]
      else:
          i, j, k = 0, 1, 2
          if M[1, 1] > M[0, 0]:
              i, j, k = 1, 2, 0
          if M[2, 2] > M[i, i]:
              i, j, k = 2, 0, 1
          t = M[i, i] - (M[j, j] + M[k, k]) + M[3, 3]
          q[i] = t
          q[j] = M[i, j] + M[j, i]
          q[k] = M[k, i] + M[i, k]
          q[3] = M[k, j] - M[j, k]
      q *= 0.5 / math.sqrt(t * M[3, 3])
      return q
docs: Return quaternion from rotation matrix.
  
      >>> R = rotation_matrix(0.123, (1, 2, 3))
      >>> q = quaternion_from_matrix(R)
      >>> numpy.allclose(q, [0.0164262, 0.0328524, 0.0492786, 0.9981095])
      True