diff --git a/skimage/filter/__init__.py b/skimage/filter/__init__.py index f1c1fd49..8894ff1a 100644 --- a/skimage/filter/__init__.py +++ b/skimage/filter/__init__.py @@ -3,7 +3,7 @@ from .ctmf import median_filter from ._canny import canny from .edges import (sobel, hsobel, vsobel, scharr, hscharr, vscharr, prewitt, hprewitt, vprewitt) -from .denoise import tv_denoise, denoise_tv -from ._denoise import denoise_bilateral +from .denoise import tv_denoise +from ._denoise import denoise_bilateral, denoise_tv from ._rank_order import rank_order from .thresholding import threshold_otsu, threshold_adaptive diff --git a/skimage/filter/_denoise.pyx b/skimage/filter/_denoise.pyx index b60a5a85..ca75ab4f 100644 --- a/skimage/filter/_denoise.pyx +++ b/skimage/filter/_denoise.pyx @@ -7,6 +7,7 @@ cimport numpy as cnp import numpy as np from libc.math cimport exp, fabs, sqrt from libc.stdlib cimport malloc, free +from libc.float cimport DBL_MAX from skimage._shared.interpolation cimport get_pixel3d from skimage.util import img_as_float @@ -174,4 +175,170 @@ def denoise_bilateral(image, int win_size=5, sigma_range=None, free(centres) free(total_values) - return out + return np.squeeze(out) + + +cdef inline double _get_elem(double* image, Py_ssize_t rows, Py_ssize_t cols, + Py_ssize_t dims, Py_ssize_t r, Py_ssize_t c, + Py_ssize_t k): + return image[r * cols * dims + c * dims + k] + + +cdef inline void _set_elem(double* image, Py_ssize_t rows, Py_ssize_t cols, + Py_ssize_t dims, Py_ssize_t r, Py_ssize_t c, + Py_ssize_t k, double value): + image[r * cols * dims + c * dims + k] = value + + +cdef inline void _incr_elem(double* image, Py_ssize_t rows, Py_ssize_t cols, + Py_ssize_t dims, Py_ssize_t r, Py_ssize_t c, + Py_ssize_t k, double value): + image[r * cols * dims + c * dims + k] += value + + +def denoise_tv(image, double weight, int max_iter=100, double eps=1e-3): + """Perform total-variation denoising using split-Bregman optimization. + + Total-variation denoising (also know as total-variation regularization) + tries to find an image with less total total-variation under the constraint + of being similar to the input image, which is controlled by the + regularization parameter. + + Parameters + ---------- + image : ndarray + Input data to be denoised (converted using img_as_float`). + weight : float, optional + Denoising weight. The smaller the `weight`, the more denoising (at + the expense of less similarity to the `input`). The regularization + parameter `lambda` is chosen as `2 * weight`. + eps : float, optional + Relative difference of the value of the cost function that determines + the stop criterion. The algorithm stops when:: + + SUM((u(n) - u(n-1))**2) < eps + + max_iter: int, optional + Maximal number of iterations used for the optimization. + + Returns + ------- + u : ndarray + Denoised image. + + References + ---------- + .. [1] http://en.wikipedia.org/wiki/Total_variation_denoising + .. [2] ftp://ftp.math.ucla.edu/pub/camreport/cam08-29.pdf + .. [3] http://www.ipol.im/pub/art/2012/g-tvd/article_lr.pdf + + """ + + image = np.atleast_3d(img_as_float(image)) + + cdef: + Py_ssize_t rows = image.shape[0] + Py_ssize_t cols = image.shape[1] + Py_ssize_t dims = image.shape[2] + Py_ssize_t rows2 = rows + 2 + Py_ssize_t cols2 = cols + 2 + Py_ssize_t r, c, k + + Py_ssize_t total = rows * cols * dims + + shape_ext = (rows2, cols2, dims) + + cnp.ndarray[dtype=cnp.double_t, ndim=3, mode='c'] cimage = \ + np.ascontiguousarray(image) + cnp.ndarray[dtype=cnp.double_t, ndim=3, mode='c'] u = \ + np.zeros(shape_ext, dtype=np.double) + + cnp.ndarray[dtype=cnp.double_t, ndim=3, mode='c'] dx = \ + np.zeros(shape_ext, dtype=np.double) + cnp.ndarray[dtype=cnp.double_t, ndim=3, mode='c'] dy = \ + np.zeros(shape_ext, dtype=np.double) + cnp.ndarray[dtype=cnp.double_t, ndim=3, mode='c'] bx = \ + np.zeros(shape_ext, dtype=np.double) + cnp.ndarray[dtype=cnp.double_t, ndim=3, mode='c'] by = \ + np.zeros(shape_ext, dtype=np.double) + + double* image_data = cimage.data + double* u_data = u.data + + double* dx_data = dx.data + double* dy_data = dy.data + double* bx_data = bx.data + double* by_data = by.data + + double ux, uy, uprev, unew, bxx, byy, dxx, dyy, s + int i = 0 + double lam = 2 * weight + double rmse = DBL_MAX + double norm = (weight + 4 * lam) + + u[1:-1, 1:-1] = image + + # reflect image + u[0, 1:-1] = image[1, :] + u[1:-1, 0] = image[:, 1] + u[-1, 1:-1] = image[-2, :] + u[1:-1, -1] = image[:, -2] + + while i < max_iter and rmse > eps: + + rmse = 0 + + for k in range(dims): + for r in range(1, rows + 1): + for c in range(1, cols + 1): + + uprev = _get_elem(u_data, rows2, cols2, dims, r, c, k) + + # forward derivatives + ux = _get_elem(u_data, rows2, cols2, dims, + r, c+1, k) - uprev + uy = _get_elem(u_data, rows2, cols2, dims, + r+1, c, k) - uprev + + # Gauss-Seidel method + unew = ( + lam * ( + + _get_elem(u_data, rows2, cols2, dims, r+1, c, k) + + _get_elem(u_data, rows2, cols2, dims, r-1, c, k) + + _get_elem(u_data, rows2, cols2, dims, r, c+1, k) + + _get_elem(u_data, rows2, cols2, dims, r, c-1, k) + + + _get_elem(dx_data, rows2, cols2, dims, r, c-1, k) + - _get_elem(dx_data, rows2, cols2, dims, r, c, k) + + _get_elem(dy_data, rows2, cols2, dims, r-1, c, k) + - _get_elem(dy_data, rows2, cols2, dims, r, c, k) + + - _get_elem(bx_data, rows2, cols2, dims, r, c-1, k) + + _get_elem(bx_data, rows2, cols2, dims, r, c, k) + - _get_elem(by_data, rows2, cols2, dims, r-1, c, k) + + _get_elem(by_data, rows2, cols2, dims, r, c, k) + ) + weight * _get_elem(image_data, rows, cols, dims, + r-1, c-1, k) + ) / norm + _set_elem(u_data, rows2, cols2, dims, r, c, k, unew) + + # update root mean square error + rmse += (unew - uprev)**2 + + bxx = _get_elem(bx_data, rows2, cols2, dims, r, c, k) + byy = _get_elem(by_data, rows2, cols2, dims, r, c, k) + + s = sqrt((ux + bxx)**2 + (uy + byy)**2) + dxx = s * lam * (ux + bxx) / (s * lam + 1) + dyy = s * lam * (uy + byy) / (s * lam + 1) + + _set_elem(dx_data, rows2, cols2, dims, r, c, k, dxx) + _set_elem(dy_data, rows2, cols2, dims, r, c, k, dyy) + + _incr_elem(bx_data, rows2, cols2, dims, r, c, k, ux - dxx) + _incr_elem(by_data, rows2, cols2, dims, r, c, k, uy - dyy) + + rmse = sqrt(rmse / total) + i += 1 + + return np.squeeze(u[1:-1, 1:-1])