Secure your code as it's written. Use Snyk Code to scan source code in minutes - no build needed - and fix issues immediately.
ind_mat = mapmri_isotropic_index_matrix(radial_order)
n_elem = ind_mat.shape[0]
n_rgrad = rgrad.shape[0]
K = np.zeros((n_rgrad, n_elem))
counter = 0
for n in range(0, radial_order + 1, 2):
for j in range(1, 2 + n // 2):
l = n + 2 - 2 * j
const = (-1) ** (j - 1) *\
(np.sqrt(2) * np.pi) ** (-1) *\
(r ** 2 / 2) ** (l / 2)
for m in range(-l, l+1):
K[:, counter] = const * real_sph_harm(m, l, theta, phi)
counter += 1
return K
qval, theta, phi = cart2sphere(q[:, 0], q[:, 1], q[:, 2])
theta[np.isnan(theta)] = 0
ind_mat = mapmri_isotropic_index_matrix(radial_order)
n_elem = ind_mat.shape[0]
n_qgrad = q.shape[0]
M = np.zeros((n_qgrad, n_elem))
counter = 0
for n in range(0, radial_order + 1, 2):
for j in range(1, 2 + n // 2):
l = n + 2 - 2 * j
const = mapmri_isotropic_radial_signal_basis(j, l, mu, qval)
for m in range(-l, l+1):
M[:, counter] = const * real_sph_harm(m, l, theta, phi)
counter += 1
return M
def angular_basis_EAP_opt(j, l, m, r, theta, phi):
angular_part = (
(-1) ** (j - 1) * (np.sqrt(2) * np.pi) ** (-1) *
(r ** 2 / 2) ** (l / 2) * real_sph_harm(m, l, theta, phi)
)
return angular_part
Returns
-------
B : ndarray
Matrix of the spherical harmonics model used to fit the data
m : int ``|m| <= n``
The order of the harmonic.
n : int ``>= 0``
The degree of the harmonic.
"""
if iso_comp < 2:
msg = ("Multi-tissue CSD requires at least 2 tissue compartments")
raise ValueError(msg)
r, theta, phi = geo.cart2sphere(*gtab.gradients.T)
m, n = shm.sph_harm_ind_list(sh_order)
B = shm.real_sph_harm(m, n, theta[:, None], phi[:, None])
B[np.ix_(gtab.b0s_mask, n > 0)] = 0.
iso = np.empty([B.shape[0], iso_comp])
iso[:] = SH_CONST
B = np.concatenate([iso, B], axis=1)
return B, m, n
def sim_response(sh_order=8, bvals=bvals, evals=evals_d, csf_md=csf_md,
gm_md=cgm_md):
bvals = np.array(bvals, copy=True)
evecs = np.zeros((3, 3))
z = np.array([0, 0, 1.])
evecs[:, 0] = z
evecs[:2, 1:] = np.eye(2)
n = np.arange(0, sh_order + 1, 2)
m = np.zeros_like(n)
big_sphere = default_sphere.subdivide()
theta, phi = big_sphere.theta, big_sphere.phi
B = shm.real_sph_harm(m, n, theta[:, None], phi[:, None])
A = shm.real_sph_harm(0, 0, 0, 0)
response = np.empty([len(bvals), len(n) + 2])
for i, bvalue in enumerate(bvals):
gtab = GradientTable(big_sphere.vertices * bvalue)
wm_response = single_tensor(gtab, 1., evals, evecs, snr=None)
response[i, 2:] = np.linalg.lstsq(B, wm_response)[0]
response[i, 0] = np.exp(-bvalue * csf_md) / A
response[i, 1] = np.exp(-bvalue * gm_md) / A
return MultiShellResponse(response, sh_order, bvals)
def multi_tissue_basis(gtab, sh_order, iso_comp):
"""Builds a basis for multi-shell CSD model"""
if iso_comp < 1:
msg = ("Multi-tissue CSD requires at least 2 tissue compartments")
raise ValueError(msg)
r, theta, phi = geo.cart2sphere(*gtab.gradients.T)
m, n = shm.sph_harm_ind_list(sh_order)
B = shm.real_sph_harm(m, n, theta[:, None], phi[:, None])
B[np.ix_(gtab.b0s_mask, n > 0)] = 0.
iso = np.empty([B.shape[0], iso_comp])
iso[:] = sh_const
B = np.concatenate([iso, B], axis=1)
return B, m, n
def multi_tissue_basis(gtab, sh_order, iso_comp):
"""Builds a basis for multi-shell CSD model"""
if iso_comp < 1:
msg = ("Multi-tissue CSD requires at least 2 tissue compartments")
raise ValueError(msg)
r, theta, phi = geo.cart2sphere(*gtab.gradients.T)
m, n = shm.sph_harm_ind_list(sh_order)
B = shm.real_sph_harm(m, n, theta[:, None], phi[:, None])
B[np.ix_(gtab.b0s_mask, n > 0)] = 0.
iso = np.empty([B.shape[0], iso_comp])
iso[:] = sh_const
B = np.concatenate([iso, B], axis=1)
return B, m, n
def basis(self, sphere):
"""A basis that maps the response coefficients onto a sphere."""
theta = sphere.theta[:, None]
phi = sphere.phi[:, None]
return real_sph_harm(self.m, self.n, theta, phi)
scale factor
sphere_vertices : array, shape (N,3)
vertices of the odf sphere
"""
r, theta, phi = cart2sphere(sphere_vertices[:, 0], sphere_vertices[:, 1], sphere_vertices[:, 2])
theta[np.isnan(theta)] = 0
counter = 0
upsilon = np.zeros(
(len(sphere_vertices), (radialOrder + 1) * ((radialOrder + 1) / 2) * (2 * radialOrder + 1)))
for n in range(radialOrder + 1):
for l in range(0, n + 1, 2):
for m in range(-l, l + 1):
upsilon[:, counter] = (-1) ** (n - l / 2.0) * __kappa_odf(zeta, n, l) * \
hyp2f1(l - n, l / 2.0 + 1.5, l + 1.5, 2.0) * \
real_sph_harm(m, l, theta, phi)
counter += 1
return upsilon[:, 0:counter]