FIX: Load arrays of linear transforms from FSL files

nitransforms/io/base.py

@@ -1,4 +1,5 @@
 """Read/write linear transforms."""
+from pathlib import Path
 import numpy as np
 from nibabel import load as loadimg
 from scipy.io.matlab.miobase import get_matfile_version
@@ -174,3 +175,9 @@ def _read_mat(byte_stream):
     else:
         raise TransformFileError("Not a Matlab file.")
     return reader.get_variables()
+
+
+def _ensure_image(img):
+    if isinstance(img, (str, Path)):
+        return loadimg(img)
+    return img

nitransforms/io/fsl.py

@@ -1,10 +1,17 @@
 """Read/write FSL's transforms."""
 import os
+import warnings
 import numpy as np
+from numpy.linalg import inv
 from pathlib import Path
 from nibabel.affines import voxel_sizes
 
-from .base import BaseLinearTransformList, LinearParameters, TransformFileError
+from .base import (
+    BaseLinearTransformList,
+    LinearParameters,
+    TransformFileError,
+    _ensure_image,
+)
 
 
 class FSLLinearTransform(LinearParameters):
@@ -24,17 +31,29 @@ def to_string(self):
 
     @classmethod
     def from_ras(cls, ras, moving=None, reference=None):
-        """Create an ITK affine from a nitransform's RAS+ matrix."""
+        """Create an FSL affine from a nitransform's RAS+ matrix."""
+        if moving is None:
+            warnings.warn(
+                "[Converting FSL to RAS] moving not provided, using reference as moving"
+            )
+            moving = reference
+
+        if reference is None:
+            raise ValueError("Cannot build FSL linear transform without a reference")
+
+        reference = _ensure_image(reference)
+        moving = _ensure_image(moving)
+
         # Adjust for reference image offset and orientation
         refswp, refspc = _fsl_aff_adapt(reference)
-        pre = reference.affine.dot(np.linalg.inv(refspc).dot(np.linalg.inv(refswp)))
+        pre = reference.affine.dot(inv(refspc).dot(inv(refswp)))
 
         # Adjust for moving image offset and orientation
         movswp, movspc = _fsl_aff_adapt(moving)
-        post = np.linalg.inv(movswp).dot(movspc.dot(np.linalg.inv(moving.affine)))
+        post = inv(movswp).dot(movspc.dot(inv(moving.affine)))
 
         # Compose FSL transform
-        mat = np.linalg.inv(np.swapaxes(post.dot(ras.dot(pre)), 0, 1))
+        mat = inv(np.swapaxes(post.dot(ras.dot(pre)), 0, 1))
 
         tf = cls()
         tf.structarr["parameters"] = mat.T
@@ -55,6 +74,29 @@ def from_string(cls, string):
         )
         return tf
 
+    def to_ras(self, moving=None, reference=None):
+        """Return a nitransforms internal RAS+ matrix."""
+        if reference is None:
+            raise ValueError("Cannot build FSL linear transform without a reference")
+
+        if moving is None:
+            warnings.warn(
+                "Converting FSL to RAS: moving image not provided, using reference."
+            )
+            moving = reference
+
+        reference = _ensure_image(reference)
+        moving = _ensure_image(moving)
+
+        refswp, refspc = _fsl_aff_adapt(reference)
+
+        pre = refswp @ refspc @ inv(reference.affine)
+        # Adjust for moving image offset and orientation
+        movswp, movspc = _fsl_aff_adapt(moving)
+        post = moving.affine @ inv(movspc) @ inv(movswp)
+        mat = self.structarr["parameters"].T
+        return post @ np.swapaxes(inv(mat), 0, 1) @ pre
+
 
 class FSLLinearTransformArray(BaseLinearTransformList):
     """A string-based structure for series of FSL linear transforms."""

nitransforms/linear.py

@@ -241,6 +241,8 @@ def from_filename(cls, filename, fmt="X5", reference=None, moving=None):
             _factory = io.itk.ITKLinearTransformArray
         elif fmt.lower() in ("lta", "fs"):
             _factory = io.LinearTransformArray
+        elif fmt.lower() == "fsl":
+            _factory = io.fsl.FSLLinearTransformArray
         else:
             raise NotImplementedError
 

nitransforms/tests/test_linear.py

@@ -82,18 +82,62 @@ def test_loadsave(tmp_path, data_path, testdata_path, fmt):
 
     fname = tmp_path / ".".join(("transform-mapping", fmt))
     xfm.to_filename(fname, fmt=fmt)
-    xfm == nitl.load(fname, fmt=fmt, reference=ref_file)
+
+    if fmt == "fsl":
+        # FSL should not read a transform without reference
+        with pytest.raises(ValueError):
+            nitl.load(fname, fmt=fmt)
+            nitl.load(fname, fmt=fmt, moving=ref_file)
+
+        with pytest.warns(UserWarning):
+            assert np.allclose(
+                xfm.matrix,
+                nitl.load(fname, fmt=fmt, reference=ref_file).matrix,
+                rtol=1e-2,  # FSL incurs into large errors due to rounding
+            )
+
+        assert np.allclose(
+            xfm.matrix,
+            nitl.load(fname, fmt=fmt, reference=ref_file, moving=ref_file).matrix,
+            rtol=1e-2,  # FSL incurs into large errors due to rounding
+        )
+    else:
+        assert xfm == nitl.load(fname, fmt=fmt, reference=ref_file)
+
     xfm.to_filename(fname, fmt=fmt, moving=ref_file)
-    xfm == nitl.load(fname, fmt=fmt, reference=ref_file)
+
+    if fmt == "fsl":
+        assert np.allclose(
+            xfm.matrix,
+            nitl.load(fname, fmt=fmt, reference=ref_file, moving=ref_file).matrix,
+            rtol=1e-2,  # FSL incurs into large errors due to rounding
+        )
+    else:
+        assert xfm == nitl.load(fname, fmt=fmt, reference=ref_file)
 
     ref_file = testdata_path / "someones_anatomy.nii.gz"
     xfm = nitl.load(data_path / "affine-LAS.itk.tfm", fmt="itk")
     xfm.reference = ref_file
     fname = tmp_path / ".".join(("single-transform", fmt))
     xfm.to_filename(fname, fmt=fmt)
-    xfm == nitl.load(fname, fmt=fmt, reference=ref_file)
+    if fmt == "fsl":
+        assert np.allclose(
+            xfm.matrix,
+            nitl.load(fname, fmt=fmt, reference=ref_file, moving=ref_file).matrix,
+            rtol=1e-2,  # FSL incurs into large errors due to rounding
+        )
+    else:
+        assert xfm == nitl.load(fname, fmt=fmt, reference=ref_file)
+
     xfm.to_filename(fname, fmt=fmt, moving=ref_file)
-    xfm == nitl.load(fname, fmt=fmt, reference=ref_file)
+    if fmt == "fsl":
+        assert np.allclose(
+            xfm.matrix,
+            nitl.load(fname, fmt=fmt, reference=ref_file, moving=ref_file).matrix,
+            rtol=1e-2,  # FSL incurs into large errors due to rounding
+        )
+    else:
+        assert xfm == nitl.load(fname, fmt=fmt, reference=ref_file)
 
 
 @pytest.mark.xfail(reason="Not fully implemented")