Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

[ENH] Vector-checking utilities #26

Merged
merged 13 commits into from
Nov 5, 2019
Merged
1 change: 1 addition & 0 deletions .travis.yml
Original file line number Diff line number Diff line change
Expand Up @@ -6,6 +6,7 @@ language: python
cache:
directories:
- $HOME/.cache/pip
- $HOME/.cache/data

python:
- 3.6
Expand Down
29 changes: 29 additions & 0 deletions dmriprep/conftest.py
Original file line number Diff line number Diff line change
Expand Up @@ -5,6 +5,28 @@
import numpy as np
import nibabel as nb
import pytest
from dipy.data.fetcher import _make_fetcher, UW_RW_URL

_dipy_datadir_root = os.getenv('DMRIPREP_TESTS_DATA') or Path.home()
dipy_datadir = Path(_dipy_datadir_root) / '.cache' / 'data'
dipy_datadir.mkdir(parents=True, exist_ok=True)

_make_fetcher(
"fetch_sherbrooke_3shell",
str(dipy_datadir),
UW_RW_URL + "1773/38475/",
['HARDI193.nii.gz', 'HARDI193.bval', 'HARDI193.bvec'],
['HARDI193.nii.gz', 'HARDI193.bval', 'HARDI193.bvec'],
['0b735e8f16695a37bfbd66aab136eb66',
'e9b9bb56252503ea49d31fb30a0ac637',
'0c83f7e8b917cd677ad58a078658ebb7'],
doc="Download a 3shell HARDI dataset with 192 gradient direction")()

_sherbrooke_data = {
'dwi_file': dipy_datadir / "HARDI193.nii.gz",
'bvecs': np.loadtxt(dipy_datadir / "HARDI193.bvec").T,
'bvals': np.loadtxt(dipy_datadir / "HARDI193.bval"),
}


@pytest.fixture(autouse=True)
Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is awesome.

Expand All @@ -15,7 +37,14 @@ def doctest_autoimport(doctest_namespace):
doctest_namespace['os'] = os
doctest_namespace['Path'] = Path
doctest_namespace['data_dir'] = Path(__file__).parent / 'data' / 'tests'
doctest_namespace['dipy_datadir'] = dipy_datadir
tmpdir = tempfile.TemporaryDirectory()
doctest_namespace['tmpdir'] = tmpdir.name
yield
tmpdir.cleanup()


@pytest.fixture()
def dipy_test_data(scope='session'):
"""Create a temporal directory shared across tests to pull data in."""
return _sherbrooke_data
99 changes: 99 additions & 0 deletions dmriprep/interfaces/vectors.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,99 @@
"""Handling the gradient table."""
from pathlib import Path
import numpy as np
from nipype.utils.filemanip import fname_presuffix
from nipype.interfaces.base import (
SimpleInterface, BaseInterfaceInputSpec, TraitedSpec,
File, traits, isdefined
)
from ..utils.vectors import DiffusionGradientTable, B0_THRESHOLD, BVEC_NORM_EPSILON


class _CheckGradientTableInputSpec(BaseInterfaceInputSpec):
dwi_file = File(exists=True, mandatory=True)
in_bvec = File(exists=True, xor=['in_rasb'])
in_bval = File(exists=True, xor=['in_rasb'])
in_rasb = File(exists=True, xor=['in_bval', 'in_bvec'])
b0_threshold = traits.Float(B0_THRESHOLD, usedefault=True)
bvec_norm_epsilon = traits.Float(BVEC_NORM_EPSILON, usedefault=True)
b_scale = traits.Bool(True, usedefault=True)


class _CheckGradientTableOutputSpec(TraitedSpec):
out_rasb = File(exists=True)
out_bval = File(exists=True)
out_bvec = File(exists=True)
full_sphere = traits.Bool()
pole = traits.Tuple(traits.Float, traits.Float, traits.Float)
b0_ixs = traits.List(traits.Int)


class CheckGradientTable(SimpleInterface):
"""
Ensure the correctness of the gradient table.

Example
-------

>>> os.chdir(tmpdir)
>>> check = CheckGradientTable(
... dwi_file=str(data_dir / 'dwi.nii.gz'),
... in_rasb=str(data_dir / 'dwi.tsv')).run()
>>> check.outputs.pole
(0.0, 0.0, 0.0)
>>> check.outputs.full_sphere
True

>>> check = CheckGradientTable(
... dwi_file=str(data_dir / 'dwi.nii.gz'),
... in_bvec=str(data_dir / 'bvec'),
... in_bval=str(data_dir / 'bval')).run()
>>> check.outputs.pole
(0.0, 0.0, 0.0)
>>> check.outputs.full_sphere
True
>>> newrasb = np.loadtxt(check.outputs.out_rasb, skiprows=1)
>>> oldrasb = np.loadtxt(str(data_dir / 'dwi.tsv'), skiprows=1)
>>> np.allclose(newrasb, oldrasb, rtol=1.e-3)
True

"""

input_spec = _CheckGradientTableInputSpec
output_spec = _CheckGradientTableOutputSpec

def _run_interface(self, runtime):
rasb_file = _undefined(self.inputs, 'in_rasb')

