ENH: testing DTI (no kurtosis)

single- and multi-shell
exact fit (no noise)

unit_tests/fdt/dtifit/feedsRun

+#!/usr/bin/env fslpython
+"""
+Generates data following the diffusion tensor model
+"""
+import sys
+import os
+from subprocess import run
+import numpy as np
+from numpy import testing
+import nibabel as nib
+
+OUTDIR = sys.argv[0]
+
+
+def gen_data():
+    """
+    Populates given directory with diffusion data
+
+    :yield: directories containing reference noiseless data with the eigen-vectors and eigen-values used to generate the data
+    """
+    directory = 'dti'
+
+    # eigen vectors vary along x-axis
+    eigen_vectors = np.array([
+        [(1., 0, 0), (0, 1, 0), (0, 0, 1)],
+        [(1., 0, 0), (0, np.sqrt(0.5), np.sqrt(0.5)), (0, np.sqrt(0.5), -np.sqrt(0.5))],
+        [(np.sqrt(0.5), -np.sqrt(0.5), 0), (np.sqrt(0.25), np.sqrt(0.25), np.sqrt(0.5)),
+         (np.sqrt(0.25), np.sqrt(0.25), -np.sqrt(0.5))],
+    ])[:, None, None, :, :]
+
+    for idx1 in range(3):
+        for idx2 in range(3):
+            testing.assert_allclose((eigen_vectors[..., idx1, :] * eigen_vectors[..., idx2, :]).sum(-1),
+                                    1 if idx1 == idx2 else 0., atol=1e-8)
+
+    # eigen-values vary along y-axis
+    # eigen-values can not be the same or eigen-vectors will be ill-defined
+    eigen_values = np.array([
+        [1.2, 1., 0.8],
+        [1., 0.5, 0.],
+        [2., 1.9, 1.8],
+    ])[None, :, None, :, None] * 1e-3
+
+    # S0 varies along the z-axis
+    S0 = np.array([1., 1.5])[None, None, :, None, None] * 1000.
+
+    eigen_vectors, eigen_values, S0 = np.broadcast_arrays(eigen_vectors, eigen_values, S0)
+    S0 = S0[:, :, :, 0, 0]
+
+    diffusion_tensor = (eigen_values[..., None] * eigen_vectors[..., None, :] * eigen_vectors[..., :, None]).sum(-3)
+
+    for multi_shell in (False, True):  # we should get the same results for single or multi-shell data
+        directory += '_multi' if multi_shell else '_single'
+        for kurt in (0, ):#1, 2):  # 0: no kurtosis; 1: single kurtosis; 2: parallel and perpendicular kurtosis
+            if kurt != 0 and not multi_shell:
+                continue
+            directory += ['', '_kurt', '_kurtdir'][kurt]
+
+            bvals = np.full(50, 1000.)
+            if multi_shell:
+                bvals[25:] = 2000.
+            bvals[::10] = 0.
+
+            run(['gps', '--ndir=50', '--out=bvecs'], check=True)
+            bvecs = np.genfromtxt('bvecs')
+            assert bvecs.shape == (50, 3), f"GPS produced bvecs-file with shape {bvecs.shape} rather than the expected (50, 3)"
+            testing.assert_allclose((bvecs ** 2).sum(-1), 1.)
+
+            kurt_val = np.array([
+                [0, 0, 0.],
+                [1., 1., 1.],
+                [1., 0.5, 0.5],
+            ][kurt])
+            data = S0[..., None] * np.exp(np.sum(
+                -(bvals[:, None] * eigen_values[..., None, :, 0] +
+                  kurt_val * eigen_values[..., None, :, 0] ** 2 * bvals[:, None] ** 2) *
+                np.sum(bvecs[:, None, :] * eigen_vectors[..., None, :, :], -1) ** 2, axis=-1
+            ))
+
+            if (data / S0[..., None]).min() < 0.01:
+                raise ValueError("dtifit rounds attenuations below 0.01 up to 0.01, "
+                                 "so these low values should not be in the test data")
+
+            for flipped in (False, True):
+                if flipped:
+                    directory += '_flipped'
+
+                if not os.path.isdir(directory):
+                    os.mkdir(directory)
+
+                if flipped:
+                    np.savetxt(f'{directory}/bvals', bvals[:, None])
+                    np.savetxt(f'{directory}/bvecs', bvecs.T)
+                else:
+                    np.savetxt(f'{directory}/bvals', bvals[None, :])
+                    np.savetxt(f'{directory}/bvecs', bvecs)
+
+                affine = np.eye(4) * 1.25
+                affine[-1, -1] = 1.
+                for idx in range(3):
+                    nib.Nifti1Image(eigen_vectors[..., idx, :], affine=affine).to_filename(f'{directory}/ref_V{idx + 1}.nii.gz')
+                    nib.Nifti1Image(eigen_values[..., idx, 0], affine=affine).to_filename(f'{directory}/ref_L{idx + 1}.nii.gz')
+                nib.Nifti1Image(S0, affine=affine).to_filename(f'{directory}/ref_S0.nii.gz')
+                nib.Nifti1Image(data, affine=affine).to_filename(f'{directory}/ref_data.nii.gz')
+                nib.Nifti1Image(np.ones(data.shape[:3], dtype=int), affine=affine).to_filename(f'{directory}/nodif_brain_mask.nii.gz')
+
+                tensor_components = diffusion_tensor[:, :, :, [0, 0, 0, 1, 1, 2], [0, 1, 2, 1, 2, 2]]
+                nib.Nifti1Image(tensor_components, affine=affine).to_filename(f'{directory}/ref_tensor.nii.gz')
+                yield directory, kurt
+
+
+def fit_data(directory, kurt=0):
+    for wls in (False, True):  # in this noise-free data the --wls flag should not matter
+        base_output = f'{directory}/dti{"_wls" if wls else ""}'
+        cmd = [
+            'dtifit',
+            '-k', f'{directory}/ref_data.nii.gz',
+            '-m', f'{directory}/nodif_brain_mask.nii.gz',
+            '-r', f'{directory}/bvecs',
+            '-b', f'{directory}/bvals',
+            '-o', base_output,
+            '--sse',
+            '--save_tensor',
+        ]
+        if wls:
+            cmd += ['--wls']
+        if kurt == 1:
+            cmd += ['--kurt']
+        if kurt == 2:
+            cmd += ['--kurt']
+        run(cmd, check=True)
+        yield base_output
+
+
+for directory, kurt in gen_data():
+    for base_output in fit_data(directory, kurt):
+        def compare(name):
+            ref = nib.load(f'{directory}/ref_{name}.nii.gz').get_fdata()
+            fit = nib.load(f'{base_output}_{name}.nii.gz').get_fdata()
+            assert ref.shape == fit.shape, f'incorrect NIFTI image shape for {name}'
+            testing.assert_allclose(ref, fit, atol=1e-8,
+                                    err_msg=f'mismatch in {name}')
+
+        compare('S0')
+        compare('tensor')
+        for idx in (3, 2, 1):
+            compare(f'L{idx}')
+
+            ref = nib.load(f'{directory}/ref_V{idx}.nii.gz').get_fdata()
+            fit = nib.load(f'{base_output}_V{idx}.nii.gz').get_fdata()
+            assert ref.shape == fit.shape
+            inner = (ref * fit).sum(-1)
+            testing.assert_allclose(abs(inner), 1.)
+
+        sse = nib.load(f'{base_output}_sse.nii.gz').get_fdata()
+        testing.assert_allclose(sse, 0., atol=1e-8)
+
+