From 69c02be409145ae5ff7befc0508b208df42e4f2e Mon Sep 17 00:00:00 2001 From: mathiasg Date: Mon, 3 Jun 2024 15:56:10 -0400 Subject: [PATCH 1/6] FIX: Ensure shape matches template The reverse of the TransformFixedParameters should match the template dimensions. This fell through the cracks because the templates tested (MNI152 variants) all have the same X and Z dimensions. Only when testing on a template with distinct values for shape did this become evident. Additionally, this removes any `FIXED_PARAMS` checks, as that is template specific and more useful for debugging the initial implementation. --- fmriprep/utils/transforms.py | 20 +------------------- 1 file changed, 1 insertion(+), 19 deletions(-) diff --git a/fmriprep/utils/transforms.py b/fmriprep/utils/transforms.py index 7b36eb359..d7f0210ad 100644 --- a/fmriprep/utils/transforms.py +++ b/fmriprep/utils/transforms.py @@ -1,6 +1,5 @@ """Utilities for loading transforms for resampling""" -import warnings from pathlib import Path import h5py @@ -38,16 +37,6 @@ def load_transforms(xfm_paths: list[Path], inverse: list[bool]) -> nt.base.Trans return chain -FIXED_PARAMS = np.array([ - 193.0, 229.0, 193.0, # Size - 96.0, 132.0, -78.0, # Origin - 1.0, 1.0, 1.0, # Spacing - -1.0, 0.0, 0.0, # Directions - 0.0, -1.0, 0.0, - 0.0, 0.0, 1.0, -]) # fmt:skip - - def load_ants_h5(filename: Path) -> nt.base.TransformBase: """Load ANTs H5 files as a nitransforms TransformChain""" # Borrowed from https://github.com/feilong/process @@ -81,15 +70,8 @@ def load_ants_h5(filename: Path) -> nt.base.TransformBase: raise ValueError(msg) fixed_params = transform2['TransformFixedParameters'][:] - if not np.array_equal(fixed_params, FIXED_PARAMS): - msg = 'Unexpected fixed parameters\n' - msg += f'Expected: {FIXED_PARAMS}\n' - msg += f'Found: {fixed_params}' - if not np.array_equal(fixed_params[6:], FIXED_PARAMS[6:]): - raise ValueError(msg) - warnings.warn(msg, stacklevel=1) - shape = tuple(fixed_params[:3].astype(int)) + shape = tuple(fixed_params[:3].astype(int)[::-1]) warp = h['TransformGroup']['2']['TransformParameters'][:] warp = warp.reshape((*shape, 3)).transpose(2, 1, 0, 3) warp *= np.array([-1, -1, 1]) From 83d4461db06c4f8e59e8078da36744cdac74e32a Mon Sep 17 00:00:00 2001 From: Mathias Goncalves Date: Tue, 4 Jun 2024 12:00:01 -0400 Subject: [PATCH 2/6] Update fmriprep/utils/transforms.py Co-authored-by: Chris Markiewicz --- fmriprep/utils/transforms.py | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/fmriprep/utils/transforms.py b/fmriprep/utils/transforms.py index d7f0210ad..cf53a1f36 100644 --- a/fmriprep/utils/transforms.py +++ b/fmriprep/utils/transforms.py @@ -71,9 +71,10 @@ def load_ants_h5(filename: Path) -> nt.base.TransformBase: fixed_params = transform2['TransformFixedParameters'][:] - shape = tuple(fixed_params[:3].astype(int)[::-1]) - warp = h['TransformGroup']['2']['TransformParameters'][:] - warp = warp.reshape((*shape, 3)).transpose(2, 1, 0, 3) + shape = tuple(fixed_params[:3].astype(int)) + # ITK stores warps in Fortran-order, where the vector components change fastest + # Nitransforms expects 3 volumes, not a volume of three-vectors, so transpose + warp = np.reshape(transform2['TransformParameters'], (3, *shape), order='F').transpose(1, 2, 3, 0) warp *= np.array([-1, -1, 1]) warp_affine = np.eye(4) From 0293a18008cee2efd61d3b2d975932b274038cfb Mon Sep 17 00:00:00 2001 From: Chris Markiewicz Date: Tue, 4 Jun 2024 12:26:47 -0400 Subject: [PATCH 3/6] STY: ruff, use trailing comma to force reshape to expand --- fmriprep/utils/transforms.py | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/fmriprep/utils/transforms.py b/fmriprep/utils/transforms.py index cf53a1f36..23105efa7 100644 --- a/fmriprep/utils/transforms.py +++ b/fmriprep/utils/transforms.py @@ -74,7 +74,11 @@ def load_ants_h5(filename: Path) -> nt.base.TransformBase: shape = tuple(fixed_params[:3].astype(int)) # ITK stores warps in Fortran-order, where the vector components change fastest # Nitransforms expects 3 volumes, not a volume of three-vectors, so transpose - warp = np.reshape(transform2['TransformParameters'], (3, *shape), order='F').transpose(1, 2, 3, 0) + warp = np.reshape( + transform2['TransformParameters'], + (3, *shape), + order='F', + ).transpose(1, 2, 3, 0) warp *= np.array([-1, -1, 1]) warp_affine = np.eye(4) From a3a952f8b8573f7fdb948949e45b2bc1cbc2c398 Mon Sep 17 00:00:00 2001 From: Chris Markiewicz Date: Wed, 5 Jun 2024 10:56:41 -0400 Subject: [PATCH 4/6] RF: Use all fixed parameters to construct warp affine Break the transformation steps down to be very explicit, with comments. --- fmriprep/utils/transforms.py | 45 ++++++++++++++++++++++++------------ 1 file changed, 30 insertions(+), 15 deletions(-) diff --git a/fmriprep/utils/transforms.py b/fmriprep/utils/transforms.py index 23105efa7..ad5458673 100644 --- a/fmriprep/utils/transforms.py +++ b/fmriprep/utils/transforms.py @@ -7,6 +7,7 @@ import nitransforms as nt import numpy as np from nitransforms.io.itk import ITKCompositeH5 +from tranforms3d.affines import compose as compose_affine def load_transforms(xfm_paths: list[Path], inverse: list[bool]) -> nt.base.TransformBase: @@ -45,7 +46,8 @@ def load_ants_h5(filename: Path) -> nt.base.TransformBase: # Changes: # * Tolerate a missing displacement field # * Return the original affine without a round-trip - # * Always return a nitransforms TransformChain + # * Always return a nitransforms TransformBase + # * Construct warp affine from fixed parameters # # This should be upstreamed into nitransforms h = h5py.File(filename) @@ -69,22 +71,35 @@ def load_ants_h5(filename: Path) -> nt.base.TransformBase: msg += f'[{i}]: {h["TransformGroup"][i]["TransformType"][:][0]}\n' raise ValueError(msg) - fixed_params = transform2['TransformFixedParameters'][:] + # Warp field fixed parameters as defined in + # https://itk.org/Doxygen/html/classitk_1_1DisplacementFieldTransform.html + shape = transform2['TransformFixedParameters'][:3] + origin = transform2['TransformFixedParameters'][3:6] + spacing = transform2['TransformFixedParameters'][6:9] + direction = transform2['TransformFixedParameters'][9:].reshape((3, 3)) + + # We are not yet confident that we handle non-unit spacing + # or direction cosine ordering correctly. + # If we confirm or fix, we can remove these checks. + if not np.allclose(spacing, 1): + raise ValueError(f'Unexpected spacing: {spacing}') + if not np.allclose(direction, direction.T): + raise ValueError(f'Asymmetric direction matrix: {direction}') + + # ITK uses LPS affines + lps_affine = compose_affine(T=origin, R=direction, Z=spacing) + ras_affine = np.diag([-1, -1, 1, 1]) @ lps_affine - shape = tuple(fixed_params[:3].astype(int)) # ITK stores warps in Fortran-order, where the vector components change fastest - # Nitransforms expects 3 volumes, not a volume of three-vectors, so transpose - warp = np.reshape( + # Vectors are in mm LPS + itk_warp = np.reshape( transform2['TransformParameters'], - (3, *shape), + (3, *shape.astype(int)), order='F', - ).transpose(1, 2, 3, 0) - warp *= np.array([-1, -1, 1]) - - warp_affine = np.eye(4) - warp_affine[:3, :3] = fixed_params[9:].reshape((3, 3)) - warp_affine[:3, 3] = fixed_params[3:6] - lps_to_ras = np.eye(4) * np.array([-1, -1, 1, 1]) - warp_affine = lps_to_ras @ warp_affine - transforms.insert(0, nt.DenseFieldTransform(nb.Nifti1Image(warp, warp_affine))) + ) + + # Nitransforms warps are in RAS, with the vector components changing slowest + nt_warp = itk_warp.transpose(1, 2, 3, 0) * np.array([-1, -1, 1]) + + transforms.insert(0, nt.DenseFieldTransform(nb.Nifti1Image(nt_warp, ras_affine))) return nt.TransformChain(transforms) From 55bea47152df32c7ff91fd45a20f55b3fea18ba7 Mon Sep 17 00:00:00 2001 From: Chris Markiewicz Date: Wed, 5 Jun 2024 11:47:41 -0400 Subject: [PATCH 5/6] MNT: Directly depend on transforms3d --- pyproject.toml | 1 + 1 file changed, 1 insertion(+) diff --git a/pyproject.toml b/pyproject.toml index c65d9d416..f707950ac 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -36,6 +36,7 @@ dependencies = [ "smriprep @ git+https://github.com/nipreps/smriprep.git@master", "tedana >= 23.0.2", "templateflow >= 24.1.0", + "transforms3d", "toml", "codecarbon", "APScheduler", From 758e5ac3c44bbcc1fbc68fd975e0305c08cee7da Mon Sep 17 00:00:00 2001 From: Chris Markiewicz Date: Wed, 5 Jun 2024 12:09:24 -0400 Subject: [PATCH 6/6] Update fmriprep/utils/transforms.py --- fmriprep/utils/transforms.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/fmriprep/utils/transforms.py b/fmriprep/utils/transforms.py index ad5458673..4c9850fa3 100644 --- a/fmriprep/utils/transforms.py +++ b/fmriprep/utils/transforms.py @@ -7,7 +7,7 @@ import nitransforms as nt import numpy as np from nitransforms.io.itk import ITKCompositeH5 -from tranforms3d.affines import compose as compose_affine +from transforms3d.affines import compose as compose_affine def load_transforms(xfm_paths: list[Path], inverse: list[bool]) -> nt.base.TransformBase: