diff --git a/nitransforms/io/base.py b/nitransforms/io/base.py index cbc0fc75..c5fc2646 100644 --- a/nitransforms/io/base.py +++ b/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 diff --git a/nitransforms/io/fsl.py b/nitransforms/io/fsl.py index 68fcfa35..a087a40f 100644 --- a/nitransforms/io/fsl.py +++ b/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.""" diff --git a/nitransforms/linear.py b/nitransforms/linear.py index db9cf71a..a03f2695 100644 --- a/nitransforms/linear.py +++ b/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 diff --git a/nitransforms/tests/test_linear.py b/nitransforms/tests/test_linear.py index c3c74ed4..400910a1 100644 --- a/nitransforms/tests/test_linear.py +++ b/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")