In [3]:
import os
from dmipy.signal_models import g_from_b
from scipy.stats import gamma
gradient_folder = '/user/rfick/home/microstruktur/microstruktur/gradient_tables/'
Generate schemefile¶
In [4]:
bvals_hcp = np.loadtxt(gradient_folder + 'bvals_hcp_wu_minn.txt') # in s/mm^2, shape (N)
bvecs_hcp = np.loadtxt(gradient_folder + 'bvecs_hcp_wu_minn.txt').T # Cartesian unit vectors, shape(N x 3)
delta, Delta = np.loadtxt(gradient_folder + 'delta_Delta_hcp_wu_minn.txt', skiprows=1) # in seconds
G = g_from_b(bvals_hcp * 1e6, delta, Delta) # gradient strength in T/m. bvals input as s/m^2 so x10^6
G[G == 5.] = 0. # set b0 measurements to 0 explicitly
TE = 2 * delta + Delta + 0.001
In [10]:
N = len(bvals_hcp) # length of acquisition
scheme = np.concatenate([bvecs_hcp, G[:, None].T,
np.tile(Delta, (1, N)),
np.tile(delta, (1, N)),
np.tile(TE, (1, N))], axis=0).T # camino schemefile of size (N x 7)
header = np.array([['VERSION: STEJSKALTANNER' ,'','','','','','']]) # clumsy way to make an array with text
scheme_with_header = np.concatenate([header, scheme], axis=0) # the to-be-saved schemefile
np.savetxt(gradient_folder + 'schemefile_hcp_wu_minn.scheme1', scheme_with_header, fmt="%s")
Generate parameters for simulation substrate¶
In [77]:
lattices, alpha, beta = np.loadtxt(gradient_folder + 'camino_substrate_parameters_alexander2010.txt',
skiprows=1, usecols=(0, 1, 2)).T
x = np.linspace(0, 4e-6, 200)
for alpha_, beta_ in zip(alpha[::2], beta[::2]):
gamma_dist = gamma(alpha_, scale=beta_)
plt.plot(x * 2e6, gamma_dist.pdf(x))
xlabel('Axon Diameter ($\mu$m)', fontsize=15)
ylabel('Probability Density', fontsize=15)
title('Simulated Gamma Distributions', fontsize=17)
Out[77]:
<matplotlib.text.Text at 0x92bc110>
In [78]:
plt.hist(alpha[::2] * beta[::2] * 2e6, bins=15)
xlabel('Mean Axon Diameter ($\mu$m)', fontsize=15)
ylabel('Frequency', fontsize=15)
Out[78]:
<matplotlib.text.Text at 0x92b4910>
In [90]:
latices_onesize = lattices[::2]
latice_surfaces = latices_onesize ** 2
original_surface_intra = latice_surfaces * 0.7 # volume fraction of maximum packing
largest_latice_surface = original_surface_intra * 5 # to make the intra-axonal surface fraction 0.2
In [130]:
latice_steps = 10
output = []
fnames = []
for i in np.arange(latices_onesize.shape[0]):
latices_ = np.sqrt(np.linspace(latice_surfaces[i], largest_latice_surface[i], latice_steps))
alphas_ = np.tile(alpha[::2][i], latice_steps)
betas_ = np.tile(beta[::2][i], latice_steps)
for latnum_ in np.arange(10):
fnames.append('gamma'+str(i)+'lat'+str(latnum_)+'.Bfloat')
output.append(np.c_[latices_, alphas_, betas_])
In [132]:
parameter_list = np.c_[np.concatenate(output), np.array(fnames)]
savetxt(gradient_folder + 'camino_substrate_parameters_review_fick2017.txt', parameter_list, fmt="%s")
Camino command to start simulation¶
In [ ]: