Merge branch 'sform-qform-getset' into main-master

* sform-qform-getset:
  ENH + TEST: enforce float image affines
  TEST - refactor Bago's tests to match new API
  TEST - added test for set_qform and set_sform
  NF - add get/set for sform/qform on image
  NF - add keyword arguments to sform/qform get/set

Author: matthew-brett

nibabel/nifti1.py

@@ -701,9 +701,29 @@ def get_qform_quaternion(self):
         # Adjust threshold to fact that source data was float32
         return fillpositive(bcd, FLOAT32_EPS_3)
 
-    def get_qform(self):
-        ''' Return 4x4 affine matrix from qform parameters in header '''
+    def get_qform(self, coded=False):
+        """ Return 4x4 affine matrix from qform parameters in header
+
+        Parameters
+        ----------
+        coded : bool, optional
+            If True, return {affine or None}, and qform code.  If False, just
+            return affine.  {affine or None} means, return None if qform code ==
+            0, and affine otherwise.
+
+        Returns
+        -------
+        affine : None or (4,4) ndarray
+            If `coded` is False, always return affine reconstructed from qform
+            quaternion.  If `coded` is True, return None if qform code is 0,
+            else return the affine.
+        code : int
+            Qform code. Only returned if `coded` is True.
+        """
         hdr = self._structarr
+        code = hdr['qform_code']
+        if code == 0 and coded:
+            return None, 0
         quat = self.get_qform_quaternion()
         R = quat2mat(quat)
         vox = hdr['pixdim'][1:4].copy()
@@ -718,30 +738,39 @@ def get_qform(self):
         out = np.eye(4)
         out[0:3, 0:3] = M
         out[0:3, 3] = [hdr['qoffset_x'], hdr['qoffset_y'], hdr['qoffset_z']]
+        if coded:
+            return out, code
         return out
 
-    def set_qform(self, affine, code=None):
+    def set_qform(self, affine, code=None, strip_shears=True):
         ''' Set qform header values from 4x4 affine
 
         Parameters
         ----------
         hdr : nifti1 header
-        affine : 4x4 array
-            affine transform to write into qform
+        affine : None or 4x4 array
+            affine transform to write into sform. If None, only set code.
         code : None, string or integer
             String or integer giving meaning of transform in *affine*.
-            The default is None.  If code is None, then {if current
-            qform code is not 0, leave code as it is in the header; else
-            set to 2 ('aligned')}.
+            The default is None.  If code is None, then:
+            * If affine is None, `code`-> 0
+            * If affine not None and sform code in header == 0, `code`-> 2
+              (aligned)
+            * If affine not None and sform code in header != 0, `code`-> sform
+              code in header
+        strip_shears : bool, optional
+            Whether to strip shears in `affine`.  If True, shears will be
+            silently stripped. If False, the presence of shears will raise a
+            ``HeaderDataError``
 
         Notes
         -----
         The qform transform only encodes translations, rotations and
-        zooms. If there are shear components to the `affine` transform,
-        the written qform gives the closest approximation where the
-        rotation matrix is orthogonal. This is to allow quaternion
-        representation. The orthogonal representation enforces
-        orthogonal axes.
+        zooms. If there are shear components to the `affine` transform, and
+        `strip_shears` is True (the default), the written qform gives the
+        closest approximation where the rotation matrix is orthogonal. This is
+        to allow quaternion representation. The orthogonal representation
+        enforces orthogonal axes.
 
         Examples
         --------
@@ -765,15 +794,25 @@ def set_qform(self, affine, code=None):
         >>> hdr.set_qform(affine, code='scanner')
         >>> int(hdr['qform_code'])
         1
+        >>> hdr.set_qform(None)
+        >>> int(hdr['qform_code'])
+        0
         '''
         hdr = self._structarr
+        old_code = hdr['qform_code']
         if code is None:
-            code = hdr['qform_code']
-            if code == 0: # default is 'aligned'
-                hdr['qform_code'] = 2
-        else:
+            if affine is None:
+                code = 0
+            elif old_code == 0:
+                code = 2 # aligned
+            else:
+                code = old_code
+        else: # code set
             code = self._field_recoders['qform_code'][code]
-            hdr['qform_code'] = code
+        hdr['qform_code'] = code
+        if affine is None:
+            return
+        affine = np.asarray(affine)
         if not affine.shape == (4, 4):
             raise TypeError('Need 4x4 affine as input')
         trans = affine[:3, 3]
@@ -793,6 +832,9 @@ def set_qform(self, affine, code=None):
         # orthogonal matrix PR, to input R
         P, S, Qs = npl.svd(R)
         PR = np.dot(P, Qs)
+        if not strip_shears and not np.allclose(PR, R):
+            raise HeaderDataError("Shears in affine and `strip_shears` is "
+                                  "False")
         # Convert to quaternion
         quat = mat2quat(PR)
         # Set into header
