From b0807739efaa87e412d2c0d4a6efe3652e7acb13 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= <hannesschoenberger@gmail.com>
Date: Fri, 5 Oct 2012 11:51:08 +0200
Subject: [PATCH 01/20] Use common file for denoising filters

---
 skimage/filter/__init__.py                    | 2 +-
 skimage/filter/{_tv_denoise.py => denoise.py} | 0
 2 files changed, 1 insertion(+), 1 deletion(-)
 rename skimage/filter/{_tv_denoise.py => denoise.py} (100%)

diff --git a/skimage/filter/__init__.py b/skimage/filter/__init__.py
index f76556df..dc5f9896 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 ._tv_denoise import tv_denoise
+from .denoise import tv_denoise
 from ._rank_order import rank_order
 from .thresholding import threshold_otsu, threshold_adaptive
diff --git a/skimage/filter/_tv_denoise.py b/skimage/filter/denoise.py
similarity index 100%
rename from skimage/filter/_tv_denoise.py
rename to skimage/filter/denoise.py

From 25867b9dd5f916ac22f60075c4206a4aa82b75e1 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= <hannesschoenberger@gmail.com>
Date: Fri, 5 Oct 2012 11:55:27 +0200
Subject: [PATCH 02/20] Improve doc string format of TV denoise functions

---
 skimage/filter/denoise.py | 55 ++++++++++++++++-----------------------
 1 file changed, 22 insertions(+), 33 deletions(-)

diff --git a/skimage/filter/denoise.py b/skimage/filter/denoise.py
index 3302319c..ae4d2f82 100644
--- a/skimage/filter/denoise.py
+++ b/skimage/filter/denoise.py
@@ -3,35 +3,32 @@ from skimage import img_as_float
 
 
 def _tv_denoise_3d(im, weight=100, eps=2.e-4, n_iter_max=200):
-    """
-    Perform total-variation denoising on 3-D arrays
+    """Perform total-variation denoising on 3-D arrays.
 
     Parameters
     ----------
     im: ndarray
-        3-D input data to be denoised
-
+        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``)
-
+        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
+        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.
+        Maximal number of iterations used for the optimization.
 
     Returns
     -------
     out: ndarray
-        denoised array of floats
+        Denoised array of floats.
 
     Notes
     -----
-    Rudin, Osher and Fatemi algorithm
+    Rudin, Osher and Fatemi algorithm.
 
     Examples
     ---------
@@ -86,43 +83,39 @@ def _tv_denoise_3d(im, weight=100, eps=2.e-4, n_iter_max=200):
 
 
 def _tv_denoise_2d(im, weight=50, eps=2.e-4, n_iter_max=200):
-    """
-    Perform total-variation denoising
+    """Perform total-variation denoising.
 
     Parameters
     ----------
     im: ndarray
-        input data to be denoised
-
+        Input data to be denoised.
     weight: float, optional
-        denoising weight. The greater ``weight``, the more denoising (at
+        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
+        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.
+        Maximal number of iterations used for the optimization.
 
     Returns
     -------
     out: ndarray
-        denoised array of floats
+        Denoised array of floats.
 
     Notes
     -----
     The principle of total variation denoising is explained in
-    http://en.wikipedia.org/wiki/Total_variation_denoising
+    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.
@@ -173,33 +166,30 @@ def _tv_denoise_2d(im, weight=50, eps=2.e-4, n_iter_max=200):
 
 
 def tv_denoise(im, weight=50, eps=2.e-4, n_iter_max=200):
-    """
-    Perform total-variation denoising on 2-d and 3-d images
+    """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,
+        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``)
-
+        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
+        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.
+        Maximal number of iterations used for the optimization.
 
     Returns
     -------
     out: ndarray
-        denoised array of floats
+        Denoised array of floats.
 
     Notes
     -----
@@ -217,7 +207,6 @@ def tv_denoise(im, weight=50, eps=2.e-4, n_iter_max=200):
 
     References
     ----------
-
     .. [1] A. Chambolle, An algorithm for total variation minimization and
            applications, Journal of Mathematical Imaging and Vision,
            Springer, 2004, 20, 89-97.

From 75d706ca2558e5b19bc7ae08ead10a6bc466ca5a Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= <hannesschoenberger@gmail.com>
Date: Fri, 5 Oct 2012 23:09:34 +0200
Subject: [PATCH 03/20] Add get_pixel function for 3D arrays

---
 skimage/_shared/interpolation.pxd |  7 ++--
 skimage/_shared/interpolation.pyx | 54 +++++++++++++++++++++++++------
 2 files changed, 49 insertions(+), 12 deletions(-)

diff --git a/skimage/_shared/interpolation.pxd b/skimage/_shared/interpolation.pxd
index f43ff25e..4e00dfb0 100644
--- a/skimage/_shared/interpolation.pxd
+++ b/skimage/_shared/interpolation.pxd
@@ -18,7 +18,10 @@ cdef double bicubic_interpolation(double* image, int rows, int cols,
                                   double r, double c, char mode,
                                   double cval)
 
-cdef double get_pixel(double* image, int rows, int cols, int r, int c,
-                      char mode, double cval)
+cdef double get_pixel2d(double* image, int rows, int cols, int r, int c,
+                        char mode, double cval)
+
+cdef double get_pixel3d(double* image, int rows, int cols, int dims, int r,
+                        int c, int d, char mode, double cval)
 
 cdef int coord_map(int dim, int coord, char mode)
