diff --git a/doc/examples/plot_denoise.py b/doc/examples/plot_denoise.py index 8f02f21b..debb2ea6 100644 --- a/doc/examples/plot_denoise.py +++ b/doc/examples/plot_denoise.py @@ -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) diff --git a/skimage/filter/__init__.py b/skimage/filter/__init__.py index f1c1fd49..f32c3b23 100644 --- a/skimage/filter/__init__.py +++ b/skimage/filter/__init__.py @@ -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 diff --git a/skimage/filter/denoise.py b/skimage/filter/_denoise.py similarity index 65% rename from skimage/filter/denoise.py rename to skimage/filter/_denoise.py index 351b2482..e26f8387 100644 --- a/skimage/filter/denoise.py +++ b/skimage/filter/_denoise.py @@ -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) diff --git a/skimage/filter/_denoise.pyx b/skimage/filter/_denoise_cy.pyx similarity index 56% rename from skimage/filter/_denoise.pyx rename to skimage/filter/_denoise_cy.pyx index b60a5a85..8ebf5c73 100644 --- a/skimage/filter/_denoise.pyx +++ b/skimage/filter/_denoise_cy.pyx @@ -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]) diff --git a/skimage/filter/setup.py b/skimage/filter/setup.py index 56c1e9e5..650a26a6 100644 --- a/skimage/filter/setup.py +++ b/skimage/filter/setup.py @@ -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()]) diff --git a/skimage/filter/tests/test_denoise.py b/skimage/filter/tests/test_denoise.py index 63a5a5e0..cb3b5da3 100644 --- a/skimage/filter/tests/test_denoise.py +++ b/skimage/filter/tests/test_denoise.py @@ -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():