From eb103a1039b4f8a30eae71dfac775b14a1a34963 Mon Sep 17 00:00:00 2001
From: Michiel Cottaar <MichielCottaar@protonmail.com>
Date: Fri, 6 Dec 2019 11:59:48 +0000
Subject: [PATCH] ENH: testing DTI (no kurtosis)
single- and multi-shell
exact fit (no noise)
---
unit_tests/fdt/dtifit/feedsRun | 158 +++++++++++++++++++++++++++++++++
1 file changed, 158 insertions(+)
create mode 100755 unit_tests/fdt/dtifit/feedsRun
diff --git a/unit_tests/fdt/dtifit/feedsRun b/unit_tests/fdt/dtifit/feedsRun
new file mode 100755
index 0000000..15f5528
--- /dev/null
+++ b/unit_tests/fdt/dtifit/feedsRun
@@ -0,0 +1,158 @@
+#!/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)
+
+
--
GitLab