wassname/scikit-image
1
0
Fork 0
mirror of https://github.com/wassname/scikit-image.git synced 2026-08-14 12:50:39 +08:00
Code Issues Packages Projects Releases Wiki Activity

Replace manual array with builting Cython buffer indexing

Browse Source
This commit is contained in:
Johannes Schönberger
2012-12-26 09:28:25 +01:00
parent ccd2f42643
commit fff33702ea
1 changed files with 23 additions and 52 deletions
Download Patch File Download Diff File Expand all files Collapse all files
skimage/filter/_denoise.pyx
+23 -52
View File
@@ -179,24 +179,6 @@ def denoise_bilateral(image, int win_size=5, sigma_range=None,
return np.squeeze(out)
cdef inline double _get_elem(double* image, Py_ssize_t rows, Py_ssize_t cols,
Py_ssize_t dims, Py_ssize_t r, Py_ssize_t c,
Py_ssize_t k):
return image[r * cols * dims + c * dims + k]
cdef inline void _set_elem(double* image, Py_ssize_t rows, Py_ssize_t cols,
Py_ssize_t dims, Py_ssize_t r, Py_ssize_t c,
Py_ssize_t k, double value):
image[r * cols * dims + c * dims + k] = value
cdef inline void _incr_elem(double* image, Py_ssize_t rows, Py_ssize_t cols,
Py_ssize_t dims, Py_ssize_t r, Py_ssize_t c,
Py_ssize_t k, double value):
image[r * cols * dims + c * dims + k] += value
def denoise_tv(image, double weight, int max_iter=100, double eps=1e-3):
"""Perform total-variation denoising using split-Bregman optimization.
@@ -263,14 +245,6 @@ def denoise_tv(image, double weight, int max_iter=100, double eps=1e-3):
cnp.ndarray[dtype=cnp.double_t, ndim=3, mode='c'] by = \
np.zeros(shape_ext, dtype=np.double)
double* image_data = <double*>cimage.data
double* u_data = <double*>u.data
double* dx_data = <double*>dx.data
double* dy_data = <double*>dy.data
double* bx_data = <double*>bx.data
double* by_data = <double*>by.data
double ux, uy, uprev, unew, bxx, byy, dxx, dyy, s
int i = 0
double lam = 2 * weight
@@ -293,51 +267,48 @@ def denoise_tv(image, double weight, int max_iter=100, double eps=1e-3):
for r in range(1, rows + 1):
for c in range(1, cols + 1):
uprev = _get_elem(u_data, rows2, cols2, dims, r, c, k)
uprev = u[r, c, k]
# forward derivatives
ux = _get_elem(u_data, rows2, cols2, dims,
r, c+1, k) - uprev
uy = _get_elem(u_data, rows2, cols2, dims,
r+1, c, k) - uprev
ux = u[r, c + 1, k] - uprev
uy = u[r + 1, c, k] - uprev
# Gauss-Seidel method
unew = (
lam * (
+ _get_elem(u_data, rows2, cols2, dims, r+1, c, k)
+ _get_elem(u_data, rows2, cols2, dims, r-1, c, k)
+ _get_elem(u_data, rows2, cols2, dims, r, c+1, k)
+ _get_elem(u_data, rows2, cols2, dims, r, c-1, k)
+ u[r + 1, c, k]
+ u[r - 1, c, k]
+ u[r, c + 1, k]
+ u[r, c - 1, k]
+ _get_elem(dx_data, rows2, cols2, dims, r, c-1, k)
- _get_elem(dx_data, rows2, cols2, dims, r, c, k)
+ _get_elem(dy_data, rows2, cols2, dims, r-1, c, k)
- _get_elem(dy_data, rows2, cols2, dims, r, c, k)
+ dx[r, c - 1, k]
- dx[r, c, k]
+ dy[r - 1, c, k]
- dy[r, c, k]
- _get_elem(bx_data, rows2, cols2, dims, r, c-1, k)
+ _get_elem(bx_data, rows2, cols2, dims, r, c, k)
- _get_elem(by_data, rows2, cols2, dims, r-1, c, k)
+ _get_elem(by_data, rows2, cols2, dims, r, c, k)
) + weight * _get_elem(image_data, rows, cols, dims,
r-1, c-1, k)
- bx[r, c - 1, k]
+ bx[r, c, k]
- by[r - 1, c, k]
+ by[r, c, k]
) + weight * cimage[r - 1, c - 1, k]
) / norm
_set_elem(u_data, rows2, cols2, dims, r, c, k, unew)
u[r, c, k] = unew
# update root mean square error
rmse += (unew - uprev)**2
bxx = _get_elem(bx_data, rows2, cols2, dims, r, c, k)
byy = _get_elem(by_data, rows2, cols2, dims, r, c, k)
bxx = bx[r, c, k]
byy = by[r, c, k]
s = sqrt((ux + bxx)**2 + (uy + byy)**2)
dxx = s * lam * (ux + bxx) / (s * lam + 1)
dyy = s * lam * (uy + byy) / (s * lam + 1)
_set_elem(dx_data, rows2, cols2, dims, r, c, k, dxx)
_set_elem(dy_data, rows2, cols2, dims, r, c, k, dyy)
dx[r, c, k] = dxx
dy[r, c, k] = dyy
_incr_elem(bx_data, rows2, cols2, dims, r, c, k, ux - dxx)
_incr_elem(by_data, rows2, cols2, dims, r, c, k, uy - dyy)
bx[r, c, k] += ux - dxx
by[r, c, k] += uy - dyy
rmse = sqrt(rmse / total)
i += 1