table = DiffusionGradientTable(
self.inputs.dwi_file,
bvecs=_undefined(self.inputs, 'in_bvec'),
bvals=_undefined(self.inputs, 'in_bval'),
rasb_file=rasb_file,
b_scale=self.inputs.b_scale,
bvec_norm_epsilon=self.inputs.bvec_norm_epsilon,
b0_threshold=self.inputs.b0_threshold,
)
pole = table.pole
self._results['pole'] = tuple(pole)
self._results['full_sphere'] = np.all(pole == 0.0)
self._results['b0_ixs'] = np.where(table.b0mask)[0].tolist()

cwd = Path(runtime.cwd).absolute()
if rasb_file is None:
rasb_file = fname_presuffix(
self.inputs.dwi_file, use_ext=False, suffix='.tsv',
newpath=str(cwd))
table.to_filename(rasb_file)
self._results['out_rasb'] = rasb_file
table.to_filename('%s/dwi' % cwd, filetype='fsl')
self._results['out_bval'] = str(cwd / 'dwi.bval')
self._results['out_bvec'] = str(cwd / 'dwi.bvec')
return runtime


def _undefined(objekt, name, default=None):
value = getattr(objekt, name)
if not isdefined(value):
return default
return value
Empty file.
91 changes: 91 additions & 0 deletions dmriprep/utils/tests/test_vectors.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,91 @@
"""Test vector utilities."""
import pytest
import numpy as np
from dmriprep.utils import vectors as v
from collections import namedtuple


def test_corruption(tmpdir, dipy_test_data, monkeypatch):
"""Check whether b-value rescaling is operational."""
tmpdir.chdir()

bvals = dipy_test_data['bvals']
bvecs = dipy_test_data['bvecs']

dgt = v.DiffusionGradientTable(**dipy_test_data)
affine = dgt.affine.copy()

# Test vector hemisphere coverage
assert np.all(dgt.pole == [0., 0., 0.])

dgt.to_filename('dwi.tsv')
dgt = v.DiffusionGradientTable(rasb_file='dwi.tsv')
assert dgt.normalized is False
with pytest.raises(TypeError):
dgt.to_filename('dwi', filetype='fsl') # You can do this iff the affine is set.

# check accessing obj.affine
dgt = v.DiffusionGradientTable(dwi_file=namedtuple('Affine', ['affine'])(affine))
assert np.all(dgt.affine == affine)
dgt = v.DiffusionGradientTable(dwi_file=affine)
assert np.all(dgt.affine == affine)

# Perform various corruption checks using synthetic corrupted bval-bvec.
dgt = v.DiffusionGradientTable()
dgt.bvecs = bvecs
with pytest.raises(ValueError):
dgt.bvals = bvals[:-1]

dgt = v.DiffusionGradientTable()
dgt.bvals = bvals
with pytest.raises(ValueError):
dgt.bvecs = bvecs[:-1]

# Missing b0
bval_no_b0 = bvals.copy()
bval_no_b0[0] = 51
with pytest.raises(ValueError):
dgt = v.DiffusionGradientTable(dwi_file=dipy_test_data['dwi_file'],
bvals=bval_no_b0, bvecs=bvecs)
bvec_no_b0 = bvecs.copy()
bvec_no_b0[0] = np.array([1.0, 0.0, 0.0])
with pytest.raises(ValueError):
dgt = v.DiffusionGradientTable(dwi_file=dipy_test_data['dwi_file'],
bvals=bvals, bvecs=bvec_no_b0)

# Corrupt b0 b-val
bval_odd_b0 = bvals.copy()
bval_odd_b0[bval_odd_b0 == 0] = 1e-8
dgt = v.DiffusionGradientTable(dwi_file=dipy_test_data['dwi_file'],
bvals=bval_odd_b0, bvecs=bvecs)
assert dgt.bvals[0] == 0

# Corrupt b0 b-vec
bvec_odd_b0 = bvecs.copy()
b0mask = np.all(bvec_odd_b0 == 0, axis=1)
bvec_odd_b0[b0mask] = [10, 10, 10]
dgt = v.DiffusionGradientTable(dwi_file=dipy_test_data['dwi_file'],
bvals=bvals, bvecs=bvec_odd_b0)
assert np.all(dgt.bvecs[b0mask] == [0., 0., 0.])

# Test normalization
bvecs_factor = 2.0 * bvecs
dgt = v.DiffusionGradientTable(dwi_file=dipy_test_data['dwi_file'],
bvals=bvals, bvecs=bvecs_factor)
assert -1.0 <= np.max(np.abs(dgt.gradients[..., :-1])) <= 1.0
assert dgt.normalized is True

def mock_func(*args, **kwargs):
return 'called!'

with monkeypatch.context() as m:
m.setattr(v, 'normalize_gradients', mock_func)
assert dgt.normalize() is None # Test nothing is executed.

with monkeypatch.context() as m:
m.setattr(v, 'bvecs2ras', mock_func)
assert dgt.generate_vecval() is None # Test nothing is executed.

# Miscellaneous tests
with pytest.raises(ValueError):
dgt.to_filename('path', filetype='mrtrix')
Loading