wassname/scikit-image
1
0
Fork 0
mirror of https://github.com/wassname/scikit-image.git synced 2026-07-03 13:11:14 +08:00
Code Issues Packages Projects Releases Wiki Activity

Add bilateral denoising filter

Browse Source
This commit is contained in:
Johannes Schönberger
2012-10-06 00:34:05 +02:00
parent 75d706ca25
commit f9ecd140d8
4 changed files with 221 additions and 1 deletions
Download Patch File Download Diff File Expand all files Collapse all files
skimage/filter/__init__.py
+1 -1
View File
@@ -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
skimage/filter/_denoise.pyx
+160
View File
@@ -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 = <double*>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 = <double*>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 = <double*>image.data
cdef double* out_data = <double*>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[<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):
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 = <double*>image.data
cdef double* out_data = <double*>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 = <double*>malloc(dims * sizeof(double))
cdef double* centres = <double*>malloc(dims * sizeof(double))
cdef double* total_values = <double*>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[<int>(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
skimage/filter/denoise.py
+57
View File
@@ -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)
skimage/filter/setup.py
+3
View File
@@ -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