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

DKI PR3 - NF: Adding standard kurtosis statistics on module dki.py #677

Merged
merged 52 commits into from Oct 5, 2015

Conversation

Projects
None yet
5 participants
@RafaelNH
Contributor

RafaelNH commented Jul 13, 2015

This works follows the work done in PR #664.

Now that the functions to estimate the kurtosis tensor are implemented, we can proceed with the estimate of standard kurtosis statistics (or also known as the standard DKI rotational invariant measures). In this PR, the analytical solutions of the following statistics are implemented:

  • Mean kurtosis
  • Radial kurtosis
  • Axial kurtosis

For the mean kurtosis analytical solution, the implementation of the Carlson's elliptical integrals was required. Comparative to the code that Prof. Maurizio Marrale kindly supplied, I followed a more recent, general and optimized implementation of the Carlson integrals which is described here: http://arxiv.org/pdf/math/9409227v1.pdf . Despite this, estimation of the integrals are still time consuming, so I am planing to convert functions carlson_rf and carlson_rd in cython.

In alternative to the analytical solution, two numerical methods for a faster estimation of the standard DKI measures are implemented. These numerical methods are expected to be less accurate that the analytical solutions, however they can provide alternatives less computationally demanding and provide a simpler mathematical framework which will be used to further validate the analytical solutions.

@RafaelNH

This comment has been minimized.

Contributor

RafaelNH commented Jul 13, 2015

PS: At the moment this pull request only have the implementations of the Carlson's integrals, however in the following hours I will add the missing code.

@RafaelNH

This comment has been minimized.

Contributor

RafaelNH commented Jul 14, 2015

A final version of the codes for the mean, radial and axial kurtosis were now added!

These are the planed future steps:

  • write a the numerical method for the radial kurtosis. This will depend on #674. This method is not an exact analytical solution, however it is useful since it does not require rotating the kurtosis tensors to a frame of reference that diagonalizes the diffusion tensor, and perpendicular kurtosis can be computed relative to other vectors besides the diffusion tensor principal axis.
  • replace the Sherbrooke's 3 shells data on example by an HCP-like dataset. Sherbrooke's data seems to have insufficient SNR for fitting the diffusion kurtosis model. For more detail on this see the discussion in the NeuroImage mailing list.
  • write carlson_rf and carlson_rd in cython (????).

Other point to discuss:

  • When diffusion parameters are zero, kurtosis values are defined as NaN. As pointed by @arokem, having this NaN can be problematic for image analyzing steps. Should I replace this NaN by zeros?
@RafaelNH

This comment has been minimized.

Contributor

RafaelNH commented Jul 17, 2015

Any comments so far? Ignore the last point of my last github comment. From 399cce7, dki module does not produce NaN anymore.

ax.flat[2].set_title('AK')
plt.show()
fig3.savefig('Kurtosis_tensor_standard_measures.png')

This comment has been minimized.

@arokem

arokem Jul 17, 2015

Member

You need to add the .. figure directive in the text block below for this figure to be included in the webpage

sph = Sphere(xyz=gtab.bvecs[gtab.bvals > 0])
MK_nm = mean_kurtosis(dkiF.model_params, sph)
assert_array_almost_equal(MK_as, MK_nm, decimal=1)

This comment has been minimized.

@arokem

arokem Jul 17, 2015

Member

Why would these two computations be so different (decimal=1 !)?

They seem to take the same code path, calling dipy.reconst.dki.mean_kurtosis on the same set of parameters here, and in line https://github.com/RafaelNH/dipy/blob/DKI_PR3_kurtosis_statistics/dipy/reconst/dki.py#L1257

This comment has been minimized.

@RafaelNH

RafaelNH Jul 18, 2015

Contributor

When a sphere is given as input, MK is estimated as the average of the directional kurtosis of the sphere's vertices instead of the MK's analytical solution (see documentation dipy.reconst.dki.mean_kurtosis). This is why computations can be so different.

An = (xn + yn + zn) / 3.0
Q = (3.*errtol) ** (-1/6.) * np.max([np.abs(An - xn), np.abs(An - yn),
np.abs(An - zn)])
# Convergence condition

This comment has been minimized.

@arokem

