Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

typos + pep stuf #308

Merged
merged 2 commits into from Jan 7, 2014
Merged

typos + pep stuf #308

Changes from 1 commit
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
749b609
typos + pep stuf
samuelstjean Jan 6, 2014
953d5f2
Reverted some change
samuelstjean Jan 7, 2014
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
60 changes: 31 additions & 29 deletions dipy/reconst/csdeconv.py
View file
Open in desktop
Expand Up @@ -15,6 +15,7 @@
from dipy.reconst.dti import TensorModel, fractional_anisotropy
from scipy.integrate import quad


class ConstrainedSphericalDeconvModel(OdfModel, Cache):

def __init__(self, gtab, response, reg_sphere=None, sh_order=8, lambda_=1, tau=0.1):
Expand Down Expand Up @@ -73,7 +74,7 @@ def __init__(self, gtab, response, reg_sphere=None, sh_order=8, lambda_=1, tau=0

no_params = ((sh_order + 1) * (sh_order + 2)) / 2

if no_params > np.sum(gtab.b0s_mask == False):
if no_params > np.sum(gtab.b0s_mask is False):
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

OK - this won't work. The comparison is an array comparison, so I don't think that we want to use is.

msg = "Number of parameters required for the fit are more "
msg += "than the actual data points"
warnings.warn(msg, UserWarning)
Expand Down Expand Up @@ -164,7 +165,7 @@ def __init__(self, gtab, ratio, reg_sphere=None, sh_order=8, lambda_=1., tau=0.1

no_params = ((sh_order + 1) * (sh_order + 2)) / 2

if no_params > np.sum(gtab.b0s_mask == False):
if no_params > np.sum(gtab.b0s_mask is False):
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Same here

msg = "Number of parameters required for the fit are more "
msg += "than the actual data points"
warnings.warn(msg, UserWarning)
Expand Down Expand Up @@ -338,17 +339,17 @@ def forward_sdt_deconv_mat(ratio, n, r2_term=False):
Parameters
----------
ratio : float
ratio = $\frac{\lambda_2}{\lambda_1}$ of the single fiber response
ratio = $\frac{\lambda_2}{\lambda_1}$ of the single fiber response
function
n : ndarray (N,)
The degree of spherical harmonic function associated with each row of
the deconvolution matrix. Only even degrees are allowed.
r2_term : bool
True if ODF comes from an ODF computed from a model using the $r^2$ term
in the integral. For example, DSI, GQI, SHORE, CSA, Tensor, Multi-tensor
ODFs. This results in using the proper analytical response function
solution solving from the single-fiber ODF with the r^2 term. This
derivation is not published anywhere but is very similar to [1]_.
ODFs. This results in using the proper analytical response function
solution solving from the single-fiber ODF with the r^2 term. This
derivation is not published anywhere but is very similar to [1]_.

Returns
-------
Expand All @@ -369,11 +370,11 @@ def forward_sdt_deconv_mat(ratio, n, r2_term=False):
frt = np.zeros(n_degrees) # FRT (Funk-Radon transform) q-ball matrix

for l in np.arange(0, n_degrees*2, 2):
if r2_term :
if r2_term :
sharp = quad(lambda z: lpn(l, z)[0][-1] * gamma(1.5) * np.sqrt( ratio / (4 * np.pi ** 3) ) /
np.power((1 - (1 - ratio) * z ** 2), 1.5), -1., 1.)
else :
sharp = quad(lambda z: lpn(l, z)[0][-1] * np.sqrt(1 / (1 - (1 - ratio) * z * z)), -1., 1.)
sharp = quad(lambda z: lpn(l, z)[0][-1] * np.sqrt(1 / (1 - (1 - ratio) * z * z)), -1., 1.)

sdt[l / 2] = sharp[0]
frt[l / 2] = 2 * np.pi * lpn(l, 0)[0][-1]
Expand Down Expand Up @@ -414,9 +415,9 @@ def csdeconv(s_sh, sh_order, R, B_reg, lambda_=1., tau=0.1):
Returns
-------
fodf_sh : ndarray (``(sh_order + 1)*(sh_order + 2)/2``,)
Spherical harmonics coefficients of the constrained-regarized fiber ODF
Spherical harmonics coefficients of the constrained-regularized fiber ODF
num_it : int
Number of iterations in the constrained-regarization used for convergence
Number of iterations in the constrained-regularization used for convergence

References
----------
Expand All @@ -426,7 +427,7 @@ def csdeconv(s_sh, sh_order, R, B_reg, lambda_=1., tau=0.1):
"""

# generate initial fODF estimate, truncated at SH order 4
fodf_sh = np.linalg.lstsq(R, s_sh)[0]
fodf_sh = np.linalg.lstsq(R, s_sh)[0]
fodf_sh[15:] = 0

fodf = np.dot(B_reg, fodf_sh)
Expand All @@ -452,25 +453,26 @@ def csdeconv(s_sh, sh_order, R, B_reg, lambda_=1., tau=0.1):
k = k2

# This is the super-resolved trick.
# Wherever there is a negative amplitude value on the fODF, it
# concatinates a value to the S vector so that the estimation can
# focus on trying to eliminate it. In a sense, this "adds" a
# measurement, which can help to better estimate the fodf_sh, even if
# Wherever there is a negative amplitude value on the fODF, it
# concatinates a value to the S vector so that the estimation can
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

"concatenates"

# focus on trying to eliminate it. In a sense, this "adds" a
# measurement, which can help to better estimate the fodf_sh, even if
# you have more SH coeffcients to estimate than actual S measurements.
M = np.concatenate((R, lambda_ * B_reg[k, :]))
S = np.concatenate((s_sh, np.zeros(k.shape)))
try:
fodf_sh = np.linalg.lstsq(M, S)[0]
except np.linalg.LinAlgError as lae:
# SVD did not converge in Linear Least Squares in current
# SVD did not converge in Linear Least Squares in current
# voxel. Proceeding with initial SH estimate for this voxel.
pass

warnings.warn('maximum number of iterations exceeded - failed to converge')
return fodf_sh, num_it


def odf_deconv(odf_sh, R, B_reg, lambda_=1., tau=0.1, r2_term=False):
r""" ODF constrained-regularized sherical deconvolution using
r""" ODF constrained-regularized spherical deconvolution using
the Sharpening Deconvolution Transform (SDT) [1]_, [2]_.

Parameters
Expand All @@ -490,12 +492,12 @@ def odf_deconv(odf_sh, R, B_reg, lambda_=1., tau=0.1, r2_term=False):
True if ODF is computed from model that uses the $r^2$ term in the integral.
Recall that Tuch's ODF (used in Q-ball Imaging [1]_) and the true normalized ODF
definition differ from a $r^2$ term in the ODF integral. The original Sharpening
Deconvolution Transform (SDT) technique [2]_ is expecting Tuch's ODF without
the $r^2$ (see [3]_ for the mathematical details).
Deconvolution Transform (SDT) technique [2]_ is expecting Tuch's ODF without
the $r^2$ (see [3]_ for the mathematical details).
Now, this function supports ODF that have been computed using the $r^2$ term because
the proper analytical response function has be derived.
the proper analytical response function has be derived.
For example, models such as DSI, GQI, SHORE, CSA, Tensor, Multi-tensor ODFs, should now
be deconvolved with the r2_term=True.
be deconvolved with the r2_term=True.

Returns
-------
Expand All @@ -519,7 +521,7 @@ def odf_deconv(odf_sh, R, B_reg, lambda_=1., tau=0.1, r2_term=False):

# if sharpening a q-ball odf (it is NOT properly normalized), we need to force normalization
# otherwise, for DSI, CSA, SHORE, Tensor odfs, they are normalized by construction
if ~r2_term :
if ~r2_term :
Z = np.linalg.norm(fodf)
fodf_sh /= Z

Expand Down Expand Up @@ -547,7 +549,7 @@ def odf_deconv(odf_sh, R, B_reg, lambda_=1., tau=0.1, r2_term=False):
try:
fodf_sh = np.linalg.lstsq(M, ODF)[0]
except np.linalg.LinAlgError as lae:
# SVD did not converge in Linear Least Squares in current
# SVD did not converge in Linear Least Squares in current
# voxel. Proceeding with initial SH estimate for this voxel.
pass

Expand Down Expand Up @@ -585,12 +587,12 @@ def odf_sh_to_sharp(odfs_sh, sphere, basis=None, ratio=3 / 15., sh_order=8, lamb
True if ODF is computed from model that uses the $r^2$ term in the integral.
Recall that Tuch's ODF (used in Q-ball Imaging [1]_) and the true normalized ODF
definition differ from a $r^2$ term in the ODF integral. The original Sharpening
Deconvolution Transform (SDT) technique [2]_ is expecting Tuch's ODF without
the $r^2$ (see [3]_ for the mathematical details).
Deconvolution Transform (SDT) technique [2]_ is expecting Tuch's ODF without
the $r^2$ (see [3]_ for the mathematical details).
Now, this function supports ODF that have been computed using the $r^2$ term because
the proper analytical response function has be derived.
the proper analytical response function has be derived.
For example, models such as DSI, GQI, SHORE, CSA, Tensor, Multi-tensor ODFs, should now
be deconvolved with the r2_term=True.
be deconvolved with the r2_term=True.

Returns
-------
Expand Down Expand Up @@ -618,7 +620,7 @@ def odf_sh_to_sharp(odfs_sh, sphere, basis=None, ratio=3 / 15., sh_order=8, lamb
fodf_sh = np.zeros(odfs_sh.shape)

for index in ndindex(odfs_sh.shape[:-1]):
fodf_sh[index], num_it = odf_deconv(odfs_sh[index], R, B_reg, lambda_=lambda_,
fodf_sh[index], num_it = odf_deconv(odfs_sh[index], R, B_reg, lambda_=lambda_,
tau=tau, r2_term=r2_term)

return fodf_sh
Expand All @@ -633,7 +635,7 @@ def auto_response(gtab, data, roi_center=None, roi_radius=10, fa_thr=0.7):
data : ndarray
diffusion data
roi_center : tuple, (3,)
Center of ROI in data. If center is None, it is assumed that is
Center of ROI in data. If center is None, it is assumed that it is
the center of the volume with shape `data.shape[:3]`.
roi_radius : int
radius of cubic ROI
Expand Down