BF: keep the sform and qform codes across a reorientation

nibabel/nifti1.py

@@ -29,6 +29,7 @@
 from .deprecated import alert_future_error
 from .filebasedimages import ImageFileError, SerializableImage
 from .optpkg import optional_package
+from .orientations import inv_ornt_aff
 from .quaternions import fillpositive, mat2quat, quat2mat
 from .spatialimages import HeaderDataError
 from .spm99analyze import SpmAnalyzeHeader
@@ -2402,6 +2403,26 @@ def as_reoriented(self, ornt):
 
         img.header.set_dim_info(*new_dim)
 
+        # Reorienting permutes and flips axes; it does not move the image into
+        # a different space.  Writing the new affine into the header goes
+        # through _affine2header, which files it under the default 'aligned'
+        # sform code and drops the qform, so a scanner anat image came back
+        # labelled as merely aligned (gh-1427).  Restore whichever codes were
+        # set before.  Each form is reoriented from its own matrix, not from
+        # img.affine, which is the sform whenever one is set: a qform and an
+        # sform may describe different spaces (scanner and MNI, say), and the
+        # qform slot must not end up holding the sform's matrix under the
+        # qform's code.  A code of 0 is left alone: the image had no form of
+        # that kind to preserve, and clearing what _affine2header just set
+        # could leave the image with no valid form at all.
+        vox_ornt = inv_ornt_aff(ornt, self.shape)
+        sform, sform_code = self.header.get_sform(coded=True)
+        qform, qform_code = self.header.get_qform(coded=True)
+        if sform_code != 0:
+            img.header.set_sform(sform @ vox_ornt, code=int(sform_code))
+        if qform_code != 0:
+            img.header.set_qform(qform @ vox_ornt, code=int(qform_code))
+
         return img
 
 

nibabel/tests/test_nifti1.py

@@ -41,7 +41,7 @@
 from nibabel.tmpdirs import InTemporaryDirectory
 
 from ..freesurfer import load as mghload
-from ..orientations import aff2axcodes
+from ..orientations import aff2axcodes, inv_ornt_aff
 from ..testing import (
     bytesio_filemap,
     bytesio_round_trip,
@@ -1573,6 +1573,62 @@ def test_rt_bias(self):
                 bias_thresh = np.max([max_miss / np.sqrt(count), eps])
                 assert np.abs(bias) < bias_thresh
 
+    def test_reoriented_preserves_xform_codes(self):
+        # gh-1427: reorienting permutes and flips axes, it does not move the
+        # image into another space, so the sform/qform codes must survive.
+        # dcm2niix writes scanner anat (code 1) for both.
+        affine = np.diag([-2.0, 2, 2, 1])
+        affine[:3, 3] = [10, -20, -30]
+        img = self.single_class(np.zeros((4, 4, 4)), affine)
+        img.header.set_sform(affine, code=1)
+        img.header.set_qform(affine, code=1)
+
+        reoriented = img.as_reoriented([[0, -1], [1, 1], [2, 1]])
+
+        assert int(reoriented.header['sform_code']) == 1
+        assert int(reoriented.header['qform_code']) == 1
+        # The header must still describe the affine the image reports.
+        assert_array_equal(reoriented.affine, reoriented.header.get_best_affine())
+
+    def test_reoriented_keeps_qform_and_sform_in_their_own_spaces(self):
+        # gh-1546 review: a scanner qform (code 1) next to an sform registered
+        # to MNI (code 4).  Each form must be reoriented from its own matrix;
+        # img.affine is the sform here, and copying it into the qform slot
+        # would label MNI coordinates as scanner coordinates.
+        shape = (4, 5, 6)
+        qform = np.diag([-2.0, 2, 2, 1])
+        qform[:3, 3] = [10, -20, -30]
+        sform = np.array([[0, 0, -2.0, 5], [0, 2.0, 0, -7], [2.0, 0, 0, 3], [0, 0, 0, 1]])
+        img = self.single_class(np.zeros(shape), sform)
+        img.header.set_sform(sform, code=4)
+        img.header.set_qform(qform, code=1)
+        ornt = [[1, -1], [0, 1], [2, -1]]
+
+        reoriented = img.as_reoriented(ornt)
+
+        vox_ornt = inv_ornt_aff(ornt, shape)
+        new_sform, new_scode = reoriented.header.get_sform(coded=True)
+        new_qform, new_qcode = reoriented.header.get_qform(coded=True)
+        assert new_scode == 4
+        assert new_qcode == 1
+        assert_almost_equal(new_sform, sform @ vox_ornt)
+        assert_almost_equal(new_qform, qform @ vox_ornt)
+        assert_array_equal(reoriented.affine, reoriented.header.get_best_affine())
+
+    def test_reoriented_leaves_unset_xform_codes_alone(self):
+        # An image with no qform has nothing to preserve there, and clearing
+        # what the affine write-back set could leave it with no valid form.
+        affine = np.diag([-2.0, 2, 2, 1])
+        img = self.single_class(np.zeros((4, 4, 4)), affine)
+        img.header.set_sform(affine, code=4)
+        img.header.set_qform(None, code=0)
+
+        reoriented = img.as_reoriented([[0, -1], [1, 1], [2, 1]])
+
+        assert int(reoriented.header['sform_code']) == 4
+        assert int(reoriented.header['qform_code']) == 0
+        assert_array_equal(reoriented.affine, reoriented.header.get_best_affine())
+
     def test_reoriented_dim_info(self):
         # Check that dim_info is reoriented correctly
         arr = np.arange(24, dtype='f4').reshape((2, 3, 4))