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

WIP: Local pca and noise estimation #1053

Closed
wants to merge 48 commits into
base: master
from

Conversation

Projects
None yet
8 participants
@riddhishb
Contributor

riddhishb commented May 20, 2016

This is the rough PR for the LocalPCA denoising which I am working on. This not at all the final version, and is for reviewing.
@Garyfallidis @RafaelNH @omarocegueda @villalonreina do have a look at this.

@arokem

This comment has been minimized.

Member

arokem commented May 20, 2016

Exciting to see the beginning of the GSoC project! I realize this is early days. For now, just a naming comment: since the localPCA_denoise module is a submodule of the denoise module, I think that it would be fine to call it localPCA.

@Garyfallidis Garyfallidis changed the title from Local pca and noise estimation to WIP: Local pca and noise estimation May 20, 2016

@Garyfallidis

This comment has been minimized.

Member

Garyfallidis commented May 20, 2016

Hi @riddhishb I changed the title to be WIP (work in progress). That allows to the other developers to see that this not ready to be merged but that people are allowed to look and give suggestions/comments.

import numpy as np
from dipy.denoise.rician_adaptation import rician_adaptation
def localPCA_denoise(arr, sigma, patch_radius=1, tou=0, rician=True):

This comment has been minimized.

@Garyfallidis

Garyfallidis May 20, 2016

Member

Please change the name of the file to localpca.py
And the name of the function to localpca. We do not use capital letters in function names.

@riddhishb

This comment has been minimized.

Contributor

riddhishb commented May 20, 2016

@Garyfallidis @arokem will do :)

Riddhish Bhalodia Riddhish Bhalodia

@arokem arokem added the gsoc2016 label May 21, 2016

Riddhish Bhalodia Riddhish Bhalodia
'''
Local PCA Based Denoising of Diffusion Datasets
References

This comment has been minimized.

@RafaelNH

RafaelNH May 22, 2016

Contributor

References should be after Parameters and Returns

Returns
-------

This comment has been minimized.

@RafaelNH

RafaelNH May 22, 2016

Contributor

Remove unnecessary white space.

denoised_arr : 4D array
this is the denoised array of the same size as that of
arr (input data)

This comment has been minimized.

@RafaelNH

RafaelNH May 22, 2016

Contributor

Same as above

Riddhish Bhalodia added some commits May 24, 2016

Riddhish Bhalodia Riddhish Bhalodia
Corrected rician adaptation and SNR rescaling, added look up table in…
… dipy.data.files and modified fetcher.py
def localpca(arr, sigma, patch_radius=1, tou=0, rician=True):
'''
Local PCA Based Denoising of Diffusion Datasets

This comment has been minimized.

@RafaelNH

RafaelNH May 28, 2016

Contributor

Add the reference in first line description of function: "Local PCA Based Denoising of Diffusion MRI Datasets [1]_"

this is the denoised array of the same size as that of
arr (input data)
'''

This comment has been minimized.

@RafaelNH

RafaelNH May 28, 2016

Contributor

Add a reference field:

References

.. [1] Manjon JV, Coupe P, Concha L, Buades A, Collins DL, 2013. Diffusion Weighted Image Denoising Using Overcomplete Local PCA. PLOS one http://dx.doi.org/10.1371/journal.pone.0073021

import scipy as sp
def localpca(arr, sigma, patch_radius=1, tou=0, rician=True):
'''

This comment has been minimized.

@RafaelNH

RafaelNH May 28, 2016

Contributor

default tou should be None. Below when checking tou you should use:
if tou is None:

@riddhishb

This comment has been minimized.

Contributor

riddhishb commented May 28, 2016

@RafaelNH will do :)

Riddhish Bhalodia added some commits May 28, 2016

Y = X.dot(W)
# When the block covers each pixel identify it into the label matrix theta
X_est = Y.dot(np.transpose(W))
X_est = X_est.dot(D_hat)

This comment has been minimized.

@RafaelNH

RafaelNH May 28, 2016

Contributor

Perhaps a good test will be comparing these lines of code with:
X_est = X.dot(D_hat) # your previous way to compute X_est

@coveralls

This comment has been minimized.

coveralls commented Jul 10, 2016

Coverage Status

Coverage increased (+0.2%) to 83.017% when pulling 9613ae7 on riddhishb:localPCA into 7c236b6 on nipy:master.

@coveralls

This comment has been minimized.

coveralls commented Jul 11, 2016

Coverage Status

Coverage increased (+0.2%) to 83.017% when pulling e4521ed on riddhishb:localPCA into 7c236b6 on nipy:master.

@cython.boundscheck(False)
@cython.wraparound(False)
@cython.cdivision(True)
def fast_lpca(double[:, :, :, :] I, int radius, double[:, :, :, :] sigma):

This comment has been minimized.

@omarocegueda

omarocegueda Jul 12, 2016

Contributor

Are we considering a different sigma field for each gradient?

This comment has been minimized.

@riddhishb

riddhishb Jul 12, 2016

Contributor

no we are not, the function is having 4D for convenience, but the sigma at all directions are same

l0norm = 0
for p in range(ndiff):
if d[p] >= sigma[i0, j0, k0, p] * 2.3 * 2.3:

This comment has been minimized.

@omarocegueda

omarocegueda Jul 12, 2016

Contributor

I don't think having duplicated sigmas is more convenient, it is actually decreasing performance because this product can be precomputed out of the for loop here

This comment has been minimized.

@omarocegueda

omarocegueda Jul 12, 2016

Contributor

Oh! and it may be using a lot of unnecessary memory =s

This comment has been minimized.

@riddhishb

