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

Gibbs removal #1771

Merged
merged 33 commits into from Jul 16, 2019

Conversation

@RafaelNH
Copy link
Contributor

commented Mar 14, 2019

In this pull request, I am submitting a new pre-processing algorithm to remove Gibbs Ringing artefacts.

This algorithm, is an adaption of the algorithm suggested by:

Kellner et al. Gibbs‐ringing artifact removal based on local subvoxel‐shifts. MRM 2016

The full adapted version of algorithm is described in chapter 3 of my phd thesis:

Neto Henriques, R. (2018). Advanced Methods for Diffusion MRI Data Analysis and their Application to the Healthy Ageing Brain (Doctoral thesis). https://doi.org/10.17863/CAM.29356

@pep8speaks

This comment has been minimized.

Copy link

commented Mar 14, 2019

Hello @RafaelNH, Thank you for updating !

Line 30:1: E402 module level import not at top of file
Line 40:1: E402 module level import not at top of file
Line 41:1: E402 module level import not at top of file
Line 144:1: E402 module level import not at top of file
Line 204:1: E402 module level import not at top of file
Line 210:1: E402 module level import not at top of file

Comment last updated at 2019-07-12 11:28:41 UTC
@arokem arokem referenced this pull request Mar 14, 2019

@RafaelNH RafaelNH changed the title WIP - Gibbs removal Gibbs removal Mar 18, 2019

@RafaelNH RafaelNH marked this pull request as ready for review Mar 18, 2019

@RafaelNH RafaelNH requested review from skoudoro and arokem Mar 18, 2019

@arokem
Copy link
Member

left a comment

Overall looks good. I had some questions and comments. In particular, I wonder about the naming of the different functions. Do you expect users to access any of the functions except gibb_removal? Maybe we make that one the only one that isn't private to the module? Or is that too much?

Gibbs suppression
~~~~~~~~~~~~~~~~~

- :ref:`example_suppress_gibbs`

This comment has been minimized.

Copy link
@arokem

arokem Mar 19, 2019

Member

Maybe put this under the "denoising" header below?

dipy/denoise/gibbs.py Outdated Show resolved Hide resolved
dipy/denoise/gibbs.py Outdated Show resolved Hide resolved
dipy/denoise/gibbs.py Outdated Show resolved Hide resolved
dipy/denoise/gibbs.py Outdated Show resolved Hide resolved
return imagec


def gibbs_removal(vol, slice_axis=2, nn=3):

This comment has been minimized.

Copy link
@arokem

arokem Mar 19, 2019

Member

So this will be the main public interface for this module? This is the function that most users will want to access? I find the naming scheme here a little bit confusing. Maybe the others (the "sub-functions") should be designated as private functions by adding an underscore prefix?

doc/examples/reconst_dki.py Show resolved Hide resolved
[1]_, [2]_.
In the following example, we show how to suppress Gibbs artefacts of MRI images
in dipy. This algorithm is based on an adapted version of the sub-voxel

This comment has been minimized.

Copy link
@arokem

arokem Mar 19, 2019

Member

You don't need "in dipy" here.

In the following example, we show how to suppress Gibbs artefacts of MRI images
in dipy. This algorithm is based on an adapted version of the sub-voxel
Gibbs suppression procedure [3]_. Full details of the implemented algorithm
can be found in the Chapter 3 of [4]_ (please cite [3]_, [4]_ if you are using

This comment has been minimized.

Copy link
@arokem

arokem Mar 19, 2019

Member

"the Chapter" => "chapter"

@@ -0,0 +1,155 @@
# -*- coding: utf-8 -*-

This comment has been minimized.

Copy link
@arokem

arokem Mar 19, 2019

Member

Maybe rename this file "denoise_gibbs.py"?

@codecov-io

This comment has been minimized.

Copy link

commented Mar 20, 2019

Codecov Report

Merging #1771 into master will increase coverage by 0.09%.
The diff coverage is 100%.

Impacted file tree graph

@@            Coverage Diff             @@
##           master    #1771      +/-   ##
==========================================
+ Coverage   85.31%   85.41%   +0.09%     
==========================================
  Files         118      119       +1     
  Lines       14218    14312      +94     
  Branches     2236     2249      +13     
==========================================
+ Hits        12130    12224      +94     
  Misses       1577     1577              
  Partials      511      511
Impacted Files Coverage Δ
dipy/denoise/gibbs.py 100% <100%> (ø)
dipy/reconst/mapmri.py 90.33% <0%> (ø) ⬆️
dipy/sims/voxel.py 90.58% <0%> (+0.12%) ⬆️

@RafaelNH RafaelNH force-pushed the RafaelNH:Gibbs_Removal branch from ec19d98 to c51e49c Mar 21, 2019

@RafaelNH

This comment has been minimized.

Copy link
Contributor Author

commented Mar 21, 2019

@arokem! Thanks for your review. I am still not sure if this should go under the denoising heather because I don't see Gibbs suppression as a denoising algorithm. But, anyway, I also agree that perhaps is not efficient to create a new heather just for this algorithm. Saying this, I think I've done all that you asked me on your review. Let's me know if I skipped something by distraction or if you have any further comments. Thanks ;), Rafael

@RafaelNH

This comment has been minimized.

Copy link
Contributor Author

commented Mar 22, 2019

@arokem - just realised that I didn't address the comment about avoiding the if statement in line 88. Bit busy at the moment but I will do it asap!

@skoudoro
Copy link
Member

left a comment

Hi @RafaelNH,

Nice Works! Sorry for the delay but here a quick first review. I need to test this code and read the paper. My main concern during this first review is the 2D images and swap_axis management.

dipy/denoise/gibbs.py Show resolved Hide resolved
raise ValueError("Different slices have to be organized along" +
"one of the 3 first matrix dimensions")
elif slice_axis < 2:
vol = np.swapaxes(vol, slice_axis, 2)

This comment has been minimized.

Copy link
@skoudoro

skoudoro Mar 25, 2019

Member

I did not test it but I think it will crash when vol.ndim == 2? What do you think ?

This comment has been minimized.

Copy link
@RafaelNH

RafaelNH Mar 28, 2019

Author Contributor

Data is already in the format assumed by the function when vol.ndim == 2, therefore you don't need that case in the elif statements. In the tests, I am testing a 2D image and all seems to run fine. I will add in the function documentation that data can also be a ([X, Y]) matrix.

This comment has been minimized.

Copy link
@skoudoro

skoudoro Mar 28, 2019

Member

In your test, you do not swap axis with the 2d data and assume that the user will do the correct move 😄

if you do the same test with swap_axis=0 (which is possible), this should crash. Since you are managing 2D data as well, it will be good to manage this case.

This comment has been minimized.

Copy link
@RafaelNH

RafaelNH Mar 29, 2019

Author Contributor

Very well spotted! I've just fixed the issue! For 2D data, swap_axis is a dummy variable, because for 2D data you need just to run _gibbs_removal_2d. See more information in the comments of my code.

This comment has been minimized.

Copy link
@arokem

arokem Jul 2, 2019

Member

Is this resolved? If yes, please hit the "resolve conversation" button

dipy/denoise/gibbs.py Outdated Show resolved Hide resolved

# Run gibbs removal of 2D images
if nd == 2:
vol = _gibbs_removal_2d(vol, n_points=n_points, G0=G0, G1=G1)

This comment has been minimized.

Copy link
@skoudoro

skoudoro Mar 25, 2019

Member

As said before, your code should crash before reaching here. I propose you to transform your 2d images to 3d by adding a fake dimension on your data when nd==2. If you do that, you do not need this part.

This comment has been minimized.

Copy link
@RafaelNH

RafaelNH Mar 28, 2019

Author Contributor

True! However, to do this change I will need to add two extra if statements, one in "check matrix dimension" and another in "Reshape data to original format". Since my code is not crashing for 2D image I will leave the code as it is.

dipy/denoise/gibbs.py Outdated Show resolved Hide resolved
dipy/denoise/gibbs.py Outdated Show resolved Hide resolved
dipy/denoise/gibbs.py Outdated Show resolved Hide resolved
gibbs_removal, _image_tv)
from numpy.testing import (assert_, assert_array_almost_equal, assert_raises)

# Produce a 2D image

This comment has been minimized.

Copy link
@skoudoro

skoudoro Mar 25, 2019

Member

Can you move this part to a setup_module() function like here. Do not forget to explicit what is global or not

This comment has been minimized.

Copy link
@RafaelNH

RafaelNH Mar 29, 2019

Author Contributor

Very nice suggestion!

assert_array_almost_equal(image4d_cor[:, :, 1, 1], image_cor)


def test_swaped_gibbs_3d():

This comment has been minimized.

Copy link
@skoudoro

skoudoro Mar 25, 2019

Member

Can you add test_swaped_gibbs_2d() ?

This comment has been minimized.

Copy link
@RafaelNH

RafaelNH Mar 29, 2019

Author Contributor

Done!

@RafaelNH

This comment has been minimized.

Copy link
Contributor Author

commented Mar 28, 2019

@arokem - I've just recalled why I am using the if statement in line 88 instead of using the axis information. Basically, with this if statement, indexing in lines 121-135 is easier to read.

@skoudoro - I've addressed most of your comments - I will do the rest asap (2 comments are missing)

@RafaelNH

This comment has been minimized.

Copy link
Contributor Author

commented Mar 29, 2019

@arokem , @skoudoro - I think I've addressed all your comments - let me know if you spot something else that I have to address. Many thanks for your inputs!

@skoudoro skoudoro added this to the 1.0 milestone Apr 23, 2019

computed using the time series values on the right boundary and ptv values
are computed using the time series values on the left boundary.
"""
if axis:

This comment has been minimized.

Copy link
@Borda

Borda Apr 27, 2019

Contributor

what about xs = x.copy() if axis else x.T.copy()

# original grid points
xs[idx] = (isp[idx] - isn[idx])/(sp[idx] + sn[idx])*sn[idx] + isn[idx]

if axis:

This comment has been minimized.

Copy link
@Borda

Borda Apr 27, 2019

Contributor

what about return xs if axis else xs.T


def _weights(shape):
""" Computes the weights necessary to combine two images processed by
the 1D gibbs removal procedure along two different axis [1]_.

This comment has been minimized.

Copy link
@Garyfallidis
K1, K0 = np.meshgrid(k[1:-1], k[1:-1])
cosk0 = 1.0 + np.cos(K0)
cosk1 = 1.0 + np.cos(K1)
G1[1:-1, 1:-1] = cosk0 / (cosk0+cosk1)

This comment has been minimized.

Copy link
@Garyfallidis

Garyfallidis May 7, 2019

Member

cosk0 + cosk1

cosk1 = 1.0 + np.cos(K1)
G1[1:-1, 1:-1] = cosk0 / (cosk0+cosk1)
G0[1:-1, 1:-1] = cosk1 / (cosk0+cosk1)

This comment has been minimized.

Copy link
@Garyfallidis
# check matrix dimension
if nd == 4:
inishap = vol.shape
vol = vol.reshape((inishap[0], inishap[1], inishap[2]*inishap[3]))

This comment has been minimized.

Copy link
@Garyfallidis
----
This function was created to deal with gibbs artefacts of MR images.
Assuming that MR images are reconstructed from estimates of their Fourier
expansion coefficients, during TV calculation matrix x can taken as and

This comment has been minimized.

Copy link
@Garyfallidis
ax.flat[1].set_title('Corrected')

plt.show()
fig1.savefig('Gibbs_suppression_b0.png')

This comment has been minimized.

Copy link
@Garyfallidis

Garyfallidis May 7, 2019

Member

Very hard to see the Gibbs artifact in the original data.

Note
----
This function decreases the effects of gibbs oscillations based on the
analysis of local total variation (TV) along the two axis of the image.

This comment has been minimized.

Copy link
@Garyfallidis

Garyfallidis May 7, 2019

Member

two axes

Although artefact correction is done based on each point primary adjanced
neighbors, total variation should be accessed in a larger range of
neigbors. If you want to adjust the number of the neigbors to be
considered in TV calculation please change parameter nn.

This comment has been minimized.

Copy link
@Garyfallidis

Garyfallidis May 7, 2019

Member

Where is that nn parameter?

@RafaelNH RafaelNH referenced this pull request Jun 9, 2019

@RafaelNH RafaelNH force-pushed the RafaelNH:Gibbs_Removal branch from 00bff7d to 895fcde Jun 20, 2019

@RafaelNH

This comment has been minimized.

Copy link
Contributor Author

commented Jun 25, 2019

@arokem , @Garyfallidis , @skoudoro, @Borda - FYI: I am done with this PR. Continuous integration issues are not related to this PR. Please let me know of there is something else that you want me to address here.

@skoudoro

This comment has been minimized.

Copy link
Member

commented Jun 25, 2019

Overall, looks good to me, I just need to play a bit more with the example but I am ok to merge it by the end of this week if there is no other comment

@arokem
Copy link
Member

left a comment

This is almost ready to go. A few suggestions from me. Let me know when you've integrated this, and we can go ahead and merge. If you could also rebase before the merge, that would be appreciated!

----
This function was created to deal with gibbs artefacts of MR images.
Assuming that MR images are reconstructed from estimates of their Fourier
expansion coefficients, during TV calculation matrix x can taken as and

This comment has been minimized.

Copy link
@arokem

arokem Jul 2, 2019

Member
Suggested change
expansion coefficients, during TV calculation matrix x can taken as and
expansion coefficients, during TV calculation matrix x can taken be taken as a


