From a508ec54ca8516b40a74aeb2a13ada9c1de39da0 Mon Sep 17 00:00:00 2001 From: Emmanuelle Gouillart Date: Sun, 2 Feb 2014 18:37:33 +0100 Subject: [PATCH 01/30] 2-D and 3-D implementation of NL-means denoising --- skimage/filter/_nl_means_denoising.pyx | 232 +++++++++++++++++++++++++ skimage/filters/setup.py | 4 + 2 files changed, 236 insertions(+) create mode 100644 skimage/filter/_nl_means_denoising.pyx diff --git a/skimage/filter/_nl_means_denoising.pyx b/skimage/filter/_nl_means_denoising.pyx new file mode 100644 index 00000000..72a45880 --- /dev/null +++ b/skimage/filter/_nl_means_denoising.pyx @@ -0,0 +1,232 @@ +import numpy as np +from skimage import util +cimport numpy as np +cimport cython +from libc.math cimport exp + +DTYPE = np.float +ctypedef np.float32_t DTYPE_t + + +@cython.boundscheck(False) +cdef inline float patch_distance2d(DTYPE_t [:, :] p1, + DTYPE_t [:, :] p2, + DTYPE_t [:, ::] w, int s): + cdef int i, j + cdef float distance = 0 + cdef float tmp_diff + for i in range(s): + 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(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): + 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 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). 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 DTYPE_t [:, ::1] padded = np.ascontiguousarray(util.pad(image, + offset, 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)). + astype(np.float32)) + cdef float distance + cdef int x, y, i, j + w = 1./ (np.sum(w) * 2 * h ** 2) * w + for x in range(offset, n_x + offset): + for y in range(offset, n_y + offset): + new_value = 0 + weight_sum = 0 + for i in range(max(- d, offset - x), + min(d + 1, n_x - x - 1)): + 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], + 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_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) + """ + if s % 2 == 0: + s += 1 # odd value for symmetric patch + cdef int n_x, n_y, n_z + n_x, n_y, n_z = image.shape + cdef int offset = s / 2 + cdef DTYPE_t [:, :, ::1] padded = np.ascontiguousarray(util.pad(image, + offset, 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, 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)). + astype(np.float32)) + cdef float distance + cdef int x, y, z, i, j, k + w = 1./ (np.sum(w) * 2 * h ** 2) * w + for x in range(offset, n_x + offset): + for y in range(offset, n_y + offset): + for z in range(offset, n_z + offset): + new_value = 0 + weight_sum = 0 + for i in range(max(- d, offset - x), + min(d + 1, n_x - x - 1)): + for j in range(max(- d, offset - y), + min(d + 1, n_y - y - 1)): + 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], + 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/filters/setup.py b/skimage/filters/setup.py index 8bfebcb7..21728902 100644 --- a/skimage/filters/setup.py +++ b/skimage/filters/setup.py @@ -14,6 +14,7 @@ def configuration(parent_package='', top_path=None): config.add_data_dir('rank/tests') cython(['_ctmf.pyx'], working_path=base_path) + cython(['_nl_means_denoising.pyx'], working_path=base_path) cython(['rank/core_cy.pyx'], working_path=base_path) cython(['rank/generic_cy.pyx'], working_path=base_path) cython(['rank/percentile_cy.pyx'], working_path=base_path) @@ -21,6 +22,9 @@ def configuration(parent_package='', top_path=None): config.add_extension('_ctmf', sources=['_ctmf.c'], include_dirs=[get_numpy_include_dirs()]) + config.add_extension('_nl_means_denoising', + sources=['_nl_means_denoising.c'], + include_dirs=[get_numpy_include_dirs()]) config.add_extension('rank.core_cy', sources=['rank/core_cy.c'], include_dirs=[get_numpy_include_dirs()]) config.add_extension('rank.generic_cy', sources=['rank/generic_cy.c'], From a5ed4acf86047a492f7e7df0c6dc6034cd8735f0 Mon Sep 17 00:00:00 2001 From: Emmanuelle Gouillart Date: Sun, 23 Feb 2014 13:28:27 +0100 Subject: [PATCH 02/30] 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() From 5f89395fcd3096fee43b3565eb9ba534b1cd1975 Mon Sep 17 00:00:00 2001 From: Emmanuelle Gouillart Date: Sun, 23 Feb 2014 18:30:58 +0100 Subject: [PATCH 03/30] Modified indentation in places (PEP8) --- skimage/filter/_nl_means_denoising.pyx | 20 ++++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/skimage/filter/_nl_means_denoising.pyx b/skimage/filter/_nl_means_denoising.pyx index cadbfd09..bdf4fd6e 100644 --- a/skimage/filter/_nl_means_denoising.pyx +++ b/skimage/filter/_nl_means_denoising.pyx @@ -11,8 +11,8 @@ cdef eps = 1.e-8 @cython.boundscheck(False) cdef inline float patch_distance_2d(DTYPE_t [:, :] p1, - DTYPE_t [:, :] p2, - DTYPE_t [:, ::] w, int s): + 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 @@ -34,8 +34,8 @@ cdef inline float patch_distance_2d(DTYPE_t [:, :] p1, @cython.boundscheck(False) cdef inline float patch_distance_2drgb(DTYPE_t [:, :, :] p1, - DTYPE_t [:, :, :] p2, - DTYPE_t [:, ::] w, int s): + DTYPE_t [:, :, :] p2, + DTYPE_t [:, ::] w, int s): cdef int i, j cdef int center = s / 2 cdef int color @@ -55,8 +55,8 @@ cdef inline float patch_distance_2drgb(DTYPE_t [:, :, :] p1, @cython.boundscheck(False) cdef inline float patch_distance_3d(DTYPE_t [:, :, :] p1, - DTYPE_t [:, :, :] p2, - DTYPE_t [:, :, ::] w, int s): + DTYPE_t [:, :, :] p2, + DTYPE_t [:, :, ::] w, int s): cdef int i, j, k cdef float distance = 0 cdef float tmp_diff @@ -139,7 +139,7 @@ def _nl_means_denoising_2d(image, int s=7, int d=13, float h=0.1): y_start: y_end], padded[x_start_i: x_end_i, y_start_j: y_end_j], - w, s) + w, s) weight_sum += weight new_value += weight * padded[x + i, y + j] result[x, y] = new_value / weight_sum @@ -214,7 +214,7 @@ def _nl_means_denoising_2drgb(image, int s=7, int d=13, float h=0.1): y_start: y_end, :], padded[x_start_i: x_end_i, y_start_j: y_end_j, :], - w, s) + w, s) weight_sum += weight for color in range(3): new_values[color] += weight * padded[x + i, y + j, @@ -295,8 +295,8 @@ def _nl_means_denoising_3d(image, int s=7, z_end_k = z_end + k weight = patch_distance_3d( padded[x_start: x_end, - y_start: y_end, - z_start: z_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], From 2e91b0db953ad7e72eb09ce6f421a66b5f2b30c1 Mon Sep 17 00:00:00 2001 From: emmanuelle Date: Sat, 6 Dec 2014 13:16:40 +0100 Subject: [PATCH 04/30] Moved non-local means denoising to restoration submodule --- doc/examples/plot_nonlocal_means.py | 2 +- skimage/filters/setup.py | 4 --- skimage/restoration/__init__.py | 5 +-- .../_nl_means_denoising.pyx | 0 .../nl_means_denoising.py | 0 skimage/restoration/setup.py | 4 +++ skimage/restoration/tests/test_denoise.py | 32 ++++++++++++++++++- 7 files changed, 39 insertions(+), 8 deletions(-) rename skimage/{filter => restoration}/_nl_means_denoising.pyx (100%) rename skimage/{filter => restoration}/nl_means_denoising.py (100%) diff --git a/doc/examples/plot_nonlocal_means.py b/doc/examples/plot_nonlocal_means.py index 790c9bdf..c6879be4 100644 --- a/doc/examples/plot_nonlocal_means.py +++ b/doc/examples/plot_nonlocal_means.py @@ -15,7 +15,7 @@ import numpy as np import matplotlib.pyplot as plt from skimage import data, img_as_float -from skimage.filter import nl_means_denoising +from skimage.restoration import nl_means_denoising lena = img_as_float(data.lena()) diff --git a/skimage/filters/setup.py b/skimage/filters/setup.py index 21728902..8bfebcb7 100644 --- a/skimage/filters/setup.py +++ b/skimage/filters/setup.py @@ -14,7 +14,6 @@ def configuration(parent_package='', top_path=None): config.add_data_dir('rank/tests') cython(['_ctmf.pyx'], working_path=base_path) - cython(['_nl_means_denoising.pyx'], working_path=base_path) cython(['rank/core_cy.pyx'], working_path=base_path) cython(['rank/generic_cy.pyx'], working_path=base_path) cython(['rank/percentile_cy.pyx'], working_path=base_path) @@ -22,9 +21,6 @@ def configuration(parent_package='', top_path=None): config.add_extension('_ctmf', sources=['_ctmf.c'], include_dirs=[get_numpy_include_dirs()]) - config.add_extension('_nl_means_denoising', - sources=['_nl_means_denoising.c'], - include_dirs=[get_numpy_include_dirs()]) config.add_extension('rank.core_cy', sources=['rank/core_cy.c'], include_dirs=[get_numpy_include_dirs()]) config.add_extension('rank.generic_cy', sources=['rank/generic_cy.c'], diff --git a/skimage/restoration/__init__.py b/skimage/restoration/__init__.py index 2b593ccf..b008c638 100644 --- a/skimage/restoration/__init__.py +++ b/skimage/restoration/__init__.py @@ -22,7 +22,7 @@ from .deconvolution import wiener, unsupervised_wiener, richardson_lucy from .unwrap import unwrap_phase from ._denoise import denoise_tv_chambolle, denoise_tv_bregman, \ denoise_bilateral - +from .nl_means_denoising import nl_means_denoising __all__ = ['wiener', 'unsupervised_wiener', @@ -30,4 +30,5 @@ __all__ = ['wiener', 'unwrap_phase', 'denoise_tv_bregman', 'denoise_tv_chambolle', - 'denoise_bilateral'] + 'denoise_bilateral', + 'nl_means_denoising'] diff --git a/skimage/filter/_nl_means_denoising.pyx b/skimage/restoration/_nl_means_denoising.pyx similarity index 100% rename from skimage/filter/_nl_means_denoising.pyx rename to skimage/restoration/_nl_means_denoising.pyx diff --git a/skimage/filter/nl_means_denoising.py b/skimage/restoration/nl_means_denoising.py similarity index 100% rename from skimage/filter/nl_means_denoising.py rename to skimage/restoration/nl_means_denoising.py diff --git a/skimage/restoration/setup.py b/skimage/restoration/setup.py index e20073e0..12638107 100644 --- a/skimage/restoration/setup.py +++ b/skimage/restoration/setup.py @@ -17,6 +17,7 @@ def configuration(parent_package='', top_path=None): cython(['_unwrap_2d.pyx'], working_path=base_path) cython(['_unwrap_3d.pyx'], working_path=base_path) cython(['_denoise_cy.pyx'], working_path=base_path) + cython(['_nl_means_denoising.pyx'], working_path=base_path) config.add_extension('_unwrap_1d', sources=['_unwrap_1d.c'], include_dirs=[get_numpy_include_dirs()]) @@ -28,6 +29,9 @@ def configuration(parent_package='', top_path=None): include_dirs=[get_numpy_include_dirs()]) config.add_extension('_denoise_cy', sources=['_denoise_cy.c'], include_dirs=[get_numpy_include_dirs(), '../_shared']) + config.add_extension('_nl_means_denoising', + sources=['_nl_means_denoising.c'], + include_dirs=[get_numpy_include_dirs()]) return config diff --git a/skimage/restoration/tests/test_denoise.py b/skimage/restoration/tests/test_denoise.py index 46ed4b35..77d38557 100644 --- a/skimage/restoration/tests/test_denoise.py +++ b/skimage/restoration/tests/test_denoise.py @@ -46,7 +46,8 @@ def test_denoise_tv_chambolle_float_result_range(): img = astro_gray int_astro = np.multiply(img, 255).astype(np.uint8) assert np.max(int_astro) > 1 - denoised_int_astro = restoration.denoise_tv_chambolle(int_astro, weight=60.0) + denoised_int_astro = restoration.denoise_tv_chambolle(int_astro, + weight=60.0) # test if the value range of output float data is within [0.0:1.0] assert denoised_int_astro.dtype == np.float assert np.max(denoised_int_astro) <= 1.0 @@ -143,5 +144,34 @@ def test_denoise_bilateral_3d(): assert out1[30:45, 5:15].std() > out2[30:45, 5:15].std() +def test_nl_means_denoising_2d(): + img = np.zeros((40, 40)) + img[10:-10, 10:-10] = 1. + img += 0.3*np.random.randn(*img.shape) + denoised = restoration.nl_means_denoising(img, 7, 5, 0.1) + # make sure noise is reduced + assert img.std() > denoised.std() + + +def test_nl_means_denoising_2drgb(): + # reduce image size because nl means is very slow + img = np.copy(astro[:50, :50]) + # add some random noise + img += 0.5 * img.std() * np.random.random(img.shape) + img = np.clip(img, 0, 1) + denoised = restoration.nl_means_denoising(img, 7, 9, 0.08) + # make sure noise is reduced + assert img.std() > denoised.std() + + +def test_nl_means_denoising_3d(): + img = np.zeros((20, 20, 10)) + img[5:-5, 5:-5, 3:-3] = 1. + img += 0.3*np.random.randn(*img.shape) + denoised = restoration.nl_means_denoising(img, 5, 4, 0.1) + # make sure noise is reduced + assert img.std() > denoised.std() + + if __name__ == "__main__": run_module_suite() From cc6882b5845ff15edff6527b12010e0c67185f2b Mon Sep 17 00:00:00 2001 From: emmanuelle Date: Sat, 6 Dec 2014 13:19:27 +0100 Subject: [PATCH 05/30] Modified bento file with non-local means extension --- bento.info | 3 +++ 1 file changed, 3 insertions(+) diff --git a/bento.info b/bento.info index 6a8a0000..5b666d9a 100644 --- a/bento.info +++ b/bento.info @@ -151,6 +151,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 From ecae7843e9757217921d905c9ff1215afc4d24bd Mon Sep 17 00:00:00 2001 From: emmanuelle Date: Sat, 6 Dec 2014 13:46:13 +0100 Subject: [PATCH 06/30] Removed test file that has been moved to restoration submodule --- skimage/filter/tests/test_denoise.py | 169 --------------------------- 1 file changed, 169 deletions(-) delete mode 100644 skimage/filter/tests/test_denoise.py diff --git a/skimage/filter/tests/test_denoise.py b/skimage/filter/tests/test_denoise.py deleted file mode 100644 index 206b4027..00000000 --- a/skimage/filter/tests/test_denoise.py +++ /dev/null @@ -1,169 +0,0 @@ -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() From e26f9028eb80fe3b9816f56fa0404bcc6184972a Mon Sep 17 00:00:00 2001 From: emmanuelle Date: Sat, 6 Dec 2014 13:58:09 +0100 Subject: [PATCH 07/30] Modified import to be compatible with Python 3 --- skimage/restoration/nl_means_denoising.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/skimage/restoration/nl_means_denoising.py b/skimage/restoration/nl_means_denoising.py index 173d6cf7..c106795e 100644 --- a/skimage/restoration/nl_means_denoising.py +++ b/skimage/restoration/nl_means_denoising.py @@ -1,5 +1,5 @@ import numpy as np -from _nl_means_denoising import _nl_means_denoising_2d, \ +from skimage.restoration._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): From 4004f048ef37e24ef6b7125a7b1ca65af54abc0b Mon Sep 17 00:00:00 2001 From: emmanuelle Date: Sun, 7 Dec 2014 12:02:05 +0100 Subject: [PATCH 08/30] Faster version of non-local means denoising for 2D greyscale images --- skimage/restoration/__init__.py | 5 +- skimage/restoration/_nl_means_denoising.pyx | 82 +++++++++++++++++++- skimage/restoration/nl_means_denoising.py | 85 ++++++++++++++++++++- skimage/restoration/tests/test_denoise.py | 9 +++ 4 files changed, 176 insertions(+), 5 deletions(-) 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]) From 1d8e02bf387fe667ffba378960a8e2bc3b44556f Mon Sep 17 00:00:00 2001 From: emmanuelle Date: Sat, 10 Jan 2015 12:15:26 +0100 Subject: [PATCH 09/30] Removed normalization of h parameters because float images will be normalized anyway. --- skimage/restoration/_nl_means_denoising.pyx | 24 ++++++++++----------- skimage/restoration/nl_means_denoising.py | 8 +++---- 2 files changed, 15 insertions(+), 17 deletions(-) diff --git a/skimage/restoration/_nl_means_denoising.pyx b/skimage/restoration/_nl_means_denoising.pyx index 7a1b866f..f91e0170 100644 --- a/skimage/restoration/_nl_means_denoising.pyx +++ b/skimage/restoration/_nl_means_denoising.pyx @@ -102,8 +102,6 @@ def _nl_means_denoising_2d(image, int s=7, int d=13, float h=0.1): 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 @@ -115,7 +113,7 @@ def _nl_means_denoising_2d(image, int s=7, int d=13, float h=0.1): 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 + w = 1. / (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 @@ -180,7 +178,6 @@ def _nl_means_denoising_2drgb(image, int s=7, int d=13, float h=0.1): ((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 @@ -189,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) * 2 * h ** 2) * w + w = 1. / (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 @@ -255,7 +252,6 @@ def _nl_means_denoising_3d(image, int s=7, 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.) cdef float new_value cdef float weight_sum, weight @@ -268,7 +264,7 @@ def _nl_means_denoising_3d(image, int s=7, cdef int x, y, z, i, j, k 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 + w = 1. / (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 @@ -337,17 +333,17 @@ def _fast_nl_means_denoising_2d(image, int s=7, int d=13, float h=0.1): 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] result = np.zeros_like(padded) 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. + 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 + # Outer loops on patch shifts for t1 in range(-d, d + 1): for t2 in range(-d, d + 1): integral = np.zeros_like(padded) @@ -371,8 +367,10 @@ 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] - \ integral[x + offset, y - offset] - distance /= s2 - weight = exp(- distance / h2) + distance /= (s2 * h2) + if distance > 4: + continue + weight = exp(- distance) weights[x, y] += weight result[x, y] += weight * padded[x + t1, y + t2] for x in range(offset, n_x - offset): diff --git a/skimage/restoration/nl_means_denoising.py b/skimage/restoration/nl_means_denoising.py index 79c67e04..5abbec5e 100644 --- a/skimage/restoration/nl_means_denoising.py +++ b/skimage/restoration/nl_means_denoising.py @@ -66,8 +66,8 @@ def nl_means_denoising(image, patch_size=7, patch_distance=11, h=0.1): >>> 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)) + 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)) @@ -150,8 +150,8 @@ def fast_nl_means_denoising(image, patch_size=7, patch_distance=11, h=0.1): >>> 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)) + return np.array(_fast_nl_means_denoising_2d(image, s=patch_size, + d=patch_distance, h=h)) else: raise ValueError("Fast non local means denoising is only possible for \ 2D grayscale images.") From ecaa9591954f70c49f005f0ea7ad5505a8784ff9 Mon Sep 17 00:00:00 2001 From: emmanuelle Date: Sat, 10 Jan 2015 16:40:43 +0100 Subject: [PATCH 10/30] Implemented asymmetric distance computation to save speed factor of 2 --- skimage/restoration/_nl_means_denoising.pyx | 31 +++++++++++---------- skimage/restoration/tests/test_denoise.py | 2 +- 2 files changed, 18 insertions(+), 15 deletions(-) diff --git a/skimage/restoration/_nl_means_denoising.pyx b/skimage/restoration/_nl_means_denoising.pyx index f91e0170..424c02a5 100644 --- a/skimage/restoration/_nl_means_denoising.pyx +++ b/skimage/restoration/_nl_means_denoising.pyx @@ -314,7 +314,7 @@ def _fast_nl_means_denoising_2d(image, int s=7, int d=13, float h=0.1): Parameters ---------- image: ndarray - input data to be denoised + 2-D input data to be denoised s: int, optional size of patches used for denoising @@ -326,11 +326,12 @@ def _fast_nl_means_denoising_2d(image, int s=7, int d=13, float h=0.1): 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 + # 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) @@ -338,21 +339,23 @@ def _fast_nl_means_denoising_2d(image, int s=7, int d=13, float h=0.1): cdef DTYPE_t [:, ::1] integral = np.zeros_like(padded) 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. 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(-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) - 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] - @@ -370,14 +373,14 @@ def _fast_nl_means_denoising_2d(image, int s=7, int d=13, float h=0.1): distance /= (s2 * h2) if distance > 4: continue - weight = exp(- distance) + weight = alpha * exp(- distance) weights[x, y] += weight + weights[x + t1, y + t2] += weight result[x, y] += weight * padded[x + t1, y + t2] + result[x + t1, y + t2] += weight * padded[x, y] for x in range(offset, n_x - offset): - for y in range(offset, n_x - offset): + for y in range(offset, n_y - 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/tests/test_denoise.py b/skimage/restoration/tests/test_denoise.py index 04fcd184..409e2d3d 100644 --- a/skimage/restoration/tests/test_denoise.py +++ b/skimage/restoration/tests/test_denoise.py @@ -154,7 +154,7 @@ def test_nl_means_denoising_2d(): def test_fast_nl_means_denoising_2d(): - img = np.zeros((40, 40)) + 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) From c5e38dc5b1f2d3ae572734c304a5920f46829351 Mon Sep 17 00:00:00 2001 From: emmanuelle Date: Sat, 10 Jan 2015 19:11:15 +0100 Subject: [PATCH 11/30] Renamed module, since tests were breaking because module and function had the same name. --- skimage/restoration/{nl_means_denoising.py => non_local_means.py} | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename skimage/restoration/{nl_means_denoising.py => non_local_means.py} (100%) diff --git a/skimage/restoration/nl_means_denoising.py b/skimage/restoration/non_local_means.py similarity index 100% rename from skimage/restoration/nl_means_denoising.py rename to skimage/restoration/non_local_means.py From 8b9a777c79526bebdc75b5bce19c64f55ee94741 Mon Sep 17 00:00:00 2001 From: emmanuelle Date: Sun, 11 Jan 2015 22:37:42 +0100 Subject: [PATCH 12/30] Bug correction: I had forgotten to add __init__.py --- skimage/restoration/__init__.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/skimage/restoration/__init__.py b/skimage/restoration/__init__.py index a5c0ce0c..9f708852 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, fast_nl_means_denoising +from .non_local_means import nl_means_denoising, fast_nl_means_denoising __all__ = ['wiener', 'unsupervised_wiener', From dd9030d44c88f55f15901c3a38012bc7821a597b Mon Sep 17 00:00:00 2001 From: emmanuelle Date: Sun, 18 Jan 2015 22:18:29 +0100 Subject: [PATCH 13/30] 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() From 60d0c81aed7767d0f721970a1c6a0cc9b2f71205 Mon Sep 17 00:00:00 2001 From: emmanuelle Date: Tue, 20 Jan 2015 22:11:24 +0100 Subject: [PATCH 14/30] Changed type of error returned for image dimension strictly greater than 3. --- skimage/restoration/non_local_means.py | 4 ++-- skimage/restoration/tests/test_denoise.py | 2 +- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/skimage/restoration/non_local_means.py b/skimage/restoration/non_local_means.py index f727a12d..6dd20af6 100644 --- a/skimage/restoration/non_local_means.py +++ b/skimage/restoration/non_local_means.py @@ -120,6 +120,6 @@ def nl_means_denoising(image, patch_size=7, patch_distance=11, h=0.1, 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.") + raise NotImplementedError("Non-local means denoising is only \ + implemented 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 35e945ec..996c53e3 100644 --- a/skimage/restoration/tests/test_denoise.py +++ b/skimage/restoration/tests/test_denoise.py @@ -184,7 +184,7 @@ def test_nl_means_denoising_3d(): def test_nl_means_denoising_wrong_dimension(): img = np.zeros((5, 5, 5, 5)) - assert_raises(ValueError, restoration.nl_means_denoising, img) + assert_raises(NotImplementedError, restoration.nl_means_denoising, img) if __name__ == "__main__": From ec6aa2a017c0ce0a739fa3c083769ecfc2ca9a2a Mon Sep 17 00:00:00 2001 From: emmanuelle Date: Sun, 25 Jan 2015 11:49:07 +0100 Subject: [PATCH 15/30] Added docstring to helper functions in cython submodule --- skimage/restoration/_nl_means_denoising.pyx | 83 ++++++++++++++++++++- 1 file changed, 79 insertions(+), 4 deletions(-) diff --git a/skimage/restoration/_nl_means_denoising.pyx b/skimage/restoration/_nl_means_denoising.pyx index 1b8afc64..839f3039 100644 --- a/skimage/restoration/_nl_means_denoising.pyx +++ b/skimage/restoration/_nl_means_denoising.pyx @@ -13,6 +13,31 @@ cdef eps = 1.e-8 cdef inline float patch_distance_2d(DTYPE_t [:, :] p1, DTYPE_t [:, :] p2, DTYPE_t [:, ::] w, int s): + """ + Compute a Gaussian distance between two image patches. + + Parameters + ---------- + p1 : 2-D array_like + first patch + p2 : 2-D array_like + first 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 + + 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 @@ -36,6 +61,31 @@ cdef inline float patch_distance_2d(DTYPE_t [:, :] p1, cdef inline float patch_distance_2drgb(DTYPE_t [:, :, :] p1, DTYPE_t [:, :, :] p2, DTYPE_t [:, ::] 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 + first patch, 2D image with last dimension corresponding to channels + 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 + + exp( -w * (p1 - p2)**2) + """ cdef int i, j cdef int center = s / 2 cdef int color @@ -57,6 +107,31 @@ cdef inline float patch_distance_2drgb(DTYPE_t [:, :, :] p1, cdef inline float patch_distance_3d(DTYPE_t [:, :, :] p1, DTYPE_t [:, :, :] p2, DTYPE_t [:, :, ::] w, int s): + """ + Compute a Gaussian distance between two image patches. + + Parameters + ---------- + p1 : 3-D array_like + first patch + p2 : 3-D array_like + first patch + w : 3-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 + + exp( -w * (p1 - p2)**2) + """ cdef int i, j, k cdef float distance = 0 cdef float tmp_diff @@ -381,8 +456,8 @@ def _fast_nl_means_denoising_2d(image, int s=7, int d=13, float h=0.1): result[x + t1, y + t2] += weight * padded[x, y] for x in range(offset, n_x - offset): for y in range(offset, n_y - offset): - # I think there is no risk of division by zero - # except in padded zone + # No risk of division by zero, since the contribution + # of a null shift is strictly positive result[x, y] /= weights[x, y] return result[pad_size: - pad_size, pad_size: - pad_size] @@ -473,8 +548,8 @@ def _fast_nl_means_denoising_2drgb(image, int s=7, int d=13, float h=0.1): 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 + # No risk of division by zero, since the contribution + # of a null shift is strictly positive result[x, y, channel] /= weights[x, y] return result[pad_size: - pad_size, pad_size: - pad_size] From fa13744d59c50ac3ce65b4537806265b4fc795ef Mon Sep 17 00:00:00 2001 From: emmanuelle Date: Sun, 25 Jan 2015 11:49:45 +0100 Subject: [PATCH 16/30] Changed API of nl_means_denoising function to have a multichannel flag --- skimage/restoration/non_local_means.py | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/skimage/restoration/non_local_means.py b/skimage/restoration/non_local_means.py index 6dd20af6..cdd1a187 100644 --- a/skimage/restoration/non_local_means.py +++ b/skimage/restoration/non_local_means.py @@ -5,29 +5,29 @@ from skimage.restoration._nl_means_denoising import _nl_means_denoising_2d, \ _fast_nl_means_denoising_2drgb def nl_means_denoising(image, patch_size=7, patch_distance=11, h=0.1, - fast_mode=True): + 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 : ndarray - input data to be denoised - + 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. 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 @@ -105,14 +105,14 @@ def nl_means_denoising(image, patch_size=7, patch_distance=11, h=0.1, 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 + elif image.ndim == 3 and not multichannel: # 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 image.ndim == 3 and multichannel: # 2-D color (RGB) images if fast_mode: return np.array(_fast_nl_means_denoising_2drgb(image, patch_size, patch_distance, h)) From 065bc18253a49506f3e51651bfccf04aec87c613 Mon Sep 17 00:00:00 2001 From: emmanuelle Date: Sun, 25 Jan 2015 16:40:07 +0100 Subject: [PATCH 17/30] Merged 2D and 2D RGB functions for non-local means denoising --- skimage/restoration/_nl_means_denoising.pyx | 218 ++++---------------- skimage/restoration/non_local_means.py | 18 +- 2 files changed, 47 insertions(+), 189 deletions(-) diff --git a/skimage/restoration/_nl_means_denoising.pyx b/skimage/restoration/_nl_means_denoising.pyx index 839f3039..d7eac0f8 100644 --- a/skimage/restoration/_nl_means_denoising.pyx +++ b/skimage/restoration/_nl_means_denoising.pyx @@ -150,78 +150,6 @@ cdef inline float patch_distance_3d(DTYPE_t [:, :, :] p1, @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 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). 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 - # 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() - 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) * 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 @@ -242,13 +170,13 @@ def _nl_means_denoising_2drgb(image, int s=7, int d=13, float h=0.1): """ 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 n_x, n_y, n_ch + n_x, n_y, n_ch = 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] new_values = np.zeros(n_ch).astype(np.float32) cdef DTYPE_t [:, :, ::1] padded = np.ascontiguousarray(util.pad(image, ((offset, offset), (offset, offset), (0, 0)), mode='reflect').astype(np.float32)) @@ -261,13 +189,13 @@ 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. / (3 * np.sum(w) * h ** 2) * w + w = 1. / (n_ch * 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 x_end = x + offset + 1 for y in range(offset, n_y + offset): - for color in range(3): + for color in range(n_ch): new_values[color] = 0 weight_sum = 0 y_start = y - offset @@ -281,17 +209,26 @@ def _nl_means_denoising_2drgb(image, int s=7, int d=13, float h=0.1): min(d + 1, n_y - y - 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) + if n_ch == 1: + weight = patch_distance_2d( + padded[x_start: x_end, + y_start: y_end, 0], + padded[x_start_i: x_end_i, + y_start_j: y_end_j, 0], + w, s) + + else: + 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 - for color in range(3): + for color in range(n_ch): new_values[color] += weight * padded[x + i, y + j, color] - for color in range(3): + for color in range(n_ch): result[x, y, color] = new_values[color] / weight_sum return result[offset:-offset, offset:-offset] @@ -389,90 +326,7 @@ def _fast_nl_means_denoising_2d(image, int s=7, int d=13, float h=0.1): Parameters ---------- image: ndarray - 2-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, t1, t2, x, y - cdef float weight, distance - cdef float alpha - 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 - # 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) - 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 * h2) - # 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 - result[x, y] += weight * padded[x + t1, y + t2] - result[x + t1, y + t2] += weight * padded[x, y] - for x in range(offset, n_x - offset): - for y in range(offset, n_y - offset): - # No risk of division by zero, since the contribution - # of a null shift is strictly positive - 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 + 2-D input data to be denoised, grayscale or RGB s: int, optional size of patches used for denoising @@ -496,13 +350,13 @@ def _fast_nl_means_denoising_2drgb(image, int s=7, int d=13, float h=0.1): 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 int n_x, n_y, n_ch, 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, n_y, n_ch = image.shape + cdef float h2s2 = n_ch * h2 * s2 n_x += 2 * pad_size n_y += 2 * pad_size # Outer loops on patch shifts @@ -518,12 +372,16 @@ def _fast_nl_means_denoising_2drgb(image, int s=7, int d=13, float h=0.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) + if n_ch == 1: + distance = (padded[x, y, 0] - + padded[x + t1, y + t2, 0])**2 + else: + 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] @@ -542,12 +400,12 @@ def _fast_nl_means_denoising_2drgb(image, int s=7, int d=13, float h=0.1): weight = alpha * exp(- distance) weights[x, y] += weight weights[x + t1, y + t2] += weight - for ch in range(3): + for ch in range(n_ch): 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): + for channel in range(n_ch): # No risk of division by zero, since the contribution # of a null shift is strictly positive result[x, y, channel] /= weights[x, y] diff --git a/skimage/restoration/non_local_means.py b/skimage/restoration/non_local_means.py index cdd1a187..835ceed2 100644 --- a/skimage/restoration/non_local_means.py +++ b/skimage/restoration/non_local_means.py @@ -1,8 +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, \ - _fast_nl_means_denoising_2d, _fast_nl_means_denoising_3d, \ - _fast_nl_means_denoising_2drgb + _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): @@ -99,12 +98,13 @@ def nl_means_denoising(image, patch_size=7, patch_distance=11, h=0.1, >>> denoised_a = nl_means_denoising(a, 7, 5, 0.1) """ if image.ndim == 2: + image = image[..., np.newaxis] if fast_mode: - return np.array(_fast_nl_means_denoising_2d(image, s=patch_size, - d=patch_distance, h=h)) + return np.squeeze(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)) + return np.squeeze(np.array(_nl_means_denoising_2d(image, + s=patch_size, d=patch_distance, h=h))) elif image.ndim == 3 and not multichannel: # only grayscale if fast_mode: return np.array(_fast_nl_means_denoising_3d(image, s=patch_size, @@ -114,10 +114,10 @@ def nl_means_denoising(image, patch_size=7, patch_distance=11, h=0.1, patch_distance, h)) if image.ndim == 3 and multichannel: # 2-D color (RGB) images if fast_mode: - return np.array(_fast_nl_means_denoising_2drgb(image, patch_size, + return np.array(_fast_nl_means_denoising_2d(image, patch_size, patch_distance, h)) else: - return np.array(_nl_means_denoising_2drgb(image, patch_size, + return np.array(_nl_means_denoising_2d(image, patch_size, patch_distance, h)) else: raise NotImplementedError("Non-local means denoising is only \ From e83873d56383835755378a23141e94bebd0922e3 Mon Sep 17 00:00:00 2001 From: emmanuelle Date: Sun, 25 Jan 2015 16:42:39 +0100 Subject: [PATCH 18/30] Removed blank lines in docstring --- skimage/restoration/_nl_means_denoising.pyx | 12 ------------ 1 file changed, 12 deletions(-) diff --git a/skimage/restoration/_nl_means_denoising.pyx b/skimage/restoration/_nl_means_denoising.pyx index d7eac0f8..b603a6a9 100644 --- a/skimage/restoration/_nl_means_denoising.pyx +++ b/skimage/restoration/_nl_means_denoising.pyx @@ -157,13 +157,10 @@ def _nl_means_denoising_2d(image, int s=7, int d=13, float h=0.1): ---------- 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. @@ -244,13 +241,10 @@ def _nl_means_denoising_3d(image, int s=7, ---------- 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) """ @@ -327,13 +321,10 @@ def _fast_nl_means_denoising_2d(image, int s=7, int d=13, float h=0.1): ---------- 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. @@ -423,13 +414,10 @@ def _fast_nl_means_denoising_3d(image, int s=5, int d=7, float h=0.1): ---------- 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. From 2e0b37b11d5d2a4873558c4e44a459b8dd4c69d3 Mon Sep 17 00:00:00 2001 From: emmanuelle Date: Sun, 25 Jan 2015 19:21:23 +0100 Subject: [PATCH 19/30] Changed variable names for coordinates for better consistency with skimage conventions. Moved some code to inline functions to improve clarity. --- skimage/restoration/_nl_means_denoising.pyx | 387 +++++++++++--------- skimage/restoration/non_local_means.py | 2 +- skimage/restoration/tests/test_denoise.py | 21 +- 3 files changed, 239 insertions(+), 171 deletions(-) diff --git a/skimage/restoration/_nl_means_denoising.pyx b/skimage/restoration/_nl_means_denoising.pyx index b603a6a9..832efa85 100644 --- a/skimage/restoration/_nl_means_denoising.pyx +++ b/skimage/restoration/_nl_means_denoising.pyx @@ -167,12 +167,12 @@ def _nl_means_denoising_2d(image, int s=7, int d=13, float h=0.1): """ if s % 2 == 0: s += 1 # odd value for symmetric patch - cdef int n_x, n_y, n_ch - n_x, n_y, n_ch = image.shape + cdef int n_row, n_col, n_ch + n_row, n_col, n_ch = 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 int x_row, x_col, i, j, color + cdef int x_row_start, x_row_end, x_col_start, x_col_end + cdef int x_row_start_i, x_row_end_i, x_col_start_j, x_col_end_j cdef DTYPE_t [::1] new_values = np.zeros(n_ch).astype(np.float32) cdef DTYPE_t [:, :, ::1] padded = np.ascontiguousarray(util.pad(image, ((offset, offset), (offset, offset), (0, 0)), @@ -181,52 +181,52 @@ def _nl_means_denoising_2d(image, int s=7, int d=13, float h=0.1): 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] + xg_row, xg_col = 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)) + - (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 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): + for x_row in range(offset, n_row + offset): + x_row_start = x_row - offset + x_row_end = x_row + offset + 1 + for x_col in range(offset, n_col + offset): for color in range(n_ch): new_values[color] = 0 weight_sum = 0 - y_start = y - offset - y_end = y + offset + 1 + x_col_start = x_col - offset + x_col_end = x_col + 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 + for i in range(max(- d, offset - x_row), + min(d + 1, n_row - x_row - 1)): + x_row_start_i = x_row_start + i + x_row_end_i = x_row_end + i + for j in range(max(- d, offset - x_col), + min(d + 1, n_col - x_col - 1)): + x_col_start_j = x_col_start + j + x_col_end_j = x_col_end + j if n_ch == 1: weight = patch_distance_2d( - padded[x_start: x_end, - y_start: y_end, 0], - padded[x_start_i: x_end_i, - y_start_j: y_end_j, 0], + padded[x_row_start: x_row_end, + x_col_start: x_col_end, 0], + padded[x_row_start_i: x_row_end_i, + x_col_start_j: x_col_end_j, 0], w, s) else: 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, :], + padded[x_row_start: x_row_end, + x_col_start: x_col_end, :], + padded[x_row_start_i: x_row_end_i, + x_col_start_j: x_col_end_j, :], w, s) weight_sum += weight for color in range(n_ch): - new_values[color] += weight * padded[x + i, y + j, - color] + new_values[color] += weight * padded[x_row + i, + x_col + j, color] for color in range(n_ch): - result[x, y, color] = new_values[color] / weight_sum + result[x_row, x_col, color] = new_values[color] / weight_sum return result[offset:-offset, offset:-offset] @@ -250,8 +250,8 @@ def _nl_means_denoising_3d(image, int s=7, """ if s % 2 == 0: s += 1 # odd value for symmetric patch - cdef int n_x, n_y, n_z - n_x, n_y, n_z = image.shape + 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 DTYPE_t [:, :, ::1] padded = np.ascontiguousarray(util.pad( @@ -261,54 +261,104 @@ 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, - -offset: offset + 1] + xg_pln, xg_row, xg_col = 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)). - astype(np.float32)) + - (xg_pln ** 2 + xg_row ** 2 + xg_col ** 2) / + (2 * A ** 2)).astype(np.float32)) cdef float distance - cdef int x, y, z, i, j, k - 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 + cdef int x_pln, x_row, x_col, i, j, k + cdef int x_pln_start, x_pln_end, x_row_start, x_row_end, \ + x_col_start, x_col_end + cdef int x_pln_start_i, x_pln_end_i, x_row_start_j, x_row_end_j, \ + x_col_start_k, x_col_end_k w = 1. / (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 - 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 + for x_pln in range(offset, n_pln + offset): + x_pln_start = x_pln - offset + x_pln_end = x_pln + offset + 1 + for x_row in range(offset, n_row + offset): + x_row_start = x_row - offset + x_row_end = x_row + offset + 1 + for x_col in range(offset, n_col + offset): + x_col_start = x_col - offset + x_col_end = x_col + 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)): - z_start_k = z_start + k - z_end_k = z_end + k + for i in range(max(- d, offset - x_pln), + min(d + 1, n_pln - x_pln - 1)): + x_pln_start_i = x_pln_start + i + x_pln_end_i = x_pln_end + i + for j in range(max(- d, offset - x_row), + min(d + 1, n_row - x_row - 1)): + x_row_start_j = x_row_start + j + x_row_end_j = x_row_end + j + for k in range(max(- d, offset - x_col), + min(d + 1, n_col - x_col - 1)): + x_col_start_k = x_col_start + k + x_col_end_k = x_col_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], + padded[x_pln_start: x_pln_end, + x_row_start: x_row_end, + x_col_start: x_col_end], + padded[x_pln_start_i: x_pln_end_i, + x_row_start_j: x_row_end_j, + x_col_start_k: x_col_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 + new_value += weight * padded[x_pln + i, + x_row + j, x_col + k] + result[x_pln, x_row, x_col] = new_value / weight_sum 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(DTYPE_t [:, ::] integral, + int x_row, int x_col, int offset, float h2s2): + """ + See + 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[x_row + offset, x_col + offset] + \ + integral[x_row - offset, x_col - offset] - \ + integral[x_row - offset, x_col + offset] - \ + integral[x_row + offset, x_col - offset] + distance /= h2s2 + return distance + +@cython.cdivision(True) +@cython.boundscheck(False) +cdef inline float _integral_to_distance_3d(DTYPE_t [:, :, ::] integral, + int x_pln, int x_row, int x_col, int offset, + float s_cube_h_square): + """ + See + 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[x_pln + offset, x_row + offset, x_col + offset] + - integral[x_pln - offset, x_row - offset, x_col - offset] + + integral[x_pln - offset, x_row - offset, x_col + offset] + + integral[x_pln - offset, x_row + offset, x_col - offset] + + integral[x_pln + offset, x_row - offset, x_col - offset] + - integral[x_pln - offset, x_row + offset, x_col + offset] + - integral[x_pln + offset, x_row - offset, x_col + offset] + - integral[x_pln + offset, x_row + offset, x_col - offset]) + distance /= s_cube_h_square + return distance + @cython.cdivision(True) @cython.boundscheck(False) @@ -341,65 +391,66 @@ def _fast_nl_means_denoising_2d(image, int s=7, int d=13, float h=0.1): 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, n_ch, t1, t2, x, y + cdef int n_row, n_col, n_ch, t_row, t_col, x_row, x_col cdef float weight, distance cdef float alpha cdef float h2 = h ** 2. cdef float s2 = s ** 2. - n_x, n_y, n_ch = image.shape + n_row, n_col, n_ch = image.shape cdef float h2s2 = n_ch * h2 * s2 - n_x += 2 * pad_size - n_y += 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 - for t1 in range(-d, d + 1): - for t2 in range(0, d + 1): + for t_row in range(-d, d + 1): + 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 t2 == 0 and t1 is not 0: + if t_col == 0 and t_row 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)): + for x_row in range(max(1, - t_row), min(n_row, n_row - t_row)): + for x_col in range(max(1, - t_col), min(n_col, n_col - t_col)): if n_ch == 1: - distance = (padded[x, y, 0] - - padded[x + t1, y + t2, 0])**2 + distance = (padded[x_row, x_col, 0] - + padded[x_row + t_row, x_col + t_col, 0])**2 else: - 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 + distance = ((padded[x_row, x_col, 0] - + padded[x_row + t_row, x_col + t_col, 0])**2 + +(padded[x_row, x_col, 1] - + padded[x_row + t_row, x_col + t_col, 1])**2 + +(padded[x_row, x_col, 2] - + padded[x_row + t_row, x_col + t_col, 2])**2) + integral[x_row, x_col] = distance + \ + integral[x_row - 1, x_col] + \ + integral[x_row, x_col - 1] \ + - integral[x_row - 1, x_col - 1] + for x_row in range(max(offset, offset - t_row), + min(n_row - offset, n_row - offset - t_row)): + for x_col in range(max(offset, offset - t_col), + min(n_col - offset, n_col - offset - t_col)): + distance = _integral_to_distance_2d(integral, x_row, x_col, + offset, 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 + weights[x_row, x_col] += weight + weights[x_row + t_row, x_col + t_col] += weight for ch in range(n_ch): - 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): + result[x_row, x_col, ch] += weight * \ + padded[x_row + t_row, x_col + t_col, ch] + result[x_row + t_row, x_col + t_col, ch] += \ + weight * padded[x_row, x_col, ch] + # Normalize pixel values using sum of weights of contributing patches + for x_row in range(offset, n_row - offset): + for x_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[x, y, channel] /= weights[x, y] + result[x_row, x_col, channel] /= weights[x_row, x_col] return result[pad_size: - pad_size, pad_size: - pad_size] @@ -433,85 +484,85 @@ def _fast_nl_means_denoising_3d(image, int s=5, int d=7, float h=0.1): 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 int n_pln, n_row, n_col, t_pln, t_row, t_col, \ + x_pln, x_row, x_col + cdef int x_pln_integral_min, x_pln_integral_max, \ + x_row_integral_min, x_row_integral_max, \ + x_col_integral_min, x_col_integral_max + cdef int x_pln_dist_min, x_pln_dist_max, x_row_dist_min, x_row_dist_max, \ + x_col_dist_min, x_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_x, n_y, n_z = image.shape - n_x += 2 * pad_size - n_y += 2 * pad_size - n_z += 2 * pad_size + 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 - 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) + for t_pln in range(-d, d + 1): + x_pln_integral_min = max(1, - t_pln) + x_pln_integral_max = min(n_pln, n_pln - t_pln) + x_pln_dist_min = max(offset, offset - t_pln) + x_pln_dist_max = min(n_pln - offset, n_pln - offset - t_pln) + for t_row in range(-d, d + 1): + x_row_integral_min = max(1, - t_row) + x_row_integral_max = min(n_row, n_row - t_row) + x_row_dist_min = max(offset, offset - t_row) + x_row_dist_max = min(n_row - offset, n_row - offset - t_row) + for t_col in range(0, d + 1): + x_col_integral_min = max(1, - t_col) + x_col_integral_max = min(n_col, n_col - t_col) + x_col_dist_min = max(offset, offset - t_col) + x_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 t3 == 0 and (t1 is not 0 or t2 is not 0): + if t_col == 0 and (t_pln is not 0 or t_row 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 + for x_pln in range(x_pln_integral_min, x_pln_integral_max): + for x_row in range(x_row_integral_min, x_row_integral_max): + for x_col in range(x_col_integral_min, + x_col_integral_max): + integral[x_pln, x_row, x_col] = \ + ((padded[x_pln, x_row, x_col] - + padded[x_pln + t_pln, x_row + t_row, + x_col + t_col])**2 + + integral[x_pln - 1, x_row, x_col] + + integral[x_pln, x_row - 1, x_col] + + integral[x_pln, x_row, x_col - 1] + + integral[x_pln - 1, x_row - 1, x_col - 1] + - integral[x_pln - 1, x_row - 1, x_col] + - integral[x_pln, x_row - 1, x_col - 1] + - integral[x_pln - 1, x_row, x_col - 1]) + for x_pln in range(x_pln_dist_min, x_pln_dist_max): + for x_row in range(x_row_dist_min, x_row_dist_max): + for x_col in range(x_col_dist_min, x_col_dist_max): + distance = _integral_to_distance_3d(integral, + x_pln, x_row, x_col, offset, + 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): + weights[x_pln, x_row, x_col] += weight + weights[x_pln + t_pln, x_row + t_row, + x_col + t_col] += weight + result[x_pln, x_row, x_col] += weight * \ + padded[x_pln + t_pln, x_row + t_row, + x_col + t_col] + result[x_pln + t_pln, x_row + t_row, + x_col + t_col] += weight * \ + padded[x_pln, x_row, x_col] + for x_pln in range(offset, n_pln - offset): + for x_row in range(offset, n_row - offset): + for x_col in range(offset, n_col - offset): # I think there is no risk of division by zero # except in padded zone - result[x, y, z] /= weights[x, y, z] + result[x_pln, x_row, x_col] /= weights[x_pln, x_row, x_col] 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 835ceed2..13d889a7 100644 --- a/skimage/restoration/non_local_means.py +++ b/skimage/restoration/non_local_means.py @@ -26,7 +26,7 @@ def nl_means_denoising(image, patch_size=7, patch_distance=11, h=0.1, 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. + 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 diff --git a/skimage/restoration/tests/test_denoise.py b/skimage/restoration/tests/test_denoise.py index 996c53e3..934e10d4 100644 --- a/skimage/restoration/tests/test_denoise.py +++ b/skimage/restoration/tests/test_denoise.py @@ -174,14 +174,31 @@ 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, fast_mode=True) + denoised = restoration.nl_means_denoising(img, 5, 4, 0.1, fast_mode=True, + multichannel=False) # make sure noise is reduced assert img.std() > denoised.std() - denoised = restoration.nl_means_denoising(img, 5, 4, 0.1, fast_mode=False) + denoised = restoration.nl_means_denoising(img, 5, 4, 0.1, 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) From 4697d66c85524de6255e6dc5ad0bc079f1d29d26 Mon Sep 17 00:00:00 2001 From: emmanuelle Date: Mon, 26 Jan 2015 22:06:38 +0100 Subject: [PATCH 20/30] [ENH] Removed eps parameter that was not needed. Added a test for the case when no denoising at all is performed when h is very small. Corrected a bug in looping range in the classical algorithm. --- skimage/restoration/_nl_means_denoising.pyx | 23 +++++++++------------ skimage/restoration/tests/test_denoise.py | 23 +++++++++++++++------ 2 files changed, 27 insertions(+), 19 deletions(-) diff --git a/skimage/restoration/_nl_means_denoising.pyx b/skimage/restoration/_nl_means_denoising.pyx index 832efa85..7d0c1fcf 100644 --- a/skimage/restoration/_nl_means_denoising.pyx +++ b/skimage/restoration/_nl_means_denoising.pyx @@ -6,8 +6,6 @@ from libc.math cimport exp ctypedef np.float32_t DTYPE_t -cdef eps = 1.e-8 - @cython.boundscheck(False) cdef inline float patch_distance_2d(DTYPE_t [:, :] p1, @@ -44,15 +42,15 @@ cdef inline float patch_distance_2d(DTYPE_t [:, :] p1, 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 + 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 > 5: - return eps + return 0. for j in range(s): tmp_diff = p1[i, j] - p2[i, j] - distance += w[i, j] * tmp_diff * tmp_diff + distance += (w[i, j] * tmp_diff * tmp_diff) distance = exp(- distance) return distance @@ -94,7 +92,7 @@ cdef inline float patch_distance_2drgb(DTYPE_t [:, :, :] p1, for i in range(s): # exp of large negative numbers will be 0, so we'd better stop if distance > 5: - return eps + return 0. for j in range(s): for color in range(3): tmp_diff = p1[i, j, color] - p2[i, j, color] @@ -138,7 +136,7 @@ cdef inline float patch_distance_3d(DTYPE_t [:, :, :] p1, for i in range(s): # exp of large negative numbers will be 0, so we'd better stop if distance > 5: - return eps + return 0. for j in range(s): for k in range(s): tmp_diff = p1[i, j, k] - p2[i, j, k] @@ -199,11 +197,11 @@ def _nl_means_denoising_2d(image, int s=7, int d=13, float h=0.1): x_col_end = x_col + offset + 1 # Coordinates of test pixel and patch bounds for i in range(max(- d, offset - x_row), - min(d + 1, n_row - x_row - 1)): + min(d + 1, n_row + offset - x_row)): x_row_start_i = x_row_start + i x_row_end_i = x_row_end + i for j in range(max(- d, offset - x_col), - min(d + 1, n_col - x_col - 1)): + min(d + 1, n_col + offset - x_col)): x_col_start_j = x_col_start + j x_col_end_j = x_col_end + j if n_ch == 1: @@ -213,7 +211,6 @@ def _nl_means_denoising_2d(image, int s=7, int d=13, float h=0.1): padded[x_row_start_i: x_row_end_i, x_col_start_j: x_col_end_j, 0], w, s) - else: weight = patch_distance_2drgb( padded[x_row_start: x_row_end, @@ -288,15 +285,15 @@ def _nl_means_denoising_3d(image, int s=7, weight_sum = 0 # Coordinates of test pixel and patch bounds for i in range(max(- d, offset - x_pln), - min(d + 1, n_pln - x_pln - 1)): + min(d + 1, n_pln + offset - x_pln)): x_pln_start_i = x_pln_start + i x_pln_end_i = x_pln_end + i for j in range(max(- d, offset - x_row), - min(d + 1, n_row - x_row - 1)): + min(d + 1, n_row + offset - x_row)): x_row_start_j = x_row_start + j x_row_end_j = x_row_end + j for k in range(max(- d, offset - x_col), - min(d + 1, n_col - x_col - 1)): + min(d + 1, n_col + offset - x_col)): x_col_start_k = x_col_start + k x_col_end_k = x_col_end + k weight = patch_distance_3d( diff --git a/skimage/restoration/tests/test_denoise.py b/skimage/restoration/tests/test_denoise.py index 934e10d4..e84e1c49 100644 --- a/skimage/restoration/tests/test_denoise.py +++ b/skimage/restoration/tests/test_denoise.py @@ -148,10 +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, fast_mode=True) + 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.1, fast_mode=False) + denoised = restoration.nl_means_denoising(img, 7, 5, 0.2, fast_mode=False) # make sure noise is reduced assert img.std() > denoised.std() @@ -162,10 +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, fast_mode=True) + 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.08, fast_mode=False) + denoised = restoration.nl_means_denoising(img, 7, 9, 0.3, fast_mode=False) # make sure noise is reduced assert img.std() > denoised.std() @@ -174,11 +174,11 @@ 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, fast_mode=True, + 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.1, fast_mode=False, + 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() @@ -204,5 +204,16 @@ def test_nl_means_denoising_wrong_dimension(): 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() From 112b0b4fcd0d8e043a32c77aee546cf6811d3e24 Mon Sep 17 00:00:00 2001 From: emmanuelle Date: Thu, 29 Jan 2015 21:36:08 +0100 Subject: [PATCH 21/30] Global constant defined for distance cutoff --- skimage/restoration/_nl_means_denoising.pyx | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/skimage/restoration/_nl_means_denoising.pyx b/skimage/restoration/_nl_means_denoising.pyx index 7d0c1fcf..70b60921 100644 --- a/skimage/restoration/_nl_means_denoising.pyx +++ b/skimage/restoration/_nl_means_denoising.pyx @@ -6,6 +6,7 @@ from libc.math cimport exp ctypedef np.float32_t DTYPE_t +cdef float DISTANCE_CUTOFF = 5. @cython.boundscheck(False) cdef inline float patch_distance_2d(DTYPE_t [:, :] p1, @@ -46,7 +47,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 > 5: + if distance > DISTANCE_CUTOFF: return 0. for j in range(s): tmp_diff = p1[i, j] - p2[i, j] @@ -91,7 +92,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 > 5: + if distance > DISTANCE_CUTOFF: return 0. for j in range(s): for color in range(3): @@ -135,7 +136,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 > 5: + if distance > DISTANCE_CUTOFF: return 0. for j in range(s): for k in range(s): @@ -431,7 +432,7 @@ def _fast_nl_means_denoising_2d(image, int s=7, int d=13, float h=0.1): distance = _integral_to_distance_2d(integral, x_row, x_col, offset, h2s2) # exp of large negative numbers is close to zero - if distance > 5: + if distance > DISTANCE_CUTOFF: continue weight = alpha * exp(- distance) weights[x_row, x_col] += weight @@ -543,7 +544,7 @@ def _fast_nl_means_denoising_3d(image, int s=5, int d=7, float h=0.1): x_pln, x_row, x_col, offset, s_cube_h_square) # exp of large negative numbers is close to zero - if distance > 5.: + if distance > DISTANCE_CUTOFF: continue weight = alpha * exp(- distance) weights[x_pln, x_row, x_col] += weight From 009e76141137f948301da7938a15a6ee2aa65419 Mon Sep 17 00:00:00 2001 From: emmanuelle Date: Thu, 29 Jan 2015 22:40:46 +0100 Subject: [PATCH 22/30] [ENH] Docstring cleaning --- skimage/restoration/_nl_means_denoising.pyx | 134 ++++++++++++-------- skimage/restoration/non_local_means.py | 12 +- 2 files changed, 84 insertions(+), 62 deletions(-) diff --git a/skimage/restoration/_nl_means_denoising.pyx b/skimage/restoration/_nl_means_denoising.pyx index 70b60921..1854fad7 100644 --- a/skimage/restoration/_nl_means_denoising.pyx +++ b/skimage/restoration/_nl_means_denoising.pyx @@ -18,13 +18,13 @@ cdef inline float patch_distance_2d(DTYPE_t [:, :] p1, Parameters ---------- p1 : 2-D array_like - first patch + First patch. p2 : 2-D array_like - first patch + Second patch. w : 2-D array_like - array of weigths for the different pixels of the patches + Array of weigths for the different pixels of the patches. s : int - linear size of the patches + Linear size of the patches. Returns ------- @@ -35,7 +35,7 @@ cdef inline float patch_distance_2d(DTYPE_t [:, :] p1, ----- The returned distance is given by - exp( -w * (p1 - p2)**2) + .. math:: \exp( -w (p1 - p2)^2) """ cdef int i, j cdef int center = s / 2 @@ -52,7 +52,7 @@ cdef inline float patch_distance_2d(DTYPE_t [:, :] p1, for j in range(s): tmp_diff = p1[i, j] - p2[i, j] distance += (w[i, j] * tmp_diff * tmp_diff) - distance = exp(- distance) + distance = exp(-distance) return distance @@ -66,13 +66,13 @@ cdef inline float patch_distance_2drgb(DTYPE_t [:, :, :] p1, Parameters ---------- p1 : 3-D array_like - first patch, 2D image with last dimension corresponding to channels + First patch, 2D image with last dimension corresponding to channels. p2 : 3-D array_like - first patch, 2D image with last dimension corresponding to channels + Second patch, 2D image with last dimension corresponding to channels. w : 2-D array_like - array of weigths for the different pixels of the patches + Array of weights for the different pixels of the patches. s : int - linear size of the patches + Linear size of the patches. Returns ------- @@ -83,7 +83,7 @@ cdef inline float patch_distance_2drgb(DTYPE_t [:, :, :] p1, ----- The returned distance is given by - exp( -w * (p1 - p2)**2) + .. math:: \exp( -w (p1 - p2)^2) """ cdef int i, j cdef int center = s / 2 @@ -98,7 +98,7 @@ cdef inline float patch_distance_2drgb(DTYPE_t [:, :, :] p1, 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) + distance = exp(-distance) return distance @@ -112,13 +112,13 @@ cdef inline float patch_distance_3d(DTYPE_t [:, :, :] p1, Parameters ---------- p1 : 3-D array_like - first patch + First patch. p2 : 3-D array_like - first patch + Second patch. w : 3-D array_like - array of weigths for the different pixels of the patches + Array of weights for the different pixels of the patches. s : int - linear size of the patches + Linear size of the patches. Returns ------- @@ -129,7 +129,7 @@ cdef inline float patch_distance_3d(DTYPE_t [:, :, :] p1, ----- The returned distance is given by - exp( -w * (p1 - p2)**2) + .. math:: \exp( -w (p1 - p2)^2) """ cdef int i, j, k cdef float distance = 0 @@ -142,7 +142,7 @@ cdef inline float patch_distance_3d(DTYPE_t [:, :, :] p1, 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) + distance = exp(-distance) return distance @@ -154,15 +154,20 @@ def _nl_means_denoising_2d(image, int s=7, int d=13, float h=0.1): 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 + 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 @@ -182,8 +187,8 @@ def _nl_means_denoising_2d(image, int s=7, int d=13, float h=0.1): cdef float weight_sum, weight xg_row, xg_col = np.mgrid[-offset:offset + 1, -offset:offset + 1] cdef DTYPE_t [:, ::1] w = np.ascontiguousarray(np.exp( - - (xg_row ** 2 + xg_col ** 2) / (2 * A ** 2)). - astype(np.float32)) + -(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 and patch bounds @@ -237,14 +242,19 @@ def _nl_means_denoising_3d(image, int s=7, 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) + 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 @@ -263,7 +273,7 @@ def _nl_means_denoising_3d(image, int s=7, -offset: offset + 1, -offset: offset + 1] cdef DTYPE_t [:, :, ::1] w = np.ascontiguousarray(np.exp( - - (xg_pln ** 2 + xg_row ** 2 + xg_col ** 2) / + -(xg_pln ** 2 + xg_row ** 2 + xg_col ** 2) / (2 * A ** 2)).astype(np.float32)) cdef float distance cdef int x_pln, x_row, x_col, i, j, k @@ -318,7 +328,8 @@ def _nl_means_denoising_3d(image, int s=7, cdef inline float _integral_to_distance_2d(DTYPE_t [:, ::] integral, int x_row, int x_col, int offset, float h2s2): """ - See + References + ---------- Jacques Froment. Parameter-Free Fast Pixelwise Non-Local Means Denoising. Image Processing On Line, 2014, vol. 4, p. 300-326. @@ -338,7 +349,8 @@ cdef inline float _integral_to_distance_3d(DTYPE_t [:, :, ::] integral, int x_pln, int x_row, int x_col, int offset, float s_cube_h_square): """ - See + References + ---------- Jacques Froment. Parameter-Free Fast Pixelwise Non-Local Means Denoising. Image Processing On Line, 2014, vol. 4, p. 300-326. @@ -367,15 +379,20 @@ def _fast_nl_means_denoising_2d(image, int s=7, int d=13, float h=0.1): 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 + 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 @@ -434,7 +451,7 @@ def _fast_nl_means_denoising_2d(image, int s=7, int d=13, float h=0.1): # exp of large negative numbers is close to zero if distance > DISTANCE_CUTOFF: continue - weight = alpha * exp(- distance) + weight = alpha * exp(-distance) weights[x_row, x_col] += weight weights[x_row + t_row, x_col + t_col] += weight for ch in range(n_ch): @@ -461,15 +478,20 @@ def _fast_nl_means_denoising_3d(image, int s=5, int d=7, float h=0.1): 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 + 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 @@ -546,7 +568,7 @@ def _fast_nl_means_denoising_3d(image, int s=5, int d=7, float h=0.1): # exp of large negative numbers is close to zero if distance > DISTANCE_CUTOFF: continue - weight = alpha * exp(- distance) + weight = alpha * exp(-distance) weights[x_pln, x_row, x_col] += weight weights[x_pln + t_pln, x_row + t_row, x_col + t_col] += weight diff --git a/skimage/restoration/non_local_means.py b/skimage/restoration/non_local_means.py index 13d889a7..a94437c5 100644 --- a/skimage/restoration/non_local_means.py +++ b/skimage/restoration/non_local_means.py @@ -12,14 +12,14 @@ def nl_means_denoising(image, patch_size=7, patch_distance=11, h=0.1, Parameters ---------- image : 2D or 3D ndarray - input image to be denoised, which can be 2D or 3D, and grayscale + 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 + Size of patches used for denoising. patch_distance : int, optional - maximal distance in pixels where to search patches used for denoising + 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 + 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 @@ -28,7 +28,7 @@ def nl_means_denoising(image, patch_size=7, patch_distance=11, h=0.1, 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 + 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. @@ -36,7 +36,7 @@ def nl_means_denoising(image, patch_size=7, patch_distance=11, h=0.1, ------- result : ndarray - denoised image, of same shape as `image`. + Denoised image, of same shape as `image`. See Also -------- From 4dbf33b9ca6470a7308b431f95df2af897523c9e Mon Sep 17 00:00:00 2001 From: emmanuelle Date: Sun, 1 Feb 2015 20:33:57 +0100 Subject: [PATCH 23/30] PEP 8: spaces and operators --- skimage/restoration/_nl_means_denoising.pyx | 102 ++++++++++---------- 1 file changed, 51 insertions(+), 51 deletions(-) diff --git a/skimage/restoration/_nl_means_denoising.pyx b/skimage/restoration/_nl_means_denoising.pyx index 1854fad7..e2db20e9 100644 --- a/skimage/restoration/_nl_means_denoising.pyx +++ b/skimage/restoration/_nl_means_denoising.pyx @@ -202,27 +202,27 @@ def _nl_means_denoising_2d(image, int s=7, int d=13, float h=0.1): x_col_start = x_col - offset x_col_end = x_col + offset + 1 # Coordinates of test pixel and patch bounds - for i in range(max(- d, offset - x_row), + for i in range(max(-d, offset - x_row), min(d + 1, n_row + offset - x_row)): x_row_start_i = x_row_start + i x_row_end_i = x_row_end + i - for j in range(max(- d, offset - x_col), + for j in range(max(-d, offset - x_col), min(d + 1, n_col + offset - x_col)): x_col_start_j = x_col_start + j x_col_end_j = x_col_end + j if n_ch == 1: weight = patch_distance_2d( - padded[x_row_start: x_row_end, - x_col_start: x_col_end, 0], - padded[x_row_start_i: x_row_end_i, - x_col_start_j: x_col_end_j, 0], + padded[x_row_start:x_row_end, + x_col_start:x_col_end, 0], + padded[x_row_start_i:x_row_end_i, + x_col_start_j:x_col_end_j, 0], w, s) else: weight = patch_distance_2drgb( - padded[x_row_start: x_row_end, - x_col_start: x_col_end, :], - padded[x_row_start_i: x_row_end_i, - x_col_start_j: x_col_end_j, :], + padded[x_row_start:x_row_end, + x_col_start:x_col_end, :], + padded[x_row_start_i:x_row_end_i, + x_col_start_j:x_col_end_j, :], w, s) weight_sum += weight for color in range(n_ch): @@ -295,25 +295,25 @@ def _nl_means_denoising_3d(image, int s=7, new_value = 0 weight_sum = 0 # Coordinates of test pixel and patch bounds - for i in range(max(- d, offset - x_pln), + for i in range(max(-d, offset - x_pln), min(d + 1, n_pln + offset - x_pln)): x_pln_start_i = x_pln_start + i x_pln_end_i = x_pln_end + i - for j in range(max(- d, offset - x_row), + for j in range(max(-d, offset - x_row), min(d + 1, n_row + offset - x_row)): x_row_start_j = x_row_start + j x_row_end_j = x_row_end + j - for k in range(max(- d, offset - x_col), + for k in range(max(-d, offset - x_col), min(d + 1, n_col + offset - x_col)): x_col_start_k = x_col_start + k x_col_end_k = x_col_end + k weight = patch_distance_3d( - padded[x_pln_start: x_pln_end, - x_row_start: x_row_end, - x_col_start: x_col_end], - padded[x_pln_start_i: x_pln_end_i, - x_row_start_j: x_row_end_j, - x_col_start_k: x_col_end_k], + padded[x_pln_start:x_pln_end, + x_row_start:x_row_end, + x_col_start:x_col_end], + padded[x_pln_start_i:x_pln_end_i, + x_row_start_j:x_row_end_j, + x_col_start_k:x_col_end_k], w, s) weight_sum += weight new_value += weight * padded[x_pln + i, @@ -323,6 +323,7 @@ def _nl_means_denoising_3d(image, int s=7, #-------------- Accelerated algorithm of Froment 2015 ------------------ + @cython.cdivision(True) @cython.boundscheck(False) cdef inline float _integral_to_distance_2d(DTYPE_t [:, ::] integral, @@ -343,6 +344,7 @@ cdef inline float _integral_to_distance_2d(DTYPE_t [:, ::] integral, distance /= h2s2 return distance + @cython.cdivision(True) @cython.boundscheck(False) cdef inline float _integral_to_distance_3d(DTYPE_t [:, :, ::] integral, @@ -356,16 +358,15 @@ cdef inline float _integral_to_distance_3d(DTYPE_t [:, :, ::] integral, Used in _fast_nl_means_denoising_3d """ - cdef float distance - distance = (integral[x_pln + offset, x_row + offset, x_col + offset] - - integral[x_pln - offset, x_row - offset, x_col - offset] - + integral[x_pln - offset, x_row - offset, x_col + offset] - + integral[x_pln - offset, x_row + offset, x_col - offset] - + integral[x_pln + offset, x_row - offset, x_col - offset] - - integral[x_pln - offset, x_row + offset, x_col + offset] - - integral[x_pln + offset, x_row - offset, x_col + offset] - - integral[x_pln + offset, x_row + offset, x_col - offset]) + distance = (integral[x_pln + offset, x_row + offset, x_col + offset] - + integral[x_pln - offset, x_row - offset, x_col - offset] + + integral[x_pln - offset, x_row - offset, x_col + offset] + + integral[x_pln - offset, x_row + offset, x_col - offset] + + integral[x_pln + offset, x_row - offset, x_col - offset] - + integral[x_pln - offset, x_row + offset, x_col + offset] - + integral[x_pln + offset, x_row - offset, x_col + offset] - + integral[x_pln + offset, x_row + offset, x_col - offset]) distance /= s_cube_h_square return distance @@ -426,22 +427,22 @@ def _fast_nl_means_denoising_2d(image, int s=7, int d=13, float h=0.1): else: alpha = 1. integral = np.zeros_like(padded[..., 0], order='C') - for x_row in range(max(1, - t_row), min(n_row, n_row - t_row)): - for x_col in range(max(1, - t_col), min(n_col, n_col - t_col)): + for x_row in range(max(1, -t_row), min(n_row, n_row - t_row)): + for x_col in range(max(1, -t_col), min(n_col, n_col - t_col)): if n_ch == 1: distance = (padded[x_row, x_col, 0] - padded[x_row + t_row, x_col + t_col, 0])**2 else: distance = ((padded[x_row, x_col, 0] - - padded[x_row + t_row, x_col + t_col, 0])**2 - +(padded[x_row, x_col, 1] - - padded[x_row + t_row, x_col + t_col, 1])**2 - +(padded[x_row, x_col, 2] - - padded[x_row + t_row, x_col + t_col, 2])**2) + padded[x_row + t_row, x_col + t_col, 0])**2 + + (padded[x_row, x_col, 1] - + padded[x_row + t_row, x_col + t_col, 1])**2 + + (padded[x_row, x_col, 2] - + padded[x_row + t_row, x_col + t_col, 2])**2) integral[x_row, x_col] = distance + \ integral[x_row - 1, x_col] + \ - integral[x_row, x_col - 1] \ - - integral[x_row - 1, x_col - 1] + integral[x_row, x_col - 1] - \ + integral[x_row - 1, x_col - 1] for x_row in range(max(offset, offset - t_row), min(n_row - offset, n_row - offset - t_row)): for x_col in range(max(offset, offset - t_col), @@ -456,7 +457,7 @@ def _fast_nl_means_denoising_2d(image, int s=7, int d=13, float h=0.1): weights[x_row + t_row, x_col + t_col] += weight for ch in range(n_ch): result[x_row, x_col, ch] += weight * \ - padded[x_row + t_row, x_col + t_col, ch] + padded[x_row + t_row, x_col + t_col, ch] result[x_row + t_row, x_col + t_col, ch] += \ weight * padded[x_row, x_col, ch] # Normalize pixel values using sum of weights of contributing patches @@ -466,7 +467,7 @@ def _fast_nl_means_denoising_2d(image, int s=7, int d=13, float h=0.1): # No risk of division by zero, since the contribution # of a null shift is strictly positive result[x_row, x_col, channel] /= weights[x_row, x_col] - return result[pad_size: - pad_size, pad_size: - pad_size] + return result[pad_size:-pad_size, pad_size:-pad_size] @cython.cdivision(True) @@ -523,17 +524,17 @@ def _fast_nl_means_denoising_3d(image, int s=5, int d=7, float h=0.1): # Outer loops on patch shifts # With t2 >= 0, reference patch is always on the left of test patch for t_pln in range(-d, d + 1): - x_pln_integral_min = max(1, - t_pln) + x_pln_integral_min = max(1, -t_pln) x_pln_integral_max = min(n_pln, n_pln - t_pln) x_pln_dist_min = max(offset, offset - t_pln) x_pln_dist_max = min(n_pln - offset, n_pln - offset - t_pln) for t_row in range(-d, d + 1): - x_row_integral_min = max(1, - t_row) + x_row_integral_min = max(1, -t_row) x_row_integral_max = min(n_row, n_row - t_row) x_row_dist_min = max(offset, offset - t_row) x_row_dist_max = min(n_row - offset, n_row - offset - t_row) for t_col in range(0, d + 1): - x_col_integral_min = max(1, - t_col) + x_col_integral_min = max(1, -t_col) x_col_integral_max = min(n_col, n_col - t_col) x_col_dist_min = max(offset, offset - t_col) x_col_dist_max = min(n_col - offset, n_col - offset - t_col) @@ -550,15 +551,15 @@ def _fast_nl_means_denoising_3d(image, int s=5, int d=7, float h=0.1): x_col_integral_max): integral[x_pln, x_row, x_col] = \ ((padded[x_pln, x_row, x_col] - - padded[x_pln + t_pln, x_row + t_row, + padded[x_pln + t_pln, x_row + t_row, x_col + t_col])**2 + integral[x_pln - 1, x_row, x_col] + integral[x_pln, x_row - 1, x_col] + integral[x_pln, x_row, x_col - 1] + - integral[x_pln - 1, x_row - 1, x_col - 1] - - integral[x_pln - 1, x_row - 1, x_col] - - integral[x_pln, x_row - 1, x_col - 1] - - integral[x_pln - 1, x_row, x_col - 1]) + integral[x_pln - 1, x_row - 1, x_col - 1] - + integral[x_pln - 1, x_row - 1, x_col] - + integral[x_pln, x_row - 1, x_col - 1] - + integral[x_pln - 1, x_row, x_col - 1]) for x_pln in range(x_pln_dist_min, x_pln_dist_max): for x_row in range(x_row_dist_min, x_row_dist_max): for x_col in range(x_col_dist_min, x_col_dist_max): @@ -581,8 +582,7 @@ def _fast_nl_means_denoising_3d(image, int s=5, int d=7, float h=0.1): for x_pln in range(offset, n_pln - offset): for x_row in range(offset, n_row - offset): for x_col in range(offset, n_col - offset): - # I think there is no risk of division by zero - # except in padded zone + # No risk of division by zero, since the contribution + # of a null shift is strictly positive result[x_pln, x_row, x_col] /= weights[x_pln, x_row, x_col] - return result[pad_size: - pad_size, pad_size: - pad_size, - pad_size: -pad_size] + return result[pad_size:-pad_size, pad_size:-pad_size, pad_size:-pad_size] From ee81e00c0efa156c9ae98109ab002f31bc5d5ef1 Mon Sep 17 00:00:00 2001 From: emmanuelle Date: Sun, 1 Feb 2015 22:30:00 +0100 Subject: [PATCH 24/30] Helper functions for computing the integral of the difference between image and shifted image. --- skimage/restoration/_nl_means_denoising.pyx | 240 +++++++++++++------- 1 file changed, 152 insertions(+), 88 deletions(-) diff --git a/skimage/restoration/_nl_means_denoising.pyx b/skimage/restoration/_nl_means_denoising.pyx index e2db20e9..a941cbbb 100644 --- a/skimage/restoration/_nl_means_denoising.pyx +++ b/skimage/restoration/_nl_means_denoising.pyx @@ -327,7 +327,7 @@ def _nl_means_denoising_3d(image, int s=7, @cython.cdivision(True) @cython.boundscheck(False) cdef inline float _integral_to_distance_2d(DTYPE_t [:, ::] integral, - int x_row, int x_col, int offset, float h2s2): + int row, int col, int offset, float h2s2): """ References ---------- @@ -337,10 +337,10 @@ cdef inline float _integral_to_distance_2d(DTYPE_t [:, ::] integral, Used in _fast_nl_means_denoising_2d """ cdef float distance - distance = integral[x_row + offset, x_col + offset] + \ - integral[x_row - offset, x_col - offset] - \ - integral[x_row - offset, x_col + offset] - \ - integral[x_row + offset, x_col - offset] + 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 @@ -348,7 +348,7 @@ cdef inline float _integral_to_distance_2d(DTYPE_t [:, ::] integral, @cython.cdivision(True) @cython.boundscheck(False) cdef inline float _integral_to_distance_3d(DTYPE_t [:, :, ::] integral, - int x_pln, int x_row, int x_col, int offset, + int pln, int row, int col, int offset, float s_cube_h_square): """ References @@ -359,18 +359,119 @@ cdef inline float _integral_to_distance_3d(DTYPE_t [:, :, ::] integral, Used in _fast_nl_means_denoising_3d """ cdef float distance - distance = (integral[x_pln + offset, x_row + offset, x_col + offset] - - integral[x_pln - offset, x_row - offset, x_col - offset] + - integral[x_pln - offset, x_row - offset, x_col + offset] + - integral[x_pln - offset, x_row + offset, x_col - offset] + - integral[x_pln + offset, x_row - offset, x_col - offset] - - integral[x_pln - offset, x_row + offset, x_col + offset] - - integral[x_pln + offset, x_row - offset, x_col + offset] - - integral[x_pln + offset, x_row + offset, x_col - offset]) + 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(DTYPE_t [:, :, ::] padded, + DTYPE_t [:, ::] 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(DTYPE_t [:, :, ::] padded, + DTYPE_t [:, :, ::] 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): @@ -407,7 +508,7 @@ def _fast_nl_means_denoising_2d(image, int s=7, int d=13, float h=0.1): 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_row, n_col, n_ch, t_row, t_col, x_row, x_col + 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. @@ -427,46 +528,32 @@ def _fast_nl_means_denoising_2d(image, int s=7, int d=13, float h=0.1): else: alpha = 1. integral = np.zeros_like(padded[..., 0], order='C') - for x_row in range(max(1, -t_row), min(n_row, n_row - t_row)): - for x_col in range(max(1, -t_col), min(n_col, n_col - t_col)): - if n_ch == 1: - distance = (padded[x_row, x_col, 0] - - padded[x_row + t_row, x_col + t_col, 0])**2 - else: - distance = ((padded[x_row, x_col, 0] - - padded[x_row + t_row, x_col + t_col, 0])**2 + - (padded[x_row, x_col, 1] - - padded[x_row + t_row, x_col + t_col, 1])**2 + - (padded[x_row, x_col, 2] - - padded[x_row + t_row, x_col + t_col, 2])**2) - integral[x_row, x_col] = distance + \ - integral[x_row - 1, x_col] + \ - integral[x_row, x_col - 1] - \ - integral[x_row - 1, x_col - 1] - for x_row in range(max(offset, offset - t_row), - min(n_row - offset, n_row - offset - t_row)): - for x_col in range(max(offset, offset - t_col), - min(n_col - offset, n_col - offset - t_col)): - distance = _integral_to_distance_2d(integral, x_row, x_col, + _integral_image_2d(padded, integral, t_row, t_col, + n_row, n_col, n_ch) + for row in range(max(offset, offset - t_row), + min(n_row - offset, n_row - offset - t_row)): + for col in range(max(offset, offset - t_col), + min(n_col - offset, n_col - offset - t_col)): + 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) - weights[x_row, x_col] += weight - weights[x_row + t_row, x_col + t_col] += weight + weights[row, col] += weight + weights[row + t_row, col + t_col] += weight for ch in range(n_ch): - result[x_row, x_col, ch] += weight * \ - padded[x_row + t_row, x_col + t_col, ch] - result[x_row + t_row, x_col + t_col, ch] += \ - weight * padded[x_row, x_col, 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 x_row in range(offset, n_row - offset): - for x_col in range(offset, n_col - offset): + 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[x_row, x_col, channel] /= weights[x_row, x_col] + result[row, col, channel] /= weights[row, col] return result[pad_size:-pad_size, pad_size:-pad_size] @@ -506,10 +593,7 @@ def _fast_nl_means_denoising_3d(image, int s=5, int d=7, float h=0.1): cdef DTYPE_t [:, :, ::1] weights = np.zeros_like(padded) cdef DTYPE_t [:, :, ::1] integral = np.zeros_like(padded) cdef int n_pln, n_row, n_col, t_pln, t_row, t_col, \ - x_pln, x_row, x_col - cdef int x_pln_integral_min, x_pln_integral_max, \ - x_row_integral_min, x_row_integral_max, \ - x_col_integral_min, x_col_integral_max + pln, row, col cdef int x_pln_dist_min, x_pln_dist_max, x_row_dist_min, x_row_dist_max, \ x_col_dist_min, x_col_dist_max cdef float weight, distance @@ -524,18 +608,12 @@ def _fast_nl_means_denoising_3d(image, int s=5, int d=7, float h=0.1): # Outer loops on patch shifts # With t2 >= 0, reference patch is always on the left of test patch for t_pln in range(-d, d + 1): - x_pln_integral_min = max(1, -t_pln) - x_pln_integral_max = min(n_pln, n_pln - t_pln) x_pln_dist_min = max(offset, offset - t_pln) x_pln_dist_max = min(n_pln - offset, n_pln - offset - t_pln) for t_row in range(-d, d + 1): - x_row_integral_min = max(1, -t_row) - x_row_integral_max = min(n_row, n_row - t_row) x_row_dist_min = max(offset, offset - t_row) x_row_dist_max = min(n_row - offset, n_row - offset - t_row) for t_col in range(0, d + 1): - x_col_integral_min = max(1, -t_col) - x_col_integral_max = min(n_col, n_col - t_col) x_col_dist_min = max(offset, offset - t_col) x_col_dist_max = min(n_col - offset, n_col - offset - t_col) # alpha is to account for patches on the same column @@ -545,44 +623,30 @@ def _fast_nl_means_denoising_3d(image, int s=5, int d=7, float h=0.1): else: alpha = 1. integral = np.zeros_like(padded) - for x_pln in range(x_pln_integral_min, x_pln_integral_max): - for x_row in range(x_row_integral_min, x_row_integral_max): - for x_col in range(x_col_integral_min, - x_col_integral_max): - integral[x_pln, x_row, x_col] = \ - ((padded[x_pln, x_row, x_col] - - padded[x_pln + t_pln, x_row + t_row, - x_col + t_col])**2 + - integral[x_pln - 1, x_row, x_col] + - integral[x_pln, x_row - 1, x_col] + - integral[x_pln, x_row, x_col - 1] + - integral[x_pln - 1, x_row - 1, x_col - 1] - - integral[x_pln - 1, x_row - 1, x_col] - - integral[x_pln, x_row - 1, x_col - 1] - - integral[x_pln - 1, x_row, x_col - 1]) - for x_pln in range(x_pln_dist_min, x_pln_dist_max): - for x_row in range(x_row_dist_min, x_row_dist_max): - for x_col in range(x_col_dist_min, x_col_dist_max): + _integral_image_3d(padded, integral, t_pln, t_row, t_col, + n_pln, n_row, n_col) + for pln in range(x_pln_dist_min, x_pln_dist_max): + for row in range(x_row_dist_min, x_row_dist_max): + for col in range(x_col_dist_min, x_col_dist_max): distance = _integral_to_distance_3d(integral, - x_pln, x_row, x_col, offset, - s_cube_h_square) + 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) - weights[x_pln, x_row, x_col] += weight - weights[x_pln + t_pln, x_row + t_row, - x_col + t_col] += weight - result[x_pln, x_row, x_col] += weight * \ - padded[x_pln + t_pln, x_row + t_row, - x_col + t_col] - result[x_pln + t_pln, x_row + t_row, - x_col + t_col] += weight * \ - padded[x_pln, x_row, x_col] - for x_pln in range(offset, n_pln - offset): - for x_row in range(offset, n_row - offset): - for x_col in range(offset, n_col - offset): + 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] + 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[x_pln, x_row, x_col] /= weights[x_pln, x_row, x_col] + result[pln, row, col] /= weights[pln, row, col] return result[pad_size:-pad_size, pad_size:-pad_size, pad_size:-pad_size] From 50907cd8e7dc6fb88e611559d89fcc4a245b7100 Mon Sep 17 00:00:00 2001 From: emmanuelle Date: Sun, 1 Feb 2015 22:43:18 +0100 Subject: [PATCH 25/30] Modified variables' names. --- skimage/restoration/_nl_means_denoising.pyx | 147 ++++++++++---------- 1 file changed, 73 insertions(+), 74 deletions(-) diff --git a/skimage/restoration/_nl_means_denoising.pyx b/skimage/restoration/_nl_means_denoising.pyx index a941cbbb..82018bff 100644 --- a/skimage/restoration/_nl_means_denoising.pyx +++ b/skimage/restoration/_nl_means_denoising.pyx @@ -174,9 +174,9 @@ def _nl_means_denoising_2d(image, int s=7, int d=13, float h=0.1): cdef int n_row, n_col, n_ch n_row, n_col, n_ch = image.shape cdef int offset = s / 2 - cdef int x_row, x_col, i, j, color - cdef int x_row_start, x_row_end, x_col_start, x_col_end - cdef int x_row_start_i, x_row_end_i, x_col_start_j, x_col_end_j + 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 DTYPE_t [::1] new_values = np.zeros(n_ch).astype(np.float32) cdef DTYPE_t [:, :, ::1] padded = np.ascontiguousarray(util.pad(image, ((offset, offset), (offset, offset), (0, 0)), @@ -192,44 +192,44 @@ def _nl_means_denoising_2d(image, int s=7, int d=13, float h=0.1): cdef float distance w = 1. / (n_ch * np.sum(w) * h ** 2) * w # Coordinates of central pixel and patch bounds - for x_row in range(offset, n_row + offset): - x_row_start = x_row - offset - x_row_end = x_row + offset + 1 - for x_col in range(offset, n_col + offset): + for row in range(offset, n_row + offset): + row_start = row - offset + row_end = row + offset + 1 + for col in range(offset, n_col + offset): for color in range(n_ch): new_values[color] = 0 weight_sum = 0 - x_col_start = x_col - offset - x_col_end = x_col + offset + 1 + col_start = col - offset + col_end = col + offset + 1 # Coordinates of test pixel and patch bounds - for i in range(max(-d, offset - x_row), - min(d + 1, n_row + offset - x_row)): - x_row_start_i = x_row_start + i - x_row_end_i = x_row_end + i - for j in range(max(-d, offset - x_col), - min(d + 1, n_col + offset - x_col)): - x_col_start_j = x_col_start + j - x_col_end_j = x_col_end + j + 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 + 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 if n_ch == 1: weight = patch_distance_2d( - padded[x_row_start:x_row_end, - x_col_start:x_col_end, 0], - padded[x_row_start_i:x_row_end_i, - x_col_start_j:x_col_end_j, 0], + 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[x_row_start:x_row_end, - x_col_start:x_col_end, :], - padded[x_row_start_i:x_row_end_i, - x_col_start_j:x_col_end_j, :], + padded[row_start:row_end, + col_start:col_end, :], + padded[row_start_i:row_end_i, + col_start_j:col_end_j, :], w, s) weight_sum += weight for color in range(n_ch): - new_values[color] += weight * padded[x_row + i, - x_col + j, color] + new_values[color] += weight * padded[row + i, + col + j, color] for color in range(n_ch): - result[x_row, x_col, color] = new_values[color] / weight_sum + result[row, col, color] = new_values[color] / weight_sum return result[offset:-offset, offset:-offset] @@ -276,49 +276,48 @@ def _nl_means_denoising_3d(image, int s=7, -(xg_pln ** 2 + xg_row ** 2 + xg_col ** 2) / (2 * A ** 2)).astype(np.float32)) cdef float distance - cdef int x_pln, x_row, x_col, i, j, k - cdef int x_pln_start, x_pln_end, x_row_start, x_row_end, \ - x_col_start, x_col_end - cdef int x_pln_start_i, x_pln_end_i, x_row_start_j, x_row_end_j, \ - x_col_start_k, x_col_end_k + 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 and patch bounds - for x_pln in range(offset, n_pln + offset): - x_pln_start = x_pln - offset - x_pln_end = x_pln + offset + 1 - for x_row in range(offset, n_row + offset): - x_row_start = x_row - offset - x_row_end = x_row + offset + 1 - for x_col in range(offset, n_col + offset): - x_col_start = x_col - offset - x_col_end = x_col + offset + 1 + for pln in range(offset, n_pln + offset): + pln_start = pln - offset + pln_end = pln + offset + 1 + for row in range(offset, n_row + offset): + row_start = row - offset + row_end = row + offset + 1 + for col in range(offset, n_col + offset): + col_start = col - offset + col_end = col + offset + 1 new_value = 0 weight_sum = 0 # Coordinates of test pixel and patch bounds - for i in range(max(-d, offset - x_pln), - min(d + 1, n_pln + offset - x_pln)): - x_pln_start_i = x_pln_start + i - x_pln_end_i = x_pln_end + i - for j in range(max(-d, offset - x_row), - min(d + 1, n_row + offset - x_row)): - x_row_start_j = x_row_start + j - x_row_end_j = x_row_end + j - for k in range(max(-d, offset - x_col), - min(d + 1, n_col + offset - x_col)): - x_col_start_k = x_col_start + k - x_col_end_k = x_col_end + k + 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 + 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 + 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[x_pln_start:x_pln_end, - x_row_start:x_row_end, - x_col_start:x_col_end], - padded[x_pln_start_i:x_pln_end_i, - x_row_start_j:x_row_end_j, - x_col_start_k:x_col_end_k], + 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) weight_sum += weight - new_value += weight * padded[x_pln + i, - x_row + j, x_col + k] - result[x_pln, x_row, x_col] = new_value / weight_sum + new_value += weight * padded[pln + i, + row + j, col + k] + result[pln, row, col] = new_value / weight_sum return result[offset:-offset, offset:-offset, offset:-offset] #-------------- Accelerated algorithm of Froment 2015 ------------------ @@ -594,8 +593,8 @@ def _fast_nl_means_denoising_3d(image, int s=5, int d=7, float h=0.1): cdef DTYPE_t [:, :, ::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 x_pln_dist_min, x_pln_dist_max, x_row_dist_min, x_row_dist_max, \ - x_col_dist_min, x_col_dist_max + 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. @@ -608,14 +607,14 @@ def _fast_nl_means_denoising_3d(image, int s=5, int d=7, float h=0.1): # Outer loops on patch shifts # With t2 >= 0, reference patch is always on the left of test patch for t_pln in range(-d, d + 1): - x_pln_dist_min = max(offset, offset - t_pln) - x_pln_dist_max = min(n_pln - offset, n_pln - offset - t_pln) + pln_dist_min = max(offset, offset - t_pln) + pln_dist_max = min(n_pln - offset, n_pln - offset - t_pln) for t_row in range(-d, d + 1): - x_row_dist_min = max(offset, offset - t_row) - x_row_dist_max = min(n_row - offset, n_row - offset - t_row) + row_dist_min = max(offset, offset - t_row) + row_dist_max = min(n_row - offset, n_row - offset - t_row) for t_col in range(0, d + 1): - x_col_dist_min = max(offset, offset - t_col) - x_col_dist_max = min(n_col - offset, n_col - offset - t_col) + 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): @@ -625,9 +624,9 @@ def _fast_nl_means_denoising_3d(image, int s=5, int d=7, float h=0.1): integral = np.zeros_like(padded) _integral_image_3d(padded, integral, t_pln, t_row, t_col, n_pln, n_row, n_col) - for pln in range(x_pln_dist_min, x_pln_dist_max): - for row in range(x_row_dist_min, x_row_dist_max): - for col in range(x_col_dist_min, x_col_dist_max): + for pln in range(pln_dist_min, pln_dist_max): + for row in range(row_dist_min, row_dist_max): + for col in range(col_dist_min, col_dist_max): distance = _integral_to_distance_3d(integral, pln, row, col, offset, s_cube_h_square) # exp of large negative numbers is close to zero From 6770e9e5574388c4ef031c70e818c533aea6ef67 Mon Sep 17 00:00:00 2001 From: emmanuelle Date: Sun, 1 Feb 2015 22:51:00 +0100 Subject: [PATCH 26/30] Modified name of image dtype. --- skimage/restoration/_nl_means_denoising.pyx | 64 ++++++++++----------- 1 file changed, 32 insertions(+), 32 deletions(-) diff --git a/skimage/restoration/_nl_means_denoising.pyx b/skimage/restoration/_nl_means_denoising.pyx index 82018bff..8f6692d2 100644 --- a/skimage/restoration/_nl_means_denoising.pyx +++ b/skimage/restoration/_nl_means_denoising.pyx @@ -4,14 +4,14 @@ cimport numpy as np cimport cython from libc.math cimport exp -ctypedef np.float32_t DTYPE_t +ctypedef np.float32_t IMGDTYPE cdef float DISTANCE_CUTOFF = 5. @cython.boundscheck(False) -cdef inline float patch_distance_2d(DTYPE_t [:, :] p1, - DTYPE_t [:, :] p2, - DTYPE_t [:, ::] w, int s): +cdef inline float patch_distance_2d(IMGDTYPE [:, :] p1, + IMGDTYPE [:, :] p2, + IMGDTYPE [:, ::] w, int s): """ Compute a Gaussian distance between two image patches. @@ -57,9 +57,9 @@ cdef inline float patch_distance_2d(DTYPE_t [:, :] p1, @cython.boundscheck(False) -cdef inline float patch_distance_2drgb(DTYPE_t [:, :, :] p1, - DTYPE_t [:, :, :] p2, - DTYPE_t [:, ::] w, int s): +cdef inline float patch_distance_2drgb(IMGDTYPE [:, :, :] p1, + IMGDTYPE [:, :, :] p2, + IMGDTYPE [:, ::] w, int s): """ Compute a Gaussian distance between two image patches. @@ -103,9 +103,9 @@ cdef inline float patch_distance_2drgb(DTYPE_t [:, :, :] p1, @cython.boundscheck(False) -cdef inline float patch_distance_3d(DTYPE_t [:, :, :] p1, - DTYPE_t [:, :, :] p2, - DTYPE_t [:, :, ::] w, int s): +cdef inline float patch_distance_3d(IMGDTYPE [:, :, :] p1, + IMGDTYPE [:, :, :] p2, + IMGDTYPE [:, :, ::] w, int s): """ Compute a Gaussian distance between two image patches. @@ -177,16 +177,16 @@ def _nl_means_denoising_2d(image, int s=7, int d=13, float h=0.1): 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 DTYPE_t [::1] new_values = np.zeros(n_ch).astype(np.float32) - cdef DTYPE_t [:, :, ::1] padded = np.ascontiguousarray(util.pad(image, + 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 DTYPE_t [:, :, ::1] result = padded.copy() + 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 DTYPE_t [:, ::1] w = np.ascontiguousarray(np.exp( + cdef IMGDTYPE [:, ::1] w = np.ascontiguousarray(np.exp( -(xg_row ** 2 + xg_col ** 2) / (2 * A ** 2)). astype(np.float32)) cdef float distance @@ -262,19 +262,19 @@ def _nl_means_denoising_3d(image, int s=7, n_pln, n_row, n_col = 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( + cdef IMGDTYPE [:, :, ::1] padded = np.ascontiguousarray(util.pad( image.astype(np.float32), offset, mode='reflect')) - cdef DTYPE_t [:, :, ::1] result = padded.copy() + 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 DTYPE_t [:, :, ::1] w = np.ascontiguousarray(np.exp( + cdef IMGDTYPE [:, :, ::1] w = np.ascontiguousarray(np.exp( -(xg_pln ** 2 + xg_row ** 2 + xg_col ** 2) / - (2 * A ** 2)).astype(np.float32)) + (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 @@ -325,7 +325,7 @@ def _nl_means_denoising_3d(image, int s=7, @cython.cdivision(True) @cython.boundscheck(False) -cdef inline float _integral_to_distance_2d(DTYPE_t [:, ::] integral, +cdef inline float _integral_to_distance_2d(IMGDTYPE [:, ::] integral, int row, int col, int offset, float h2s2): """ References @@ -346,7 +346,7 @@ cdef inline float _integral_to_distance_2d(DTYPE_t [:, ::] integral, @cython.cdivision(True) @cython.boundscheck(False) -cdef inline float _integral_to_distance_3d(DTYPE_t [:, :, ::] integral, +cdef inline float _integral_to_distance_3d(IMGDTYPE [:, :, ::] integral, int pln, int row, int col, int offset, float s_cube_h_square): """ @@ -372,8 +372,8 @@ cdef inline float _integral_to_distance_3d(DTYPE_t [:, :, ::] integral, @cython.cdivision(True) @cython.boundscheck(False) -cdef inline _integral_image_2d(DTYPE_t [:, :, ::] padded, - DTYPE_t [:, ::] integral, int t_row, +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`` @@ -422,8 +422,8 @@ cdef inline _integral_image_2d(DTYPE_t [:, :, ::] padded, @cython.cdivision(True) @cython.boundscheck(False) -cdef inline _integral_image_3d(DTYPE_t [:, :, ::] padded, - DTYPE_t [:, :, ::] integral, int t_pln, +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): """ @@ -501,12 +501,12 @@ def _fast_nl_means_denoising_2d(image, int s=7, int d=13, float h=0.1): # 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, + 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 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 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 @@ -586,11 +586,11 @@ def _fast_nl_means_denoising_3d(image, int s=5, int d=7, float h=0.1): # 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, + cdef IMGDTYPE [:, :, ::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 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, \ From bcd24eab862d76a3bf29fad22f1df3b65a688fb5 Mon Sep 17 00:00:00 2001 From: emmanuelle Date: Sun, 1 Feb 2015 23:11:34 +0100 Subject: [PATCH 27/30] Removed a few redundant lines. --- skimage/restoration/non_local_means.py | 23 +++++++++-------------- 1 file changed, 9 insertions(+), 14 deletions(-) diff --git a/skimage/restoration/non_local_means.py b/skimage/restoration/non_local_means.py index a94437c5..2ba98a11 100644 --- a/skimage/restoration/non_local_means.py +++ b/skimage/restoration/non_local_means.py @@ -59,7 +59,7 @@ def nl_means_denoising(image, patch_size=7, patch_distance=11, h=0.1, 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=True``, for which + 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 @@ -99,27 +99,22 @@ def nl_means_denoising(image, patch_size=7, patch_distance=11, h=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, - s=patch_size, d=patch_distance, h=h))) + patch_size, patch_distance, h))) else: return np.squeeze(np.array(_nl_means_denoising_2d(image, - s=patch_size, d=patch_distance, h=h))) - elif image.ndim == 3 and not multichannel: # only grayscale + 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)) - if image.ndim == 3 and multichannel: # 2-D color (RGB) images - if fast_mode: - return np.array(_fast_nl_means_denoising_2d(image, patch_size, - patch_distance, h)) - else: - return np.array(_nl_means_denoising_2d(image, patch_size, - patch_distance, h)) - else: - raise NotImplementedError("Non-local means denoising is only \ - implemented for 2D grayscale and RGB images or 3-D grayscale images.") From 0f4b09257843d0662736671564e806055bc1601a Mon Sep 17 00:00:00 2001 From: emmanuelle Date: Sun, 1 Feb 2015 23:31:09 +0100 Subject: [PATCH 28/30] Added comment for inner loops. --- skimage/restoration/_nl_means_denoising.pyx | 2 ++ 1 file changed, 2 insertions(+) diff --git a/skimage/restoration/_nl_means_denoising.pyx b/skimage/restoration/_nl_means_denoising.pyx index 8f6692d2..1461f485 100644 --- a/skimage/restoration/_nl_means_denoising.pyx +++ b/skimage/restoration/_nl_means_denoising.pyx @@ -529,6 +529,7 @@ def _fast_nl_means_denoising_2d(image, int s=7, int d=13, float h=0.1): 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 for row in range(max(offset, offset - t_row), min(n_row - offset, n_row - offset - t_row)): for col in range(max(offset, offset - t_col), @@ -624,6 +625,7 @@ def _fast_nl_means_denoising_3d(image, int s=5, int d=7, float h=0.1): 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 for pln in range(pln_dist_min, pln_dist_max): for row in range(row_dist_min, row_dist_max): for col in range(col_dist_min, col_dist_max): From 4c9dbe1e18b5fdbf5399184b0034f09a4747067e Mon Sep 17 00:00:00 2001 From: emmanuelle Date: Mon, 2 Feb 2015 22:55:49 +0100 Subject: [PATCH 29/30] Added a comment about normalization. --- skimage/restoration/_nl_means_denoising.pyx | 1 + 1 file changed, 1 insertion(+) diff --git a/skimage/restoration/_nl_means_denoising.pyx b/skimage/restoration/_nl_means_denoising.pyx index 1461f485..93ab51c5 100644 --- a/skimage/restoration/_nl_means_denoising.pyx +++ b/skimage/restoration/_nl_means_denoising.pyx @@ -644,6 +644,7 @@ def _fast_nl_means_denoising_3d(image, int s=5, int d=7, float h=0.1): 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): From 3eb775542f414016f51fbb278312fb8b79518e27 Mon Sep 17 00:00:00 2001 From: emmanuelle Date: Tue, 10 Feb 2015 23:22:42 +0100 Subject: [PATCH 30/30] Added some comments inside Cython functions. --- skimage/restoration/_nl_means_denoising.pyx | 66 +++++++++++++++++++-- 1 file changed, 62 insertions(+), 4 deletions(-) diff --git a/skimage/restoration/_nl_means_denoising.pyx b/skimage/restoration/_nl_means_denoising.pyx index 93ab51c5..089914fc 100644 --- a/skimage/restoration/_nl_means_denoising.pyx +++ b/skimage/restoration/_nl_means_denoising.pyx @@ -191,25 +191,34 @@ def _nl_means_denoising_2d(image, int s=7, int d=13, float h=0.1): astype(np.float32)) cdef float distance w = 1. / (n_ch * np.sum(w) * h ** 2) * w - # Coordinates of central pixel and patch bounds + + # 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 - # Coordinates of test pixel and patch bounds + + # 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, @@ -224,12 +233,19 @@ def _nl_means_denoising_2d(image, int s=7, int d=13, float h=0.1): 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] @@ -281,27 +297,35 @@ def _nl_means_denoising_3d(image, int s=7, 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 and patch bounds + + # 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 - # Coordinates of test pixel and patch bounds + + # 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 @@ -314,10 +338,15 @@ def _nl_means_denoising_3d(image, int s=7, 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 ------------------ @@ -516,9 +545,12 @@ def _fast_nl_means_denoising_2d(image, int s=7, int d=13, float h=0.1): 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 @@ -526,27 +558,36 @@ def _fast_nl_means_denoising_2d(image, int s=7, int d=13, float h=0.1): 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): @@ -554,6 +595,8 @@ def _fast_nl_means_denoising_2d(image, int s=7, int d=13, float h=0.1): # 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] @@ -605,14 +648,18 @@ def _fast_nl_means_denoising_3d(image, int s=5, int d=7, float h=0.1): 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) @@ -622,19 +669,27 @@ def _fast_nl_means_denoising_3d(image, int s=5, int d=7, float h=0.1): 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 @@ -644,6 +699,7 @@ def _fast_nl_means_denoising_3d(image, int s=5, int d=7, float h=0.1): 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): @@ -651,4 +707,6 @@ def _fast_nl_means_denoising_3d(image, int s=5, int d=7, float h=0.1): # 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]