diff --git a/skimage/restoration/_denoise.py b/skimage/restoration/_denoise.py index 97fab246..62985cc7 100644 --- a/skimage/restoration/_denoise.py +++ b/skimage/restoration/_denoise.py @@ -116,13 +116,13 @@ def denoise_tv_bregman(image, weight, max_iter=100, eps=1e-3, isotropic=True): return _denoise_tv_bregman(image, weight, max_iter, eps, isotropic) -def _denoise_tv_chambolle_3d(im, weight=0.2, eps=2.e-4, n_iter_max=200): - """Perform total-variation denoising on 3D images. +def _denoise_tv_chambolle_nd(im, weight=0.1, eps=2.e-4, n_iter_max=200): + """Perform total-variation denoising on n-dimensional images. Parameters ---------- im : ndarray - 3-D input data to be denoised. + n-D input data to be denoised. weight : float, optional Denoising weight. The greater `weight`, the more denoising (at the expense of fidelity to `input`). @@ -146,36 +146,45 @@ def _denoise_tv_chambolle_3d(im, weight=0.2, eps=2.e-4, n_iter_max=200): """ - px = np.zeros_like(im) - py = np.zeros_like(im) - pz = np.zeros_like(im) - gx = np.zeros_like(im) - gy = np.zeros_like(im) - gz = np.zeros_like(im) + ndim = im.ndim + p = np.zeros((im.ndim, ) + im.shape, dtype=im.dtype) + g = np.zeros_like(p) d = np.zeros_like(im) i = 0 while i < n_iter_max: - d = - px - py - pz - d[1:] += px[:-1] - d[:, 1:] += py[:, :-1] - d[:, :, 1:] += pz[:, :, :-1] - - out = im + d + if i > 0: + # d will be the (negative) divergence of p + d = -p.sum(0) + slices_d = [slice(None), ] * ndim + slices_p = [slice(None), ] * (ndim + 1) + for ax in range(ndim): + slices_d[ax] = slice(1, None) + slices_p[ax+1] = slice(0, -1) + slices_p[0] = ax + d[slices_d] += p[slices_p] + slices_d[ax] = slice(None) + slices_p[ax+1] = slice(None) + out = im + d + else: + out = im E = (d ** 2).sum() - gx[:-1] = np.diff(out, axis=0) - gy[:, :-1] = np.diff(out, axis=1) - gz[:, :, :-1] = np.diff(out, axis=2) - norm = np.sqrt(gx ** 2 + gy ** 2 + gz ** 2) + # g stores the gradients of out along each axis + # e.g. g[0] is the first order finite difference along axis 0 + slices_g = [slice(None), ] * (ndim + 1) + for ax in range(ndim): + slices_g[ax+1] = slice(0, -1) + slices_g[0] = ax + g[slices_g] = np.diff(out, axis=ax) + slices_g[ax+1] = slice(None) + + norm = np.sqrt((g ** 2).sum(axis=0))[np.newaxis, ...] E += weight * norm.sum() - norm *= 0.5 / weight + tau = 1. / (2.*ndim) + norm *= tau / weight norm += 1. - px -= 1. / 6. * gx - px /= norm - py -= 1. / 6. * gy - py /= norm - pz -= 1 / 6. * gz - pz /= norm + p -= tau * g + p /= norm E /= float(im.size) if i == 0: E_init = E @@ -189,89 +198,13 @@ def _denoise_tv_chambolle_3d(im, weight=0.2, eps=2.e-4, n_iter_max=200): return out -def _denoise_tv_chambolle_2d(im, weight=0.2, eps=2.e-4, n_iter_max=200): - """Perform total-variation denoising on 2D images. - - Parameters - ---------- - 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 - 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. - - Returns - ------- - out : ndarray - Denoised array of floats. - - Notes - ----- - The principle of total variation denoising is explained in - 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. - - """ - - px = np.zeros_like(im) - py = np.zeros_like(im) - gx = np.zeros_like(im) - gy = np.zeros_like(im) - d = np.zeros_like(im) - i = 0 - while i < n_iter_max: - d = -px - py - d[1:] += px[:-1] - d[:, 1:] += py[:, :-1] - - out = im + d - E = (d ** 2).sum() - gx[:-1] = np.diff(out, axis=0) - gy[:, :-1] = np.diff(out, axis=1) - norm = np.sqrt(gx ** 2 + gy ** 2) - E += weight * norm.sum() - norm *= 0.5 / weight - norm += 1 - px -= 0.25 * gx - px /= norm - py -= 0.25 * gy - py /= norm - E /= float(im.size) - if i == 0: - E_init = E - E_previous = E - else: - if np.abs(E_previous - E) < eps * E_init: - break - else: - E_previous = E - i += 1 - return out - - -def denoise_tv_chambolle(im, weight=0.2, eps=2.e-4, n_iter_max=200, +def denoise_tv_chambolle(im, weight=0.1, eps=2.e-4, n_iter_max=200, multichannel=False): - """Perform total-variation denoising on 2D and 3D images. + """Perform total-variation denoising on n-dimensional images. Parameters ---------- - im : ndarray (2d or 3d) of ints, uints or floats + im : ndarray 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. @@ -289,7 +222,7 @@ def denoise_tv_chambolle(im, weight=0.2, eps=2.e-4, n_iter_max=200, 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. + also applied in the channels dimension. Returns ------- @@ -341,17 +274,11 @@ def denoise_tv_chambolle(im, weight=0.2, eps=2.e-4, n_iter_max=200, if not im_type.kind == 'f': im = img_as_float(im) - if im.ndim == 2: - out = _denoise_tv_chambolle_2d(im, weight, eps, n_iter_max) - elif im.ndim == 3: - 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) + if multichannel: + out = np.zeros_like(im) + for c in range(im.shape[-1]): + out[..., c] = _denoise_tv_chambolle_nd(im[..., c], weight, eps, + n_iter_max) else: - raise ValueError('only 2-d and 3-d images may be denoised with this ' - 'function') + out = _denoise_tv_chambolle_nd(im, weight, eps, n_iter_max) return out diff --git a/skimage/restoration/tests/test_denoise.py b/skimage/restoration/tests/test_denoise.py index fb199afd..58f0a261 100644 --- a/skimage/restoration/tests/test_denoise.py +++ b/skimage/restoration/tests/test_denoise.py @@ -1,7 +1,7 @@ import numpy as np from numpy.testing import run_module_suite, assert_raises, assert_equal -from skimage import restoration, data, color, img_as_float +from skimage import restoration, data, color, img_as_float, measure np.random.seed(1234) @@ -21,7 +21,7 @@ def test_denoise_tv_chambolle_2d(): # clip noise so that it does not exceed allowed range for float images. img = np.clip(img, 0, 1) # denoise - denoised_astro = restoration.denoise_tv_chambolle(img, weight=0.25) + denoised_astro = restoration.denoise_tv_chambolle(img, weight=0.1) # which dtype? assert denoised_astro.dtype in [np.float, np.float32, np.float64] from scipy import ndimage as ndi @@ -34,8 +34,17 @@ def test_denoise_tv_chambolle_2d(): def test_denoise_tv_chambolle_multichannel(): - denoised0 = restoration.denoise_tv_chambolle(astro[..., 0], weight=0.25) - denoised = restoration.denoise_tv_chambolle(astro, weight=0.25, + denoised0 = restoration.denoise_tv_chambolle(astro[..., 0], weight=0.1) + denoised = restoration.denoise_tv_chambolle(astro, weight=0.1, + multichannel=True) + assert_equal(denoised[..., 0], denoised0) + + # tile astronaut subset to generate 3D+channels data + astro3 = np.tile(astro[:64, :64, np.newaxis, :], [1, 1, 2, 1]) + # modify along tiled dimension to give non-zero gradient on 3rd axis + astro3[:, :, 0, :] = 2*astro3[:, :, 0, :] + denoised0 = restoration.denoise_tv_chambolle(astro3[..., 0], weight=0.1) + denoised = restoration.denoise_tv_chambolle(astro3, weight=0.1, multichannel=True) assert_equal(denoised[..., 0], denoised0) @@ -46,7 +55,7 @@ def test_denoise_tv_chambolle_float_result_range(): int_astro = np.multiply(img, 255).astype(np.uint8) assert np.max(int_astro) > 1 denoised_int_astro = restoration.denoise_tv_chambolle(int_astro, - weight=0.25) + weight=0.1) # test if the value range of output float data is within [0.0:1.0] assert denoised_int_astro.dtype == np.float assert np.max(denoised_int_astro) <= 1.0 @@ -62,13 +71,45 @@ def test_denoise_tv_chambolle_3d(): mask += 20 * np.random.rand(*mask.shape) mask[mask < 0] = 0 mask[mask > 255] = 255 - res = restoration.denoise_tv_chambolle(mask.astype(np.uint8), weight=0.4) + res = restoration.denoise_tv_chambolle(mask.astype(np.uint8), weight=0.1) assert res.dtype == np.float assert res.std() * 255 < mask.std() - # test wrong number of dimensions - assert_raises(ValueError, restoration.denoise_tv_chambolle, - np.random.rand(8, 8, 8, 8)) + +def test_denoise_tv_chambolle_1d(): + """Apply the TV denoising algorithm on a 1D sinusoid.""" + x = 125 + 100*np.sin(np.linspace(0, 8*np.pi, 1000)) + x += 20 * np.random.rand(x.size) + x = np.clip(x, 0, 255) + res = restoration.denoise_tv_chambolle(x.astype(np.uint8), weight=0.1) + assert res.dtype == np.float + assert res.std() * 255 < x.std() + + +def test_denoise_tv_chambolle_4d(): + """ TV denoising for a 4D input.""" + im = 255 * np.random.rand(8, 8, 8, 8) + res = restoration.denoise_tv_chambolle(im.astype(np.uint8), weight=0.1) + assert res.dtype == np.float + assert res.std() * 255 < im.std() + + +def test_denoise_tv_chambolle_weighting(): + # make sure a specified weight gives consistent results regardless of + # the number of input image dimensions + rstate = np.random.RandomState(1234) + img2d = astro_gray.copy() + img2d += 0.15 * rstate.standard_normal(img2d.shape) + img2d = np.clip(img2d, 0, 1) + + # generate 4D image by tiling + img4d = np.tile(img2d[..., None, None], (1, 1, 2, 2)) + + w = 0.2 + denoised_2d = restoration.denoise_tv_chambolle(img2d, weight=w) + denoised_4d = restoration.denoise_tv_chambolle(img4d, weight=w) + assert measure.structural_similarity(denoised_2d, + denoised_4d[:, :, 0, 0]) > 0.99 def test_denoise_tv_bregman_2d():