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

Added eigh version of localpca to svd version #1297

Merged
merged 24 commits into from Aug 14, 2017

Conversation

Projects
None yet
7 participants
@ssheybani
Contributor

ssheybani commented Jul 10, 2017

Denoising with LocalPCA:
The bottleneck of LocalPCA is decomposition of the data into the eigenvalues/singular values and the corresponding vectors (basis). For this purpose, Eigenvalue decomposition was used but was previously replaced with Singular Values Decomposition, for more accuracy. However, SVD is many times slower than EVD. In this pull request, the code is changed to include both of them, as an optional parameter, "pca_method".
Also for more speed up, the function accepts the brain mask as an optional parameter and the denoising will be only applied to the voxels inside the mask.

@arokem

Looks good! Please rewrite the tests to make sure that they run both methods (and maybe that they give the same answer?).

if out_dtype is None:
out_dtype = arr.dtype
# We retain float64 precision, iff the input is in this precision:
# We retain float64 precision, if the input is in this precision:

This comment has been minimized.

@arokem

arokem Jul 10, 2017

Member

This was intentional! To denote that we do this "if and only if".

@@ -97,6 +105,8 @@ def localpca(arr, sigma, patch_radius=2, tau_factor=2.3, out_dtype=None):
for j in range(patch_radius, arr.shape[1] - patch_radius):
for i in range(patch_radius, arr.shape[0] - patch_radius):
# Shorthand for indexing variables:
if mask[i, j, k] == False:

This comment has been minimized.

@arokem

arokem Jul 10, 2017

Member

Nice!

@@ -133,4 +152,5 @@ def localpca(arr, sigma, patch_radius=2, tau_factor=2.3, out_dtype=None):
denoised_arr = thetax / theta
denoised_arr.clip(min=0, out=denoised_arr)
denoised_arr[~mask] = 0

This comment has been minimized.

@arokem

arokem Jul 10, 2017

Member

Maybe put back the original values?

This comment has been minimized.

@Garyfallidis

Garyfallidis Jul 12, 2017

Member

Nope, this is an application of the mask. If we put back the original values we will have noise in the background which is undesired.

This comment has been minimized.

@arokem

arokem Jul 13, 2017

Member

That makes sense. Might be worth being explicit about this in the docstring. Some users might be surprised to see non-brain parts of the image (e.g., skull) disappear.

@Garyfallidis

This comment has been minimized.

Member

Garyfallidis commented Jul 10, 2017

Please also @ssheybani include more information in the description and show how eigh compares with svd in some slices (not in tests - just in the description for reference and helping the reviewers understand). We may have discussed some things in person but that doesn't mean that the others here know what we have discussed.

Jon Haitz Legarreta and others added some commits Jul 14, 2017

Jon Haitz Legarreta
DOC: Fix typos in dti.py reconstruction file doc.
Fix typos in method documentation for the dti.py reconstruction file.
Jon Haitz Legarreta
DOC: Add missing label to reciprocal space eq.
Add missing label to reciprocal space equation in b_and_q.rst file in
the doc/theory folder.

Although according to
http://www.sphinx-doc.org/en/stable/ext/math.html
cross-referencing equations via their label will not work across
different documents.

The equation at issue is referenced from doc/faq.rst.
@coveralls

This comment has been minimized.

coveralls commented Jul 14, 2017

Coverage Status

Coverage increased (+0.001%) to 85.436% when pulling 4414ba8 on ssheybani:master into 736cbb5 on nipy:master.

1 similar comment
@coveralls

This comment has been minimized.

coveralls commented Jul 14, 2017

Coverage Status

Coverage increased (+0.001%) to 85.436% when pulling 4414ba8 on ssheybani:master into 736cbb5 on nipy:master.

@codecov-io

This comment has been minimized.

codecov-io commented Jul 14, 2017

Codecov Report

Merging #1297 into master will decrease coverage by 0.09%.
The diff coverage is 94.59%.

Impacted file tree graph

@@            Coverage Diff            @@
##           master    #1297     +/-   ##
=========================================
- Coverage   87.13%   87.03%   -0.1%     
=========================================
  Files         228      228             
  Lines       28800    28834     +34     
  Branches     3093     3099      +6     
=========================================
+ Hits        25094    25095      +1     
- Misses       3003     3035     +32     
- Partials      703      704      +1
Impacted Files Coverage Δ
dipy/reconst/dti.py 96.5% <ø> (ø) ⬆️
dipy/denoise/tests/test_lpca.py 98.7% <100%> (+0.09%) ⬆️
dipy/denoise/localpca.py 93.05% <90%> (-2.03%) ⬇️
dipy/io/tests/test_io_peaks.py 89.81% <0%> (-7.66%) ⬇️
dipy/reconst/shore.py 85.15% <0%> (-4.49%) ⬇️
dipy/io/peaks.py 89.79% <0%> (-2.66%) ⬇️
dipy/core/graph.py 73.8% <0%> (-1.2%) ⬇️
dipy/reconst/mapmri.py 89.39% <0%> (-0.87%) ⬇️
dipy/io/dpy.py 94.02% <0%> (+1.38%) ⬆️

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update 736cbb5...be0d7de. Read the comment docs.

@ssheybani

This comment has been minimized.

Contributor

ssheybani commented Jul 14, 2017

Here are the results of comparison between the two methods (eig vs svd). The test data is a small portion of a Ferret dataset, with the shape [160, 120, 7, 210] (28224000 voxels in total)

A summary of line profiling:

With svd:

Total time: 450.4 s
Line #      Hits         Time  Per Hit   % Time  Line Contents
==============================================================
   126      7731    419366666  54244.8     93.1                      U, S, Vt = svd(X, *svd_args)[:3]
   144      7731     20806916   2691.4      4.6                  Xest = X.dot(W).dot(W.T) + M

With eig:

Total time: 83.4777 s
Line #      Hits         Time  Per Hit   % Time  Line Contents
==============================================================
137      7731      7992914   1033.9      9.6                      C = np.transpose(X).dot(X)
139      7731     47239798   6110.4     56.6                      [d, W] = eigh(C, turbo=True)
144      7731     18560852   2400.8     22.2                  Xest = X.dot(W).dot(W.T) + M

The denoised matrices resulted from the two methods were exactly the same (according to numpy.testing.assert_equal).
Here's one of the slices:
eigh_svd_compare_ferret
eigh_svd_compare_ferret_res

@coveralls

This comment has been minimized.

coveralls commented Jul 14, 2017

Coverage Status

Coverage increased (+0.001%) to 85.436% when pulling 96ceccc on ssheybani:master into 736cbb5 on nipy:master.

Merge pull request #1303 from jhlegarreta/CrossRefReciprocalSpaceForm…
…ulainDoc

DOC: Add missing label to reciprocal space eq.
# If mask is not specified, use the whole volume
mask = np.ones_like(arr, dtype=bool)[..., 0]
if pca_method == 'svd':

This comment has been minimized.

@Garyfallidis

Garyfallidis Jul 17, 2017

Member

Please add the imports at the top of the file.

@@ -223,13 +228,13 @@ def test_phantom():
rmse_noisy_wrc = np.sum(np.abs(DWI_clean_wrc - DWI)) / \
np.sum(np.abs(DWI_clean_wrc))

This comment has been minimized.

@Garyfallidis

Garyfallidis Jul 17, 2017

Member

Add a test with the mask.

# order:
W = Vt[::-1].T
if pca_method == 'svd':

This comment has been minimized.

@Garyfallidis

Garyfallidis Jul 17, 2017

Member

No need to check the string each time. Add a boolean variable.

@coveralls

This comment has been minimized.

coveralls commented Jul 21, 2017

Coverage Status

Coverage decreased (-0.0004%) to 85.434% when pulling 71b5061 on ssheybani:master into 736cbb5 on nipy:master.

1 similar comment
@coveralls

This comment has been minimized.

coveralls commented Jul 21, 2017

Coverage Status

Coverage decreased (-0.0004%) to 85.434% when pulling 71b5061 on ssheybani:master into 736cbb5 on nipy:master.

@coveralls

This comment has been minimized.

coveralls commented Jul 21, 2017

Coverage Status

Coverage decreased (-0.0004%) to 85.434% when pulling 71b5061 on ssheybani:master into 736cbb5 on nipy:master.

1 similar comment
@coveralls

This comment has been minimized.

coveralls commented Jul 21, 2017

Coverage Status

Coverage decreased (-0.0004%) to 85.434% when pulling 71b5061 on ssheybani:master into 736cbb5 on nipy:master.

matthieudumont and others added some commits Jul 26, 2017

Merge pull request #1308 from matthieudumont/FixDTIModeDoc
Fix inversion in the dti mode doc
@arokem

This comment has been minimized.

Member

arokem commented Jul 28, 2017

Where does this stand? Looks like we are still missing a test with the mask, but that this is otherwise ready to go?

@coveralls

This comment has been minimized.

coveralls commented Jul 28, 2017

Coverage Status

Coverage increased (+0.0001%) to 85.435% when pulling e710241 on ssheybani:master into 736cbb5 on nipy:master.

1 similar comment
@coveralls

This comment has been minimized.

coveralls commented Jul 28, 2017

Coverage Status

