From d865ec56a5905fd0b3cd1089ce345fe30f1dbcae Mon Sep 17 00:00:00 2001 From: Julien Marabotto <166002186+jmarabotto@users.noreply.github.com> Date: Fri, 5 Apr 2024 10:00:33 +0200 Subject: [PATCH 01/14] enh: outsource the apply function --- nitransforms/base.py | 95 ------------------------------- nitransforms/resampling.py | 114 +++++++++++++++++++++++++++++++++++++ 2 files changed, 114 insertions(+), 95 deletions(-) create mode 100644 nitransforms/resampling.py diff --git a/nitransforms/base.py b/nitransforms/base.py index 68b97f75..3f22990f 100644 --- a/nitransforms/base.py +++ b/nitransforms/base.py @@ -222,101 +222,6 @@ def ndim(self): """Access the dimensions of the reference space.""" return self.reference.ndim - def apply( - self, - spatialimage, - reference=None, - order=3, - mode="constant", - cval=0.0, - prefilter=True, - output_dtype=None, - ): - """ - Apply a transformation to an image, resampling on the reference spatial object. - - Parameters - ---------- - spatialimage : `spatialimage` - The image object containing the data to be resampled in reference - space - reference : spatial object, optional - The image, surface, or combination thereof containing the coordinates - of samples that will be sampled. - order : 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 - Determines how the input image is extended when the resamplings overflows - a border. Default is 'constant'. - cval : float, optional - Constant value for ``mode='constant'``. Default is 0.0. - prefilter: 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. - output_dtype: dtype specifier, 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 - dtype is float64, otherwise it is equivalent to the input image's - ``get_data_dtype()`` (on-disk type). - 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. - - Returns - ------- - resampled : `spatialimage` or ndarray - The data imaged after resampling to reference space. - - """ - if reference is not None and isinstance(reference, (str, Path)): - reference = _nbload(str(reference)) - - _ref = ( - self.reference if reference is None else SpatialReference.factory(reference) - ) - - if _ref is None: - raise TransformError("Cannot apply transform without reference") - - if isinstance(spatialimage, (str, Path)): - spatialimage = _nbload(str(spatialimage)) - - data = np.asanyarray(spatialimage.dataobj) - targets = ImageGrid(spatialimage).index( # data should be an image - _as_homogeneous(self.map(_ref.ndcoords.T), dim=_ref.ndim) - ) - - resampled = ndi.map_coordinates( - data, - targets.T, - output=output_dtype, - order=order, - mode=mode, - cval=cval, - prefilter=prefilter, - ) - - 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()) - moved = spatialimage.__class__( - resampled.reshape(_ref.shape), - _ref.affine, - hdr, - ) - return moved - - return resampled - def map(self, x, inverse=False): r""" Apply :math:`y = f(x)`. diff --git a/nitransforms/resampling.py b/nitransforms/resampling.py new file mode 100644 index 00000000..7d2765ac --- /dev/null +++ b/nitransforms/resampling.py @@ -0,0 +1,114 @@ +# emacs: -*- mode: python-mode; py-indent-offset: 4; indent-tabs-mode: nil -*- +# vi: set ft=python sts=4 ts=4 sw=4 et: +### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ## +# +# See COPYING file distributed along with the NiBabel package for the +# copyright and license terms. +# +### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ## +"""Resampling utilities.""" +from pathlib import Path +import numpy as np +import h5py +import warnings +from nibabel.loadsave import load as _nbload +from nibabel import funcs as _nbfuncs +from nibabel.nifti1 import intent_codes as INTENT_CODES +from nibabel.cifti2 import Cifti2Image +from scipy import ndimage as ndi + + +def apply( + transform, + spatialimage, + reference=None, + order=3, + mode="constant", + cval=0.0, + prefilter=True, + output_dtype=None, +): + """ + Apply a transformation to an image, resampling on the reference spatial object. + + Parameters + ---------- + spatialimage : `spatialimage` + The image object containing the data to be resampled in reference + space + reference : spatial object, optional + The image, surface, or combination thereof containing the coordinates + of samples that will be sampled. + order : 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 + Determines how the input image is extended when the resamplings overflows + a border. Default is 'constant'. + cval : float, optional + Constant value for ``mode='constant'``. Default is 0.0. + prefilter: 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. + output_dtype: dtype specifier, 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 + dtype is float64, otherwise it is equivalent to the input image's + ``get_data_dtype()`` (on-disk type). + 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. + + Returns + ------- + resampled : `spatialimage` or ndarray + The data imaged after resampling to reference space. + + """ + if reference is not None and isinstance(reference, (str, Path)): + reference = _nbload(str(reference)) + + _ref = ( + transform.reference if reference is None else SpatialReference.factory(reference) + ) + + if _ref is None: + raise TransformError("Cannot apply transform without reference") + + if isinstance(spatialimage, (str, Path)): + spatialimage = _nbload(str(spatialimage)) + + data = np.asanyarray(spatialimage.dataobj) + targets = ImageGrid(spatialimage).index( # data should be an image + _as_homogeneous(transform.map(_ref.ndcoords.T), dim=_ref.ndim) + ) + + resampled = ndi.map_coordinates( + data, + targets.T, + output=output_dtype, + order=order, + mode=mode, + cval=cval, + prefilter=prefilter, + ) + + 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()) + moved = spatialimage.__class__( + resampled.reshape(_ref.shape), + _ref.affine, + hdr, + ) + return moved + + return resampled From d6dca85e0e01e0c1046dc7f15f7f45c6ac66e591 Mon Sep 17 00:00:00 2001 From: Julien Marabotto <166002186+jmarabotto@users.noreply.github.com> Date: Fri, 5 Apr 2024 10:06:05 +0200 Subject: [PATCH 02/14] sty: pacify flake8 --- nitransforms/resampling.py | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/nitransforms/resampling.py b/nitransforms/resampling.py index 7d2765ac..a876bb12 100644 --- a/nitransforms/resampling.py +++ b/nitransforms/resampling.py @@ -9,14 +9,14 @@ """Resampling utilities.""" from pathlib import Path import numpy as np -import h5py -import warnings from nibabel.loadsave import load as _nbload -from nibabel import funcs as _nbfuncs -from nibabel.nifti1 import intent_codes as INTENT_CODES -from nibabel.cifti2 import Cifti2Image -from scipy import ndimage as ndi +from nitransforms.base import ( + ImageGrid, + TransformError, + SpatialReference, + _as_homogeneous, +) def apply( transform, From e299dc182087c37bf70ae1c0044b4e1285b285dd Mon Sep 17 00:00:00 2001 From: Julien Marabotto <166002186+jmarabotto@users.noreply.github.com> Date: Fri, 5 Apr 2024 10:09:33 +0200 Subject: [PATCH 03/14] sty: fix imports --- nitransforms/base.py | 1 - nitransforms/resampling.py | 2 ++ 2 files changed, 2 insertions(+), 1 deletion(-) diff --git a/nitransforms/base.py b/nitransforms/base.py index 3f22990f..3d9d6895 100644 --- a/nitransforms/base.py +++ b/nitransforms/base.py @@ -15,7 +15,6 @@ from nibabel import funcs as _nbfuncs from nibabel.nifti1 import intent_codes as INTENT_CODES from nibabel.cifti2 import Cifti2Image -from scipy import ndimage as ndi EQUALITY_TOL = 1e-5 diff --git a/nitransforms/resampling.py b/nitransforms/resampling.py index a876bb12..e89b081a 100644 --- a/nitransforms/resampling.py +++ b/nitransforms/resampling.py @@ -10,6 +10,7 @@ from pathlib import Path import numpy as np from nibabel.loadsave import load as _nbload +from scipy import ndimage as ndi from nitransforms.base import ( ImageGrid, @@ -18,6 +19,7 @@ _as_homogeneous, ) + def apply( transform, spatialimage, From ffd8a6d404e23008725bf9e162ec4549bf12d465 Mon Sep 17 00:00:00 2001 From: Oscar Esteban Date: Fri, 5 Apr 2024 10:35:24 +0200 Subject: [PATCH 04/14] fix: update many apply() calls --- nitransforms/cli.py | 4 ++- nitransforms/tests/test_base.py | 7 +++--- nitransforms/tests/test_io.py | 41 ++++++++++++++++++++----------- nitransforms/tests/test_linear.py | 16 ++++++------ 4 files changed, 43 insertions(+), 25 deletions(-) diff --git a/nitransforms/cli.py b/nitransforms/cli.py index 63b8bed4..8f8f5ce0 100644 --- a/nitransforms/cli.py +++ b/nitransforms/cli.py @@ -5,6 +5,7 @@ from .linear import load as linload from .nonlinear import load as nlinload +from .resampling import apply def cli_apply(pargs): @@ -38,7 +39,8 @@ def cli_apply(pargs): # ensure a reference is set xfm.reference = pargs.ref or pargs.moving - moved = xfm.apply( + moved = apply( + xfm, pargs.moving, order=pargs.order, mode=pargs.mode, diff --git a/nitransforms/tests/test_base.py b/nitransforms/tests/test_base.py index 8422ca10..d09096f7 100644 --- a/nitransforms/tests/test_base.py +++ b/nitransforms/tests/test_base.py @@ -6,6 +6,7 @@ from ..base import SpatialReference, SampledSpatialData, ImageGrid, TransformBase from .. import linear as nitl +from ..resampling import apply def test_SpatialReference(testdata_path): @@ -95,7 +96,7 @@ def _to_hdf5(klass, x5_root): xfm = TransformBase() xfm.reference = fname assert xfm.ndim == 3 - moved = xfm.apply(fname, order=0) + moved = apply(xfm, fname, order=0) assert np.all( imgdata == np.asanyarray(moved.dataobj, dtype=moved.get_data_dtype()) ) @@ -104,7 +105,7 @@ def _to_hdf5(klass, x5_root): xfm = TransformBase() xfm.reference = fname assert xfm.ndim == 3 - moved = xfm.apply(str(fname), reference=fname, order=0) + moved = apply(xfm, str(fname), reference=fname, order=0) assert np.all( imgdata == np.asanyarray(moved.dataobj, dtype=moved.get_data_dtype()) ) @@ -118,7 +119,7 @@ def _to_hdf5(klass, x5_root): ) ] ) - giimoved = xfm.apply(fname, reference=gii, order=0) + giimoved = apply(xfm, fname, reference=gii, order=0) assert np.allclose(giimoved.reshape(xfm.reference.shape), moved.get_fdata()) # Test to_filename diff --git a/nitransforms/tests/test_io.py b/nitransforms/tests/test_io.py index bcee9198..0cc79d15 100644 --- a/nitransforms/tests/test_io.py +++ b/nitransforms/tests/test_io.py @@ -28,6 +28,8 @@ ) from nitransforms.io.base import LinearParameters, TransformIOError, TransformFileError from nitransforms.conftest import _datadir, _testdir +from nitransforms.resampling import apply + LPS = np.diag([-1, -1, 1, 1]) ITK_MAT = LPS.dot(np.ones((4, 4)).dot(LPS)) @@ -497,10 +499,13 @@ def test_afni_oblique(tmpdir, parameters, swapaxes, testdata_path, dir_x, dir_y, assert np.allclose(card_aff, nb.load("deob_3drefit.nii.gz").affine) # Check that nitransforms can emulate 3drefit -deoblique - nt3drefit = Affine( - afni._cardinal_rotation(img.affine, False), - reference="deob_3drefit.nii.gz", - ).apply("orig.nii.gz") + nt3drefit = apply( + Affine( + afni._cardinal_rotation(img.affine, False), + reference="deob_3drefit.nii.gz", + ), + "orig.nii.gz", + ) diff = ( np.asanyarray(img.dataobj, dtype="uint8") @@ -509,10 +514,13 @@ def test_afni_oblique(tmpdir, parameters, swapaxes, testdata_path, dir_x, dir_y, assert np.sqrt((diff[10:-10, 10:-10, 10:-10] ** 2).mean()) < 0.1 # Check that nitransforms can revert 3drefit -deoblique - nt_undo3drefit = Affine( - afni._cardinal_rotation(img.affine, True), - reference="orig.nii.gz", - ).apply("deob_3drefit.nii.gz") + nt_undo3drefit = apply( + Affine( + afni._cardinal_rotation(img.affine, True), + reference="orig.nii.gz", + ), + "deob_3drefit.nii.gz", + ) diff = ( np.asanyarray(img.dataobj, dtype="uint8") @@ -531,16 +539,21 @@ def test_afni_oblique(tmpdir, parameters, swapaxes, testdata_path, dir_x, dir_y, assert np.allclose(deobaff, deobnii.affine) # Check resampling in deobliqued grid - ntdeobnii = Affine(np.eye(4), reference=deobnii.__class__( - np.zeros(deobshape, dtype="uint8"), - deobaff, - deobnii.header - )).apply(img, order=0) + ntdeobnii = apply( + Affine(np.eye(4), reference=deobnii.__class__( + np.zeros(deobshape, dtype="uint8"), + deobaff, + deobnii.header + )), + img, + order=0, + ) # Generate an internal box to exclude border effects box = np.zeros(img.shape, dtype="uint8") box[10:-10, 10:-10, 10:-10] = 1 - ntdeobmask = Affine(np.eye(4), reference=ntdeobnii).apply( + ntdeobmask = apply( + Affine(np.eye(4), reference=ntdeobnii), nb.Nifti1Image(box, img.affine, img.header), order=0, ) diff --git a/nitransforms/tests/test_linear.py b/nitransforms/tests/test_linear.py index 2957f59c..9a06fe32 100644 --- a/nitransforms/tests/test_linear.py +++ b/nitransforms/tests/test_linear.py @@ -13,6 +13,7 @@ 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 @@ -285,7 +286,7 @@ def test_apply_linear_transform(tmpdir, get_testdata, get_testmask, image_orient assert exit_code == 0 sw_moved_mask = nb.load("resampled_brainmask.nii.gz") - nt_moved_mask = xfm.apply(msk, order=0) + 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) @@ -305,7 +306,7 @@ def test_apply_linear_transform(tmpdir, get_testdata, get_testmask, image_orient sw_moved = nb.load("resampled.nii.gz") sw_moved.set_data_dtype(img.get_data_dtype()) - nt_moved = xfm.apply(img, order=0) + 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()) @@ -314,7 +315,7 @@ def test_apply_linear_transform(tmpdir, get_testdata, get_testmask, image_orient # A certain tolerance is necessary because of resampling at borders assert np.sqrt((diff[brainmask] ** 2).mean()) < RMSE_TOL - nt_moved = xfm.apply("img.nii.gz", order=0) + 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()) @@ -343,8 +344,8 @@ def test_LinearTransformsMapping_apply(tmp_path, data_path, testdata_path): assert isinstance(hmc, nitl.LinearTransformsMapping) # Test-case: realign functional data on to sbref - nii = hmc.apply( - testdata_path / "func.nii.gz", order=1, reference=testdata_path / "sbref.nii.gz" + nii = apply( + hmc, testdata_path / "func.nii.gz", order=1, reference=testdata_path / "sbref.nii.gz" ) assert nii.dataobj.shape[-1] == len(hmc) @@ -352,13 +353,14 @@ def test_LinearTransformsMapping_apply(tmp_path, data_path, testdata_path): hmcinv = nitl.LinearTransformsMapping( np.linalg.inv(hmc.matrix), reference=testdata_path / "func.nii.gz" ) - nii = hmcinv.apply(testdata_path / "fmap.nii.gz", order=1) + 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): - hmc.apply( + apply( + hmc, testdata_path / "func.nii.gz", order=1, reference=testdata_path / "sbref.nii.gz", From 150c14579c66ddb776dc8cacf718ddffb9c117cd Mon Sep 17 00:00:00 2001 From: Julien Marabotto Date: Thu, 18 Apr 2024 10:43:53 +0200 Subject: [PATCH 05/14] Implement ndim Tranformbase -> now raises TypeError; Affine -> overshadows TransformBase's property and returns self._matrix.ndim + 1. To be applied to LinearTransformMapping (Currently TypeError) --- commandline_for_tests.txt | 7 +++++++ nitransforms/.DS_Store | Bin 0 -> 6148 bytes nitransforms/base.py | 16 ++++++++++++++-- nitransforms/io/.DS_Store | Bin 0 -> 6148 bytes nitransforms/linear.py | 9 +++++++-- nitransforms/tests/test_base.py | 11 +++++++++-- 6 files changed, 37 insertions(+), 6 deletions(-) create mode 100644 commandline_for_tests.txt create mode 100644 nitransforms/.DS_Store create mode 100644 nitransforms/io/.DS_Store diff --git a/commandline_for_tests.txt b/commandline_for_tests.txt new file mode 100644 index 00000000..87c85ea2 --- /dev/null +++ b/commandline_for_tests.txt @@ -0,0 +1,7 @@ +Terminal command line to test code: + +>>> TEST_DATA_HOME=$HOME/Desktop/nitransforms-tests python3 -m pytest + +Also works for a path to a specific test, eg: + +>>> TEST_DATA_HOME=$HOME/Desktop/nitransforms-tests python3 -m pytest ~/workspace/nitransforms/nitransforms/tests/test_base.py::test_TransformBase \ No newline at end of file diff --git a/nitransforms/.DS_Store b/nitransforms/.DS_Store new file mode 100644 index 0000000000000000000000000000000000000000..5008ddfcf53c02e82d7eee2e57c38e5672ef89f6 GIT binary patch literal 6148 zcmeH~Jr2S!425mzP>H1@V-^m;4Wg<&0T*E43hX&L&p$$qDprKhvt+--jT7}7np#A3 zem<@ulZcFPQ@L2!n>{z**++&mCkOWA81W14cNZlEfg7;MkzE(HCqgga^y>{tEnwC%0;vJ&^%eQ zLs35+`xjp>T0VePX7uvh|E6Wc&!U>a0lP&G#k4Lb5A>uO>f7WZy#UK86urz7ulAb$o-7aA4#Zv`F!AQhhg literal 0 HcmV?d00001 diff --git a/nitransforms/linear.py b/nitransforms/linear.py index eb4a95d7..569ac9c2 100644 --- a/nitransforms/linear.py +++ b/nitransforms/linear.py @@ -143,6 +143,11 @@ def matrix(self): """Access the internal representation of this affine.""" return self._matrix + @property + def ndim(self): + """Access the internal representation of this affine.""" + return self._matrix.ndim + 1 + def map(self, x, inverse=False): r""" Apply :math:`y = f(x)`. @@ -270,7 +275,7 @@ def __repr__(self): """ return repr(self.matrix) - + class LinearTransformsMapping(Affine): """Represents a series of linear transforms.""" @@ -533,4 +538,4 @@ def load(filename, fmt=None, reference=None, moving=None): if isinstance(xfm, LinearTransformsMapping) and len(xfm) == 1: xfm = xfm[0] - return xfm + return xfm \ No newline at end of file diff --git a/nitransforms/tests/test_base.py b/nitransforms/tests/test_base.py index d09096f7..893fe80c 100644 --- a/nitransforms/tests/test_base.py +++ b/nitransforms/tests/test_base.py @@ -95,7 +95,8 @@ def _to_hdf5(klass, x5_root): # Test identity transform xfm = TransformBase() xfm.reference = fname - assert xfm.ndim == 3 + with pytest.raises(TypeError): + _ = xfm.ndim moved = apply(xfm, fname, order=0) assert np.all( imgdata == np.asanyarray(moved.dataobj, dtype=moved.get_data_dtype()) @@ -104,12 +105,18 @@ def _to_hdf5(klass, x5_root): # Test identity transform - setting reference xfm = TransformBase() xfm.reference = fname - assert xfm.ndim == 3 + with pytest.raises(TypeError): + _ = xfm.ndim moved = apply(xfm, str(fname), reference=fname, order=0) assert np.all( imgdata == np.asanyarray(moved.dataobj, dtype=moved.get_data_dtype()) ) + #Test ndim returned by affine + #breakpoint() + assert nitl.Affine().ndim == 3 + assert nitl.LinearTransformsMapping(xfm.matrix).ndim == 4 + # Test applying to Gifti gii = nb.gifti.GiftiImage( darrays=[ From f3529ced803ed1fef95034629d8de9c295f171db Mon Sep 17 00:00:00 2001 From: Julien Marabotto Date: Thu, 18 Apr 2024 15:57:33 +0200 Subject: [PATCH 06/14] ENH: Implement_ndim Added two tests for ndim: assert nitl.Affine().ndim == 3; assert nitl.LinearTransformsMapping([nitl.Affine()]).ndim == 4 . --- nitransforms/tests/test_base.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/nitransforms/tests/test_base.py b/nitransforms/tests/test_base.py index 893fe80c..97c452c6 100644 --- a/nitransforms/tests/test_base.py +++ b/nitransforms/tests/test_base.py @@ -113,9 +113,8 @@ def _to_hdf5(klass, x5_root): ) #Test ndim returned by affine - #breakpoint() assert nitl.Affine().ndim == 3 - assert nitl.LinearTransformsMapping(xfm.matrix).ndim == 4 + assert nitl.LinearTransformsMapping([nitl.Affine()]).ndim == 4 # Test applying to Gifti gii = nb.gifti.GiftiImage( From 3d40b93d1f2bc9ac85d25475189b95fb8704ac7a Mon Sep 17 00:00:00 2001 From: Julien Marabotto Date: Thu, 18 Apr 2024 16:44:33 +0200 Subject: [PATCH 07/14] enh: add a test on ndim --- nitransforms/base.py | 1 - nitransforms/tests/test_base.py | 4 +++- 2 files changed, 3 insertions(+), 2 deletions(-) diff --git a/nitransforms/base.py b/nitransforms/base.py index 42119d9c..05f72d7a 100644 --- a/nitransforms/base.py +++ b/nitransforms/base.py @@ -15,7 +15,6 @@ from nibabel import funcs as _nbfuncs from nibabel.nifti1 import intent_codes as INTENT_CODES from nibabel.cifti2 import Cifti2Image -from nitransforms import linear as nitl EQUALITY_TOL = 1e-5 diff --git a/nitransforms/tests/test_base.py b/nitransforms/tests/test_base.py index 97c452c6..5b50e656 100644 --- a/nitransforms/tests/test_base.py +++ b/nitransforms/tests/test_base.py @@ -114,7 +114,9 @@ def _to_hdf5(klass, x5_root): #Test ndim returned by affine assert nitl.Affine().ndim == 3 - assert nitl.LinearTransformsMapping([nitl.Affine()]).ndim == 4 + assert nitl.LinearTransformsMapping( + [nitl.Affine(), nitl.Affine()] + ).ndim == 4 # Test applying to Gifti gii = nb.gifti.GiftiImage( From abe016ef9b6c9c30a51692a5da6a17ae89b5aba7 Mon Sep 17 00:00:00 2001 From: Julien Marabotto Date: Thu, 18 Apr 2024 16:55:30 +0200 Subject: [PATCH 08/14] fix: remove straneous files, revert changes from other branch --- commandline_for_tests.txt | 7 --- nitransforms/.DS_Store | Bin 6148 -> 0 bytes nitransforms/base.py | 92 +++++++++++++++++++++++++++++ nitransforms/resampling.py | 116 ------------------------------------- 4 files changed, 92 insertions(+), 123 deletions(-) delete mode 100644 commandline_for_tests.txt delete mode 100644 nitransforms/.DS_Store delete mode 100644 nitransforms/resampling.py diff --git a/commandline_for_tests.txt b/commandline_for_tests.txt deleted file mode 100644 index 87c85ea2..00000000 --- a/commandline_for_tests.txt +++ /dev/null @@ -1,7 +0,0 @@ -Terminal command line to test code: - ->>> TEST_DATA_HOME=$HOME/Desktop/nitransforms-tests python3 -m pytest - -Also works for a path to a specific test, eg: - ->>> TEST_DATA_HOME=$HOME/Desktop/nitransforms-tests python3 -m pytest ~/workspace/nitransforms/nitransforms/tests/test_base.py::test_TransformBase \ No newline at end of file diff --git a/nitransforms/.DS_Store b/nitransforms/.DS_Store deleted file mode 100644 index 5008ddfcf53c02e82d7eee2e57c38e5672ef89f6..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 6148 zcmeH~Jr2S!425mzP>H1@V-^m;4Wg<&0T*E43hX&L&p$$qDprKhvt+--jT7}7np#A3 zem<@ulZcFPQ@L2!n>{z**++&mCkOWA81W14cNZlEfg7;MkzE(HCqgga^y>{tEnwC%0;vJ&^%eQ zLs35+`xjp>T0 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. + output_dtype: dtype specifier, 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 + dtype is float64, otherwise it is equivalent to the input image's + ``get_data_dtype()`` (on-disk type). + 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. + Returns + ------- + resampled : `spatialimage` or ndarray + The data imaged after resampling to reference space. + """ + if reference is not None and isinstance(reference, (str, Path)): + reference = _nbload(str(reference)) + + _ref = ( + self.reference if reference is None else SpatialReference.factory(reference) + ) + + if _ref is None: + raise TransformError("Cannot apply transform without reference") + + if isinstance(spatialimage, (str, Path)): + spatialimage = _nbload(str(spatialimage)) + + data = np.asanyarray(spatialimage.dataobj) + targets = ImageGrid(spatialimage).index( # data should be an image + _as_homogeneous(self.map(_ref.ndcoords.T), dim=_ref.ndim) + ) + + resampled = ndi.map_coordinates( + data, + targets.T, + output=output_dtype, + order=order, + mode=mode, + cval=cval, + prefilter=prefilter, + ) + + 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()) + moved = spatialimage.__class__( + resampled.reshape(_ref.shape), + _ref.affine, + hdr, + ) + return moved + + return resampled + def map(self, x, inverse=False): r""" Apply :math:`y = f(x)`. diff --git a/nitransforms/resampling.py b/nitransforms/resampling.py deleted file mode 100644 index e89b081a..00000000 --- a/nitransforms/resampling.py +++ /dev/null @@ -1,116 +0,0 @@ -# emacs: -*- mode: python-mode; py-indent-offset: 4; indent-tabs-mode: nil -*- -# vi: set ft=python sts=4 ts=4 sw=4 et: -### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ## -# -# See COPYING file distributed along with the NiBabel package for the -# copyright and license terms. -# -### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ## -"""Resampling utilities.""" -from pathlib import Path -import numpy as np -from nibabel.loadsave import load as _nbload -from scipy import ndimage as ndi - -from nitransforms.base import ( - ImageGrid, - TransformError, - SpatialReference, - _as_homogeneous, -) - - -def apply( - transform, - spatialimage, - reference=None, - order=3, - mode="constant", - cval=0.0, - prefilter=True, - output_dtype=None, -): - """ - Apply a transformation to an image, resampling on the reference spatial object. - - Parameters - ---------- - spatialimage : `spatialimage` - The image object containing the data to be resampled in reference - space - reference : spatial object, optional - The image, surface, or combination thereof containing the coordinates - of samples that will be sampled. - order : 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 - Determines how the input image is extended when the resamplings overflows - a border. Default is 'constant'. - cval : float, optional - Constant value for ``mode='constant'``. Default is 0.0. - prefilter: 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. - output_dtype: dtype specifier, 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 - dtype is float64, otherwise it is equivalent to the input image's - ``get_data_dtype()`` (on-disk type). - 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. - - Returns - ------- - resampled : `spatialimage` or ndarray - The data imaged after resampling to reference space. - - """ - if reference is not None and isinstance(reference, (str, Path)): - reference = _nbload(str(reference)) - - _ref = ( - transform.reference if reference is None else SpatialReference.factory(reference) - ) - - if _ref is None: - raise TransformError("Cannot apply transform without reference") - - if isinstance(spatialimage, (str, Path)): - spatialimage = _nbload(str(spatialimage)) - - data = np.asanyarray(spatialimage.dataobj) - targets = ImageGrid(spatialimage).index( # data should be an image - _as_homogeneous(transform.map(_ref.ndcoords.T), dim=_ref.ndim) - ) - - resampled = ndi.map_coordinates( - data, - targets.T, - output=output_dtype, - order=order, - mode=mode, - cval=cval, - prefilter=prefilter, - ) - - 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()) - moved = spatialimage.__class__( - resampled.reshape(_ref.shape), - _ref.affine, - hdr, - ) - return moved - - return resampled From 180ec18fdb0d8a2c1a18e6777d3807530130a78a Mon Sep 17 00:00:00 2001 From: Julien Marabotto Date: Thu, 18 Apr 2024 16:59:39 +0200 Subject: [PATCH 09/14] enh: revert unrelated changes --- nitransforms/tests/test_base.py | 7 +++--- nitransforms/tests/test_io.py | 41 +++++++++++-------------------- nitransforms/tests/test_linear.py | 16 ++++++------ 3 files changed, 24 insertions(+), 40 deletions(-) diff --git a/nitransforms/tests/test_base.py b/nitransforms/tests/test_base.py index 5b50e656..f315cc6c 100644 --- a/nitransforms/tests/test_base.py +++ b/nitransforms/tests/test_base.py @@ -6,7 +6,6 @@ from ..base import SpatialReference, SampledSpatialData, ImageGrid, TransformBase from .. import linear as nitl -from ..resampling import apply def test_SpatialReference(testdata_path): @@ -97,7 +96,7 @@ def _to_hdf5(klass, x5_root): xfm.reference = fname with pytest.raises(TypeError): _ = xfm.ndim - moved = apply(xfm, fname, order=0) + moved = xfm.apply(fname, order=0) assert np.all( imgdata == np.asanyarray(moved.dataobj, dtype=moved.get_data_dtype()) ) @@ -107,7 +106,7 @@ def _to_hdf5(klass, x5_root): xfm.reference = fname with pytest.raises(TypeError): _ = xfm.ndim - moved = apply(xfm, str(fname), reference=fname, order=0) + moved = xfm.apply(str(fname), reference=fname, order=0) assert np.all( imgdata == np.asanyarray(moved.dataobj, dtype=moved.get_data_dtype()) ) @@ -127,7 +126,7 @@ def _to_hdf5(klass, x5_root): ) ] ) - giimoved = apply(xfm, fname, reference=gii, order=0) + giimoved = xfm.apply(fname, reference=gii, order=0) assert np.allclose(giimoved.reshape(xfm.reference.shape), moved.get_fdata()) # Test to_filename diff --git a/nitransforms/tests/test_io.py b/nitransforms/tests/test_io.py index 0cc79d15..bcee9198 100644 --- a/nitransforms/tests/test_io.py +++ b/nitransforms/tests/test_io.py @@ -28,8 +28,6 @@ ) from nitransforms.io.base import LinearParameters, TransformIOError, TransformFileError from nitransforms.conftest import _datadir, _testdir -from nitransforms.resampling import apply - LPS = np.diag([-1, -1, 1, 1]) ITK_MAT = LPS.dot(np.ones((4, 4)).dot(LPS)) @@ -499,13 +497,10 @@ def test_afni_oblique(tmpdir, parameters, swapaxes, testdata_path, dir_x, dir_y, assert np.allclose(card_aff, nb.load("deob_3drefit.nii.gz").affine) # Check that nitransforms can emulate 3drefit -deoblique - nt3drefit = apply( - Affine( - afni._cardinal_rotation(img.affine, False), - reference="deob_3drefit.nii.gz", - ), - "orig.nii.gz", - ) + nt3drefit = Affine( + afni._cardinal_rotation(img.affine, False), + reference="deob_3drefit.nii.gz", + ).apply("orig.nii.gz") diff = ( np.asanyarray(img.dataobj, dtype="uint8") @@ -514,13 +509,10 @@ def test_afni_oblique(tmpdir, parameters, swapaxes, testdata_path, dir_x, dir_y, assert np.sqrt((diff[10:-10, 10:-10, 10:-10] ** 2).mean()) < 0.1 # Check that nitransforms can revert 3drefit -deoblique - nt_undo3drefit = apply( - Affine( - afni._cardinal_rotation(img.affine, True), - reference="orig.nii.gz", - ), - "deob_3drefit.nii.gz", - ) + nt_undo3drefit = Affine( + afni._cardinal_rotation(img.affine, True), + reference="orig.nii.gz", + ).apply("deob_3drefit.nii.gz") diff = ( np.asanyarray(img.dataobj, dtype="uint8") @@ -539,21 +531,16 @@ def test_afni_oblique(tmpdir, parameters, swapaxes, testdata_path, dir_x, dir_y, assert np.allclose(deobaff, deobnii.affine) # Check resampling in deobliqued grid - ntdeobnii = apply( - Affine(np.eye(4), reference=deobnii.__class__( - np.zeros(deobshape, dtype="uint8"), - deobaff, - deobnii.header - )), - img, - order=0, - ) + ntdeobnii = Affine(np.eye(4), reference=deobnii.__class__( + np.zeros(deobshape, dtype="uint8"), + deobaff, + deobnii.header + )).apply(img, order=0) # Generate an internal box to exclude border effects box = np.zeros(img.shape, dtype="uint8") box[10:-10, 10:-10, 10:-10] = 1 - ntdeobmask = apply( - Affine(np.eye(4), reference=ntdeobnii), + ntdeobmask = Affine(np.eye(4), reference=ntdeobnii).apply( nb.Nifti1Image(box, img.affine, img.header), order=0, ) diff --git a/nitransforms/tests/test_linear.py b/nitransforms/tests/test_linear.py index 9a06fe32..2957f59c 100644 --- a/nitransforms/tests/test_linear.py +++ b/nitransforms/tests/test_linear.py @@ -13,7 +13,6 @@ 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 @@ -286,7 +285,7 @@ def test_apply_linear_transform(tmpdir, get_testdata, get_testmask, image_orient assert exit_code == 0 sw_moved_mask = nb.load("resampled_brainmask.nii.gz") - nt_moved_mask = apply(xfm, msk, order=0) + nt_moved_mask = xfm.apply(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) @@ -306,7 +305,7 @@ def test_apply_linear_transform(tmpdir, get_testdata, get_testmask, image_orient sw_moved = nb.load("resampled.nii.gz") sw_moved.set_data_dtype(img.get_data_dtype()) - nt_moved = apply(xfm, img, order=0) + nt_moved = xfm.apply(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()) @@ -315,7 +314,7 @@ def test_apply_linear_transform(tmpdir, get_testdata, get_testmask, image_orient # 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) + nt_moved = xfm.apply("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()) @@ -344,8 +343,8 @@ def test_LinearTransformsMapping_apply(tmp_path, data_path, testdata_path): 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" + nii = hmc.apply( + testdata_path / "func.nii.gz", order=1, reference=testdata_path / "sbref.nii.gz" ) assert nii.dataobj.shape[-1] == len(hmc) @@ -353,14 +352,13 @@ def test_LinearTransformsMapping_apply(tmp_path, data_path, testdata_path): hmcinv = nitl.LinearTransformsMapping( np.linalg.inv(hmc.matrix), reference=testdata_path / "func.nii.gz" ) - nii = apply(hmcinv, testdata_path / "fmap.nii.gz", order=1) + nii = hmcinv.apply(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, + hmc.apply( testdata_path / "func.nii.gz", order=1, reference=testdata_path / "sbref.nii.gz", From c0e2a98fe3459243ecb3ac6aaba18fd91c209b42 Mon Sep 17 00:00:00 2001 From: Julien Marabotto Date: Thu, 18 Apr 2024 17:00:55 +0200 Subject: [PATCH 10/14] enh: revert unrelated changes --- nitransforms/cli.py | 4 +--- nitransforms/io/.DS_Store | Bin 6148 -> 0 bytes 2 files changed, 1 insertion(+), 3 deletions(-) delete mode 100644 nitransforms/io/.DS_Store diff --git a/nitransforms/cli.py b/nitransforms/cli.py index 8f8f5ce0..63b8bed4 100644 --- a/nitransforms/cli.py +++ b/nitransforms/cli.py @@ -5,7 +5,6 @@ from .linear import load as linload from .nonlinear import load as nlinload -from .resampling import apply def cli_apply(pargs): @@ -39,8 +38,7 @@ def cli_apply(pargs): # ensure a reference is set xfm.reference = pargs.ref or pargs.moving - moved = apply( - xfm, + moved = xfm.apply( pargs.moving, order=pargs.order, mode=pargs.mode, diff --git a/nitransforms/io/.DS_Store b/nitransforms/io/.DS_Store deleted file mode 100644 index db6d5fb599944be892b217b8202a87b66720ca2f..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 6148 zcmeHKJBk895Uo-J!Nkz8uiyVePX7uvh|E6Wc&!U>a0lP&G#k4Lb5A>uO>f7WZy#UK86urz7ulAb$o-7aA4#Zv`F!AQhhg From 90e529ba664e54d52a6fae914dc3b9fc309e2190 Mon Sep 17 00:00:00 2001 From: Julien Marabotto Date: Thu, 18 Apr 2024 17:03:52 +0200 Subject: [PATCH 11/14] enh: revert unrelated changes --- nitransforms/base.py | 14 ++++---------- 1 file changed, 4 insertions(+), 10 deletions(-) diff --git a/nitransforms/base.py b/nitransforms/base.py index 66b3808f..117d4064 100644 --- a/nitransforms/base.py +++ b/nitransforms/base.py @@ -15,6 +15,7 @@ from nibabel import funcs as _nbfuncs from nibabel.nifti1 import intent_codes as INTENT_CODES from nibabel.cifti2 import Cifti2Image +from scipy import ndimage as ndi EQUALITY_TOL = 1e-5 @@ -212,13 +213,6 @@ def reference(self): if self._reference is None: warnings.warn("Reference space not set") return self._reference - - @property - def reference(_ndim): - """Access a reference space where data will be resampled onto.""" - if _ndim._reference is None: - warnings.warn("Reference space not set") - return _ndim._reference @reference.setter def reference(self, image): @@ -227,7 +221,6 @@ def reference(self, image): @property def ndim(self): """Access the dimensions of the reference space.""" - #return self.reference.ndim raise TypeError("TransformBase has no dimensions") def apply( @@ -242,6 +235,7 @@ def apply( ): """ Apply a transformation to an image, resampling on the reference spatial object. + Parameters ---------- spatialimage : `spatialimage` @@ -275,10 +269,12 @@ 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. + Returns ------- resampled : `spatialimage` or ndarray The data imaged after resampling to reference space. + """ if reference is not None and isinstance(reference, (str, Path)): reference = _nbload(str(reference)) @@ -388,5 +384,3 @@ def _as_homogeneous(xyz, dtype="float32", dim=3): def _apply_affine(x, affine, dim): """Get the image array's indexes corresponding to coordinates.""" return affine.dot(_as_homogeneous(x, dim=dim).T)[:dim, ...].T - -#import pdb; pdb.set_trace() From bafd2ed9a6248779c5c682a603d17baf9585457f Mon Sep 17 00:00:00 2001 From: Julien Marabotto Date: Thu, 18 Apr 2024 17:06:23 +0200 Subject: [PATCH 12/14] enh: revert unrelated changes --- nitransforms/base.py | 1 - nitransforms/linear.py | 4 ++-- 2 files changed, 2 insertions(+), 3 deletions(-) diff --git a/nitransforms/base.py b/nitransforms/base.py index 117d4064..96f00edb 100644 --- a/nitransforms/base.py +++ b/nitransforms/base.py @@ -180,7 +180,6 @@ class TransformBase: __slots__ = ("_reference", "_ndim",) - def __init__(self, reference=None): """Instantiate a transform.""" self._reference = None diff --git a/nitransforms/linear.py b/nitransforms/linear.py index 569ac9c2..c7eb608a 100644 --- a/nitransforms/linear.py +++ b/nitransforms/linear.py @@ -275,7 +275,7 @@ def __repr__(self): """ return repr(self.matrix) - + class LinearTransformsMapping(Affine): """Represents a series of linear transforms.""" @@ -538,4 +538,4 @@ def load(filename, fmt=None, reference=None, moving=None): if isinstance(xfm, LinearTransformsMapping) and len(xfm) == 1: xfm = xfm[0] - return xfm \ No newline at end of file + return xfm From acb0736291a17dd3af46a25e2ca7efc65762b68e Mon Sep 17 00:00:00 2001 From: Oscar Esteban Date: Fri, 19 Apr 2024 10:15:32 +0200 Subject: [PATCH 13/14] Update nitransforms/tests/test_base.py --- nitransforms/tests/test_base.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/nitransforms/tests/test_base.py b/nitransforms/tests/test_base.py index f315cc6c..07a7e4ec 100644 --- a/nitransforms/tests/test_base.py +++ b/nitransforms/tests/test_base.py @@ -111,7 +111,7 @@ def _to_hdf5(klass, x5_root): imgdata == np.asanyarray(moved.dataobj, dtype=moved.get_data_dtype()) ) - #Test ndim returned by affine + # Test ndim returned by affine assert nitl.Affine().ndim == 3 assert nitl.LinearTransformsMapping( [nitl.Affine(), nitl.Affine()] From 4a182000449165c45e3a8617ba15a8c079a576c3 Mon Sep 17 00:00:00 2001 From: Oscar Esteban Date: Fri, 19 Apr 2024 10:17:27 +0200 Subject: [PATCH 14/14] Update nitransforms/linear.py --- nitransforms/linear.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/nitransforms/linear.py b/nitransforms/linear.py index c7eb608a..af14f396 100644 --- a/nitransforms/linear.py +++ b/nitransforms/linear.py @@ -147,7 +147,7 @@ def matrix(self): def ndim(self): """Access the internal representation of this affine.""" return self._matrix.ndim + 1 - + def map(self, x, inverse=False): r""" Apply :math:`y = f(x)`.