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

Implemented fast algorithm also for 3-D and 2D-RGB images. Changed API so

Browse Source 
that there is only one function for fast and classic algorithms.
This commit is contained in:
emmanuelle
2015-01-18 22:18:29 +01:00
parent 8b9a777c79
commit dd9030d44c
5 changed files with 289 additions and 107 deletions
Download Patch File Download Diff File Expand all files Collapse all files
doc/examples/plot_nonlocal_means.py
+5 -5
View File
@@ -3,7 +3,7 @@
Non-local means denoising for preserving textures
=================================================
In this example, we denoise a detail of the Lena image using the non-local
In this example, we denoise a detail of the astronaut image using the non-local
means filter. The non-local means algorithm replaces the value of a pixel by an
average of a selection of other pixels values: small patches centered on the
other pixels are compared to the patch centered on the pixel of interest, and
@@ -18,13 +18,13 @@ from skimage import data, img_as_float
from skimage.restoration import nl_means_denoising
lena = img_as_float(data.lena())
lena = lena[200:300, 100:200]
astro = img_as_float(data.astronaut())
astro = astro[30:180, 150:300]
noisy = lena + 0.6 * lena.std() * np.random.random(lena.shape)
noisy = astro + 0.3 * np.random.random(astro.shape)
noisy = np.clip(noisy, 0, 1)
denoise = nl_means_denoising(noisy, 7, 9, 0.06)
denoise = nl_means_denoising(noisy, 7, 9, 0.08)
fig, ax = plt.subplots(ncols=2, figsize=(8, 4))
skimage/restoration/__init__.py
+2 -3
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 .non_local_means import nl_means_denoising, fast_nl_means_denoising
from .non_local_means import nl_means_denoising
__all__ = ['wiener',
'unsupervised_wiener',
@@ -31,5 +31,4 @@ __all__ = ['wiener',
'denoise_tv_bregman',
'denoise_tv_chambolle',
'denoise_bilateral',
'nl_means_denoising',
'fast_nl_means_denoising']
'nl_means_denoising']
skimage/restoration/_nl_means_denoising.pyx
+215 -5
View File
@@ -23,7 +23,7 @@ cdef inline float patch_distance_2d(DTYPE_t [:, :] p1,
cdef float distance = 0
for i in range(s):
# exp of large negative numbers will be 0, so we'd better stop
if distance > 4:
if distance > 5:
return eps
for j in range(s):
tmp_diff = p1[i, j] - p2[i, j]
@@ -43,7 +43,7 @@ cdef inline float patch_distance_2drgb(DTYPE_t [:, :, :] p1,
cdef float distance = 0
for i in range(s):
# exp of large negative numbers will be 0, so we'd better stop
if distance > 4:
if distance > 5:
return eps
for j in range(s):
for color in range(3):
@@ -62,7 +62,7 @@ cdef inline float patch_distance_3d(DTYPE_t [:, :, :] p1,
cdef float tmp_diff
for i in range(s):
# exp of large negative numbers will be 0, so we'd better stop
if distance > 4:
if distance > 5:
return eps
for j in range(s):
for k in range(s):
@@ -186,7 +186,7 @@ def _nl_means_denoising_2drgb(image, int s=7, int d=13, float h=0.1):
- (xg ** 2 + yg ** 2) / (2 * A ** 2)).
astype(np.float32))
cdef float distance
w = 1. / (np.sum(w) * h ** 2) * w
w = 1. / (3 * np.sum(w) * h ** 2) * w
# Coordinates of central pixel and patch bounds
for x in range(offset, n_x + offset):
x_start = x - offset
@@ -371,7 +371,8 @@ def _fast_nl_means_denoising_2d(image, int s=7, int d=13, float h=0.1):
integral[x - offset, y + offset] - \
integral[x + offset, y - offset]
distance /= (s2 * h2)
if distance > 4:
# exp of large negative numbers is close to zero
if distance > 5:
continue
weight = alpha * exp(- distance)
weights[x, y] += weight
@@ -384,3 +385,212 @@ def _fast_nl_means_denoising_2d(image, int s=7, int d=13, float h=0.1):
# except in padded zone
result[x, y] /= weights[x, y]
return result[pad_size: - pad_size, pad_size: - pad_size]
@cython.cdivision(True)
@cython.boundscheck(False)
def _fast_nl_means_denoising_2drgb(image, int s=7, int d=13, float h=0.1):
"""
Perform fast non-local means denoising on 2-D RGB array, with the outer
loop on patch shifts in order to reduce the number of operations.
Parameters
----------
image: ndarray
2-D RGB 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
# Image padding: we need to account for patch size, possible shift,
# + 1 for the boundary effects in finite differences
cdef int pad_size = offset + d + 1
cdef DTYPE_t [:, :, ::1] padded = np.ascontiguousarray(util.pad(image,
((pad_size, pad_size), (pad_size, pad_size), (0, 0)),
mode='reflect').astype(np.float32))
cdef DTYPE_t [:, :, ::1] result = np.zeros_like(padded)
cdef DTYPE_t [:, ::1] weights = np.zeros_like(padded[..., 0], order='C')
cdef DTYPE_t [:, ::1] integral = np.zeros_like(padded[..., 0], order='C')
cdef int n_x, n_y, t1, t2, x, y
cdef float weight, distance
cdef float alpha
cdef float h2 = h ** 2.
cdef float s2 = s ** 2.
cdef float h2s2 = 3 * h2 * s2
n_x, n_y, _ = image.shape
n_x += 2 * pad_size
n_y += 2 * pad_size
# Outer loops on patch shifts
# With t2 >= 0, reference patch is always on the left of test patch
for t1 in range(-d, d + 1):
for t2 in range(0, d + 1):
# alpha is to account for patches on the same column
# distance is computed twice in this case
if t2 == 0 and t1 is not 0:
alpha = 0.5
else:
alpha = 1.
integral = np.zeros_like(padded[..., 0], order='C')
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)):
distance = ((padded[x, y, 0] -
padded[x + t1, y + t2, 0])**2
+(padded[x, y, 1] -
padded[x + t1, y + t2, 1])**2
+(padded[x, y, 2] -
padded[x + t1, y + t2, 2])**2)
integral[x, y] = distance + \
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 /= h2s2
# exp of large negative numbers is close to zero
if distance > 5:
continue
weight = alpha * exp(- distance)
weights[x, y] += weight
weights[x + t1, y + t2] += weight
for ch in range(3):
result[x, y, ch] += weight * padded[x + t1, y + t2, ch]
result[x + t1, y + t2, ch] += weight * padded[x, y, ch]
for x in range(offset, n_x - offset):
for y in range(offset, n_y - offset):
for channel in range(3):
# no risk of division by zero
# except in padded zone
result[x, y, channel] /= weights[x, y]
return result[pad_size: - pad_size, pad_size: - pad_size]
@cython.cdivision(True)
@cython.boundscheck(False)
def _fast_nl_means_denoising_3d(image, int s=5, int d=7, float h=0.1):
"""
Perform fast non-local means denoising on 3-D array, with the outer
loop on patch shifts in order to reduce the number of operations.
Parameters
----------
image: ndarray
3-D 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
# Image padding: we need to account for patch size, possible shift,
# + 1 for the boundary effects in finite differences
cdef int pad_size = offset + d + 1
cdef DTYPE_t [:, :, ::1] padded = np.ascontiguousarray(util.pad(image,
pad_size, mode='reflect').astype(np.float32))
cdef DTYPE_t [:, :, ::1] result = np.zeros_like(padded)
cdef DTYPE_t [:, :, ::1] weights = np.zeros_like(padded)
cdef DTYPE_t [:, :, ::1] integral = np.zeros_like(padded)
cdef int n_x, n_y, n_z, t1, t2, t3, x, y, z
cdef int x_integral_min, x_integral_max, y__integral_min, y_integral_max, \
z_integral_min, z_integral_max
cdef int x_dist_min, x_dist_max, y_dist_min, y_dist_max, \
z_dist_min, z_dist_max
cdef float weight, distance
cdef float alpha
cdef float h_square = h ** 2.
cdef float s_cube = s ** 3.
cdef float s_cube_h_square = h_square * s_cube
n_x, n_y, n_z = image.shape
n_x += 2 * pad_size
n_y += 2 * pad_size
n_z += 2 * pad_size
# Outer loops on patch shifts
# With t2 >= 0, reference patch is always on the left of test patch
for t1 in range(-d, d + 1):
x_integral_min = max(1, - t1)
x_integral_max = min(n_x, n_x - t1)
x_dist_min = max(offset, offset - t1)
x_dist_max = min(n_x - offset, n_x - offset - t1)
for t2 in range(-d, d + 1):
y_integral_min = max(1, - t2)
y_integral_max = min(n_y, n_y - t2)
y_dist_min = max(offset, offset - t2)
y_dist_max = min(n_y - offset, n_y - offset - t2)
for t3 in range(0, d + 1):
z_integral_min = max(1, - t3)
z_integral_max = min(n_z, n_z - t3)
z_dist_min = max(offset, offset - t3)
z_dist_max = min(n_z - offset, n_z - offset - t3)
# alpha is to account for patches on the same column
# distance is computed twice in this case
if t3 == 0 and (t1 is not 0 or t2 is not 0):
alpha = 0.5
else:
alpha = 1.
integral = np.zeros_like(padded)
for x in range(x_integral_min, x_integral_max):
for y in range(y_integral_min, y_integral_max):
for z in range(z_integral_min, z_integral_max):
integral[x, y, z] = ((padded[x, y, z] -
padded[x + t1, y + t2, z + t3])**2 +
integral[x - 1, y, z] +
integral[x, y - 1, z] +
integral[x, y, z - 1] +
integral[x - 1, y - 1, z - 1]
- integral[x - 1, y - 1, z]
- integral[x, y - 1, z - 1]
- integral[x - 1, y, z - 1])
for x in range(x_dist_min, x_dist_max):
for y in range(y_dist_min, y_dist_max):
for z in range(z_dist_min, z_dist_max):
distance = (integral[x + offset,
y + offset,
z + offset]
- integral[x - offset, y - offset, z - offset]
+ integral[x - offset, y - offset, z + offset]
+ integral[x - offset, y + offset, z - offset]
+ integral[x + offset, y - offset, z - offset]
- integral[x - offset, y + offset, z + offset]
- integral[x + offset, y - offset, z + offset]
- integral[x + offset, y + offset, z - offset])
distance /= s_cube_h_square
# exp of large negative numbers is close to zero
if distance > 5.:
continue
weight = alpha * exp(- distance)
weights[x, y, z] += weight
weights[x + t1, y + t2, z + t3] += weight
result[x, y, z] += weight * padded[x + t1, y + t2,
z + t3]
result[x + t1, y + t2, z + t3] += weight * \
padded[x, y, z]
for x in range(offset, n_x - offset):
for y in range(offset, n_y - offset):
for z in range(offset, n_z - offset):
# I think there is no risk of division by zero
# except in padded zone
result[x, y, z] /= weights[x, y, z]
return result[pad_size: - pad_size, pad_size: - pad_size,
pad_size: -pad_size]
skimage/restoration/non_local_means.py
+52 -84
View File
@@ -1,33 +1,42 @@
import numpy as np
from skimage.restoration._nl_means_denoising import _nl_means_denoising_2d, \
_nl_means_denoising_2drgb, _nl_means_denoising_3d, \
_fast_nl_means_denoising_2d
_fast_nl_means_denoising_2d, _fast_nl_means_denoising_3d, \
_fast_nl_means_denoising_2drgb
def nl_means_denoising(image, patch_size=7, patch_distance=11, h=0.1):
def nl_means_denoising(image, patch_size=7, patch_distance=11, h=0.1,
fast_mode=True):
"""
Perform non-local means denoising on 2-D or 3-D grayscale images, and
2-D RGB images.
Parameters
----------
image: ndarray
image : ndarray
input data to be denoised
patch_size: int, optional
patch_size : int, optional
size of patches used for denoising
patch_distance: int, optional
patch_distance : int, optional
maximal distance in pixels where to search patches used for denoising
h: float, optional
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.
at the expense of blurring features. For a Gaussian noise of standard
deviation sigma, a rule of thumb is to choose the value of h to be
sigma of slightly less.
fast_mode : bool, optional
if True (default value), a fast version of the non-local means
algorithm is used. If False, the original version of non-local means is
used. See the Notes section for more details about the algorithms.
Returns
-------
result: ndarray
result : ndarray
denoised image, of same shape as `image`.
See Also
@@ -43,82 +52,17 @@ def nl_means_denoising(image, patch_size=7, patch_distance=11, h=0.1):
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
In the original version of the algorithm [1]_, corresponding to
``fast=False``, the computational complexity is
image.size * patch_size ** image.ndim * patch_distance ** image.ndim
Hence, changing the size of patches or their maximal distance has a
strong effect on computing times, especially for 3-D images.
The image is padded using the `reflect` mode of `skimage.util.pad`
before denoising.
References
----------
.. [1] Buades, A., Coll, B., & Morel, J. M. (2005, June). A non-local
algorithm for image denoising. In CVPR 2005, Vol. 2, pp. 60-65, IEEE.
Examples
--------
>>> a = np.zeros((40, 40))
>>> a[10:-10, 10:-10] = 1.
>>> a += 0.3*np.random.randn(*a.shape)
>>> denoised_a = nl_means_denoising(a, 7, 5, 0.1)
"""
if image.ndim == 2:
return np.array(_nl_means_denoising_2d(image, s=patch_size,
d=patch_distance, h=h))
if image.ndim == 3 and image.shape[-1] > 4: # only grayscale
return np.array(_nl_means_denoising_3d(image, patch_size,
patch_distance, h))
if image.ndim == 3 and image.shape[-1] == 3: # 2-D color (RGB) images
return np.array(_nl_means_denoising_2drgb(image, patch_size,
patch_distance, h))
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
However, the default behavior corresponds to ``fast=True``, for which
another version of non-local means [2]_ is used, corresponding to a
complexity of
image.size * patch_distance ** image.ndim
@@ -132,14 +76,19 @@ def fast_nl_means_denoising(image, patch_size=7, patch_distance=11, h=0.1):
`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.
slightly poorer denoising performance. Moreover, for small images (images
with a linear size that is only a few times the patch size), the classic
algorithm can be faster due to boundary effects.
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
.. [1] Buades, A., Coll, B., & Morel, J. M. (2005, June). A non-local
algorithm for image denoising. In CVPR 2005, Vol. 2, pp. 60-65, IEEE.
.. [2] Jacques Froment. Parameter-Free Fast Pixelwise Non-Local Means
Denoising. Image Processing On Line, 2014, vol. 4, p. 300-326.
Examples
@@ -147,11 +96,30 @@ def fast_nl_means_denoising(image, patch_size=7, patch_distance=11, h=0.1):
>>> 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)
>>> denoised_a = nl_means_denoising(a, 7, 5, 0.1)
"""
if image.ndim == 2:
return np.array(_fast_nl_means_denoising_2d(image, s=patch_size,
if fast_mode:
return np.array(_fast_nl_means_denoising_2d(image, s=patch_size,
d=patch_distance, h=h))
else:
return np.array(_nl_means_denoising_2d(image, s=patch_size,
d=patch_distance, h=h))
if image.ndim == 3 and image.shape[-1] > 4: # only grayscale
if fast_mode:
return np.array(_fast_nl_means_denoising_3d(image, s=patch_size,
d=patch_distance, h=h))
else:
return np.array(_nl_means_denoising_3d(image, patch_size,
patch_distance, h))
if image.ndim == 3 and image.shape[-1] == 3: # 2-D color (RGB) images
if fast_mode:
return np.array(_fast_nl_means_denoising_2drgb(image, patch_size,
patch_distance, h))
else:
return np.array(_nl_means_denoising_2drgb(image, patch_size,
patch_distance, h))
else:
raise ValueError("Fast non local means denoising is only possible for \
2D grayscale images.")
raise ValueError("Non local means denoising is only possible for \
2D grayscale and RGB images or 3-D grayscale images.")
skimage/restoration/tests/test_denoise.py
+15 -10
View File
@@ -148,16 +148,10 @@ def test_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.nl_means_denoising(img, 7, 5, 0.1)
denoised = restoration.nl_means_denoising(img, 7, 5, 0.1, fast_mode=True)
# make sure noise is reduced
assert img.std() > denoised.std()
def test_fast_nl_means_denoising_2d():
img = np.zeros((40, 50))
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)
denoised = restoration.nl_means_denoising(img, 7, 5, 0.1, fast_mode=False)
# make sure noise is reduced
assert img.std() > denoised.std()
@@ -168,7 +162,10 @@ def test_nl_means_denoising_2drgb():
# add some random noise
img += 0.5 * img.std() * np.random.random(img.shape)
img = np.clip(img, 0, 1)
denoised = restoration.nl_means_denoising(img, 7, 9, 0.08)
denoised = restoration.nl_means_denoising(img, 7, 9, 0.08, fast_mode=True)
# make sure noise is reduced
assert img.std() > denoised.std()
denoised = restoration.nl_means_denoising(img, 7, 9, 0.08, fast_mode=False)
# make sure noise is reduced
assert img.std() > denoised.std()
@@ -177,10 +174,18 @@ def test_nl_means_denoising_3d():
img = np.zeros((20, 20, 10))
img[5:-5, 5:-5, 3:-3] = 1.
img += 0.3*np.random.randn(*img.shape)
denoised = restoration.nl_means_denoising(img, 5, 4, 0.1)
denoised = restoration.nl_means_denoising(img, 5, 4, 0.1, fast_mode=True)
# make sure noise is reduced
assert img.std() > denoised.std()
denoised = restoration.nl_means_denoising(img, 5, 4, 0.1, fast_mode=False)
# make sure noise is reduced
assert img.std() > denoised.std()
def test_nl_means_denoising_wrong_dimension():
img = np.zeros((5, 5, 5, 5))
assert_raises(ValueError, restoration.nl_means_denoising, img)
if __name__ == "__main__":
run_module_suite()