From 8683d675832b68d626b1452a7c6744d52b2017a0 Mon Sep 17 00:00:00 2001 From: Oscar Esteban Date: Tue, 15 Feb 2022 16:10:13 +0100 Subject: [PATCH 1/3] fix(FSL): load arrays of linear transforms - one less XFail (#40) --- nitransforms/io/fsl.py | 51 +++++++++++++++++++++++++++++++++++++++++- nitransforms/linear.py | 2 ++ 2 files changed, 52 insertions(+), 1 deletion(-) diff --git a/nitransforms/io/fsl.py b/nitransforms/io/fsl.py index 68fcfa35..07f77c6c 100644 --- a/nitransforms/io/fsl.py +++ b/nitransforms/io/fsl.py @@ -1,7 +1,9 @@ """Read/write FSL's transforms.""" import os +import warnings import numpy as np from pathlib import Path +from nibabel import load as _nbload from nibabel.affines import voxel_sizes from .base import BaseLinearTransformList, LinearParameters, TransformFileError @@ -24,7 +26,19 @@ 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))) @@ -55,6 +69,35 @@ def from_string(cls, string): ) return tf + def to_ras(self, moving=None, reference=None): + """Return a nitransforms internal RAS+ matrix.""" + if moving is None: + warnings.warn( + "Converting FSL to RAS: moving image not provided, using reference." + ) + moving = reference + + if reference is None: + raise ValueError("Cannot build FSL linear transform without a reference") + + reference = _ensure_image(reference) + moving = _ensure_image(moving) + + refswp, refspc = _fsl_aff_adapt(reference) + pre = reference.affine.dot(np.linalg.inv(refspc).dot(np.linalg.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))) + + mat = self.structarr["parameters"].T + + return ( + np.linalg.inv(post) + @ np.swapaxes(np.linalg.inv(mat), 0, 1) + @ np.linalg.inv(pre) + ) + class FSLLinearTransformArray(BaseLinearTransformList): """A string-based structure for series of FSL linear transforms.""" @@ -144,3 +187,9 @@ def _fsl_aff_adapt(space): swp[0, 0] = -1.0 swp[0, 3] = (space.shape[0] - 1) * zooms[0] return swp, np.diag(zooms) + + +def _ensure_image(img): + if isinstance(img, (str, Path)): + return _nbload(img) + return img 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 From 1112ac37aaccbfdf1b84fdfbf8d1d634946e7ef1 Mon Sep 17 00:00:00 2001 From: Oscar Esteban Date: Tue, 15 Feb 2022 18:23:00 +0100 Subject: [PATCH 2/3] enh: review comments Co-authored-by: Chris Markiewicz --- nitransforms/io/base.py | 7 +++++++ nitransforms/io/fsl.py | 33 +++++++++++++-------------------- 2 files changed, 20 insertions(+), 20 deletions(-) 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 07f77c6c..6b3dffd7 100644 --- a/nitransforms/io/fsl.py +++ b/nitransforms/io/fsl.py @@ -2,11 +2,16 @@ import os import warnings import numpy as np +from numpy.linalg import inv from pathlib import Path -from nibabel import load as _nbload from nibabel.affines import voxel_sizes -from .base import BaseLinearTransformList, LinearParameters, TransformFileError +from .base import ( + BaseLinearTransformList, + LinearParameters, + TransformFileError, + _ensure_image, +) class FSLLinearTransform(LinearParameters): @@ -41,14 +46,14 @@ def from_ras(cls, ras, moving=None, reference=None): # 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 @@ -84,19 +89,13 @@ def to_ras(self, moving=None, reference=None): moving = _ensure_image(moving) refswp, refspc = _fsl_aff_adapt(reference) - pre = reference.affine.dot(np.linalg.inv(refspc).dot(np.linalg.inv(refswp))) + pre = refswp @ refspc @ inv(reference.affine) # 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 = moving.affine @ inv(movspc) @ inv(movswp) mat = self.structarr["parameters"].T - - return ( - np.linalg.inv(post) - @ np.swapaxes(np.linalg.inv(mat), 0, 1) - @ np.linalg.inv(pre) - ) + return post @ np.swapaxes(inv(mat), 0, 1) @ pre class FSLLinearTransformArray(BaseLinearTransformList): @@ -187,9 +186,3 @@ def _fsl_aff_adapt(space): swp[0, 0] = -1.0 swp[0, 3] = (space.shape[0] - 1) * zooms[0] return swp, np.diag(zooms) - - -def _ensure_image(img): - if isinstance(img, (str, Path)): - return _nbload(img) - return img From 67d56622fe77a4b6510d82138e1a3536fe39753f Mon Sep 17 00:00:00 2001 From: Oscar Esteban Date: Tue, 15 Feb 2022 22:16:39 +0100 Subject: [PATCH 3/3] fix: actually enable the tests --- nitransforms/io/fsl.py | 6 ++-- nitransforms/tests/test_linear.py | 52 ++++++++++++++++++++++++++++--- 2 files changed, 51 insertions(+), 7 deletions(-) diff --git a/nitransforms/io/fsl.py b/nitransforms/io/fsl.py index 6b3dffd7..a087a40f 100644 --- a/nitransforms/io/fsl.py +++ b/nitransforms/io/fsl.py @@ -76,15 +76,15 @@ def from_string(cls, string): 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 - if reference is None: - raise ValueError("Cannot build FSL linear transform without a reference") - reference = _ensure_image(reference) moving = _ensure_image(moving) 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")