From 24b49fc8ee7693f2a6c2f3648a3dd66ab39c54c5 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= <hannesschoenberger@gmail.com>
Date: Fri, 2 Nov 2012 20:10:45 +0100
Subject: [PATCH 01/16] Improve total-varation denoising algorithm

Implementation of fast split-Bregman optimization algorithm in Cython. This
implementation also fixes the previously broken 3D version, which darkened
the images.
---
 skimage/filter/__init__.py  |   4 +-
 skimage/filter/_denoise.pyx | 169 +++++++++++++++++++++++++++++++++++-
 2 files changed, 170 insertions(+), 3 deletions(-)

diff --git a/skimage/filter/__init__.py b/skimage/filter/__init__.py
index f1c1fd49..8894ff1a 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 tv_denoise
+from ._denoise import denoise_bilateral, denoise_tv
 from ._rank_order import rank_order
 from .thresholding import threshold_otsu, threshold_adaptive
diff --git a/skimage/filter/_denoise.pyx b/skimage/filter/_denoise.pyx
index b60a5a85..ca75ab4f 100644
--- a/skimage/filter/_denoise.pyx
+++ b/skimage/filter/_denoise.pyx
@@ -7,6 +7,7 @@ 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
 
@@ -174,4 +175,170 @@ def denoise_bilateral(image, int win_size=5, sigma_range=None,
     free(centres)
     free(total_values)
 
-    return out
+    return np.squeeze(out)
+
+
+cdef inline double _get_elem(double* image, Py_ssize_t rows, Py_ssize_t cols,
+                             Py_ssize_t dims, Py_ssize_t r, Py_ssize_t c,
+                             Py_ssize_t k):
+    return image[r * cols * dims + c * dims + k]
+
+
+cdef inline void _set_elem(double* image, Py_ssize_t rows, Py_ssize_t cols,
+                             Py_ssize_t dims, Py_ssize_t r, Py_ssize_t c,
+                             Py_ssize_t k, double value):
+    image[r * cols * dims + c * dims + k] = value
+
+
+cdef inline void _incr_elem(double* image, Py_ssize_t rows, Py_ssize_t cols,
+                              Py_ssize_t dims, Py_ssize_t r, Py_ssize_t c,
+                              Py_ssize_t k, double value):
+    image[r * cols * dims + c * dims + k] += value
+
+
+def denoise_tv(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 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] ftp://ftp.math.ucla.edu/pub/camreport/cam08-29.pdf
+    .. [3] 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* image_data = <double*>cimage.data
+        double* u_data = <double*>u.data
+
+        double* dx_data = <double*>dx.data
+        double* dy_data = <double*>dy.data
+        double* bx_data = <double*>bx.data
+        double* by_data = <double*>by.data
+
+        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 = _get_elem(u_data, rows2, cols2, dims, r, c, k)
+
+                    # forward derivatives
+                    ux = _get_elem(u_data, rows2, cols2, dims,
+                                   r, c+1, k) - uprev
+                    uy = _get_elem(u_data, rows2, cols2, dims,
+                                   r+1, c, k) - uprev
+
+                    # Gauss-Seidel method
+                    unew = (
+                        lam * (
+                            + _get_elem(u_data, rows2, cols2, dims, r+1, c, k)
+                            + _get_elem(u_data, rows2, cols2, dims, r-1, c, k)
+                            + _get_elem(u_data, rows2, cols2, dims, r, c+1, k)
+                            + _get_elem(u_data, rows2, cols2, dims, r, c-1, k)
+
+                            + _get_elem(dx_data, rows2, cols2, dims, r, c-1, k)
+                            - _get_elem(dx_data, rows2, cols2, dims, r, c, k)
+                            + _get_elem(dy_data, rows2, cols2, dims, r-1, c, k)
+                            - _get_elem(dy_data, rows2, cols2, dims, r, c, k)
+
+                            - _get_elem(bx_data, rows2, cols2, dims, r, c-1, k)
+                            + _get_elem(bx_data, rows2, cols2, dims, r, c, k)
+                            - _get_elem(by_data, rows2, cols2, dims, r-1, c, k)
+                            + _get_elem(by_data, rows2, cols2, dims, r, c, k)
+                        ) + weight * _get_elem(image_data, rows, cols, dims,
+                                               r-1, c-1, k)
+                    ) / norm
+                    _set_elem(u_data, rows2, cols2, dims, r, c, k, unew)
+
+                    # update root mean square error
+                    rmse += (unew - uprev)**2
+
+                    bxx = _get_elem(bx_data, rows2, cols2, dims, r, c, k)
+                    byy = _get_elem(by_data, rows2, cols2, dims, 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)
+
+                    _set_elem(dx_data, rows2, cols2, dims, r, c, k, dxx)
+                    _set_elem(dy_data, rows2, cols2, dims, r, c, k, dyy)
+
+                    _incr_elem(bx_data, rows2, cols2, dims, r, c, k, ux - dxx)
+                    _incr_elem(by_data, rows2, cols2, dims, r, c, k, uy - dyy)
+
+        rmse = sqrt(rmse / total)
+        i += 1
+
+    return np.squeeze(u[1:-1, 1:-1])

From cda03cfba4920a3c526424614305ae7e6a421e82 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= <hannesschoenberger@gmail.com>
Date: Fri, 2 Nov 2012 20:17:48 +0100
Subject: [PATCH 02/16] Fix test cases for TV denoising

---
 skimage/filter/tests/test_denoise.py | 46 +++++++++++-----------------
 1 file changed, 18 insertions(+), 28 deletions(-)

diff --git a/skimage/filter/tests/test_denoise.py b/skimage/filter/tests/test_denoise.py
index 63a5a5e0..03c8c58a 100644
--- a/skimage/filter/tests/test_denoise.py
+++ b/skimage/filter/tests/test_denoise.py
@@ -9,24 +9,17 @@ lena_gray = color.rgb2gray(lena)
 
 
 def test_denoise_tv_2d():
-    # lena image
     img = lena_gray
-    # add noise to lena
+    # add some random noise
     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)
+
+    out1 = filter.denoise_tv(img, weight=10)
+    out2 = filter.denoise_tv(img, weight=5)
+
+    # make sure noise is reduced
+    assert img.std() > out1.std()
+    assert out1.std() > out2.std()
 
 
 def test_denoise_tv_float_result_range():
@@ -42,20 +35,17 @@ def test_denoise_tv_float_result_range():
 
 
 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()
+    img = lena
+    # add some random noise
+    img += 0.5 * img.std() * np.random.random(img.shape)
+    img = np.clip(img, 0, 1)
 
-    # test wrong number of dimensions
-    assert_raises(ValueError, filter.denoise_tv, np.random.random((8, 8, 8, 8)))
+    out1 = filter.denoise_tv(img, weight=10)
+    out2 = filter.denoise_tv(img, weight=5)
+
+    # make sure noise is reduced
+    assert img.std() > out1.std()
+    assert out1.std() > out2.std()
 
 
 def test_denoise_bilateral_2d():

From dd34097fe4d65b24a83fafa7c8348672a688f474 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= <hannesschoenberger@gmail.com>
Date: Fri, 2 Nov 2012 20:18:28 +0100
Subject: [PATCH 03/16] Remove old implementation if TV filter

---
 skimage/filter/__init__.py  |   3 +-
 skimage/filter/_denoise.pyx |   3 +
 skimage/filter/denoise.py   | 245 ------------------------------------
 3 files changed, 4 insertions(+), 247 deletions(-)
 delete mode 100644 skimage/filter/denoise.py

diff --git a/skimage/filter/__init__.py b/skimage/filter/__init__.py
index 8894ff1a..4f8b129a 100644
--- a/skimage/filter/__init__.py
+++ b/skimage/filter/__init__.py
@@ -3,7 +3,6 @@ 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
-from ._denoise import denoise_bilateral, denoise_tv
+from ._denoise import denoise_bilateral, denoise_tv, tv_denoise
 from ._rank_order import rank_order
 from .thresholding import threshold_otsu, threshold_adaptive
diff --git a/skimage/filter/_denoise.pyx b/skimage/filter/_denoise.pyx
index ca75ab4f..2e519b8f 100644
--- a/skimage/filter/_denoise.pyx
+++ b/skimage/filter/_denoise.pyx
@@ -10,6 +10,7 @@ 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):
@@ -342,3 +343,5 @@ def denoise_tv(image, double weight, int max_iter=100, double eps=1e-3):
         i += 1
 
     return np.squeeze(u[1:-1, 1:-1])
