diff --git a/bento.info b/bento.info
index bfc607a8..9987ebca 100644
--- a/bento.info
+++ b/bento.info
@@ -149,6 +149,9 @@ Library:
     Extension: skimage.restoration._denoise_cy
         Sources:
             skimage/restoration/_denoise_cy.pyx
+    Extension: skimage.restoration._nl_means_denoising
+        Sources:
+            skimage/restoration/_nl_means_denoising.pyx
     Extension: skimage.feature._hessian_det_appx
         Sources:
             skimage/exposure/_hessian_det_appx.pyx
diff --git a/doc/examples/plot_nonlocal_means.py b/doc/examples/plot_nonlocal_means.py
new file mode 100644
index 00000000..8349caa9
--- /dev/null
+++ b/doc/examples/plot_nonlocal_means.py
@@ -0,0 +1,41 @@
+"""
+=================================================
+Non-local means denoising for preserving textures
+=================================================
+
+In this example, we denoise a detail of the astronaut image using the non-local
+means filter. The non-local means algorithm replaces the value of a pixel by an
+average of a selection of other pixels values: small patches centered on the
+other pixels are compared to the patch centered on the pixel of interest, and
+the average is performed only for pixels that have patches close to the current
+patch. As a result, this algorithm can restore well textures, that would be
+blurred by other denoising algoritm.
+"""
+import numpy as np
+import matplotlib.pyplot as plt
+
+from skimage import data, img_as_float
+from skimage.restoration import nl_means_denoising
+
+
+astro = img_as_float(data.astronaut())
+astro = astro[30:180, 150:300]
+
+noisy = astro + 0.3 * np.random.random(astro.shape)
+noisy = np.clip(noisy, 0, 1)
+
+denoise = nl_means_denoising(noisy, 7, 9, 0.08)
+
+fig, ax = plt.subplots(ncols=2, figsize=(8, 4))
+
+ax[0].imshow(noisy)
+ax[0].axis('off')
+ax[0].set_title('noisy')
+ax[1].imshow(denoise)
+ax[1].axis('off')
+ax[1].set_title('non-local means')
+
+fig.subplots_adjust(wspace=0.02, hspace=0.2,
+                    top=0.9, bottom=0.05, left=0, right=1)
+
+plt.show()
diff --git a/skimage/restoration/__init__.py b/skimage/restoration/__init__.py
index 2b593ccf..66beecfa 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 .non_local_means 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/restoration/_nl_means_denoising.pyx b/skimage/restoration/_nl_means_denoising.pyx
new file mode 100644
index 00000000..089914fc
--- /dev/null
+++ b/skimage/restoration/_nl_means_denoising.pyx
@@ -0,0 +1,712 @@
+import numpy as np
+from skimage import util
+cimport numpy as np
+cimport cython
+from libc.math cimport exp
+
+ctypedef np.float32_t IMGDTYPE
+
+cdef float DISTANCE_CUTOFF = 5.
+
+@cython.boundscheck(False)
+cdef inline float patch_distance_2d(IMGDTYPE [:, :] p1,
+                                    IMGDTYPE [:, :] p2,
+                                    IMGDTYPE [:, ::] w, int s):
+    """
+    Compute a Gaussian distance between two image patches.
+
+    Parameters
+    ----------
+    p1 : 2-D array_like
+        First patch.
+    p2 : 2-D array_like
+        Second patch.
+    w : 2-D array_like
+        Array of weigths for the different pixels of the patches.
+    s : int
+        Linear size of the patches.
+
+    Returns
+    -------
+    distance : float
+        Gaussian distance between the two patches
+
+    Notes
+    -----
+    The returned distance is given by
+
+    .. math::  \exp( -w (p1 - p2)^2)
+    """
+    cdef int i, j
+    cdef int center = s / 2
+    # Check if central pixel is too different in the 2 patches
+    cdef float tmp_diff = p1[center, center] - p2[center, center]
+    cdef float init = w[center, center] * tmp_diff * tmp_diff
+    if init > 1:
+        return 0.
+    cdef float distance = 0
+    for i in range(s):
+        # exp of large negative numbers will be 0, so we'd better stop
+        if distance > DISTANCE_CUTOFF:
+            return 0.
+        for j in range(s):
+            tmp_diff = p1[i, j] - p2[i, j]
+            distance += (w[i, j] * tmp_diff * tmp_diff)
+    distance = exp(-distance)
+    return distance
+
+
+@cython.boundscheck(False)
+cdef inline float patch_distance_2drgb(IMGDTYPE [:, :, :] p1,
+                                       IMGDTYPE [:, :, :] p2,
+                                       IMGDTYPE [:, ::] w, int s):
+    """
+    Compute a Gaussian distance between two image patches.
+
+    Parameters
+    ----------
+    p1 : 3-D array_like
+        First patch, 2D image with last dimension corresponding to channels.
+    p2 : 3-D array_like
+        Second patch, 2D image with last dimension corresponding to channels.
+    w : 2-D array_like
+        Array of weights for the different pixels of the patches.
+    s : int
+        Linear size of the patches.
+
+    Returns
+    -------
+    distance : float
+        Gaussian distance between the two patches
+
+    Notes
+    -----
+    The returned distance is given by
+
+    .. math::  \exp( -w (p1 - p2)^2)
+    """
+    cdef int i, j
+    cdef int center = s / 2
+    cdef int color
+    cdef float tmp_diff = 0
+    cdef float distance = 0
+    for i in range(s):
+        # exp of large negative numbers will be 0, so we'd better stop
+        if distance > DISTANCE_CUTOFF:
+            return 0.
+        for j in range(s):
+            for color in range(3):
+                tmp_diff = p1[i, j, color] - p2[i, j, color]
+                distance += w[i, j] * tmp_diff * tmp_diff
+    distance = exp(-distance)
+    return distance
+
+
+@cython.boundscheck(False)
+cdef inline float patch_distance_3d(IMGDTYPE [:, :, :] p1,
+                                    IMGDTYPE [:, :, :] p2,
+                                    IMGDTYPE [:, :, ::] w, int s):
+    """
+    Compute a Gaussian distance between two image patches.
+
+    Parameters
+    ----------
+    p1 : 3-D array_like
+        First patch.
+    p2 : 3-D array_like
+        Second patch.
+    w : 3-D array_like
+        Array of weights for the different pixels of the patches.
+    s : int
+        Linear size of the patches.
+
+    Returns
+    -------
+    distance : float
+        Gaussian distance between the two patches
+
+    Notes
+    -----
+    The returned distance is given by
+
+    .. math::  \exp( -w (p1 - p2)^2)
+    """
+    cdef int i, j, k
+    cdef float distance = 0
+    cdef float tmp_diff
+    for i in range(s):
+        # exp of large negative numbers will be 0, so we'd better stop
+        if distance > DISTANCE_CUTOFF:
+            return 0.
+        for j in range(s):
+            for k in range(s):
+                tmp_diff = p1[i, j, k] - p2[i, j, k]
+                distance += w[i, j, k] * tmp_diff * tmp_diff
+    distance = exp(-distance)
+    return distance
+
+
+@cython.cdivision(True)
+@cython.boundscheck(False)
+def _nl_means_denoising_2d(image, int s=7, int d=13, float h=0.1):
+    """
+    Perform non-local means denoising on 2-D RGB image
+
+    Parameters
+    ----------
+    image : ndarray
+        Input RGB image to be denoised
+    s : int, optional
+        Size of patches used for denoising
+    d : int, optional
+        Maximal distance in pixels where to search patches used for denoising
+    h : float, optional
+        Cut-off distance (in gray levels). The higher h, the more permissive
+        one is in accepting patches.
+
+    Returns
+    -------
+    result : ndarray
+        Denoised image, of same shape as input image.
+    """
+    if s % 2 == 0:
+        s += 1  # odd value for symmetric patch
+    cdef int n_row, n_col, n_ch
+    n_row, n_col, n_ch = image.shape
+    cdef int offset = s / 2
+    cdef int row, col, i, j, color
+    cdef int row_start, row_end, col_start, col_end
+    cdef int row_start_i, row_end_i, col_start_j, col_end_j
+    cdef IMGDTYPE [::1] new_values = np.zeros(n_ch).astype(np.float32)
+    cdef IMGDTYPE [:, :, ::1] padded = np.ascontiguousarray(util.pad(image,
+                       ((offset, offset), (offset, offset), (0, 0)),
+                        mode='reflect').astype(np.float32))
+    cdef IMGDTYPE [:, :, ::1] result = padded.copy()
+    cdef float A = ((s - 1.) / 4.)
+    cdef float new_value
+    cdef float weight_sum, weight
+    xg_row, xg_col = np.mgrid[-offset:offset + 1, -offset:offset + 1]
+    cdef IMGDTYPE [:, ::1] w = np.ascontiguousarray(np.exp(
+                             -(xg_row ** 2 + xg_col ** 2) / (2 * A ** 2)).
+                             astype(np.float32))
+    cdef float distance
+    w = 1. / (n_ch * np.sum(w) * h ** 2) * w
+
+    # Coordinates of central pixel
+    # Iterate over rows, taking padding into account
+    for row in range(offset, n_row + offset):
+        row_start = row - offset
+        row_end = row + offset + 1
+        # Iterate over columns, taking padding into account
+        for col in range(offset, n_col + offset):
+            # Initialize per-channel bins
+            for color in range(n_ch):
+                new_values[color] = 0
+            # Reset weights for each local region
+            weight_sum = 0
+            col_start = col - offset
+            col_end = col + offset + 1
+
+            # Iterate over local 2d patch for each pixel
+            # First rows
+            for i in range(max(-d, offset - row),
+                           min(d + 1, n_row + offset - row)):
+                row_start_i = row_start + i
+                row_end_i = row_end + i
+                # Local patch columns
+                for j in range(max(-d, offset - col),
+                               min(d + 1, n_col + offset - col)):
+                    col_start_j = col_start + j
+                    col_end_j = col_end + j
+                    # Shortcut for grayscale, else assume RGB
+                    if n_ch == 1:
+                        weight = patch_distance_2d(
+                                 padded[row_start:row_end,
+                                        col_start:col_end, 0],
+                                 padded[row_start_i:row_end_i,
+                                        col_start_j:col_end_j, 0],
+                                 w, s)
+                    else:
+                        weight = patch_distance_2drgb(
+                                 padded[row_start:row_end,
+                                        col_start:col_end, :],
+                                 padded[row_start_i:row_end_i,
+                                        col_start_j:col_end_j, :],
+                                        w, s)
+
+                    # Collect results in weight sum
+                    weight_sum += weight
+                    # Apply to each channel multiplicatively
+                    for color in range(n_ch):
+                        new_values[color] += weight * padded[row + i,
+                                                             col + j, color]
+
+            # Normalize the result
+            for color in range(n_ch):
+                result[row, col, color] = new_values[color] / weight_sum
+
+    # Return cropped result, undoing padding
+    return result[offset:-offset, offset:-offset]
+
+
+@cython.cdivision(True)
+@cython.boundscheck(False)
+def _nl_means_denoising_3d(image, int s=7,
+            int d=13, float h=0.1):
+    """
+    Perform non-local means denoising on 3-D array
+
+    Parameters
+    ----------
+    image : ndarray
+        Input data to be denoised.
+    s : int, optional
+        Size of patches used for denoising.
+    d : int, optional
+        Maximal distance in pixels where to search patches used for denoising.
+    h : float, optional
+        Cut-off distance (in gray levels).
+
+    Returns
+    -------
+    result : ndarray
+        Denoised image, of same shape as input image.
+    """
+    if s % 2 == 0:
+        s += 1  # odd value for symmetric patch
+    cdef int n_pln, n_row, n_col
+    n_pln, n_row, n_col = image.shape
+    cdef int offset = s / 2
+    # padd the image so that boundaries are denoised as well
+    cdef IMGDTYPE [:, :, ::1] padded = np.ascontiguousarray(util.pad(
+                                        image.astype(np.float32),
+                                        offset, mode='reflect'))
+    cdef IMGDTYPE [:, :, ::1] result = padded.copy()
+    cdef float A = ((s - 1.) / 4.)
+    cdef float new_value
+    cdef float weight_sum, weight
+    xg_pln, xg_row, xg_col = np.mgrid[-offset: offset + 1,
+                                      -offset: offset + 1,
+                                      -offset: offset + 1]
+    cdef IMGDTYPE [:, :, ::1] w = np.ascontiguousarray(np.exp(
+                            -(xg_pln ** 2 + xg_row ** 2 + xg_col ** 2) /
+                             (2 * A ** 2)).astype(np.float32))
+    cdef float distance
+    cdef int pln, row, col, i, j, k
+    cdef int pln_start, pln_end, row_start, row_end, col_start, col_end
+    cdef int pln_start_i, pln_end_i, row_start_j, row_end_j, \
+             col_start_k, col_end_k
+    w = 1. / (np.sum(w) * h ** 2) * w
+
+    # Coordinates of central pixel
+    # Iterate over planes, taking padding into account
+    for pln in range(offset, n_pln + offset):
+        pln_start = pln - offset
+        pln_end = pln + offset + 1
+        # Iterate over rows, taking padding into account
+        for row in range(offset, n_row + offset):
+            row_start = row - offset
+            row_end = row + offset + 1
+            # Iterate over columns, taking padding into account
+            for col in range(offset, n_col + offset):
+                col_start = col - offset
+                col_end = col + offset + 1
+                new_value = 0
+                weight_sum = 0
+
+                # Iterate over local 3d patch for each pixel
+                # First planes
+                for i in range(max(-d, offset - pln),
+                               min(d + 1, n_pln + offset - pln)):
+                    pln_start_i = pln_start + i
+                    pln_end_i = pln_end + i
+                    # Rows
+                    for j in range(max(-d, offset - row),
+                                   min(d + 1, n_row + offset - row)):
+                        row_start_j = row_start + j
+                        row_end_j = row_end + j
+                        # Columns
+                        for k in range(max(-d, offset - col),
+                                       min(d + 1, n_col + offset - col)):
+                            col_start_k = col_start + k
+                            col_end_k = col_end + k
+                            weight = patch_distance_3d(
+                                    padded[pln_start:pln_end,
+                                           row_start:row_end,
+                                           col_start:col_end],
+                                    padded[pln_start_i:pln_end_i,
+                                           row_start_j:row_end_j,
+                                           col_start_k:col_end_k],
+                                    w, s)
+                            # Collect results in weight sum
+                            weight_sum += weight
+                            new_value += weight * padded[pln + i,
+                                                         row + j, col + k]
+
+                # Normalize the result
+                result[pln, row, col] = new_value / weight_sum
+
+    # Return cropped result, undoing padding
+    return result[offset:-offset, offset:-offset, offset:-offset]
+
+#-------------- Accelerated algorithm of Froment 2015 ------------------
+
+
+@cython.cdivision(True)
+@cython.boundscheck(False)
+cdef inline float _integral_to_distance_2d(IMGDTYPE [:, ::] integral,
+                        int row, int col, int offset, float h2s2):
+    """
+    References
+    ----------
+    Jacques Froment. Parameter-Free Fast Pixelwise Non-Local Means
+    Denoising. Image Processing On Line, 2014, vol. 4, p. 300-326.
+
+    Used in _fast_nl_means_denoising_2d
+    """
+    cdef float distance
+    distance =  integral[row + offset, col + offset] + \
+                integral[row - offset, col - offset] - \
+                integral[row - offset, col + offset] - \
+                integral[row + offset, col - offset]
+    distance /= h2s2
+    return distance
+
+
+@cython.cdivision(True)
+@cython.boundscheck(False)
+cdef inline float _integral_to_distance_3d(IMGDTYPE [:, :, ::] integral,
+                    int pln, int row, int col, int offset,
+                    float s_cube_h_square):
+    """
+    References
+    ----------
+    Jacques Froment. Parameter-Free Fast Pixelwise Non-Local Means
+    Denoising. Image Processing On Line, 2014, vol. 4, p. 300-326.
+
+    Used in _fast_nl_means_denoising_3d
+    """
+    cdef float distance
+    distance = (integral[pln + offset, row + offset, col + offset] -
+                integral[pln - offset, row - offset, col - offset] +
+                integral[pln - offset, row - offset, col + offset] +
+                integral[pln - offset, row + offset, col - offset] +
+                integral[pln + offset, row - offset, col - offset] -
+                integral[pln - offset, row + offset, col + offset] -
+                integral[pln + offset, row - offset, col + offset] -
+                integral[pln + offset, row + offset, col - offset])
+    distance /= s_cube_h_square
+    return distance
+
+
+@cython.cdivision(True)
+@cython.boundscheck(False)
+cdef inline _integral_image_2d(IMGDTYPE [:, :, ::] padded,
+                               IMGDTYPE [:, ::] integral, int t_row,
+                               int t_col, int n_row, int n_col, int n_ch):
+    """
+    Computes the integral of the squared difference between an image ``padded``
+    and the same image shifted by ``(t_row, t_col)``.
+
+    Parameters
+    ----------
+    padded : ndarray of shape (n_row, n_col, n_ch)
+        Image of interest.
+    integral : ndarray
+        Output of the function. The array is filled with integral values.
+        ``integral`` should have the same shape as ``padded``.
+    t_row : int
+        Shift along the row axis.
+    t_col : int
+        Shift along the column axis.
+    n_row : int
+    n_col : int
+    n_ch : int
+
+    Notes
+    -----
+
+    The integral computation could be performed using
+    ``transform.integral_image``, but this helper function saves memory
+    by avoiding copies of ``padded``.
+    """
+    cdef int row, col
+    cdef float distance
+    for row in range(max(1, -t_row), min(n_row, n_row - t_row)):
+        for col in range(max(1, -t_col), min(n_col, n_col - t_col)):
+            if n_ch == 1:
+                distance = (padded[row, col, 0] -
+                            padded[row + t_row, col + t_col, 0])**2
+            else:
+                distance = ((padded[row, col, 0] -
+                             padded[row + t_row, col + t_col, 0])**2 +
+                            (padded[row, col, 1] -
+                             padded[row + t_row, col + t_col, 1])**2 +
+                            (padded[row, col, 2] -
+                             padded[row + t_row, col + t_col, 2])**2)
+            integral[row, col] = distance + \
+                                 integral[row - 1, col] + \
+                                 integral[row, col - 1] - \
+                                 integral[row - 1, col - 1]
+
+@cython.cdivision(True)
+@cython.boundscheck(False)
+cdef inline _integral_image_3d(IMGDTYPE [:, :, ::] padded,
+                               IMGDTYPE [:, :, ::] integral, int t_pln,
+                               int t_row, int t_col, int n_pln, int n_row,
+                               int n_col):
+    """
+    Computes the integral of the squared difference between an image ``padded``
+    and the same image shifted by ``(t_pln, t_row, t_col)``.
+
+    Parameters
+    ----------
+    padded : ndarray of shape (n_pln, n_row, n_col)
+        Image of interest.
+    integral : ndarray
+        Output of the function. The array is filled with integral values.
+        ``integral`` should have the same shape as ``padded``.
+    t_pln : int
+        Shift along the plane axis.
+    t_row : int
+        Shift along the row axis.
+    t_col : int
+        Shift along the column axis.
+    n_pln : int
+    n_row : int
+    n_col : int
+
+    Notes
+    -----
+
+    The integral computation could be performed using
+    ``transform.integral_image``, but this helper function saves memory
+    by avoiding copies of ``padded``.
+    """
+    cdef int pln, row, col
+    cdef float distance
+    for pln in range(max(1, -t_pln), min(n_pln, n_pln - t_pln)):
+        for row in range(max(1, -t_row), min(n_row, n_row - t_row)):
+            for col in range(max(1, -t_col), min(n_col, n_col - t_col)):
+                integral[pln, row, col] = \
+                    ((padded[pln, row, col] -
+                      padded[pln + t_pln, row + t_row, col + t_col])**2 +
+                    integral[pln - 1, row, col] +
+                    integral[pln, row - 1, col] +
+                    integral[pln, row, col - 1] +
+                    integral[pln - 1, row - 1, col - 1] -
+                    integral[pln - 1, row - 1, col] -
+                    integral[pln, row - 1, col - 1] -
+                    integral[pln - 1, row, col - 1])
+
+
+@cython.cdivision(True)
+@cython.boundscheck(False)
+def _fast_nl_means_denoising_2d(image, int s=7, int d=13, float h=0.1):
+    """
+    Perform fast non-local means denoising on 2-D array, with the outer
+    loop on patch shifts in order to reduce the number of operations.
+
+    Parameters
+    ----------
+    image : ndarray
+        2-D input data to be denoised, grayscale or RGB.
+    s : int, optional
+        Size of patches used for denoising.
+    d : int, optional
+        Maximal distance in pixels where to search patches used for denoising.
+    h : float, optional
+        Cut-off distance (in gray levels). The higher h, the more permissive
+        one is in accepting patches.
+
+    Returns
+    -------
+    result : ndarray
+        Denoised image, of same shape as input image.
+    """
+    if s % 2 == 0:
+        s += 1  # odd value for symmetric patch
+    cdef int offset = s / 2
+    # Image padding: we need to account for patch size, possible shift,
+    # + 1 for the boundary effects in finite differences
+    cdef int pad_size = offset + d + 1
+    cdef IMGDTYPE [:, :, ::1] padded = np.ascontiguousarray(util.pad(image,
+                          ((pad_size, pad_size), (pad_size, pad_size), (0, 0)),
+                          mode='reflect').astype(np.float32))
+    cdef IMGDTYPE [:, :, ::1] result = np.zeros_like(padded)
+    cdef IMGDTYPE [:, ::1] weights = np.zeros_like(padded[..., 0], order='C')
+    cdef IMGDTYPE [:, ::1] integral = np.zeros_like(padded[..., 0], order='C')
+    cdef int n_row, n_col, n_ch, t_row, t_col, row, col
+    cdef float weight, distance
+    cdef float alpha
+    cdef float h2 = h ** 2.
+    cdef float s2 = s ** 2.
+    n_row, n_col, n_ch = image.shape
+    cdef float h2s2 = n_ch * h2 * s2
+    n_row += 2 * pad_size
+    n_col += 2 * pad_size
+
+    # Outer loops on patch shifts
+    # With t2 >= 0, reference patch is always on the left of test patch
+    # Iterate over shifts along the row axis
+    for t_row in range(-d, d + 1):
+        # Iterate over shifts along the column axis
+        for t_col in range(0, d + 1):
+            # alpha is to account for patches on the same column
+            # distance is computed twice in this case
+            if t_col == 0 and t_row is not 0:
+                alpha = 0.5
+            else:
+                alpha = 1.
+            # Compute integral image of the squared difference between
+            # padded and the same image shifted by (t_row, t_col)
+            integral = np.zeros_like(padded[..., 0], order='C')
+            _integral_image_2d(padded, integral, t_row, t_col,
+                               n_row, n_col, n_ch)
+
+            # Inner loops on pixel coordinates
+            # Iterate over rows, taking offset and shift into account
+            for row in range(max(offset, offset - t_row),
+                             min(n_row - offset, n_row - offset - t_row)):
+                # Iterate over columns, taking offset and shift into account
+                for col in range(max(offset, offset - t_col),
+                                 min(n_col - offset, n_col - offset - t_col)):
+                    # Compute squared distance between shifted patches
+                    distance = _integral_to_distance_2d(integral, row, col,
+                                                     offset, h2s2)
+                    # exp of large negative numbers is close to zero
+                    if distance > DISTANCE_CUTOFF:
+                        continue
+                    weight = alpha * exp(-distance)
+                    # Accumulate weights corresponding to different shifts
+                    weights[row, col] += weight
+                    weights[row + t_row, col + t_col] += weight
+                    # Iterate over channels
+                    for ch in range(n_ch):
+                        result[row, col, ch] += weight * \
+                                    padded[row + t_row, col + t_col, ch]
+                        result[row + t_row, col + t_col, ch] += \
+                                        weight * padded[row, col, ch]
+
+    # Normalize pixel values using sum of weights of contributing patches
+    for row in range(offset, n_row - offset):
+        for col in range(offset, n_col - offset):
+            for channel in range(n_ch):
+                # No risk of division by zero, since the contribution
+                # of a null shift is strictly positive
+                result[row, col, channel] /= weights[row, col]
+
+    # Return cropped result, undoing padding
+    return result[pad_size:-pad_size, pad_size:-pad_size]
+
+
+@cython.cdivision(True)
+@cython.boundscheck(False)
+def _fast_nl_means_denoising_3d(image, int s=5, int d=7, float h=0.1):
+    """
+    Perform fast non-local means denoising on 3-D array, with the outer
+    loop on patch shifts in order to reduce the number of operations.
+
+    Parameters
+    ----------
+    image : ndarray
+        3-D input data to be denoised.
+    s : int, optional
+        Size of patches used for denoising.
+    d : int, optional
+        Maximal distance in pixels where to search patches used for denoising.
+    h : float, optional
+        cut-off distance (in gray levels). The higher h, the more permissive
+        one is in accepting patches.
+
+    Returns
+    -------
+    result : ndarray
+        Denoised image, of same shape as input image.
+    """
+    if s % 2 == 0:
+        s += 1  # odd value for symmetric patch
+    cdef int offset = s / 2
+    # Image padding: we need to account for patch size, possible shift,
+    # + 1 for the boundary effects in finite differences
+    cdef int pad_size = offset + d + 1
+    cdef IMGDTYPE [:, :, ::1] padded = np.ascontiguousarray(util.pad(image,
+                                pad_size, mode='reflect').astype(np.float32))
+    cdef IMGDTYPE [:, :, ::1] result = np.zeros_like(padded)
+    cdef IMGDTYPE [:, :, ::1] weights = np.zeros_like(padded)
+    cdef IMGDTYPE [:, :, ::1] integral = np.zeros_like(padded)
+    cdef int n_pln, n_row, n_col, t_pln, t_row, t_col, \
+             pln, row, col
+    cdef int pln_dist_min, pln_dist_max, row_dist_min, row_dist_max, \
+             col_dist_min, col_dist_max
+    cdef float weight, distance
+    cdef float alpha
+    cdef float h_square = h ** 2.
+    cdef float s_cube = s ** 3.
+    cdef float s_cube_h_square = h_square * s_cube
+    n_pln, n_row, n_col = image.shape
+    n_pln += 2 * pad_size
+    n_row += 2 * pad_size
+    n_col += 2 * pad_size
+
+    # Outer loops on patch shifts
+    # With t2 >= 0, reference patch is always on the left of test patch
+    # Iterate over shifts along the plane axis
+    for t_pln in range(-d, d + 1):
+        pln_dist_min = max(offset, offset - t_pln)
+        pln_dist_max = min(n_pln - offset, n_pln - offset - t_pln)
+        # Iterate over shifts along the row axis
+        for t_row in range(-d, d + 1):
+            row_dist_min = max(offset, offset - t_row)
+            row_dist_max = min(n_row - offset, n_row - offset - t_row)
+            # Iterate over shifts along the column axis
+            for t_col in range(0, d + 1):
+                col_dist_min = max(offset, offset - t_col)
+                col_dist_max = min(n_col - offset, n_col - offset - t_col)
+                # alpha is to account for patches on the same column
+                # distance is computed twice in this case
+                if t_col == 0 and (t_pln is not 0 or t_row is not 0):
+                    alpha = 0.5
+                else:
+                    alpha = 1.
+                # Compute integral image of the squared difference between
+                # padded and the same image shifted by (t_pln, t_row, t_col)
+                integral = np.zeros_like(padded)
+                _integral_image_3d(padded, integral, t_pln, t_row, t_col,
+                                   n_pln, n_row, n_col)
+
+                # Inner loops on pixel coordinates
+                # Iterate over planes, taking offset and shift into account
+                for pln in range(pln_dist_min, pln_dist_max):
+                    # Iterate over rows, taking offset and shift into account
+                    for row in range(row_dist_min, row_dist_max):
+                        # Iterate over columns
+                        for col in range(col_dist_min, col_dist_max):
+                            # Compute squared distance between shifted patches
+                            distance = _integral_to_distance_3d(integral,
+                                        pln, row, col, offset, s_cube_h_square)
+                            # exp of large negative numbers is close to zero
+                            if distance > DISTANCE_CUTOFF:
+                                continue
+                            weight = alpha * exp(-distance)
+                            # Accumulate weights for the different shifts
+                            weights[pln, row, col] += weight
+                            weights[pln + t_pln, row + t_row,
+                                                 col + t_col] += weight
+                            result[pln, row, col] += weight * \
+                                    padded[pln + t_pln, row + t_row,
+                                                        col + t_col]
+                            result[pln + t_pln, row + t_row,
+                                                col + t_col] += weight * \
+                                                  padded[pln, row, col]
+
+    # Normalize pixel values using sum of weights of contributing patches
+    for pln in range(offset, n_pln - offset):
+        for row in range(offset, n_row - offset):
+            for col in range(offset, n_col - offset):
+                # No risk of division by zero, since the contribution
+                # of a null shift is strictly positive
+                result[pln, row, col] /= weights[pln, row, col]
+
+    # Return cropped result, undoing padding
+    return result[pad_size:-pad_size, pad_size:-pad_size, pad_size:-pad_size]
diff --git a/skimage/restoration/non_local_means.py b/skimage/restoration/non_local_means.py
new file mode 100644
index 00000000..2ba98a11
--- /dev/null
+++ b/skimage/restoration/non_local_means.py
@@ -0,0 +1,120 @@
+import numpy as np
+from skimage.restoration._nl_means_denoising import _nl_means_denoising_2d, \
+                    _nl_means_denoising_3d, \
+                    _fast_nl_means_denoising_2d, _fast_nl_means_denoising_3d
+
+def nl_means_denoising(image, patch_size=7, patch_distance=11, h=0.1,
+                       multichannel=True, fast_mode=True):
+    """
+    Perform non-local means denoising on 2-D or 3-D grayscale images, and
+    2-D RGB images.
+
+    Parameters
+    ----------
+    image : 2D or 3D ndarray
+        Input image to be denoised, which can be 2D or 3D, and grayscale
+        or RGB (for 2D images only, see ``multichannel`` parameter).
+    patch_size : int, optional
+        Size of patches used for denoising.
+    patch_distance : int, optional
+        Maximal distance in pixels where to search patches used for denoising.
+    h : float, optional
+        Cut-off distance (in gray levels). The higher h, the more permissive
+        one is in accepting patches. A higher h results in a smoother image,
+        at the expense of blurring features. For a Gaussian noise of standard
+        deviation sigma, a rule of thumb is to choose the value of h to be
+        sigma of slightly less.
+    multichannel : bool, optional
+        Whether the last axis of the image is to be interpreted as multiple
+        channels or another spatial dimension. Set to ``False`` for 3-D images.
+    fast_mode : bool, optional
+        If True (default value), a fast version of the non-local means
+        algorithm is used. If False, the original version of non-local means is
+        used. See the Notes section for more details about the algorithms.
+
+    Returns
+    -------
+
+    result : ndarray
+        Denoised image, of same shape as `image`.
+
+    See Also
+    --------
+    fast_nl_means_denoising
+
+    Notes
+    -----
+
+    The non-local means algorithm is well suited for denoising images with
+    specific textures. The principle of the algorithm is to average the value
+    of a given pixel with values of other pixels in a limited neighbourhood,
+    provided that the *patches* centered on the other pixels are similar enough
+    to the patch centered on the pixel of interest.
+
+    In the original version of the algorithm [1]_, corresponding to
+    ``fast=False``, the computational complexity is
+
+    image.size * patch_size ** image.ndim * patch_distance ** image.ndim
+
+    Hence, changing the size of patches or their maximal distance has a
+    strong effect on computing times, especially for 3-D images.
+
+    However, the default behavior corresponds to ``fast_mode=True``, for which
+    another version of non-local means [2]_ is used, corresponding to a
+    complexity of
+
+    image.size * patch_distance ** image.ndim
+
+    The computing time depends only weakly on the patch size, thanks to the
+    computation of the integral of patches distances for a given shift, that
+    reduces the number of operations [1]_. Therefore, this algorithm executes
+    faster than `nl_means_denoising`, at the expense of using twice as much
+    memory.
+
+    Compared to the classic non-local means algorithm implemented in
+    `nl_means_denoising`, all pixels of a patch contribute to the distance to
+    another patch with the same weight, no matter their distance to the center
+    of the patch. This coarser computation of the distance can result in a
+    slightly poorer denoising performance. Moreover, for small images (images
+    with a linear size that is only a few times the patch size), the classic
+    algorithm can be faster due to boundary effects.
+
+    The image is padded using the `reflect` mode of `skimage.util.pad`
+    before denoising.
+
+    References
+    ----------
+    .. [1] Buades, A., Coll, B., & Morel, J. M. (2005, June). A non-local
+        algorithm for image denoising. In CVPR 2005, Vol. 2, pp. 60-65, IEEE.
+
+    .. [2] Jacques Froment. Parameter-Free Fast Pixelwise Non-Local Means
+           Denoising. Image Processing On Line, 2014, vol. 4, p. 300-326.
+
+    Examples
+    --------
+    >>> a = np.zeros((40, 40))
+    >>> a[10:-10, 10:-10] = 1.
+    >>> a += 0.3*np.random.randn(*a.shape)
+    >>> denoised_a = nl_means_denoising(a, 7, 5, 0.1)
+    """
+    if image.ndim == 2:
+        image = image[..., np.newaxis]
+        multichannel = True
+    if image.ndim != 3:
+        raise NotImplementedError("Non-local means denoising is only \
+        implemented for 2D grayscale and RGB images or 3-D grayscale images.")
+    if multichannel:  # 2-D images
+        if fast_mode:
+            return np.squeeze(np.array(_fast_nl_means_denoising_2d(image,
+                                       patch_size, patch_distance, h)))
+        else:
+            return np.squeeze(np.array(_nl_means_denoising_2d(image,
+                                       patch_size, patch_distance, h)))
+    else:  # 3-D grayscale
+        if fast_mode:
+            return np.array(_fast_nl_means_denoising_3d(image, s=patch_size,
+                                              d=patch_distance, h=h))
+        else:
+            return np.array(_nl_means_denoising_3d(image, patch_size,
+                                patch_distance, h))
+
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..e84e1c49 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,76 @@ 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.2, fast_mode=True)
+    # make sure noise is reduced
+    assert img.std() > denoised.std()
+    denoised = restoration.nl_means_denoising(img, 7, 5, 0.2, fast_mode=False)
+    # 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.3, fast_mode=True)
+    # make sure noise is reduced
+    assert img.std() > denoised.std()
+    denoised = restoration.nl_means_denoising(img, 7, 9, 0.3, fast_mode=False)
+    # 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.2, fast_mode=True,
+                                              multichannel=False)
+    # make sure noise is reduced
+    assert img.std() > denoised.std()
+    denoised = restoration.nl_means_denoising(img, 5, 4, 0.2, fast_mode=False,
+                                              multichannel=False)
+    # make sure noise is reduced
+    assert img.std() > denoised.std()
+
+
+def test_nl_means_denoising_multichannel():
+    img = np.zeros((21, 20, 10))
+    img[10, 9:11, 2:-2] = 1.
+    img += 0.3*np.random.randn(*img.shape)
+    denoised_wrong_multichannel = restoration.nl_means_denoising(img,
+                    5, 4, 0.1, fast_mode=True, multichannel=True)
+    denoised_ok_multichannel = restoration.nl_means_denoising(img,
+                    5, 4, 0.1, fast_mode=True, multichannel=False)
+    snr_wrong = 10 * np.log10(1. /
+                            ((denoised_wrong_multichannel - img)**2).mean())
+    snr_ok = 10 * np.log10(1. /
+                            ((denoised_ok_multichannel - img)**2).mean())
+    assert snr_ok > snr_wrong
+
+
+def test_nl_means_denoising_wrong_dimension():
+    img = np.zeros((5, 5, 5, 5))
+    assert_raises(NotImplementedError, restoration.nl_means_denoising, img)
+
+
+def test_no_denoising_for_small_h():
+    img = np.zeros((40, 40))
+    img[10:-10, 10:-10] = 1.
+    img += 0.3*np.random.randn(*img.shape)
+    # very small h should result in no averaging with other patches
+    denoised = restoration.nl_means_denoising(img, 7, 5, 0.01, fast_mode=True)
+    assert np.allclose(denoised, img)
+    denoised = restoration.nl_means_denoising(img, 7, 5, 0.01, fast_mode=False)
+    assert np.allclose(denoised, img)
+
+
 if __name__ == "__main__":
     run_module_suite()