mirror of
https://github.com/wassname/scikit-image.git
synced 2026-07-10 07:00:31 +08:00
Improve total-varation denoising algorithm
Implementation of fast split-Bregman optimization algorithm in Cython. This implementation also fixes the previously broken 3D version, which darkened the images.
This commit is contained in:
Johannes Schönberger
@@ -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
|
||||
|
||||
+168
-1
@@ -7,6 +7,7 @@ cimport numpy as cnp
|
||||
import numpy as np
|
||||
from libc.math cimport exp, fabs, sqrt
|
||||
from libc.stdlib cimport malloc, free
|
||||
from libc.float cimport DBL_MAX
|
||||
from skimage._shared.interpolation cimport get_pixel3d
|
||||
from skimage.util import img_as_float
|
||||
|
||||
@@ -174,4 +175,170 @@ def denoise_bilateral(image, int win_size=5, sigma_range=None,
|
||||
free(centres)
|
||||
free(total_values)
|
||||
|
||||
return out
|
||||
return np.squeeze(out)
|
||||
|
||||
|
||||
cdef inline double _get_elem(double* image, Py_ssize_t rows, Py_ssize_t cols,
|
||||
Py_ssize_t dims, Py_ssize_t r, Py_ssize_t c,
|
||||
Py_ssize_t k):
|
||||
return image[r * cols * dims + c * dims + k]
|
||||
|
||||
|
||||
cdef inline void _set_elem(double* image, Py_ssize_t rows, Py_ssize_t cols,
|
||||
Py_ssize_t dims, Py_ssize_t r, Py_ssize_t c,
|
||||
Py_ssize_t k, double value):
|
||||
image[r * cols * dims + c * dims + k] = value
|
||||
|
||||
|
||||
cdef inline void _incr_elem(double* image, Py_ssize_t rows, Py_ssize_t cols,
|
||||
Py_ssize_t dims, Py_ssize_t r, Py_ssize_t c,
|
||||
Py_ssize_t k, double value):
|
||||
image[r * cols * dims + c * dims + k] += value
|
||||
|
||||
|
||||
def denoise_tv(image, double weight, int max_iter=100, double eps=1e-3):
|
||||
"""Perform total-variation denoising using split-Bregman optimization.
|
||||
|
||||
Total-variation denoising (also know as total-variation regularization)
|
||||
tries to find an image with less total total-variation under the constraint
|
||||
of being similar to the input image, which is controlled by the
|
||||
regularization parameter.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
image : ndarray
|
||||
Input data to be denoised (converted using img_as_float`).
|
||||
weight : float, optional
|
||||
Denoising weight. The smaller the `weight`, the more denoising (at
|
||||
the expense of less similarity to the `input`). The regularization
|
||||
parameter `lambda` is chosen as `2 * weight`.
|
||||
eps : float, optional
|
||||
Relative difference of the value of the cost function that determines
|
||||
the stop criterion. The algorithm stops when::
|
||||
|
||||
SUM((u(n) - u(n-1))**2) < eps
|
||||
|
||||
max_iter: int, optional
|
||||
Maximal number of iterations used for the optimization.
|
||||
|
||||
Returns
|
||||
-------
|
||||
u : ndarray
|
||||
Denoised image.
|
||||
|
||||
References
|
||||
----------
|
||||
.. [1] http://en.wikipedia.org/wiki/Total_variation_denoising
|
||||
.. [2] ftp://ftp.math.ucla.edu/pub/camreport/cam08-29.pdf
|
||||
.. [3] http://www.ipol.im/pub/art/2012/g-tvd/article_lr.pdf
|
||||
|
||||
"""
|
||||
|
||||
image = np.atleast_3d(img_as_float(image))
|
||||
|
||||
cdef:
|
||||
Py_ssize_t rows = image.shape[0]
|
||||
Py_ssize_t cols = image.shape[1]
|
||||
Py_ssize_t dims = image.shape[2]
|
||||
Py_ssize_t rows2 = rows + 2
|
||||
Py_ssize_t cols2 = cols + 2
|
||||
Py_ssize_t r, c, k
|
||||
|
||||
Py_ssize_t total = rows * cols * dims
|
||||
|
||||
shape_ext = (rows2, cols2, dims)
|
||||
|
||||
cnp.ndarray[dtype=cnp.double_t, ndim=3, mode='c'] cimage = \
|
||||
np.ascontiguousarray(image)
|
||||
cnp.ndarray[dtype=cnp.double_t, ndim=3, mode='c'] u = \
|
||||
np.zeros(shape_ext, dtype=np.double)
|
||||
|
||||
cnp.ndarray[dtype=cnp.double_t, ndim=3, mode='c'] dx = \
|
||||
np.zeros(shape_ext, dtype=np.double)
|
||||
cnp.ndarray[dtype=cnp.double_t, ndim=3, mode='c'] dy = \
|
||||
np.zeros(shape_ext, dtype=np.double)
|
||||
cnp.ndarray[dtype=cnp.double_t, ndim=3, mode='c'] bx = \
|
||||
np.zeros(shape_ext, dtype=np.double)
|
||||
cnp.ndarray[dtype=cnp.double_t, ndim=3, mode='c'] by = \
|
||||
np.zeros(shape_ext, dtype=np.double)
|
||||
|
||||
double* image_data = <double*>cimage.data
|
||||
double* u_data = <double*>u.data
|
||||
|
||||
double* dx_data = <double*>dx.data
|
||||
double* dy_data = <double*>dy.data
|
||||
double* bx_data = <double*>bx.data
|
||||
double* by_data = <double*>by.data
|
||||
|
||||
double ux, uy, uprev, unew, bxx, byy, dxx, dyy, s
|
||||
int i = 0
|
||||
double lam = 2 * weight
|
||||
double rmse = DBL_MAX
|
||||
double norm = (weight + 4 * lam)
|
||||
|
||||
u[1:-1, 1:-1] = image
|
||||
|
||||
# reflect image
|
||||
u[0, 1:-1] = image[1, :]
|
||||
u[1:-1, 0] = image[:, 1]
|
||||
u[-1, 1:-1] = image[-2, :]
|
||||
u[1:-1, -1] = image[:, -2]
|
||||
|
||||
while i < max_iter and rmse > eps:
|
||||
|
||||
rmse = 0
|
||||
|
||||
for k in range(dims):
|
||||
for r in range(1, rows + 1):
|
||||
for c in range(1, cols + 1):
|
||||
|
||||
uprev = _get_elem(u_data, rows2, cols2, dims, r, c, k)
|
||||
|
||||
# forward derivatives
|
||||
ux = _get_elem(u_data, rows2, cols2, dims,
|
||||
r, c+1, k) - uprev
|
||||
uy = _get_elem(u_data, rows2, cols2, dims,
|
||||
r+1, c, k) - uprev
|
||||
|
||||
# Gauss-Seidel method
|
||||
unew = (
|
||||
lam * (
|
||||
+ _get_elem(u_data, rows2, cols2, dims, r+1, c, k)
|
||||
+ _get_elem(u_data, rows2, cols2, dims, r-1, c, k)
|
||||
+ _get_elem(u_data, rows2, cols2, dims, r, c+1, k)
|
||||
+ _get_elem(u_data, rows2, cols2, dims, r, c-1, k)
|
||||
|
||||
+ _get_elem(dx_data, rows2, cols2, dims, r, c-1, k)
|
||||
- _get_elem(dx_data, rows2, cols2, dims, r, c, k)
|
||||
+ _get_elem(dy_data, rows2, cols2, dims, r-1, c, k)
|
||||
- _get_elem(dy_data, rows2, cols2, dims, r, c, k)
|
||||
|
||||
- _get_elem(bx_data, rows2, cols2, dims, r, c-1, k)
|
||||
+ _get_elem(bx_data, rows2, cols2, dims, r, c, k)
|
||||
- _get_elem(by_data, rows2, cols2, dims, r-1, c, k)
|
||||
+ _get_elem(by_data, rows2, cols2, dims, r, c, k)
|
||||
) + weight * _get_elem(image_data, rows, cols, dims,
|
||||
r-1, c-1, k)
|
||||
) / norm
|
||||
_set_elem(u_data, rows2, cols2, dims, r, c, k, unew)
|
||||
|
||||
# update root mean square error
|
||||
rmse += (unew - uprev)**2
|
||||
|
||||
bxx = _get_elem(bx_data, rows2, cols2, dims, r, c, k)
|
||||
byy = _get_elem(by_data, rows2, cols2, dims, r, c, k)
|
||||
|
||||
s = sqrt((ux + bxx)**2 + (uy + byy)**2)
|
||||
dxx = s * lam * (ux + bxx) / (s * lam + 1)
|
||||
dyy = s * lam * (uy + byy) / (s * lam + 1)
|
||||
|
||||
_set_elem(dx_data, rows2, cols2, dims, r, c, k, dxx)
|
||||
_set_elem(dy_data, rows2, cols2, dims, r, c, k, dyy)
|
||||
|
||||
_incr_elem(bx_data, rows2, cols2, dims, r, c, k, ux - dxx)
|
||||
_incr_elem(by_data, rows2, cols2, dims, r, c, k, uy - dyy)
|
||||
|
||||
rmse = sqrt(rmse / total)
|
||||
i += 1
|
||||
|
||||
return np.squeeze(u[1:-1, 1:-1])
|
||||
|
||||
Reference in New Issue
Block a user