+
+tv_denoise = deprecated('skimage.filter.denoise_tv')(denoise_tv)
diff --git a/skimage/filter/denoise.py b/skimage/filter/denoise.py
deleted file mode 100644
index 351b2482..00000000
--- a/skimage/filter/denoise.py
+++ /dev/null
@@ -1,245 +0,0 @@
-import numpy as np
-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.
-
-    Parameters
-    ----------
-    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
-        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
-    -----
-    Rudin, Osher and Fatemi algorithm.
-
-    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 = mask.astype(np.float)
-    >>> mask += 0.2*np.random.randn(*mask.shape)
-    >>> res = denoise_tv_3d(mask, weight=100)
-    """
-    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)
-    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
-        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)
-        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_2d(im, weight=50, eps=2.e-4, n_iter_max=200):
-    """Perform total-variation denoising.
-
-    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.
-
-    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)
-    """
-    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(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,
-        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
-        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
-
-    The principle of total variation denoising is to minimize the
-    total variation of the image, which can be roughly described as
-    the integral of the norm of the image gradient. Total variation
-    denoising tends to produce "cartoon-like" images, that is,
-    piecewise-constant images.
-
-    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.
-
-    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)
-    >>> 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 = 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 = _denoise_tv_2d(im, weight, eps, n_iter_max)
-    elif im.ndim == 3:
-        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)

From ccd2f4264308979b56867388e7b5cc7a1beea26f Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= <hannesschoenberger@gmail.com>
Date: Fri, 2 Nov 2012 20:18:49 +0100
Subject: [PATCH 04/16] Update denoising example for new TV implementation

---
 doc/examples/plot_denoise.py | 8 +++++---
 1 file changed, 5 insertions(+), 3 deletions(-)

diff --git a/doc/examples/plot_denoise.py b/doc/examples/plot_denoise.py
index 8f02f21b..cfcec5e4 100644
--- a/doc/examples/plot_denoise.py
+++ b/doc/examples/plot_denoise.py
@@ -32,7 +32,7 @@ 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 = color.rgb2gray(img_as_float(data.lena()))
 lena = lena[220:300, 220:320]
 
 noisy = lena + 0.5 * lena.std() * np.random.random(lena.shape)
@@ -40,17 +40,19 @@ 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(noisy, weight=40))
 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].imshow(denoise_tv(noisy, weight=20))
 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))

From fff33702eaf1d12f14f967d4c3b15164664a99a3 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= <hannesschoenberger@gmail.com>
Date: Sun, 4 Nov 2012 13:39:24 +0100
Subject: [PATCH 05/16] Replace manual array with builting Cython buffer
 indexing

---
 skimage/filter/_denoise.pyx | 75 ++++++++++++-------------------------
 1 file changed, 23 insertions(+), 52 deletions(-)

diff --git a/skimage/filter/_denoise.pyx b/skimage/filter/_denoise.pyx
index 2e519b8f..436097f6 100644
--- a/skimage/filter/_denoise.pyx
+++ b/skimage/filter/_denoise.pyx
@@ -179,24 +179,6 @@ def denoise_bilateral(image, int win_size=5, sigma_range=None,
     return np.squeeze(out)
 
 
-cdef inline double _get_elem(double* image, Py_ssize_t rows, Py_ssize_t cols,
-                             Py_ssize_t dims, Py_ssize_t r, Py_ssize_t c,
-                             Py_ssize_t k):
-    return image[r * cols * dims + c * dims + k]
-
-
-cdef inline void _set_elem(double* image, Py_ssize_t rows, Py_ssize_t cols,
-                             Py_ssize_t dims, Py_ssize_t r, Py_ssize_t c,
-                             Py_ssize_t k, double value):
-    image[r * cols * dims + c * dims + k] = value
-
-
-cdef inline void _incr_elem(double* image, Py_ssize_t rows, Py_ssize_t cols,
-                              Py_ssize_t dims, Py_ssize_t r, Py_ssize_t c,
-                              Py_ssize_t k, double value):
-    image[r * cols * dims + c * dims + k] += value
-
-
 def denoise_tv(image, double weight, int max_iter=100, double eps=1e-3):
     """Perform total-variation denoising using split-Bregman optimization.
 
@@ -263,14 +245,6 @@ def denoise_tv(image, double weight, int max_iter=100, double eps=1e-3):
         cnp.ndarray[dtype=cnp.double_t, ndim=3, mode='c'] by = \
             np.zeros(shape_ext, dtype=np.double)
 