arokem Jul 17, 2015

Member

Looks to me like you can do all the computations up to this point on the entire vector at once. So no need to have them inside the for loop.

An = (An + lamda)*0.250
n = n + 1
# post convergence calculation

This comment has been minimized.

@arokem

arokem Jul 17, 2015

Member

Same goes for all the computations below this point

Q = (errtol/4.) ** (-1/6.) * np.max([np.abs(An - xn), np.abs(An - yn),
np.abs(An - zn)])
sum_term = 0
# Convergence condition

This comment has been minimized.

@arokem

arokem Jul 17, 2015

Member

I think that the same is true in this function

Parameters
----------
L1 : ndarray (n,)

This comment has been minimized.

@arokem

arokem Jul 17, 2015

Member

I think this would work on ndarrays regardless of their shape. This seems to imply that it would only work with 1D arrays.

zero.
"""
if er is None:
er = np.finfo(L1[0]).eps * 1e3

This comment has been minimized.

@arokem

arokem Jul 17, 2015

Member

Except for this line. You would have write this as:

np.finfo(L1.ravel()[0]).eps * 1e3

This comment has been minimized.

@arokem

arokem Jul 17, 2015

Member

Once you do that. I believe that all the other functions below could also be written to operate on arbitrary-shape ndarrays. In particular, check out from dipy.core.ndindex import ndindex to iterate over arbitray-shape arrays

Q = (3.*errtol) ** (-1/6.) * np.max([np.abs(An - xn), np.abs(An - yn),
np.abs(An - zn)], axis=0)
# Convergence has to be done voxel by voxel
for v in range(len(x)):

This comment has been minimized.

@arokem

arokem Jul 18, 2015

Member

Could you use ndindex here, so as to make these functions work on arbitrary-dimensional arrays?

This comment has been minimized.

@arokem

arokem Jul 20, 2015

Member

I believe this should work here (above, line 102):

np.max(np.abs(An - np.concatenate([xn, yn, zn], axis=-1)), -1)
n = np.zeros(len(x))
# Convergence has to be done voxel by voxel
for v in range(len(x)):

This comment has been minimized.

@arokem

arokem Jul 18, 2015

Member

Here as well

@arokem

This comment has been minimized.

Member

arokem commented Jul 18, 2015

Have you had a chance to profile this code? Might be worth seeing what the performance bottle-necks are, before considering cythonizing some of it (e.g. the internal loop in the Carlson integrals). Some notes on profiling here: https://github.com/arokem/learning-python/tree/master/profiling

@RafaelNH RafaelNH referenced this pull request Jul 20, 2015

Closed

(WIP) DKI PR5 - NF: DKI-ODF estimation #685

3 of 3 tasks complete
Wzzzz[vox] = Wrotate(kt[vox], evecs[vox], [2, 2, 2, 2])
Wyyzz[vox] = Wrotate(kt[vox], evecs[vox], [1, 1, 2, 2])
# Compute MK

This comment has been minimized.

@arokem

arokem Jul 20, 2015

Member

You mean "compute RK?"

Returns
-------
W : array(4,4,4,4)
Full 4D kutosis tensor

This comment has been minimized.

@arokem

arokem Jul 20, 2015

Member

typo: "kutosis" => "kurtosis"

while 4.**(-n[v]) * Q[v] > abs(An[v]):
xnroot = np.sqrt(xn[v])
ynroot = np.sqrt(yn[v])
znroot = np.sqrt(zn[v])

This comment has been minimized.

@arokem

arokem Jul 20, 2015

Member

math.sqrt might be faster

@RafaelNH RafaelNH force-pushed the RafaelNH:DKI_PR3_kurtosis_statistics branch 2 times, most recently from 67d611f to b8f683d Jul 22, 2015

@RafaelNH

This comment has been minimized.

Contributor

RafaelNH commented Jul 22, 2015

A review of what I have done on the last couple of days:

  • According to a suggestion made by @arokem, Carlson integrals are now not limited to 'flat vectors'.
  • @arokem, I also try to replace the np.sqrt inside the Carlson integral by the math.sqrt, however at the end, I decided to keep np.sqrt because it automatically deals better with complex numbers. I could had use cmath instead, however when working with real numbers the output variables from cmath will have a 0j imaginary part.
  • By profiling the mean_kurtosis function, I figured that Carlson integrals are not the unique computational demanding steps in MK estimation. Another computational demanding step is the kurtosis rotation done by Wrotate. I tried to use the kurtosis tensor symmetry to improve Wrotate performance by avoiding the computational demanding loops, however after analysis better the code, I noticed that this is not that trivial. The calculations inside the loops also depend on the non-symmetric tensor R. Instead, I factorize Wrotate by removing its dependency to Wcons. In addition, some unnecessary steps on MK and RK functions were removed by calling directly the helpers function _Wrotate_element. Re-processing the data of the mid-term summary, procedures are now 1.6 times faster (MK takes around 14 mins to compute, RK 7 mins and AK 1min).
  • I add two test to better evaluate the statistical kurtosis measures based on less mathematical demanding approaches. Tests are covering all the lines of dki.py as you can see from nosetest reports.

Please let me know if you have suggestions for the next steps of this PR. Do you think that we could improve performance at this point by cythonizing Wrotate and Carlson integrals? Rather than this, do you have any other comments that I could address?

RafaelNH added a commit to RafaelNH/dipy that referenced this pull request Jul 24, 2015

The test for evaluating the dki-odf estimate depends on function dki.…
…Wcons which is being revised in PR nipy#677. I am adding here this function to carry with thr dki-odf test. However when the 2 pull requests are merged we have to insure that duplicates of this function are removed
MK = MK.reshape(outshape)

This comment has been minimized.

@arokem

arokem Aug 3, 2015

Member

Please remove these lines - too much whitespace!

@RafaelNH RafaelNH force-pushed the RafaelNH:DKI_PR3_kurtosis_statistics branch from 4ef1e33 to a433bdb Aug 5, 2015

@RafaelNH

This comment has been minimized.

Contributor

RafaelNH commented Aug 6, 2015

I think all the work here is ready to be merged. I already finish all points that were remaining on this PR. In particular:

  • I replace the Sherbrooke's 3 shells data on the example by the HCP-like dataset added in @arokem 's PR #691. Considering that we are not using any noise removal approach, kurtosis metrics of the example are now looking nice.
  • I further refactored the kurtosis statistics metrics and now procedures take less that 1min to process the data of the mid-term summary. I also tried to accelerate more the codes by cythonizing the Carlson integrals, however, at the end, I decided to keep the codes in python since the gains weren't anymore significant.

Please let me know if you have any other comment to address.

RafaelNH added a commit to RafaelNH/dipy that referenced this pull request Aug 10, 2015

The test for evaluating the dki-odf estimate depends on function dki.…
…Wcons which is being revised in PR nipy#677. I am adding here this function to carry with thr dki-odf test. However when the 2 pull requests are merged we have to insure that duplicates of this function are removed
biologic tissues is non-Gaussian using the kurtosis tensor (KT) [Jensen2005]_.
Measurements of non-Gaussian diffusion from the diffusion kurtosis model are of
interest because they can be used to charaterize tissue microstructural
heterogeneity [Jensen2010]_ and to derive concrete biophysical parameters as

This comment has been minimized.

@arokem

arokem Aug 10, 2015

Member

=> "...concrete biophysical parameters, such as..."

signal in the absence of diffusion gradient sensitisation, $S_0$, and the
values of diffusion, $\mathbf{D(n)}$, and diffusion kurtosis, $\mathbf{K(n)}$,
along the spatial direction $\mathbf{n}$ [NetoHe2015]_:
The diffusion kurtosis model express the diffusion-weighted signal as:

This comment has been minimized.

@arokem

arokem Aug 10, 2015

Member

=> "expresses"

$\mathbf{D(n)}$ and $\mathbf{K(n)}$ can be computed from the DT and KT using
the following equations:
where $\mathbf{b}$ is the applied diffusion weighting (which is dependent on
the measuremenent parameters measurement), $S_0$ the signal in the absence of

This comment has been minimized.

@arokem

arokem Aug 10, 2015

Member

the word "measurement" appears twice here.

the following equations:
where $\mathbf{b}$ is the applied diffusion weighting (which is dependent on
the measuremenent parameters measurement), $S_0$ the signal in the absence of
diffusion gradient sensitisation, $\mathbf{D(n)}$ the value of diffusion along

This comment has been minimized.

@arokem

arokem Aug 10, 2015

Member

=> "sensitization"

RafaelNH added some commits Aug 13, 2015

Tests: add tests to cover cases that functions gives ValueError - 1) …
…when unknow fit is given; 2) when min_signal is smaller than 0; 3) when given mask as a different shape to data
Tests: add test_split_dki_params to check that outputs of function sp…
…lit_dki_params match the parameters extracted by dkifit.kt, dkifit.evals, and dkifit.evecs
TESTS, BF: Add lines of code in test_dki_predict to also test the sig…
…nal predictions of the DiffusionKurtosisFit object. From this tests, an error on this prediction procedure was found and corrected
RF: Redefine the creteria of eigenvalue comparison to correctproblema…
…tic performance in regions near to the MK analytical solution. For more information see post #13 of http://gsoc2015dipydki.blogspot.co.uk

@RafaelNH RafaelNH force-pushed the RafaelNH:DKI_PR3_kurtosis_statistics branch from 76c81fe to d2bf028 Sep 25, 2015

@RafaelNH

This comment has been minimized.

Contributor

RafaelNH commented Sep 25, 2015

@arokem - I already rebase this PR, and yup the new version of function nlmean cannot accept anymore an array of sigmas for each diffusion-weighted volume!

@samuelstjean - I give a look on the code of the functions nlmean. So basically the new version of nlmeans can accept two types of sigma: 1) single value; 2) sigma profile. Regardless the type of sigma, nlmean will act the same according to the format of the data. For 3D data nlmean calls nlmean_3D once, while for 4D data nlmean calls nlmean_3D for each volume.

For my understanding there isn't no function in dipy that generates this sigma profiles, right? Even piesno will not be compatible with the new version of nlmean since this produces an array of sigmas per slice and not an 3D noise profile.

@samuelstjean

This comment has been minimized.

Contributor

samuelstjean commented Sep 25, 2015

That's correct, so far no publicly, easily available function will do 3D noise estimation. You can still use piesno and broadcast the 2D result to a 3D array though. The default behavior should still be (and should also always be) to call the nlmeans3D version on each volume separately.

According to your first comment, does something not work as expected with the new format of the nlmeans function?

@RafaelNH

This comment has been minimized.

Contributor

RafaelNH commented Sep 25, 2015

Yes. Previously if an array of sigmas per volume was given (output of estimate_sigma), nlmean was able to process each volume with the respective sigma on the array.

Indeed, we still can use piesno and estimate_sigma by converting their outputs to 3D arrays. My only concern is if these steps are also trivial for new Dipy users or users that are new to python. I think we should make nlmean compatible with both piesno and estimate_sigma outputs without changing the new version of nlmeans3D.

@samuelstjean

This comment has been minimized.

Contributor

samuelstjean commented Sep 25, 2015

That's weird, it should work as-is with the output of estimate_sigma if you call nlmeans directly (or I copy pasted the wrong function, but the tests are passing...), as it wraps nlmeans_3D in the proper way if I did not make any mistake. Afaik, nlmeans3D should not be really called directly unless you know how it works, as the function simply named nlmeans should set everything for you.

@arokem

This comment has been minimized.

Member

arokem commented Sep 25, 2015

To be concrete, the test implemented in #721 should pass?

@arokem

This comment has been minimized.

Member

arokem commented Sep 25, 2015

I also propose a BF in #721 - should be a quick fix that will make this merge finally possible

@RafaelNH

This comment has been minimized.

Contributor

RafaelNH commented Sep 26, 2015

Hi @arokem! Yup, #721 seems to do the magic!

arokem added a commit that referenced this pull request Oct 5, 2015

Merge pull request #677 from RafaelNH/DKI_PR3_kurtosis_statistics
DKI PR3 - NF: Adding standard kurtosis statistics on module dki.py

@arokem arokem merged commit 7fb36d9 into nipy:master Oct 5, 2015

1 check passed

continuous-integration/travis-ci/pr The Travis CI build passed
Details
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment