diff --git a/doc/examples/plot_nonlocal_means.py b/doc/examples/plot_nonlocal_means.py
index 790c9bdf..c6879be4 100644
--- a/doc/examples/plot_nonlocal_means.py
+++ b/doc/examples/plot_nonlocal_means.py
@@ -15,7 +15,7 @@ import numpy as np
 import matplotlib.pyplot as plt
 
 from skimage import data, img_as_float
-from skimage.filter import nl_means_denoising
+from skimage.restoration import nl_means_denoising
 
 
 lena = img_as_float(data.lena())
diff --git a/skimage/filters/setup.py b/skimage/filters/setup.py
index 21728902..8bfebcb7 100644
--- a/skimage/filters/setup.py
+++ b/skimage/filters/setup.py
@@ -14,7 +14,6 @@ def configuration(parent_package='', top_path=None):
     config.add_data_dir('rank/tests')
 
     cython(['_ctmf.pyx'], working_path=base_path)
-    cython(['_nl_means_denoising.pyx'], working_path=base_path)
     cython(['rank/core_cy.pyx'], working_path=base_path)
     cython(['rank/generic_cy.pyx'], working_path=base_path)
     cython(['rank/percentile_cy.pyx'], working_path=base_path)
@@ -22,9 +21,6 @@ def configuration(parent_package='', top_path=None):
 
     config.add_extension('_ctmf', sources=['_ctmf.c'],
                          include_dirs=[get_numpy_include_dirs()])
-    config.add_extension('_nl_means_denoising',
-        sources=['_nl_means_denoising.c'],
-        include_dirs=[get_numpy_include_dirs()])
     config.add_extension('rank.core_cy', sources=['rank/core_cy.c'],
         include_dirs=[get_numpy_include_dirs()])
     config.add_extension('rank.generic_cy', sources=['rank/generic_cy.c'],
diff --git a/skimage/restoration/__init__.py b/skimage/restoration/__init__.py
index 2b593ccf..b008c638 100644
--- a/skimage/restoration/__init__.py
+++ b/skimage/restoration/__init__.py
@@ -22,7 +22,7 @@ from .deconvolution import wiener, unsupervised_wiener, richardson_lucy
 from .unwrap import unwrap_phase
 from ._denoise import denoise_tv_chambolle, denoise_tv_bregman, \
                       denoise_bilateral
-
+from .nl_means_denoising import nl_means_denoising
 
 __all__ = ['wiener',
            'unsupervised_wiener',
@@ -30,4 +30,5 @@ __all__ = ['wiener',
            'unwrap_phase',
            'denoise_tv_bregman',
            'denoise_tv_chambolle',
-           'denoise_bilateral']
+           'denoise_bilateral',
+           'nl_means_denoising']
diff --git a/skimage/filter/_nl_means_denoising.pyx b/skimage/restoration/_nl_means_denoising.pyx
similarity index 100%
rename from skimage/filter/_nl_means_denoising.pyx
rename to skimage/restoration/_nl_means_denoising.pyx
diff --git a/skimage/filter/nl_means_denoising.py b/skimage/restoration/nl_means_denoising.py
similarity index 100%
rename from skimage/filter/nl_means_denoising.py
rename to skimage/restoration/nl_means_denoising.py
diff --git a/skimage/restoration/setup.py b/skimage/restoration/setup.py
index e20073e0..12638107 100644
--- a/skimage/restoration/setup.py
+++ b/skimage/restoration/setup.py
@@ -17,6 +17,7 @@ def configuration(parent_package='', top_path=None):
     cython(['_unwrap_2d.pyx'], working_path=base_path)
     cython(['_unwrap_3d.pyx'], working_path=base_path)
     cython(['_denoise_cy.pyx'], working_path=base_path)
+    cython(['_nl_means_denoising.pyx'], working_path=base_path)
 
     config.add_extension('_unwrap_1d', sources=['_unwrap_1d.c'],
                          include_dirs=[get_numpy_include_dirs()])
@@ -28,6 +29,9 @@ def configuration(parent_package='', top_path=None):
                          include_dirs=[get_numpy_include_dirs()])
     config.add_extension('_denoise_cy', sources=['_denoise_cy.c'],
         include_dirs=[get_numpy_include_dirs(), '../_shared'])
+    config.add_extension('_nl_means_denoising',
+                         sources=['_nl_means_denoising.c'],
+                         include_dirs=[get_numpy_include_dirs()])
 
     return config
 
diff --git a/skimage/restoration/tests/test_denoise.py b/skimage/restoration/tests/test_denoise.py
index 46ed4b35..77d38557 100644
--- a/skimage/restoration/tests/test_denoise.py
+++ b/skimage/restoration/tests/test_denoise.py
@@ -46,7 +46,8 @@ def test_denoise_tv_chambolle_float_result_range():
     img = astro_gray
     int_astro = np.multiply(img, 255).astype(np.uint8)
     assert np.max(int_astro) > 1
-    denoised_int_astro = restoration.denoise_tv_chambolle(int_astro, weight=60.0)
+    denoised_int_astro = restoration.denoise_tv_chambolle(int_astro,
+                                                          weight=60.0)
     # test if the value range of output float data is within [0.0:1.0]
     assert denoised_int_astro.dtype == np.float
     assert np.max(denoised_int_astro) <= 1.0
@@ -143,5 +144,34 @@ def test_denoise_bilateral_3d():
     assert out1[30:45, 5:15].std() > out2[30:45, 5:15].std()
 
 
+def test_nl_means_denoising_2d():
+    img = np.zeros((40, 40))
+    img[10:-10, 10:-10] = 1.
+    img += 0.3*np.random.randn(*img.shape)
+    denoised = restoration.nl_means_denoising(img, 7, 5, 0.1)
+    # make sure noise is reduced
+    assert img.std() > denoised.std()
+
+
+def test_nl_means_denoising_2drgb():
+    # reduce image size because nl means is very slow
+    img = np.copy(astro[:50, :50])
+    # add some random noise
+    img += 0.5 * img.std() * np.random.random(img.shape)
+    img = np.clip(img, 0, 1)
+    denoised = restoration.nl_means_denoising(img, 7, 9, 0.08)
+    # make sure noise is reduced
+    assert img.std() > denoised.std()
+
+
+def test_nl_means_denoising_3d():
+    img = np.zeros((20, 20, 10))
+    img[5:-5, 5:-5, 3:-3] = 1.
+    img += 0.3*np.random.randn(*img.shape)
+    denoised = restoration.nl_means_denoising(img, 5, 4, 0.1)
+    # make sure noise is reduced
+    assert img.std() > denoised.std()
+
+
 if __name__ == "__main__":
     run_module_suite()