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

Some improvements of non-local means denoising:

Browse Source 
- denoising RGB is now possible, and "colored patches" are then compared

- the main function is now in a pure Python file so that default values
  of kw arguments are visible in the help

- reduced the number of computations of patches bound (but this doesn't
  change much the total speed).

- added an example for the gallery

I also played with functions that could replace the exponential by a
faster and less precise function, but it turns out that most of the time
is spent in additions and multiplications when computing the distance
between two patches.
This commit is contained in:
Emmanuelle Gouillart
2014-02-23 13:28:27 +01:00
committed by emmanuelle
parent a508ec54ca
commit a5ed4acf86
4 changed files with 451 additions and 92 deletions
Download Patch File Download Diff File Expand all files Collapse all files
doc/examples/plot_nonlocal_means.py
+41
View File
@@ -0,0 +1,41 @@
"""
=================================================
Non-local means denoising for preserving textures
=================================================
In this example, we denoise a detail of the Lena 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
the average is performed only for pixels that have patches close to the current
patch. As a result, this algorithm can restore well textures, that would be
blurred by other denoising algoritm.
"""
import numpy as np
import matplotlib.pyplot as plt
from skimage import data, img_as_float
from skimage.filter import nl_means_denoising
lena = img_as_float(data.lena())
lena = lena[200:300, 100:200]
noisy = lena + 0.6 * lena.std() * np.random.random(lena.shape)
noisy = np.clip(noisy, 0, 1)
denoise = nl_means_denoising(noisy, 7, 9, 0.06)
fig, ax = plt.subplots(ncols=2, figsize=(8, 4))
ax[0].imshow(noisy)
ax[0].axis('off')
ax[0].set_title('noisy')
ax[1].imshow(denoise)
ax[1].axis('off')
ax[1].set_title('non-local means')
fig.subplots_adjust(wspace=0.02, hspace=0.2,
top=0.9, bottom=0.05, left=0, right=1)
plt.show()
skimage/filter/_nl_means_denoising.pyx
+167 -92
View File
@@ -4,18 +4,27 @@ cimport numpy as np
cimport cython
from libc.math cimport exp
DTYPE = np.float
ctypedef np.float32_t DTYPE_t
cdef eps = 1.e-8
@cython.boundscheck(False)
cdef inline float patch_distance2d(DTYPE_t [:, :] p1,
cdef inline float patch_distance_2d(DTYPE_t [:, :] p1,
DTYPE_t [:, :] p2,
DTYPE_t [:, ::] w, int s):
cdef int i, j
cdef int center = s / 2
# Check if central pixel is too different in the 2 patches
cdef float tmp_diff = p1[center, center] - p2[center, center]
cdef float init = w[center, center] * tmp_diff * tmp_diff
if init > 1:
return eps
cdef float distance = 0
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:
return eps
for j in range(s):
tmp_diff = p1[i, j] - p2[i, j]
distance += w[i, j] * tmp_diff * tmp_diff
@@ -24,13 +33,37 @@ cdef inline float patch_distance2d(DTYPE_t [:, :] p1,
@cython.boundscheck(False)
cdef inline float patch_distance(DTYPE_t [:, :, :] p1,
cdef inline float patch_distance_2drgb(DTYPE_t [:, :, :] p1,
DTYPE_t [:, :, :] p2,
DTYPE_t [:, ::] w, int s):
cdef int i, j
cdef int center = s / 2
cdef int color
cdef float tmp_diff = 0
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:
return eps
for j in range(s):
for color in range(3):
tmp_diff = p1[i, j, color] - p2[i, j, color]
distance += w[i, j] * tmp_diff * tmp_diff
distance = exp(- distance)
return distance
@cython.boundscheck(False)
cdef inline float patch_distance_3d(DTYPE_t [:, :, :] p1,
DTYPE_t [:, :, :] p2,
DTYPE_t [:, :, ::] w, int s):
cdef int i, j, k
cdef float distance = 0
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:
return eps
for j in range(s):
for k in range(s):
tmp_diff = p1[i, j, k] - p2[i, j, k]
@@ -65,37 +98,129 @@ def _nl_means_denoising_2d(image, int s=7, int d=13, float h=0.1):
cdef int n_x, n_y
n_x, n_y = 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] result = padded.copy()
# We normalize by the image contrast, and divide by 3 because of 3 channels
h *= (np.max(padded) - np.min(padded)) / 3.
cdef float A = ((s - 1.) / 4.)
cdef float new_value
cdef float weight_sum, weight
xg, yg = np.mgrid[-offset:offset + 1, -offset:offset + 1]
cdef DTYPE_t [:, ::1] w = np.ascontiguousarray(np.exp(
- (xg ** 2 + yg ** 2) / (2 * A ** 2)).
astype(np.float32))
cdef float distance
cdef int x, y, i, j
cdef int x_start, x_end, y_start, y_end
cdef int x_start_i, x_end_i, y_start_j, y_end_j
w = 1. / (np.sum(w) * 2 * h ** 2) * w
# Coordinates of central pixel and patch bounds
for x in range(offset, n_x + offset):
x_start = x - offset
x_end = x + offset + 1
for y in range(offset, n_y + offset):
new_value = 0
weight_sum = 0
y_start = y - offset
y_end = y + offset + 1
# Coordinates of test pixel and patch bounds
for i in range(max(- d, offset - x),
min(d + 1, n_x - x - 1)):
x_start_i = x_start + i
x_end_i = x_end + i
for j in range(max(- d, offset - y),
min(d + 1, n_y - y - 1)):
y_start_j = y_start + j
y_end_j = y_end + j
weight = patch_distance_2d(
padded[x_start: x_end,
y_start: y_end],
padded[x_start_i: x_end_i,
y_start_j: y_end_j],
w, s)
weight_sum += weight
new_value += weight * padded[x + i, y + j]
result[x, y] = new_value / weight_sum
return result[offset:-offset, offset:-offset]
@cython.cdivision(True)
@cython.boundscheck(False)
def _nl_means_denoising_2drgb(image, int s=7, int d=13, float h=0.1):
"""
Perform non-local means denoising on 2-D RGB image
Parameters
----------
image: ndarray
input RGB image 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 n_x, n_y
n_x, n_y, _ = image.shape
cdef int offset = s / 2
cdef int x, y, i, j, color
cdef int x_start, x_end, y_start, y_end
cdef int x_start_i, x_end_i, y_start_j, y_end_j
cdef DTYPE_t [::1] new_values = np.zeros(3).astype(np.float32)
cdef DTYPE_t [:, :, ::1] padded = np.ascontiguousarray(util.pad(image,
((offset, offset), (offset, offset), (0, 0)),
mode='reflect').astype(np.float32))
cdef DTYPE_t [:, :, ::1] result = padded.copy()
h *= (np.max(padded) - np.min(padded))
cdef float A = ((s - 1.) / 4.)
cdef float new_value
cdef float weight_sum, weight
xg, yg = np.mgrid[-offset:offset + 1, -offset:offset + 1]
cdef DTYPE_t [:, ::1] w = np.ascontiguousarray(np.exp(
- (xg ** 2 + yg ** 2)/(2 * A ** 2)).
- (xg ** 2 + yg ** 2) / (2 * A ** 2)).
astype(np.float32))
cdef float distance
cdef int x, y, i, j
w = 1./ (np.sum(w) * 2 * h ** 2) * w
w = 1. / (np.sum(w) * 2 * h ** 2) * w
# Coordinates of central pixel and patch bounds
for x in range(offset, n_x + offset):
x_start = x - offset
x_end = x + offset + 1
for y in range(offset, n_y + offset):
new_value = 0
for color in range(3):
new_values[color] = 0
weight_sum = 0
y_start = y - offset
y_end = y + offset + 1
# Coordinates of test pixel and patch bounds
for i in range(max(- d, offset - x),
min(d + 1, n_x - x - 1)):
x_start_i = x_start + i
x_end_i = x_end + i
for j in range(max(- d, offset - y),
min(d + 1, n_y - y - 1)):
weight = patch_distance2d(
padded[x - offset: x + offset + 1,
y - offset: y + offset + 1],
padded[x + i - offset: x + i + offset + 1,
y + j - offset: y + j + offset + 1],
y_start_j = y_start + j
y_end_j = y_end + j
weight = patch_distance_2drgb(
padded[x_start: x_end,
y_start: y_end, :],
padded[x_start_i: x_end_i,
y_start_j: y_end_j, :],
w, s)
weight_sum += weight
new_value += weight * padded[x + i, y + j]
result[x, y] = new_value / weight_sum
for color in range(3):
new_values[color] += weight * padded[x + i, y + j,
color]
for color in range(3):
result[x, y, color] = new_values[color] / weight_sum
return result[offset:-offset, offset:-offset]
@@ -125,6 +250,7 @@ def _nl_means_denoising_3d(image, int s=7,
cdef int n_x, n_y, n_z
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] result = padded.copy()
@@ -132,101 +258,50 @@ def _nl_means_denoising_3d(image, int s=7,
cdef float A = ((s - 1.) / 4.)
cdef float new_value
cdef float weight_sum, weight
xg, yg, zg = np.mgrid[-offset: offset + 1, -offset: offset+1,
xg, yg, zg = np.mgrid[-offset: offset + 1, -offset: offset + 1,
-offset: offset + 1]
cdef DTYPE_t [:, :, ::1] w = np.ascontiguousarray(np.exp(
- (xg ** 2 + yg ** 2 + zg ** 2)/(2 * A ** 2)).
- (xg ** 2 + yg ** 2 + zg ** 2) / (2 * A ** 2)).
astype(np.float32))
cdef float distance
cdef int x, y, z, i, j, k
w = 1./ (np.sum(w) * 2 * h ** 2) * w
cdef int x_start, x_end, y_start, y_end, z_start, z_end
cdef int x_start_i, x_end_i, y_start_j, y_end_j, z_start_k, z_end_k
w = 1. / (np.sum(w) * 2 * h ** 2) * w
# Coordinates of central pixel and patch bounds
for x in range(offset, n_x + offset):
x_start = x - offset
x_end = x + offset + 1
for y in range(offset, n_y + offset):
y_start = y - offset
y_end = y + offset + 1
for z in range(offset, n_z + offset):
z_start = z - offset
z_end = z + offset + 1
new_value = 0
weight_sum = 0
# Coordinates of test pixel and patch bounds
for i in range(max(- d, offset - x),
min(d + 1, n_x - x - 1)):
x_start_i = x_start + i
x_end_i = x_end + i
for j in range(max(- d, offset - y),
min(d + 1, n_y - y - 1)):
y_start_j = y_start + j
y_end_j = y_end + j
for k in range(max(- d, offset - z),
min(d + 1, n_z - z - 1)):
weight = patch_distance(
padded[x - offset: x + offset +1,
y - offset: y + offset +1,
z - offset: z + offset +1],
padded[x + i - offset: x + i + offset +1,
y + j - offset: y + j + offset +1,
z + k - offset: z + k + offset +1],
z_start_k = z_start + k
z_end_k = z_end + k
weight = patch_distance_3d(
padded[x_start: x_end,
y_start: y_end,
z_start: z_end],
padded[x_start_i: x_end_i,
y_start_j: y_end_j,
z_start_k: z_end_k],
w, s)
weight_sum += weight
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]
def nl_means_denoising(image, patch_size=7, patch_distance=11, h=0.1):
"""
Perform non-local means denoising on 2-D or 3-D grayscale arrays
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.
Returns
-------
result: ndarray
denoised image, of same shape as `image`.
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_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, patch_size,
patch_distance, 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))
else:
raise ValueError("Non local means denoising is only possible for \
2D and 3-D grayscale images.")
skimage/filter/nl_means_denoising.py
+74
View File
@@ -0,0 +1,74 @@
import numpy as np
from _nl_means_denoising import _nl_means_denoising_2d, \
_nl_means_denoising_2drgb, _nl_means_denoising_3d
def nl_means_denoising(image, patch_size=7, patch_distance=11, h=0.1):
"""
Perform non-local means denoising on 2-D or 3-D grayscale images, and
2-D RGB 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`.
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_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, patch_size,
patch_distance, 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.")
skimage/filter/tests/test_denoise.py
+169
View File
@@ -0,0 +1,169 @@
import numpy as np
from numpy.testing import run_module_suite, assert_raises, assert_equal
from skimage import filter, data, color, img_as_float
np.random.seed(1234)
lena = img_as_float(data.lena()[:256, :256])
lena_gray = color.rgb2gray(lena)
def test_denoise_tv_chambolle_2d():
# lena image
img = lena_gray
# add noise to lena
img += 0.5 * img.std() * np.random.random(img.shape)
# clip noise so that it does not exceed allowed range for float images.
img = np.clip(img, 0, 1)
# denoise
denoised_lena = filter.denoise_tv_chambolle(img, weight=60.0)
# which dtype?
assert denoised_lena.dtype in [np.float, np.float32, np.float64]
from scipy import ndimage
grad = ndimage.morphological_gradient(img, size=((3, 3)))
grad_denoised = ndimage.morphological_gradient(
denoised_lena, size=((3, 3)))
# test if the total variation has decreased
assert grad_denoised.dtype == np.float
assert (np.sqrt((grad_denoised**2).sum())
< np.sqrt((grad**2).sum()) / 2)
def test_denoise_tv_chambolle_multichannel():
denoised0 = filter.denoise_tv_chambolle(lena[..., 0], weight=60.0)
denoised = filter.denoise_tv_chambolle(lena, weight=60.0, multichannel=True)
assert_equal(denoised[..., 0], denoised0)
def test_denoise_tv_chambolle_float_result_range():
# lena image
img = lena_gray
int_lena = np.multiply(img, 255).astype(np.uint8)
assert np.max(int_lena) > 1
denoised_int_lena = filter.denoise_tv_chambolle(int_lena, weight=60.0)
# test if the value range of output float data is within [0.0:1.0]
assert denoised_int_lena.dtype == np.float
assert np.max(denoised_int_lena) <= 1.0
assert np.min(denoised_int_lena) >= 0.0
def test_denoise_tv_chambolle_3d():
"""Apply the TV denoising algorithm on a 3D image representing a sphere."""
x, y, z = np.ogrid[0:40, 0:40, 0:40]
mask = (x - 22)**2 + (y - 20)**2 + (z - 17)**2 < 8**2
mask = 100 * mask.astype(np.float)
mask += 60
mask += 20 * np.random.random(mask.shape)
mask[mask < 0] = 0
mask[mask > 255] = 255
res = filter.denoise_tv_chambolle(mask.astype(np.uint8), weight=100)
assert res.dtype == np.float
assert res.std() * 255 < mask.std()
# test wrong number of dimensions
assert_raises(ValueError, filter.denoise_tv_chambolle,
np.random.random((8, 8, 8, 8)))
def test_denoise_tv_bregman_2d():
img = lena_gray
# add some random noise
img += 0.5 * img.std() * np.random.random(img.shape)
img = np.clip(img, 0, 1)
out1 = filter.denoise_tv_bregman(img, weight=10)
out2 = filter.denoise_tv_bregman(img, weight=5)
# make sure noise is reduced
assert img.std() > out1.std()
assert out1.std() > out2.std()
def test_denoise_tv_bregman_float_result_range():
# lena image
img = lena_gray
int_lena = np.multiply(img, 255).astype(np.uint8)
assert np.max(int_lena) > 1
denoised_int_lena = filter.denoise_tv_bregman(int_lena, weight=60.0)
# test if the value range of output float data is within [0.0:1.0]
assert denoised_int_lena.dtype == np.float
assert np.max(denoised_int_lena) <= 1.0
assert np.min(denoised_int_lena) >= 0.0
def test_denoise_tv_bregman_3d():
img = lena
# add some random noise
img += 0.5 * img.std() * np.random.random(img.shape)
img = np.clip(img, 0, 1)
out1 = filter.denoise_tv_bregman(img, weight=10)
out2 = filter.denoise_tv_bregman(img, weight=5)
# make sure noise is reduced
assert img.std() > out1.std()
assert out1.std() > out2.std()
def test_denoise_bilateral_2d():
img = lena_gray
# add some random noise
img += 0.5 * img.std() * np.random.random(img.shape)
img = np.clip(img, 0, 1)
out1 = filter.denoise_bilateral(img, sigma_range=0.1, sigma_spatial=20)
out2 = filter.denoise_bilateral(img, sigma_range=0.2, sigma_spatial=30)
# make sure noise is reduced
assert img.std() > out1.std()
assert out1.std() > out2.std()
def test_denoise_bilateral_3d():
img = lena
# add some random noise
img += 0.5 * img.std() * np.random.random(img.shape)
img = np.clip(img, 0, 1)
out1 = filter.denoise_bilateral(img, sigma_range=0.1, sigma_spatial=20)
out2 = filter.denoise_bilateral(img, sigma_range=0.2, sigma_spatial=30)
# make sure noise is reduced
assert img.std() > out1.std()
assert out1.std() > out2.std()
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 = filter.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 = lena[-100:, -100:]
# add some random noise
img += 0.5 * img.std() * np.random.random(img.shape)
img = np.clip(img, 0, 1)
denoised = filter.nl_means_denoising(img, 7, 9, 0.08)
# make sure noise is reduced
assert img.std() > denoised.std()
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 = filter.nl_means_denoising(img, 5, 4, 0.1)
# make sure noise is reduced
assert img.std() > denoised.std()
if __name__ == "__main__":
run_module_suite()