@@ -801,13 +843,35 @@ def set_qform(self, affine, code=None):
         hdr['pixdim'][1:4] = zooms
         hdr['quatern_b'], hdr['quatern_c'], hdr['quatern_d'] = quat[1:]
 
-    def get_sform(self):
-        ''' Return sform 4x4 affine matrix from header '''
+    def get_sform(self, coded=False):
+        """ Return 4x4 affine matrix from sform parameters in header
+
+        Parameters
+        ----------
+        coded : bool, optional
+            If True, return {affine or None}, and sform code.  If False, just
+            return affine.  {affine or None} means, return None if sform code ==
+            0, and affine otherwise.
+
+        Returns
+        -------
+        affine : None or (4,4) ndarray
+            If `coded` is False, always return affine from sform fields. If
+            `coded` is True, return None if sform code is 0, else return the
+            affine.
+        code : int
+            Sform code. Only returned if `coded` is True.
+        """
         hdr = self._structarr
+        code = hdr['sform_code']
+        if code == 0 and coded:
+            return None, 0
         out = np.eye(4)
         out[0, :] = hdr['srow_x'][:]
         out[1, :] = hdr['srow_y'][:]
         out[2, :] = hdr['srow_z'][:]
+        if coded:
+            return out, code
         return out
 
     def set_sform(self, affine, code=None):
@@ -816,13 +880,16 @@ def set_sform(self, affine, code=None):
         Parameters
         ----------
         hdr : nifti1 header
-        affine : 4x4 array
-            affine transform to write into sform
+        affine : None or 4x4 array
+            affine transform to write into sform.  If None, only set `code`
         code : None, string or integer
             String or integer giving meaning of transform in *affine*.
-            The default is None.  If code is None, then {if current
-            sform code is not 0, leave code as it is in the header; else
-            set to 2 ('aligned')}.
+            The default is None.  If code is None, then:
+            * If affine is None, `code`-> 0
+            * If affine not None and sform code in header == 0, `code`-> 2
+              (aligned)
+            * If affine not None and sform code in header != 0, `code`-> sform
+              code in header
 
         Examples
         --------
@@ -846,15 +913,25 @@ def set_sform(self, affine, code=None):
         >>> hdr.set_sform(affine, code='scanner')
         >>> int(hdr['sform_code'])
         1
