From b922fa5fe473d43d03f56afe2aff75fbe52a4f55 Mon Sep 17 00:00:00 2001 From: Oscar Esteban Date: Tue, 23 Jul 2024 12:56:36 +0200 Subject: [PATCH 01/20] wip: initiate implementation --- nitransforms/resampling.py | 18 +++++++++++++++--- 1 file changed, 15 insertions(+), 3 deletions(-) diff --git a/nitransforms/resampling.py b/nitransforms/resampling.py index 9de0d2d6..bc343231 100644 --- a/nitransforms/resampling.py +++ b/nitransforms/resampling.py @@ -7,6 +7,7 @@ # ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ## """Resampling utilities.""" +from warnings import warn from pathlib import Path import numpy as np from nibabel.loadsave import load as _nbload @@ -19,6 +20,9 @@ _as_homogeneous, ) +SERIALIZE_VOLUME_WINDOW_WIDTH : int = 8 +"""Minimum number of volumes to automatically serialize 4D transforms.""" + def apply( transform, @@ -29,6 +33,8 @@ def apply( cval=0.0, prefilter=True, output_dtype=None, + serialize_nvols=SERIALIZE_VOLUME_WINDOW_WIDTH, + njobs=None, ): """ Apply a transformation to an image, resampling on the reference spatial object. @@ -89,14 +95,20 @@ def apply( spatialimage = _nbload(str(spatialimage)) data = np.asanyarray(spatialimage.dataobj) + data_nvols = 1 if data.ndim < 4 else data.shape[-1] + xfm_nvols = len(transforms) - if data.ndim == 4 and data.shape[-1] != len(transform): + if data_nvols == 1 and xfm_nvols > 1: + data = data[..., np.newaxis] + elif data_nvols != xfm_nvols: raise ValueError( "The fourth dimension of the data does not match the tranform's shape." ) - if data.ndim < transform.ndim: - data = data[..., np.newaxis] + serialize_nvols = serialize_nvols if serialize_nvols and serialize_nvols > 1 else np.inf + serialize_4d = max(data_nvols, xfm_nvols) > serialize_nvols + if serialize_4d: + warn("4D transforms serialization into 3D+t not implemented") # For model-based nonlinear transforms, generate the corresponding dense field if hasattr(transform, "to_field") and callable(transform.to_field): From 6064b8c056c2797b1d6dad3ab4a4365054291982 Mon Sep 17 00:00:00 2001 From: Julien Marabotto Date: Wed, 24 Jul 2024 11:19:56 +0200 Subject: [PATCH 02/20] enh: draft implementation of serialize 4d --- nitransforms/resampling.py | 87 ++++++++++++++++++++++++++------------ 1 file changed, 60 insertions(+), 27 deletions(-) diff --git a/nitransforms/resampling.py b/nitransforms/resampling.py index bc343231..ad37c768 100644 --- a/nitransforms/resampling.py +++ b/nitransforms/resampling.py @@ -7,12 +7,13 @@ # ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ## """Resampling utilities.""" -from warnings import warn from pathlib import Path import numpy as np from nibabel.loadsave import load as _nbload +from nibabel.arrayproxy import get_obj_dtype from scipy import ndimage as ndi +from nitransforms.linear import Affine, get from nitransforms.base import ( ImageGrid, TransformError, @@ -96,45 +97,77 @@ def apply( data = np.asanyarray(spatialimage.dataobj) data_nvols = 1 if data.ndim < 4 else data.shape[-1] - xfm_nvols = len(transforms) + xfm_nvols = len(transform) + assert xfm_nvols == transform.ndim == _ref.ndim if data_nvols == 1 and xfm_nvols > 1: data = data[..., np.newaxis] elif data_nvols != xfm_nvols: raise ValueError( - "The fourth dimension of the data does not match the tranform's shape." + "The fourth dimension of the data does not match the transform's shape." ) serialize_nvols = serialize_nvols if serialize_nvols and serialize_nvols > 1 else np.inf serialize_4d = max(data_nvols, xfm_nvols) > serialize_nvols if serialize_4d: - warn("4D transforms serialization into 3D+t not implemented") - - # For model-based nonlinear transforms, generate the corresponding dense field - if hasattr(transform, "to_field") and callable(transform.to_field): - targets = ImageGrid(spatialimage).index( - _as_homogeneous( - transform.to_field(reference=reference).map(_ref.ndcoords.T), - dim=_ref.ndim, + for t, xfm_t in enumerate(transform): + ras2vox = ~Affine(spatialimage.affine) + input_dtype = get_obj_dtype(spatialimage.dataobj) + output_dtype = output_dtype or input_dtype + + # Map the input coordinates on to timepoint t of the target (moving) + xcoords = _ref.ndcoords.astype("f4").T + ycoords = xfm_t.map(xcoords)[..., : _ref.ndim] + + # Calculate corresponding voxel coordinates + yvoxels = ras2vox.map(ycoords)[..., : _ref.ndim] + + # Interpolate + dataobj = ( + np.asanyarray(spatialimage.dataobj, dtype=input_dtype) + if spatialimage.ndim in (2, 3) + else None ) - ) + resampled[..., t] = ndi.map_coordinates( + ( + dataobj + if dataobj is not None + else spatialimage.dataobj[..., t].astype(input_dtype, copy=False) + ), + yvoxels.T, + output=output_dtype, + order=order, + mode=mode, + cval=cval, + prefilter=prefilter, + ) + else: - targets = ImageGrid(spatialimage).index( # data should be an image - _as_homogeneous(transform.map(_ref.ndcoords.T), dim=_ref.ndim) - ) + # For model-based nonlinear transforms, generate the corresponding dense field + if hasattr(transform, "to_field") and callable(transform.to_field): + targets = ImageGrid(spatialimage).index( + _as_homogeneous( + transform.to_field(reference=reference).map(_ref.ndcoords.T), + dim=_ref.ndim, + ) + ) + else: + targets = ImageGrid(spatialimage).index( # data should be an image + _as_homogeneous(transform.map(_ref.ndcoords.T), dim=_ref.ndim) + ) - if transform.ndim == 4: - targets = _as_homogeneous(targets.reshape(-2, targets.shape[0])).T - - resampled = ndi.map_coordinates( - data, - targets, - output=output_dtype, - order=order, - mode=mode, - cval=cval, - prefilter=prefilter, - ) + if transform.ndim == 4: + targets = _as_homogeneous(targets.reshape(-2, targets.shape[0])).T + + resampled = ndi.map_coordinates( + data, + targets, + output=output_dtype, + order=order, + mode=mode, + cval=cval, + prefilter=prefilter, + ) if isinstance(_ref, ImageGrid): # If reference is grid, reshape hdr = None From e47a4769b03c351a8e907e380e3dffd74e3a2955 Mon Sep 17 00:00:00 2001 From: Julien Marabotto Date: Thu, 25 Jul 2024 09:34:44 +0200 Subject: [PATCH 03/20] fix: passes more tests, more suggestions in progress --- nitransforms/resampling.py | 18 +++++++++++++++--- nitransforms/tests/test_base.py | 3 ++- 2 files changed, 17 insertions(+), 4 deletions(-) diff --git a/nitransforms/resampling.py b/nitransforms/resampling.py index ad37c768..b9ca65b8 100644 --- a/nitransforms/resampling.py +++ b/nitransforms/resampling.py @@ -13,7 +13,7 @@ from nibabel.arrayproxy import get_obj_dtype from scipy import ndimage as ndi -from nitransforms.linear import Affine, get +from nitransforms.linear import Affine, LinearTransformsMapping from nitransforms.base import ( ImageGrid, TransformError, @@ -97,15 +97,27 @@ def apply( data = np.asanyarray(spatialimage.dataobj) data_nvols = 1 if data.ndim < 4 else data.shape[-1] - xfm_nvols = len(transform) - assert xfm_nvols == transform.ndim == _ref.ndim + if type(transform) == Affine or type(transform) == LinearTransformsMapping: + xfm_nvols = len(transform) + else: + xfm_nvols = transform.ndim + """ if data_nvols == 1 and xfm_nvols > 1: data = data[..., np.newaxis] elif data_nvols != xfm_nvols: raise ValueError( "The fourth dimension of the data does not match the transform's shape." ) + RESAMPLING FAILS. SUGGEST: + """ + if data.ndim < transform.ndim: + data = data[..., np.newaxis] + elif data_nvols > 1 and data_nvols != xfm_nvols: + import pdb; pdb.set_trace() + raise ValueError( + "The fourth dimension of the data does not match the transform's shape." + ) serialize_nvols = serialize_nvols if serialize_nvols and serialize_nvols > 1 else np.inf serialize_4d = max(data_nvols, xfm_nvols) > serialize_nvols diff --git a/nitransforms/tests/test_base.py b/nitransforms/tests/test_base.py index fb4be8d8..74bc3358 100644 --- a/nitransforms/tests/test_base.py +++ b/nitransforms/tests/test_base.py @@ -186,6 +186,7 @@ def test_SurfaceMesh(testdata_path): with pytest.raises(ValueError): SurfaceMesh(nb.load(img_path)) - + """ with pytest.raises(TypeError): SurfaceMesh(nb.load(shape_path)) + """ \ No newline at end of file From 1616a35bf454898a6ff95b4d2925b4496da5be81 Mon Sep 17 00:00:00 2001 From: Julien Marabotto Date: Thu, 25 Jul 2024 11:37:32 +0200 Subject: [PATCH 04/20] fix: pass tests --- nitransforms/resampling.py | 1 - nitransforms/tests/test_base.py | 3 +-- 2 files changed, 1 insertion(+), 3 deletions(-) diff --git a/nitransforms/resampling.py b/nitransforms/resampling.py index b9ca65b8..c36750ef 100644 --- a/nitransforms/resampling.py +++ b/nitransforms/resampling.py @@ -114,7 +114,6 @@ def apply( if data.ndim < transform.ndim: data = data[..., np.newaxis] elif data_nvols > 1 and data_nvols != xfm_nvols: - import pdb; pdb.set_trace() raise ValueError( "The fourth dimension of the data does not match the transform's shape." ) diff --git a/nitransforms/tests/test_base.py b/nitransforms/tests/test_base.py index 74bc3358..fb4be8d8 100644 --- a/nitransforms/tests/test_base.py +++ b/nitransforms/tests/test_base.py @@ -186,7 +186,6 @@ def test_SurfaceMesh(testdata_path): with pytest.raises(ValueError): SurfaceMesh(nb.load(img_path)) - """ + with pytest.raises(TypeError): SurfaceMesh(nb.load(shape_path)) - """ \ No newline at end of file From 6292daf1d0f7dc56ae51d1d87a83fe827f72dd5c Mon Sep 17 00:00:00 2001 From: Julien Marabotto Date: Thu, 25 Jul 2024 13:44:11 +0200 Subject: [PATCH 05/20] fix: pass tests, serialization implemented --- nitransforms/resampling.py | 44 +++++++++++++++++++++----------------- 1 file changed, 24 insertions(+), 20 deletions(-) diff --git a/nitransforms/resampling.py b/nitransforms/resampling.py index c36750ef..52f831ef 100644 --- a/nitransforms/resampling.py +++ b/nitransforms/resampling.py @@ -102,15 +102,7 @@ def apply( xfm_nvols = len(transform) else: xfm_nvols = transform.ndim - """ - if data_nvols == 1 and xfm_nvols > 1: - data = data[..., np.newaxis] - elif data_nvols != xfm_nvols: - raise ValueError( - "The fourth dimension of the data does not match the transform's shape." - ) - RESAMPLING FAILS. SUGGEST: - """ + if data.ndim < transform.ndim: data = data[..., np.newaxis] elif data_nvols > 1 and data_nvols != xfm_nvols: @@ -119,26 +111,38 @@ def apply( ) serialize_nvols = serialize_nvols if serialize_nvols and serialize_nvols > 1 else np.inf - serialize_4d = max(data_nvols, xfm_nvols) > serialize_nvols + serialize_4d = max(data_nvols, xfm_nvols) >= serialize_nvols + if serialize_4d: - for t, xfm_t in enumerate(transform): - ras2vox = ~Affine(spatialimage.affine) - input_dtype = get_obj_dtype(spatialimage.dataobj) - output_dtype = output_dtype or input_dtype + # Avoid opening the data array just yet + input_dtype = get_obj_dtype(spatialimage.dataobj) + output_dtype = output_dtype or input_dtype + + # Prepare physical coordinates of input (grid, points) + xcoords = _ref.ndcoords.astype("f4").T + + # Invert target's (moving) affine once + ras2vox = ~Affine(spatialimage.affine) + dataobj = ( + np.asanyarray(spatialimage.dataobj, dtype=input_dtype) + if spatialimage.ndim in (2, 3) + else None + ) + # Order F ensures individual volumes are contiguous in memory + # Also matches NIfTI, making final save more efficient + resampled = np.zeros( + (xcoords.shape[0], len(transform)), dtype=output_dtype, order="F" + ) + + for t, xfm_t in enumerate(transform): # Map the input coordinates on to timepoint t of the target (moving) - xcoords = _ref.ndcoords.astype("f4").T ycoords = xfm_t.map(xcoords)[..., : _ref.ndim] # Calculate corresponding voxel coordinates yvoxels = ras2vox.map(ycoords)[..., : _ref.ndim] # Interpolate - dataobj = ( - np.asanyarray(spatialimage.dataobj, dtype=input_dtype) - if spatialimage.ndim in (2, 3) - else None - ) resampled[..., t] = ndi.map_coordinates( ( dataobj From 86b3d111f7635a04b403f4eb5c39000fd7637e20 Mon Sep 17 00:00:00 2001 From: Oscar Esteban Date: Tue, 23 Jul 2024 12:56:36 +0200 Subject: [PATCH 06/20] wip: initiate implementation --- nitransforms/resampling.py | 25 ++++++++++++++++++++++--- 1 file changed, 22 insertions(+), 3 deletions(-) diff --git a/nitransforms/resampling.py b/nitransforms/resampling.py index 9de0d2d6..b54329a5 100644 --- a/nitransforms/resampling.py +++ b/nitransforms/resampling.py @@ -7,6 +7,7 @@ # ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ## """Resampling utilities.""" +from warnings import warn from pathlib import Path import numpy as np from nibabel.loadsave import load as _nbload @@ -19,6 +20,12 @@ _as_homogeneous, ) +SERIALIZE_VOLUME_WINDOW_WIDTH : int = 8 +"""Minimum number of volumes to automatically serialize 4D transforms.""" + +class NotImplementedWarning(UserWarning): + """A custom class for warnings.""" + def apply( transform, @@ -29,6 +36,8 @@ def apply( cval=0.0, prefilter=True, output_dtype=None, + serialize_nvols=SERIALIZE_VOLUME_WINDOW_WIDTH, + njobs=None, ): """ Apply a transformation to an image, resampling on the reference spatial object. @@ -89,14 +98,24 @@ def apply( spatialimage = _nbload(str(spatialimage)) data = np.asanyarray(spatialimage.dataobj) + data_nvols = 1 if data.ndim < 4 else data.shape[-1] + xfm_nvols = len(transform) - if data.ndim == 4 and data.shape[-1] != len(transform): + if data_nvols == 1 and xfm_nvols > 1: + data = data[..., np.newaxis] + elif data_nvols != xfm_nvols: raise ValueError( "The fourth dimension of the data does not match the tranform's shape." ) - if data.ndim < transform.ndim: - data = data[..., np.newaxis] + serialize_nvols = serialize_nvols if serialize_nvols and serialize_nvols > 1 else np.inf + serialize_4d = max(data_nvols, xfm_nvols) > serialize_nvols + if serialize_4d: + warn( + "4D transforms serialization into 3D+t not implemented", + NotImplementedWarning, + stacklevel=2, + ) # For model-based nonlinear transforms, generate the corresponding dense field if hasattr(transform, "to_field") and callable(transform.to_field): From 79e5cadc6bfec93982dd70e62f3e916c0a28ab78 Mon Sep 17 00:00:00 2001 From: Oscar Esteban Date: Tue, 30 Jul 2024 16:42:05 +0200 Subject: [PATCH 07/20] enh: integrating @jmarabotto's code --- nitransforms/resampling.py | 84 ++++++++++++++++++++------------------ 1 file changed, 45 insertions(+), 39 deletions(-) diff --git a/nitransforms/resampling.py b/nitransforms/resampling.py index 09801b8d..1d6e7f76 100644 --- a/nitransforms/resampling.py +++ b/nitransforms/resampling.py @@ -7,14 +7,13 @@ # ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ## """Resampling utilities.""" -from warnings import warn + from pathlib import Path import numpy as np from nibabel.loadsave import load as _nbload from nibabel.arrayproxy import get_obj_dtype from scipy import ndimage as ndi -from nitransforms.linear import Affine, LinearTransformsMapping from nitransforms.base import ( ImageGrid, TransformError, @@ -22,7 +21,7 @@ _as_homogeneous, ) -SERIALIZE_VOLUME_WINDOW_WIDTH : int = 8 +SERIALIZE_VOLUME_WINDOW_WIDTH: int = 8 """Minimum number of volumes to automatically serialize 4D transforms.""" @@ -96,58 +95,67 @@ def apply( if isinstance(spatialimage, (str, Path)): spatialimage = _nbload(str(spatialimage)) - data = np.asanyarray(spatialimage.dataobj) - data_nvols = 1 if data.ndim < 4 else data.shape[-1] + # Avoid opening the data array just yet + input_dtype = get_obj_dtype(spatialimage.dataobj) + output_dtype = output_dtype or input_dtype + # Number of transformations + data_nvols = 1 if spatialimage.ndim < 4 else spatialimage.shape[-1] xfm_nvols = len(transform) - if data_nvols == 1 and xfm_nvols > 1: - data = data[..., np.newaxis] - elif data_nvols != xfm_nvols: + if data_nvols != xfm_nvols and min(data_nvols, xfm_nvols) > 1: raise ValueError( "The fourth dimension of the data does not match the transform's shape." ) - serialize_nvols = serialize_nvols if serialize_nvols and serialize_nvols > 1 else np.inf - serialize_4d = max(data_nvols, xfm_nvols) >= serialize_nvols + serialize_nvols = ( + serialize_nvols if serialize_nvols and serialize_nvols > 1 else np.inf + ) + n_resamplings = max(data_nvols, xfm_nvols) + serialize_4d = n_resamplings >= serialize_nvols + + targets = None + if hasattr(transform, "to_field") and callable(transform.to_field): + targets = ImageGrid(spatialimage).index( + _as_homogeneous( + transform.to_field(reference=reference).map(_ref.ndcoords.T), + dim=_ref.ndim, + ) + ) + elif xfm_nvols == 1: + targets = ImageGrid(spatialimage).index( # data should be an image + _as_homogeneous(transform.map(_ref.ndcoords.T), dim=_ref.ndim) + ) if serialize_4d: - # Avoid opening the data array just yet - input_dtype = get_obj_dtype(spatialimage.dataobj) - output_dtype = output_dtype or input_dtype - - # Prepare physical coordinates of input (grid, points) - xcoords = _ref.ndcoords.astype("f4").T - - # Invert target's (moving) affine once - ras2vox = ~Affine(spatialimage.affine) - dataobj = ( + data = ( np.asanyarray(spatialimage.dataobj, dtype=input_dtype) - if spatialimage.ndim in (2, 3) + if data_nvols == 1 else None ) # Order F ensures individual volumes are contiguous in memory # Also matches NIfTI, making final save more efficient resampled = np.zeros( - (xcoords.shape[0], len(transform)), dtype=output_dtype, order="F" + (spatialimage.size, len(transform)), dtype=output_dtype, order="F" ) - for t, xfm_t in enumerate(transform): - # Map the input coordinates on to timepoint t of the target (moving) - ycoords = xfm_t.map(xcoords)[..., : _ref.ndim] + for t in range(n_resamplings): + xfm_t = transform if n_resamplings == 1 else transform[t] - # Calculate corresponding voxel coordinates - yvoxels = ras2vox.map(ycoords)[..., : _ref.ndim] + if targets is None: + targets = ImageGrid(spatialimage).index( # data should be an image + _as_homogeneous(xfm_t.map(_ref.ndcoords.T), dim=_ref.ndim) + ) # Interpolate resampled[..., t] = ndi.map_coordinates( ( - dataobj - if dataobj is not None + data + if data is not None else spatialimage.dataobj[..., t].astype(input_dtype, copy=False) ), - yvoxels.T, + targets, output=output_dtype, order=order, mode=mode, @@ -156,19 +164,17 @@ def apply( ) else: - # For model-based nonlinear transforms, generate the corresponding dense field - if hasattr(transform, "to_field") and callable(transform.to_field): - targets = ImageGrid(spatialimage).index( - _as_homogeneous( - transform.to_field(reference=reference).map(_ref.ndcoords.T), - dim=_ref.ndim, - ) - ) - else: + data = np.asanyarray(spatialimage.dataobj, dtype=input_dtype) + + if targets is None: targets = ImageGrid(spatialimage).index( # data should be an image _as_homogeneous(transform.map(_ref.ndcoords.T), dim=_ref.ndim) ) + # Cast 3D data into 4D if 4D nonsequential transform + if data_nvols == 1 and xfm_nvols > 1: + data = data[..., np.newaxis] + if transform.ndim == 4: targets = _as_homogeneous(targets.reshape(-2, targets.shape[0])).T From e0bde092d14e67491078823e4218aa2afcd35144 Mon Sep 17 00:00:00 2001 From: Oscar Esteban Date: Tue, 30 Jul 2024 16:57:55 +0200 Subject: [PATCH 08/20] fix: ensure output dtype when resampling --- nitransforms/resampling.py | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/nitransforms/resampling.py b/nitransforms/resampling.py index 1d6e7f76..45474008 100644 --- a/nitransforms/resampling.py +++ b/nitransforms/resampling.py @@ -189,10 +189,9 @@ def apply( ) if isinstance(_ref, ImageGrid): # If reference is grid, reshape - hdr = None - if _ref.header is not None: - hdr = _ref.header.copy() - hdr.set_data_dtype(output_dtype or spatialimage.get_data_dtype()) + hdr = _ref.header.copy() if _ref.header is not None else spatialimage.header.__class__() + hdr.set_data_dtype(output_dtype) + moved = spatialimage.__class__( resampled.reshape(_ref.shape if data.ndim < 4 else _ref.shape + (-1,)), _ref.affine, From fbb04511dd210df8a08101fe883e9ab140807e8b Mon Sep 17 00:00:00 2001 From: Oscar Esteban Date: Tue, 30 Jul 2024 18:59:34 +0200 Subject: [PATCH 09/20] fix: resolve some failing tests --- nitransforms/nonlinear.py | 4 ++++ nitransforms/resampling.py | 21 ++++++++++----------- nitransforms/tests/test_base.py | 7 +++++-- 3 files changed, 19 insertions(+), 13 deletions(-) diff --git a/nitransforms/nonlinear.py b/nitransforms/nonlinear.py index 9c29c53c..ced348a2 100644 --- a/nitransforms/nonlinear.py +++ b/nitransforms/nonlinear.py @@ -94,6 +94,10 @@ def __repr__(self): """Beautify the python representation.""" return f"<{self.__class__.__name__}[{self._field.shape[-1]}D] {self._field.shape[:3]}>" + def __len__(self): + """Enable len() -- for compatibility, only len == 1 is supported.""" + return 1 + @property def ndim(self): """Get the dimensions of the transform.""" diff --git a/nitransforms/resampling.py b/nitransforms/resampling.py index 45474008..eb3f9ad0 100644 --- a/nitransforms/resampling.py +++ b/nitransforms/resampling.py @@ -97,7 +97,6 @@ def apply( # Avoid opening the data array just yet input_dtype = get_obj_dtype(spatialimage.dataobj) - output_dtype = output_dtype or input_dtype # Number of transformations data_nvols = 1 if spatialimage.ndim < 4 else spatialimage.shape[-1] @@ -115,16 +114,17 @@ def apply( serialize_4d = n_resamplings >= serialize_nvols targets = None + ref_ndcoords = _ref.ndcoords.T if hasattr(transform, "to_field") and callable(transform.to_field): targets = ImageGrid(spatialimage).index( _as_homogeneous( - transform.to_field(reference=reference).map(_ref.ndcoords.T), + transform.to_field(reference=reference).map(ref_ndcoords), dim=_ref.ndim, ) ) elif xfm_nvols == 1: targets = ImageGrid(spatialimage).index( # data should be an image - _as_homogeneous(transform.map(_ref.ndcoords.T), dim=_ref.ndim) + _as_homogeneous(transform.map(ref_ndcoords), dim=_ref.ndim) ) if serialize_4d: @@ -137,7 +137,7 @@ def apply( # Order F ensures individual volumes are contiguous in memory # Also matches NIfTI, making final save more efficient resampled = np.zeros( - (spatialimage.size, len(transform)), dtype=output_dtype, order="F" + (len(ref_ndcoords), len(transform)), dtype=input_dtype, order="F" ) for t in range(n_resamplings): @@ -145,7 +145,7 @@ def apply( if targets is None: targets = ImageGrid(spatialimage).index( # data should be an image - _as_homogeneous(xfm_t.map(_ref.ndcoords.T), dim=_ref.ndim) + _as_homogeneous(xfm_t.map(ref_ndcoords), dim=_ref.ndim) ) # Interpolate @@ -156,7 +156,6 @@ def apply( else spatialimage.dataobj[..., t].astype(input_dtype, copy=False) ), targets, - output=output_dtype, order=order, mode=mode, cval=cval, @@ -168,7 +167,7 @@ def apply( if targets is None: targets = ImageGrid(spatialimage).index( # data should be an image - _as_homogeneous(transform.map(_ref.ndcoords.T), dim=_ref.ndim) + _as_homogeneous(transform.map(ref_ndcoords), dim=_ref.ndim) ) # Cast 3D data into 4D if 4D nonsequential transform @@ -181,7 +180,6 @@ def apply( resampled = ndi.map_coordinates( data, targets, - output=output_dtype, order=order, mode=mode, cval=cval, @@ -190,13 +188,14 @@ def apply( if isinstance(_ref, ImageGrid): # If reference is grid, reshape hdr = _ref.header.copy() if _ref.header is not None else spatialimage.header.__class__() - hdr.set_data_dtype(output_dtype) + hdr.set_data_dtype(output_dtype or spatialimage.header.get_data_dtype()) moved = spatialimage.__class__( - resampled.reshape(_ref.shape if data.ndim < 4 else _ref.shape + (-1,)), + resampled.reshape(_ref.shape if n_resamplings == 1 else _ref.shape + (-1,)), _ref.affine, hdr, ) return moved - return resampled + output_dtype = output_dtype or input_dtype + return resampled.astype(output_dtype) diff --git a/nitransforms/tests/test_base.py b/nitransforms/tests/test_base.py index fb4be8d8..c85ac2e2 100644 --- a/nitransforms/tests/test_base.py +++ b/nitransforms/tests/test_base.py @@ -1,6 +1,8 @@ """Tests of the base module.""" import numpy as np import nibabel as nb +from nibabel.arrayproxy import get_obj_dtype + import pytest import h5py @@ -97,7 +99,7 @@ def _to_hdf5(klass, x5_root): fname = testdata_path / "someones_anatomy.nii.gz" img = nb.load(fname) - imgdata = np.asanyarray(img.dataobj, dtype=img.get_data_dtype()) + imgdata = np.asanyarray(img.dataobj, dtype=get_obj_dtype(img.dataobj)) # Test identity transform - setting reference xfm = TransformBase() @@ -111,7 +113,8 @@ def _to_hdf5(klass, x5_root): xfm = nitl.Affine() xfm.reference = fname moved = apply(xfm, fname, order=0) - assert np.all(imgdata == np.asanyarray(moved.dataobj, dtype=moved.get_data_dtype())) + + assert np.all(imgdata == np.asanyarray(moved.dataobj, dtype=get_obj_dtype(moved.dataobj))) # Test ndim returned by affine assert nitl.Affine().ndim == 3 From 015347272558798f53b627d3b93e40159adba7b9 Mon Sep 17 00:00:00 2001 From: Oscar Esteban Date: Wed, 31 Jul 2024 07:25:47 +0200 Subject: [PATCH 10/20] fix: ensure ``__len__`` is defined for all transforms`` --- nitransforms/base.py | 10 ++++++++++ nitransforms/nonlinear.py | 4 ---- 2 files changed, 10 insertions(+), 4 deletions(-) diff --git a/nitransforms/base.py b/nitransforms/base.py index 81ed1a5e..a40998c5 100644 --- a/nitransforms/base.py +++ b/nitransforms/base.py @@ -279,6 +279,16 @@ def __add__(self, b): return TransformChain(transforms=[self, b]) + def __len__(self): + """ + Enable ``len()``. + + By default, all transforms are of length one. + This must be overriden by transforms arrays and chains. + + """ + return 1 + @property def reference(self): """Access a reference space where data will be resampled onto.""" diff --git a/nitransforms/nonlinear.py b/nitransforms/nonlinear.py index ced348a2..9c29c53c 100644 --- a/nitransforms/nonlinear.py +++ b/nitransforms/nonlinear.py @@ -94,10 +94,6 @@ def __repr__(self): """Beautify the python representation.""" return f"<{self.__class__.__name__}[{self._field.shape[-1]}D] {self._field.shape[:3]}>" - def __len__(self): - """Enable len() -- for compatibility, only len == 1 is supported.""" - return 1 - @property def ndim(self): """Get the dimensions of the transform.""" From 85d03b426d9987ceec535a337a647d719a77d298 Mon Sep 17 00:00:00 2001 From: Oscar Esteban Date: Wed, 31 Jul 2024 07:55:07 +0200 Subject: [PATCH 11/20] fix: clarify what is a 3D transform chain and a 4D transform 3D transform chains resulting of composing several transformations (e.g., affine and deformation fields in spatial normalization) should not be split into its components. This is in contrast to lists of 3D transforms such as head-motion correcting affines, where each applies to one timepoint. These should be considered 4D and in some future they may integrate slice timing correction in them. --- nitransforms/resampling.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/nitransforms/resampling.py b/nitransforms/resampling.py index eb3f9ad0..e2de9a2c 100644 --- a/nitransforms/resampling.py +++ b/nitransforms/resampling.py @@ -98,9 +98,10 @@ def apply( # Avoid opening the data array just yet input_dtype = get_obj_dtype(spatialimage.dataobj) - # Number of transformations + # Number of data volumes data_nvols = 1 if spatialimage.ndim < 4 else spatialimage.shape[-1] - xfm_nvols = len(transform) + # Number of transforms: transforms chains (e.g., affine + field, are a single transform) + xfm_nvols = 1 if transform.ndim < 4 else len(transform) if data_nvols != xfm_nvols and min(data_nvols, xfm_nvols) > 1: raise ValueError( From 06a1c01ba8492c9b0e6f08fe4e0f7758075c790f Mon Sep 17 00:00:00 2001 From: Oscar Esteban Date: Wed, 31 Jul 2024 08:43:59 +0200 Subject: [PATCH 12/20] maint: reorganize tests around the spun-off apply --- nitransforms/tests/test_linear.py | 146 ----------- nitransforms/tests/test_manip.py | 54 +--- nitransforms/tests/test_nonlinear.py | 142 ----------- nitransforms/tests/test_resampling.py | 352 ++++++++++++++++++++++++++ 4 files changed, 353 insertions(+), 341 deletions(-) create mode 100644 nitransforms/tests/test_resampling.py diff --git a/nitransforms/tests/test_linear.py b/nitransforms/tests/test_linear.py index 50cc5371..5746d5f7 100644 --- a/nitransforms/tests/test_linear.py +++ b/nitransforms/tests/test_linear.py @@ -4,8 +4,6 @@ import os import pytest import numpy as np -from subprocess import check_call -import shutil import h5py import nibabel as nb @@ -13,28 +11,8 @@ from nibabel.affines import from_matvec from nitransforms import linear as nitl from nitransforms import io -from nitransforms.resampling import apply from .utils import assert_affines_by_filename -RMSE_TOL = 0.1 -APPLY_LINEAR_CMD = { - "fsl": """\ -flirt -setbackground 0 -interp nearestneighbour -in {moving} -ref {reference} \ --applyxfm -init {transform} -out {resampled}\ -""".format, - "itk": """\ -antsApplyTransforms -d 3 -r {reference} -i {moving} \ --o {resampled} -n NearestNeighbor -t {transform} --float\ -""".format, - "afni": """\ -3dAllineate -base {reference} -input {moving} \ --prefix {resampled} -1Dmatrix_apply {transform} -final NN\ -""".format, - "fs": """\ -mri_vol2vol --mov {moving} --targ {reference} --lta {transform} \ ---o {resampled} --nearest""".format, -} - @pytest.mark.parametrize("matrix", [[0.0], np.ones((3, 3, 3)), np.ones((3, 4)), ]) def test_linear_typeerrors1(matrix): @@ -234,96 +212,6 @@ def test_linear_save(tmpdir, data_path, get_testdata, image_orientation, sw_tool assert_affines_by_filename(xfm_fname1, xfm_fname2) -@pytest.mark.parametrize("image_orientation", ["RAS", "LAS", "LPS", 'oblique', ]) -@pytest.mark.parametrize("sw_tool", ["itk", "fsl", "afni", "fs"]) -def test_apply_linear_transform(tmpdir, get_testdata, get_testmask, image_orientation, sw_tool): - """Check implementation of exporting affines to formats.""" - tmpdir.chdir() - - img = get_testdata[image_orientation] - msk = get_testmask[image_orientation] - - # Generate test transform - T = from_matvec(euler2mat(x=0.9, y=0.001, z=0.001), [4.0, 2.0, -1.0]) - xfm = nitl.Affine(T) - xfm.reference = img - - ext = "" - if sw_tool == "itk": - ext = ".tfm" - elif sw_tool == "fs": - ext = ".lta" - - img.to_filename("img.nii.gz") - msk.to_filename("mask.nii.gz") - - # Write out transform file (software-dependent) - xfm_fname = f"M.{sw_tool}{ext}" - # Change reference dataset for AFNI & oblique - if (sw_tool, image_orientation) == ("afni", "oblique"): - io.afni.AFNILinearTransform.from_ras( - T, - moving=img, - reference=img, - ).to_filename(xfm_fname) - else: - xfm.to_filename(xfm_fname, fmt=sw_tool) - - cmd = APPLY_LINEAR_CMD[sw_tool]( - transform=os.path.abspath(xfm_fname), - reference=os.path.abspath("mask.nii.gz"), - moving=os.path.abspath("mask.nii.gz"), - resampled=os.path.abspath("resampled_brainmask.nii.gz"), - ) - - # skip test if command is not available on host - exe = cmd.split(" ", 1)[0] - if not shutil.which(exe): - pytest.skip(f"Command {exe} not found on host") - - # resample mask - exit_code = check_call([cmd], shell=True) - assert exit_code == 0 - sw_moved_mask = nb.load("resampled_brainmask.nii.gz") - - nt_moved_mask = apply(xfm, msk, order=0) - nt_moved_mask.set_data_dtype(msk.get_data_dtype()) - nt_moved_mask.to_filename("ntmask.nii.gz") - diff = np.asanyarray(sw_moved_mask.dataobj) - np.asanyarray(nt_moved_mask.dataobj) - - assert np.sqrt((diff ** 2).mean()) < RMSE_TOL - brainmask = np.asanyarray(nt_moved_mask.dataobj, dtype=bool) - - cmd = APPLY_LINEAR_CMD[sw_tool]( - transform=os.path.abspath(xfm_fname), - reference=os.path.abspath("img.nii.gz"), - moving=os.path.abspath("img.nii.gz"), - resampled=os.path.abspath("resampled.nii.gz"), - ) - - exit_code = check_call([cmd], shell=True) - assert exit_code == 0 - sw_moved = nb.load("resampled.nii.gz") - sw_moved.set_data_dtype(img.get_data_dtype()) - - nt_moved = apply(xfm, img, order=0) - diff = ( - np.asanyarray(sw_moved.dataobj, dtype=sw_moved.get_data_dtype()) - - np.asanyarray(nt_moved.dataobj, dtype=nt_moved.get_data_dtype()) - ) - - # A certain tolerance is necessary because of resampling at borders - assert np.sqrt((diff[brainmask] ** 2).mean()) < RMSE_TOL - - nt_moved = apply(xfm, "img.nii.gz", order=0) - diff = ( - np.asanyarray(sw_moved.dataobj, dtype=sw_moved.get_data_dtype()) - - np.asanyarray(nt_moved.dataobj, dtype=nt_moved.get_data_dtype()) - ) - # A certain tolerance is necessary because of resampling at borders - assert np.sqrt((diff[brainmask] ** 2).mean()) < RMSE_TOL - - def test_Affine_to_x5(tmpdir, testdata_path): """Test affine's operations.""" tmpdir.chdir() @@ -336,40 +224,6 @@ def test_Affine_to_x5(tmpdir, testdata_path): aff._to_hdf5(f.create_group("Affine")) -def test_LinearTransformsMapping_apply(tmp_path, data_path, testdata_path): - """Apply transform mappings.""" - hmc = nitl.load( - data_path / "hmc-itk.tfm", fmt="itk", reference=testdata_path / "sbref.nii.gz" - ) - assert isinstance(hmc, nitl.LinearTransformsMapping) - - # Test-case: realign functional data on to sbref - nii = apply( - hmc, testdata_path / "func.nii.gz", order=1, reference=testdata_path / "sbref.nii.gz" - ) - assert nii.dataobj.shape[-1] == len(hmc) - - # Test-case: write out a fieldmap moved with head - hmcinv = nitl.LinearTransformsMapping( - np.linalg.inv(hmc.matrix), reference=testdata_path / "func.nii.gz" - ) - - nii = apply( - hmcinv, testdata_path / "fmap.nii.gz", order=1 - ) - assert nii.dataobj.shape[-1] == len(hmc) - - # Ensure a ValueError is issued when trying to do weird stuff - hmc = nitl.LinearTransformsMapping(hmc.matrix[:1, ...]) - with pytest.raises(ValueError): - apply( - hmc, - testdata_path / "func.nii.gz", - order=1, - reference=testdata_path / "sbref.nii.gz", - ) - - def test_mulmat_operator(testdata_path): """Check the @ operator.""" ref = testdata_path / "someones_anatomy.nii.gz" diff --git a/nitransforms/tests/test_manip.py b/nitransforms/tests/test_manip.py index b7f6a6e4..2a2d6ffb 100644 --- a/nitransforms/tests/test_manip.py +++ b/nitransforms/tests/test_manip.py @@ -1,67 +1,15 @@ # emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*- # vi: set ft=python sts=4 ts=4 sw=4 et: """Tests of nonlinear transforms.""" -import os -import shutil -from subprocess import check_call import pytest import numpy as np -import nibabel as nb -from ..manip import load as _load, TransformChain +from ..manip import TransformChain from ..linear import Affine -from .test_nonlinear import ( - RMSE_TOL, - APPLY_NONLINEAR_CMD, -) -from nitransforms.resampling import apply FMT = {"lta": "fs", "tfm": "itk"} -def test_itk_h5(tmp_path, testdata_path): - """Check a translation-only field on one or more axes, different image orientations.""" - os.chdir(str(tmp_path)) - img_fname = testdata_path / "T1w_scanner.nii.gz" - xfm_fname = ( - testdata_path - / "ds-005_sub-01_from-T1w_to-MNI152NLin2009cAsym_mode-image_xfm.h5" - ) - - xfm = _load(xfm_fname) - - assert len(xfm) == 2 - - ref_fname = tmp_path / "reference.nii.gz" - nb.Nifti1Image( - np.zeros(xfm.reference.shape, dtype="uint16"), xfm.reference.affine, - ).to_filename(str(ref_fname)) - - # Then apply the transform and cross-check with software - cmd = APPLY_NONLINEAR_CMD["itk"]( - transform=xfm_fname, - reference=ref_fname, - moving=img_fname, - output="resampled.nii.gz", - extra="", - ) - - # skip test if command is not available on host - exe = cmd.split(" ", 1)[0] - if not shutil.which(exe): - pytest.skip(f"Command {exe} not found on host") - - exit_code = check_call([cmd], shell=True) - assert exit_code == 0 - sw_moved = nb.load("resampled.nii.gz") - - nt_moved = apply(xfm, img_fname, order=0) - nt_moved.to_filename("nt_resampled.nii.gz") - diff = sw_moved.get_fdata() - nt_moved.get_fdata() - # A certain tolerance is necessary because of resampling at borders - assert (np.abs(diff) > 1e-3).sum() / diff.size < RMSE_TOL - - @pytest.mark.parametrize("ext0", ["lta", "tfm"]) @pytest.mark.parametrize("ext1", ["lta", "tfm"]) @pytest.mark.parametrize("ext2", ["lta", "tfm"]) diff --git a/nitransforms/tests/test_nonlinear.py b/nitransforms/tests/test_nonlinear.py index 24d1f83e..43b4584f 100644 --- a/nitransforms/tests/test_nonlinear.py +++ b/nitransforms/tests/test_nonlinear.py @@ -19,22 +19,6 @@ from ..io.itk import ITKDisplacementsField -RMSE_TOL = 0.05 -APPLY_NONLINEAR_CMD = { - "itk": """\ -antsApplyTransforms -d 3 -r {reference} -i {moving} \ --o {output} -n NearestNeighbor -t {transform} {extra}\ -""".format, - "afni": """\ -3dNwarpApply -nwarp {transform} -source {moving} \ --master {reference} -interp NN -prefix {output} {extra}\ -""".format, - "fsl": """\ -applywarp -i {moving} -r {reference} -o {output} {extra}\ --w {transform} --interp=nn""".format, -} - - @pytest.mark.parametrize("size", [(20, 20, 20), (20, 20, 20, 3)]) def test_itk_disp_load(size): """Checks field sizes.""" @@ -113,132 +97,6 @@ def test_bsplines_references(testdata_path): ) -@pytest.mark.parametrize("image_orientation", ["RAS", "LAS", "LPS", "oblique"]) -@pytest.mark.parametrize("sw_tool", ["itk", "afni"]) -@pytest.mark.parametrize("axis", [0, 1, 2, (0, 1), (1, 2), (0, 1, 2)]) -def test_displacements_field1( - tmp_path, - get_testdata, - get_testmask, - image_orientation, - sw_tool, - axis, -): - """Check a translation-only field on one or more axes, different image orientations.""" - if (image_orientation, sw_tool) == ("oblique", "afni"): - pytest.skip("AFNI obliques are not yet implemented for displacements fields") - - os.chdir(str(tmp_path)) - nii = get_testdata[image_orientation] - msk = get_testmask[image_orientation] - nii.to_filename("reference.nii.gz") - msk.to_filename("mask.nii.gz") - - fieldmap = np.zeros( - (*nii.shape[:3], 1, 3) if sw_tool != "fsl" else (*nii.shape[:3], 3), - dtype="float32", - ) - fieldmap[..., axis] = -10.0 - - _hdr = nii.header.copy() - if sw_tool in ("itk",): - _hdr.set_intent("vector") - _hdr.set_data_dtype("float32") - - xfm_fname = "warp.nii.gz" - field = nb.Nifti1Image(fieldmap, nii.affine, _hdr) - field.to_filename(xfm_fname) - - xfm = nlload(xfm_fname, fmt=sw_tool) - - # Then apply the transform and cross-check with software - cmd = APPLY_NONLINEAR_CMD[sw_tool]( - transform=os.path.abspath(xfm_fname), - reference=tmp_path / "mask.nii.gz", - moving=tmp_path / "mask.nii.gz", - output=tmp_path / "resampled_brainmask.nii.gz", - extra="--output-data-type uchar" if sw_tool == "itk" else "", - ) - - # skip test if command is not available on host - exe = cmd.split(" ", 1)[0] - if not shutil.which(exe): - pytest.skip(f"Command {exe} not found on host") - - # resample mask - exit_code = check_call([cmd], shell=True) - assert exit_code == 0 - sw_moved_mask = nb.load("resampled_brainmask.nii.gz") - nt_moved_mask = apply(xfm, msk, order=0) - nt_moved_mask.set_data_dtype(msk.get_data_dtype()) - diff = np.asanyarray(sw_moved_mask.dataobj) - np.asanyarray(nt_moved_mask.dataobj) - - assert np.sqrt((diff**2).mean()) < RMSE_TOL - brainmask = np.asanyarray(nt_moved_mask.dataobj, dtype=bool) - - # Then apply the transform and cross-check with software - cmd = APPLY_NONLINEAR_CMD[sw_tool]( - transform=os.path.abspath(xfm_fname), - reference=tmp_path / "reference.nii.gz", - moving=tmp_path / "reference.nii.gz", - output=tmp_path / "resampled.nii.gz", - extra="--output-data-type uchar" if sw_tool == "itk" else "", - ) - - exit_code = check_call([cmd], shell=True) - assert exit_code == 0 - sw_moved = nb.load("resampled.nii.gz") - - nt_moved = apply(xfm, nii, order=0) - nt_moved.set_data_dtype(nii.get_data_dtype()) - nt_moved.to_filename("nt_resampled.nii.gz") - sw_moved.set_data_dtype(nt_moved.get_data_dtype()) - diff = np.asanyarray( - sw_moved.dataobj, dtype=sw_moved.get_data_dtype() - ) - np.asanyarray(nt_moved.dataobj, dtype=nt_moved.get_data_dtype()) - # A certain tolerance is necessary because of resampling at borders - assert np.sqrt((diff[brainmask] ** 2).mean()) < RMSE_TOL - - -@pytest.mark.parametrize("sw_tool", ["itk", "afni"]) -def test_displacements_field2(tmp_path, testdata_path, sw_tool): - """Check a translation-only field on one or more axes, different image orientations.""" - os.chdir(str(tmp_path)) - img_fname = testdata_path / "tpl-OASIS30ANTs_T1w.nii.gz" - xfm_fname = testdata_path / "ds-005_sub-01_from-OASIS_to-T1_warp_{}.nii.gz".format( - sw_tool - ) - - xfm = nlload(xfm_fname, fmt=sw_tool) - - # Then apply the transform and cross-check with software - cmd = APPLY_NONLINEAR_CMD[sw_tool]( - transform=xfm_fname, - reference=img_fname, - moving=img_fname, - output="resampled.nii.gz", - extra="", - ) - - # skip test if command is not available on host - exe = cmd.split(" ", 1)[0] - if not shutil.which(exe): - pytest.skip(f"Command {exe} not found on host") - - exit_code = check_call([cmd], shell=True) - assert exit_code == 0 - sw_moved = nb.load("resampled.nii.gz") - - nt_moved = apply(xfm, img_fname, order=0) - nt_moved.to_filename("nt_resampled.nii.gz") - sw_moved.set_data_dtype(nt_moved.get_data_dtype()) - diff = np.asanyarray( - sw_moved.dataobj, dtype=sw_moved.get_data_dtype() - ) - np.asanyarray(nt_moved.dataobj, dtype=nt_moved.get_data_dtype()) - # A certain tolerance is necessary because of resampling at borders - assert np.sqrt((diff**2).mean()) < RMSE_TOL - - def test_bspline(tmp_path, testdata_path): """Cross-check B-Splines and deformation field.""" os.chdir(str(tmp_path)) diff --git a/nitransforms/tests/test_resampling.py b/nitransforms/tests/test_resampling.py new file mode 100644 index 00000000..3dd9aff4 --- /dev/null +++ b/nitransforms/tests/test_resampling.py @@ -0,0 +1,352 @@ +# emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*- +# vi: set ft=python sts=4 ts=4 sw=4 et: +"""Exercise the standalone ``apply()`` implementation.""" +import os +import pytest +import numpy as np +from subprocess import check_call +import shutil + +import nibabel as nb +from nibabel.eulerangles import euler2mat +from nibabel.affines import from_matvec +from nitransforms import linear as nitl +from nitransforms import nonlinear as nitnl +from nitransforms import manip as nitm +from nitransforms import io +from nitransforms.resampling import apply + +RMSE_TOL_LINEAR = 0.09 +RMSE_TOL_NONLINEAR = 0.05 +APPLY_LINEAR_CMD = { + "fsl": """\ +flirt -setbackground 0 -interp nearestneighbour -in {moving} -ref {reference} \ +-applyxfm -init {transform} -out {resampled}\ +""".format, + "itk": """\ +antsApplyTransforms -d 3 -r {reference} -i {moving} \ +-o {resampled} -n NearestNeighbor -t {transform} --float\ +""".format, + "afni": """\ +3dAllineate -base {reference} -input {moving} \ +-prefix {resampled} -1Dmatrix_apply {transform} -final NN\ +""".format, + "fs": """\ +mri_vol2vol --mov {moving} --targ {reference} --lta {transform} \ +--o {resampled} --nearest""".format, +} +APPLY_NONLINEAR_CMD = { + "itk": """\ +antsApplyTransforms -d 3 -r {reference} -i {moving} \ +-o {output} -n NearestNeighbor -t {transform} {extra}\ +""".format, + "afni": """\ +3dNwarpApply -nwarp {transform} -source {moving} \ +-master {reference} -interp NN -prefix {output} {extra}\ +""".format, + "fsl": """\ +applywarp -i {moving} -r {reference} -o {output} {extra}\ +-w {transform} --interp=nn""".format, +} + + +@pytest.mark.parametrize("image_orientation", ["RAS", "LAS", "LPS", 'oblique', ]) +@pytest.mark.parametrize("sw_tool", ["itk", "fsl", "afni", "fs"]) +def test_apply_linear_transform(tmpdir, get_testdata, get_testmask, image_orientation, sw_tool): + """Check implementation of exporting affines to formats.""" + tmpdir.chdir() + + img = get_testdata[image_orientation] + msk = get_testmask[image_orientation] + + # Generate test transform + T = from_matvec(euler2mat(x=0.9, y=0.001, z=0.001), [4.0, 2.0, -1.0]) + xfm = nitl.Affine(T) + xfm.reference = img + + ext = "" + if sw_tool == "itk": + ext = ".tfm" + elif sw_tool == "fs": + ext = ".lta" + + img.to_filename("img.nii.gz") + msk.to_filename("mask.nii.gz") + + # Write out transform file (software-dependent) + xfm_fname = f"M.{sw_tool}{ext}" + # Change reference dataset for AFNI & oblique + if (sw_tool, image_orientation) == ("afni", "oblique"): + io.afni.AFNILinearTransform.from_ras( + T, + moving=img, + reference=img, + ).to_filename(xfm_fname) + else: + xfm.to_filename(xfm_fname, fmt=sw_tool) + + cmd = APPLY_LINEAR_CMD[sw_tool]( + transform=os.path.abspath(xfm_fname), + reference=os.path.abspath("mask.nii.gz"), + moving=os.path.abspath("mask.nii.gz"), + resampled=os.path.abspath("resampled_brainmask.nii.gz"), + ) + + # skip test if command is not available on host + exe = cmd.split(" ", 1)[0] + if not shutil.which(exe): + pytest.skip(f"Command {exe} not found on host") + + # resample mask + exit_code = check_call([cmd], shell=True) + assert exit_code == 0 + sw_moved_mask = nb.load("resampled_brainmask.nii.gz") + + nt_moved_mask = apply(xfm, msk, order=0) + nt_moved_mask.set_data_dtype(msk.get_data_dtype()) + nt_moved_mask.to_filename("ntmask.nii.gz") + diff = np.asanyarray(sw_moved_mask.dataobj) - np.asanyarray(nt_moved_mask.dataobj) + + assert np.sqrt((diff ** 2).mean()) < RMSE_TOL_LINEAR + brainmask = np.asanyarray(nt_moved_mask.dataobj, dtype=bool) + + cmd = APPLY_LINEAR_CMD[sw_tool]( + transform=os.path.abspath(xfm_fname), + reference=os.path.abspath("img.nii.gz"), + moving=os.path.abspath("img.nii.gz"), + resampled=os.path.abspath("resampled.nii.gz"), + ) + + exit_code = check_call([cmd], shell=True) + assert exit_code == 0 + sw_moved = nb.load("resampled.nii.gz") + sw_moved.set_data_dtype(img.get_data_dtype()) + + nt_moved = apply(xfm, img, order=0) + diff = ( + np.asanyarray(sw_moved.dataobj, dtype=sw_moved.get_data_dtype()) + - np.asanyarray(nt_moved.dataobj, dtype=nt_moved.get_data_dtype()) + ) + + # A certain tolerance is necessary because of resampling at borders + assert np.sqrt((diff[brainmask] ** 2).mean()) < RMSE_TOL_LINEAR + + nt_moved = apply(xfm, "img.nii.gz", order=0) + diff = ( + np.asanyarray(sw_moved.dataobj, dtype=sw_moved.get_data_dtype()) + - np.asanyarray(nt_moved.dataobj, dtype=nt_moved.get_data_dtype()) + ) + # A certain tolerance is necessary because of resampling at borders + assert np.sqrt((diff[brainmask] ** 2).mean()) < RMSE_TOL_LINEAR + + +@pytest.mark.parametrize("image_orientation", ["RAS", "LAS", "LPS", "oblique"]) +@pytest.mark.parametrize("sw_tool", ["itk", "afni"]) +@pytest.mark.parametrize("axis", [0, 1, 2, (0, 1), (1, 2), (0, 1, 2)]) +def test_displacements_field1( + tmp_path, + get_testdata, + get_testmask, + image_orientation, + sw_tool, + axis, +): + """Check a translation-only field on one or more axes, different image orientations.""" + if (image_orientation, sw_tool) == ("oblique", "afni"): + pytest.skip("AFNI obliques are not yet implemented for displacements fields") + + os.chdir(str(tmp_path)) + nii = get_testdata[image_orientation] + msk = get_testmask[image_orientation] + nii.to_filename("reference.nii.gz") + msk.to_filename("mask.nii.gz") + + fieldmap = np.zeros( + (*nii.shape[:3], 1, 3) if sw_tool != "fsl" else (*nii.shape[:3], 3), + dtype="float32", + ) + fieldmap[..., axis] = -10.0 + + _hdr = nii.header.copy() + if sw_tool in ("itk",): + _hdr.set_intent("vector") + _hdr.set_data_dtype("float32") + + xfm_fname = "warp.nii.gz" + field = nb.Nifti1Image(fieldmap, nii.affine, _hdr) + field.to_filename(xfm_fname) + + xfm = nitnl.load(xfm_fname, fmt=sw_tool) + + # Then apply the transform and cross-check with software + cmd = APPLY_NONLINEAR_CMD[sw_tool]( + transform=os.path.abspath(xfm_fname), + reference=tmp_path / "mask.nii.gz", + moving=tmp_path / "mask.nii.gz", + output=tmp_path / "resampled_brainmask.nii.gz", + extra="--output-data-type uchar" if sw_tool == "itk" else "", + ) + + # skip test if command is not available on host + exe = cmd.split(" ", 1)[0] + if not shutil.which(exe): + pytest.skip(f"Command {exe} not found on host") + + # resample mask + exit_code = check_call([cmd], shell=True) + assert exit_code == 0 + sw_moved_mask = nb.load("resampled_brainmask.nii.gz") + nt_moved_mask = apply(xfm, msk, order=0) + nt_moved_mask.set_data_dtype(msk.get_data_dtype()) + diff = np.asanyarray(sw_moved_mask.dataobj) - np.asanyarray(nt_moved_mask.dataobj) + + assert np.sqrt((diff**2).mean()) < RMSE_TOL_LINEAR + brainmask = np.asanyarray(nt_moved_mask.dataobj, dtype=bool) + + # Then apply the transform and cross-check with software + cmd = APPLY_NONLINEAR_CMD[sw_tool]( + transform=os.path.abspath(xfm_fname), + reference=tmp_path / "reference.nii.gz", + moving=tmp_path / "reference.nii.gz", + output=tmp_path / "resampled.nii.gz", + extra="--output-data-type uchar" if sw_tool == "itk" else "", + ) + + exit_code = check_call([cmd], shell=True) + assert exit_code == 0 + sw_moved = nb.load("resampled.nii.gz") + + nt_moved = apply(xfm, nii, order=0) + nt_moved.set_data_dtype(nii.get_data_dtype()) + nt_moved.to_filename("nt_resampled.nii.gz") + sw_moved.set_data_dtype(nt_moved.get_data_dtype()) + diff = np.asanyarray( + sw_moved.dataobj, dtype=sw_moved.get_data_dtype() + ) - np.asanyarray(nt_moved.dataobj, dtype=nt_moved.get_data_dtype()) + # A certain tolerance is necessary because of resampling at borders + assert np.sqrt((diff[brainmask] ** 2).mean()) < RMSE_TOL_LINEAR + + +@pytest.mark.parametrize("sw_tool", ["itk", "afni"]) +def test_displacements_field2(tmp_path, testdata_path, sw_tool): + """Check a translation-only field on one or more axes, different image orientations.""" + os.chdir(str(tmp_path)) + img_fname = testdata_path / "tpl-OASIS30ANTs_T1w.nii.gz" + xfm_fname = testdata_path / "ds-005_sub-01_from-OASIS_to-T1_warp_{}.nii.gz".format( + sw_tool + ) + + xfm = nitnl.load(xfm_fname, fmt=sw_tool) + + # Then apply the transform and cross-check with software + cmd = APPLY_NONLINEAR_CMD[sw_tool]( + transform=xfm_fname, + reference=img_fname, + moving=img_fname, + output="resampled.nii.gz", + extra="", + ) + + # skip test if command is not available on host + exe = cmd.split(" ", 1)[0] + if not shutil.which(exe): + pytest.skip(f"Command {exe} not found on host") + + exit_code = check_call([cmd], shell=True) + assert exit_code == 0 + sw_moved = nb.load("resampled.nii.gz") + + nt_moved = apply(xfm, img_fname, order=0) + nt_moved.to_filename("nt_resampled.nii.gz") + sw_moved.set_data_dtype(nt_moved.get_data_dtype()) + diff = np.asanyarray( + sw_moved.dataobj, dtype=sw_moved.get_data_dtype() + ) - np.asanyarray(nt_moved.dataobj, dtype=nt_moved.get_data_dtype()) + # A certain tolerance is necessary because of resampling at borders + assert np.sqrt((diff**2).mean()) < RMSE_TOL_LINEAR + + +def test_apply_transformchain(tmp_path, testdata_path): + """Check a translation-only field on one or more axes, different image orientations.""" + os.chdir(str(tmp_path)) + img_fname = testdata_path / "T1w_scanner.nii.gz" + xfm_fname = ( + testdata_path + / "ds-005_sub-01_from-T1w_to-MNI152NLin2009cAsym_mode-image_xfm.h5" + ) + + xfm = nitm.load(xfm_fname) + + assert len(xfm) == 2 + + ref_fname = tmp_path / "reference.nii.gz" + nb.Nifti1Image( + np.zeros(xfm.reference.shape, dtype="uint16"), xfm.reference.affine, + ).to_filename(str(ref_fname)) + + # Then apply the transform and cross-check with software + cmd = APPLY_NONLINEAR_CMD["itk"]( + transform=xfm_fname, + reference=ref_fname, + moving=img_fname, + output="resampled.nii.gz", + extra="", + ) + + # skip test if command is not available on host + exe = cmd.split(" ", 1)[0] + if not shutil.which(exe): + pytest.skip(f"Command {exe} not found on host") + + exit_code = check_call([cmd], shell=True) + assert exit_code == 0 + sw_moved = nb.load("resampled.nii.gz") + + nt_moved = apply(xfm, img_fname, order=0) + nt_moved.to_filename("nt_resampled.nii.gz") + diff = sw_moved.get_fdata() - nt_moved.get_fdata() + # A certain tolerance is necessary because of resampling at borders + assert (np.abs(diff) > 1e-3).sum() / diff.size < RMSE_TOL_LINEAR + + +@pytest.mark.parametrize("serialize_4d", [True, False]) +def test_LinearTransformsMapping_apply(tmp_path, data_path, testdata_path, serialize_4d): + """Apply transform mappings.""" + hmc = nitl.load( + data_path / "hmc-itk.tfm", fmt="itk", reference=testdata_path / "sbref.nii.gz" + ) + assert isinstance(hmc, nitl.LinearTransformsMapping) + + # Test-case: realign functional data on to sbref + nii = apply( + hmc, + testdata_path / "func.nii.gz", + order=1, + reference=testdata_path / "sbref.nii.gz", + serialize_nvols=2 if serialize_4d else np.inf, + ) + assert nii.dataobj.shape[-1] == len(hmc) + + # Test-case: write out a fieldmap moved with head + hmcinv = nitl.LinearTransformsMapping( + np.linalg.inv(hmc.matrix), reference=testdata_path / "func.nii.gz" + ) + + nii = apply( + hmcinv, testdata_path / "fmap.nii.gz", + order=1, + serialize_nvols=2 if serialize_4d else np.inf, + ) + assert nii.dataobj.shape[-1] == len(hmc) + + # Ensure a ValueError is issued when trying to apply mismatched transforms + # (e.g., in this case, two transforms while the functional has 8 volumes) + hmc = nitl.LinearTransformsMapping(hmc.matrix[:2, ...]) + with pytest.raises(ValueError): + apply( + hmc, + testdata_path / "func.nii.gz", + order=1, + reference=testdata_path / "sbref.nii.gz", + serialize_nvols=2 if serialize_4d else np.inf, + ) From 8dd883dcd6fb91bd2dfc101620e9711301e3dc5f Mon Sep 17 00:00:00 2001 From: Oscar Esteban Date: Wed, 31 Jul 2024 09:15:44 +0200 Subject: [PATCH 13/20] sty: format changed files --- nitransforms/base.py | 17 ++++++----- nitransforms/resampling.py | 6 +++- nitransforms/tests/test_base.py | 10 +++++-- nitransforms/tests/test_linear.py | 41 +++++++++++++++++++-------- nitransforms/tests/test_manip.py | 1 + nitransforms/tests/test_nonlinear.py | 4 +-- nitransforms/tests/test_resampling.py | 41 ++++++++++++++++++--------- 7 files changed, 79 insertions(+), 41 deletions(-) diff --git a/nitransforms/base.py b/nitransforms/base.py index a40998c5..67acc073 100644 --- a/nitransforms/base.py +++ b/nitransforms/base.py @@ -7,6 +7,7 @@ # ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ## """Common interface for transforms.""" + from pathlib import Path import numpy as np import h5py @@ -146,13 +147,13 @@ def from_arrays(cls, coordinates, triangles): darrays = [ nb.gifti.GiftiDataArray( coordinates.astype(np.float32), - intent=nb.nifti1.intent_codes['NIFTI_INTENT_POINTSET'], - datatype=nb.nifti1.data_type_codes['NIFTI_TYPE_FLOAT32'], + intent=nb.nifti1.intent_codes["NIFTI_INTENT_POINTSET"], + datatype=nb.nifti1.data_type_codes["NIFTI_TYPE_FLOAT32"], ), nb.gifti.GiftiDataArray( triangles.astype(np.int32), - intent=nb.nifti1.intent_codes['NIFTI_INTENT_TRIANGLE'], - datatype=nb.nifti1.data_type_codes['NIFTI_TYPE_INT32'], + intent=nb.nifti1.intent_codes["NIFTI_INTENT_TRIANGLE"], + datatype=nb.nifti1.data_type_codes["NIFTI_TYPE_INT32"], ), ] gii = nb.gifti.GiftiImage(darrays=darrays) @@ -282,7 +283,7 @@ def __add__(self, b): def __len__(self): """ Enable ``len()``. - + By default, all transforms are of length one. This must be overriden by transforms arrays and chains. @@ -345,10 +346,8 @@ def apply(self, *args, **kwargs): Deprecated. Please use ``nitransforms.resampling.apply`` instead. """ - message = ( - "The `apply` method is deprecated. Please use `nitransforms.resampling.apply` instead." - ) - warnings.warn(message, DeprecationWarning, stacklevel=2) + _msg = "This method is deprecated. Please use `nitransforms.resampling.apply` instead." + warnings.warn(_msg, DeprecationWarning, stacklevel=2) from .resampling import apply return apply(self, *args, **kwargs) diff --git a/nitransforms/resampling.py b/nitransforms/resampling.py index e2de9a2c..abfe2b71 100644 --- a/nitransforms/resampling.py +++ b/nitransforms/resampling.py @@ -188,7 +188,11 @@ def apply( ) if isinstance(_ref, ImageGrid): # If reference is grid, reshape - hdr = _ref.header.copy() if _ref.header is not None else spatialimage.header.__class__() + hdr = ( + _ref.header.copy() + if _ref.header is not None + else spatialimage.header.__class__() + ) hdr.set_data_dtype(output_dtype or spatialimage.header.get_data_dtype()) moved = spatialimage.__class__( diff --git a/nitransforms/tests/test_base.py b/nitransforms/tests/test_base.py index c85ac2e2..4bb147fd 100644 --- a/nitransforms/tests/test_base.py +++ b/nitransforms/tests/test_base.py @@ -1,4 +1,5 @@ """Tests of the base module.""" + import numpy as np import nibabel as nb from nibabel.arrayproxy import get_obj_dtype @@ -114,7 +115,9 @@ def _to_hdf5(klass, x5_root): xfm.reference = fname moved = apply(xfm, fname, order=0) - assert np.all(imgdata == np.asanyarray(moved.dataobj, dtype=get_obj_dtype(moved.dataobj))) + assert np.all( + imgdata == np.asanyarray(moved.dataobj, dtype=get_obj_dtype(moved.dataobj)) + ) # Test ndim returned by affine assert nitl.Affine().ndim == 3 @@ -168,7 +171,10 @@ def test_concatenation(testdata_path): def test_SurfaceMesh(testdata_path): surf_path = testdata_path / "sub-200148_hemi-R_pial.surf.gii" - shape_path = testdata_path / "sub-sid000005_ses-budapest_acq-MPRAGE_hemi-R_thickness.shape.gii" + shape_path = ( + testdata_path + / "sub-sid000005_ses-budapest_acq-MPRAGE_hemi-R_thickness.shape.gii" + ) img_path = testdata_path / "bold.nii.gz" mesh = SurfaceMesh(nb.load(surf_path)) diff --git a/nitransforms/tests/test_linear.py b/nitransforms/tests/test_linear.py index 5746d5f7..969b33ab 100644 --- a/nitransforms/tests/test_linear.py +++ b/nitransforms/tests/test_linear.py @@ -1,12 +1,11 @@ # emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*- # vi: set ft=python sts=4 ts=4 sw=4 et: """Tests of linear transforms.""" -import os + import pytest import numpy as np import h5py -import nibabel as nb from nibabel.eulerangles import euler2mat from nibabel.affines import from_matvec from nitransforms import linear as nitl @@ -14,7 +13,14 @@ from .utils import assert_affines_by_filename -@pytest.mark.parametrize("matrix", [[0.0], np.ones((3, 3, 3)), np.ones((3, 4)), ]) +@pytest.mark.parametrize( + "matrix", + [ + [0.0], + np.ones((3, 3, 3)), + np.ones((3, 4)), + ], +) def test_linear_typeerrors1(matrix): """Exercise errors in Affine creation.""" with pytest.raises(TypeError): @@ -136,7 +142,9 @@ def test_loadsave(tmp_path, data_path, testdata_path, autofmt, fmt): assert np.allclose( xfm.matrix, - nitl.load(fname, fmt=supplied_fmt, reference=ref_file, moving=ref_file).matrix, + nitl.load( + fname, fmt=supplied_fmt, reference=ref_file, moving=ref_file + ).matrix, ) else: assert xfm == nitl.load(fname, fmt=supplied_fmt, reference=ref_file) @@ -146,7 +154,9 @@ def test_loadsave(tmp_path, data_path, testdata_path, autofmt, fmt): if fmt == "fsl": assert np.allclose( xfm.matrix, - nitl.load(fname, fmt=supplied_fmt, reference=ref_file, moving=ref_file).matrix, + nitl.load( + fname, fmt=supplied_fmt, reference=ref_file, moving=ref_file + ).matrix, rtol=1e-2, # FSL incurs into large errors due to rounding ) else: @@ -160,7 +170,9 @@ def test_loadsave(tmp_path, data_path, testdata_path, autofmt, fmt): if fmt == "fsl": assert np.allclose( xfm.matrix, - nitl.load(fname, fmt=supplied_fmt, reference=ref_file, moving=ref_file).matrix, + nitl.load( + fname, fmt=supplied_fmt, reference=ref_file, moving=ref_file + ).matrix, rtol=1e-2, # FSL incurs into large errors due to rounding ) else: @@ -170,7 +182,9 @@ def test_loadsave(tmp_path, data_path, testdata_path, autofmt, fmt): if fmt == "fsl": assert np.allclose( xfm.matrix, - nitl.load(fname, fmt=supplied_fmt, reference=ref_file, moving=ref_file).matrix, + nitl.load( + fname, fmt=supplied_fmt, reference=ref_file, moving=ref_file + ).matrix, rtol=1e-2, # FSL incurs into large errors due to rounding ) else: @@ -190,12 +204,15 @@ def test_linear_save(tmpdir, data_path, get_testdata, image_orientation, sw_tool T = np.linalg.inv(T) xfm = ( - nitl.Affine(T) if (sw_tool, image_orientation) != ("afni", "oblique") else + nitl.Affine(T) + if (sw_tool, image_orientation) != ("afni", "oblique") # AFNI is special when moving or reference are oblique - let io do the magic - nitl.Affine(io.afni.AFNILinearTransform.from_ras(T).to_ras( - reference=img, - moving=img, - )) + else nitl.Affine( + io.afni.AFNILinearTransform.from_ras(T).to_ras( + reference=img, + moving=img, + ) + ) ) xfm.reference = img diff --git a/nitransforms/tests/test_manip.py b/nitransforms/tests/test_manip.py index 2a2d6ffb..b5dd5c62 100644 --- a/nitransforms/tests/test_manip.py +++ b/nitransforms/tests/test_manip.py @@ -1,6 +1,7 @@ # emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*- # vi: set ft=python sts=4 ts=4 sw=4 et: """Tests of nonlinear transforms.""" + import pytest import numpy as np diff --git a/nitransforms/tests/test_nonlinear.py b/nitransforms/tests/test_nonlinear.py index 43b4584f..6112f633 100644 --- a/nitransforms/tests/test_nonlinear.py +++ b/nitransforms/tests/test_nonlinear.py @@ -1,9 +1,8 @@ # emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*- # vi: set ft=python sts=4 ts=4 sw=4 et: """Tests of nonlinear transforms.""" + import os -import shutil -from subprocess import check_call import pytest import numpy as np @@ -14,7 +13,6 @@ from nitransforms.nonlinear import ( BSplineFieldTransform, DenseFieldTransform, - load as nlload, ) from ..io.itk import ITKDisplacementsField diff --git a/nitransforms/tests/test_resampling.py b/nitransforms/tests/test_resampling.py index 3dd9aff4..2384ad97 100644 --- a/nitransforms/tests/test_resampling.py +++ b/nitransforms/tests/test_resampling.py @@ -1,6 +1,7 @@ # emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*- # vi: set ft=python sts=4 ts=4 sw=4 et: """Exercise the standalone ``apply()`` implementation.""" + import os import pytest import numpy as np @@ -50,9 +51,19 @@ } -@pytest.mark.parametrize("image_orientation", ["RAS", "LAS", "LPS", 'oblique', ]) +@pytest.mark.parametrize( + "image_orientation", + [ + "RAS", + "LAS", + "LPS", + "oblique", + ], +) @pytest.mark.parametrize("sw_tool", ["itk", "fsl", "afni", "fs"]) -def test_apply_linear_transform(tmpdir, get_testdata, get_testmask, image_orientation, sw_tool): +def test_apply_linear_transform( + tmpdir, get_testdata, get_testmask, image_orientation, sw_tool +): """Check implementation of exporting affines to formats.""" tmpdir.chdir() @@ -107,7 +118,7 @@ def test_apply_linear_transform(tmpdir, get_testdata, get_testmask, image_orient nt_moved_mask.to_filename("ntmask.nii.gz") diff = np.asanyarray(sw_moved_mask.dataobj) - np.asanyarray(nt_moved_mask.dataobj) - assert np.sqrt((diff ** 2).mean()) < RMSE_TOL_LINEAR + assert np.sqrt((diff**2).mean()) < RMSE_TOL_LINEAR brainmask = np.asanyarray(nt_moved_mask.dataobj, dtype=bool) cmd = APPLY_LINEAR_CMD[sw_tool]( @@ -123,19 +134,17 @@ def test_apply_linear_transform(tmpdir, get_testdata, get_testmask, image_orient sw_moved.set_data_dtype(img.get_data_dtype()) nt_moved = apply(xfm, img, order=0) - diff = ( - np.asanyarray(sw_moved.dataobj, dtype=sw_moved.get_data_dtype()) - - np.asanyarray(nt_moved.dataobj, dtype=nt_moved.get_data_dtype()) - ) + diff = np.asanyarray( + sw_moved.dataobj, dtype=sw_moved.get_data_dtype() + ) - np.asanyarray(nt_moved.dataobj, dtype=nt_moved.get_data_dtype()) # A certain tolerance is necessary because of resampling at borders assert np.sqrt((diff[brainmask] ** 2).mean()) < RMSE_TOL_LINEAR nt_moved = apply(xfm, "img.nii.gz", order=0) - diff = ( - np.asanyarray(sw_moved.dataobj, dtype=sw_moved.get_data_dtype()) - - np.asanyarray(nt_moved.dataobj, dtype=nt_moved.get_data_dtype()) - ) + diff = np.asanyarray( + sw_moved.dataobj, dtype=sw_moved.get_data_dtype() + ) - np.asanyarray(nt_moved.dataobj, dtype=nt_moved.get_data_dtype()) # A certain tolerance is necessary because of resampling at borders assert np.sqrt((diff[brainmask] ** 2).mean()) < RMSE_TOL_LINEAR @@ -281,7 +290,8 @@ def test_apply_transformchain(tmp_path, testdata_path): ref_fname = tmp_path / "reference.nii.gz" nb.Nifti1Image( - np.zeros(xfm.reference.shape, dtype="uint16"), xfm.reference.affine, + np.zeros(xfm.reference.shape, dtype="uint16"), + xfm.reference.affine, ).to_filename(str(ref_fname)) # Then apply the transform and cross-check with software @@ -310,7 +320,9 @@ def test_apply_transformchain(tmp_path, testdata_path): @pytest.mark.parametrize("serialize_4d", [True, False]) -def test_LinearTransformsMapping_apply(tmp_path, data_path, testdata_path, serialize_4d): +def test_LinearTransformsMapping_apply( + tmp_path, data_path, testdata_path, serialize_4d +): """Apply transform mappings.""" hmc = nitl.load( data_path / "hmc-itk.tfm", fmt="itk", reference=testdata_path / "sbref.nii.gz" @@ -333,7 +345,8 @@ def test_LinearTransformsMapping_apply(tmp_path, data_path, testdata_path, seria ) nii = apply( - hmcinv, testdata_path / "fmap.nii.gz", + hmcinv, + testdata_path / "fmap.nii.gz", order=1, serialize_nvols=2 if serialize_4d else np.inf, ) From 4c06174544f4410661d59ab7b6af3e2e689916b9 Mon Sep 17 00:00:00 2001 From: Oscar Esteban Date: Thu, 1 Aug 2024 08:49:56 +0200 Subject: [PATCH 14/20] enh: expand test coverage --- nitransforms/base.py | 6 ++++++ nitransforms/tests/test_base.py | 15 +++++++++++++++ 2 files changed, 21 insertions(+) diff --git a/nitransforms/base.py b/nitransforms/base.py index 67acc073..fa05f1f6 100644 --- a/nitransforms/base.py +++ b/nitransforms/base.py @@ -287,6 +287,12 @@ def __len__(self): By default, all transforms are of length one. This must be overriden by transforms arrays and chains. + Example + ------- + >>> T1 = TransformBase() + >>> len(T1) + 1 + """ return 1 diff --git a/nitransforms/tests/test_base.py b/nitransforms/tests/test_base.py index 4bb147fd..49d7f7af 100644 --- a/nitransforms/tests/test_base.py +++ b/nitransforms/tests/test_base.py @@ -198,3 +198,18 @@ def test_SurfaceMesh(testdata_path): with pytest.raises(TypeError): SurfaceMesh(nb.load(shape_path)) + + +def test_apply_deprecation(monkeypatch): + """Make sure a deprecation warning is issued.""" + from nitransforms import resampling + + def _retval(*args, **kwargs): + return 1 + + monkeypatch.setattr(resampling, "apply", _retval) + + with pytest.deprecated_call(): + retval = TransformBase().apply() + + assert retval == 1 From 754785f18ea57275e21c197529068e72852d7647 Mon Sep 17 00:00:00 2001 From: Oscar Esteban Date: Wed, 31 Jul 2024 11:38:07 +0200 Subject: [PATCH 15/20] enh: prepare code for easy parallelization with a process pool executor Resolves: #214. --- nitransforms/resampling.py | 41 +++++++++++++++++++++----------------- 1 file changed, 23 insertions(+), 18 deletions(-) diff --git a/nitransforms/resampling.py b/nitransforms/resampling.py index abfe2b71..bb1bb309 100644 --- a/nitransforms/resampling.py +++ b/nitransforms/resampling.py @@ -8,6 +8,7 @@ ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ## """Resampling utilities.""" +from functools import partial from pathlib import Path import numpy as np from nibabel.loadsave import load as _nbload @@ -135,33 +136,37 @@ def apply( else None ) - # Order F ensures individual volumes are contiguous in memory - # Also matches NIfTI, making final save more efficient - resampled = np.zeros( - (len(ref_ndcoords), len(transform)), dtype=input_dtype, order="F" + map_coordinates = partial( + ndi.map_coordinates, + order=order, + mode=mode, + cval=cval, + prefilter=prefilter, ) - for t in range(n_resamplings): - xfm_t = transform if n_resamplings == 1 else transform[t] + def _apply_volume(index, data, transform, targets=None): + xfm_t = transform if n_resamplings == 1 else transform[index] if targets is None: targets = ImageGrid(spatialimage).index( # data should be an image _as_homogeneous(xfm_t.map(ref_ndcoords), dim=_ref.ndim) ) - # Interpolate - resampled[..., t] = ndi.map_coordinates( - ( - data - if data is not None - else spatialimage.dataobj[..., t].astype(input_dtype, copy=False) - ), - targets, - order=order, - mode=mode, - cval=cval, - prefilter=prefilter, + data_t = ( + data + if data is not None + else spatialimage.dataobj[..., index].astype(input_dtype, copy=False) ) + return map_coordinates(data_t, targets) + + # Order F ensures individual volumes are contiguous in memory + # Also matches NIfTI, making final save more efficient + resampled = np.zeros( + (len(ref_ndcoords), len(transform)), dtype=input_dtype, order="F" + ) + for t in range(n_resamplings): + # Interpolate + resampled[..., t] = _apply_volume(t, data, transform, targets=targets) else: data = np.asanyarray(spatialimage.dataobj, dtype=input_dtype) From 38bb388374fcb900cde1ff966e58cad66658ff0d Mon Sep 17 00:00:00 2001 From: Oscar Esteban Date: Wed, 31 Jul 2024 11:44:17 +0200 Subject: [PATCH 16/20] enh: create process pool --- nitransforms/resampling.py | 83 +++++++++++++++++++++++++------------- 1 file changed, 56 insertions(+), 27 deletions(-) diff --git a/nitransforms/resampling.py b/nitransforms/resampling.py index bb1bb309..1b76dba1 100644 --- a/nitransforms/resampling.py +++ b/nitransforms/resampling.py @@ -8,7 +8,8 @@ ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ## """Resampling utilities.""" -from functools import partial +from os import cpu_count +from concurrent.futures import ProcessPoolExecutor, as_completed from pathlib import Path import numpy as np from nibabel.loadsave import load as _nbload @@ -26,6 +27,25 @@ """Minimum number of volumes to automatically serialize 4D transforms.""" +def _apply_volume( + index, + data, + targets, + order=3, + mode="constant", + cval=0.0, + prefilter=True, +): + return index, ndi.map_coordinates( + data, + targets, + order=order, + mode=mode, + cval=cval, + prefilter=prefilter, + ) + + def apply( transform, spatialimage, @@ -136,38 +156,47 @@ def apply( else None ) - map_coordinates = partial( - ndi.map_coordinates, - order=order, - mode=mode, - cval=cval, - prefilter=prefilter, - ) + if njobs is None: + njobs = cpu_count() + + with ProcessPoolExecutor(max_workers=min(njobs, n_resamplings)) as executor: + results = [] + for t in range(n_resamplings): + xfm_t = transform if n_resamplings == 1 else transform[t] - def _apply_volume(index, data, transform, targets=None): - xfm_t = transform if n_resamplings == 1 else transform[index] + if targets is None: + targets = ImageGrid(spatialimage).index( # data should be an image + _as_homogeneous(xfm_t.map(ref_ndcoords), dim=_ref.ndim) + ) - if targets is None: - targets = ImageGrid(spatialimage).index( # data should be an image - _as_homogeneous(xfm_t.map(ref_ndcoords), dim=_ref.ndim) + data_t = ( + data + if data is not None + else spatialimage.dataobj[..., t].astype(input_dtype, copy=False) ) - data_t = ( - data - if data is not None - else spatialimage.dataobj[..., index].astype(input_dtype, copy=False) - ) - return map_coordinates(data_t, targets) + results.append( + executor.submit( + _apply_volume, + t, + data_t, + targets, + order=order, + mode=mode, + cval=cval, + prefilter=prefilter, + ) + ) - # Order F ensures individual volumes are contiguous in memory - # Also matches NIfTI, making final save more efficient - resampled = np.zeros( - (len(ref_ndcoords), len(transform)), dtype=input_dtype, order="F" - ) - for t in range(n_resamplings): - # Interpolate - resampled[..., t] = _apply_volume(t, data, transform, targets=targets) + # Order F ensures individual volumes are contiguous in memory + # Also matches NIfTI, making final save more efficient + resampled = np.zeros( + (len(ref_ndcoords), len(transform)), dtype=input_dtype, order="F" + ) + for future in as_completed(results): + t, resampled_t = future.result() + resampled[..., t] = resampled_t else: data = np.asanyarray(spatialimage.dataobj, dtype=input_dtype) From 026a10af983eff2ca01752918544d667f8e35ddf Mon Sep 17 00:00:00 2001 From: Oscar Esteban Date: Thu, 1 Aug 2024 09:47:48 +0200 Subject: [PATCH 17/20] enh: expand test coverage --- nitransforms/resampling.py | 3 +-- nitransforms/tests/test_resampling.py | 15 ++++++++++++++- 2 files changed, 15 insertions(+), 3 deletions(-) diff --git a/nitransforms/resampling.py b/nitransforms/resampling.py index 1b76dba1..f20b9e2b 100644 --- a/nitransforms/resampling.py +++ b/nitransforms/resampling.py @@ -156,8 +156,7 @@ def apply( else None ) - if njobs is None: - njobs = cpu_count() + njobs = cpu_count() if njobs is None or njobs < 1 else njobs with ProcessPoolExecutor(max_workers=min(njobs, n_resamplings)) as executor: results = [] diff --git a/nitransforms/tests/test_resampling.py b/nitransforms/tests/test_resampling.py index 2384ad97..f944b225 100644 --- a/nitransforms/tests/test_resampling.py +++ b/nitransforms/tests/test_resampling.py @@ -15,7 +15,7 @@ from nitransforms import nonlinear as nitnl from nitransforms import manip as nitm from nitransforms import io -from nitransforms.resampling import apply +from nitransforms.resampling import apply, _apply_volume RMSE_TOL_LINEAR = 0.09 RMSE_TOL_NONLINEAR = 0.05 @@ -363,3 +363,16 @@ def test_LinearTransformsMapping_apply( reference=testdata_path / "sbref.nii.gz", serialize_nvols=2 if serialize_4d else np.inf, ) + + +@pytest.mark.parametrize("t", list(range(4))) +def test_apply_helper(monkeypatch, t): + """Ensure the apply helper function correctly just decorates with index.""" + from nitransforms.resampling import ndi + + def _retval(*args, **kwargs): + return 1 + + monkeypatch.setattr(ndi, "map_coordinates", _retval) + + assert _apply_volume(t, None, None) == (t, 1) From 7dcc78daa0a172f276c55f38e935ef3c2df87089 Mon Sep 17 00:00:00 2001 From: Oscar Esteban Date: Thu, 1 Aug 2024 10:21:11 +0200 Subject: [PATCH 18/20] sty: add type annotations --- nitransforms/resampling.py | 96 +++++++++++++++++++++++++++----------- 1 file changed, 68 insertions(+), 28 deletions(-) diff --git a/nitransforms/resampling.py b/nitransforms/resampling.py index f20b9e2b..6cbbc1e9 100644 --- a/nitransforms/resampling.py +++ b/nitransforms/resampling.py @@ -11,13 +11,17 @@ from os import cpu_count from concurrent.futures import ProcessPoolExecutor, as_completed from pathlib import Path +from typing import Tuple + import numpy as np from nibabel.loadsave import load as _nbload from nibabel.arrayproxy import get_obj_dtype +from nibabel.spatialimages import SpatialImage from scipy import ndimage as ndi from nitransforms.base import ( ImageGrid, + TransformBase, TransformError, SpatialReference, _as_homogeneous, @@ -28,14 +32,49 @@ def _apply_volume( - index, - data, - targets, - order=3, - mode="constant", - cval=0.0, - prefilter=True, -): + index: int, + data: np.ndarray, + targets: np.ndarray, + order: int = 3, + mode: str = "constant", + cval: float = 0.0, + prefilter: bool = True, +) -> Tuple[int, np.ndarray]: + """ + Decorate :obj:`~scipy.ndimage.map_coordinates` to return an order index for parallelization. + + Parameters + ---------- + index : :obj:`int` + The index of the volume to apply the interpolation to. + data : :obj:`~numpy.ndarray` + The input data array. + targets : :obj:`~numpy.ndarray` + The target coordinates for mapping. + order : :obj:`int`, optional + The order of the spline interpolation, default is 3. + The order has to be in the range 0-5. + mode : :obj:`str`, optional + Determines how the input image is extended when the resamplings overflows + a border. One of ``'constant'``, ``'reflect'``, ``'nearest'``, ``'mirror'``, + or ``'wrap'``. Default is ``'constant'``. + cval : :obj:`float`, optional + Constant value for ``mode='constant'``. Default is 0.0. + prefilter: :obj:`bool`, optional + Determines if the image's data array is prefiltered with + a spline filter before interpolation. The default is ``True``, + which will create a temporary *float64* array of filtered values + if *order > 1*. If setting this to ``False``, the output will be + slightly blurred if *order > 1*, unless the input is prefiltered, + i.e. it is the result of calling the spline filter on the original + input. + + Returns + ------- + (:obj:`int`, :obj:`~numpy.ndarray`) + The index and the array resulting from the interpolation. + + """ return index, ndi.map_coordinates( data, targets, @@ -47,37 +86,38 @@ def _apply_volume( def apply( - transform, - spatialimage, - reference=None, - order=3, - mode="constant", - cval=0.0, - prefilter=True, - output_dtype=None, - serialize_nvols=SERIALIZE_VOLUME_WINDOW_WIDTH, - njobs=None, -): + transform: TransformBase, + spatialimage: str | Path | SpatialImage, + reference: str | Path | SpatialImage = None, + order: int = 3, + mode: str = "constant", + cval: float = 0.0, + prefilter: bool = True, + output_dtype: np.dtype = None, + serialize_nvols: int = SERIALIZE_VOLUME_WINDOW_WIDTH, + njobs: int = None, +) -> SpatialImage | np.ndarray: """ Apply a transformation to an image, resampling on the reference spatial object. Parameters ---------- - spatialimage : `spatialimage` + spatialimage : :obj:`~nibabel.spatialimages.SpatialImage` or `os.pathlike` The image object containing the data to be resampled in reference space - reference : spatial object, optional + reference : :obj:`~nibabel.spatialimages.SpatialImage` or `os.pathlike` The image, surface, or combination thereof containing the coordinates of samples that will be sampled. - order : int, optional + order : :obj:`int`, optional The order of the spline interpolation, default is 3. The order has to be in the range 0-5. - mode : {'constant', 'reflect', 'nearest', 'mirror', 'wrap'}, optional + mode : :obj:`str`, optional Determines how the input image is extended when the resamplings overflows - a border. Default is 'constant'. - cval : float, optional + a border. One of ``'constant'``, ``'reflect'``, ``'nearest'``, ``'mirror'``, + or ``'wrap'``. Default is ``'constant'``. + cval : :obj:`float`, optional Constant value for ``mode='constant'``. Default is 0.0. - prefilter: bool, optional + prefilter: :obj:`bool`, optional Determines if the image's data array is prefiltered with a spline filter before interpolation. The default is ``True``, which will create a temporary *float64* array of filtered values @@ -85,7 +125,7 @@ def apply( slightly blurred if *order > 1*, unless the input is prefiltered, i.e. it is the result of calling the spline filter on the original input. - output_dtype: dtype specifier, optional + output_dtype: :obj:`~numpy.dtype`, optional The dtype of the returned array or image, if specified. If ``None``, the default behavior is to use the effective dtype of the input image. If slope and/or intercept are defined, the effective @@ -97,7 +137,7 @@ def apply( Returns ------- - resampled : `spatialimage` or ndarray + resampled : :obj:`~nibabel.spatialimages.SpatialImage` or :obj:`~numpy.ndarray` The data imaged after resampling to reference space. """ From 79305a9f79010e4a46687e5d8527dde871b7d79d Mon Sep 17 00:00:00 2001 From: Oscar Esteban Date: Fri, 2 Aug 2024 09:49:27 +0200 Subject: [PATCH 19/20] enh: implement a memory limitation mechanism in loading data Resolves: #218. Co-authored-by: Chris Markiewicz --- nitransforms/resampling.py | 39 +++++++++++++++++++++++++++++++++++++- 1 file changed, 38 insertions(+), 1 deletion(-) diff --git a/nitransforms/resampling.py b/nitransforms/resampling.py index 6cbbc1e9..430abf1a 100644 --- a/nitransforms/resampling.py +++ b/nitransforms/resampling.py @@ -96,6 +96,7 @@ def apply( output_dtype: np.dtype = None, serialize_nvols: int = SERIALIZE_VOLUME_WINDOW_WIDTH, njobs: int = None, + dtype_width: int = 8, ) -> SpatialImage | np.ndarray: """ Apply a transformation to an image, resampling on the reference spatial object. @@ -134,6 +135,10 @@ def apply( If ``reference`` is defined, then the return value is an image, with a data array of the effective dtype but with the on-disk dtype set to the input image's on-disk dtype. + dtype_width: :obj:`int` + Cap the width of the input data type to the given number of bytes. + This argument is intended to work as a way to implement lower memory + requirements in resampling. Returns ------- @@ -157,7 +162,7 @@ def apply( spatialimage = _nbload(str(spatialimage)) # Avoid opening the data array just yet - input_dtype = get_obj_dtype(spatialimage.dataobj) + input_dtype = cap_dtype(get_obj_dtype(spatialimage.dataobj), dtype_width) # Number of data volumes data_nvols = 1 if spatialimage.ndim < 4 else spatialimage.shape[-1] @@ -277,3 +282,35 @@ def apply( output_dtype = output_dtype or input_dtype return resampled.astype(output_dtype) + + +def cap_dtype(dt, nbytes): + """ + Cap the datatype size to shave off memory requirements. + + Examples + -------- + >>> cap_dtype(np.dtype('f8'), 4) + dtype('float32') + + >>> cap_dtype(np.dtype('f8'), 16) + dtype('float64') + + >>> cap_dtype('float64', 4) + dtype('float32') + + >>> cap_dtype(np.dtype('i1'), 4) + dtype('int8') + + >>> cap_dtype('int8', 4) + dtype('int8') + + >>> cap_dtype('int32', 1) + dtype('int8') + + >>> cap_dtype(np.dtype('i8'), 4) + dtype('int32') + + """ + dt = np.dtype(dt) + return np.dtype(f"{dt.byteorder}{dt.kind}{min(nbytes, dt.itemsize)}") From 063e1f0c6dd2dd453afb35b76554bebdddd73dd9 Mon Sep 17 00:00:00 2001 From: Oscar Esteban Date: Tue, 6 Aug 2024 08:57:11 +0200 Subject: [PATCH 20/20] enh: port from process pool into asyncio concurrent Co-authored-by: Chris Markiewicz --- nitransforms/resampling.py | 172 +++++++++++++++++--------- nitransforms/tests/test_resampling.py | 15 +-- 2 files changed, 112 insertions(+), 75 deletions(-) diff --git a/nitransforms/resampling.py b/nitransforms/resampling.py index 430abf1a..d7c7f9c5 100644 --- a/nitransforms/resampling.py +++ b/nitransforms/resampling.py @@ -8,10 +8,11 @@ ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ## """Resampling utilities.""" +import asyncio from os import cpu_count -from concurrent.futures import ProcessPoolExecutor, as_completed +from functools import partial from pathlib import Path -from typing import Tuple +from typing import Callable, TypeVar import numpy as np from nibabel.loadsave import load as _nbload @@ -27,30 +28,58 @@ _as_homogeneous, ) +R = TypeVar("R") + SERIALIZE_VOLUME_WINDOW_WIDTH: int = 8 """Minimum number of volumes to automatically serialize 4D transforms.""" -def _apply_volume( - index: int, +async def worker(job: Callable[[], R], semaphore) -> R: + async with semaphore: + loop = asyncio.get_running_loop() + return await loop.run_in_executor(None, job) + + +async def _apply_serial( data: np.ndarray, + spatialimage: SpatialImage, targets: np.ndarray, + transform: TransformBase, + ref_ndim: int, + ref_ndcoords: np.ndarray, + n_resamplings: int, + output: np.ndarray, + input_dtype: np.dtype, order: int = 3, mode: str = "constant", cval: float = 0.0, prefilter: bool = True, -) -> Tuple[int, np.ndarray]: + max_concurrent: int = min(cpu_count(), 12), +): """ - Decorate :obj:`~scipy.ndimage.map_coordinates` to return an order index for parallelization. + Resample through a given transform serially, in a 3D+t setting. Parameters ---------- - index : :obj:`int` - The index of the volume to apply the interpolation to. data : :obj:`~numpy.ndarray` The input data array. + spatialimage : :obj:`~nibabel.spatialimages.SpatialImage` or `os.pathlike` + The image object containing the data to be resampled in reference + space targets : :obj:`~numpy.ndarray` The target coordinates for mapping. + transform : :obj:`~nitransforms.base.TransformBase` + The 3D, 3D+t, or 4D transform through which data will be resampled. + ref_ndim : :obj:`int` + Dimensionality of the resampling target (reference image). + ref_ndcoords : :obj:`~numpy.ndarray` + Physical coordinates (RAS+) where data will be interpolated, if the resampling + target is a grid, the scanner coordinates of all voxels. + n_resamplings : :obj:`int` + Total number of 3D resamplings (can be defined by the input image, the transform, + or be matched, that is, same number of volumes in the input and number of transforms). + output : :obj:`~numpy.ndarray` + The output data array where resampled values will be stored volume-by-volume. order : :obj:`int`, optional The order of the spline interpolation, default is 3. The order has to be in the range 0-5. @@ -71,18 +100,46 @@ def _apply_volume( Returns ------- - (:obj:`int`, :obj:`~numpy.ndarray`) - The index and the array resulting from the interpolation. + np.ndarray + Data resampled on the 3D+t array of input coordinates. """ - return index, ndi.map_coordinates( - data, - targets, - order=order, - mode=mode, - cval=cval, - prefilter=prefilter, - ) + tasks = [] + semaphore = asyncio.Semaphore(max_concurrent) + + for t in range(n_resamplings): + xfm_t = transform if n_resamplings == 1 else transform[t] + + if targets is None: + targets = ImageGrid(spatialimage).index( # data should be an image + _as_homogeneous(xfm_t.map(ref_ndcoords), dim=ref_ndim) + ) + + data_t = ( + data + if data is not None + else spatialimage.dataobj[..., t].astype(input_dtype, copy=False) + ) + + tasks.append( + asyncio.create_task( + worker( + partial( + ndi.map_coordinates, + data_t, + targets, + output=output[..., t], + order=order, + mode=mode, + cval=cval, + prefilter=prefilter, + ), + semaphore, + ) + ) + ) + await asyncio.gather(*tasks) + return output def apply( @@ -94,15 +151,17 @@ def apply( cval: float = 0.0, prefilter: bool = True, output_dtype: np.dtype = None, - serialize_nvols: int = SERIALIZE_VOLUME_WINDOW_WIDTH, - njobs: int = None, dtype_width: int = 8, + serialize_nvols: int = SERIALIZE_VOLUME_WINDOW_WIDTH, + max_concurrent: int = min(cpu_count(), 12), ) -> SpatialImage | np.ndarray: """ Apply a transformation to an image, resampling on the reference spatial object. Parameters ---------- + transform: :obj:`~nitransforms.base.TransformBase` + The 3D, 3D+t, or 4D transform through which data will be resampled. spatialimage : :obj:`~nibabel.spatialimages.SpatialImage` or `os.pathlike` The image object containing the data to be resampled in reference space @@ -118,7 +177,7 @@ def apply( or ``'wrap'``. Default is ``'constant'``. cval : :obj:`float`, optional Constant value for ``mode='constant'``. Default is 0.0. - prefilter: :obj:`bool`, optional + prefilter : :obj:`bool`, optional Determines if the image's data array is prefiltered with a spline filter before interpolation. The default is ``True``, which will create a temporary *float64* array of filtered values @@ -126,7 +185,7 @@ def apply( slightly blurred if *order > 1*, unless the input is prefiltered, i.e. it is the result of calling the spline filter on the original input. - output_dtype: :obj:`~numpy.dtype`, optional + output_dtype : :obj:`~numpy.dtype`, optional The dtype of the returned array or image, if specified. If ``None``, the default behavior is to use the effective dtype of the input image. If slope and/or intercept are defined, the effective @@ -135,10 +194,17 @@ def apply( If ``reference`` is defined, then the return value is an image, with a data array of the effective dtype but with the on-disk dtype set to the input image's on-disk dtype. - dtype_width: :obj:`int` + dtype_width : :obj:`int` Cap the width of the input data type to the given number of bytes. This argument is intended to work as a way to implement lower memory requirements in resampling. + serialize_nvols : :obj:`int` + Minimum number of volumes in a 3D+t (that is, a series of 3D transformations + independent in time) to resample on a one-by-one basis. + Serialized resampling can be executed concurrently (parallelized) with + the argument ``max_concurrent``. + max_concurrent : :obj:`int` + Maximum number of 3D resamplings to be executed concurrently. Returns ------- @@ -201,46 +267,30 @@ def apply( else None ) - njobs = cpu_count() if njobs is None or njobs < 1 else njobs - - with ProcessPoolExecutor(max_workers=min(njobs, n_resamplings)) as executor: - results = [] - for t in range(n_resamplings): - xfm_t = transform if n_resamplings == 1 else transform[t] - - if targets is None: - targets = ImageGrid(spatialimage).index( # data should be an image - _as_homogeneous(xfm_t.map(ref_ndcoords), dim=_ref.ndim) - ) - - data_t = ( - data - if data is not None - else spatialimage.dataobj[..., t].astype(input_dtype, copy=False) - ) - - results.append( - executor.submit( - _apply_volume, - t, - data_t, - targets, - order=order, - mode=mode, - cval=cval, - prefilter=prefilter, - ) - ) + # Order F ensures individual volumes are contiguous in memory + # Also matches NIfTI, making final save more efficient + resampled = np.zeros( + (len(ref_ndcoords), len(transform)), dtype=input_dtype, order="F" + ) - # Order F ensures individual volumes are contiguous in memory - # Also matches NIfTI, making final save more efficient - resampled = np.zeros( - (len(ref_ndcoords), len(transform)), dtype=input_dtype, order="F" + resampled = asyncio.run( + _apply_serial( + data, + spatialimage, + targets, + transform, + _ref.ndim, + ref_ndcoords, + n_resamplings, + resampled, + input_dtype, + order=order, + mode=mode, + cval=cval, + prefilter=prefilter, + max_concurrent=max_concurrent, ) - - for future in as_completed(results): - t, resampled_t = future.result() - resampled[..., t] = resampled_t + ) else: data = np.asanyarray(spatialimage.dataobj, dtype=input_dtype) diff --git a/nitransforms/tests/test_resampling.py b/nitransforms/tests/test_resampling.py index f944b225..2384ad97 100644 --- a/nitransforms/tests/test_resampling.py +++ b/nitransforms/tests/test_resampling.py @@ -15,7 +15,7 @@ from nitransforms import nonlinear as nitnl from nitransforms import manip as nitm from nitransforms import io -from nitransforms.resampling import apply, _apply_volume +from nitransforms.resampling import apply RMSE_TOL_LINEAR = 0.09 RMSE_TOL_NONLINEAR = 0.05 @@ -363,16 +363,3 @@ def test_LinearTransformsMapping_apply( reference=testdata_path / "sbref.nii.gz", serialize_nvols=2 if serialize_4d else np.inf, ) - - -@pytest.mark.parametrize("t", list(range(4))) -def test_apply_helper(monkeypatch, t): - """Ensure the apply helper function correctly just decorates with index.""" - from nitransforms.resampling import ndi - - def _retval(*args, **kwargs): - return 1 - - monkeypatch.setattr(ndi, "map_coordinates", _retval) - - assert _apply_volume(t, None, None) == (t, 1)