code: def rotation_vec2mat(r):
      """
      R = rotation_vec2mat(r)
  
      The rotation matrix is given by the Rodrigues formula:
  
      R = Id + sin(theta)*Sn + (1-cos(theta))*Sn^2
  
      with:
  
             0  -nz  ny
      Sn =   nz   0 -nx
            -ny  nx   0
  
      where n = r / ||r||
  
      In case the angle ||r|| is very small, the above formula may lead
      to numerical instabilities. We instead use a Taylor expansion
      around theta=0:
  
      R = I + sin(theta)/tetha Sr + (1-cos(theta))/teta2 Sr^2
  
      leading to:
  
      R = I + (1-theta2/6)*Sr + (1/2-theta2/24)*Sr^2
  
      To avoid numerical instabilities, an upper threshold is applied to
      the angle. It is chosen to be a multiple of 2*pi, hence the
      resulting rotation is then the identity matrix. This strategy warrants
      that the output matrix is a continuous function of the input vector.
      """
      theta = np.sqrt(np.sum(r**2))
      if theta > MAX_ANGLE:
          return np.eye(3)
      elif theta > SMALL_ANGLE:
          n = r / theta
          Sn = np.array(
              [[0, -n[2], n[1]], [n[2], 0, -n[0]], [-n[1], n[0], 0]]
          )
          R = (
              np.eye(3)
              + np.sin(theta) * Sn
              + (1 - np.cos(theta)) * np.dot(Sn, Sn)
          )
      else:
          Sr = np.array(
              [[0, -r[2], r[1]], [r[2], 0, -r[0]], [-r[1], r[0], 0]]
          )
          theta2 = theta * theta
          R = (
              np.eye(3)
              + (1 - theta2 / 6.0) * Sr
              + (0.5 - theta2 / 24.0) * np.dot(Sr, Sr)
          )
      return R
  
imports: import numpy
  import scipy.linalg
  from nibabel.affines import apply_affine
  from transforms3d.quaternions import mat2quat, quat2axangle
  from .transform import Transform
label_desc: Use image node_id 3 for calling a global function with example usage: rotation_vec2mat(r) and returns: R, np