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

Merge pull request #372 from ahojnnes/tv-filter

Browse Source 
ENH: New implementation of TV denoising.
This commit is contained in:
Stefan van der Walt
2012-12-28 12:16:17 -08:00
parent eff906bfb3 81d88b8cc3
commit c3b1984190
6 changed files with 279 additions and 69 deletions
Download Patch File Download Diff File Expand all files Collapse all files
doc/examples/plot_denoise.py
+8 -6
View File
@@ -30,30 +30,32 @@ 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
from skimage.filter import denoise_tv_chambolle, 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 = lena + 0.6 * lena.std() * np.random.random(lena.shape)
noisy = np.clip(noisy, 0, 1)
fig, ax = plt.subplots(nrows=2, ncols=3, figsize=(8, 5))
plt.gray()
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].imshow(denoise_tv_chambolle(noisy, weight=0.1, multichannel=True))
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].imshow(denoise_bilateral(noisy, sigma_range=0.05, 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].imshow(denoise_tv_chambolle(noisy, weight=0.2, multichannel=True))
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].imshow(denoise_bilateral(noisy, sigma_range=0.1, sigma_spatial=15))
ax[1, 1].axis('off')
ax[1, 1].set_title('(more) Bilateral')
ax[1, 2].imshow(lena)
skimage/filter/__init__.py
+2 -2
View File
@@ -3,7 +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
from ._denoise import denoise_bilateral
from ._denoise import denoise_tv_chambolle, tv_denoise
from ._denoise_cy import denoise_bilateral, denoise_tv_bregman
from ._rank_order import rank_order
from .thresholding import threshold_otsu, threshold_adaptive
skimage/filter/denoise.py → skimage/filter/_denoise.py
+68 -50
View File
@@ -3,28 +3,28 @@ from skimage import img_as_float
from skimage._shared.utils import deprecated
def _denoise_tv_3d(im, weight=100, eps=2.e-4, n_iter_max=200):
"""Perform total-variation denoising on 3-D arrays.
def _denoise_tv_chambolle_3d(im, weight=100, eps=2.e-4, n_iter_max=200):
"""Perform total-variation denoising on 3D images.
Parameters
----------
im: ndarray
im : ndarray
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``).
eps: float, optional
weight : float, optional
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
the stop criterion. The algorithm stops when:
(E_(n-1) - E_n) < eps * E_0
n_iter_max: int, optional
n_iter_max : int, optional
Maximal number of iterations used for the optimization.
Returns
-------
out: ndarray
out : ndarray
Denoised array of floats.
Notes
@@ -33,13 +33,14 @@ def _denoise_tv_3d(im, weight=100, eps=2.e-4, n_iter_max=200):
Examples
---------
First build synthetic noisy 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 = (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 = denoise_tv_3d(mask, weight=100)
>>> mask += 0.2 * np.random.randn(*mask.shape)
>>> res = denoise_tv(mask, weight=100)
"""
px = np.zeros_like(im)
py = np.zeros_like(im)
pz = np.zeros_like(im)
@@ -83,28 +84,28 @@ def _denoise_tv_3d(im, weight=100, eps=2.e-4, n_iter_max=200):
return out
def _denoise_tv_2d(im, weight=50, eps=2.e-4, n_iter_max=200):
"""Perform total-variation denoising.
def _denoise_tv_chambolle_2d(im, weight=50, eps=2.e-4, n_iter_max=200):
"""Perform total-variation denoising on 2D images.
Parameters
----------
im: ndarray
im : ndarray
Input data to be denoised.
weight: float, optional
Denoising weight. The greater ``weight``, the more denoising (at
the expense of fidelity to ``input``)
eps: float, optional
weight : float, optional
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
the stop criterion. The algorithm stops when:
(E_(n-1) - E_n) < eps * E_0
n_iter_max: int, optional
n_iter_max : int, optional
Maximal number of iterations used for the optimization.
Returns
-------
out: ndarray
out : ndarray
Denoised array of floats.
Notes
@@ -123,13 +124,13 @@ def _denoise_tv_2d(im, weight=50, eps=2.e-4, n_iter_max=200):
Examples
---------
>>> import scipy
>>> lena = scipy.lena()
>>> import scipy
>>> lena = scipy.lena().astype(np.float)
>>> lena += 0.5 * lena.std()*np.random.randn(*lena.shape)
>>> denoised_lena = denoise_tv(lena, weight=60.0)
>>> from skimage import color, data
>>> lena = color.rgb2gray(data.lena())
>>> lena += 0.5 * lena.std() * np.random.randn(*lena.shape)
>>> denoised_lena = denoise_tv(lena, weight=60)
"""
px = np.zeros_like(im)
py = np.zeros_like(im)
gx = np.zeros_like(im)
@@ -166,34 +167,41 @@ def _denoise_tv_2d(im, weight=50, eps=2.e-4, n_iter_max=200):
return out
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.
def denoise_tv_chambolle(im, weight=50, eps=2.e-4, n_iter_max=200,
multichannel=False):
"""Perform total-variation denoising on 2D and 3D images.
Parameters
----------
im: ndarray (2d or 3d) of ints, uints or floats
im : ndarray (2d or 3d) of ints, uints or floats
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``).
eps: float, optional
weight : float, optional
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 the stop criterion. The algorithm stops when:
(E_(n-1) - E_n) < eps * E_0
n_iter_max: int, optional
n_iter_max : int, optional
Maximal number of iterations used for the optimization.
multichannel : bool, optional
Apply total-variation denoising separately for each channel. This
option should be true for color images, otherwise the denoising is
also applied in the 3rd dimension.
Returns
-------
out: ndarray
Denoised array of floats.
out : ndarray
Denoised image.
Notes
-----
Make sure to set the multichannel parameter appropriately for color images.
The principle of total variation denoising is explained in
http://en.wikipedia.org/wiki/Total_variation_denoising
@@ -214,32 +222,42 @@ def denoise_tv(im, weight=50, eps=2.e-4, n_iter_max=200):
Examples
---------
>>> import scipy
>>> # 2D example using lena
>>> lena = scipy.lena()
>>> import scipy
>>> lena = scipy.lena().astype(np.float)
>>> lena += 0.5 * lena.std()*np.random.randn(*lena.shape)
2D example on Lena image:
>>> from skimage import color, data
>>> lena = color.rgb2gray(data.lena())
>>> lena += 0.5 * lena.std() * np.random.randn(*lena.shape)
>>> denoised_lena = denoise_tv(lena, weight=60)
>>> # 3D example on synthetic data
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 = (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 = denoise_tv_3d(mask, weight=100)
>>> res = denoise_tv(mask, weight=100)
"""
im_type = im.dtype
if not im_type.kind == 'f':
im = img_as_float(im)
if im.ndim == 2:
out = _denoise_tv_2d(im, weight, eps, n_iter_max)
out = _denoise_tv_chambolle_2d(im, weight, eps, n_iter_max)
elif im.ndim == 3:
out = _denoise_tv_3d(im, weight, eps, n_iter_max)
if multichannel:
out = np.zeros_like(im)
for c in range(im.shape[2]):
out[..., c] = _denoise_tv_chambolle_2d(im[..., c], weight, eps,
n_iter_max)
else:
out = _denoise_tv_chambolle_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)
tv_denoise = deprecated('skimage.filter.denoise_tv_chambolle')\
(denoise_tv_chambolle)
skimage/filter/_denoise.pyx → skimage/filter/_denoise_cy.pyx
+144 -1
View File
@@ -7,8 +7,10 @@ cimport numpy as cnp
import numpy as np
from libc.math cimport exp, fabs, sqrt
from libc.stdlib cimport malloc, free
from libc.float cimport DBL_MAX
from skimage._shared.interpolation cimport get_pixel3d
from skimage.util import img_as_float
from skimage._shared.utils import deprecated
cdef inline double _gaussian_weight(double sigma, double value):
@@ -174,4 +176,145 @@ def denoise_bilateral(image, int win_size=5, sigma_range=None,
free(centres)
free(total_values)
return out
return np.squeeze(out)
def denoise_tv_bregman(image, double weight, int max_iter=100, double eps=1e-3):
"""Perform total-variation denoising using split-Bregman optimization.
Total-variation denoising (also know as total-variation regularization)
tries to find an image with less total-variation under the constraint
of being similar to the input image, which is controlled by the
regularization parameter.
Parameters
----------
image : ndarray
Input data to be denoised (converted using img_as_float`).
weight : float, optional
Denoising weight. The smaller the `weight`, the more denoising (at
the expense of less similarity to the `input`). The regularization
parameter `lambda` is chosen as `2 * weight`.
eps : float, optional
Relative difference of the value of the cost function that determines
the stop criterion. The algorithm stops when::
SUM((u(n) - u(n-1))**2) < eps
max_iter: int, optional
Maximal number of iterations used for the optimization.
Returns
-------
u : ndarray
Denoised image.
References
----------
.. [1] http://en.wikipedia.org/wiki/Total_variation_denoising
.. [2] Tom Goldstein and Stanley Osher, "The Split Bregman Method For L1
Regularized Problems",
ftp://ftp.math.ucla.edu/pub/camreport/cam08-29.pdf
.. [3] Pascal Getreuer, "Rudin–Osher–Fatemi Total Variation Denoising
using Split Bregman" in Image Processing On Line on 2012–05–19,
http://www.ipol.im/pub/art/2012/g-tvd/article_lr.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 rows2 = rows + 2
Py_ssize_t cols2 = cols + 2
Py_ssize_t r, c, k
Py_ssize_t total = rows * cols * dims
shape_ext = (rows2, cols2, dims)
cnp.ndarray[dtype=cnp.double_t, ndim=3, mode='c'] cimage = \
np.ascontiguousarray(image)
cnp.ndarray[dtype=cnp.double_t, ndim=3, mode='c'] u = \
np.zeros(shape_ext, dtype=np.double)
cnp.ndarray[dtype=cnp.double_t, ndim=3, mode='c'] dx = \
np.zeros(shape_ext, dtype=np.double)
cnp.ndarray[dtype=cnp.double_t, ndim=3, mode='c'] dy = \
np.zeros(shape_ext, dtype=np.double)
cnp.ndarray[dtype=cnp.double_t, ndim=3, mode='c'] bx = \
np.zeros(shape_ext, dtype=np.double)
cnp.ndarray[dtype=cnp.double_t, ndim=3, mode='c'] by = \
np.zeros(shape_ext, dtype=np.double)
double ux, uy, uprev, unew, bxx, byy, dxx, dyy, s
int i = 0
double lam = 2 * weight
double rmse = DBL_MAX
double norm = (weight + 4 * lam)
u[1:-1, 1:-1] = image
# reflect image
u[0, 1:-1] = image[1, :]
u[1:-1, 0] = image[:, 1]
u[-1, 1:-1] = image[-2, :]
u[1:-1, -1] = image[:, -2]
while i < max_iter and rmse > eps:
rmse = 0
for k in range(dims):
for r in range(1, rows + 1):
for c in range(1, cols + 1):
uprev = u[r, c, k]
# forward derivatives
ux = u[r, c + 1, k] - uprev
uy = u[r + 1, c, k] - uprev
# Gauss-Seidel method
unew = (
lam * (
+ u[r + 1, c, k]
+ u[r - 1, c, k]
+ u[r, c + 1, k]
+ u[r, c - 1, k]
+ dx[r, c - 1, k]
- dx[r, c, k]
+ dy[r - 1, c, k]
- dy[r, c, k]
- bx[r, c - 1, k]
+ bx[r, c, k]
- by[r - 1, c, k]
+ by[r, c, k]
) + weight * cimage[r - 1, c - 1, k]
) / norm
u[r, c, k] = unew
# update root mean square error
rmse += (unew - uprev)**2
bxx = bx[r, c, k]
byy = by[r, c, k]
s = sqrt((ux + bxx)**2 + (uy + byy)**2)
dxx = s * lam * (ux + bxx) / (s * lam + 1)
dyy = s * lam * (uy + byy) / (s * lam + 1)
dx[r, c, k] = dxx
dy[r, c, k] = dyy
bx[r, c, k] += ux - dxx
by[r, c, k] += uy - dyy
rmse = sqrt(rmse / total)
i += 1
return np.squeeze(u[1:-1, 1:-1])
skimage/filter/setup.py
+2 -2
View File
@@ -13,7 +13,7 @@ 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)
cython(['_denoise_cy.pyx'], working_path=base_path)
cython(['rank/_core8.pyx'], working_path=base_path)
cython(['rank/_core16.pyx'], working_path=base_path)
cython(['rank/_crank8.pyx'], working_path=base_path)
@@ -27,7 +27,7 @@ def configuration(parent_package='', top_path=None):
config.add_extension('_ctmf', sources=['_ctmf.c'],
include_dirs=[get_numpy_include_dirs()])
config.add_extension('_denoise', sources=['_denoise.c'],
config.add_extension('_denoise_cy', sources=['_denoise_cy.c'],
include_dirs=[get_numpy_include_dirs(), '../_shared'])
config.add_extension('rank/_core8', sources=['rank/_core8.c'],
include_dirs=[get_numpy_include_dirs()])
skimage/filter/tests/test_denoise.py
+55 -8
View File
@@ -1,5 +1,5 @@
import numpy as np
from numpy.testing import run_module_suite, assert_raises
from numpy.testing import run_module_suite, assert_raises, assert_equal
from skimage import filter, data, color, img_as_float
@@ -8,7 +8,7 @@ lena = img_as_float(data.lena()[:256, :256])
lena_gray = color.rgb2gray(lena)
def test_denoise_tv_2d():
def test_denoise_tv_chambolle_2d():
# lena image
img = lena_gray
# add noise to lena
@@ -16,7 +16,7 @@ def test_denoise_tv_2d():
# 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)
denoised_lena = filter.denoise_tv_chambolle(img, weight=60.0)
# which dtype?
assert denoised_lena.dtype in [np.float, np.float32, np.float64]
from scipy import ndimage
@@ -29,19 +29,25 @@ def test_denoise_tv_2d():
< np.sqrt((grad**2).sum()) / 2)
def test_denoise_tv_float_result_range():
def test_denoise_tv_chambolle_multichannel():
denoised0 = filter.denoise_tv_chambolle(lena[..., 0], weight=60.0)
denoised = filter.denoise_tv_chambolle(lena, weight=60.0, multichannel=True)
assert_equal(denoised[..., 0], denoised0)
def test_denoise_tv_chambolle_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)
denoised_int_lena = filter.denoise_tv_chambolle(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():
def test_denoise_tv_chambolle_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
@@ -50,12 +56,53 @@ def test_denoise_tv_3d():
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)
res = filter.denoise_tv_chambolle(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)))
assert_raises(ValueError, filter.denoise_tv_chambolle,
np.random.random((8, 8, 8, 8)))
def test_denoise_tv_bregman_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_tv_bregman(img, weight=10)
out2 = filter.denoise_tv_bregman(img, weight=5)
# make sure noise is reduced
assert img.std() > out1.std()
assert out1.std() > out2.std()
def test_denoise_tv_bregman_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_bregman(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_bregman_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_tv_bregman(img, weight=10)
out2 = filter.denoise_tv_bregman(img, weight=5)
# make sure noise is reduced
assert img.std() > out1.std()
assert out1.std() > out2.std()
def test_denoise_bilateral_2d():