diff --git a/skimage/_shared/interpolation.pyx b/skimage/_shared/interpolation.pyx
index e0fd0067..b34ba507 100644
--- a/skimage/_shared/interpolation.pyx
+++ b/skimage/_shared/interpolation.pyx
@@ -35,8 +35,8 @@ cdef inline double nearest_neighbour_interpolation(double* image, int rows,
 
     """
 
-    return get_pixel(image, rows, cols, <int>round(r), <int>round(c),
-                     mode, cval)
+    return get_pixel2d(image, rows, cols, <int>round(r), <int>round(c),
+                       mode, cval)
 
 
 cdef inline double bilinear_interpolation(double* image, int rows, int cols,
@@ -72,10 +72,10 @@ cdef inline double bilinear_interpolation(double* image, int rows, int cols,
     maxc = <int>ceil(c)
     dr = r - minr
     dc = c - minc
-    top = (1 - dc) * get_pixel(image, rows, cols, minr, minc, mode, cval) \
-          + dc * get_pixel(image, rows, cols, minr, maxc, mode, cval)
-    bottom = (1 - dc) * get_pixel(image, rows, cols, maxr, minc, mode, cval) \
-             + dc * get_pixel(image, rows, cols, maxr, maxc, mode, cval)
+    top = (1 - dc) * get_pixel2d(image, rows, cols, minr, minc, mode, cval) \
+          + dc * get_pixel2d(image, rows, cols, minr, maxc, mode, cval)
+    bottom = (1 - dc) * get_pixel2d(image, rows, cols, maxr, minc, mode, cval) \
+             + dc * get_pixel2d(image, rows, cols, maxr, maxc, mode, cval)
     return (1 - dr) * top + dr * bottom
 
 
@@ -144,7 +144,7 @@ cdef inline double biquadratic_interpolation(double* image, int rows, int cols,
     # row-wise cubic interpolation
     for pr in range(r0, r0 + 3):
         for pc in range(c0, c0 + 3):
-            fc[pc - c0] = get_pixel(image, rows, cols, pr, pc, mode, cval)
+            fc[pc - c0] = get_pixel2d(image, rows, cols, pr, pc, mode, cval)
         fr[pr - r0] = quadratic_interpolation(xc, fc)
 
     # cubic interpolation for interpolated values of each row
@@ -216,15 +216,15 @@ cdef inline double bicubic_interpolation(double* image, int rows, int cols,
     # row-wise cubic interpolation
     for pr in range(r0, r0 + 4):
         for pc in range(c0, c0 + 4):
-            fc[pc - c0] = get_pixel(image, rows, cols, pr, pc, mode, cval)
+            fc[pc - c0] = get_pixel2d(image, rows, cols, pr, pc, mode, cval)
         fr[pr - r0] = cubic_interpolation(xc, fc)
 
     # cubic interpolation for interpolated values of each row
     return cubic_interpolation(xr, fr)
 
 
-cdef inline double get_pixel(double* image, int rows, int cols, int r, int c,
-                             char mode, double cval):
+cdef inline double get_pixel2d(double* image, int rows, int cols, int r, int c,
+                               char mode, double cval):
     """Get a pixel from the image, taking wrapping mode into consideration.
 
     Parameters
@@ -255,6 +255,40 @@ cdef inline double get_pixel(double* image, int rows, int cols, int r, int c,
         return image[coord_map(rows, r, mode) * cols + coord_map(cols, c, mode)]
 
 
+cdef inline double get_pixel3d(double* image, int rows, int cols, int dims, int r,
+                               int c, int d, char mode, double cval):
+    """Get a pixel from the image, taking wrapping mode into consideration.
+
+    Parameters
+    ----------
+    image : double array
+        Input image.
+    rows, cols, dims : int
+        Shape of image.
+    r, c, d : int
+        Position at which to get the pixel.
+    mode : {'C', 'W', 'R', 'N'}
+        Wrapping mode. Constant, Wrap, Reflect or Nearest.
+    cval : double
+        Constant value to use for constant mode.
+
+    Returns
+    -------
+    value : double
+        Pixel value at given position.
+
+    """
+    if mode == 'C':
+        if (r < 0) or (r > rows - 1) or (c < 0) or (c > cols - 1):
+            return cval
+        else:
+            return image[r * cols * dims + c * dims + d]
+    else:
+        return image[coord_map(rows, r, mode) * cols * dims
+                     + coord_map(cols, c, mode) * dims
+                     + d]
+
+
 cdef inline int coord_map(int dim, int coord, char mode):
     """
     Wrap a coordinate, according to a given mode.

From f9ecd140d8942f1db6226df4e7a15d56d79d922c Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= <hannesschoenberger@gmail.com>
Date: Sat, 6 Oct 2012 00:34:05 +0200
Subject: [PATCH 04/20] Add bilateral denoising filter

---
 skimage/filter/__init__.py  |   2 +-
 skimage/filter/_denoise.pyx | 160 ++++++++++++++++++++++++++++++++++++
 skimage/filter/denoise.py   |  57 +++++++++++++
 skimage/filter/setup.py     |   3 +
 4 files changed, 221 insertions(+), 1 deletion(-)
 create mode 100644 skimage/filter/_denoise.pyx

diff --git a/skimage/filter/__init__.py b/skimage/filter/__init__.py
index dc5f9896..eec88038 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 tv_denoise
+from .denoise import tv_denoise, denoise_bilateral
 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
new file mode 100644
index 00000000..bf6c372a
--- /dev/null
+++ b/skimage/filter/_denoise.pyx
@@ -0,0 +1,160 @@
+#cython: cdivision=True
+#cython: boundscheck=False
+#cython: nonecheck=False
+#cython: wraparound=False
+
+cimport numpy as cnp
+import numpy as np
+from libc.math cimport exp, fabs, sqrt
+from libc.stdlib cimport malloc, free
+from skimage._shared.interpolation cimport get_pixel2d, get_pixel3d
+
+
+cdef inline double _gaussian_weight(double sigma, double value):
+    return exp(- 0.5 * (value / sigma)**2)
+
+
+cdef double* _compute_color_lut(int bins, double sigma, double max_value):
+    cdef double* color_lut = <double*>malloc(bins * sizeof(double))
+    cdef Py_ssize_t b
+    for b in range(bins):
+        color_lut[b] = _gaussian_weight(sigma, b * max_value / bins)
+    return color_lut
+
+
+cdef double* _compute_range_lut(int win_size, double sigma):
+    cdef double* range_lut = <double*>malloc(win_size**2 * sizeof(double))
+    cdef Py_ssize_t kr, kc
+    cdef Py_ssize_t window_ext = (win_size - 1) / 2
+    cdef double dist
+    for kr in range(win_size):
+        for kc in range(win_size):
+            dist = sqrt((kr - window_ext)**2 + (kc - window_ext)**2)
+            range_lut[kr * win_size + kc] = _gaussian_weight(sigma, dist)
+    return range_lut
+
+
+def _denoise_bilateral2d(cnp.ndarray[dtype=cnp.double_t, ndim=2, mode='c'] image,
+                         int win_size, double sigma_color,
+                         double sigma_range, int bins, char mode,
+                         double cval):
+
+    cdef Py_ssize_t rows = image.shape[0]
+    cdef Py_ssize_t cols = image.shape[1]
+    cdef Py_ssize_t window_ext = (win_size - 1) / 2
+    cdef double max_value = image.max()
+
+    cdef cnp.ndarray[dtype=cnp.double_t, ndim=2, mode='c'] out = \
+         np.zeros((rows, cols), dtype=np.double)
+
+    cdef double* image_data = <double*>image.data
+    cdef double* out_data = <double*>out.data
+
+    cdef double* color_lut = _compute_color_lut(bins, sigma_color, max_value)
+    cdef double* range_lut = _compute_range_lut(win_size, sigma_range)
+
+    cdef Py_ssize_t r, c, wr, wc, kr, kc, rr, cc, pixel_addr
+    cdef double centre, value, weight, total_value, total_weight, \
+                color_weight, range_weight, diff
+    cdef double dist_scale = bins / max_value
+
+    for r in range(rows):
+        for c in range(cols):
+            pixel_addr = r * cols + c
+            total_value = 0
+            total_weight = 0
+            centre = image_data[pixel_addr]
+            for wr in range(- window_ext, window_ext + 1):
+                rr = wr + r
+                kr = wr + window_ext
+                for wc in range(- window_ext, window_ext + 1):
+                    cc = wc + c
+                    kc = wc + window_ext
+
+                    value = get_pixel2d(image_data, rows, cols,
+                                        rr, cc, mode, cval)
+                    diff = fabs(centre - value)
+
+                    range_weight = range_lut[kr * win_size + kc]
+                    color_weight = color_lut[<int>(diff * dist_scale)]
+
+                    weight = range_weight * color_weight
+                    total_value += value * weight
+                    total_weight += weight
+
+            out_data[pixel_addr] = total_value / total_weight
+
+    free(color_lut)
+    free(range_lut)
+
+    return out
+
+
+def _denoise_bilateral3d(cnp.ndarray[dtype=cnp.double_t, ndim=3, mode='c'] image,
+                         int win_size, double sigma_color,
+                         double sigma_range, int bins, char mode,
+                         double cval):
+
+    cdef Py_ssize_t rows = image.shape[0]
+    cdef Py_ssize_t cols = image.shape[1]
+    cdef Py_ssize_t dims = image.shape[2]
+    cdef Py_ssize_t window_ext = (win_size - 1) / 2
+    cdef double max_value = image.max()
+
+    cdef cnp.ndarray[dtype=cnp.double_t, ndim=3, mode='c'] out = \
+         np.zeros((rows, cols, dims), dtype=np.double)
+
+    cdef double* image_data = <double*>image.data
+    cdef double* out_data = <double*>out.data
+
+    cdef double* color_lut = _compute_color_lut(bins, sigma_color, max_value)
+    cdef double* range_lut = _compute_range_lut(win_size, sigma_range)
+
+    cdef Py_ssize_t r, c, d, wr, wc, kr, kc, rr, cc, pixel_addr
+    cdef double value, weight, dist, total_weight, color_weight, range_weight
+    cdef double dist_scale = bins / dims / max_value
+    cdef double* values = <double*>malloc(dims * sizeof(double))
+    cdef double* centres = <double*>malloc(dims * sizeof(double))
+    cdef double* total_values = <double*>malloc(dims * sizeof(double))
+
+    for r in range(rows):
+        for c in range(cols):
+            pixel_addr = r * cols * dims + c * dims
+            total_weight = 0
+            for d in range(dims):
+                total_values[d] = 0
+                centres[d] = image_data[pixel_addr + d]
+            for wr in range(- window_ext, window_ext + 1):
+                rr = wr + r
+                kr = wr + window_ext
+                for wc in range(- window_ext, window_ext + 1):
+                    cc = wc + c
+                    kc = wc + window_ext
+
+                    # save pixel values for all dims and compute euclidian
+                    # distance between centre stack and current position
+                    dist = 0
+                    for d in range(dims):
+                        value = get_pixel3d(image_data, rows, cols, dims,
+                                            rr, cc, d, mode, cval)
+                        values[d] = value
+                        dist += (centres[d] - value)**2
+                    dist = sqrt(dist)
+
+                    range_weight = range_lut[kr * win_size + kc]
+                    color_weight = color_lut[<int>(dist * dist_scale)]
+
+                    weight = range_weight * color_weight
+                    for d in range(dims):
+                        total_values[d] += values[d] * weight
+                    total_weight += weight
+            for d in range(dims):
+                out_data[pixel_addr + d] = total_values[d] / total_weight
+
+    free(color_lut)
+    free(range_lut)
+    free(values)
+    free(centres)
+    free(total_values)
+
+    return out
diff --git a/skimage/filter/denoise.py b/skimage/filter/denoise.py
index ae4d2f82..f6b30f6f 100644
--- a/skimage/filter/denoise.py
+++ b/skimage/filter/denoise.py
@@ -1,5 +1,6 @@
 import numpy as np
 from skimage import img_as_float
+import _denoise
 
 
 def _tv_denoise_3d(im, weight=100, eps=2.e-4, n_iter_max=200):
@@ -239,3 +240,59 @@ def tv_denoise(im, weight=50, eps=2.e-4, n_iter_max=200):
         raise ValueError('only 2-d and 3-d images may be denoised with this '
                          'function')
     return out
+
+
+def denoise_bilateral(image, win_size=5, sigma_color=1, sigma_range=1, bins=1e4,
+                      mode='constant', cval=0):
+    """Denoise image using bilateral filter.
+
+    Parameters
+    ----------
+    image : ndarray
+        Input image.
+    win_size : int
+        Window size for filtering.
+    sigma_color : float
+        Standard deviation for color distance. A larger value results in
+        averaging of pixels with larger color differences.
+    sigma_range : float
+        Standard deviation for range distance. A larger value results in
+        averaging of pixels with larger spatial differences.
+    bins : int
+        Number of discrete values for gaussian weights of color filtering.
+        A larger value results in improved accuracy.
+    mode : string
+        How to handle values outside the image borders. See
+        `scipy.ndimage.map_coordinates` for detail.
+    cval : string
+        Used in conjunction with mode 'constant', the value outside
+        the image boundaries.
+
+    Returns
+    -------
+    denoised : ndarray
+        Denoised image.
+
+    References
+    ----------
+    .. [1] http://users.soe.ucsc.edu/~manduchi/Papers/ICCV98.pdf
+
+    """
+
+    # not using img_as_float to preserve original range of values, which is
+    # necessary so sigma_color is applied as user desires
+    image = np.array(image, dtype=np.double)
+
+    if mode not in ('constant', 'wrap', 'reflect', 'nearest'):
+        raise ValueError("Invalid mode specified.  Please use "
+                         "`constant`, `nearest`, `wrap` or `reflect`.")
+    mode = ord(mode[0].upper())
+
+    if image.ndim == 2 or (image.ndim == 3 and image.shape[2] == 1):
+        if image.ndim == 3 and image.shape[2] == 1:
+            image = np.squeeze(image)
+        func = _denoise._denoise_bilateral2d
+    else:
+        func = _denoise._denoise_bilateral3d
+    image = np.ascontiguousarray(image)
+    return func(image, win_size, sigma_color, sigma_range, bins, mode, cval)
diff --git a/skimage/filter/setup.py b/skimage/filter/setup.py
index 3e573aec..b996055b 100644
--- a/skimage/filter/setup.py
+++ b/skimage/filter/setup.py
@@ -13,9 +13,12 @@ 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)
 
     config.add_extension('_ctmf', sources=['_ctmf.c'],
                          include_dirs=[get_numpy_include_dirs()])
+    config.add_extension('_denoise', sources=['_denoise.c'],
+                         include_dirs=[get_numpy_include_dirs(), '../_shared'])
 
     return config
 

From cec7b6055d21ae0a9ac65b10d9fbc73feac70979 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= <hannesschoenberger@gmail.com>
Date: Sat, 6 Oct 2012 08:32:08 +0200
Subject: [PATCH 05/20] Restructure cdef declarations

---
 skimage/filter/_denoise.pyx | 85 +++++++++++++++++++++----------------
 1 file changed, 48 insertions(+), 37 deletions(-)

diff --git a/skimage/filter/_denoise.pyx b/skimage/filter/_denoise.pyx
index bf6c372a..b0da5ddb 100644
--- a/skimage/filter/_denoise.pyx
+++ b/skimage/filter/_denoise.pyx
@@ -15,22 +15,30 @@ cdef inline double _gaussian_weight(double sigma, double value):
 
 
 cdef double* _compute_color_lut(int bins, double sigma, double max_value):
-    cdef double* color_lut = <double*>malloc(bins * sizeof(double))
-    cdef Py_ssize_t b
+
+    cdef:
+        double* color_lut = <double*>malloc(bins * sizeof(double))
+        Py_ssize_t b
+
     for b in range(bins):
         color_lut[b] = _gaussian_weight(sigma, b * max_value / bins)
+
     return color_lut
 
 
 cdef double* _compute_range_lut(int win_size, double sigma):
-    cdef double* range_lut = <double*>malloc(win_size**2 * sizeof(double))
-    cdef Py_ssize_t kr, kc
-    cdef Py_ssize_t window_ext = (win_size - 1) / 2
-    cdef double dist
+
+    cdef:
+        double* range_lut = <double*>malloc(win_size**2 * sizeof(double))
+        Py_ssize_t kr, kc
+        Py_ssize_t window_ext = (win_size - 1) / 2
+        double dist
+
     for kr in range(win_size):
         for kc in range(win_size):
             dist = sqrt((kr - window_ext)**2 + (kc - window_ext)**2)
             range_lut[kr * win_size + kc] = _gaussian_weight(sigma, dist)
+
     return range_lut
 
 
@@ -39,24 +47,25 @@ def _denoise_bilateral2d(cnp.ndarray[dtype=cnp.double_t, ndim=2, mode='c'] image
                          double sigma_range, int bins, char mode,
                          double cval):
 
-    cdef Py_ssize_t rows = image.shape[0]
-    cdef Py_ssize_t cols = image.shape[1]
-    cdef Py_ssize_t window_ext = (win_size - 1) / 2
-    cdef double max_value = image.max()
+    cdef:
+        Py_ssize_t rows = image.shape[0]
+        Py_ssize_t cols = image.shape[1]
+        Py_ssize_t window_ext = (win_size - 1) / 2
+        double max_value = image.max()
 
-    cdef cnp.ndarray[dtype=cnp.double_t, ndim=2, mode='c'] out = \
-         np.zeros((rows, cols), dtype=np.double)
+        cnp.ndarray[dtype=cnp.double_t, ndim=2, mode='c'] out = \
+        np.zeros((rows, cols), dtype=np.double)
 
-    cdef double* image_data = <double*>image.data
-    cdef double* out_data = <double*>out.data
+        double* image_data = <double*>image.data
+        double* out_data = <double*>out.data
 
-    cdef double* color_lut = _compute_color_lut(bins, sigma_color, max_value)
-    cdef double* range_lut = _compute_range_lut(win_size, sigma_range)
+        double* color_lut = _compute_color_lut(bins, sigma_color, max_value)
+        double* range_lut = _compute_range_lut(win_size, sigma_range)
 
-    cdef Py_ssize_t r, c, wr, wc, kr, kc, rr, cc, pixel_addr
-    cdef double centre, value, weight, total_value, total_weight, \
-                color_weight, range_weight, diff
-    cdef double dist_scale = bins / max_value
+        Py_ssize_t r, c, wr, wc, kr, kc, rr, cc, pixel_addr
+        double centre, value, weight, total_value, total_weight, \
+               color_weight, range_weight, diff
+        double dist_scale = bins / max_value
 
     for r in range(rows):
         for c in range(cols):
@@ -95,27 +104,29 @@ def _denoise_bilateral3d(cnp.ndarray[dtype=cnp.double_t, ndim=3, mode='c'] image
                          double sigma_range, int bins, char mode,
                          double cval):
 
-    cdef Py_ssize_t rows = image.shape[0]
-    cdef Py_ssize_t cols = image.shape[1]
-    cdef Py_ssize_t dims = image.shape[2]
-    cdef Py_ssize_t window_ext = (win_size - 1) / 2
-    cdef double max_value = image.max()
+    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 window_ext = (win_size - 1) / 2
 
-    cdef cnp.ndarray[dtype=cnp.double_t, ndim=3, mode='c'] out = \
-         np.zeros((rows, cols, dims), dtype=np.double)
+        double max_value = image.max()
 
-    cdef double* image_data = <double*>image.data
-    cdef double* out_data = <double*>out.data
+        cnp.ndarray[dtype=cnp.double_t, ndim=3, mode='c'] out = \
+            np.zeros((rows, cols, dims), dtype=np.double)
 
-    cdef double* color_lut = _compute_color_lut(bins, sigma_color, max_value)
-    cdef double* range_lut = _compute_range_lut(win_size, sigma_range)
+        double* image_data = <double*>image.data
+        double* out_data = <double*>out.data
 
-    cdef Py_ssize_t r, c, d, wr, wc, kr, kc, rr, cc, pixel_addr
-    cdef double value, weight, dist, total_weight, color_weight, range_weight
-    cdef double dist_scale = bins / dims / max_value
-    cdef double* values = <double*>malloc(dims * sizeof(double))
-    cdef double* centres = <double*>malloc(dims * sizeof(double))
-    cdef double* total_values = <double*>malloc(dims * sizeof(double))
+        double* color_lut = _compute_color_lut(bins, sigma_color, max_value)
+        double* range_lut = _compute_range_lut(win_size, sigma_range)
+
+        Py_ssize_t r, c, d, wr, wc, kr, kc, rr, cc, pixel_addr
+        double value, weight, dist, total_weight, color_weight, range_weight
+        double dist_scale = bins / dims / max_value
+        double* values = <double*>malloc(dims * sizeof(double))
+        double* centres = <double*>malloc(dims * sizeof(double))
+        double* total_values = <double*>malloc(dims * sizeof(double))
 
     for r in range(rows):
         for c in range(cols):

From 3757e5c5eb0f142d6f23dcab3cfae98811c014f5 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= <hannesschoenberger@gmail.com>
Date: Sat, 6 Oct 2012 22:03:04 +0200
Subject: [PATCH 06/20] Use common file for denoising filters

---
 skimage/filter/tests/{test_tv_denoise.py => test_denoise.py} | 0
 1 file changed, 0 insertions(+), 0 deletions(-)
 rename skimage/filter/tests/{test_tv_denoise.py => test_denoise.py} (100%)

diff --git a/skimage/filter/tests/test_tv_denoise.py b/skimage/filter/tests/test_denoise.py
similarity index 100%
rename from skimage/filter/tests/test_tv_denoise.py
rename to skimage/filter/tests/test_denoise.py

From 7967d5fb49cd1bb0b1ed8d2417c3ace36f47600d Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= <hannesschoenberger@gmail.com>
Date: Sat, 6 Oct 2012 22:18:41 +0200
Subject: [PATCH 07/20] Refactor denoise tests and add tests for bilateral
 filter

---
 skimage/filter/tests/test_denoise.py | 135 ++++++++++++++++-----------
 1 file changed, 78 insertions(+), 57 deletions(-)

diff --git a/skimage/filter/tests/test_denoise.py b/skimage/filter/tests/test_denoise.py
index cc4fae7e..2efbba84 100644
--- a/skimage/filter/tests/test_denoise.py
+++ b/skimage/filter/tests/test_denoise.py
@@ -1,68 +1,89 @@
 import numpy as np
-from numpy.testing import run_module_suite
+from numpy.testing import run_module_suite, assert_raises
 
-from skimage import filter, data, color
+from skimage import filter, data, color, img_as_float
 
 
-class TestTvDenoise():
+lena = img_as_float(data.lena()[:256, :256])
+lena_gray = color.rgb2gray(lena)
 
-    def test_tv_denoise_2d(self):
-        """
-        Apply the TV denoising algorithm on the lena image provided
-        by scipy
-        """
-        # lena image
-        lena = color.rgb2gray(data.lena())[:256, :256]
-        # add noise to lena
-        lena += 0.5 * lena.std() * np.random.randn(*lena.shape)
-        # clip noise so that it does not exceed allowed range for float images.
-        lena = np.clip(lena, 0, 1)
-        # denoise
-        denoised_lena = filter.tv_denoise(lena, weight=60.0)
-        # which dtype?
-        assert denoised_lena.dtype in [np.float, np.float32, np.float64]
-        from scipy import ndimage
-        grad = ndimage.morphological_gradient(lena, 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_tv_denoise_float_result_range(self):
-        # lena image
-        lena = color.rgb2gray(data.lena())[:256, :256]
-        int_lena = np.multiply(lena, 255).astype(np.uint8)
-        assert np.max(int_lena) > 1
-        denoised_int_lena = filter.tv_denoise(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_tv_denoise_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.tv_denoise(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_tv_denoise_3d(self):
-        """
-        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.randn(*mask.shape)
-        mask[mask < 0] = 0
-        mask[mask > 255] = 255
-        res = filter.tv_denoise(mask.astype(np.uint8), weight=100)
-        assert res.dtype == np.float
-        assert res.std() * 255 < mask.std()
 
-        # test wrong number of dimensions
-        a = np.random.random((8, 8, 8, 8))
-        try:
-            res = filter.tv_denoise(a)
-        except ValueError:
-            pass
+def test_tv_denoise_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.tv_denoise(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_tv_denoise_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.tv_denoise(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.tv_denoise, np.random.random((8, 8, 8, 8)))
+
+
+def test_denoise_bilateral_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_bilateral(img, sigma_color=0.1, sigma_range=20)
+    out2 = filter.denoise_bilateral(img, sigma_color=0.2, sigma_range=30)
+
+    # make sure noise is reduced
+    assert img.std() > out1.std()
+    assert out1.std() > out2.std()
+
+
+def test_denoise_bilateral_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_bilateral(img, sigma_color=0.1, sigma_range=20)
+    out2 = filter.denoise_bilateral(img, sigma_color=0.2, sigma_range=30)
+
+    # make sure noise is reduced
+    assert img.std() > out1.std()
+    assert out1.std() > out2.std()
 
 
 if __name__ == "__main__":

From e39012aa240736fcc570416f2e0ece8f15c9744a Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= <hannesschoenberger@gmail.com>
Date: Sat, 6 Oct 2012 22:31:34 +0200
Subject: [PATCH 08/20] Add more detailed description of bilateral filter

---
 skimage/filter/denoise.py | 11 +++++++++++
 1 file changed, 11 insertions(+)

diff --git a/skimage/filter/denoise.py b/skimage/filter/denoise.py
index f6b30f6f..50ca0247 100644
--- a/skimage/filter/denoise.py
+++ b/skimage/filter/denoise.py
@@ -246,6 +246,17 @@ def denoise_bilateral(image, win_size=5, sigma_color=1, sigma_range=1, bins=1e4,
                       mode='constant', cval=0):
     """Denoise image using bilateral filter.
 
+    This is an edge-preserving and noise reducing denoising filter. It averages
+    pixel based on their spatial closeness and radiometric similarity.
+
+    Spatial closeness is measured by the gaussian function of the euclidian
+    distance between two pixels and a certain standard deviation
+    (`sigma_range`).
+
+    Radiometric similarity is measured by the gaussian function of the euclidian
+    distance between two color values and a certain standard deviation
+    (`sigma_color`).
+
     Parameters
     ----------
     image : ndarray

From 8de3ec988de42d9ceb92e6d436f7d332eb5d337f Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= <hannesschoenberger@gmail.com>
Date: Sat, 6 Oct 2012 22:51:16 +0200
Subject: [PATCH 09/20] Add example script for bilateral denoising

---
 doc/examples/plot_denoise.py         | 64 ++++++++++++++++++++++++++++
 doc/examples/plot_lena_tv_denoise.py | 51 ----------------------
 2 files changed, 64 insertions(+), 51 deletions(-)
 create mode 100644 doc/examples/plot_denoise.py
 delete mode 100644 doc/examples/plot_lena_tv_denoise.py

diff --git a/doc/examples/plot_denoise.py b/doc/examples/plot_denoise.py
new file mode 100644
index 00000000..30dfc10c
--- /dev/null
+++ b/doc/examples/plot_denoise.py
@@ -0,0 +1,64 @@
+"""
+=============================
+Denoising the picture of Lena
+=============================
+
+In this example, we denoise a noisy version of the picture of Lena using the
+total variation and bilateral denoising filter.
+
+These algorithms typically produce "posterized" images with flat domains
+separated by sharp edges. It is possible to change the degree of posterization
+by controlling the tradeoff between denoising and faithfulness to the original
+image.
+
+Total variation filter
+----------------------
+
+The result of this filter is an image that has a minimal total variation norm,
+while being as close to the initial image as possible. The total variation is
+the L1 norm of the gradient of the image, and minimizing the total variation.
+
+Bilateral filter
+----------------
+A bilateral filter is an edge-preserving and noise reducing denoising filter.
+It averages pixel based on their spatial closeness and radiometric similarity.
+
+"""
+
+import numpy as np
+import matplotlib.pyplot as plt
+
+from skimage import data, color, img_as_float
+from skimage.filter import tv_denoise, 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 = np.clip(noisy, 0, 1)
+
+fig, ax = plt.subplots(nrows=2, ncols=3, figsize=(8, 5))
+
+ax[0, 0].imshow(lena)
+ax[0, 0].axis('off')
+ax[0, 0].set_title('original')
+ax[0, 1].imshow(tv_denoise(noisy, weight=0.02))
+ax[0, 1].axis('off')
+ax[0, 1].set_title('TV')
+ax[0, 2].imshow(tv_denoise(noisy, weight=0.05))
+ax[0, 2].axis('off')
+ax[0, 2].set_title('(more) TV')
+
+ax[1, 0].imshow(noisy)
+ax[1, 0].axis('off')
+ax[1, 0].set_title('original')
+ax[1, 1].imshow(denoise_bilateral(noisy, sigma_color=0.02, sigma_range=15))
+ax[1, 1].axis('off')
+ax[1, 1].set_title('Bilateral')
+ax[1, 2].imshow(denoise_bilateral(noisy, sigma_color=0.05, sigma_range=15))
+ax[1, 2].axis('off')
+ax[1, 2].set_title('(more) Bilateral')
+
+fig.subplots_adjust(wspace=0.02, hspace=0.2, top=0.9, bottom=0.05, left=0, right=1)
+
+plt.show()
diff --git a/doc/examples/plot_lena_tv_denoise.py b/doc/examples/plot_lena_tv_denoise.py
deleted file mode 100644
index bce07845..00000000
--- a/doc/examples/plot_lena_tv_denoise.py
+++ /dev/null
@@ -1,51 +0,0 @@
-"""
-====================================================
-Denoising the picture of Lena using total variation
-====================================================
-
-In this example, we denoise a noisy version of the picture of Lena
-using the total variation denoising filter. The result of this filter
-is an image that has a minimal total variation norm, while being as
-close to the initial image as possible. The total variation is the L1
-norm of the gradient of the image, and minimizing the total variation
-typically produces "posterized" images with flat domains separated by
-sharp edges.
-
-It is possible to change the degree of posterization by controlling
-the tradeoff between denoising and faithfulness to the original image.
-
-"""
-
-import numpy as np
-import matplotlib.pyplot as plt
-
-from skimage import data, color, img_as_ubyte
-from skimage.filter import tv_denoise
-
-l = img_as_ubyte(color.rgb2gray(data.lena()))
-l = l[230:290, 220:320]
-
-noisy = l + 0.4 * l.std() * np.random.random(l.shape)
-
-tv_denoised = tv_denoise(noisy, weight=10)
-
-plt.figure(figsize=(8, 2))
-
-plt.subplot(131)
-plt.imshow(noisy, cmap=plt.cm.gray, vmin=40, vmax=220)
-plt.axis('off')
-plt.title('noisy', fontsize=20)
-plt.subplot(132)
-plt.imshow(tv_denoised, cmap=plt.cm.gray, vmin=40, vmax=220)
-plt.axis('off')
-plt.title('TV denoising', fontsize=20)
-
-tv_denoised = tv_denoise(noisy, weight=50)
-plt.subplot(133)
-plt.imshow(tv_denoised, cmap=plt.cm.gray, vmin=40, vmax=220)
-plt.axis('off')
-plt.title('(more) TV denoising', fontsize=20)
-
-plt.subplots_adjust(wspace=0.02, hspace=0.02, top=0.9, bottom=0, left=0,
-                    right=1)
-plt.show()

From 415a3c69b368e587aac158ea852785dba472f88b Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= <hannesschoenberger@gmail.com>
Date: Sat, 6 Oct 2012 22:52:05 +0200
Subject: [PATCH 10/20] Wrap column

---
 doc/examples/plot_denoise.py | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/doc/examples/plot_denoise.py b/doc/examples/plot_denoise.py
index 30dfc10c..4a05e339 100644
--- a/doc/examples/plot_denoise.py
+++ b/doc/examples/plot_denoise.py
@@ -59,6 +59,7 @@ ax[1, 2].imshow(denoise_bilateral(noisy, sigma_color=0.05, sigma_range=15))
 ax[1, 2].axis('off')
 ax[1, 2].set_title('(more) Bilateral')
 
-fig.subplots_adjust(wspace=0.02, hspace=0.2, top=0.9, bottom=0.05, left=0, right=1)
+fig.subplots_adjust(wspace=0.02, hspace=0.2,
+                    top=0.9, bottom=0.05, left=0, right=1)
 
 plt.show()

From 52962f30f1215b705bee0beb70b40819ddf0164e Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= <hannesschoenberger@gmail.com>
Date: Tue, 16 Oct 2012 16:50:54 +0200
Subject: [PATCH 11/20] Reorder subplots

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

diff --git a/doc/examples/plot_denoise.py b/doc/examples/plot_denoise.py
index 4a05e339..45e5e4b5 100644
--- a/doc/examples/plot_denoise.py
+++ b/doc/examples/plot_denoise.py
@@ -39,25 +39,25 @@ noisy = np.clip(noisy, 0, 1)
 
 fig, ax = plt.subplots(nrows=2, ncols=3, figsize=(8, 5))
 
-ax[0, 0].imshow(lena)
+ax[0, 0].imshow(noisy)
 ax[0, 0].axis('off')
-ax[0, 0].set_title('original')
-ax[0, 1].imshow(tv_denoise(noisy, weight=0.02))
+ax[0, 0].set_title('noisy')
+ax[0, 1].imshow(tv_denoise(noisy, weight=0.1))
 ax[0, 1].axis('off')
 ax[0, 1].set_title('TV')
-ax[0, 2].imshow(tv_denoise(noisy, weight=0.05))
+ax[0, 2].imshow(denoise_bilateral(noisy, sigma_color=0.03, sigma_range=15))
 ax[0, 2].axis('off')
-ax[0, 2].set_title('(more) TV')
+ax[0, 2].set_title('Bilateral')
 
-ax[1, 0].imshow(noisy)
+ax[1, 0].imshow(tv_denoise(noisy, weight=0.2))
 ax[1, 0].axis('off')
-ax[1, 0].set_title('original')
-ax[1, 1].imshow(denoise_bilateral(noisy, sigma_color=0.02, sigma_range=15))
+ax[1, 0].set_title('(more) TV')
+ax[1, 1].imshow(denoise_bilateral(noisy, sigma_color=0.06, sigma_range=15))
 ax[1, 1].axis('off')
-ax[1, 1].set_title('Bilateral')
-ax[1, 2].imshow(denoise_bilateral(noisy, sigma_color=0.05, sigma_range=15))
+ax[1, 1].set_title('(more) Bilateral')
+ax[1, 2].imshow(lena)
 ax[1, 2].axis('off')
-ax[1, 2].set_title('(more) Bilateral')
+ax[1, 2].set_title('original')
 
 fig.subplots_adjust(wspace=0.02, hspace=0.2,
                     top=0.9, bottom=0.05, left=0, right=1)

From 62c1beccb59d714b96e9597fbbadd0bddab55972 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= <hannesschoenberger@gmail.com>
Date: Tue, 16 Oct 2012 17:28:42 +0200
Subject: [PATCH 12/20] Fix description of denoising filters

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

diff --git a/doc/examples/plot_denoise.py b/doc/examples/plot_denoise.py
index 45e5e4b5..12e8b9b6 100644
--- a/doc/examples/plot_denoise.py
+++ b/doc/examples/plot_denoise.py
@@ -16,12 +16,13 @@ Total variation filter
 
 The result of this filter is an image that has a minimal total variation norm,
 while being as close to the initial image as possible. The total variation is
-the L1 norm of the gradient of the image, and minimizing the total variation.
+the L1 norm of the gradient of the image.
 
 Bilateral filter
 ----------------
-A bilateral filter is an edge-preserving and noise reducing denoising filter.
-It averages pixel based on their spatial closeness and radiometric similarity.
+
+A bilateral filter is an edge-preserving and noise reducing filter. It averages
+pixels based on their spatial closeness and radiometric similarity.
 
 """
 

From 1a1560f69679240507668878a1c69e794554f42a Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= <hannesschoenberger@gmail.com>
Date: Tue, 16 Oct 2012 17:29:39 +0200
Subject: [PATCH 13/20] Fix PEP8 issues

---
 skimage/filter/_denoise.pyx | 10 +++++-----
 1 file changed, 5 insertions(+), 5 deletions(-)

diff --git a/skimage/filter/_denoise.pyx b/skimage/filter/_denoise.pyx
index b0da5ddb..ad25d9d0 100644
--- a/skimage/filter/_denoise.pyx
+++ b/skimage/filter/_denoise.pyx
@@ -11,7 +11,7 @@ from skimage._shared.interpolation cimport get_pixel2d, get_pixel3d
 
 
 cdef inline double _gaussian_weight(double sigma, double value):
-    return exp(- 0.5 * (value / sigma)**2)
+    return exp(-0.5 * (value / sigma)**2)
 
 
 cdef double* _compute_color_lut(int bins, double sigma, double max_value):
@@ -73,10 +73,10 @@ def _denoise_bilateral2d(cnp.ndarray[dtype=cnp.double_t, ndim=2, mode='c'] image
             total_value = 0
             total_weight = 0
             centre = image_data[pixel_addr]
-            for wr in range(- window_ext, window_ext + 1):
+            for wr in range(-window_ext, window_ext + 1):
                 rr = wr + r
                 kr = wr + window_ext
-                for wc in range(- window_ext, window_ext + 1):
+                for wc in range(-window_ext, window_ext + 1):
                     cc = wc + c
                     kc = wc + window_ext
 
@@ -135,10 +135,10 @@ def _denoise_bilateral3d(cnp.ndarray[dtype=cnp.double_t, ndim=3, mode='c'] image
             for d in range(dims):
                 total_values[d] = 0
                 centres[d] = image_data[pixel_addr + d]
-            for wr in range(- window_ext, window_ext + 1):
+            for wr in range(-window_ext, window_ext + 1):
                 rr = wr + r
                 kr = wr + window_ext
-                for wc in range(- window_ext, window_ext + 1):
+                for wc in range(-window_ext, window_ext + 1):
                     cc = wc + c
                     kc = wc + window_ext
 

From 108cbf90f8cbc48d5dc98e2c8d2be88f02c89735 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= <hannesschoenberger@gmail.com>
Date: Tue, 16 Oct 2012 17:34:49 +0200
Subject: [PATCH 14/20] Fix typo in doc string of biliteral filter

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

diff --git a/skimage/filter/denoise.py b/skimage/filter/denoise.py
index 50ca0247..f5c9e0ee 100644
--- a/skimage/filter/denoise.py
+++ b/skimage/filter/denoise.py
@@ -247,7 +247,7 @@ def denoise_bilateral(image, win_size=5, sigma_color=1, sigma_range=1, bins=1e4,
     """Denoise image using bilateral filter.
 
     This is an edge-preserving and noise reducing denoising filter. It averages
-    pixel based on their spatial closeness and radiometric similarity.
+    pixels based on their spatial closeness and radiometric similarity.
 
     Spatial closeness is measured by the gaussian function of the euclidian
     distance between two pixels and a certain standard deviation

From 5b8f554a2e86f0450301844e9ee93a26ab00ee3b Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= <hannesschoenberger@gmail.com>
Date: Wed, 17 Oct 2012 12:09:09 +0200
Subject: [PATCH 15/20] Rename tv_denoise to denoise_tv and deprecate

---
 skimage/filter/__init__.py           |  2 +-
 skimage/filter/denoise.py            | 22 +++++++++++++---------
 skimage/filter/tests/test_denoise.py | 14 +++++++-------
 3 files changed, 21 insertions(+), 17 deletions(-)

diff --git a/skimage/filter/__init__.py b/skimage/filter/__init__.py
index eec88038..02fddf23 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 tv_denoise, denoise_bilateral
+from .denoise import tv_denoise, denoise_tv, denoise_bilateral
 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
index f5c9e0ee..d7c96f90 100644
--- a/skimage/filter/denoise.py
+++ b/skimage/filter/denoise.py
@@ -1,9 +1,10 @@
 import numpy as np
 from skimage import img_as_float
+from skimage._shared.utils import deprecated
 import _denoise
 
 
-def _tv_denoise_3d(im, weight=100, eps=2.e-4, n_iter_max=200):
+def _denoise_tv_3d(im, weight=100, eps=2.e-4, n_iter_max=200):
     """Perform total-variation denoising on 3-D arrays.
 
     Parameters
@@ -38,7 +39,7 @@ def _tv_denoise_3d(im, weight=100, eps=2.e-4, n_iter_max=200):
     >>> 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 = tv_denoise_3d(mask, weight=100)
+    >>> res = denoise_tv_3d(mask, weight=100)
     """
     px = np.zeros_like(im)
     py = np.zeros_like(im)
@@ -83,7 +84,7 @@ def _tv_denoise_3d(im, weight=100, eps=2.e-4, n_iter_max=200):
     return out
 
 
-def _tv_denoise_2d(im, weight=50, eps=2.e-4, n_iter_max=200):
+def _denoise_tv_2d(im, weight=50, eps=2.e-4, n_iter_max=200):
     """Perform total-variation denoising.
 
     Parameters
@@ -128,7 +129,7 @@ def _tv_denoise_2d(im, weight=50, eps=2.e-4, n_iter_max=200):
     >>> import scipy
     >>> lena = scipy.lena().astype(np.float)
     >>> lena += 0.5 * lena.std()*np.random.randn(*lena.shape)
-    >>> denoised_lena = tv_denoise(lena, weight=60.0)
+    >>> denoised_lena = denoise_tv(lena, weight=60.0)
     """
     px = np.zeros_like(im)
     py = np.zeros_like(im)
@@ -166,7 +167,7 @@ def _tv_denoise_2d(im, weight=50, eps=2.e-4, n_iter_max=200):
     return out
 
 
-def tv_denoise(im, weight=50, eps=2.e-4, n_iter_max=200):
+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
@@ -220,28 +221,31 @@ def tv_denoise(im, weight=50, eps=2.e-4, n_iter_max=200):
     >>> import scipy
     >>> lena = scipy.lena().astype(np.float)
     >>> lena += 0.5 * lena.std()*np.random.randn(*lena.shape)
-    >>> denoised_lena = tv_denoise(lena, weight=60)
+    >>> 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 = tv_denoise_3d(mask, weight=100)
+    >>> 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 = _tv_denoise_2d(im, weight, eps, n_iter_max)
+        out = _denoise_tv_2d(im, weight, eps, n_iter_max)
     elif im.ndim == 3:
-        out = _tv_denoise_3d(im, weight, eps, n_iter_max)
+        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)
+
+
 def denoise_bilateral(image, win_size=5, sigma_color=1, sigma_range=1, bins=1e4,
                       mode='constant', cval=0):
     """Denoise image using bilateral filter.
diff --git a/skimage/filter/tests/test_denoise.py b/skimage/filter/tests/test_denoise.py
index 2efbba84..b8e6c9ef 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_tv_denoise_2d():
+def test_denoise_tv_2d():
     # lena image
     img = lena_gray
     # add noise to lena
@@ -16,7 +16,7 @@ def test_tv_denoise_2d():
     # clip noise so that it does not exceed allowed range for float images.
     img = np.clip(img, 0, 1)
     # denoise
-    denoised_lena = filter.tv_denoise(img, weight=60.0)
+    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
@@ -29,19 +29,19 @@ def test_tv_denoise_2d():
             < np.sqrt((grad**2).sum()) / 2)
 
 
-def test_tv_denoise_float_result_range():
+def test_denoise_tv_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.tv_denoise(int_lena, weight=60.0)
+    denoised_int_lena = filter.denoise_tv(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_tv_denoise_3d():
+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
@@ -50,12 +50,12 @@ def test_tv_denoise_3d():
     mask += 20 * np.random.random(mask.shape)
     mask[mask < 0] = 0
     mask[mask > 255] = 255
-    res = filter.tv_denoise(mask.astype(np.uint8), weight=100)
+    res = filter.denoise_tv(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.tv_denoise, np.random.random((8, 8, 8, 8)))
+    assert_raises(ValueError, filter.denoise_tv, np.random.random((8, 8, 8, 8)))
 
 
 def test_denoise_bilateral_2d():

From c04f14c12941572c61b3046a8e058dbd9cc16c8f Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= <hannesschoenberger@gmail.com>
Date: Wed, 17 Oct 2012 16:28:44 +0200
Subject: [PATCH 16/20] Combine 2D and 3D bilateral filter in one function and
 rename sigma parameters

---
 doc/examples/plot_denoise.py         |   4 +-
 skimage/filter/__init__.py           |   3 +-
 skimage/filter/_denoise.pyx          | 122 ++++++++++++++-------------
 skimage/filter/denoise.py            |  67 ---------------
 skimage/filter/tests/test_denoise.py |   8 +-
 5 files changed, 71 insertions(+), 133 deletions(-)

diff --git a/doc/examples/plot_denoise.py b/doc/examples/plot_denoise.py
index 12e8b9b6..a3e83a65 100644
--- a/doc/examples/plot_denoise.py
+++ b/doc/examples/plot_denoise.py
@@ -46,14 +46,14 @@ ax[0, 0].set_title('noisy')
 ax[0, 1].imshow(tv_denoise(noisy, weight=0.1))
 ax[0, 1].axis('off')
 ax[0, 1].set_title('TV')
-ax[0, 2].imshow(denoise_bilateral(noisy, sigma_color=0.03, sigma_range=15))
+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(tv_denoise(noisy, weight=0.2))
 ax[1, 0].axis('off')
 ax[1, 0].set_title('(more) TV')
-ax[1, 1].imshow(denoise_bilateral(noisy, sigma_color=0.06, sigma_range=15))
+ax[1, 1].imshow(denoise_bilateral(noisy, sigma_range=0.06, 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 02fddf23..f1c1fd49 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 tv_denoise, denoise_tv, denoise_bilateral
+from .denoise import tv_denoise, denoise_tv
+from ._denoise import denoise_bilateral
 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 ad25d9d0..bda46b0d 100644
--- a/skimage/filter/_denoise.pyx
+++ b/skimage/filter/_denoise.pyx
@@ -7,7 +7,8 @@ cimport numpy as cnp
 import numpy as np
 from libc.math cimport exp, fabs, sqrt
 from libc.stdlib cimport malloc, free
-from skimage._shared.interpolation cimport get_pixel2d, get_pixel3d
+from skimage._shared.interpolation cimport get_pixel3d
+from skimage.util import img_as_float
 
 
 cdef inline double _gaussian_weight(double sigma, double value):
@@ -42,67 +43,57 @@ cdef double* _compute_range_lut(int win_size, double sigma):
     return range_lut
 
 
-def _denoise_bilateral2d(cnp.ndarray[dtype=cnp.double_t, ndim=2, mode='c'] image,
-                         int win_size, double sigma_color,
-                         double sigma_range, int bins, char mode,
-                         double cval):
+def denoise_bilateral(image, int win_size=5, sigma_range=None,
+                      double sigma_spatial=1, int bins=10000, mode='constant',
+                      double cval=0):
+    """Denoise image using bilateral filter.
 
-    cdef:
-        Py_ssize_t rows = image.shape[0]
-        Py_ssize_t cols = image.shape[1]
-        Py_ssize_t window_ext = (win_size - 1) / 2
-        double max_value = image.max()
+    This is an edge-preserving and noise reducing denoising filter. It averages
+    pixels based on their spatial closeness and radiometric similarity.
 
-        cnp.ndarray[dtype=cnp.double_t, ndim=2, mode='c'] out = \
-        np.zeros((rows, cols), dtype=np.double)
+    Spatial closeness is measured by the gaussian function of the euclidian
+    distance between two pixels and a certain standard deviation
+    (`sigma_spatial`).
 
-        double* image_data = <double*>image.data
-        double* out_data = <double*>out.data
+    Radiometric similarity is measured by the gaussian function of the euclidian
+    distance between two color values and a certain standard deviation
+    (`sigma_range`).
 
-        double* color_lut = _compute_color_lut(bins, sigma_color, max_value)
-        double* range_lut = _compute_range_lut(win_size, sigma_range)
+    Parameters
+    ----------
+    image : ndarray
+        Input image.
+    win_size : int
+        Window size for filtering.
+    sigma_range : float
+        Standard deviation for grayvalue/color distance (radiometric
+        similarity). A larger value results in averaging of pixels with larger
+        radiometric differences.
+    sigma_spatial : float
+        Standard deviation for range distance. A larger value results in
+        averaging of pixels with larger spatial differences.
+    bins : int
+        Number of discrete values for gaussian weights of color filtering.
+        A larger value results in improved accuracy.
+    mode : string
+        How to handle values outside the image borders. See
+        `scipy.ndimage.map_coordinates` for detail.
+    cval : string
+        Used in conjunction with mode 'constant', the value outside
+        the image boundaries.
 
-        Py_ssize_t r, c, wr, wc, kr, kc, rr, cc, pixel_addr
-        double centre, value, weight, total_value, total_weight, \
-               color_weight, range_weight, diff
-        double dist_scale = bins / max_value
+    Returns
+    -------
+    denoised : ndarray
+        Denoised image.
 
-    for r in range(rows):
-        for c in range(cols):
-            pixel_addr = r * cols + c
-            total_value = 0
-            total_weight = 0
-            centre = image_data[pixel_addr]
-            for wr in range(-window_ext, window_ext + 1):
-                rr = wr + r
-                kr = wr + window_ext
-                for wc in range(-window_ext, window_ext + 1):
-                    cc = wc + c
-                    kc = wc + window_ext
+    References
+    ----------
+    .. [1] http://users.soe.ucsc.edu/~manduchi/Papers/ICCV98.pdf
 
-                    value = get_pixel2d(image_data, rows, cols,
-                                        rr, cc, mode, cval)
-                    diff = fabs(centre - value)
+    """
 
-                    range_weight = range_lut[kr * win_size + kc]
-                    color_weight = color_lut[<int>(diff * dist_scale)]
-
-                    weight = range_weight * color_weight
-                    total_value += value * weight
-                    total_weight += weight
-
-            out_data[pixel_addr] = total_value / total_weight
-
-    free(color_lut)
-    free(range_lut)
-
-    return out
-
-
-def _denoise_bilateral3d(cnp.ndarray[dtype=cnp.double_t, ndim=3, mode='c'] image,
-                         int win_size, double sigma_color,
-                         double sigma_range, int bins, char mode,
-                         double cval):
+    image = np.atleast_3d(img_as_float(image))
 
     cdef:
         Py_ssize_t rows = image.shape[0]
@@ -112,22 +103,35 @@ def _denoise_bilateral3d(cnp.ndarray[dtype=cnp.double_t, ndim=3, mode='c'] image
 
         double max_value = image.max()
 
+        cnp.ndarray[dtype=cnp.double_t, ndim=3, mode='c'] cimage = \
+            np.ascontiguousarray(image)
         cnp.ndarray[dtype=cnp.double_t, ndim=3, mode='c'] out = \
             np.zeros((rows, cols, dims), dtype=np.double)
 
-        double* image_data = <double*>image.data
+        double* image_data = <double*>cimage.data
         double* out_data = <double*>out.data
 
-        double* color_lut = _compute_color_lut(bins, sigma_color, max_value)
-        double* range_lut = _compute_range_lut(win_size, sigma_range)
+        double* color_lut = _compute_color_lut(bins, sigma_range, max_value)
+        double* range_lut = _compute_range_lut(win_size, sigma_spatial)
 
         Py_ssize_t r, c, d, wr, wc, kr, kc, rr, cc, pixel_addr
-        double value, weight, dist, total_weight, color_weight, range_weight
+        double value, weight, dist, total_weight, csigma_range, color_weight, \
+               range_weight
         double dist_scale = bins / dims / max_value
         double* values = <double*>malloc(dims * sizeof(double))
         double* centres = <double*>malloc(dims * sizeof(double))
         double* total_values = <double*>malloc(dims * sizeof(double))
 
+    if sigma_range is None:
+        csigma_range = image.std()
+    else:
+        csigma_range = sigma_range
+
+    if mode not in ('constant', 'wrap', 'reflect', 'nearest'):
+        raise ValueError("Invalid mode specified.  Please use "
+                         "`constant`, `nearest`, `wrap` or `reflect`.")
+    cdef char cmode = ord(mode[0].upper())
+
     for r in range(rows):
         for c in range(cols):
             pixel_addr = r * cols * dims + c * dims
@@ -147,7 +151,7 @@ def _denoise_bilateral3d(cnp.ndarray[dtype=cnp.double_t, ndim=3, mode='c'] image
                     dist = 0
                     for d in range(dims):
                         value = get_pixel3d(image_data, rows, cols, dims,
-                                            rr, cc, d, mode, cval)
+                                            rr, cc, d, cmode, cval)
                         values[d] = value
                         dist += (centres[d] - value)**2
                     dist = sqrt(dist)
diff --git a/skimage/filter/denoise.py b/skimage/filter/denoise.py
index d7c96f90..d62f8b27 100644
--- a/skimage/filter/denoise.py
+++ b/skimage/filter/denoise.py
@@ -244,70 +244,3 @@ def denoise_tv(im, weight=50, eps=2.e-4, n_iter_max=200):
 
 
 tv_denoise = deprecated('skimage.filter.denoise_tv')(denoise_tv)
-
-
-def denoise_bilateral(image, win_size=5, sigma_color=1, sigma_range=1, bins=1e4,
-                      mode='constant', cval=0):
-    """Denoise image using bilateral filter.
-
-    This is an edge-preserving and noise reducing denoising filter. It averages
-    pixels based on their spatial closeness and radiometric similarity.
-
-    Spatial closeness is measured by the gaussian function of the euclidian
-    distance between two pixels and a certain standard deviation
-    (`sigma_range`).
-
-    Radiometric similarity is measured by the gaussian function of the euclidian
-    distance between two color values and a certain standard deviation
-    (`sigma_color`).
-
-    Parameters
-    ----------
-    image : ndarray
-        Input image.
-    win_size : int
-        Window size for filtering.
-    sigma_color : float
-        Standard deviation for color distance. A larger value results in
-        averaging of pixels with larger color differences.
-    sigma_range : float
-        Standard deviation for range distance. A larger value results in
-        averaging of pixels with larger spatial differences.
-    bins : int
-        Number of discrete values for gaussian weights of color filtering.
-        A larger value results in improved accuracy.
-    mode : string
-        How to handle values outside the image borders. See
-        `scipy.ndimage.map_coordinates` for detail.
-    cval : string
-        Used in conjunction with mode 'constant', the value outside
-        the image boundaries.
-
-    Returns
-    -------
-    denoised : ndarray
-        Denoised image.
-
-    References
-    ----------
-    .. [1] http://users.soe.ucsc.edu/~manduchi/Papers/ICCV98.pdf
-
-    """
-
-    # not using img_as_float to preserve original range of values, which is
-    # necessary so sigma_color is applied as user desires
-    image = np.array(image, dtype=np.double)
-
-    if mode not in ('constant', 'wrap', 'reflect', 'nearest'):
-        raise ValueError("Invalid mode specified.  Please use "
-                         "`constant`, `nearest`, `wrap` or `reflect`.")
-    mode = ord(mode[0].upper())
-
-    if image.ndim == 2 or (image.ndim == 3 and image.shape[2] == 1):
-        if image.ndim == 3 and image.shape[2] == 1:
-            image = np.squeeze(image)
-        func = _denoise._denoise_bilateral2d
-    else:
-        func = _denoise._denoise_bilateral3d
-    image = np.ascontiguousarray(image)
-    return func(image, win_size, sigma_color, sigma_range, bins, mode, cval)
diff --git a/skimage/filter/tests/test_denoise.py b/skimage/filter/tests/test_denoise.py
index b8e6c9ef..63a5a5e0 100644
--- a/skimage/filter/tests/test_denoise.py
+++ b/skimage/filter/tests/test_denoise.py
@@ -64,8 +64,8 @@ def test_denoise_bilateral_2d():
     img += 0.5 * img.std() * np.random.random(img.shape)
     img = np.clip(img, 0, 1)
 
-    out1 = filter.denoise_bilateral(img, sigma_color=0.1, sigma_range=20)
-    out2 = filter.denoise_bilateral(img, sigma_color=0.2, sigma_range=30)
+    out1 = filter.denoise_bilateral(img, sigma_range=0.1, sigma_spatial=20)
+    out2 = filter.denoise_bilateral(img, sigma_range=0.2, sigma_spatial=30)
 
     # make sure noise is reduced
     assert img.std() > out1.std()
@@ -78,8 +78,8 @@ def test_denoise_bilateral_3d():
     img += 0.5 * img.std() * np.random.random(img.shape)
     img = np.clip(img, 0, 1)
 
-    out1 = filter.denoise_bilateral(img, sigma_color=0.1, sigma_range=20)
-    out2 = filter.denoise_bilateral(img, sigma_color=0.2, sigma_range=30)
+    out1 = filter.denoise_bilateral(img, sigma_range=0.1, sigma_spatial=20)
+    out2 = filter.denoise_bilateral(img, sigma_range=0.2, sigma_spatial=30)
 
     # make sure noise is reduced
     assert img.std() > out1.std()

From 2cdb3896800a4466e76a29a31138ca026d7a0db6 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= <hannesschoenberger@gmail.com>
Date: Wed, 17 Oct 2012 16:29:32 +0200
Subject: [PATCH 17/20] Remove unused import

---
 skimage/filter/denoise.py | 1 -
 1 file changed, 1 deletion(-)

diff --git a/skimage/filter/denoise.py b/skimage/filter/denoise.py
index d62f8b27..351b2482 100644
--- a/skimage/filter/denoise.py
+++ b/skimage/filter/denoise.py
@@ -1,7 +1,6 @@
 import numpy as np
 from skimage import img_as_float
 from skimage._shared.utils import deprecated
-import _denoise
 
 
 def _denoise_tv_3d(im, weight=100, eps=2.e-4, n_iter_max=200):

From dd421dcc80e14157a0a4e8f6f557201b425713ec Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= <hannesschoenberger@gmail.com>
Date: Wed, 17 Oct 2012 16:32:31 +0200
Subject: [PATCH 18/20] Add note about range of values of sigma_range parameter

---
 skimage/filter/_denoise.pyx | 4 +++-
 1 file changed, 3 insertions(+), 1 deletion(-)

diff --git a/skimage/filter/_denoise.pyx b/skimage/filter/_denoise.pyx
index bda46b0d..b60a5a85 100644
--- a/skimage/filter/_denoise.pyx
+++ b/skimage/filter/_denoise.pyx
@@ -68,7 +68,9 @@ def denoise_bilateral(image, int win_size=5, sigma_range=None,
     sigma_range : float
         Standard deviation for grayvalue/color distance (radiometric
         similarity). A larger value results in averaging of pixels with larger
-        radiometric differences.
+        radiometric differences. Note, that the image will be converted using
+        the `img_as_float` function and thus the standard deviation is in
+        respect to the range `[0, 1]`.
     sigma_spatial : float
         Standard deviation for range distance. A larger value results in
         averaging of pixels with larger spatial differences.

From a41163beb0bd63194f58a03faf5563ae8fe60365 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= <hannesschoenberger@gmail.com>
Date: Wed, 17 Oct 2012 16:45:25 +0200
Subject: [PATCH 19/20] Add Cython denoise source file to bento file

---
 bento.info | 5 ++++-
 1 file changed, 4 insertions(+), 1 deletion(-)

diff --git a/bento.info b/bento.info
index 50a2743a..55b3690d 100644
--- a/bento.info
+++ b/bento.info
@@ -64,7 +64,10 @@ Library:
     Extension: skimage.filter._ctmf
         Sources:
             skimage/filter/_ctmf.pyx
-    Extension: skimage.morphology.ccomp
+    Extension: skimage.filter._denoise
+        Sources:
+            skimage/filter/_ctmf.pyx
+    Extension: skimage.morphology.denoise
         Sources:
             skimage/morphology/ccomp.pyx
     Extension: skimage.morphology._watershed

From 8347bd449fc137cbbc986461b3f4aa5c68c2137e Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= <hannesschoenberger@gmail.com>
Date: Thu, 18 Oct 2012 17:10:07 +0200
Subject: [PATCH 20/20] Replace deprecated function call

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

diff --git a/doc/examples/plot_denoise.py b/doc/examples/plot_denoise.py
index a3e83a65..8f02f21b 100644
--- a/doc/examples/plot_denoise.py
+++ b/doc/examples/plot_denoise.py
@@ -30,7 +30,7 @@ import numpy as np
 import matplotlib.pyplot as plt
 
 from skimage import data, color, img_as_float
-from skimage.filter import tv_denoise, denoise_bilateral
+from skimage.filter import denoise_tv, denoise_bilateral
 
 lena = img_as_float(data.lena())
 lena = lena[220:300, 220:320]
@@ -43,14 +43,14 @@ fig, ax = plt.subplots(nrows=2, ncols=3, figsize=(8, 5))
 ax[0, 0].imshow(noisy)
 ax[0, 0].axis('off')
 ax[0, 0].set_title('noisy')
-ax[0, 1].imshow(tv_denoise(noisy, weight=0.1))
+ax[0, 1].imshow(denoise_tv(noisy, weight=0.1))
 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(tv_denoise(noisy, weight=0.2))
+ax[1, 0].imshow(denoise_tv(noisy, weight=0.2))
 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))