Coverage increased (+0.0001%) to 85.435% when pulling e710241 on ssheybani:master into 736cbb5 on nipy:master.

@arokem

Small additional comments about mask documentation and test.

# Try this with a sigma volume, instead of a scalar
sigma_vol = sigma * np.ones(DWI.shape[:-1])
DWI_den = localpca(DWI, sigma_vol, patch_radius=3)
mask = np.ones_like(DWI, dtype=bool)[..., 0]
DWI_den = localpca(DWI, sigma_vol, mask, patch_radius=3)

This comment has been minimized.

@arokem

arokem Jul 30, 2017

Member

Right now, this doesn't really test that the mask does anything. Maybe make a mask with some ones and some zeros and check that the part that has ones in it gets denoised, while the part that has zeros returns as all zeros.

This comment has been minimized.

@ssheybani

ssheybani Aug 1, 2017

Contributor

I did that. But assert_(rmse_den < rmse_noisy) gave me an error. I think that's most likely because setting the out-of-mask areas to zero has caused the denoised image to have more error than the original noisy one. I checked the DWI_clean matrix and it looks like that for a good number of slices, the background is not actually black (0) but white (1)! Do you have any suggestions on what is the preferred way to deal with it?

This comment has been minimized.

@arokem

arokem Aug 1, 2017

Member

Sounds to me like the mask is doing something funky! I am not sure where that arises. Could you implement this and push it here? I might be able to fetch your branch and play around with it.

This comment has been minimized.

@ssheybani

ssheybani Aug 3, 2017

Contributor

I just added this for instance:
mask[0, ...] = False

This comment has been minimized.

@arokem

arokem Aug 4, 2017

Member

I get this failure only for the second comparison:

assert_(rmse_den_wrc < rmse_noisy_wrc)

not for the first one:

assert_(rmse_den < rmse_noisy)

This comment has been minimized.

@ssheybani

ssheybani Aug 4, 2017

Contributor

I just corrected it . What I've done is when it comes to calculating rmse values, I applied the mask on the original data as well, to avoid the problem with the background. Thus it didn't raise any error.

r"""Local PCA-based denoising of diffusion datasets.
Parameters
----------
arr : 4D array
Array of data to be denoised. The dimensions are (X, Y, Z, N), where N
are the diffusion gradient directions.
mask : 3D boolean array
A mask with voxels that are true inside the brain and false outside of
it.

This comment has been minimized.

@arokem

arokem Jul 30, 2017

Member

I'm still missing an explanation of what the mask does here: it denoises within the True part and returns zeros outside of those voxels, right?

This comment has been minimized.

@ssheybani

ssheybani Aug 1, 2017

Contributor

Yes.

This comment has been minimized.

@arokem

arokem Aug 1, 2017

Member

Please add that to the docstring.

This comment has been minimized.

@ssheybani

ssheybani Aug 3, 2017

Contributor

Done

@Garyfallidis

This comment has been minimized.

Member

Garyfallidis commented Aug 1, 2017

Ping @ssheybani there are some final requests here. It would be great to merge this asap. Thanks.

@coveralls

This comment has been minimized.

coveralls commented Aug 4, 2017

Coverage Status

Coverage decreased (-0.09%) to 85.344% when pulling be0d7de on ssheybani:master into 736cbb5 on nipy:master.

1 similar comment
@coveralls

This comment has been minimized.

coveralls commented Aug 4, 2017

Coverage Status

Coverage decreased (-0.09%) to 85.344% when pulling be0d7de on ssheybani:master into 736cbb5 on nipy:master.

@arokem

This comment has been minimized.

Member

arokem commented Aug 4, 2017

Yup. I think that was it.

This is ready to merge from my point of view. @Garyfallidis : feel free to merge if all your comments have been addressed.

@arokem

arokem approved these changes Aug 9, 2017

@Garyfallidis

This comment has been minimized.

Member

Garyfallidis commented Aug 14, 2017

Great! Thank you @ssheybani :)

@Garyfallidis Garyfallidis merged commit c268bd2 into nipy:master Aug 14, 2017

3 of 4 checks passed

coverage/coveralls Coverage decreased (-0.09%) to 85.344%
Details
codecov/patch 94.59% of diff hit (target 87.13%)
Details
codecov/project Absolute coverage decreased by -0.09% but relative coverage increased by +7.46% compared to 736cbb5
Details
continuous-integration/travis-ci/pr The Travis CI build passed
Details

ShreyasFadnavis pushed a commit to ShreyasFadnavis/dipy that referenced this pull request Sep 20, 2018

Merge pull request nipy#1297 from ssheybani/master
Added eigh version of localpca to svd version
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment