From 2d326a464ed7c7cdf9b14a445333405dd7408958 Mon Sep 17 00:00:00 2001 From: "Gregory R. Lee" Date: Tue, 22 Dec 2015 16:36:14 -0500 Subject: [PATCH 01/11] ENH: n-dimensional refactor of TV denoising. _denoise_tv_chambolle_2d and _denoise_tv_chambolle_3d are replaced by _denoise_tv_chambolle_nd The restriction to 2D and 3D in denoies_tv_chambolle was removed. --- skimage/restoration/_denoise.py | 149 ++++++++------------------------ 1 file changed, 36 insertions(+), 113 deletions(-) diff --git a/skimage/restoration/_denoise.py b/skimage/restoration/_denoise.py index 97fab246..f86e4e0f 100644 --- a/skimage/restoration/_denoise.py +++ b/skimage/restoration/_denoise.py @@ -116,7 +116,7 @@ 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): +def _denoise_tv_chambolle_nd(im, weight=0.2, eps=2.e-4, n_iter_max=200): """Perform total-variation denoising on 3D images. Parameters @@ -146,112 +146,41 @@ 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.shape + (im.ndim, ), 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 = -p.sum(-1) + 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] = slice(0, -1) + slices_p[-1] = ax + d[slices_d] += p[slices_p] + slices_d[ax] = slice(None) + slices_p[ax] = 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) + slices_g = [slice(None), ] * (ndim + 1) + for ax in range(ndim): + slices_g[ax] = slice(0, -1) + slices_g[-1] = ax + g[slices_g] = np.diff(out, axis=ax) + slices_g[ax] = slice(None) + + norm = np.sqrt((g*g).sum(-1, keepdims=True)) E += weight * norm.sum() norm *= 0.5 / weight norm += 1. - px -= 1. / 6. * gx - px /= norm - py -= 1. / 6. * gy - py /= norm - pz -= 1 / 6. * gz - pz /= 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_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 + p -= 1. / (2.*ndim) * g + p /= norm E /= float(im.size) if i == 0: E_init = E @@ -271,7 +200,7 @@ def denoise_tv_chambolle(im, weight=0.2, eps=2.e-4, n_iter_max=200, 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 +218,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 +270,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[2]): + 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 From c2b86d8a02e63c1771d2961d53b1d46dc75380cf Mon Sep 17 00:00:00 2001 From: "Gregory R. Lee" Date: Tue, 22 Dec 2015 16:45:51 -0500 Subject: [PATCH 02/11] fix numpy 1.6 compatibility --- skimage/restoration/_denoise.py | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/skimage/restoration/_denoise.py b/skimage/restoration/_denoise.py index f86e4e0f..50eed95b 100644 --- a/skimage/restoration/_denoise.py +++ b/skimage/restoration/_denoise.py @@ -175,7 +175,12 @@ def _denoise_tv_chambolle_nd(im, weight=0.2, eps=2.e-4, n_iter_max=200): g[slices_g] = np.diff(out, axis=ax) slices_g[ax] = slice(None) - norm = np.sqrt((g*g).sum(-1, keepdims=True)) + try: + norm = np.sqrt((g*g).sum(-1, keepdims=True)) + except: + # no keepdims argument available prior to numpy 1.7 + norm = np.sqrt((g*g).sum(-1))[..., np.newaxis] + E += weight * norm.sum() norm *= 0.5 / weight norm += 1. From 4a0ad777806d9eaac4820640533a712e6697ac35 Mon Sep 17 00:00:00 2001 From: "Gregory R. Lee" Date: Tue, 22 Dec 2015 17:18:50 -0500 Subject: [PATCH 03/11] update tests for denoise_tv_chambolle --- skimage/restoration/tests/test_denoise.py | 20 +++++++++++++++++--- 1 file changed, 17 insertions(+), 3 deletions(-) diff --git a/skimage/restoration/tests/test_denoise.py b/skimage/restoration/tests/test_denoise.py index fb199afd..07d04b3c 100644 --- a/skimage/restoration/tests/test_denoise.py +++ b/skimage/restoration/tests/test_denoise.py @@ -66,9 +66,23 @@ def test_denoise_tv_chambolle_3d(): 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.5) + 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.5) + assert res.dtype == np.float + assert res.std() * 255 < im.std() def test_denoise_tv_bregman_2d(): From 4bc4f0a16903eef80c3fda11d4263764b3f4cd72 Mon Sep 17 00:00:00 2001 From: "Gregory R. Lee" Date: Tue, 22 Dec 2015 20:22:47 -0500 Subject: [PATCH 04/11] remove try/except in _denoise_chambolle_nd. just use the numpy 1.6 compatible case --- skimage/restoration/_denoise.py | 7 +------ 1 file changed, 1 insertion(+), 6 deletions(-) diff --git a/skimage/restoration/_denoise.py b/skimage/restoration/_denoise.py index 50eed95b..c4d64ee6 100644 --- a/skimage/restoration/_denoise.py +++ b/skimage/restoration/_denoise.py @@ -175,12 +175,7 @@ def _denoise_tv_chambolle_nd(im, weight=0.2, eps=2.e-4, n_iter_max=200): g[slices_g] = np.diff(out, axis=ax) slices_g[ax] = slice(None) - try: - norm = np.sqrt((g*g).sum(-1, keepdims=True)) - except: - # no keepdims argument available prior to numpy 1.7 - norm = np.sqrt((g*g).sum(-1))[..., np.newaxis] - + norm = np.sqrt((g*g).sum(-1))[..., np.newaxis] E += weight * norm.sum() norm *= 0.5 / weight norm += 1. From 4a8a2189f25bc5f626e877918052e6277a786a80 Mon Sep 17 00:00:00 2001 From: "Gregory R. Lee" Date: Wed, 23 Dec 2015 12:08:28 -0500 Subject: [PATCH 05/11] use Fortran order for g and p in _denoise_chambolle_nd to avoid a performance regression --- skimage/restoration/_denoise.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/skimage/restoration/_denoise.py b/skimage/restoration/_denoise.py index c4d64ee6..a676fd58 100644 --- a/skimage/restoration/_denoise.py +++ b/skimage/restoration/_denoise.py @@ -147,7 +147,7 @@ def _denoise_tv_chambolle_nd(im, weight=0.2, eps=2.e-4, n_iter_max=200): """ ndim = im.ndim - p = np.zeros(im.shape + (im.ndim, ), dtype=im.dtype) + p = np.zeros(im.shape + (im.ndim, ), dtype=im.dtype, order='F') g = np.zeros_like(p) d = np.zeros_like(im) i = 0 From 9f419024ad44b922e1299f532e4627616ca3b5bd Mon Sep 17 00:00:00 2001 From: "Gregory R. Lee" Date: Thu, 24 Dec 2015 09:35:34 -0500 Subject: [PATCH 06/11] Switch _denoise_tv_chambolle to use C-ordering --- skimage/restoration/_denoise.py | 24 ++++++++++++------------ 1 file changed, 12 insertions(+), 12 deletions(-) diff --git a/skimage/restoration/_denoise.py b/skimage/restoration/_denoise.py index a676fd58..aa1807fe 100644 --- a/skimage/restoration/_denoise.py +++ b/skimage/restoration/_denoise.py @@ -117,12 +117,12 @@ def denoise_tv_bregman(image, weight, max_iter=100, eps=1e-3, isotropic=True): def _denoise_tv_chambolle_nd(im, weight=0.2, eps=2.e-4, n_iter_max=200): - """Perform total-variation denoising on 3D images. + """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`). @@ -147,22 +147,22 @@ def _denoise_tv_chambolle_nd(im, weight=0.2, eps=2.e-4, n_iter_max=200): """ ndim = im.ndim - p = np.zeros(im.shape + (im.ndim, ), dtype=im.dtype, order='F') + 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: if i > 0: - d = -p.sum(-1) + 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] = slice(0, -1) - slices_p[-1] = ax + 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] = slice(None) + slices_p[ax+1] = slice(None) out = im + d else: out = im @@ -170,12 +170,12 @@ def _denoise_tv_chambolle_nd(im, weight=0.2, eps=2.e-4, n_iter_max=200): slices_g = [slice(None), ] * (ndim + 1) for ax in range(ndim): - slices_g[ax] = slice(0, -1) - slices_g[-1] = ax + slices_g[ax+1] = slice(0, -1) + slices_g[0] = ax g[slices_g] = np.diff(out, axis=ax) - slices_g[ax] = slice(None) + slices_g[ax+1] = slice(None) - norm = np.sqrt((g*g).sum(-1))[..., np.newaxis] + norm = np.sqrt((g ** 2).sum(0))[np.newaxis, ...] E += weight * norm.sum() norm *= 0.5 / weight norm += 1. @@ -196,7 +196,7 @@ def _denoise_tv_chambolle_nd(im, weight=0.2, eps=2.e-4, n_iter_max=200): def denoise_tv_chambolle(im, weight=0.2, 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 ---------- From dd5a708a33a25f0903c11cf3873f37873393af53 Mon Sep 17 00:00:00 2001 From: "Gregory R. Lee" Date: Thu, 24 Dec 2015 12:28:41 -0500 Subject: [PATCH 07/11] minor bugfix so that weight gives a consistent amount of smoothing regardless of the number of dimensions. As in the original reference, the stepsize should be in both the numerator and denominator of p. --- skimage/restoration/_denoise.py | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/skimage/restoration/_denoise.py b/skimage/restoration/_denoise.py index aa1807fe..1fb3d89a 100644 --- a/skimage/restoration/_denoise.py +++ b/skimage/restoration/_denoise.py @@ -153,6 +153,7 @@ def _denoise_tv_chambolle_nd(im, weight=0.2, eps=2.e-4, n_iter_max=200): i = 0 while i < n_iter_max: 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) @@ -168,6 +169,8 @@ def _denoise_tv_chambolle_nd(im, weight=0.2, eps=2.e-4, n_iter_max=200): out = im E = (d ** 2).sum() + # 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) @@ -177,9 +180,10 @@ def _denoise_tv_chambolle_nd(im, weight=0.2, eps=2.e-4, n_iter_max=200): norm = np.sqrt((g ** 2).sum(0))[np.newaxis, ...] E += weight * norm.sum() - norm *= 0.5 / weight + tau = 1. / (2.*ndim) + norm *= tau / weight norm += 1. - p -= 1. / (2.*ndim) * g + p -= tau * g p /= norm E /= float(im.size) if i == 0: From e6a1a1933704f20f0fa41659dfad1f86385a39d1 Mon Sep 17 00:00:00 2001 From: "Gregory R. Lee" Date: Thu, 24 Dec 2015 12:34:11 -0500 Subject: [PATCH 08/11] update default weight to reflect the bugfix from hardcoded 0.5 to 1/(2*ndim) use the default weight in all the tests as well. --- skimage/restoration/_denoise.py | 4 ++-- skimage/restoration/tests/test_denoise.py | 14 +++++++------- 2 files changed, 9 insertions(+), 9 deletions(-) diff --git a/skimage/restoration/_denoise.py b/skimage/restoration/_denoise.py index 1fb3d89a..9e06fbd8 100644 --- a/skimage/restoration/_denoise.py +++ b/skimage/restoration/_denoise.py @@ -116,7 +116,7 @@ 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_nd(im, weight=0.2, eps=2.e-4, n_iter_max=200): +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 @@ -198,7 +198,7 @@ def _denoise_tv_chambolle_nd(im, weight=0.2, eps=2.e-4, n_iter_max=200): 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 n-dimensional images. diff --git a/skimage/restoration/tests/test_denoise.py b/skimage/restoration/tests/test_denoise.py index 07d04b3c..9aa91526 100644 --- a/skimage/restoration/tests/test_denoise.py +++ b/skimage/restoration/tests/test_denoise.py @@ -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,8 @@ 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) @@ -46,7 +46,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,7 +62,7 @@ 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() @@ -72,7 +72,7 @@ def test_denoise_tv_chambolle_1d(): 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.5) + res = restoration.denoise_tv_chambolle(x.astype(np.uint8), weight=0.1) assert res.dtype == np.float assert res.std() * 255 < x.std() @@ -80,7 +80,7 @@ def test_denoise_tv_chambolle_1d(): 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.5) + res = restoration.denoise_tv_chambolle(im.astype(np.uint8), weight=0.1) assert res.dtype == np.float assert res.std() * 255 < im.std() From a30592a22aa693723b1d0efe74f7baceba7abbe9 Mon Sep 17 00:00:00 2001 From: "Gregory R. Lee" Date: Sat, 23 Jan 2016 11:23:21 -0500 Subject: [PATCH 09/11] MAINT: explicitly specify axis argument --- skimage/restoration/_denoise.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/skimage/restoration/_denoise.py b/skimage/restoration/_denoise.py index 9e06fbd8..acd25cb0 100644 --- a/skimage/restoration/_denoise.py +++ b/skimage/restoration/_denoise.py @@ -178,7 +178,7 @@ def _denoise_tv_chambolle_nd(im, weight=0.1, eps=2.e-4, n_iter_max=200): g[slices_g] = np.diff(out, axis=ax) slices_g[ax+1] = slice(None) - norm = np.sqrt((g ** 2).sum(0))[np.newaxis, ...] + norm = np.sqrt((g ** 2).sum(axis=0))[np.newaxis, ...] E += weight * norm.sum() tau = 1. / (2.*ndim) norm *= tau / weight From 53742be9c632c668bae171d080bf4c75007aa0f9 Mon Sep 17 00:00:00 2001 From: "Gregory R. Lee" Date: Sun, 24 Jan 2016 02:24:37 -0500 Subject: [PATCH 10/11] fix axis bug in multichannel case and add a corresponding test --- skimage/restoration/_denoise.py | 2 +- skimage/restoration/tests/test_denoise.py | 9 +++++++++ 2 files changed, 10 insertions(+), 1 deletion(-) diff --git a/skimage/restoration/_denoise.py b/skimage/restoration/_denoise.py index acd25cb0..62985cc7 100644 --- a/skimage/restoration/_denoise.py +++ b/skimage/restoration/_denoise.py @@ -276,7 +276,7 @@ def denoise_tv_chambolle(im, weight=0.1, eps=2.e-4, n_iter_max=200, if multichannel: out = np.zeros_like(im) - for c in range(im.shape[2]): + for c in range(im.shape[-1]): out[..., c] = _denoise_tv_chambolle_nd(im[..., c], weight, eps, n_iter_max) else: diff --git a/skimage/restoration/tests/test_denoise.py b/skimage/restoration/tests/test_denoise.py index 9aa91526..5f4ce178 100644 --- a/skimage/restoration/tests/test_denoise.py +++ b/skimage/restoration/tests/test_denoise.py @@ -39,6 +39,15 @@ def test_denoise_tv_chambolle_multichannel(): 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) + def test_denoise_tv_chambolle_float_result_range(): # astronaut image From 682f0895e7f045ddbdebc1450f7660966212dd4f Mon Sep 17 00:00:00 2001 From: "Gregory R. Lee" Date: Sun, 24 Jan 2016 12:27:10 -0500 Subject: [PATCH 11/11] TST: add test for constent TV weighting across number of image dimensions for Chambolle algorithm --- skimage/restoration/tests/test_denoise.py | 20 +++++++++++++++++++- 1 file changed, 19 insertions(+), 1 deletion(-) diff --git a/skimage/restoration/tests/test_denoise.py b/skimage/restoration/tests/test_denoise.py index 5f4ce178..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) @@ -94,6 +94,24 @@ def test_denoise_tv_chambolle_4d(): 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(): img = checkerboard_gray.copy() # add some random noise