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

Denoising images using non-local means #320

Merged
merged 43 commits into from Mar 2, 2014

Conversation

Projects
None yet
5 participants
@Garyfallidis
Member

Garyfallidis commented Jan 23, 2014

In this PR I am adding a new general way to apply denoising in 3D and 4D images using NLMEANS. I have rewrote and optimized the algorithm. For the first time we use OpenMP in dipy through cython. OpenMP increases considerable execution time when you want to use multi-threading with shared memory.

A few things that are missing and they should come soon in another PR.

  1. estimate SNR for each of the 3D volumes of a DWI dataset
  2. estimate SNR for T1 and add an example of T1 denoising in the existing example.
  3. explain that we first estimate the noise and then apply the mask
  4. add any other references.
  5. extend the method to its adaptive version which is what we used for the ISBI competition
  6. make sure that the results correspond with Pierrick's matlab code.

@samuelstjean will work on these steps as I need to focus on a different project and denoising is his main interest.

GGT!

@stefanv

This comment has been minimized.

Contributor

stefanv commented Jan 23, 2014

Unless you have awesome models, of course :)

# calculate weights between the central patch and the moving patch in block
# (m, n, o) coordinates are the center of the moving patch
# (a, b, c) run incide both patches

This comment has been minimized.

@samuelstjean

samuelstjean Feb 15, 2014

Contributor

inside

block_radius : int
block size is ``2 x block_radius + 1``. Default is 5.
rician : boolean
If True the noise is estimate as Rician, otherwise Gaussian noise

This comment has been minimized.

@samuelstjean

samuelstjean Feb 15, 2014

Contributor

is estimated

Using the non-local means filter [Coupe2008]_ you can denoise 3D or 4D images and
boost the SNR of your datasets. You can also decide between modeling the noise
as Gaussian or Rician (default).

This comment has been minimized.

@samuelstjean

samuelstjean Feb 15, 2014

Contributor

Actually, the only difference is the substraction of 2_sigma_*2 at the end, so it should be something like removing the bias introducedby the rician noise instead of really modelling stuff.

This comment has been minimized.

@Garyfallidis

Garyfallidis Feb 16, 2014

Member

I don't see why what you suggest makes better sense. Please be more specific.

This comment has been minimized.

@samuelstjean

samuelstjean Feb 16, 2014

Contributor

Because it's not noise modelling, it removes the rician introduced bias, but in no way does any modelling of any sort. Or it depends on your definition of modelling :P

@samuelstjean

This comment has been minimized.

Contributor

samuelstjean commented Feb 15, 2014

Seems good, it change the code and example afterward I guess.

@mdesco

This comment has been minimized.

Contributor

mdesco commented Feb 16, 2014

Another question/comment. I just denoised a 500Mb full brain dataset and my RAM consumption went up to 3 Giga. So, a factor of 6 increase.

I see that I inputed an float32 file and the output was a float64. So, this might be contributing to the issue.

Do you think that if the input is int16 or float32, the output should as well? My old binary used to do that.

Still fast and good results!

@samuelstjean

This comment has been minimized.

Contributor

samuelstjean commented Feb 16, 2014

Should be the responsability of the user to save it himself as another dtype I'd say. Dunno for the ram usage though.

@Garyfallidis

This comment has been minimized.

Member

Garyfallidis commented Feb 16, 2014

@mdesco good catch. Please give it a try with the a new fix. It should be okay now. Please let me know otherwise.

@mdesco

This comment has been minimized.

Contributor

mdesco commented Feb 16, 2014

I still have a large RAM consumption. 251Mb input and my RAM goes up to 3Gig during the processing.

@mdesco mdesco closed this Feb 16, 2014

@mdesco mdesco reopened this Feb 16, 2014

@Garyfallidis

This comment has been minimized.

Member

Garyfallidis commented Feb 16, 2014

What type of dataset is the one that you use? Is it a 4D image of dMRI data?

@Garyfallidis

This comment has been minimized.

Member

Garyfallidis commented Feb 16, 2014

