Merge branch 'issue-23-tensor-compartment' into 'master'

Issue 23 tensor compartment

Closes #23

See merge request neichert/DIVE!12

dive/compartments/tensor.py

+from .base import Parameter, BaseCompartment
+import numpy as np
+
+
+class Tensor(BaseCompartment):
+    """
+    Represents water with gaussian diffusion in 3D.
+    The tensor is assumed to be aligned along the z-axis ([0,0,1]).
+
+    """
+    name = 'tensor'
+    required_acquisition = ('bval', 'bvec')
+    parameters = (
+        Parameter("diffusivity_ax", 0., 3.1, 1.7),
+        Parameter("diffusivity_perp1", 0., 3.1, 1.),
+        Parameter("diffusivity_perp2", 0., 3.1, 0.5),
+        Parameter("S0", 0., 2., 1.),
+    )
+
+    @staticmethod
+    def predict(S0,
+                diffusivity_ax,
+                diffusivity_perp1,
+                diffusivity_perp2,
+                bval,
+                bvec):
+        """
+        Predicts signal of water in a tensor
+
+        Args:
+            S0: base signal at b=0
+            diffusivity_ax: diffusivity of the water parallel to the fibre direction
+            diffusivity_perp1: diffusivity of the water orthogonal to the fibre direction
+            diffusivity_perp2: diffusivity of the water along the other orthogonal 
+            bval: b-value of the acquisition
+            bvec: direction of the applied gradient array. (..., 3)
+
+            Here, we assume the axial diffusivity to be along the z-axis
+        """
+
+        D = np.array([[diffusivity_perp2, 0, 0],
+                      [0, diffusivity_perp1, 0],
+                      [0, 0, diffusivity_ax]])
+
+        decay = - bval * np.sum(bvec * (D @ bvec.T).T, axis=-1)
+
+        return S0 * np.exp(decay)

test/compartments/test_tensor.py

+from dive.compartments.tensor import Tensor
+from dive.compartments.stick import Stick
+from dive.compartments.ball import Ball
+import pytest
+import numpy as np
+
+
+def test_valid_tensors():
+    # No diffusion
+    dot_tensor = Tensor(S0=1.3,
+                        diffusivity_ax=0,
+                        diffusivity_perp1=0,
+                        diffusivity_perp2=0)
+    z_direction = np.array([0., 0., 1.]).T  # (N,3)
+
+    assert dot_tensor(bval=0., bvec=z_direction) == 1.3
+    assert dot_tensor(bval=1., bvec=z_direction) == 1.3
+    assert dot_tensor.values['S0'] == 1.3
+    assert dot_tensor.values['diffusivity_ax'] == 0
+    assert dot_tensor.values['diffusivity_perp1'] == 0
+    assert dot_tensor.values['diffusivity_perp1'] == 0
+
+    # Only axial diffusion
+    stick_tensor = Tensor(S0=1.3,
+                          diffusivity_ax=3.,
+                          diffusivity_perp1=0,
+                          diffusivity_perp2=0)
+    stick = Stick(S0=1.3,
+                  diffusivity=3)
+
+    assert stick_tensor(bval=0., bvec=z_direction) == 1.3
+    assert stick_tensor(bval=2., bvec=z_direction) == stick(bval=2., bvec=z_direction)
+
+    # Isotropic diffusion
+    ball_tensor = Tensor(S0=1.3,
+                         diffusivity_ax=3.,
+                         diffusivity_perp1=3,
+                         diffusivity_perp2=3)
+    ball = Ball(S0=1.3,
+                diffusivity=3)
+    assert ball_tensor(bval=0., bvec=z_direction) == 1.3
+    assert ball_tensor(bval=2., bvec=z_direction) == ball(bval=2.)
+
+    # Fast isotropic diffusion
+    fast_ball_tensor = Tensor(S0=1.3,
+                              diffusivity_ax=1e5,
+                              diffusivity_perp1=1e5,
+                              diffusivity_perp2=1e5)
+    assert fast_ball_tensor(bval=0., bvec=z_direction) == 1.3
+    assert fast_ball_tensor(bval=1., bvec=z_direction) < 1e-20
+
+    # Default tensor
+    default_tensor = Tensor()
+    assert default_tensor.values['S0'] == 1.
+    assert default_tensor.values['diffusivity_ax'] == 1.7
+    assert default_tensor.values['diffusivity_perp1'] == 1.
+    assert default_tensor.values['diffusivity_perp2'] == 0.5
+
+    # Test dimensions of signal based on multiple bvecs
+    mutliple_direction_tensor = Tensor()
+    bvals = np.array([100., 300., 500., 500., ]).T
+    bvecs = np.array([[1., 0., 0.],
+                      [0., 1., 0.],
+                      [0., 0., 1.],
+                      [1 / 3, 1 / 3, 1 / 3]])
+
+    assert np.squeeze(mutliple_direction_tensor(bval=bvals, bvec=bvecs)).shape[0] == 4
+
+
+def test_invalid_tensor():
+    with pytest.raises(ValueError):
+        Tensor(unkown_parameter=20)