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

Merge pull request #349 from ahojnnes/bilateral-filter

Browse Source 
ENH: Implement bilateral filter.
This commit is contained in:
Stefan van der Walt
2012-10-18 12:34:51 -07:00
parent c3a37a87ea 8347bd449f
commit fc998310ef
11 changed files with 426 additions and 176 deletions
Download Patch File Download Diff File Expand all files Collapse all files
bento.info
+4 -1
View File
@@ -64,7 +64,10 @@ Library:
Extension: skimage.filter._ctmf
Sources:
skimage/filter/_ctmf.pyx
Extension: skimage.morphology.ccomp
Extension: skimage.filter._denoise
Sources:
skimage/filter/_ctmf.pyx
Extension: skimage.morphology.denoise
Sources:
skimage/morphology/ccomp.pyx
Extension: skimage.morphology._watershed
doc/examples/plot_denoise.py
+66
View File
@@ -0,0 +1,66 @@
"""
=============================
Denoising the picture of Lena
=============================
In this example, we denoise a noisy version of the picture of Lena using the
total variation and bilateral denoising filter.
These algorithms typically produce "posterized" images with flat domains
separated by sharp edges. It is possible to change the degree of posterization
by controlling the tradeoff between denoising and faithfulness to the original
image.
Total variation filter
----------------------
The result of this filter is an image that has a minimal total variation norm,
while being as close to the initial image as possible. The total variation is
the L1 norm of the gradient of the image.
Bilateral filter
----------------
A bilateral filter is an edge-preserving and noise reducing filter. It averages
pixels based on their spatial closeness and radiometric similarity.
"""
import numpy as np
import matplotlib.pyplot as plt
from skimage import data, color, img_as_float
from skimage.filter import denoise_tv, denoise_bilateral
lena = img_as_float(data.lena())
lena = lena[220:300, 220:320]
noisy = lena + 0.5 * lena.std() * np.random.random(lena.shape)
noisy = np.clip(noisy, 0, 1)
fig, ax = plt.subplots(nrows=2, ncols=3, figsize=(8, 5))
ax[0, 0].imshow(noisy)
ax[0, 0].axis('off')
ax[0, 0].set_title('noisy')
ax[0, 1].imshow(denoise_tv(noisy, weight=0.1))
ax[0, 1].axis('off')
ax[0, 1].set_title('TV')
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(denoise_tv(noisy, weight=0.2))
ax[1, 0].axis('off')
ax[1, 0].set_title('(more) TV')
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)
ax[1, 2].axis('off')
ax[1, 2].set_title('original')
fig.subplots_adjust(wspace=0.02, hspace=0.2,
top=0.9, bottom=0.05, left=0, right=1)
plt.show()
doc/examples/plot_lena_tv_denoise.py
-51
View File
@@ -1,51 +0,0 @@
"""
====================================================
Denoising the picture of Lena using total variation
====================================================
In this example, we denoise a noisy version of the picture of Lena
using the total variation denoising filter. The result of this filter
is an image that has a minimal total variation norm, while being as
close to the initial image as possible. The total variation is the L1
norm of the gradient of the image, and minimizing the total variation
typically produces "posterized" images with flat domains separated by
sharp edges.
It is possible to change the degree of posterization by controlling
the tradeoff between denoising and faithfulness to the original image.
"""
import numpy as np
import matplotlib.pyplot as plt
from skimage import data, color, img_as_ubyte
from skimage.filter import tv_denoise
l = img_as_ubyte(color.rgb2gray(data.lena()))
l = l[230:290, 220:320]
noisy = l + 0.4 * l.std() * np.random.random(l.shape)
tv_denoised = tv_denoise(noisy, weight=10)
plt.figure(figsize=(8, 2))
plt.subplot(131)
plt.imshow(noisy, cmap=plt.cm.gray, vmin=40, vmax=220)
plt.axis('off')
plt.title('noisy', fontsize=20)
plt.subplot(132)
plt.imshow(tv_denoised, cmap=plt.cm.gray, vmin=40, vmax=220)
plt.axis('off')
plt.title('TV denoising', fontsize=20)
tv_denoised = tv_denoise(noisy, weight=50)
plt.subplot(133)
plt.imshow(tv_denoised, cmap=plt.cm.gray, vmin=40, vmax=220)
plt.axis('off')
plt.title('(more) TV denoising', fontsize=20)
plt.subplots_adjust(wspace=0.02, hspace=0.02, top=0.9, bottom=0, left=0,
right=1)
plt.show()
skimage/_shared/interpolation.pxd
+5 -2
View File
@@ -18,7 +18,10 @@ cdef double bicubic_interpolation(double* image, int rows, int cols,
double r, double c, char mode,
double cval)
cdef double get_pixel(double* image, int rows, int cols, int r, int c,
char mode, double cval)
cdef double get_pixel2d(double* image, int rows, int cols, int r, int c,
char mode, double cval)
cdef double get_pixel3d(double* image, int rows, int cols, int dims, int r,
int c, int d, char mode, double cval)
cdef int coord_map(int dim, int coord, char mode)
skimage/_shared/interpolation.pyx
+44 -10
View File
@@ -35,8 +35,8 @@ cdef inline double nearest_neighbour_interpolation(double* image, int rows,
"""
return get_pixel(image, rows, cols, <int>round(r), <int>round(c),
mode, cval)
return get_pixel2d(image, rows, cols, <int>round(r), <int>round(c),
mode, cval)
cdef inline double bilinear_interpolation(double* image, int rows, int cols,
@@ -72,10 +72,10 @@ cdef inline double bilinear_interpolation(double* image, int rows, int cols,
maxc = <int>ceil(c)
dr = r - minr
dc = c - minc
top = (1 - dc) * get_pixel(image, rows, cols, minr, minc, mode, cval) \
+ dc * get_pixel(image, rows, cols, minr, maxc, mode, cval)
bottom = (1 - dc) * get_pixel(image, rows, cols, maxr, minc, mode, cval) \
+ dc * get_pixel(image, rows, cols, maxr, maxc, mode, cval)
top = (1 - dc) * get_pixel2d(image, rows, cols, minr, minc, mode, cval) \
+ dc * get_pixel2d(image, rows, cols, minr, maxc, mode, cval)
bottom = (1 - dc) * get_pixel2d(image, rows, cols, maxr, minc, mode, cval) \
+ dc * get_pixel2d(image, rows, cols, maxr, maxc, mode, cval)
return (1 - dr) * top + dr * bottom
@@ -144,7 +144,7 @@ cdef inline double biquadratic_interpolation(double* image, int rows, int cols,
# row-wise cubic interpolation
for pr in range(r0, r0 + 3):
for pc in range(c0, c0 + 3):
fc[pc - c0] = get_pixel(image, rows, cols, pr, pc, mode, cval)
fc[pc - c0] = get_pixel2d(image, rows, cols, pr, pc, mode, cval)
fr[pr - r0] = quadratic_interpolation(xc, fc)
# cubic interpolation for interpolated values of each row
@@ -216,15 +216,15 @@ cdef inline double bicubic_interpolation(double* image, int rows, int cols,
# row-wise cubic interpolation
for pr in range(r0, r0 + 4):
for pc in range(c0, c0 + 4):
fc[pc - c0] = get_pixel(image, rows, cols, pr, pc, mode, cval)
fc[pc - c0] = get_pixel2d(image, rows, cols, pr, pc, mode, cval)
fr[pr - r0] = cubic_interpolation(xc, fc)
# cubic interpolation for interpolated values of each row
return cubic_interpolation(xr, fr)
cdef inline double get_pixel(double* image, int rows, int cols, int r, int c,
char mode, double cval):
cdef inline double get_pixel2d(double* image, int rows, int cols, int r, int c,
char mode, double cval):
"""Get a pixel from the image, taking wrapping mode into consideration.
Parameters
@@ -255,6 +255,40 @@ cdef inline double get_pixel(double* image, int rows, int cols, int r, int c,
return image[coord_map(rows, r, mode) * cols + coord_map(cols, c, mode)]
cdef inline double get_pixel3d(double* image, int rows, int cols, int dims, int r,
int c, int d, char mode, double cval):
"""Get a pixel from the image, taking wrapping mode into consideration.
Parameters
----------
image : double array
Input image.
rows, cols, dims : int
Shape of image.
r, c, d : int
Position at which to get the pixel.
mode : {'C', 'W', 'R', 'N'}
Wrapping mode. Constant, Wrap, Reflect or Nearest.
cval : double
Constant value to use for constant mode.
Returns
-------
value : double
Pixel value at given position.
"""
if mode == 'C':
if (r < 0) or (r > rows - 1) or (c < 0) or (c > cols - 1):
return cval
else:
return image[r * cols * dims + c * dims + d]
else:
return image[coord_map(rows, r, mode) * cols * dims
+ coord_map(cols, c, mode) * dims
+ d]
cdef inline int coord_map(int dim, int coord, char mode):
"""
Wrap a coordinate, according to a given mode.
skimage/filter/__init__.py
+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 ._tv_denoise import tv_denoise
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
skimage/filter/_denoise.pyx
+177
View File
@@ -0,0 +1,177 @@
#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_pixel3d
from skimage.util import img_as_float
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))
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))
Py_ssize_t kr, kc
Py_ssize_t window_ext = (win_size - 1) / 2
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_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.
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_spatial`).
Radiometric similarity is measured by the gaussian function of the euclidian
distance between two color values and a certain standard deviation
(`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. Note, that the image will be converted using
the `img_as_float` function and thus the standard deviation is in
respect to the range `[0, 1]`.
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.
Returns
-------
denoised : ndarray
Denoised image.
References
----------
.. [1] http://users.soe.ucsc.edu/~manduchi/Papers/ICCV98.pdf
"""
image = np.atleast_3d(img_as_float(image))
cdef:
Py_ssize_t rows = image.shape[0]
Py_ssize_t cols = image.shape[1]
Py_ssize_t dims = image.shape[2]
Py_ssize_t window_ext = (win_size - 1) / 2
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*>cimage.data
double* out_data = <double*>out.data
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, 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
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, cmode, 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/_tv_denoise.py → skimage/filter/denoise.py
+35 -42
View File
@@ -1,37 +1,35 @@
import numpy as np
from skimage import img_as_float
from skimage._shared.utils import deprecated
def _tv_denoise_3d(im, weight=100, eps=2.e-4, n_iter_max=200):
"""
Perform total-variation denoising on 3-D arrays
def _denoise_tv_3d(im, weight=100, eps=2.e-4, n_iter_max=200):
"""Perform total-variation denoising on 3-D arrays.
Parameters
----------
im: ndarray
3-D input data to be denoised
3-D input data to be denoised.
weight: float, optional
denoising weight. The greater ``weight``, the more denoising (at
the expense of fidelity to ``input``)
Denoising weight. The greater ``weight``, the more denoising (at
the expense of fidelity to ``input``).
eps: float, optional
relative difference of the value of the cost function that determines
Relative difference of the value of the cost function that determines
the stop criterion. The algorithm stops when:
(E_(n-1) - E_n) < eps * E_0
n_iter_max: int, optional
maximal number of iterations used for the optimization.
Maximal number of iterations used for the optimization.
Returns
-------
out: ndarray
denoised array of floats
Denoised array of floats.
Notes
-----
Rudin, Osher and Fatemi algorithm
Rudin, Osher and Fatemi algorithm.
Examples
---------
@@ -40,7 +38,7 @@ def _tv_denoise_3d(im, weight=100, eps=2.e-4, n_iter_max=200):
>>> mask = (x -22)**2 + (y - 20)**2 + (z - 17)**2 < 8**2
>>> mask = mask.astype(np.float)
>>> mask += 0.2*np.random.randn(*mask.shape)
>>> res = tv_denoise_3d(mask, weight=100)
>>> res = denoise_tv_3d(mask, weight=100)
"""
px = np.zeros_like(im)
py = np.zeros_like(im)
@@ -85,44 +83,40 @@ def _tv_denoise_3d(im, weight=100, eps=2.e-4, n_iter_max=200):
return out
def _tv_denoise_2d(im, weight=50, eps=2.e-4, n_iter_max=200):
"""
Perform total-variation denoising
def _denoise_tv_2d(im, weight=50, eps=2.e-4, n_iter_max=200):
"""Perform total-variation denoising.
Parameters
----------
im: ndarray
input data to be denoised
Input data to be denoised.
weight: float, optional
denoising weight. The greater ``weight``, the more denoising (at
Denoising weight. The greater ``weight``, the more denoising (at
the expense of fidelity to ``input``)
eps: float, optional
relative difference of the value of the cost function that determines
Relative difference of the value of the cost function that determines
the stop criterion. The algorithm stops when:
(E_(n-1) - E_n) < eps * E_0
n_iter_max: int, optional
maximal number of iterations used for the optimization.
Maximal number of iterations used for the optimization.
Returns
-------
out: ndarray
denoised array of floats
Denoised array of floats.
Notes
-----
The principle of total variation denoising is explained in
http://en.wikipedia.org/wiki/Total_variation_denoising
http://en.wikipedia.org/wiki/Total_variation_denoising.
This code is an implementation of the algorithm of Rudin, Fatemi and Osher
that was proposed by Chambolle in [1]_.
References
----------
.. [1] A. Chambolle, An algorithm for total variation minimization and
applications, Journal of Mathematical Imaging and Vision,
Springer, 2004, 20, 89-97.
@@ -134,7 +128,7 @@ def _tv_denoise_2d(im, weight=50, eps=2.e-4, n_iter_max=200):
>>> import scipy
>>> lena = scipy.lena().astype(np.float)
>>> lena += 0.5 * lena.std()*np.random.randn(*lena.shape)
>>> denoised_lena = tv_denoise(lena, weight=60.0)
>>> denoised_lena = denoise_tv(lena, weight=60.0)
"""
px = np.zeros_like(im)
py = np.zeros_like(im)
@@ -172,34 +166,31 @@ def _tv_denoise_2d(im, weight=50, eps=2.e-4, n_iter_max=200):
return out
def tv_denoise(im, weight=50, eps=2.e-4, n_iter_max=200):
"""
Perform total-variation denoising on 2-d and 3-d images
def denoise_tv(im, weight=50, eps=2.e-4, n_iter_max=200):
"""Perform total-variation denoising on 2-d and 3-d images.
Parameters
----------
im: ndarray (2d or 3d) of ints, uints or floats
input data to be denoised. `im` can be of any numeric type,
Input data to be denoised. `im` can be of any numeric type,
but it is cast into an ndarray of floats for the computation
of the denoised image.
weight: float, optional
denoising weight. The greater ``weight``, the more denoising (at
the expense of fidelity to ``input``)
Denoising weight. The greater ``weight``, the more denoising (at
the expense of fidelity to ``input``).
eps: float, optional
relative difference of the value of the cost function that
Relative difference of the value of the cost function that
determines the stop criterion. The algorithm stops when:
(E_(n-1) - E_n) < eps * E_0
n_iter_max: int, optional
maximal number of iterations used for the optimization.
Maximal number of iterations used for the optimization.
Returns
-------
out: ndarray
denoised array of floats
Denoised array of floats.
Notes
-----
@@ -217,7 +208,6 @@ def tv_denoise(im, weight=50, eps=2.e-4, n_iter_max=200):
References
----------
.. [1] A. Chambolle, An algorithm for total variation minimization and
applications, Journal of Mathematical Imaging and Vision,
Springer, 2004, 20, 89-97.
@@ -230,23 +220,26 @@ def tv_denoise(im, weight=50, eps=2.e-4, n_iter_max=200):
>>> import scipy
>>> lena = scipy.lena().astype(np.float)
>>> lena += 0.5 * lena.std()*np.random.randn(*lena.shape)
>>> denoised_lena = tv_denoise(lena, weight=60)
>>> denoised_lena = denoise_tv(lena, weight=60)
>>> # 3D example on synthetic data
>>> x, y, z = np.ogrid[0:40, 0:40, 0:40]
>>> mask = (x -22)**2 + (y - 20)**2 + (z - 17)**2 < 8**2
>>> mask = mask.astype(np.float)
>>> mask += 0.2*np.random.randn(*mask.shape)
>>> res = tv_denoise_3d(mask, weight=100)
>>> res = denoise_tv_3d(mask, weight=100)
"""
im_type = im.dtype
if not im_type.kind == 'f':
im = img_as_float(im)
if im.ndim == 2:
out = _tv_denoise_2d(im, weight, eps, n_iter_max)
out = _denoise_tv_2d(im, weight, eps, n_iter_max)
elif im.ndim == 3:
out = _tv_denoise_3d(im, weight, eps, n_iter_max)
out = _denoise_tv_3d(im, weight, eps, n_iter_max)
else:
raise ValueError('only 2-d and 3-d images may be denoised with this '
'function')
return out
tv_denoise = deprecated('skimage.filter.denoise_tv')(denoise_tv)
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
skimage/filter/tests/test_denoise.py
+90
View File
@@ -0,0 +1,90 @@
import numpy as np
from numpy.testing import run_module_suite, assert_raises
from skimage import filter, data, color, img_as_float
lena = img_as_float(data.lena()[:256, :256])
lena_gray = color.rgb2gray(lena)
def test_denoise_tv_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(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_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(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_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(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, np.random.random((8, 8, 8, 8)))
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()
if __name__ == "__main__":
run_module_suite()
skimage/filter/tests/test_tv_denoise.py
-69
View File
@@ -1,69 +0,0 @@
import numpy as np
from numpy.testing import run_module_suite
from skimage import filter, data, color
class TestTvDenoise():
def test_tv_denoise_2d(self):
"""
Apply the TV denoising algorithm on the lena image provided
by scipy
"""
# lena image
lena = color.rgb2gray(data.lena())[:256, :256]
# add noise to lena
lena += 0.5 * lena.std() * np.random.randn(*lena.shape)
# clip noise so that it does not exceed allowed range for float images.
lena = np.clip(lena, 0, 1)
# denoise
denoised_lena = filter.tv_denoise(lena, weight=60.0)
# which dtype?
assert denoised_lena.dtype in [np.float, np.float32, np.float64]
from scipy import ndimage
grad = ndimage.morphological_gradient(lena, 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_tv_denoise_float_result_range(self):
# lena image
lena = color.rgb2gray(data.lena())[:256, :256]
int_lena = np.multiply(lena, 255).astype(np.uint8)
assert np.max(int_lena) > 1
denoised_int_lena = filter.tv_denoise(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_tv_denoise_3d(self):
"""
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.randn(*mask.shape)
mask[mask < 0] = 0
mask[mask > 255] = 255
res = filter.tv_denoise(mask.astype(np.uint8), weight=100)
assert res.dtype == np.float
assert res.std() * 255 < mask.std()
# test wrong number of dimensions
a = np.random.random((8, 8, 8, 8))
try:
res = filter.tv_denoise(a)
except ValueError:
pass
if __name__ == "__main__":
run_module_suite()