diff --git a/doc/examples/plot_denoise.py b/doc/examples/plot_denoise.py index 12e8b9b6..a3e83a65 100644 --- a/doc/examples/plot_denoise.py +++ b/doc/examples/plot_denoise.py @@ -46,14 +46,14 @@ ax[0, 0].set_title('noisy') ax[0, 1].imshow(tv_denoise(noisy, weight=0.1)) ax[0, 1].axis('off') ax[0, 1].set_title('TV') -ax[0, 2].imshow(denoise_bilateral(noisy, sigma_color=0.03, sigma_range=15)) +ax[0, 2].imshow(denoise_bilateral(noisy, sigma_range=0.03, sigma_spatial=15)) ax[0, 2].axis('off') ax[0, 2].set_title('Bilateral') ax[1, 0].imshow(tv_denoise(noisy, weight=0.2)) ax[1, 0].axis('off') ax[1, 0].set_title('(more) TV') -ax[1, 1].imshow(denoise_bilateral(noisy, sigma_color=0.06, sigma_range=15)) +ax[1, 1].imshow(denoise_bilateral(noisy, sigma_range=0.06, sigma_spatial=15)) ax[1, 1].axis('off') ax[1, 1].set_title('(more) Bilateral') ax[1, 2].imshow(lena) diff --git a/skimage/filter/__init__.py b/skimage/filter/__init__.py index 02fddf23..f1c1fd49 100644 --- a/skimage/filter/__init__.py +++ b/skimage/filter/__init__.py @@ -3,6 +3,7 @@ from .ctmf import median_filter from ._canny import canny from .edges import (sobel, hsobel, vsobel, scharr, hscharr, vscharr, prewitt, hprewitt, vprewitt) -from .denoise import tv_denoise, denoise_tv, denoise_bilateral +from .denoise import tv_denoise, denoise_tv +from ._denoise import denoise_bilateral from ._rank_order import rank_order from .thresholding import threshold_otsu, threshold_adaptive diff --git a/skimage/filter/_denoise.pyx b/skimage/filter/_denoise.pyx index ad25d9d0..bda46b0d 100644 --- a/skimage/filter/_denoise.pyx +++ b/skimage/filter/_denoise.pyx @@ -7,7 +7,8 @@ cimport numpy as cnp import numpy as np from libc.math cimport exp, fabs, sqrt from libc.stdlib cimport malloc, free -from skimage._shared.interpolation cimport get_pixel2d, get_pixel3d +from skimage._shared.interpolation cimport get_pixel3d +from skimage.util import img_as_float cdef inline double _gaussian_weight(double sigma, double value): @@ -42,67 +43,57 @@ cdef double* _compute_range_lut(int win_size, double sigma): return range_lut -def _denoise_bilateral2d(cnp.ndarray[dtype=cnp.double_t, ndim=2, mode='c'] image, - int win_size, double sigma_color, - double sigma_range, int bins, char mode, - double cval): +def denoise_bilateral(image, int win_size=5, sigma_range=None, + double sigma_spatial=1, int bins=10000, mode='constant', + double cval=0): + """Denoise image using bilateral filter. - cdef: - Py_ssize_t rows = image.shape[0] - Py_ssize_t cols = image.shape[1] - Py_ssize_t window_ext = (win_size - 1) / 2 - double max_value = image.max() + This is an edge-preserving and noise reducing denoising filter. It averages + pixels based on their spatial closeness and radiometric similarity. - cnp.ndarray[dtype=cnp.double_t, ndim=2, mode='c'] out = \ - np.zeros((rows, cols), dtype=np.double) + Spatial closeness is measured by the gaussian function of the euclidian + distance between two pixels and a certain standard deviation + (`sigma_spatial`). - double* image_data = image.data - double* out_data = out.data + Radiometric similarity is measured by the gaussian function of the euclidian + distance between two color values and a certain standard deviation + (`sigma_range`). - double* color_lut = _compute_color_lut(bins, sigma_color, max_value) - double* range_lut = _compute_range_lut(win_size, sigma_range) + Parameters + ---------- + image : ndarray + Input image. + win_size : int + Window size for filtering. + sigma_range : float + Standard deviation for grayvalue/color distance (radiometric + similarity). A larger value results in averaging of pixels with larger + radiometric differences. + sigma_spatial : float + Standard deviation for range distance. A larger value results in + averaging of pixels with larger spatial differences. + bins : int + Number of discrete values for gaussian weights of color filtering. + A larger value results in improved accuracy. + mode : string + How to handle values outside the image borders. See + `scipy.ndimage.map_coordinates` for detail. + cval : string + Used in conjunction with mode 'constant', the value outside + the image boundaries. - Py_ssize_t r, c, wr, wc, kr, kc, rr, cc, pixel_addr - double centre, value, weight, total_value, total_weight, \ - color_weight, range_weight, diff - double dist_scale = bins / max_value + Returns + ------- + denoised : ndarray + Denoised image. - for r in range(rows): - for c in range(cols): - pixel_addr = r * cols + c - total_value = 0 - total_weight = 0 - centre = image_data[pixel_addr] - for wr in range(-window_ext, window_ext + 1): - rr = wr + r - kr = wr + window_ext - for wc in range(-window_ext, window_ext + 1): - cc = wc + c - kc = wc + window_ext + References + ---------- + .. [1] http://users.soe.ucsc.edu/~manduchi/Papers/ICCV98.pdf - value = get_pixel2d(image_data, rows, cols, - rr, cc, mode, cval) - diff = fabs(centre - value) + """ - range_weight = range_lut[kr * win_size + kc] - color_weight = color_lut[(diff * dist_scale)] - - weight = range_weight * color_weight - total_value += value * weight - total_weight += weight - - out_data[pixel_addr] = total_value / total_weight - - free(color_lut) - free(range_lut) - - return out - - -def _denoise_bilateral3d(cnp.ndarray[dtype=cnp.double_t, ndim=3, mode='c'] image, - int win_size, double sigma_color, - double sigma_range, int bins, char mode, - double cval): + image = np.atleast_3d(img_as_float(image)) cdef: Py_ssize_t rows = image.shape[0] @@ -112,22 +103,35 @@ def _denoise_bilateral3d(cnp.ndarray[dtype=cnp.double_t, ndim=3, mode='c'] image double max_value = image.max() + cnp.ndarray[dtype=cnp.double_t, ndim=3, mode='c'] cimage = \ + np.ascontiguousarray(image) cnp.ndarray[dtype=cnp.double_t, ndim=3, mode='c'] out = \ np.zeros((rows, cols, dims), dtype=np.double) - double* image_data = image.data + double* image_data = cimage.data double* out_data = out.data - double* color_lut = _compute_color_lut(bins, sigma_color, max_value) - double* range_lut = _compute_range_lut(win_size, sigma_range) + double* color_lut = _compute_color_lut(bins, sigma_range, max_value) + double* range_lut = _compute_range_lut(win_size, sigma_spatial) Py_ssize_t r, c, d, wr, wc, kr, kc, rr, cc, pixel_addr - double value, weight, dist, total_weight, color_weight, range_weight + double value, weight, dist, total_weight, csigma_range, color_weight, \ + range_weight double dist_scale = bins / dims / max_value double* values = malloc(dims * sizeof(double)) double* centres = malloc(dims * sizeof(double)) double* total_values = malloc(dims * sizeof(double)) + if sigma_range is None: + csigma_range = image.std() + else: + csigma_range = sigma_range + + if mode not in ('constant', 'wrap', 'reflect', 'nearest'): + raise ValueError("Invalid mode specified. Please use " + "`constant`, `nearest`, `wrap` or `reflect`.") + cdef char cmode = ord(mode[0].upper()) + for r in range(rows): for c in range(cols): pixel_addr = r * cols * dims + c * dims @@ -147,7 +151,7 @@ def _denoise_bilateral3d(cnp.ndarray[dtype=cnp.double_t, ndim=3, mode='c'] image dist = 0 for d in range(dims): value = get_pixel3d(image_data, rows, cols, dims, - rr, cc, d, mode, cval) + rr, cc, d, cmode, cval) values[d] = value dist += (centres[d] - value)**2 dist = sqrt(dist) diff --git a/skimage/filter/denoise.py b/skimage/filter/denoise.py index d7c96f90..d62f8b27 100644 --- a/skimage/filter/denoise.py +++ b/skimage/filter/denoise.py @@ -244,70 +244,3 @@ def denoise_tv(im, weight=50, eps=2.e-4, n_iter_max=200): tv_denoise = deprecated('skimage.filter.denoise_tv')(denoise_tv) - - -def denoise_bilateral(image, win_size=5, sigma_color=1, sigma_range=1, bins=1e4, - mode='constant', cval=0): - """Denoise image using bilateral filter. - - This is an edge-preserving and noise reducing denoising filter. It averages - pixels based on their spatial closeness and radiometric similarity. - - Spatial closeness is measured by the gaussian function of the euclidian - distance between two pixels and a certain standard deviation - (`sigma_range`). - - Radiometric similarity is measured by the gaussian function of the euclidian - distance between two color values and a certain standard deviation - (`sigma_color`). - - Parameters - ---------- - image : ndarray - Input image. - win_size : int - Window size for filtering. - sigma_color : float - Standard deviation for color distance. A larger value results in - averaging of pixels with larger color differences. - sigma_range : float - Standard deviation for range distance. A larger value results in - averaging of pixels with larger spatial differences. - bins : int - Number of discrete values for gaussian weights of color filtering. - A larger value results in improved accuracy. - mode : string - How to handle values outside the image borders. See - `scipy.ndimage.map_coordinates` for detail. - cval : string - Used in conjunction with mode 'constant', the value outside - the image boundaries. - - Returns - ------- - denoised : ndarray - Denoised image. - - References - ---------- - .. [1] http://users.soe.ucsc.edu/~manduchi/Papers/ICCV98.pdf - - """ - - # not using img_as_float to preserve original range of values, which is - # necessary so sigma_color is applied as user desires - image = np.array(image, dtype=np.double) - - if mode not in ('constant', 'wrap', 'reflect', 'nearest'): - raise ValueError("Invalid mode specified. Please use " - "`constant`, `nearest`, `wrap` or `reflect`.") - mode = ord(mode[0].upper()) - - if image.ndim == 2 or (image.ndim == 3 and image.shape[2] == 1): - if image.ndim == 3 and image.shape[2] == 1: - image = np.squeeze(image) - func = _denoise._denoise_bilateral2d - else: - func = _denoise._denoise_bilateral3d - image = np.ascontiguousarray(image) - return func(image, win_size, sigma_color, sigma_range, bins, mode, cval) diff --git a/skimage/filter/tests/test_denoise.py b/skimage/filter/tests/test_denoise.py index b8e6c9ef..63a5a5e0 100644 --- a/skimage/filter/tests/test_denoise.py +++ b/skimage/filter/tests/test_denoise.py @@ -64,8 +64,8 @@ def test_denoise_bilateral_2d(): img += 0.5 * img.std() * np.random.random(img.shape) img = np.clip(img, 0, 1) - out1 = filter.denoise_bilateral(img, sigma_color=0.1, sigma_range=20) - out2 = filter.denoise_bilateral(img, sigma_color=0.2, sigma_range=30) + 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() @@ -78,8 +78,8 @@ def test_denoise_bilateral_3d(): img += 0.5 * img.std() * np.random.random(img.shape) img = np.clip(img, 0, 1) - out1 = filter.denoise_bilateral(img, sigma_color=0.1, sigma_range=20) - out2 = filter.denoise_bilateral(img, sigma_color=0.2, sigma_range=30) + 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()