+        >>> hdr.set_sform(None)
+        >>> int(hdr['sform_code'])
+        0
         '''
         hdr = self._structarr
+        old_code = hdr['sform_code']
         if code is None:
-            code = hdr['sform_code']
-            if code == 0:
-                hdr['sform_code'] = 2 # aligned
-        else:
+            if affine is None:
+                code = 0
+            elif old_code == 0:
+                code = 2 # aligned
+            else:
+                code = old_code
+        else: # code set
             code = self._field_recoders['sform_code'][code]
-            hdr['sform_code'] = code
+        hdr['sform_code'] = code
+        if affine is None:
+            return
+        affine = np.asarray(affine)
         hdr['srow_x'][:] = affine[0, :]
         hdr['srow_y'][:] = affine[1, :]
         hdr['srow_z'][:] = affine[2, :]
@@ -1512,6 +1589,170 @@ def update_header(self):
         # Make qform 'unknown'
         hdr.set_qform(self._affine, code='unknown')
 
+    def get_qform(self, coded=False):
+        """ Return 4x4 affine matrix from qform parameters in header
+
+        Parameters
+        ----------
+        coded : bool, optional
+            If True, return {affine or None}, and qform code.  If False, just
+            return affine.  {affine or None} means, return None if qform code ==
+            0, and affine otherwise.
+
+        Returns
+        -------
+        affine : None or (4,4) ndarray
+            If `coded` is False, always return affine reconstructed from qform
+            quaternion.  If `coded` is True, return None if qform code is 0,
+            else return the affine.
+        code : int
+            Qform code. Only returned if `coded` is True.
+
+        See also
+        --------
+        Nifti1Header.set_qform
+        """
+        return self._header.get_qform(coded)
+
+    def set_qform(self, affine, code=None, strip_shears=True, **kwargs):
+        ''' Set qform header values from 4x4 affine
+
+        Parameters
+        ----------
+        hdr : nifti1 header
+        affine : None or 4x4 array
+            affine transform to write into sform. If None, only set code.
+        code : None, string or integer
+            String or integer giving meaning of transform in *affine*.
+            The default is None.  If code is None, then:
+            * If affine is None, `code`-> 0
+            * If affine not None and sform code in header == 0, `code`-> 2
+              (aligned)
+            * If affine not None and sform code in header != 0, `code`-> sform
+              code in header
+        strip_shears : bool, optional
+            Whether to strip shears in `affine`.  If True, shears will be
+            silently stripped. If False, the presence of shears will raise a
+            ``HeaderDataError``
+        update_affine : bool, optional
+            Whether to update the image affine from the header best affine after
+            setting the qform. Must be keyword argumemt (because of different
+            position in `set_qform`). Default is True
+
+        See also
+        --------
+        Nifti1Header.set_qform
+
+        Examples
+        --------
+        >>> data = np.arange(24).reshape((2,3,4))
+        >>> aff = np.diag([2, 3, 4, 1])
+        >>> img = Nifti1Pair(data, aff)
+        >>> img.get_qform()
+        array([[ 2.,  0.,  0.,  0.],
+               [ 0.,  3.,  0.,  0.],
+               [ 0.,  0.,  4.,  0.],
+               [ 0.,  0.,  0.,  1.]])
+        >>> img.get_qform(coded=True)
+        (None, 0)
+        >>> aff2 = np.diag([3, 4, 5, 1])
+        >>> img.set_qform(aff2, 'talairach')
+        >>> qaff, code = img.get_qform(coded=True)
+        >>> np.all(qaff == aff2)
+        True
+        >>> int(code)
+        3
+        '''
+        update_affine = kwargs.pop('update_affine', True)
+        if kwargs:
+            raise TypeError('Unexpected keyword argument(s) %s' % kwargs)
+        self._header.set_qform(affine, code, strip_shears)
+        if update_affine:
+            self._affine[:] = self._header.get_best_affine()
+
+    def get_sform(self, coded=False):
+        """ Return 4x4 affine matrix from sform parameters in header
+
+        Parameters
+        ----------
+        coded : bool, optional
+            If True, return {affine or None}, and sform code.  If False, just
+            return affine.  {affine or None} means, return None if sform code ==
+            0, and affine otherwise.
+
+        Returns
+        -------
+        affine : None or (4,4) ndarray
+            If `coded` is False, always return affine from sform fields. If
+            `coded` is True, return None if sform code is 0, else return the
+            affine.
+        code : int
+            Sform code. Only returned if `coded` is True.
+
+        See also
+        --------
+        Nifti1Header.get_sform
+        """
+        return self._header.get_sform(coded)
+
+    def set_sform(self, affine, code=None, **kwargs):
+        ''' Set sform transform from 4x4 affine
+
+        Parameters
+        ----------
+        hdr : nifti1 header
+        affine : None or 4x4 array
+            affine transform to write into sform.  If None, only set `code`
+        code : None, string or integer
+            String or integer giving meaning of transform in *affine*.
+            The default is None.  If code is None, then:
+            * If affine is None, `code`-> 0
+            * If affine not None and sform code in header == 0, `code`-> 2
+              (aligned)
+            * If affine not None and sform code in header != 0, `code`-> sform
+              code in header
+        update_affine : bool, optional
+            Whether to update the image affine from the header best affine after
+            setting the qform.  Must be keyword argumemt (because of different
+            position in `set_qform`). Default is True
+
+        See also
+        --------
+        Nifti1Header.set_sform
+
+        Examples
+        --------
+        >>> data = np.arange(24).reshape((2,3,4))
+        >>> aff = np.diag([2, 3, 4, 1])
+        >>> img = Nifti1Pair(data, aff)
+        >>> img.get_sform()
+        array([[ 2.,  0.,  0.,  0.],
+               [ 0.,  3.,  0.,  0.],
+               [ 0.,  0.,  4.,  0.],
+               [ 0.,  0.,  0.,  1.]])
+        >>> saff, code = img.get_sform(coded=True)
+        >>> saff
+        array([[ 2.,  0.,  0.,  0.],
+               [ 0.,  3.,  0.,  0.],
+               [ 0.,  0.,  4.,  0.],
+               [ 0.,  0.,  0.,  1.]])
+        >>> int(code)
+        2
+        >>> aff2 = np.diag([3, 4, 5, 1])
+        >>> img.set_sform(aff2, 'talairach')
+        >>> saff, code = img.get_sform(coded=True)
+        >>> np.all(saff == aff2)
+        True
+        >>> int(code)
+        3
+        '''
+        update_affine = kwargs.pop('update_affine', True)
+        if kwargs:
+            raise TypeError('Unexpected keyword argument(s) %s' % kwargs)
+        self._header.set_sform(affine, code)
+        if update_affine:
+            self._affine[:] = self._header.get_best_affine()
+
 
 class Nifti1Image(Nifti1Pair):
     header_class = Nifti1Header

nibabel/spatialimages.py

@@ -301,8 +301,9 @@ def __init__(self, data, affine, header=None,
             # Check that affine is array-like 4,4.  Maybe this is too strict at
             # this abstract level, but so far I think all image formats we know
             # do need 4,4.
-            # Copy affine to  isolate from environment
-            affine = np.array(affine, copy=True)
+            # Copy affine to  isolate from environment.  Specify float type to
+            # avoid surprising integer rounding when setting values into affine
+            affine = np.array(affine, dtype=np.float64, copy=True)
             if not affine.shape == (4,4):
                 raise ValueError('Affine should be shape 4,4')
         self._affine = affine