Merge pull request nipreps#534 from fliem/conform-fix-zooms

FIX: Error conforming T1w images with differing zooms before recon-all

Resolves: nipreps/smriprep#70

Author: oesteban

niworkflows/interfaces/images.py

@@ -305,31 +305,37 @@ def _run_interface(self, runtime):
             xyz_unit = "mm"
 
         # Set a 0.05mm threshold to performing rescaling
-        atol = {"meter": 1e-5, "mm": 0.01, "micron": 10}[xyz_unit]
-
-        # Rescale => change zooms
-        # Resize => update image dimensions
-        rescale = not np.allclose(zooms, target_zooms, atol=atol)
+        atol_gross = {"meter": 5e-5, "mm": 0.05, "micron": 50}[xyz_unit]
+        # if 0.01 > difference > 0.001mm, freesurfer won't be able to merge the images
+        atol_fine = {"meter": 1e-6, "mm": 0.001, "micron": 1}[xyz_unit]
+
+        # Update zooms => Modify affine
+        # Rescale => Resample to resized voxels
+        # Resize => Resample to new image dimensions
+        update_zooms = not np.allclose(zooms, target_zooms, atol=atol_fine, rtol=0)
+        rescale = not np.allclose(zooms, target_zooms, atol=atol_gross, rtol=0)
         resize = not np.all(shape == target_shape)
-        if rescale or resize:
-            if rescale:
+        resample = rescale or resize
+        if resample or update_zooms:
+            # Use an affine with the corrected zooms, whether or not we resample
+            if update_zooms:
                 scale_factor = target_zooms / zooms
-                target_affine[:3, :3] = reoriented.affine[:3, :3].dot(
-                    np.diag(scale_factor)
-                )
+                target_affine[:3, :3] = reoriented.affine[:3, :3] @ np.diag(scale_factor)
 
             if resize:
                 # The shift is applied after scaling.
                 # Use a proportional shift to maintain relative position in dataset
                 size_factor = target_span / (zooms * shape)
                 # Use integer shifts to avoid unnecessary interpolation
-                offset = (
-                    reoriented.affine[:3, 3] * size_factor - reoriented.affine[:3, 3]
-                )
+                offset = reoriented.affine[:3, 3] * size_factor - reoriented.affine[:3, 3]
                 target_affine[:3, 3] = reoriented.affine[:3, 3] + offset.astype(int)
 
-            data = nli.resample_img(reoriented, target_affine, target_shape).dataobj
-            conform_xfm = np.linalg.inv(reoriented.affine).dot(target_affine)
+            conform_xfm = np.linalg.inv(reoriented.affine) @ target_affine
+
+            # Create new image
+            data = reoriented.dataobj
+            if resample:
+                data = nli.resample_img(reoriented, target_affine, target_shape).dataobj
             reoriented = reoriented.__class__(data, target_affine, reoriented.header)
 
         # Image may be reoriented, rescaled, and/or resized

niworkflows/interfaces/tests/test_images.py

@@ -180,3 +180,35 @@ def test_IntraModalMerge(tmpdir, shape, mshape):
     merged_data = nb.load(merged).get_fdata(dtype="float32")
     new_mshape = (*mshape[:3], 2 if len(mshape) == 3 else mshape[3] * 2)
     assert merged_data.shape == new_mshape
+
+
+def test_conform_resize(tmpdir):
+    fname = str(tmpdir / "test.nii")
+
+    random_data = np.random.random(size=(5, 5, 5))
+    img = nb.Nifti1Image(random_data, np.eye(4))
+    img.to_filename(fname)
+    conform = pe.Node(im.Conform(), name="conform", base_dir=str(tmpdir))
+    conform.inputs.in_file = fname
+    conform.inputs.target_zooms = (1, 1, 1.5)
+    conform.inputs.target_shape = (5, 5, 5)
+    res = conform.run()
+
+    out_img = nb.load(res.outputs.out_file)
+    assert out_img.header.get_zooms() == conform.inputs.target_zooms
+
+
+def test_conform_set_zooms(tmpdir):
+    fname = str(tmpdir / "test.nii")
+
+    random_data = np.random.random(size=(5, 5, 5))
+    img = nb.Nifti1Image(random_data, np.eye(4))
+    img.to_filename(fname)
+    conform = pe.Node(im.Conform(), name="conform", base_dir=str(tmpdir))
+    conform.inputs.in_file = fname
+    conform.inputs.target_zooms = (1, 1, 1.002)
+    conform.inputs.target_shape = (5, 5, 5)
+    res = conform.run()
+
+    out_img = nb.load(res.outputs.out_file)
+    assert np.allclose(out_img.header.get_zooms(), conform.inputs.target_zooms)