riddhishb Jul 12, 2016

Contributor

Yeah :) great point will change asap.

This comment has been minimized.

@omarocegueda

omarocegueda Jul 12, 2016

Contributor

Let's squeeze up to the last bit of memory =P

double[:, :] cov = np.empty((ndiff, ndiff))
double[:] mu = np.empty(ndiff)
double[:] x
double[:, :, :, :] theta = np.zeros((n0, n1, n2, ndiff))

This comment has been minimized.

@omarocegueda

omarocegueda Jul 12, 2016

Contributor

We also have several copies of theta here, one (n0, n1, n2) volume would be enough

This comment has been minimized.

@riddhishb

riddhishb Jul 12, 2016

Contributor

@omarocegueda I changed that in the code which I will commit soon :), memory saving mode on!

riddhishb added some commits Jul 12, 2016

@coveralls

This comment has been minimized.

coveralls commented Jul 12, 2016

Coverage Status

Coverage increased (+0.2%) to 83.019% when pulling 07782cb on riddhishb:localPCA into 7c236b6 on nipy:master.

@coveralls

This comment has been minimized.

coveralls commented Jul 25, 2016

Coverage Status

Coverage increased (+0.2%) to 83.019% when pulling 2b92916 on riddhishb:localPCA into 7c236b6 on nipy:master.

@cython.boundscheck(False)
@cython.wraparound(False)
@cython.cdivision(True)
cdef void fast_vecmat_mul(double[:] x, double[:, :] A, double[:] out) nogil:

This comment has been minimized.

@omarocegueda

omarocegueda Jul 26, 2016

Contributor

Is it possible to pass raw pointers instead of memory views here?:
fast_vecmat_mul(double *x, double ** A, double* out)

temp[s] = I[p, q, r, s] - mu[s]
fast_vecmat_mul(temp, P, temp1)
update_vector(temp1, mu, 1)

This comment has been minimized.

@omarocegueda

omarocegueda Jul 26, 2016

Contributor

If we refactor fast_vecmat_mult then temp and temp1 will be plain arrays, so we may need another version of update_vector receiving double* too

for p in range(ndiff):
if d[p] >= sigma[i0, j0, k0] * 2.3 * 2.3:
l0norm += 1
update_outer_prod(P, W_hat[:, p], 1)

This comment has been minimized.

@omarocegueda

omarocegueda Jul 26, 2016

Contributor

If we refactor fast_vecmat_mul then P is now double ** so we may need a new update_outer_prod as well

cnp.npy_intp cur_i, prev_i, i, j, k, ii
cnp.npy_intp i0, j0, k0, p, q, r, s
double l0norm
double[:, :, :, :, :] T = np.zeros((2, n1, n2, ndiff, ndiff))

This comment has been minimized.

@omarocegueda

omarocegueda Jul 26, 2016

Contributor

It may be tempting to refactor everything to plain arrays. I just tried but allocating memory resulted in an extreme overhead, so unless we reserve the temporary memory as one single memory block it won't improve performance.

@coveralls

This comment has been minimized.

coveralls commented Jul 27, 2016

Coverage Status

Coverage increased (+0.2%) to 83.019% when pulling 203bb82 on riddhishb:localPCA into 7c236b6 on nipy:master.

@coveralls

This comment has been minimized.

coveralls commented Aug 3, 2016

Coverage Status

Coverage increased (+0.2%) to 83.019% when pulling dae5958 on riddhishb:localPCA into 7c236b6 on nipy:master.

@coveralls

This comment has been minimized.

coveralls commented Aug 4, 2016

Coverage Status

Coverage increased (+0.2%) to 83.018% when pulling 3d8252f on riddhishb:localPCA into 7c236b6 on nipy:master.

Denoise images using Local PCA
===============================
Using the local PCA based denoising for diffusion

This comment has been minimized.

@RafaelNH

RafaelNH Aug 5, 2016

Contributor

Can you rewrite this paragraph? A bit confusing - perhaps it will be nice to give some more information of the PCA denoising techniques.

This comment has been minimized.

@riddhishb

riddhishb Aug 5, 2016

Contributor

@RafaelNH On it!

@coveralls

This comment has been minimized.

coveralls commented Aug 5, 2016

Coverage Status

Coverage increased (+0.2%) to 83.018% when pulling 2cb9df0 on riddhishb:localPCA into 7c236b6 on nipy:master.

@coveralls

This comment has been minimized.

coveralls commented Aug 9, 2016

Coverage Status

Coverage increased (+0.2%) to 83.018% when pulling 6e6b780 on riddhishb:localPCA into 7c236b6 on nipy:master.

Using the local PCA based denoising for diffusion
images [Manjon2013]_ we can denoise diffusion images significantly. This algorithm is
popular as it takes into account the directional information in diffusion data.

This comment has been minimized.

@RafaelNH

RafaelNH Aug 9, 2016

Contributor

Last sentence: 'For diffusion-weighted data this algorithm is advantageous for taking into account the directional information in diffusion data.'

if arr.ndim == 4:
tou = 2.3 * 2.3 * sigma

This comment has been minimized.

@samuelstjean

samuelstjean Sep 29, 2016

Contributor

This should be the variance, not the standard deviation, so your docstring is misnommed or you missed the square. See see http://journals.plos.org/plosone/article?id=10.1371/journal.pone.0073021#s2

After an exhaustive search for the optimum value of τ parameter, it was set to (2.3σ)2, where σ2 is the estimated local noise variance

This comment has been minimized.

@riddhishb

riddhishb Sep 29, 2016

Contributor

I checked it's not named correctly, and I will change that.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment