From 18a28efe090caabd0aad947e22f2954ba39d3472 Mon Sep 17 00:00:00 2001 From: Chris Markiewicz Date: Wed, 9 Oct 2024 20:44:52 -0400 Subject: [PATCH 1/5] ref: Drop custom load_ants_h5 for fixed nt.manip.load --- fmriprep/utils/transforms.py | 75 +----------------------------------- pyproject.toml | 2 +- 2 files changed, 3 insertions(+), 74 deletions(-) diff --git a/fmriprep/utils/transforms.py b/fmriprep/utils/transforms.py index 4c9850fa3..e15010954 100644 --- a/fmriprep/utils/transforms.py +++ b/fmriprep/utils/transforms.py @@ -2,12 +2,7 @@ from pathlib import Path -import h5py -import nibabel as nb import nitransforms as nt -import numpy as np -from nitransforms.io.itk import ITKCompositeH5 -from transforms3d.affines import compose as compose_affine def load_transforms(xfm_paths: list[Path], inverse: list[bool]) -> nt.base.TransformBase: @@ -24,7 +19,8 @@ def load_transforms(xfm_paths: list[Path], inverse: list[bool]) -> nt.base.Trans for path, inv in zip(xfm_paths[::-1], inverse[::-1], strict=False): path = Path(path) if path.suffix == '.h5': - xfm = load_ants_h5(path) + # Load as a TransformChain + xfm = nt.manip.load(path) else: xfm = nt.linear.load(path) if inv: @@ -36,70 +32,3 @@ def load_transforms(xfm_paths: list[Path], inverse: list[bool]) -> nt.base.Trans if chain is None: chain = nt.base.TransformBase() return chain - - -def load_ants_h5(filename: Path) -> nt.base.TransformBase: - """Load ANTs H5 files as a nitransforms TransformChain""" - # Borrowed from https://github.com/feilong/process - # process.resample.parse_combined_hdf5() - # - # Changes: - # * Tolerate a missing displacement field - # * Return the original affine without a round-trip - # * Always return a nitransforms TransformBase - # * Construct warp affine from fixed parameters - # - # This should be upstreamed into nitransforms - h = h5py.File(filename) - xform = ITKCompositeH5.from_h5obj(h) - - # nt.Affine - transforms = [nt.Affine(xform[0].to_ras())] - - if '2' not in h['TransformGroup']: - return transforms[0] - - transform2 = h['TransformGroup']['2'] - - # Confirm these transformations are applicable - if transform2['TransformType'][:][0] not in ( - b'DisplacementFieldTransform_float_3_3', - b'DisplacementFieldTransform_double_3_3', - ): - msg = 'Unknown transform type [2]\n' - for i in h['TransformGroup'].keys(): - msg += f'[{i}]: {h["TransformGroup"][i]["TransformType"][:][0]}\n' - raise ValueError(msg) - - # 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 - - # ITK stores warps in Fortran-order, where the vector components change fastest - # Vectors are in mm LPS - itk_warp = np.reshape( - transform2['TransformParameters'], - (3, *shape.astype(int)), - order='F', - ) - - # 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) diff --git a/pyproject.toml b/pyproject.toml index 59434ca0e..6af2e831d 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -25,7 +25,7 @@ dependencies = [ "nipype >= 1.8.5", "nireports >= 23.2.2", "nitime", - "nitransforms >= 21.0.0, != 24.0.0", + "nitransforms >= 24.0.2", "niworkflows >= 1.11.0", "numpy >= 1.22", "packaging", From 4db1127d3f32b8615b30a451027ff9adc64a30ed Mon Sep 17 00:00:00 2001 From: Chris Markiewicz Date: Wed, 9 Oct 2024 20:53:31 -0400 Subject: [PATCH 2/5] ref: Switch to nitransforms.resampling.apply --- fmriprep/interfaces/resampling.py | 3 ++- fmriprep/workflows/bold/registration.py | 7 +++++-- 2 files changed, 7 insertions(+), 3 deletions(-) diff --git a/fmriprep/interfaces/resampling.py b/fmriprep/interfaces/resampling.py index 788739bfc..de5c7fcf7 100644 --- a/fmriprep/interfaces/resampling.py +++ b/fmriprep/interfaces/resampling.py @@ -6,6 +6,7 @@ import nibabel as nb import nitransforms as nt +import nitransforms.resampling import numpy as np from nipype.interfaces.base import ( File, @@ -689,7 +690,7 @@ def reconstruct_fieldmap( ) if not direct: - fmap_img = transforms.apply(fmap_img, reference=target) + fmap_img = nt.resampling.apply(transforms, fmap_img, reference=target) fmap_img.header.set_intent('estimate', name='fieldmap Hz') fmap_img.header.set_data_dtype('float32') diff --git a/fmriprep/workflows/bold/registration.py b/fmriprep/workflows/bold/registration.py index 3f7b5e1d6..8455b1182 100644 --- a/fmriprep/workflows/bold/registration.py +++ b/fmriprep/workflows/bold/registration.py @@ -737,6 +737,7 @@ def _conditional_downsampling(in_file, in_mask, zoom_th=4.0): import nibabel as nb import nitransforms as nt import numpy as np + from nitransforms.resampling import apply as transform from scipy.ndimage.filters import gaussian_filter img = nb.load(in_file) @@ -756,14 +757,16 @@ def _conditional_downsampling(in_file, in_mask, zoom_th=4.0): offset = old_center - newrot.dot((newshape - 1) * 0.5) newaffine = nb.affines.from_matvec(newrot, offset) + identity = nt.base.TransformBase() + newref = nb.Nifti1Image(np.zeros(newshape, dtype=np.uint8), newaffine) - nt.Affine(reference=newref).apply(img).to_filename(out_file) + transform(identity, img, reference=newref).to_filename(out_file) mask = nb.load(in_mask) mask.set_data_dtype(float) mdata = gaussian_filter(mask.get_fdata(dtype=float), scaling) floatmask = nb.Nifti1Image(mdata, mask.affine, mask.header) - newmask = nt.Affine(reference=newref).apply(floatmask) + newmask = transform(identity, floatmask, reference=newref) hdr = newmask.header.copy() hdr.set_data_dtype(np.uint8) newmaskdata = (newmask.get_fdata(dtype=float) > 0.5).astype(np.uint8) From bd523b987b3ba38f653efdeee85342565e14b7c2 Mon Sep 17 00:00:00 2001 From: Chris Markiewicz Date: Wed, 9 Oct 2024 20:55:24 -0400 Subject: [PATCH 3/5] fix: Update requirements.txt --- requirements.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/requirements.txt b/requirements.txt index 77ca17846..6f7130ea3 100644 --- a/requirements.txt +++ b/requirements.txt @@ -232,7 +232,7 @@ nireports==23.2.2 # via fmriprep (pyproject.toml) nitime==0.10.2 # via fmriprep (pyproject.toml) -nitransforms==23.0.1 +nitransforms==24.0.2 # via # fmriprep (pyproject.toml) # niworkflows From 02d617ac2865d13d682a7900b0bf147288478098 Mon Sep 17 00:00:00 2001 From: Chris Markiewicz Date: Thu, 10 Oct 2024 09:49:15 -0400 Subject: [PATCH 4/5] fix: Use nt.Affine() for identity, not nt.base.TransformBase() --- fmriprep/utils/transforms.py | 2 +- fmriprep/workflows/bold/registration.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/fmriprep/utils/transforms.py b/fmriprep/utils/transforms.py index e15010954..b1049ebd2 100644 --- a/fmriprep/utils/transforms.py +++ b/fmriprep/utils/transforms.py @@ -30,5 +30,5 @@ def load_transforms(xfm_paths: list[Path], inverse: list[bool]) -> nt.base.Trans else: chain += xfm if chain is None: - chain = nt.base.TransformBase() + chain = nt.Affine() # Identity return chain diff --git a/fmriprep/workflows/bold/registration.py b/fmriprep/workflows/bold/registration.py index 8455b1182..8d9c52664 100644 --- a/fmriprep/workflows/bold/registration.py +++ b/fmriprep/workflows/bold/registration.py @@ -757,7 +757,7 @@ def _conditional_downsampling(in_file, in_mask, zoom_th=4.0): offset = old_center - newrot.dot((newshape - 1) * 0.5) newaffine = nb.affines.from_matvec(newrot, offset) - identity = nt.base.TransformBase() + identity = nt.Affine() newref = nb.Nifti1Image(np.zeros(newshape, dtype=np.uint8), newaffine) transform(identity, img, reference=newref).to_filename(out_file) From 7a7a623362460df93f6047e34cdaa30cdbf0ad53 Mon Sep 17 00:00:00 2001 From: Mathias Goncalves Date: Thu, 10 Oct 2024 10:51:57 -0400 Subject: [PATCH 5/5] Apply suggestions from code review Co-authored-by: Chris Markiewicz --- fmriprep/workflows/bold/registration.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/fmriprep/workflows/bold/registration.py b/fmriprep/workflows/bold/registration.py index 8d9c52664..c145edc2f 100644 --- a/fmriprep/workflows/bold/registration.py +++ b/fmriprep/workflows/bold/registration.py @@ -737,7 +737,7 @@ def _conditional_downsampling(in_file, in_mask, zoom_th=4.0): import nibabel as nb import nitransforms as nt import numpy as np - from nitransforms.resampling import apply as transform + from nitransforms.resampling import apply as applyxfm from scipy.ndimage.filters import gaussian_filter img = nb.load(in_file) @@ -760,13 +760,13 @@ def _conditional_downsampling(in_file, in_mask, zoom_th=4.0): identity = nt.Affine() newref = nb.Nifti1Image(np.zeros(newshape, dtype=np.uint8), newaffine) - transform(identity, img, reference=newref).to_filename(out_file) + applyxfm(identity, img, reference=newref).to_filename(out_file) mask = nb.load(in_mask) mask.set_data_dtype(float) mdata = gaussian_filter(mask.get_fdata(dtype=float), scaling) floatmask = nb.Nifti1Image(mdata, mask.affine, mask.header) - newmask = transform(identity, floatmask, reference=newref) + newmask = applyxfm(identity, floatmask, reference=newref) hdr = newmask.header.copy() hdr.set_data_dtype(np.uint8) newmaskdata = (newmask.get_fdata(dtype=float) > 0.5).astype(np.uint8)