Dirty example to display signal for individual shell of multi-shell as spherical function.

signal_multi_shell.py

"""
Dirty example to display signal for individual shell of multi-shell as spherical function.
No background maps, just signal spheres.
Option to assume or not antipodal symmetry.
Assume the gradient table as at least 1 b0.
Uses so not very robust ad-hoc function to detect which sample belong to which shell.
"""

import numpy as np

import nibabel as nib

from dipy.viz import fvtk
from dipy.core.sphere import HemiSphere, Sphere
from dipy.io.gradients import read_bvals_bvecs
# from dipy.core.gradients import gradient_table

from dipy.data import fetch_isbi2013_2shell, read_isbi2013_2shell



def main():
	# # CHANGE THIS
	# path_dwi = '/media/Data/data_MRI/slider_kawin/'
	# path_bvecs = '/media/Data/data_MRI/slider_kawin/'
	# path_bvals = '/media/Data/data_MRI/slider_kawin/'
	apply_antipodal_symmetry = True
	shell_idx_to_plot = 0 # excluding b0s


	# # Get data / gradient info
	# data = nib.load(path_dwi).get_data()
	# bvals, bvecs = read_bvals_bvecs(path_bvals, path_bvecs)



	fetch_isbi2013_2shell()
	img, gtab = read_isbi2013_2shell()
	data = img.get_data()
	data= data[10:40, 22, None, 10:40]
	bvecs = gtab.bvecs
	bvals = gtab.bvals
	del gtab, img


	# At this point, you need to have bvecs/bvals/data, 
	# comment the fetching and uncomment + change path for your data



	# Find b-values in the sampling scheme
	target_bvalues = _find_target_bvalues(bvals, shell_th = 50)

	# Drop the first (assumes it's a b0)
	target_bvalues = target_bvalues[1:]

	# Assign shell idx to bvecs based on target_bvalues
	shells = _find_shells(bvals, target_bvalues, shell_th = 50)

	# Extract desired shell
	bvecs_to_use = (shells == shell_idx_to_plot)
	bvecs_shell = bvecs[bvecs_to_use]
	data_shell = data[..., bvecs_to_use]


	# Making a sphere to get vertices and faces
	if apply_antipodal_symmetry:
		signal_sphere = HemiSphere(xyz = bvecs_shell)
		signal_sphere = signal_sphere.mirror()
		data_shell = np.concatenate((data_shell, data_shell), axis = -1)
	else:
		signal_sphere = Sphere(xyz = bvecs_shell)


	r = fvtk.ren()
	sfu = fvtk.sphere_funcs(data_shell, signal_sphere, scale=2.2, norm=True)
	sfu.RotateX(90)
	sfu.RotateY(180)
	fvtk.add(r, sfu)
	# outname = 'screenshot_signal.png'
	# fvtk.record(r, n_frames=1, out_path = outname, size=(3000, 1500), magnification = 2)
	fvtk.show(r)



## SOON: from scilpy.viz.sampling_scheme import ...
def _find_target_bvalues(bvals, shell_th = 50):
	# Not robust
	target_bvalues = []

	bvalues = np.sort(np.array(list(set(bvals))))

	for bval in bvalues:
		add_bval = True
		for target_bval in target_bvalues:
			if (bval <= target_bval + shell_th) & (bval >= target_bval - shell_th):
				add_bval = False
		if add_bval:
			target_bvalues.append(bval)

	return target_bvalues


# Assign bvecs to a target shell
def _find_shells(bvals, target_bvalues, shell_th = 50):
	# Not robust
	# shell -1 means nbvecs not part of target_bvalues
	shells = -1 * np.ones_like(bvals)

	for shell_id, bval in enumerate(target_bvalues):
		shells[(bvals <= bval + shell_th) & (bvals >= bval - shell_th)] = shell_id

	return shells


if __name__ == "__main__":
    main()