ENH, WIP: Compare CSD on multishell data to DTI on downsampled data #3157
Conversation
|
Hello @alexrockhill, Thank you for submitting the Pull Request !
Do see the DIPY coding Style guideline |
There was a problem hiding this comment.
Hi @alexrockhill,
Thank you for starting this. See below my comments.
Furthermore, the file need to be added in doc/examples/_valid_examples.toml.
I suppose the tutorial is not finished because I do not see any comparison.
EDIT: ok I just saw the WIP + draft, so it make sense.
| bvals = gtab.bvals | ||
| bvecs = gtab.bvecs |
There was a problem hiding this comment.
I do not think you need this, it is deffine above
|
|
||
| ############################################################################### | ||
| # For the sake of simplicity, we only select two non-zero b-values for this | ||
| # example. |
There was a problem hiding this comment.
Can you add the number of b-value we have on this data and which one we select. it might not be obvious for some people when reading the code
| sel_b = np.logical_or(np.logical_or(bvals == 0, bvals == 1000), bvals == 2000) | ||
| data = data[..., sel_b] | ||
|
|
||
| gtab = gradient_table(bvals[sel_b], bvecs[sel_b]) |
There was a problem hiding this comment.
Note to myself: we might need an explicity function for that in DIPY
| plt.imshow(np.rot90(maskdata[:, :, 9, 0]), cmap=plt.cm.bone) | ||
| plt.show() | ||
|
|
||
| denoised = mppca(data, mask=mask, patch_radius=2) | ||
|
|
||
| axial_slice = 10 | ||
| plt.imshow(np.rot90(denoised[:, :, axial_slice, 0]), cmap=plt.cm.bone) | ||
| plt.imshow(np.rot90(data[:, :, axial_slice, 0]), cmap=plt.cm.hot, | ||
| alpha=0.25) | ||
| plt.show() |
There was a problem hiding this comment.
can you check if we do not have already a function in dipy.viz.plottingfor that. if yes, use it, if no, create it and just call it in the tutorial
| plt.imshow(np.rot90(ap[:, :, axial_slice]), cmap=plt.cm.bone) | ||
| plt.show() |
There was a problem hiding this comment.
same as above and explain what we see
| hmrf = TissueClassifierHMRF() | ||
| initial_segmentation, final_segmentation, PVE = hmrf.classify( | ||
| ap, nclasses=3, beta=0.1) | ||
|
|
||
| csf = np.where(final_segmentation == 1, 1, 0) | ||
| gm = np.where(final_segmentation == 2, 1, 0) | ||
| wm = np.where(final_segmentation == 3, 1, 0) | ||
|
|
||
| plt.imshow(np.rot90(denoised[:, :, axial_slice, 0]), cmap=plt.cm.bone) | ||
| plt.imshow(np.rot90(np.where(gm[:, :, axial_slice], 0.5, np.nan)), | ||
| cmap=plt.cm.bone_r, vmin=0, vmax=1, alpha=0.5) | ||
| plt.imshow(np.rot90(np.where(wm[:, :, axial_slice], 0.9, np.nan)), | ||
| cmap=plt.cm.bone, vmin=0, vmax=1, alpha=0.5) | ||
| plt.imshow(np.rot90(np.where(csf[:, :, axial_slice], 1, np.nan)), | ||
| cmap=plt.cm.cool, alpha=0.5) | ||
| plt.show() |
There was a problem hiding this comment.
same as above and explain what we see
|
I think I might just add this to the end of the existing tutorial instead |
|
I'm having a bit of issues with this deciding what's the best way for comparison, I feel like just visualizing the corpus callosum like the other tutorials isn't really that telling about how much better having multishell data is. Is there phantom data to use to compare? Just thinking out loud and sharing in case anyone had advice, I haven't searched the tutorials but I will soon! |
Fixes #3123