From a5ed4acf86047a492f7e7df0c6dc6034cd8735f0 Mon Sep 17 00:00:00 2001 From: Emmanuelle Gouillart Date: Sun, 23 Feb 2014 13:28:27 +0100 Subject: [PATCH] Some improvements of non-local means denoising: - denoising RGB is now possible, and "colored patches" are then compared - the main function is now in a pure Python file so that default values of kw arguments are visible in the help - reduced the number of computations of patches bound (but this doesn't change much the total speed). - added an example for the gallery I also played with functions that could replace the exponential by a faster and less precise function, but it turns out that most of the time is spent in additions and multiplications when computing the distance between two patches. --- doc/examples/plot_nonlocal_means.py | 41 ++++ skimage/filter/_nl_means_denoising.pyx | 259 ++++++++++++++++--------- skimage/filter/nl_means_denoising.py | 74 +++++++ skimage/filter/tests/test_denoise.py | 169 ++++++++++++++++ 4 files changed, 451 insertions(+), 92 deletions(-) create mode 100644 doc/examples/plot_nonlocal_means.py create mode 100644 skimage/filter/nl_means_denoising.py create mode 100644 skimage/filter/tests/test_denoise.py diff --git a/doc/examples/plot_nonlocal_means.py b/doc/examples/plot_nonlocal_means.py new file mode 100644 index 00000000..790c9bdf --- /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 Lena 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.filter import nl_means_denoising + + +lena = img_as_float(data.lena()) +lena = lena[200:300, 100:200] + +noisy = lena + 0.6 * lena.std() * np.random.random(lena.shape) +noisy = np.clip(noisy, 0, 1) + +denoise = nl_means_denoising(noisy, 7, 9, 0.06) + +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/filter/_nl_means_denoising.pyx b/skimage/filter/_nl_means_denoising.pyx index 72a45880..cadbfd09 100644 --- a/skimage/filter/_nl_means_denoising.pyx +++ b/skimage/filter/_nl_means_denoising.pyx @@ -4,18 +4,27 @@ cimport numpy as np cimport cython from libc.math cimport exp -DTYPE = np.float ctypedef np.float32_t DTYPE_t +cdef eps = 1.e-8 + @cython.boundscheck(False) -cdef inline float patch_distance2d(DTYPE_t [:, :] p1, +cdef inline float patch_distance_2d(DTYPE_t [:, :] p1, DTYPE_t [:, :] p2, DTYPE_t [:, ::] w, int s): 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 eps 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 > 4: + return eps for j in range(s): tmp_diff = p1[i, j] - p2[i, j] distance += w[i, j] * tmp_diff * tmp_diff @@ -24,13 +33,37 @@ cdef inline float patch_distance2d(DTYPE_t [:, :] p1, @cython.boundscheck(False) -cdef inline float patch_distance(DTYPE_t [:, :, :] p1, +cdef inline float patch_distance_2drgb(DTYPE_t [:, :, :] p1, + DTYPE_t [:, :, :] p2, + DTYPE_t [:, ::] w, int s): + 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 > 4: + return eps + 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(DTYPE_t [:, :, :] p1, DTYPE_t [:, :, :] p2, DTYPE_t [:, :, ::] w, int s): 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 > 4: + return eps for j in range(s): for k in range(s): tmp_diff = p1[i, j, k] - p2[i, j, k] @@ -65,37 +98,129 @@ def _nl_means_denoising_2d(image, int s=7, int d=13, float h=0.1): cdef int n_x, n_y n_x, n_y = 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] result = padded.copy() + # We normalize by the image contrast, and divide by 3 because of 3 channels + h *= (np.max(padded) - np.min(padded)) / 3. + cdef float A = ((s - 1.) / 4.) + cdef float new_value + cdef float weight_sum, weight + xg, yg = np.mgrid[-offset:offset + 1, -offset:offset + 1] + cdef DTYPE_t [:, ::1] w = np.ascontiguousarray(np.exp( + - (xg ** 2 + yg ** 2) / (2 * A ** 2)). + astype(np.float32)) + cdef float distance + cdef int x, y, i, j + cdef int x_start, x_end, y_start, y_end + cdef int x_start_i, x_end_i, y_start_j, y_end_j + w = 1. / (np.sum(w) * 2 * h ** 2) * w + # Coordinates of central pixel and patch bounds + for x in range(offset, n_x + offset): + x_start = x - offset + x_end = x + offset + 1 + for y in range(offset, n_y + offset): + new_value = 0 + weight_sum = 0 + y_start = y - offset + y_end = y + offset + 1 + # Coordinates of test pixel and patch bounds + for i in range(max(- d, offset - x), + min(d + 1, n_x - x - 1)): + x_start_i = x_start + i + x_end_i = x_end + i + for j in range(max(- d, offset - y), + min(d + 1, n_y - y - 1)): + y_start_j = y_start + j + y_end_j = y_end + j + weight = patch_distance_2d( + padded[x_start: x_end, + y_start: y_end], + padded[x_start_i: x_end_i, + y_start_j: y_end_j], + w, s) + weight_sum += weight + new_value += weight * padded[x + i, y + j] + result[x, y] = new_value / weight_sum + return result[offset:-offset, offset:-offset] + + +@cython.cdivision(True) +@cython.boundscheck(False) +def _nl_means_denoising_2drgb(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. + """ + if s % 2 == 0: + s += 1 # odd value for symmetric patch + cdef int n_x, n_y + n_x, n_y, _ = image.shape + cdef int offset = s / 2 + cdef int x, y, i, j, color + cdef int x_start, x_end, y_start, y_end + cdef int x_start_i, x_end_i, y_start_j, y_end_j + cdef DTYPE_t [::1] new_values = np.zeros(3).astype(np.float32) + cdef DTYPE_t [:, :, ::1] padded = np.ascontiguousarray(util.pad(image, + ((offset, offset), (offset, offset), (0, 0)), + mode='reflect').astype(np.float32)) + cdef DTYPE_t [:, :, ::1] result = padded.copy() h *= (np.max(padded) - np.min(padded)) cdef float A = ((s - 1.) / 4.) cdef float new_value cdef float weight_sum, weight xg, yg = np.mgrid[-offset:offset + 1, -offset:offset + 1] cdef DTYPE_t [:, ::1] w = np.ascontiguousarray(np.exp( - - (xg ** 2 + yg ** 2)/(2 * A ** 2)). + - (xg ** 2 + yg ** 2) / (2 * A ** 2)). astype(np.float32)) cdef float distance - cdef int x, y, i, j - w = 1./ (np.sum(w) * 2 * h ** 2) * w + w = 1. / (np.sum(w) * 2 * h ** 2) * w + # Coordinates of central pixel and patch bounds for x in range(offset, n_x + offset): + x_start = x - offset + x_end = x + offset + 1 for y in range(offset, n_y + offset): - new_value = 0 + for color in range(3): + new_values[color] = 0 weight_sum = 0 + y_start = y - offset + y_end = y + offset + 1 + # Coordinates of test pixel and patch bounds for i in range(max(- d, offset - x), min(d + 1, n_x - x - 1)): + x_start_i = x_start + i + x_end_i = x_end + i for j in range(max(- d, offset - y), min(d + 1, n_y - y - 1)): - weight = patch_distance2d( - padded[x - offset: x + offset + 1, - y - offset: y + offset + 1], - padded[x + i - offset: x + i + offset + 1, - y + j - offset: y + j + offset + 1], + y_start_j = y_start + j + y_end_j = y_end + j + weight = patch_distance_2drgb( + padded[x_start: x_end, + y_start: y_end, :], + padded[x_start_i: x_end_i, + y_start_j: y_end_j, :], w, s) weight_sum += weight - new_value += weight * padded[x + i, y + j] - result[x, y] = new_value / weight_sum + for color in range(3): + new_values[color] += weight * padded[x + i, y + j, + color] + for color in range(3): + result[x, y, color] = new_values[color] / weight_sum return result[offset:-offset, offset:-offset] @@ -125,6 +250,7 @@ def _nl_means_denoising_3d(image, int s=7, cdef int n_x, n_y, n_z 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] result = padded.copy() @@ -132,101 +258,50 @@ def _nl_means_denoising_3d(image, int s=7, cdef float A = ((s - 1.) / 4.) cdef float new_value cdef float weight_sum, weight - xg, yg, zg = np.mgrid[-offset: offset + 1, -offset: offset+1, + xg, yg, zg = np.mgrid[-offset: offset + 1, -offset: offset + 1, -offset: offset + 1] cdef DTYPE_t [:, :, ::1] w = np.ascontiguousarray(np.exp( - - (xg ** 2 + yg ** 2 + zg ** 2)/(2 * A ** 2)). + - (xg ** 2 + yg ** 2 + zg ** 2) / (2 * A ** 2)). astype(np.float32)) cdef float distance cdef int x, y, z, i, j, k - w = 1./ (np.sum(w) * 2 * h ** 2) * w + cdef int x_start, x_end, y_start, y_end, z_start, z_end + cdef int x_start_i, x_end_i, y_start_j, y_end_j, z_start_k, z_end_k + w = 1. / (np.sum(w) * 2 * h ** 2) * w + # Coordinates of central pixel and patch bounds for x in range(offset, n_x + offset): + x_start = x - offset + x_end = x + offset + 1 for y in range(offset, n_y + offset): + y_start = y - offset + y_end = y + offset + 1 for z in range(offset, n_z + offset): + z_start = z - offset + z_end = z + offset + 1 new_value = 0 weight_sum = 0 + # Coordinates of test pixel and patch bounds for i in range(max(- d, offset - x), min(d + 1, n_x - x - 1)): + x_start_i = x_start + i + x_end_i = x_end + i for j in range(max(- d, offset - y), min(d + 1, n_y - y - 1)): + y_start_j = y_start + j + y_end_j = y_end + j for k in range(max(- d, offset - z), min(d + 1, n_z - z - 1)): - weight = patch_distance( - padded[x - offset: x + offset +1, - y - offset: y + offset +1, - z - offset: z + offset +1], - padded[x + i - offset: x + i + offset +1, - y + j - offset: y + j + offset +1, - z + k - offset: z + k + offset +1], + z_start_k = z_start + k + z_end_k = z_end + k + weight = patch_distance_3d( + padded[x_start: x_end, + y_start: y_end, + z_start: z_end], + padded[x_start_i: x_end_i, + y_start_j: y_end_j, + z_start_k: z_end_k], w, s) weight_sum += weight 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] - - -def nl_means_denoising(image, patch_size=7, patch_distance=11, h=0.1): - """ - Perform non-local means denoising on 2-D or 3-D grayscale arrays - - 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. - - Returns - ------- - - result: ndarray - denoised image, of same shape as `image`. - - 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_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, patch_size, - patch_distance, 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)) - else: - raise ValueError("Non local means denoising is only possible for \ - 2D and 3-D grayscale images.") diff --git a/skimage/filter/nl_means_denoising.py b/skimage/filter/nl_means_denoising.py new file mode 100644 index 00000000..173d6cf7 --- /dev/null +++ b/skimage/filter/nl_means_denoising.py @@ -0,0 +1,74 @@ +import numpy as np +from _nl_means_denoising import _nl_means_denoising_2d, \ + _nl_means_denoising_2drgb, _nl_means_denoising_3d + +def nl_means_denoising(image, patch_size=7, patch_distance=11, h=0.1): + """ + Perform non-local means denoising on 2-D or 3-D grayscale images, and + 2-D RGB 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`. + + 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_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, patch_size, + patch_distance, 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.") diff --git a/skimage/filter/tests/test_denoise.py b/skimage/filter/tests/test_denoise.py new file mode 100644 index 00000000..206b4027 --- /dev/null +++ b/skimage/filter/tests/test_denoise.py @@ -0,0 +1,169 @@ +import numpy as np +from numpy.testing import run_module_suite, assert_raises, assert_equal + +from skimage import filter, data, color, img_as_float + + +np.random.seed(1234) + + +lena = img_as_float(data.lena()[:256, :256]) +lena_gray = color.rgb2gray(lena) + + +def test_denoise_tv_chambolle_2d(): + # lena image + img = lena_gray + # add noise to lena + img += 0.5 * img.std() * np.random.random(img.shape) + # clip noise so that it does not exceed allowed range for float images. + img = np.clip(img, 0, 1) + # denoise + denoised_lena = filter.denoise_tv_chambolle(img, weight=60.0) + # which dtype? + assert denoised_lena.dtype in [np.float, np.float32, np.float64] + from scipy import ndimage + grad = ndimage.morphological_gradient(img, size=((3, 3))) + grad_denoised = ndimage.morphological_gradient( + denoised_lena, size=((3, 3))) + # test if the total variation has decreased + assert grad_denoised.dtype == np.float + assert (np.sqrt((grad_denoised**2).sum()) + < np.sqrt((grad**2).sum()) / 2) + + +def test_denoise_tv_chambolle_multichannel(): + denoised0 = filter.denoise_tv_chambolle(lena[..., 0], weight=60.0) + denoised = filter.denoise_tv_chambolle(lena, weight=60.0, multichannel=True) + assert_equal(denoised[..., 0], denoised0) + + +def test_denoise_tv_chambolle_float_result_range(): + # lena image + img = lena_gray + int_lena = np.multiply(img, 255).astype(np.uint8) + assert np.max(int_lena) > 1 + denoised_int_lena = filter.denoise_tv_chambolle(int_lena, weight=60.0) + # test if the value range of output float data is within [0.0:1.0] + assert denoised_int_lena.dtype == np.float + assert np.max(denoised_int_lena) <= 1.0 + assert np.min(denoised_int_lena) >= 0.0 + + +def test_denoise_tv_chambolle_3d(): + """Apply the TV denoising algorithm on a 3D image representing a sphere.""" + x, y, z = np.ogrid[0:40, 0:40, 0:40] + mask = (x - 22)**2 + (y - 20)**2 + (z - 17)**2 < 8**2 + mask = 100 * mask.astype(np.float) + mask += 60 + mask += 20 * np.random.random(mask.shape) + mask[mask < 0] = 0 + mask[mask > 255] = 255 + res = filter.denoise_tv_chambolle(mask.astype(np.uint8), weight=100) + assert res.dtype == np.float + assert res.std() * 255 < mask.std() + + # test wrong number of dimensions + assert_raises(ValueError, filter.denoise_tv_chambolle, + np.random.random((8, 8, 8, 8))) + + +def test_denoise_tv_bregman_2d(): + img = lena_gray + # add some random noise + img += 0.5 * img.std() * np.random.random(img.shape) + img = np.clip(img, 0, 1) + + out1 = filter.denoise_tv_bregman(img, weight=10) + out2 = filter.denoise_tv_bregman(img, weight=5) + + # make sure noise is reduced + assert img.std() > out1.std() + assert out1.std() > out2.std() + + +def test_denoise_tv_bregman_float_result_range(): + # lena image + img = lena_gray + int_lena = np.multiply(img, 255).astype(np.uint8) + assert np.max(int_lena) > 1 + denoised_int_lena = filter.denoise_tv_bregman(int_lena, weight=60.0) + # test if the value range of output float data is within [0.0:1.0] + assert denoised_int_lena.dtype == np.float + assert np.max(denoised_int_lena) <= 1.0 + assert np.min(denoised_int_lena) >= 0.0 + + +def test_denoise_tv_bregman_3d(): + img = lena + # add some random noise + img += 0.5 * img.std() * np.random.random(img.shape) + img = np.clip(img, 0, 1) + + out1 = filter.denoise_tv_bregman(img, weight=10) + out2 = filter.denoise_tv_bregman(img, weight=5) + + # make sure noise is reduced + assert img.std() > out1.std() + assert out1.std() > out2.std() + + +def test_denoise_bilateral_2d(): + img = lena_gray + # add some random noise + img += 0.5 * img.std() * np.random.random(img.shape) + img = np.clip(img, 0, 1) + + out1 = filter.denoise_bilateral(img, sigma_range=0.1, sigma_spatial=20) + out2 = filter.denoise_bilateral(img, sigma_range=0.2, sigma_spatial=30) + + # make sure noise is reduced + assert img.std() > out1.std() + assert out1.std() > out2.std() + + +def test_denoise_bilateral_3d(): + img = lena + # add some random noise + img += 0.5 * img.std() * np.random.random(img.shape) + img = np.clip(img, 0, 1) + + out1 = filter.denoise_bilateral(img, sigma_range=0.1, sigma_spatial=20) + out2 = filter.denoise_bilateral(img, sigma_range=0.2, sigma_spatial=30) + + # make sure noise is reduced + assert img.std() > out1.std() + assert out1.std() > out2.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 = filter.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 = lena[-100:, -100:] + # add some random noise + img += 0.5 * img.std() * np.random.random(img.shape) + img = np.clip(img, 0, 1) + denoised = filter.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 = filter.nl_means_denoising(img, 5, 4, 0.1) + # make sure noise is reduced + assert img.std() > denoised.std() + + +if __name__ == "__main__": + run_module_suite()