-        double* image_data = <double*>cimage.data
-        double* u_data = <double*>u.data
-
-        double* dx_data = <double*>dx.data
-        double* dy_data = <double*>dy.data
-        double* bx_data = <double*>bx.data
-        double* by_data = <double*>by.data
-
         double ux, uy, uprev, unew, bxx, byy, dxx, dyy, s
         int i = 0
         double lam = 2 * weight
@@ -293,51 +267,48 @@ def denoise_tv(image, double weight, int max_iter=100, double eps=1e-3):
             for r in range(1, rows + 1):
                 for c in range(1, cols + 1):
 
-                    uprev = _get_elem(u_data, rows2, cols2, dims, r, c, k)
+                    uprev = u[r, c, k]
 
                     # forward derivatives
-                    ux = _get_elem(u_data, rows2, cols2, dims,
-                                   r, c+1, k) - uprev
-                    uy = _get_elem(u_data, rows2, cols2, dims,
-                                   r+1, c, k) - uprev
+                    ux = u[r, c + 1, k] - uprev
+                    uy = u[r + 1, c, k] - uprev
 
                     # Gauss-Seidel method
                     unew = (
                         lam * (
-                            + _get_elem(u_data, rows2, cols2, dims, r+1, c, k)
-                            + _get_elem(u_data, rows2, cols2, dims, r-1, c, k)
-                            + _get_elem(u_data, rows2, cols2, dims, r, c+1, k)
-                            + _get_elem(u_data, rows2, cols2, dims, r, c-1, k)
+                            + u[r + 1, c, k]
+                            + u[r - 1, c, k]
+                            + u[r, c + 1, k]
+                            + u[r, c - 1, k]
 
-                            + _get_elem(dx_data, rows2, cols2, dims, r, c-1, k)
-                            - _get_elem(dx_data, rows2, cols2, dims, r, c, k)
-                            + _get_elem(dy_data, rows2, cols2, dims, r-1, c, k)
-                            - _get_elem(dy_data, rows2, cols2, dims, r, c, k)
+                            + dx[r, c - 1, k]
+                            - dx[r, c, k]
+                            + dy[r - 1, c, k]
+                            - dy[r, c, k]
 
-                            - _get_elem(bx_data, rows2, cols2, dims, r, c-1, k)
-                            + _get_elem(bx_data, rows2, cols2, dims, r, c, k)
-                            - _get_elem(by_data, rows2, cols2, dims, r-1, c, k)
-                            + _get_elem(by_data, rows2, cols2, dims, r, c, k)
-                        ) + weight * _get_elem(image_data, rows, cols, dims,
-                                               r-1, c-1, 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
-                    _set_elem(u_data, rows2, cols2, dims, r, c, k, unew)
+                    u[r, c, k] = unew
 
                     # update root mean square error
                     rmse += (unew - uprev)**2
 
-                    bxx = _get_elem(bx_data, rows2, cols2, dims, r, c, k)
-                    byy = _get_elem(by_data, rows2, cols2, dims, r, c, k)
+                    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)
 
-                    _set_elem(dx_data, rows2, cols2, dims, r, c, k, dxx)
-                    _set_elem(dy_data, rows2, cols2, dims, r, c, k, dyy)
+                    dx[r, c, k] = dxx
+                    dy[r, c, k] = dyy
 
-                    _incr_elem(bx_data, rows2, cols2, dims, r, c, k, ux - dxx)
-                    _incr_elem(by_data, rows2, cols2, dims, r, c, k, uy - dyy)
+                    bx[r, c, k] += ux - dxx
+                    by[r, c, k] += uy - dyy
 
         rmse = sqrt(rmse / total)
         i += 1

From 780d886db95f45ff17c2535f30bd1ec026357792 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= <hannesschoenberger@gmail.com>
Date: Tue, 25 Dec 2012 10:44:19 +0100
Subject: [PATCH 06/16] Rename denoise_tv to denoise_tv_bregman

---
 skimage/filter/__init__.py           |  2 +-
 skimage/filter/_denoise.pyx          |  4 +---
 skimage/filter/tests/test_denoise.py | 10 +++++-----
 3 files changed, 7 insertions(+), 9 deletions(-)

diff --git a/skimage/filter/__init__.py b/skimage/filter/__init__.py
index 4f8b129a..39d50376 100644
--- a/skimage/filter/__init__.py
+++ b/skimage/filter/__init__.py
@@ -3,6 +3,6 @@ from .ctmf import median_filter
 from ._canny import canny
 from .edges import (sobel, hsobel, vsobel, scharr, hscharr, vscharr, prewitt,
                     hprewitt, vprewitt)
-from ._denoise import denoise_bilateral, denoise_tv, tv_denoise
+from ._denoise 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.pyx b/skimage/filter/_denoise.pyx
index 436097f6..ea046a0a 100644
--- a/skimage/filter/_denoise.pyx
+++ b/skimage/filter/_denoise.pyx
@@ -179,7 +179,7 @@ def denoise_bilateral(image, int win_size=5, sigma_range=None,
     return np.squeeze(out)
 
 
-def denoise_tv(image, double weight, int max_iter=100, double eps=1e-3):
+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)
@@ -314,5 +314,3 @@ def denoise_tv(image, double weight, int max_iter=100, double eps=1e-3):
         i += 1
 
     return np.squeeze(u[1:-1, 1:-1])
-
-tv_denoise = deprecated('skimage.filter.denoise_tv')(denoise_tv)
diff --git a/skimage/filter/tests/test_denoise.py b/skimage/filter/tests/test_denoise.py
index 03c8c58a..6f77fdc7 100644
--- a/skimage/filter/tests/test_denoise.py
+++ b/skimage/filter/tests/test_denoise.py
@@ -14,8 +14,8 @@ def test_denoise_tv_2d():
     img += 0.5 * img.std() * np.random.random(img.shape)
     img = np.clip(img, 0, 1)
 
-    out1 = filter.denoise_tv(img, weight=10)
-    out2 = filter.denoise_tv(img, weight=5)
+    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()
@@ -27,7 +27,7 @@ def test_denoise_tv_float_result_range():
     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_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
@@ -40,8 +40,8 @@ def test_denoise_tv_3d():
     img += 0.5 * img.std() * np.random.random(img.shape)
     img = np.clip(img, 0, 1)
 
-    out1 = filter.denoise_tv(img, weight=10)
-    out2 = filter.denoise_tv(img, weight=5)
+    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()

From f2840ba9067872f7f217d145d52ef91ed9ba1576 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= <hannesschoenberger@gmail.com>
Date: Tue, 25 Dec 2012 10:55:46 +0100
Subject: [PATCH 07/16] Rename denoising cython source file

---
 skimage/filter/__init__.py                       | 2 +-
 skimage/filter/{_denoise.pyx => _denoise_cy.pyx} | 0
 skimage/filter/setup.py                          | 4 ++--
 3 files changed, 3 insertions(+), 3 deletions(-)
 rename skimage/filter/{_denoise.pyx => _denoise_cy.pyx} (100%)

diff --git a/skimage/filter/__init__.py b/skimage/filter/__init__.py
index 39d50376..80d637cf 100644
--- a/skimage/filter/__init__.py
+++ b/skimage/filter/__init__.py
@@ -3,6 +3,6 @@ from .ctmf import median_filter
 from ._canny import canny
 from .edges import (sobel, hsobel, vsobel, scharr, hscharr, vscharr, prewitt,
                     hprewitt, vprewitt)
-from ._denoise import denoise_bilateral, denoise_tv_bregman
+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.pyx b/skimage/filter/_denoise_cy.pyx
similarity index 100%
rename from skimage/filter/_denoise.pyx
rename to skimage/filter/_denoise_cy.pyx
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()])

From dfd8cb4fc203594df6ed14dede882381041e5930 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= <hannesschoenberger@gmail.com>
Date: Tue, 25 Dec 2012 11:00:38 +0100
Subject: [PATCH 08/16] Fix typo in denoise_tv_bregman doc string

---
 skimage/filter/_denoise_cy.pyx | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/skimage/filter/_denoise_cy.pyx b/skimage/filter/_denoise_cy.pyx
index ea046a0a..47e2a6d5 100644
--- a/skimage/filter/_denoise_cy.pyx
+++ b/skimage/filter/_denoise_cy.pyx
@@ -183,7 +183,7 @@ 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 total-variation under the constraint
+    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.
 

From ae5488a2b189af807d008bf4e23ccba167ac8907 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= <hannesschoenberger@gmail.com>
Date: Tue, 25 Dec 2012 11:04:39 +0100
Subject: [PATCH 09/16] Rename test functions for bregman algorithm

---
 skimage/filter/tests/test_denoise.py | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/skimage/filter/tests/test_denoise.py b/skimage/filter/tests/test_denoise.py
index 6f77fdc7..10745b72 100644
--- a/skimage/filter/tests/test_denoise.py
+++ b/skimage/filter/tests/test_denoise.py
@@ -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_bregman_2d():
     img = lena_gray
     # add some random noise
     img += 0.5 * img.std() * np.random.random(img.shape)
