From f9ecd140d8942f1db6226df4e7a15d56d79d922c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Sat, 6 Oct 2012 00:34:05 +0200 Subject: [PATCH] Add bilateral denoising filter --- skimage/filter/__init__.py | 2 +- skimage/filter/_denoise.pyx | 160 ++++++++++++++++++++++++++++++++++++ skimage/filter/denoise.py | 57 +++++++++++++ skimage/filter/setup.py | 3 + 4 files changed, 221 insertions(+), 1 deletion(-) create mode 100644 skimage/filter/_denoise.pyx diff --git a/skimage/filter/__init__.py b/skimage/filter/__init__.py index dc5f9896..eec88038 100644 --- a/skimage/filter/__init__.py +++ b/skimage/filter/__init__.py @@ -3,6 +3,6 @@ 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 +from .denoise import tv_denoise, 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 new file mode 100644 index 00000000..bf6c372a --- /dev/null +++ b/skimage/filter/_denoise.pyx @@ -0,0 +1,160 @@ +#cython: cdivision=True +#cython: boundscheck=False +#cython: nonecheck=False +#cython: wraparound=False + +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 + + +cdef inline double _gaussian_weight(double sigma, double value): + return exp(- 0.5 * (value / sigma)**2) + + +cdef double* _compute_color_lut(int bins, double sigma, double max_value): + cdef double* color_lut = malloc(bins * sizeof(double)) + cdef Py_ssize_t b + for b in range(bins): + color_lut[b] = _gaussian_weight(sigma, b * max_value / bins) + return color_lut + + +cdef double* _compute_range_lut(int win_size, double sigma): + cdef double* range_lut = malloc(win_size**2 * sizeof(double)) + cdef Py_ssize_t kr, kc + cdef Py_ssize_t window_ext = (win_size - 1) / 2 + cdef double dist + for kr in range(win_size): + for kc in range(win_size): + dist = sqrt((kr - window_ext)**2 + (kc - window_ext)**2) + range_lut[kr * win_size + kc] = _gaussian_weight(sigma, dist) + 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): + + cdef Py_ssize_t rows = image.shape[0] + cdef Py_ssize_t cols = image.shape[1] + cdef Py_ssize_t window_ext = (win_size - 1) / 2 + cdef double max_value = image.max() + + cdef cnp.ndarray[dtype=cnp.double_t, ndim=2, mode='c'] out = \ + np.zeros((rows, cols), dtype=np.double) + + cdef double* image_data = image.data + cdef double* out_data = out.data + + cdef double* color_lut = _compute_color_lut(bins, sigma_color, max_value) + cdef double* range_lut = _compute_range_lut(win_size, sigma_range) + + cdef Py_ssize_t r, c, wr, wc, kr, kc, rr, cc, pixel_addr + cdef double centre, value, weight, total_value, total_weight, \ + color_weight, range_weight, diff + cdef double dist_scale = bins / max_value + + 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 + + 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): + + cdef Py_ssize_t rows = image.shape[0] + cdef Py_ssize_t cols = image.shape[1] + cdef Py_ssize_t dims = image.shape[2] + cdef Py_ssize_t window_ext = (win_size - 1) / 2 + cdef double max_value = image.max() + + cdef cnp.ndarray[dtype=cnp.double_t, ndim=3, mode='c'] out = \ + np.zeros((rows, cols, dims), dtype=np.double) + + cdef double* image_data = image.data + cdef double* out_data = out.data + + cdef double* color_lut = _compute_color_lut(bins, sigma_color, max_value) + cdef double* range_lut = _compute_range_lut(win_size, sigma_range) + + cdef Py_ssize_t r, c, d, wr, wc, kr, kc, rr, cc, pixel_addr + cdef double value, weight, dist, total_weight, color_weight, range_weight + cdef double dist_scale = bins / dims / max_value + cdef double* values = malloc(dims * sizeof(double)) + cdef double* centres = malloc(dims * sizeof(double)) + cdef double* total_values = malloc(dims * sizeof(double)) + + for r in range(rows): + for c in range(cols): + pixel_addr = r * cols * dims + c * dims + total_weight = 0 + for d in range(dims): + total_values[d] = 0 + centres[d] = image_data[pixel_addr + d] + 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 + + # save pixel values for all dims and compute euclidian + # distance between centre stack and current position + dist = 0 + for d in range(dims): + value = get_pixel3d(image_data, rows, cols, dims, + rr, cc, d, mode, cval) + values[d] = value + dist += (centres[d] - value)**2 + dist = sqrt(dist) + + range_weight = range_lut[kr * win_size + kc] + color_weight = color_lut[(dist * dist_scale)] + + weight = range_weight * color_weight + for d in range(dims): + total_values[d] += values[d] * weight + total_weight += weight + for d in range(dims): + out_data[pixel_addr + d] = total_values[d] / total_weight + + free(color_lut) + free(range_lut) + free(values) + free(centres) + free(total_values) + + return out diff --git a/skimage/filter/denoise.py b/skimage/filter/denoise.py index ae4d2f82..f6b30f6f 100644 --- a/skimage/filter/denoise.py +++ b/skimage/filter/denoise.py @@ -1,5 +1,6 @@ import numpy as np from skimage import img_as_float +import _denoise def _tv_denoise_3d(im, weight=100, eps=2.e-4, n_iter_max=200): @@ -239,3 +240,59 @@ def tv_denoise(im, weight=50, eps=2.e-4, n_iter_max=200): raise ValueError('only 2-d and 3-d images may be denoised with this ' 'function') return out + + +def denoise_bilateral(image, win_size=5, sigma_color=1, sigma_range=1, bins=1e4, + mode='constant', cval=0): + """Denoise image using bilateral filter. + + 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/setup.py b/skimage/filter/setup.py index 3e573aec..b996055b 100644 --- a/skimage/filter/setup.py +++ b/skimage/filter/setup.py @@ -13,9 +13,12 @@ def configuration(parent_package='', top_path=None): config.add_data_dir('tests') cython(['_ctmf.pyx'], working_path=base_path) + cython(['_denoise.pyx'], working_path=base_path) config.add_extension('_ctmf', sources=['_ctmf.c'], include_dirs=[get_numpy_include_dirs()]) + config.add_extension('_denoise', sources=['_denoise.c'], + include_dirs=[get_numpy_include_dirs(), '../_shared']) return config