Hi Max, I am not able to replicate your problem. I monitored the memory usage with a 4D dataset of 400MB and the memory needed at the end and during execution time was an extra 400MB as expected. So, the total amount of memory was ~800MB. I would suggest maybe to recompile your dipy installation after having fetch my latest changes for this brunch. If that doesn't work then let's have a look together at IMEKA.

@mdesco

This comment has been minimized.

Contributor

mdesco commented Feb 23, 2014

All good on my side. It works fine and nlmeans supports 4D. Thanks.

@Garyfallidis

This comment has been minimized.

Member

Garyfallidis commented Feb 24, 2014

Okay, guys. This is ready on my side. Can someone who is not from the SCIL give it a last review and merge this baby?

@@ -0,0 +1,8 @@
#init for denoisa aka the denoising module

This comment has been minimized.

@arokem

arokem Feb 24, 2014

Member

Typo : "denoise"

@arokem

This comment has been minimized.

Member

arokem commented Feb 24, 2014

I will take a look.

On Mon, Feb 24, 2014 at 6:19 AM, Eleftherios Garyfallidis <
notifications@github.com> wrote:

Okay, guys. This is ready on my side. Can someone who is not from the SCIL
give it a last review and merge this baby?


Reply to this email directly or view it on GitHubhttps://github.com//pull/320#issuecomment-35889891
.

the denoised ``arr`` which has the same shape as ``arr``.
"""
if arr.ndim != 3:

This comment has been minimized.

@arokem

arokem Feb 24, 2014

Member

The docstring suggests that either 4D or 3D arrays are acceptable inputs. Which is it?

This comment has been minimized.

@Garyfallidis

Garyfallidis Mar 2, 2014

Member

Good catch. Thx!

cdef double process_block(double [:, :, ::1] arr,
cnp.npy_intp i, cnp.npy_intp j, cnp.npy_intp k,
cnp.npy_intp B, cnp.npy_intp P, double sigma) nogil:

This comment has been minimized.

@arokem

arokem Feb 24, 2014

Member

docstring?

patch_radius, block_radius,
rician).astype(arr.dtype)
if arr.ndim == 4:

This comment has been minimized.

@arokem

arokem Feb 24, 2014

Member

OK - this makes sense

assert_array_almost_equal)
from dipy.denoise.nlmeans import nlmeans
from dipy.denoise.denspeed import add_padding_reflection, remove_padding
from matplotlib.pyplot import *

This comment has been minimized.

@arokem

arokem Feb 24, 2014

Member

Why do you need matplotlib to test this? I assume that you did some debugging with plotting at some point, but this import might not be necessary any more. If it is, please set it up with the tripwire mechanism. Also from foo import * is a no no. Consider importing what you need or a namespace.

This comment has been minimized.

@Garyfallidis

Garyfallidis Mar 2, 2014

Member

True, sorry my bad. Will clean it asap!

"""
In order to call ``nlmeans`` first you need to estimate the standard deviation
of the noise. TODO

This comment has been minimized.

@arokem

arokem Feb 24, 2014

Member

I don't think we can merge an example that has a TODO in it. Please do this. There's some code that does something like this in the RESTORE example.

This comment has been minimized.

@samuelstjean

samuelstjean Feb 24, 2014

Contributor

For this one, I'm gonna replace it with a better example using a robust estimation of the noise instead, so either we remove it for nom and I'll do the changes on my part or it's just gonna be replaced in the next month.

This comment has been minimized.

@arokem

arokem Feb 24, 2014

Member

OK - but please don't forget to circle back to this as soon as you get a
chance.

On Mon, Feb 24, 2014 at 9:36 AM, Samuel St-Jean notifications@github.comwrote:

In doc/examples/denoise_nlmeans.py:

