mirror of
https://github.com/wassname/scikit-image.git
synced 2026-07-04 07:01:16 +08:00
2-D and 3-D implementation of NL-means denoising
This commit is contained in:
Emmanuelle Gouillart
committed by
emmanuelle
emmanuelle
parent
f82243904f
commit
a508ec54ca
@@ -0,0 +1,232 @@
|
||||
import numpy as np
|
||||
from skimage import util
|
||||
cimport numpy as np
|
||||
cimport cython
|
||||
from libc.math cimport exp
|
||||
|
||||
DTYPE = np.float
|
||||
ctypedef np.float32_t DTYPE_t
|
||||
|
||||
|
||||
@cython.boundscheck(False)
|
||||
cdef inline float patch_distance2d(DTYPE_t [:, :] p1,
|
||||
DTYPE_t [:, :] p2,
|
||||
DTYPE_t [:, ::] w, int s):
|
||||
cdef int i, j
|
||||
cdef float distance = 0
|
||||
cdef float tmp_diff
|
||||
for i in range(s):
|
||||
for j in range(s):
|
||||
tmp_diff = p1[i, j] - p2[i, j]
|
||||
distance += w[i, j] * tmp_diff * tmp_diff
|
||||
distance = exp(- distance)
|
||||
return distance
|
||||
|
||||
|
||||
@cython.boundscheck(False)
|
||||
cdef inline float patch_distance(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):
|
||||
for j in range(s):
|
||||
for k in range(s):
|
||||
tmp_diff = p1[i, j, k] - p2[i, j, k]
|
||||
distance += w[i, j, k] * tmp_diff * tmp_diff
|
||||
distance = exp(- distance)
|
||||
return distance
|
||||
|
||||
|
||||
@cython.cdivision(True)
|
||||
@cython.boundscheck(False)
|
||||
def _nl_means_denoising_2d(image, int s=7, int d=13, float h=0.1):
|
||||
"""
|
||||
Perform non-local means denoising on 2-D array
|
||||
|
||||
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 n_x, n_y
|
||||
n_x, n_y = image.shape
|
||||
cdef int offset = s / 2
|
||||
cdef DTYPE_t [:, ::1] padded = np.ascontiguousarray(util.pad(image,
|
||||
offset, 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)).
|
||||
astype(np.float32))
|
||||
cdef float distance
|
||||
cdef int x, y, i, j
|
||||
w = 1./ (np.sum(w) * 2 * h ** 2) * w
|
||||
for x in range(offset, n_x + offset):
|
||||
for y in range(offset, n_y + offset):
|
||||
new_value = 0
|
||||
weight_sum = 0
|
||||
for i in range(max(- d, offset - x),
|
||||
min(d + 1, n_x - x - 1)):
|
||||
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],
|
||||
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_3d(image, int s=7,
|
||||
int d=13, float h=0.1):
|
||||
"""
|
||||
Perform non-local means denoising on 3-D array
|
||||
|
||||
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)
|
||||
"""
|
||||
if s % 2 == 0:
|
||||
s += 1 # odd value for symmetric patch
|
||||
cdef int n_x, n_y, n_z
|
||||
n_x, n_y, n_z = image.shape
|
||||
cdef int offset = s / 2
|
||||
cdef DTYPE_t [:, :, ::1] padded = np.ascontiguousarray(util.pad(image,
|
||||
offset, 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, 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)).
|
||||
astype(np.float32))
|
||||
cdef float distance
|
||||
cdef int x, y, z, i, j, k
|
||||
w = 1./ (np.sum(w) * 2 * h ** 2) * w
|
||||
for x in range(offset, n_x + offset):
|
||||
for y in range(offset, n_y + offset):
|
||||
for z in range(offset, n_z + offset):
|
||||
new_value = 0
|
||||
weight_sum = 0
|
||||
for i in range(max(- d, offset - x),
|
||||
min(d + 1, n_x - x - 1)):
|
||||
for j in range(max(- d, offset - y),
|
||||
min(d + 1, n_y - y - 1)):
|
||||
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],
|
||||
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.")
|
||||
@@ -14,6 +14,7 @@ def configuration(parent_package='', top_path=None):
|
||||
config.add_data_dir('rank/tests')
|
||||
|
||||
cython(['_ctmf.pyx'], working_path=base_path)
|
||||
cython(['_nl_means_denoising.pyx'], working_path=base_path)
|
||||
cython(['rank/core_cy.pyx'], working_path=base_path)
|
||||
cython(['rank/generic_cy.pyx'], working_path=base_path)
|
||||
cython(['rank/percentile_cy.pyx'], working_path=base_path)
|
||||
@@ -21,6 +22,9 @@ def configuration(parent_package='', top_path=None):
|
||||
|
||||
config.add_extension('_ctmf', sources=['_ctmf.c'],
|
||||
include_dirs=[get_numpy_include_dirs()])
|
||||
config.add_extension('_nl_means_denoising',
|
||||
sources=['_nl_means_denoising.c'],
|
||||
include_dirs=[get_numpy_include_dirs()])
|
||||
config.add_extension('rank.core_cy', sources=['rank/core_cy.c'],
|
||||
include_dirs=[get_numpy_include_dirs()])
|
||||
config.add_extension('rank.generic_cy', sources=['rank/generic_cy.c'],
|
||||
|
||||
Reference in New Issue
Block a user