code: def mat2quat(M):
      """Calculate quaternion corresponding to given rotation matrix
  
      Parameters
      ----------
      M : array-like
        3x3 rotation matrix
  
      Returns
      -------
      q : (4,) array
        closest quaternion to input matrix, having positive q[0]
  
      Notes
      -----
      Method claimed to be robust to numerical errors in M
  
      Constructs quaternion by calculating maximum eigenvector for matrix
      K (constructed from input `M`).  Although this is not tested, a
      maximum eigenvalue of 1 corresponds to a valid rotation.
  
      A quaternion q*-1 corresponds to the same rotation as q; thus the
      sign of the reconstructed quaternion is arbitrary, and we return
      quaternions with positive w (q[0]).
  
      References
      ----------
      * https://en.wikipedia.org/wiki/Rotation_matrix#Quaternion
      * Bar-Itzhack, Itzhack Y. (2000), "New method for extracting the
        quaternion from a rotation matrix", AIAA Journal of Guidance,
        Control and Dynamics 23(6):1085-1087 (Engineering Note), ISSN
        0731-5090
  
      Examples
      --------
      >>> import numpy as np
      >>> q = mat2quat(np.eye(3)) # Identity rotation
      >>> np.allclose(q, [1, 0, 0, 0])
      True
      >>> q = mat2quat(np.diag([1, -1, -1]))
      >>> np.allclose(q, [0, 1, 0, 0]) # 180 degree rotn around axis 0
      True
  
      """
      # Qyx refers to the contribution of the y input vector component to
      # the x output vector component.  Qyx is therefore the same as
      # M[0,1].  The notation is from the Wikipedia article.
      Qxx, Qyx, Qzx, Qxy, Qyy, Qzy, Qxz, Qyz, Qzz = M.flat
      # Fill only lower half of symmetric matrix
      K = (
          np.array(
              [
                  [Qxx - Qyy - Qzz, 0, 0, 0],
                  [Qyx + Qxy, Qyy - Qxx - Qzz, 0, 0],
                  [Qzx + Qxz, Qzy + Qyz, Qzz - Qxx - Qyy, 0],
                  [Qyz - Qzy, Qzx - Qxz, Qxy - Qyx, Qxx + Qyy + Qzz],
              ]
          )
          / 3.0
      )
      # Use Hermitian eigenvectors, values for speed
      vals, vecs = np.linalg.eigh(K)
      # Select largest eigenvector, reorder to w,x,y,z quaternion
      q = vecs[[3, 0, 1, 2], np.argmax(vals)]
      # Prefer quaternion with positive w
      # (q * -1 corresponds to same rotation as q)
      if q[0] < 0:
          q *= -1
      return q
  
file_path: nibabel/nibabel/quaternions.py
imports: import math
  import numpy
  from .casting import sctypes
label_desc: Use image node_id 3 for calling a global function with example usage: mat2quat(M) and returns: q