ENH: More nifti tests and fixes

Author: cookpa

ants/utils/nifti_utils.py

@@ -266,7 +266,9 @@ def get_nifti_sform_spatial_info(metadata, shear_threshold=1e-6, max_angle_devia
 
 def get_nifti_qform_spatial_info(metadata):
     """
-    Extract qform-derived spacing, origin, direction. Uses 'qto_xyz' from the metadata dict.
+    Extract qform-derived spacing, origin, direction. Uses the 4x4 'qto_xyz' from the metadata dict. This is the rotation
+    matrix derived from quaternion parameters in the NIfTI header, multiplied by the pixdim scales and with the qoffset
+    translation.
 
     Note: output is in ITK LPS coordinates
 
@@ -358,9 +360,9 @@ def get_nifti_spatial_transform_from_metadata(metadata, prefer_sform=True, shear
             if verbose:
                 print(f"[sform] spacing={info_s['spacing']} origin={info_s['origin']} (desheared={info_s['desheared']})")
             # Verify spacing vs pixdim
-            if not _spacing_matches(info_s["spacing"], pixdims):
+            if not _spacing_matches(info_s['transform_spacing'], pixdims):
                 warnings.warn(
-                    f"sform-derived spacing {info_s['spacing']} does not match NIfTI pixdim {pixdims}; "
+                    f"sform-derived spacing {info_s['transform_spacing']} does not match NIfTI pixdim {pixdims}; "
                     "ignoring sform and trying qform.",
                     RuntimeWarning,
                 )
@@ -376,22 +378,23 @@ def get_nifti_spatial_transform_from_metadata(metadata, prefer_sform=True, shear
             warnings.warn("No usable sform or qform present in metadata; image not modified", RuntimeWarning)
             return
         info_q = get_nifti_qform_spatial_info(metadata)
-        if not _spacing_matches(info_q["spacing"], pixdims):
-                raise ValueError(f"qform-derived spacing {info_q['spacing']} does not match NIfTI pixdim {pixdims}")
+        if not _spacing_matches(info_q['transform_spacing'], pixdims):
+                raise ValueError(f"qform-derived spacing {info_q['transform_spacing']} does not match NIfTI pixdim {pixdims}")
         if verbose:
-            print(f"[qform] spacing={info_q['spacing']} origin={info_q['origin']}")
+            print(f"[qform] spacing={info_q['transform_spacing']} origin={info_q['origin']}")
+            print(f"Setting spacing to pixdims={pixdims}")
         return dict(
-            origin=info_q["origin"],
+            origin=info_q['origin'],
             spacing=pixdims,
-            direction=info_q["direction"],
+            direction=info_q['direction'],
             transform_source="qform",
         )
     else:
         # Use sform
         return dict(
-            origin=info_s["origin"],
+            origin=info_s['origin'],
             spacing=pixdims,
-            direction=info_s["direction"],
+            direction=info_s['direction'],
             transform_source="sform",
         )
 
@@ -439,11 +442,13 @@ def set_nifti_spatial_transform_from_metadata(image, metadata, prefer_sform=True
 
     if image.dimension == 4:
         # Use original definition of time spacing and origin
-        image.set_origin(info["origin"].append(image.origin[-1]))
+        origin_4d = info['origin'].copy()
+        origin_4d.append(image.origin[-1])
+        image.set_origin(origin_4d)
         # direction is 4x4 with identity time
         dir4 = np.eye(4)
-        dir4[:3, :3] = np.array(info["direction"])
+        dir4[:3, :3] = np.array(info['direction'])
         image.set_direction(dir4.tolist())
     else:
-        image.set_origin(info["origin"])
-        image.set_direction(info["direction"])
+        image.set_origin(info['origin'])
+        image.set_direction(info['direction'])

tests/test_utils.py

@@ -7,6 +7,7 @@
 from tempfile import TemporaryDirectory
 
 from common import run_tests
+from copy import deepcopy
 
 import math
 import numpy.testing as nptest
@@ -1351,19 +1352,21 @@ def get_nifti_spatial_transform_from_metadata(self):
         self.assertTrue(xform['transform_source'] == 'sform')
 
         # Flip sform just to be different
-        mni_metadata_flip = self.mni_metadata.deepcopy()
+        mni_metadata_flip = deepcopy(self.mni_metadata)
         srow_x = self.mni_metadata['srow_x'].split()
         srow_x = [str(-float(v)) for v in srow_x]
         mni_metadata_flip['srow_x'] = ' '.join(srow_x)
 
         xform = ants.get_nifti_spatial_transform_from_metadata(mni_metadata_flip)
         self.assertTrue(np.allclose(xform['direction'][0], -self.mni.direction[0], atol=1e-5))
         self.assertTrue(np.allclose(xform['origin'][0], -self.mni.origin[0], atol=1e-5))
+
+        # Prefer qform should give original
         xform_q = ants.get_nifti_spatial_transform_from_metadata(mni_metadata_flip, prefer_sform=False)
         self.assertTrue(np.allclose(xform_q['direction'], self.mni.direction, atol=1e-5))
 
         # Check that shear gets removed
-        mni_metadata_shear = self.mni_metadata.deepcopy()
+        mni_metadata_shear = deepcopy(self.mni_metadata)
         srow_x = self.mni_metadata['srow_x'].split()
         srow_x[1] = str('-0.0001')
         mni_metadata_shear['srow_x'] = ' '.join(srow_x)
@@ -1383,5 +1386,44 @@ def get_nifti_spatial_transform_from_metadata(self):
         self.assertTrue(xform['transform_source'] == 'qform')
 
 
+    def test_set_nifti_spatial_transform_from_metadata(self):
+        # Flip sform
+        mni_metadata_flip = deepcopy(self.mni_metadata)
+        srow_x = self.mni_metadata['srow_x'].split()
+        srow_x = [str(-float(v)) for v in srow_x]
+        mni_metadata_flip['srow_x'] = ' '.join(srow_x)
+
+        mni_copy = ants.image_clone(self.mni)
+        ants.set_nifti_spatial_transform_from_metadata(mni_copy, mni_metadata_flip)
+        self.assertTrue(np.allclose(mni_copy.direction[0], -self.mni.direction[0], atol=1e-5))
+        self.assertTrue(np.allclose(mni_copy.direction[1], self.mni.direction[1], atol=1e-5))
+        self.assertTrue(np.allclose(mni_copy.direction[2], self.mni.direction[2], atol=1e-5))
+        self.assertTrue(np.allclose(mni_copy.origin[0], -self.mni.origin[0], atol=1e-5))
+        self.assertTrue(np.allclose(mni_copy.origin[1], self.mni.origin[1], atol=1e-5))
+        self.assertTrue(np.allclose(mni_copy.origin[2], self.mni.origin[2], atol=1e-5))
+
+        # Make a time series image
+        mni_metadata_ts = deepcopy(mni_metadata_flip)
+        mni_metadata_ts['dim[0]'] = '4'
+        mni_metadata_ts['dim[4]'] = '3'
+        mni_metadata_ts['pixdim[4]'] = '2.0'
+        array_list = [img.numpy() for img in [self.mni, self.mni, self.mni]]
+        stacked_array = np.stack(array_list, axis=-1)
+        mni_ts = ants.from_numpy(
+            stacked_array,
+            has_components=False,
+            spacing=self.mni.spacing + (2.0,),
+            origin=self.mni.origin + (0.0,),
+            direction=np.eye(4)
+        )
+        ants.set_nifti_spatial_transform_from_metadata(mni_ts, mni_metadata_ts)
+
+        mni_ts_expected_direction = np.eye(4)
+        mni_ts_expected_direction[:3,:3] = mni_copy.direction
+
+        self.assertTrue(np.allclose(mni_ts.direction, mni_ts_expected_direction, atol=1e-5))
+        self.assertTrue(np.allclose(mni_ts.origin, mni_copy.origin + (0.0,), atol=1e-5))
+        self.assertTrue(np.allclose(mni_ts.spacing, mni_copy.spacing + (2.0,), atol=1e-5))
+
 if __name__ == "__main__":
     run_tests()