wassname/scikit-image
1
0
Fork 0
mirror of https://github.com/wassname/scikit-image.git synced 2026-07-03 17:15:46 +08:00
Code Issues Packages Projects Releases Wiki Activity

Faster version of non-local means denoising for 2D greyscale images

Browse Source
This commit is contained in:
emmanuelle
2014-12-07 12:02:05 +01:00
parent e26f9028eb
commit 4004f048ef
4 changed files with 176 additions and 5 deletions
Download Patch File Download Diff File Expand all files Collapse all files
skimage/restoration/__init__.py
+3 -2
View File
@@ -22,7 +22,7 @@ from .deconvolution import wiener, unsupervised_wiener, richardson_lucy
from .unwrap import unwrap_phase
from ._denoise import denoise_tv_chambolle, denoise_tv_bregman, \
denoise_bilateral
from .nl_means_denoising import nl_means_denoising
from .nl_means_denoising import nl_means_denoising, fast_nl_means_denoising
__all__ = ['wiener',
'unsupervised_wiener',
@@ -31,4 +31,5 @@ __all__ = ['wiener',
'denoise_tv_bregman',
'denoise_tv_chambolle',
'denoise_bilateral',
'nl_means_denoising']
'nl_means_denoising',
'fast_nl_means_denoising']
skimage/restoration/_nl_means_denoising.pyx
+80 -2
View File
@@ -251,8 +251,9 @@ def _nl_means_denoising_3d(image, int s=7,
n_x, n_y, n_z = image.shape
cdef int offset = s / 2
# padd the image so that boundaries are denoised as well
cdef DTYPE_t [:, :, ::1] padded = np.ascontiguousarray(util.pad(image,
offset, mode='reflect').astype(np.float32))
cdef DTYPE_t [:, :, ::1] padded = np.ascontiguousarray(util.pad(
image.astype(np.float32),
offset, mode='reflect'))
cdef DTYPE_t [:, :, ::1] result = padded.copy()
h *= (np.max(padded) - np.min(padded))
cdef float A = ((s - 1.) / 4.)
@@ -305,3 +306,80 @@ def _nl_means_denoising_3d(image, int s=7,
new_value += weight * padded[x + i, y + j, z + k]
result[x, y, z] = new_value / weight_sum
return result[offset:-offset, offset:-offset, offset:-offset]
@cython.cdivision(True)
@cython.boundscheck(False)
def _fast_nl_means_denoising_2d(image, int s=7, int d=13, float h=0.1):
"""
Perform fast non-local means denoising on 2-D array, with the outer
loop on patch shifts in order to reduce the number of operations.
Parameters
----------
image: ndarray
input data to be denoised
s: int, optional
size of patches used for denoising
d: int, optional
maximal distance in pixels where to search patches used for denoising
h: float, optional
cut-off distance (in gray levels). The higher h, the more permissive
one is in accepting patches.
"""
if s % 2 == 0:
s += 1 # odd value for symmetric patch
cdef int offset = s / 2
cdef int pad_size = offset + d
cdef DTYPE_t [:, ::1] padded = np.ascontiguousarray(util.pad(image,
pad_size, mode='reflect').astype(np.float32))
cdef DTYPE_t [:, ::1] result = padded.copy()
cdef DTYPE_t [:, ::1] weights = np.zeros_like(padded)
cdef DTYPE_t [:, ::1] integral = np.zeros_like(padded)
h *= (np.max(padded) - np.min(padded))
cdef int n_x, n_y, t1, t2, x, y
cdef float weight, distance
cdef float h2 = h**2
cdef float s2 = s**2.
n_x, n_y = image.shape
n_x += 2 * pad_size
n_y += 2 * pad_size
for t1 in range(-d, d + 1):
for t2 in range(-d, d + 1):
integral = np.zeros_like(padded)
for x in range(max(1, - t1), min(n_x, n_x - t1)):
integral[x, 0] = (padded[x, 0] - padded[x + t1, t2])**2 + \
integral[x - 1, 0]
for y in range(max(1, - t2), min(n_y, n_y - t2)):
integral[0, y] = (padded[0, y] - padded[t1, y + t2])**2 + \
integral[0, y - 1]
for x in range(max(1, - t1), min(n_x, n_x - t1)):
for y in range(max(1, - t2), min(n_y, n_y - t2)):
integral[x, y] = (padded[x, y] -
padded[x + t1, y + t2])**2 + \
integral[x - 1, y] + integral[x, y - 1] \
- integral[x - 1, y - 1]
for x in range(max(offset, offset - t1),
min(n_x - offset, n_x - offset - t1)):
for y in range(max(offset, offset - t2),
min(n_y - offset, n_y - offset - t2)):
distance = integral[x + offset, y + offset] + \
integral[x - offset, y - offset] - \
integral[x - offset, y + offset] - \
integral[x + offset, y - offset]
distance /= s2
weight = exp(- distance / h2)
weights[x, y] += weight
result[x, y] += weight * padded[x + t1, y + t2]
for x in range(offset, n_x - offset):
for y in range(offset, n_x - offset):
# I think there is no risk of division by zero
# except in padded zone
result[x, y] /= weights[x, y]
return result[pad_size: - pad_size, pad_size: - pad_size]
skimage/restoration/nl_means_denoising.py
+84 -1
View File
@@ -1,6 +1,7 @@
import numpy as np
from skimage.restoration._nl_means_denoising import _nl_means_denoising_2d, \
_nl_means_denoising_2drgb, _nl_means_denoising_3d
_nl_means_denoising_2drgb, _nl_means_denoising_3d, \
_fast_nl_means_denoising_2d
def nl_means_denoising(image, patch_size=7, patch_distance=11, h=0.1):
"""
@@ -29,6 +30,10 @@ def nl_means_denoising(image, patch_size=7, patch_distance=11, h=0.1):
result: ndarray
denoised image, of same shape as `image`.
See Also
--------
fast_nl_means_denoising
Notes
-----
@@ -72,3 +77,81 @@ def nl_means_denoising(image, patch_size=7, patch_distance=11, h=0.1):
else:
raise ValueError("Non local means denoising is only possible for \
2D grayscale and RGB images or 3-D grayscale images.")
def fast_nl_means_denoising(image, patch_size=7, patch_distance=11, h=0.1):
"""
Performs fast non-local means denoising on 2-D grayscale images.
Parameters
----------
image: ndarray
input data to be denoised
patch_size: int, optional
size of patches used for denoising
patch_distance: int, optional
maximal distance in pixels where to search patches used for denoising
h: float, optional
cut-off distance (in gray levels). The higher h, the more permissive
one is in accepting patches. A higher h results in a smoother image,
at the expense of blurring features.
Returns
-------
result: ndarray
denoised image, of same shape as `image`.
See Also
--------
nl_means_denoising
Notes
-----
The non-local means algorithm is well suited for denoising images with
specific textures. The principle of the algorithm is to average the value
of a given pixel with values of other pixels in a limited neighbourhood,
provided that the *patches* centered on the other pixels are similar enough
to the patch centered on the pixel of interest.
The complexity of the algorithm is
image.size * patch_distance ** image.ndim
The computing time depends only weakly on the patch size, thanks to the
computation of the integral of patches distances for a given shift, that
reduces the number of operations [1]_. Therefore, this algorithm executes
faster than `nl_means_denoising`, at the expense of using twice as much
memory.
Compared to the classic non-local means algorithm implemented in
`nl_means_denoising`, all pixels of a patch contribute to the distance to
another patch with the same weight, no matter their distance to the center
of the patch. This coarser computation of the distance can result in a
slightly poorer denoising performance.
The image is padded using the `reflect` mode of `skimage.util.pad`
before denoising.
References
----------
.. [1] Jacques Froment. Parameter-Free Fast Pixelwise Non-Local Means
Denoising. Image Processing On Line, 2014, vol. 4, p. 300-326.
Examples
--------
>>> a = np.zeros((40, 40))
>>> a[10:-10, 10:-10] = 1.
>>> a += 0.3*np.random.randn(*a.shape)
>>> denoised_a = fast_nl_means_denoising(a, 7, 5, 0.1)
"""
if image.ndim == 2:
return np.array(_fast_nl_means_denoising_2d(image, patch_size,
patch_distance, h))
else:
raise ValueError("Fast non local means denoising is only possible for \
2D grayscale images.")
skimage/restoration/tests/test_denoise.py
+9
View File
@@ -153,6 +153,15 @@ def test_nl_means_denoising_2d():
assert img.std() > denoised.std()
def test_fast_nl_means_denoising_2d():
img = np.zeros((40, 40))
img[10:-10, 10:-10] = 1.
img += 0.3*np.random.randn(*img.shape)
denoised = restoration.fast_nl_means_denoising(img, 7, 5, 0.1)
# make sure noise is reduced
assert img.std() > denoised.std()
def test_nl_means_denoising_2drgb():
# reduce image size because nl means is very slow
img = np.copy(astro[:50, :50])