def _image_tv(x, axis=0, n_points=3):
""" Computes total variation (TV) of matrix x accross a given axis and

This comment has been minimized.

Copy link
@arokem

arokem Jul 2, 2019

Member
Suggested change
""" Computes total variation (TV) of matrix x accross a given axis and
""" Computes total variation (TV) of matrix x across a given axis and
----
This function suppresses the effects of Gibbs oscillations based on the
analysis of local total variation (TV). Although artefact correction is
done based on two adjanced points for each voxel, total variation should be

This comment has been minimized.

Copy link
@arokem

arokem Jul 2, 2019

Member
Suggested change
done based on two adjanced points for each voxel, total variation should be
done based on two adjacent points for each voxel, total variation should be
This function suppresses the effects of Gibbs oscillations based on the
analysis of local total variation (TV). Although artefact correction is
done based on two adjanced points for each voxel, total variation should be
accessed in a larger range of neighbours. The number of the neighbours to

This comment has been minimized.

Copy link
@arokem

arokem Jul 2, 2019

Member
Suggested change
accessed in a larger range of neighbours. The number of the neighbours to
accessed in a larger range of neighbours. The number of neighbours to
cosk1 = 1.0 + np.cos(K1)
G1[1:-1, 1:-1] = cosk0 / (cosk0 + cosk1)
G0[1:-1, 1:-1] = cosk1 / (cosk0 + cosk1)

This comment has been minimized.

Copy link
@arokem

arokem Jul 2, 2019

Member
Suggested change

data_corrected = gibbs_removal(data_slices, slice_axis=2)

""" Due to the high dimensionality of diffusion-weighted data, we recommend you

This comment has been minimized.

Copy link
@arokem

arokem Jul 2, 2019

Member
Suggested change
""" Due to the high dimensionality of diffusion-weighted data, we recommend you
""" Due to the high dimensionality of diffusion-weighted data, we recommend that
data_corrected = gibbs_removal(data_slices, slice_axis=2)

""" Due to the high dimensionality of diffusion-weighted data, we recommend you
to specify which is the axis of data matrix that corresponds to different

This comment has been minimized.

Copy link
@arokem

arokem Jul 2, 2019

Member
Suggested change
to specify which is the axis of data matrix that corresponds to different
you specify which is the axis of data matrix that corresponds to different
:align: center
Uncorrected (left panel) and corrected (middle panel) b-value=0 images. For
a reference, the difference between uncorrected and corrected images are

This comment has been minimized.

Copy link
@arokem

arokem Jul 2, 2019

Member
Suggested change
a reference, the difference between uncorrected and corrected images are
reference, the difference between uncorrected and corrected images is
a reference, the difference between uncorrected and corrected images are
shown in the right panel.
The above figure show that the benefits of suppressing Gibbs artefacts is hard

This comment has been minimized.

Copy link
@arokem

arokem Jul 2, 2019

Member
Suggested change
The above figure show that the benefits of suppressing Gibbs artefacts is hard
The above figure shows that the benefits of suppressing Gibbs artefacts are hard
shown in the right panel.
The above figure show that the benefits of suppressing Gibbs artefacts is hard
to observed on b-value=0 data. Therefore, diffusion derived metrics for both

This comment has been minimized.

Copy link
@arokem

arokem Jul 2, 2019

Member
Suggested change
to observed on b-value=0 data. Therefore, diffusion derived metrics for both
to observe on b-value=0 data. Therefore, diffusion derived metrics for both

@RafaelNH RafaelNH force-pushed the RafaelNH:Gibbs_Removal branch from 2e12577 to a71c467 Jul 12, 2019

RafaelNH added 24 commits Mar 17, 2019
BF, TEST: Case of swap axis in 2D images
In this case, swap axis is a dummy variable

@RafaelNH RafaelNH force-pushed the RafaelNH:Gibbs_Removal branch from a71c467 to 4c1d8d5 Jul 15, 2019

@RafaelNH

This comment has been minimized.

Copy link
Contributor Author

commented Jul 15, 2019

Done! All suggestion addressed, rebased, and example updated according to latest version of dipy.

@arokem

This comment has been minimized.

Copy link
Member

commented Jul 15, 2019

OK. Once the CI goes green, I will merge this.

@arokem arokem merged commit 17b68f0 into nipy:master Jul 16, 2019

3 checks passed

Codacy/PR Quality Review Up to standards. A positive pull request.
Details
continuous-integration/appveyor/pr AppVeyor build succeeded
Details
continuous-integration/travis-ci/pr The Travis CI build passed
Details
@arokem

This comment has been minimized.

Copy link
Member

commented Jul 16, 2019

Thanks @RafaelNH for the strong work!

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Projects
None yet
7 participants
You can’t perform that action at this time.