+img, gtab = read_sherbrooke_3shell()
+
+data = img.get_data()
+aff = img.get_affine()
+
+mask = data[..., 0] > 100
+
+data = data[..., 0]
+
+print("vol size", data.shape)
+
+t = time()
+
+"""
+In order to call nlmeans first you need to estimate the standard deviation
+of the noise. TODO

For this one, I'm gonna replace it with a better example using a robust
estimation of the noise instead, so either we remove it for nom and I'll do
the changes on my part or it's just gonna be replaced in the next month.


Reply to this email directly or view it on GitHubhttps://github.com//pull/320/files#r10000657
.

import matplotlib.pyplot as plt
axial_middle = data.shape[2] / 2
plt.figure('Showing the datasets')

This comment has been minimized.

@arokem

arokem Feb 24, 2014

Member

One more comment here: I think that it is generally preferable to use the following subplots API in matplotlib:

fig, ax = plt.subplots(2)
ax[0].imshow(data[:, :, axial_middle].T, cmap='gray', origin='lower')
ax[0].set_axis_off()
ax[1].imshow(den[:, :, axial_middle].T, cmap='gray', origin='lower')
ax[1].set_axis_off()
fig.savefig('denoised_S0.png')
@arokem

This comment has been minimized.

Member

arokem commented Feb 24, 2014

I am getting this error when running the example. Are you sure the input to the bbox_inches kwarg makes sense?

AssertionError Traceback (most recent call last)
/Users/arokem/Desktop/denoise_nlmeans.py in ()
49 plt.imshow(den[:, :, axial_middle].T, cmap='gray', origin='lower')
50 plt.show()
---> 51 plt.savefig('denoised_S0.png', bbox_inches='tight')
52
53 """

/Users/arokem/Library/Enthought/Canopy_64bit/User/lib/python2.7/site-packages/matplotlib/pyplot.pyc in savefig(_args, *_kwargs)
559 def savefig(_args, *_kwargs):
560 fig = gcf()
--> 561 return fig.savefig(_args, *_kwargs)
562
563

/Users/arokem/Library/Enthought/Canopy_64bit/User/lib/python2.7/site-packages/matplotlib/figure.pyc in savefig(self, _args, *_kwargs)
1419 self.set_frameon(frameon)
1420
-> 1421 self.canvas.print_figure(_args, *_kwargs)
1422
1423 if frameon:

/Users/arokem/Library/Enthought/Canopy_64bit/User/lib/python2.7/site-packages/matplotlib/backend_bases.pyc in print_figure(self, filename, dpi, facecolor, edgecolor, orientation, format, *_kwargs)
2167 *_kwargs)
2168 renderer = self.figure._cachedRenderer
-> 2169 bbox_inches = self.figure.get_tightbbox(renderer)
2170
2171 bbox_artists = kwargs.pop("bbox_extra_artists", None)

/Users/arokem/Library/Enthought/Canopy_64bit/User/lib/python2.7/site-packages/matplotlib/figure.pyc in get_tightbbox(self, renderer)
1562 bb.append(ax.get_tightbbox(renderer))
1563
-> 1564 _bbox = Bbox.union([b for b in bb if b.width != 0 or b.height != 0])
1565
1566 bbox_inches = TransformedBbox(_bbox,

/Users/arokem/Library/Enthought/Canopy_64bit/User/lib/python2.7/site-packages/matplotlib/transforms.pyc in union(bboxes)
712 Return a :class:Bbox that contains all of the given bboxes.
713 """
--> 714 assert(len(bboxes))
715
716 if len(bboxes) == 1:

@Garyfallidis

This comment has been minimized.

Member

Garyfallidis commented Mar 2, 2014

Hi @arokem I added all your suggestions except from the bbox_inches='tight'. I think this problem is on your side. This option is recommended in many places. Let me know if the problem persists. Could it be that you have some problems with your matplotlib installation or version? My version is 1.2.1.

Also if you know a better way to remove the padding from the figure let me know.

arokem added a commit that referenced this pull request Mar 2, 2014

Merge pull request #320 from Garyfallidis/nlmeans_opt
Denoising images using non-local means

@arokem arokem merged commit c068cb2 into nipy:master Mar 2, 2014

1 check passed

default The Travis CI build passed
Details
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment