From 553d7e058787bf0d798f205bace0bed5edcd5ba6 Mon Sep 17 00:00:00 2001 From: mathiasg Date: Mon, 3 Jun 2024 15:56:10 -0400 Subject: [PATCH 1/5] 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 7b36eb35..d7f0210a 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 44bf3708e78bd4a619233a91baa55686a0047782 Mon Sep 17 00:00:00 2001 From: Mathias Goncalves Date: Tue, 4 Jun 2024 12:00:01 -0400 Subject: [PATCH 2/5] 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 d7f0210a..cf53a1f3 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 a28212dfb688c065739d7ca9aa4aac6db7287ed7 Mon Sep 17 00:00:00 2001 From: Chris Markiewicz Date: Tue, 4 Jun 2024 12:26:47 -0400 Subject: [PATCH 3/5] 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 cf53a1f3..23105efa 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 bcc7d3e2c3bcf8635febafdf97ce5727d67b2cec Mon Sep 17 00:00:00 2001 From: mathiasg Date: Wed, 5 Jun 2024 15:08:05 -0400 Subject: [PATCH 4/5] FIX: Catch cases of unexpected resampling --- fmriprep/utils/transforms.py | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/fmriprep/utils/transforms.py b/fmriprep/utils/transforms.py index 23105efa..8937f5b2 100644 --- a/fmriprep/utils/transforms.py +++ b/fmriprep/utils/transforms.py @@ -71,6 +71,16 @@ def load_ants_h5(filename: Path) -> nt.base.TransformBase: fixed_params = transform2['TransformFixedParameters'][:] + spacing = fixed_params[6:9] + direction = fixed_params[9:] + + # Supported spacing + if not np.array_equal(spacing, np.array([1.0, 1.0, 1.0])): + raise ValueError(f'Unexpected spacing: {spacing}') + + if not np.array_equal(direction, np.diag([-1, -1, -1])): + raise ValueError(f'Asymmetric direction matrix: {direction}') + 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 From c9619a81bc6cf45eac05bc8171307ca189c453e0 Mon Sep 17 00:00:00 2001 From: Mathias Goncalves Date: Wed, 5 Jun 2024 15:19:06 -0400 Subject: [PATCH 5/5] Apply suggestions from code review Co-authored-by: Chris Markiewicz --- fmriprep/utils/transforms.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/fmriprep/utils/transforms.py b/fmriprep/utils/transforms.py index 8937f5b2..bcd96e5f 100644 --- a/fmriprep/utils/transforms.py +++ b/fmriprep/utils/transforms.py @@ -72,13 +72,13 @@ def load_ants_h5(filename: Path) -> nt.base.TransformBase: fixed_params = transform2['TransformFixedParameters'][:] spacing = fixed_params[6:9] - direction = fixed_params[9:] + direction = fixed_params[9:].reshape((3, 3)) # Supported spacing if not np.array_equal(spacing, np.array([1.0, 1.0, 1.0])): raise ValueError(f'Unexpected spacing: {spacing}') - if not np.array_equal(direction, np.diag([-1, -1, -1])): + if not np.array_equal(direction, direction.T): raise ValueError(f'Asymmetric direction matrix: {direction}') shape = tuple(fixed_params[:3].astype(int))