@@ -22,7 +22,7 @@ def test_denoise_tv_2d():
     assert out1.std() > out2.std()
 
 
-def test_denoise_tv_float_result_range():
+def test_denoise_tv_bregman_float_result_range():
     # lena image
     img = lena_gray
     int_lena = np.multiply(img, 255).astype(np.uint8)
@@ -34,7 +34,7 @@ def test_denoise_tv_float_result_range():
     assert np.min(denoised_int_lena) >= 0.0
 
 
-def test_denoise_tv_3d():
+def test_denoise_tv_bregman_3d():
     img = lena
     # add some random noise
     img += 0.5 * img.std() * np.random.random(img.shape)

From 236ef48ec97fe389aeaf5bfeb195863d59dd928f Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= <hannesschoenberger@gmail.com>
Date: Tue, 25 Dec 2012 11:15:41 +0100
Subject: [PATCH 10/16] Add original (Chambolle) implementation of TV denoising

---
 skimage/filter/__init__.py |   1 +
 skimage/filter/_denoise.py | 261 +++++++++++++++++++++++++++++++++++++
 2 files changed, 262 insertions(+)
 create mode 100644 skimage/filter/_denoise.py

diff --git a/skimage/filter/__init__.py b/skimage/filter/__init__.py
index 80d637cf..f32c3b23 100644
--- a/skimage/filter/__init__.py
+++ b/skimage/filter/__init__.py
@@ -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 ._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
new file mode 100644
index 00000000..315efc05
--- /dev/null
+++ b/skimage/filter/_denoise.py
@@ -0,0 +1,261 @@
+import numpy as np
+from skimage import img_as_float
+from skimage._shared.utils import deprecated
+
+
+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
+        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
+        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
+    -----
+    Rudin, Osher and Fatemi algorithm.
+
+    Examples
+    ---------
+    >>> 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 = denoise_tv(mask, weight=100)
+
+    """
+
+    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)
+    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
+        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)
+        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=50, 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.
+
+    Examples
+    ---------
+    >>> 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)
+    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=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
+        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
+        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.
+    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 image.
+
+    Notes
+    -----
+    The principle of total variation denoising is explained in
+    http://en.wikipedia.org/wiki/Total_variation_denoising
+
+    The principle of total variation denoising is to minimize the
+    total variation of the image, which can be roughly described as
+    the integral of the norm of the image gradient. Total variation
+    denoising tends to produce "cartoon-like" images, that is,
+    piecewise-constant images.
+
+    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.
+
+    Examples
+    ---------
+    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:
+
+    >>> 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 = 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)
+    elif im.ndim == 3:
+        if multichannel:
+            out = np.zeros_like(im)
+            for c in range(im.shape[2]):
+                out[..., c] = _denoise_tv_2d(im[..., c], weight, eps,
+                                             n_iter_max)
+        else:
+            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_chambolle')\
+                       (denoise_tv_chambolle)

From b7254f225ce10c371310363f5af611d3d3ac0f48 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= <hannesschoenberger@gmail.com>
Date: Tue, 25 Dec 2012 11:26:43 +0100
Subject: [PATCH 11/16] Update denoising example

---
 doc/examples/plot_denoise.py | 14 +++++++-------
 1 file changed, 7 insertions(+), 7 deletions(-)

diff --git a/doc/examples/plot_denoise.py b/doc/examples/plot_denoise.py
index cfcec5e4..debb2ea6 100644
--- a/doc/examples/plot_denoise.py
+++ b/doc/examples/plot_denoise.py
@@ -30,12 +30,12 @@ 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 = color.rgb2gray(img_as_float(data.lena()))
+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))
@@ -45,17 +45,17 @@ 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=40))
+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=20))
+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)

From aadc45cd5b42f5d15b6fb041197246a77bf2dd7a Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= <hannesschoenberger@gmail.com>
Date: Tue, 25 Dec 2012 11:27:08 +0100
Subject: [PATCH 12/16] Fix bugs caused by renaming functions

---
 skimage/filter/_denoise.py | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/skimage/filter/_denoise.py b/skimage/filter/_denoise.py
index 315efc05..07f3c58f 100644
--- a/skimage/filter/_denoise.py
+++ b/skimage/filter/_denoise.py
@@ -242,15 +242,15 @@ def denoise_tv_chambolle(im, weight=50, eps=2.e-4, n_iter_max=200,
         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:
         if multichannel:
             out = np.zeros_like(im)
             for c in range(im.shape[2]):
-                out[..., c] = _denoise_tv_2d(im[..., c], weight, eps,
+                out[..., c] = _denoise_tv_chambolle_2d(im[..., c], weight, eps,
                                              n_iter_max)
         else:
-            out = _denoise_tv_3d(im, weight, eps, n_iter_max)
+            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')

From 5078da0aed06e4b512a0bc880ad407745e128d35 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= <hannesschoenberger@gmail.com>
Date: Tue, 25 Dec 2012 11:30:21 +0100
Subject: [PATCH 13/16] Add test cases for chambolle tv denoising
 implementation

---
 skimage/filter/tests/test_denoise.py | 51 ++++++++++++++++++++++++++++
 1 file changed, 51 insertions(+)

diff --git a/skimage/filter/tests/test_denoise.py b/skimage/filter/tests/test_denoise.py
index 10745b72..c8e7bb02 100644
--- a/skimage/filter/tests/test_denoise.py
+++ b/skimage/filter/tests/test_denoise.py
@@ -8,6 +8,57 @@ lena = img_as_float(data.lena()[:256, :256])
 lena_gray = color.rgb2gray(lena)
 
 
+def test_denoise_tv_chambolle_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_chambolle(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_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_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_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
+    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_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_chambolle,
+                  np.random.random((8, 8, 8, 8)))
+
+
 def test_denoise_tv_bregman_2d():
     img = lena_gray
     # add some random noise

From db5460591abecf55bd09e2d76721b66efd213557 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= <hannesschoenberger@gmail.com>
Date: Tue, 25 Dec 2012 11:33:16 +0100
Subject: [PATCH 14/16] Add test case for multi-channel chambolle denoising

---
 skimage/filter/tests/test_denoise.py | 8 +++++++-
 1 file changed, 7 insertions(+), 1 deletion(-)

diff --git a/skimage/filter/tests/test_denoise.py b/skimage/filter/tests/test_denoise.py
index c8e7bb02..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
 
@@ -29,6 +29,12 @@ def test_denoise_tv_chambolle_2d():
             < np.sqrt((grad**2).sum()) / 2)
 
 
+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

From d82ae4dd95b3184cb09617c04969d7e5483ffdc2 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= <hannesschoenberger@gmail.com>
Date: Tue, 25 Dec 2012 11:36:48 +0100
Subject: [PATCH 15/16] Add note about multichannel parameter

---
 skimage/filter/_denoise.py | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/skimage/filter/_denoise.py b/skimage/filter/_denoise.py
index 07f3c58f..e26f8387 100644
--- a/skimage/filter/_denoise.py
+++ b/skimage/filter/_denoise.py
@@ -200,6 +200,8 @@ def denoise_tv_chambolle(im, weight=50, eps=2.e-4, n_iter_max=200,
 
     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
 

From 81d88b8cc3a4a1b35ad3e6c29cf02ef95cd94331 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= <hannesschoenberger@gmail.com>
Date: Thu, 27 Dec 2012 08:44:50 +0100
Subject: [PATCH 16/16] Include more information about reference papers

---
 skimage/filter/_denoise_cy.pyx | 8 ++++++--
 1 file changed, 6 insertions(+), 2 deletions(-)

diff --git a/skimage/filter/_denoise_cy.pyx b/skimage/filter/_denoise_cy.pyx
index 47e2a6d5..8ebf5c73 100644
--- a/skimage/filter/_denoise_cy.pyx
+++ b/skimage/filter/_denoise_cy.pyx
@@ -212,8 +212,12 @@ def denoise_tv_bregman(image, double weight, int max_iter=100, double eps=1e-3):
     References
     ----------
     .. [1] http://en.wikipedia.org/wiki/Total_variation_denoising
-    .. [2] ftp://ftp.math.ucla.edu/pub/camreport/cam08-29.pdf
-    .. [3] http://www.ipol.im/pub/art/2012/g-tvd/article_lr.pdf
+    .. [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
 
     """