From dd9030d44c88f55f15901c3a38012bc7821a597b Mon Sep 17 00:00:00 2001 From: emmanuelle Date: Sun, 18 Jan 2015 22:18:29 +0100 Subject: [PATCH] Implemented fast algorithm also for 3-D and 2D-RGB images. Changed API so that there is only one function for fast and classic algorithms. --- doc/examples/plot_nonlocal_means.py | 10 +- skimage/restoration/__init__.py | 5 +- skimage/restoration/_nl_means_denoising.pyx | 220 +++++++++++++++++++- skimage/restoration/non_local_means.py | 136 +++++------- skimage/restoration/tests/test_denoise.py | 25 ++- 5 files changed, 289 insertions(+), 107 deletions(-) diff --git a/doc/examples/plot_nonlocal_means.py b/doc/examples/plot_nonlocal_means.py index c6879be4..8349caa9 100644 --- a/doc/examples/plot_nonlocal_means.py +++ b/doc/examples/plot_nonlocal_means.py @@ -3,7 +3,7 @@ Non-local means denoising for preserving textures ================================================= -In this example, we denoise a detail of the Lena image using the non-local +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 @@ -18,13 +18,13 @@ from skimage import data, img_as_float from skimage.restoration import nl_means_denoising -lena = img_as_float(data.lena()) -lena = lena[200:300, 100:200] +astro = img_as_float(data.astronaut()) +astro = astro[30:180, 150:300] -noisy = lena + 0.6 * lena.std() * np.random.random(lena.shape) +noisy = astro + 0.3 * np.random.random(astro.shape) noisy = np.clip(noisy, 0, 1) -denoise = nl_means_denoising(noisy, 7, 9, 0.06) +denoise = nl_means_denoising(noisy, 7, 9, 0.08) fig, ax = plt.subplots(ncols=2, figsize=(8, 4)) diff --git a/skimage/restoration/__init__.py b/skimage/restoration/__init__.py index 9f708852..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, fast_nl_means_denoising +from .non_local_means import nl_means_denoising __all__ = ['wiener', 'unsupervised_wiener', @@ -31,5 +31,4 @@ __all__ = ['wiener', 'denoise_tv_bregman', 'denoise_tv_chambolle', 'denoise_bilateral', - 'nl_means_denoising', - 'fast_nl_means_denoising'] + 'nl_means_denoising'] diff --git a/skimage/restoration/_nl_means_denoising.pyx b/skimage/restoration/_nl_means_denoising.pyx index 424c02a5..1b8afc64 100644 --- a/skimage/restoration/_nl_means_denoising.pyx +++ b/skimage/restoration/_nl_means_denoising.pyx @@ -23,7 +23,7 @@ cdef inline float patch_distance_2d(DTYPE_t [:, :] p1, cdef float distance = 0 for i in range(s): # exp of large negative numbers will be 0, so we'd better stop - if distance > 4: + if distance > 5: return eps for j in range(s): tmp_diff = p1[i, j] - p2[i, j] @@ -43,7 +43,7 @@ cdef inline float patch_distance_2drgb(DTYPE_t [:, :, :] p1, cdef float distance = 0 for i in range(s): # exp of large negative numbers will be 0, so we'd better stop - if distance > 4: + if distance > 5: return eps for j in range(s): for color in range(3): @@ -62,7 +62,7 @@ cdef inline float patch_distance_3d(DTYPE_t [:, :, :] p1, cdef float tmp_diff for i in range(s): # exp of large negative numbers will be 0, so we'd better stop - if distance > 4: + if distance > 5: return eps for j in range(s): for k in range(s): @@ -186,7 +186,7 @@ def _nl_means_denoising_2drgb(image, int s=7, int d=13, float h=0.1): - (xg ** 2 + yg ** 2) / (2 * A ** 2)). astype(np.float32)) cdef float distance - w = 1. / (np.sum(w) * h ** 2) * w + w = 1. / (3 * np.sum(w) * h ** 2) * w # Coordinates of central pixel and patch bounds for x in range(offset, n_x + offset): x_start = x - offset @@ -371,7 +371,8 @@ def _fast_nl_means_denoising_2d(image, int s=7, int d=13, float h=0.1): integral[x - offset, y + offset] - \ integral[x + offset, y - offset] distance /= (s2 * h2) - if distance > 4: + # exp of large negative numbers is close to zero + if distance > 5: continue weight = alpha * exp(- distance) weights[x, y] += weight @@ -384,3 +385,212 @@ def _fast_nl_means_denoising_2d(image, int s=7, int d=13, float h=0.1): # except in padded zone result[x, y] /= weights[x, y] return result[pad_size: - pad_size, pad_size: - pad_size] + + +@cython.cdivision(True) +@cython.boundscheck(False) +def _fast_nl_means_denoising_2drgb(image, int s=7, int d=13, float h=0.1): + """ + Perform fast non-local means denoising on 2-D RGB array, with the outer + loop on patch shifts in order to reduce the number of operations. + + Parameters + ---------- + image: ndarray + 2-D RGB 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 + # 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 DTYPE_t [:, :, ::1] padded = np.ascontiguousarray(util.pad(image, + ((pad_size, pad_size), (pad_size, pad_size), (0, 0)), + mode='reflect').astype(np.float32)) + cdef DTYPE_t [:, :, ::1] result = np.zeros_like(padded) + cdef DTYPE_t [:, ::1] weights = np.zeros_like(padded[..., 0], order='C') + cdef DTYPE_t [:, ::1] integral = np.zeros_like(padded[..., 0], order='C') + cdef int n_x, n_y, t1, t2, x, y + cdef float weight, distance + cdef float alpha + cdef float h2 = h ** 2. + cdef float s2 = s ** 2. + cdef float h2s2 = 3 * h2 * s2 + n_x, n_y, _ = image.shape + n_x += 2 * pad_size + n_y += 2 * pad_size + # Outer loops on patch shifts + # With t2 >= 0, reference patch is always on the left of test patch + for t1 in range(-d, d + 1): + for t2 in range(0, d + 1): + # alpha is to account for patches on the same column + # distance is computed twice in this case + if t2 == 0 and t1 is not 0: + alpha = 0.5 + else: + alpha = 1. + integral = np.zeros_like(padded[..., 0], order='C') + 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)): + distance = ((padded[x, y, 0] - + padded[x + t1, y + t2, 0])**2 + +(padded[x, y, 1] - + padded[x + t1, y + t2, 1])**2 + +(padded[x, y, 2] - + padded[x + t1, y + t2, 2])**2) + integral[x, y] = distance + \ + 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 /= h2s2 + # exp of large negative numbers is close to zero + if distance > 5: + continue + weight = alpha * exp(- distance) + weights[x, y] += weight + weights[x + t1, y + t2] += weight + for ch in range(3): + result[x, y, ch] += weight * padded[x + t1, y + t2, ch] + result[x + t1, y + t2, ch] += weight * padded[x, y, ch] + for x in range(offset, n_x - offset): + for y in range(offset, n_y - offset): + for channel in range(3): + # no risk of division by zero + # except in padded zone + result[x, y, channel] /= weights[x, y] + 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. + """ + 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 DTYPE_t [:, :, ::1] padded = np.ascontiguousarray(util.pad(image, + pad_size, mode='reflect').astype(np.float32)) + cdef DTYPE_t [:, :, ::1] result = np.zeros_like(padded) + cdef DTYPE_t [:, :, ::1] weights = np.zeros_like(padded) + cdef DTYPE_t [:, :, ::1] integral = np.zeros_like(padded) + cdef int n_x, n_y, n_z, t1, t2, t3, x, y, z + cdef int x_integral_min, x_integral_max, y__integral_min, y_integral_max, \ + z_integral_min, z_integral_max + cdef int x_dist_min, x_dist_max, y_dist_min, y_dist_max, \ + z_dist_min, z_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_x, n_y, n_z = image.shape + n_x += 2 * pad_size + n_y += 2 * pad_size + n_z += 2 * pad_size + # Outer loops on patch shifts + # With t2 >= 0, reference patch is always on the left of test patch + for t1 in range(-d, d + 1): + x_integral_min = max(1, - t1) + x_integral_max = min(n_x, n_x - t1) + x_dist_min = max(offset, offset - t1) + x_dist_max = min(n_x - offset, n_x - offset - t1) + for t2 in range(-d, d + 1): + y_integral_min = max(1, - t2) + y_integral_max = min(n_y, n_y - t2) + y_dist_min = max(offset, offset - t2) + y_dist_max = min(n_y - offset, n_y - offset - t2) + for t3 in range(0, d + 1): + z_integral_min = max(1, - t3) + z_integral_max = min(n_z, n_z - t3) + z_dist_min = max(offset, offset - t3) + z_dist_max = min(n_z - offset, n_z - offset - t3) + # alpha is to account for patches on the same column + # distance is computed twice in this case + if t3 == 0 and (t1 is not 0 or t2 is not 0): + alpha = 0.5 + else: + alpha = 1. + integral = np.zeros_like(padded) + for x in range(x_integral_min, x_integral_max): + for y in range(y_integral_min, y_integral_max): + for z in range(z_integral_min, z_integral_max): + integral[x, y, z] = ((padded[x, y, z] - + padded[x + t1, y + t2, z + t3])**2 + + integral[x - 1, y, z] + + integral[x, y - 1, z] + + integral[x, y, z - 1] + + integral[x - 1, y - 1, z - 1] + - integral[x - 1, y - 1, z] + - integral[x, y - 1, z - 1] + - integral[x - 1, y, z - 1]) + for x in range(x_dist_min, x_dist_max): + for y in range(y_dist_min, y_dist_max): + for z in range(z_dist_min, z_dist_max): + distance = (integral[x + offset, + y + offset, + z + offset] + - integral[x - offset, y - offset, z - offset] + + integral[x - offset, y - offset, z + offset] + + integral[x - offset, y + offset, z - offset] + + integral[x + offset, y - offset, z - offset] + - integral[x - offset, y + offset, z + offset] + - integral[x + offset, y - offset, z + offset] + - integral[x + offset, y + offset, z - offset]) + distance /= s_cube_h_square + # exp of large negative numbers is close to zero + if distance > 5.: + continue + weight = alpha * exp(- distance) + weights[x, y, z] += weight + weights[x + t1, y + t2, z + t3] += weight + result[x, y, z] += weight * padded[x + t1, y + t2, + z + t3] + result[x + t1, y + t2, z + t3] += weight * \ + padded[x, y, z] + for x in range(offset, n_x - offset): + for y in range(offset, n_y - offset): + for z in range(offset, n_z - offset): + # I think there is no risk of division by zero + # except in padded zone + result[x, y, z] /= weights[x, y, z] + 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 index 5abbec5e..f727a12d 100644 --- a/skimage/restoration/non_local_means.py +++ b/skimage/restoration/non_local_means.py @@ -1,33 +1,42 @@ import numpy as np from skimage.restoration._nl_means_denoising import _nl_means_denoising_2d, \ _nl_means_denoising_2drgb, _nl_means_denoising_3d, \ - _fast_nl_means_denoising_2d + _fast_nl_means_denoising_2d, _fast_nl_means_denoising_3d, \ + _fast_nl_means_denoising_2drgb -def nl_means_denoising(image, patch_size=7, patch_distance=11, h=0.1): +def nl_means_denoising(image, patch_size=7, patch_distance=11, h=0.1, + fast_mode=True): """ Perform non-local means denoising on 2-D or 3-D grayscale images, and 2-D RGB images. Parameters ---------- - image: ndarray + image : ndarray input data to be denoised - patch_size: int, optional + patch_size : int, optional size of patches used for denoising - patch_distance: int, optional + patch_distance : int, optional maximal distance in pixels where to search patches used for denoising - h: float, optional + 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. + 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. + + 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 + result : ndarray denoised image, of same shape as `image`. See Also @@ -43,82 +52,17 @@ def nl_means_denoising(image, patch_size=7, patch_distance=11, h=0.1): 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 + 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. - 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. - - 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: - return np.array(_nl_means_denoising_2d(image, s=patch_size, - d=patch_distance, h=h)) - if image.ndim == 3 and image.shape[-1] > 4: # only grayscale - return np.array(_nl_means_denoising_3d(image, patch_size, - patch_distance, h)) - if image.ndim == 3 and image.shape[-1] == 3: # 2-D color (RGB) images - return np.array(_nl_means_denoising_2drgb(image, patch_size, - patch_distance, h)) - 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 + However, the default behavior corresponds to ``fast=True``, for which + another version of non-local means [2]_ is used, corresponding to a + complexity of image.size * patch_distance ** image.ndim @@ -132,14 +76,19 @@ def fast_nl_means_denoising(image, patch_size=7, patch_distance=11, h=0.1): `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. + 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] Jacques Froment. Parameter-Free Fast Pixelwise Non-Local Means + .. [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 @@ -147,11 +96,30 @@ def fast_nl_means_denoising(image, patch_size=7, patch_distance=11, h=0.1): >>> 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) + >>> denoised_a = nl_means_denoising(a, 7, 5, 0.1) """ if image.ndim == 2: - return np.array(_fast_nl_means_denoising_2d(image, s=patch_size, + if fast_mode: + return np.array(_fast_nl_means_denoising_2d(image, s=patch_size, + d=patch_distance, h=h)) + else: + return np.array(_nl_means_denoising_2d(image, s=patch_size, d=patch_distance, h=h)) + if image.ndim == 3 and image.shape[-1] > 4: # only 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)) + if image.ndim == 3 and image.shape[-1] == 3: # 2-D color (RGB) images + if fast_mode: + return np.array(_fast_nl_means_denoising_2drgb(image, patch_size, + patch_distance, h)) + else: + return np.array(_nl_means_denoising_2drgb(image, patch_size, + patch_distance, h)) else: - raise ValueError("Fast non local means denoising is only possible for \ - 2D grayscale images.") + raise ValueError("Non local means denoising is only possible for \ + 2D grayscale and RGB images or 3-D grayscale images.") + diff --git a/skimage/restoration/tests/test_denoise.py b/skimage/restoration/tests/test_denoise.py index 409e2d3d..35e945ec 100644 --- a/skimage/restoration/tests/test_denoise.py +++ b/skimage/restoration/tests/test_denoise.py @@ -148,16 +148,10 @@ 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) + denoised = restoration.nl_means_denoising(img, 7, 5, 0.1, fast_mode=True) # make sure noise is reduced assert img.std() > denoised.std() - - -def test_fast_nl_means_denoising_2d(): - img = np.zeros((40, 50)) - 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) + denoised = restoration.nl_means_denoising(img, 7, 5, 0.1, fast_mode=False) # make sure noise is reduced assert img.std() > denoised.std() @@ -168,7 +162,10 @@ def test_nl_means_denoising_2drgb(): # 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) + denoised = restoration.nl_means_denoising(img, 7, 9, 0.08, fast_mode=True) + # make sure noise is reduced + assert img.std() > denoised.std() + denoised = restoration.nl_means_denoising(img, 7, 9, 0.08, fast_mode=False) # make sure noise is reduced assert img.std() > denoised.std() @@ -177,10 +174,18 @@ 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) + denoised = restoration.nl_means_denoising(img, 5, 4, 0.1, fast_mode=True) + # make sure noise is reduced + assert img.std() > denoised.std() + denoised = restoration.nl_means_denoising(img, 5, 4, 0.1, fast_mode=False) # make sure noise is reduced assert img.std() > denoised.std() +def test_nl_means_denoising_wrong_dimension(): + img = np.zeros((5, 5, 5, 5)) + assert_raises(ValueError, restoration.nl_means_denoising, img) + + if __name__ == "__main__": run_module_suite()