crossing_bval_fideling.py

# Crossing tissues sensitivity estimation
import numpy as np

# D1, D2, D3, the tissues diffusivities
# assumes D1 >= D2=D3


# comment and uncomment as needed
# um^2 / ms
MD = 0.14
D1 = 0.25
# D2 = 0.
# compute D2 from MD and D1
D2 = (3.*MD - D1)/2.
# compute D1 from MD and D2
# D1 = 3.*MD - 2.*D2

# ms / um^2
bval = 4.

print('D1 = {}, D2=D3 = {}, b = {}\n'.format(D1, D2, bval))

# single tensor min vs max difference
S_single_min = np.exp(-bval*D1)
S_single_max = np.exp(-bval*D2)
print('Single peak signal min(at peak direction), max(perpendicular to peak)')
print('min = {}  max = {}  ratio max/min = {}, diff max-min {}\n'.format(S_single_min, S_single_max, S_single_max/S_single_min, S_single_max - S_single_min))


ds = np.arange(20,91,10.)

ratios = []
diffs = []

# crossing

for d in ds:
	S_at_peak = 0.5 * np.exp(-bval*(D1*np.cos(d*np.pi/180.)**2 + D2*(np.sin(d*np.pi/180.)**2))) + 0.5 * S_single_min

	S_at_valley = np.exp(-bval*(D1*np.cos((d/2.)*np.pi/180.)**2 + D2*(np.sin((d/2.)*np.pi/180.)**2)))

	print('Crossing signal at {} degrees'.format(d))
	print('at_peak = {}  in_valley = {}  ratio peak/valley = {}, diff peak-valley {}\n'.format(S_at_peak, S_at_valley, S_at_peak/S_at_valley, S_at_peak - S_at_valley))

	ratios.append(S_at_peak/S_at_valley)
	diffs.append(S_at_peak - S_at_valley)




import pylab as pl 

pl.figure()
pl.plot(ds, ratios)
pl.title('ratio peak/valley')

pl.figure()
pl.plot(ds, diffs)
pl.title('diff peak - valley')

pl.show()