Resample DWI signal single-shell

resample_dwi_signal.py

import numpy as np
import nibabel as nib


# shitty way to get 1 voxel of single-shell signal
from dipy.data import fetch_isbi2013_2shell, read_isbi2013_2shell
fetch_isbi2013_2shell()
img, gtab = read_isbi2013_2shell()
data = img.get_data()
del img

# find an interesting voxel
vox = np.unravel_index(np.argmax(data.std(axis = 3)), data.shape[:3])
shell_1 = (gtab.bvals <= 1501) & (gtab.bvals >= 1499)
data = data[vox]
pts = gtab.bvecs[shell_1]
data = data[shell_1]


# plot signal on "native" bvecs
from dipy.core.sphere import HemiSphere
signal_native_pts = HemiSphere(xyz = pts)
signal_native_pts_mirror = signal_native_pts.mirror()
data_mirror = np.concatenate((data, data), axis = -1)

from dipy.viz import fvtk
r = fvtk.ren()
sfu = fvtk.sphere_funcs(data_mirror, signal_native_pts_mirror, 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)



# Fit SH to signal
from dipy.reconst.shm import sph_harm_lookup, smooth_pinv

sh_order = 8
smooth=0.006
sph_harm_basis = sph_harm_lookup.get('mrtrix')
Ba, m, n = sph_harm_basis(sh_order, signal_native_pts.theta, signal_native_pts.phi)
L = -n * (n + 1)
invB = smooth_pinv(Ba, np.sqrt(smooth) * L)
data_sh = np.dot(data, invB.T)


# Resample
from dipy.data import get_sphere
sphere = get_sphere('repulsion724')

sph_harm_basis = sph_harm_lookup.get('mrtrix')
Ba, m, n = sph_harm_basis(sh_order, sphere.theta, sphere.phi)
data_resampled = np.dot(Ba, data_sh)


# viz resampled signal
from dipy.viz import fvtk
r = fvtk.ren()
sfu = fvtk.sphere_funcs(data_resampled, sphere, scale=2.2, norm=True)
sfu.RotateX(90)
sfu.RotateY(180)
fvtk.add(r, sfu)
# outname = 'screenshot_signal_r.png'
# fvtk.record(r, n_frames=1, out_path = outname, size=(3000, 1500), magnification = 2)
fvtk.show(r)