diff --git a/skimage/restoration/__init__.py b/skimage/restoration/__init__.py index b008c638..a5c0ce0c 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 +from .nl_means_denoising import nl_means_denoising, fast_nl_means_denoising __all__ = ['wiener', 'unsupervised_wiener', @@ -31,4 +31,5 @@ __all__ = ['wiener', 'denoise_tv_bregman', 'denoise_tv_chambolle', 'denoise_bilateral', - 'nl_means_denoising'] + 'nl_means_denoising', + 'fast_nl_means_denoising'] diff --git a/skimage/restoration/_nl_means_denoising.pyx b/skimage/restoration/_nl_means_denoising.pyx index bdf4fd6e..7a1b866f 100644 --- a/skimage/restoration/_nl_means_denoising.pyx +++ b/skimage/restoration/_nl_means_denoising.pyx @@ -251,8 +251,9 @@ def _nl_means_denoising_3d(image, int s=7, n_x, n_y, n_z = image.shape cdef int offset = s / 2 # padd the image so that boundaries are denoised as well - cdef DTYPE_t [:, :, ::1] padded = np.ascontiguousarray(util.pad(image, - offset, mode='reflect').astype(np.float32)) + cdef DTYPE_t [:, :, ::1] padded = np.ascontiguousarray(util.pad( + image.astype(np.float32), + offset, mode='reflect')) cdef DTYPE_t [:, :, ::1] result = padded.copy() h *= (np.max(padded) - np.min(padded)) cdef float A = ((s - 1.) / 4.) @@ -305,3 +306,80 @@ def _nl_means_denoising_3d(image, int s=7, new_value += weight * padded[x + i, y + j, z + k] result[x, y, z] = new_value / weight_sum return result[offset:-offset, offset:-offset, offset:-offset] + + +@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 + 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. + """ + + if s % 2 == 0: + s += 1 # odd value for symmetric patch + cdef int offset = s / 2 + cdef int pad_size = offset + d + cdef DTYPE_t [:, ::1] padded = np.ascontiguousarray(util.pad(image, + pad_size, mode='reflect').astype(np.float32)) + cdef DTYPE_t [:, ::1] result = padded.copy() + cdef DTYPE_t [:, ::1] weights = np.zeros_like(padded) + cdef DTYPE_t [:, ::1] integral = np.zeros_like(padded) + h *= (np.max(padded) - np.min(padded)) + cdef int n_x, n_y, t1, t2, x, y + cdef float weight, distance + cdef float h2 = h**2 + cdef float s2 = s**2. + n_x, n_y = image.shape + n_x += 2 * pad_size + n_y += 2 * pad_size + for t1 in range(-d, d + 1): + for t2 in range(-d, d + 1): + integral = np.zeros_like(padded) + for x in range(max(1, - t1), min(n_x, n_x - t1)): + integral[x, 0] = (padded[x, 0] - padded[x + t1, t2])**2 + \ + integral[x - 1, 0] + for y in range(max(1, - t2), min(n_y, n_y - t2)): + integral[0, y] = (padded[0, y] - padded[t1, y + t2])**2 + \ + integral[0, y - 1] + for x in range(max(1, - t1), min(n_x, n_x - t1)): + for y in range(max(1, - t2), min(n_y, n_y - t2)): + integral[x, y] = (padded[x, y] - + padded[x + t1, y + t2])**2 + \ + integral[x - 1, y] + integral[x, y - 1] \ + - integral[x - 1, y - 1] + for x in range(max(offset, offset - t1), + min(n_x - offset, n_x - offset - t1)): + for y in range(max(offset, offset - t2), + min(n_y - offset, n_y - offset - t2)): + distance = integral[x + offset, y + offset] + \ + integral[x - offset, y - offset] - \ + integral[x - offset, y + offset] - \ + integral[x + offset, y - offset] + distance /= s2 + weight = exp(- distance / h2) + weights[x, y] += weight + result[x, y] += weight * padded[x + t1, y + t2] + for x in range(offset, n_x - offset): + for y in range(offset, n_x - offset): + # I think there is no risk of division by zero + # except in padded zone + result[x, y] /= weights[x, y] + return result[pad_size: - pad_size, pad_size: - pad_size] + + diff --git a/skimage/restoration/nl_means_denoising.py b/skimage/restoration/nl_means_denoising.py index c106795e..79c67e04 100644 --- a/skimage/restoration/nl_means_denoising.py +++ b/skimage/restoration/nl_means_denoising.py @@ -1,6 +1,7 @@ import numpy as np from skimage.restoration._nl_means_denoising import _nl_means_denoising_2d, \ - _nl_means_denoising_2drgb, _nl_means_denoising_3d + _nl_means_denoising_2drgb, _nl_means_denoising_3d, \ + _fast_nl_means_denoising_2d def nl_means_denoising(image, patch_size=7, patch_distance=11, h=0.1): """ @@ -29,6 +30,10 @@ def nl_means_denoising(image, patch_size=7, patch_distance=11, h=0.1): result: ndarray denoised image, of same shape as `image`. + See Also + -------- + fast_nl_means_denoising + Notes ----- @@ -72,3 +77,81 @@ def nl_means_denoising(image, patch_size=7, patch_distance=11, h=0.1): else: raise ValueError("Non local means denoising is only possible for \ 2D grayscale and RGB images or 3-D grayscale images.") + + +def fast_nl_means_denoising(image, patch_size=7, patch_distance=11, h=0.1): + """ + Performs fast non-local means denoising on 2-D grayscale images. + + Parameters + ---------- + image: ndarray + input data to be denoised + + 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. + + Returns + ------- + + result: ndarray + denoised image, of same shape as `image`. + + See Also + -------- + 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. + + The complexity of the algorithm is + + 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. + + The image is padded using the `reflect` mode of `skimage.util.pad` + before denoising. + + References + ---------- + .. [1] 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 = fast_nl_means_denoising(a, 7, 5, 0.1) + """ + if image.ndim == 2: + return np.array(_fast_nl_means_denoising_2d(image, patch_size, + patch_distance, h)) + else: + raise ValueError("Fast non local means denoising is only possible for \ + 2D grayscale images.") diff --git a/skimage/restoration/tests/test_denoise.py b/skimage/restoration/tests/test_denoise.py index 77d38557..04fcd184 100644 --- a/skimage/restoration/tests/test_denoise.py +++ b/skimage/restoration/tests/test_denoise.py @@ -153,6 +153,15 @@ def test_nl_means_denoising_2d(): assert img.std() > denoised.std() +def test_fast_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.fast_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])