Combine 2D and 3D bilateral filter in one function and rename sigma parameters

This commit is contained in:
Johannes Schönberger
2012-10-17 16:28:44 +02:00
parent 5b8f554a2e
commit c04f14c129
5 changed files with 71 additions and 133 deletions
+2 -2
View File
@@ -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)
+2 -1
View File
@@ -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
+63 -59
View File
@@ -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 = <double*>image.data
double* out_data = <double*>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[<int>(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 = <double*>image.data
double* image_data = <double*>cimage.data
double* out_data = <double*>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 = <double*>malloc(dims * sizeof(double))
double* centres = <double*>malloc(dims * sizeof(double))
double* total_values = <double*>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)
-67
View File
@@ -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)
+4 -4
View File
@@ -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()