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 diff --git a/doc/examples/plot_denoise.py b/doc/examples/plot_denoise.py new file mode 100644 index 00000000..8f02f21b --- /dev/null +++ b/doc/examples/plot_denoise.py @@ -0,0 +1,66 @@ +""" +============================= +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. + +Bilateral filter +---------------- + +A bilateral filter is an edge-preserving and noise reducing filter. It averages +pixels 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 denoise_tv, 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(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].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].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].axis('off') +ax[1, 1].set_title('(more) Bilateral') +ax[1, 2].imshow(lena) +ax[1, 2].axis('off') +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) + +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() 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, round(r), round(c), - mode, cval) + return get_pixel2d(image, rows, cols, round(r), 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 = 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. diff --git a/skimage/filter/__init__.py b/skimage/filter/__init__.py index f76556df..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 ._tv_denoise import tv_denoise +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 new file mode 100644 index 00000000..b60a5a85 --- /dev/null +++ b/skimage/filter/_denoise.pyx @@ -0,0 +1,177 @@ +#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_pixel3d +from skimage.util import img_as_float + + +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 = 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 = 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 + + +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. + + 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_spatial`). + + Radiometric similarity is measured by the gaussian function of the euclidian + distance between two color values and a certain standard deviation + (`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. 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. + 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 + + """ + + 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 window_ext = (win_size - 1) / 2 + + 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 = cimage.data + double* out_data = out.data + + 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, csigma_range, color_weight, \ + range_weight + double dist_scale = bins / dims / max_value + double* values = malloc(dims * sizeof(double)) + double* centres = malloc(dims * sizeof(double)) + double* total_values = 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 + 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, cmode, cval) + values[d] = value + dist += (centres[d] - value)**2 + dist = sqrt(dist) + + range_weight = range_lut[kr * win_size + kc] + color_weight = color_lut[(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/_tv_denoise.py b/skimage/filter/denoise.py similarity index 75% rename from skimage/filter/_tv_denoise.py rename to skimage/filter/denoise.py index 3302319c..351b2482 100644 --- a/skimage/filter/_tv_denoise.py +++ b/skimage/filter/denoise.py @@ -1,37 +1,35 @@ import numpy as np from skimage import img_as_float +from skimage._shared.utils import deprecated -def _tv_denoise_3d(im, weight=100, eps=2.e-4, n_iter_max=200): - """ - Perform total-variation denoising on 3-D arrays +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 - + 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 --------- @@ -40,7 +38,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) @@ -85,44 +83,40 @@ 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): - """ - Perform total-variation denoising +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 - + 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. @@ -134,7 +128,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) @@ -172,34 +166,31 @@ 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): - """ - Perform total-variation denoising on 2-d and 3-d images +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, + 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 +208,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. @@ -230,23 +220,26 @@ 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) 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 diff --git a/skimage/filter/tests/test_denoise.py b/skimage/filter/tests/test_denoise.py new file mode 100644 index 00000000..63a5a5e0 --- /dev/null +++ b/skimage/filter/tests/test_denoise.py @@ -0,0 +1,90 @@ +import numpy as np +from numpy.testing import run_module_suite, assert_raises + +from skimage import filter, data, color, img_as_float + + +lena = img_as_float(data.lena()[:256, :256]) +lena_gray = color.rgb2gray(lena) + + +def test_denoise_tv_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(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_float_result_range(): + # lena image + img = lena_gray + int_lena = np.multiply(img, 255).astype(np.uint8) + assert np.max(int_lena) > 1 + denoised_int_lena = filter.denoise_tv(int_lena, weight=60.0) + # test if the value range of output float data is within [0.0:1.0] + assert denoised_int_lena.dtype == np.float + assert np.max(denoised_int_lena) <= 1.0 + assert np.min(denoised_int_lena) >= 0.0 + + +def test_denoise_tv_3d(): + """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() + + # test wrong number of dimensions + assert_raises(ValueError, filter.denoise_tv, 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_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() + 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_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() + assert out1.std() > out2.std() + + +if __name__ == "__main__": + run_module_suite() diff --git a/skimage/filter/tests/test_tv_denoise.py b/skimage/filter/tests/test_tv_denoise.py deleted file mode 100644 index cc4fae7e..00000000 --- a/skimage/filter/tests/test_tv_denoise.py +++ /dev/null @@ -1,69 +0,0 @@ -import numpy as np -from numpy.testing import run_module_suite - -from skimage import filter, data, color - - -class TestTvDenoise(): - - 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_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 - - -if __name__ == "__main__": - run_module_suite()