Added eigh version of localpca to svd version #1297

ssheybani · 2017-07-10T19:35:41Z

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?).

arokem · 2017-07-10T19:46:52Z

dipy/denoise/localpca.py

    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 was intentional! To denote that we do this "if and only if".

arokem · 2017-07-10T19:47:34Z

dipy/denoise/localpca.py

@@ -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:


arokem · 2017-07-10T19:48:05Z

dipy/denoise/localpca.py

@@ -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


Maybe put back the original values?

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.

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 · 2017-07-10T21:09:34Z

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.

Fix typos in method documentation for the dti.py reconstruction file.

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 · 2017-07-14T20:47:34Z

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

coveralls · 2017-07-14T20:47:34Z

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

codecov-io · 2017-07-14T20:56:22Z

Codecov Report

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

@@            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 · 2017-07-14T22:00:38Z

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:

coveralls · 2017-07-14T23:17:17Z

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

…FormulainDoc DOC: Add missing label to reciprocal space eq.

Garyfallidis · 2017-07-17T18:18:33Z

dipy/denoise/localpca.py

+        # If mask is not specified, use the whole volume
+        mask = np.ones_like(arr, dtype=bool)[..., 0]
+
+    if pca_method == 'svd':


Please add the imports at the top of the file.

Garyfallidis · 2017-07-17T18:18:56Z

dipy/denoise/tests/test_lpca.py

@@ -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))

+


Add a test with the mask.

Garyfallidis · 2017-07-17T18:23:52Z

dipy/denoise/localpca.py

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


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

DOC: Fix typos in dti.py reconstruction file doc.

coveralls · 2017-07-21T21:44:01Z

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

coveralls · 2017-07-21T21:44:01Z

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

coveralls · 2017-07-21T22:59:50Z

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

coveralls · 2017-07-21T22:59:50Z

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

Fix inversion in the dti mode doc

arokem · 2017-07-28T04:43:11Z

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

coveralls · 2017-07-28T21:30:52Z

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

coveralls · 2017-07-28T21:30:52Z

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.

arokem · 2017-07-30T19:19:04Z

dipy/denoise/tests/test_lpca.py

    # 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)


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.

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?

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.

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

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)

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.

arokem · 2017-07-30T19:20:04Z

dipy/denoise/localpca.py

    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.


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?

Please add that to the docstring.

Garyfallidis · 2017-08-01T20:43:13Z

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

coveralls · 2017-08-04T21:30:20Z

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

coveralls · 2017-08-04T21:30:20Z

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

arokem · 2017-08-04T23:34:30Z

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.

Garyfallidis · 2017-08-14T21:07:54Z

Great! Thank you @ssheybani :)

Added eigh version of localpca to svd version

Added eigh version of localpca to svd version

64aff62

arokem requested changes Jul 10, 2017

View reviewed changes

Jon Haitz Legarreta and others added 3 commits July 14, 2017 12:34

DOC: Fix typos in dti.py reconstruction file doc.

25489b2

Fix typos in method documentation for the dti.py reconstruction file.

Made the mask optional, Corrected the tests

4414ba8

Corrected if back to iff

96ceccc

Merge pull request dipy#1303 from jhlegarreta/CrossRefReciprocalSpace…

c155950

…FormulainDoc DOC: Add missing label to reciprocal space eq.

Garyfallidis reviewed Jul 17, 2017

View reviewed changes

RafaelNH and others added 6 commits July 17, 2017 23:09

Merge pull request dipy#1304 from jhlegarreta/FixDTIReconstDocTypos

08d8d90

DOC: Fix typos in dti.py reconstruction file doc.

Added eigh version of localpca to svd version

cdebd37

Made the mask optional, Corrected the tests

3279e2b

Corrected if back to iff

a28620b

Merge branch 'master' of github.com:ssheybani/dipy

fa6ac5a

A few small fixes for PEP8 consistency

71b5061

matthieudumont and others added 2 commits July 26, 2017 08:52

Fix inversion in the dti mode doc

cf1c232

Merge pull request dipy#1308 from matthieudumont/FixDTIModeDoc

59f596a

Fix inversion in the dti mode doc

Added eigh version of localpca to svd version

b17ab32

ssheybani added 6 commits July 28, 2017 14:53

Made the mask optional, Corrected the tests

d512cb8

Corrected if back to iff

2cb81f4

A few small fixes for PEP8 consistency

3f0bd02

Included mask in the test

b064938

Merge branch 'master' of github.com:ssheybani/dipy

6fdd37e

fix:Included mask in the test

e710241

arokem reviewed Jul 30, 2017

View reviewed changes

ssheybani added 3 commits August 3, 2017 10:35

Added False values to the mask of the test

57c0cf6

Added False values to the mask of the test

f42dd46

Corrected the test to avoid the rmse error

be0d7de

arokem approved these changes Aug 9, 2017

View reviewed changes

Garyfallidis merged commit c268bd2 into dipy:master Aug 14, 2017

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

Merge pull request dipy#1297 from ssheybani/master

e26fb9f

Added eigh version of localpca to svd version

		@@ -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))

Added eigh version of localpca to svd version #1297

Added eigh version of localpca to svd version #1297

Conversation

ssheybani commented Jul 10, 2017 • edited

arokem left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Garyfallidis commented Jul 10, 2017

coveralls commented Jul 14, 2017

coveralls commented Jul 14, 2017

codecov-io commented Jul 14, 2017 • edited

Codecov Report

ssheybani commented Jul 14, 2017 • edited

A summary of line profiling:

With svd:

With eig:

coveralls commented Jul 14, 2017 • edited

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

coveralls commented Jul 21, 2017

coveralls commented Jul 21, 2017

coveralls commented Jul 21, 2017

coveralls commented Jul 21, 2017

arokem commented Jul 28, 2017

coveralls commented Jul 28, 2017

coveralls commented Jul 28, 2017

arokem left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

ssheybani Aug 4, 2017 • edited

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Garyfallidis commented Aug 1, 2017

coveralls commented Aug 4, 2017

coveralls commented Aug 4, 2017

arokem commented Aug 4, 2017

Garyfallidis commented Aug 14, 2017

ssheybani commented Jul 10, 2017 •

edited

codecov-io commented Jul 14, 2017 •

edited

ssheybani commented Jul 14, 2017 •

edited

coveralls commented Jul 14, 2017 •

edited

ssheybani Aug 4, 2017 •

edited