Merge pull request #26 from dPys/enh/vectors

ENH: Vector representation and checking utilities

Author: oesteban

.travis.yml

@@ -6,6 +6,7 @@ language: python
 cache:
   directories:
     - $HOME/.cache/pip
+    - $HOME/.cache/data
 
 python:
   - 3.6

dmriprep/conftest.py

@@ -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)
@@ -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

dmriprep/interfaces/vectors.py

@@ -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

dmriprep/utils/tests/__init__.py

@@ -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')