mirror of
https://github.com/wassname/scikit-image.git
synced 2026-06-27 19:48:43 +08:00
Merge pull request #874 from emmanuelle/nlm_denoise
FEAT: NL-means denoising
This commit is contained in:
Josh Warner
@@ -149,6 +149,9 @@ Library:
|
||||
Extension: skimage.restoration._denoise_cy
|
||||
Sources:
|
||||
skimage/restoration/_denoise_cy.pyx
|
||||
Extension: skimage.restoration._nl_means_denoising
|
||||
Sources:
|
||||
skimage/restoration/_nl_means_denoising.pyx
|
||||
Extension: skimage.feature._hessian_det_appx
|
||||
Sources:
|
||||
skimage/exposure/_hessian_det_appx.pyx
|
||||
|
||||
@@ -0,0 +1,41 @@
|
||||
"""
|
||||
=================================================
|
||||
Non-local means denoising for preserving textures
|
||||
=================================================
|
||||
|
||||
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
|
||||
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.restoration import nl_means_denoising
|
||||
|
||||
|
||||
astro = img_as_float(data.astronaut())
|
||||
astro = astro[30:180, 150:300]
|
||||
|
||||
noisy = astro + 0.3 * np.random.random(astro.shape)
|
||||
noisy = np.clip(noisy, 0, 1)
|
||||
|
||||
denoise = nl_means_denoising(noisy, 7, 9, 0.08)
|
||||
|
||||
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()
|
||||
@@ -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
|
||||
|
||||
__all__ = ['wiener',
|
||||
'unsupervised_wiener',
|
||||
@@ -30,4 +30,5 @@ __all__ = ['wiener',
|
||||
'unwrap_phase',
|
||||
'denoise_tv_bregman',
|
||||
'denoise_tv_chambolle',
|
||||
'denoise_bilateral']
|
||||
'denoise_bilateral',
|
||||
'nl_means_denoising']
|
||||
|
||||
@@ -0,0 +1,712 @@
|
||||
import numpy as np
|
||||
from skimage import util
|
||||
cimport numpy as np
|
||||
cimport cython
|
||||
from libc.math cimport exp
|
||||
|
||||
ctypedef np.float32_t IMGDTYPE
|
||||
|
||||
cdef float DISTANCE_CUTOFF = 5.
|
||||
|
||||
@cython.boundscheck(False)
|
||||
cdef inline float patch_distance_2d(IMGDTYPE [:, :] p1,
|
||||
IMGDTYPE [:, :] p2,
|
||||
IMGDTYPE [:, ::] w, int s):
|
||||
"""
|
||||
Compute a Gaussian distance between two image patches.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
p1 : 2-D array_like
|
||||
First patch.
|
||||
p2 : 2-D array_like
|
||||
Second patch.
|
||||
w : 2-D array_like
|
||||
Array of weigths for the different pixels of the patches.
|
||||
s : int
|
||||
Linear size of the patches.
|
||||
|
||||
Returns
|
||||
-------
|
||||
distance : float
|
||||
Gaussian distance between the two patches
|
||||
|
||||
Notes
|
||||
-----
|
||||
The returned distance is given by
|
||||
|
||||
.. math:: \exp( -w (p1 - p2)^2)
|
||||
"""
|
||||
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 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 > DISTANCE_CUTOFF:
|
||||
return 0.
|
||||
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_2drgb(IMGDTYPE [:, :, :] p1,
|
||||
IMGDTYPE [:, :, :] p2,
|
||||
IMGDTYPE [:, ::] w, int s):
|
||||
"""
|
||||
Compute a Gaussian distance between two image patches.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
p1 : 3-D array_like
|
||||
First patch, 2D image with last dimension corresponding to channels.
|
||||
p2 : 3-D array_like
|
||||
Second patch, 2D image with last dimension corresponding to channels.
|
||||
w : 2-D array_like
|
||||
Array of weights for the different pixels of the patches.
|
||||
s : int
|
||||
Linear size of the patches.
|
||||
|
||||
Returns
|
||||
-------
|
||||
distance : float
|
||||
Gaussian distance between the two patches
|
||||
|
||||
Notes
|
||||
-----
|
||||
The returned distance is given by
|
||||
|
||||
.. math:: \exp( -w (p1 - p2)^2)
|
||||
"""
|
||||
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 > DISTANCE_CUTOFF:
|
||||
return 0.
|
||||
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(IMGDTYPE [:, :, :] p1,
|
||||
IMGDTYPE [:, :, :] p2,
|
||||
IMGDTYPE [:, :, ::] w, int s):
|
||||
"""
|
||||
Compute a Gaussian distance between two image patches.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
p1 : 3-D array_like
|
||||
First patch.
|
||||
p2 : 3-D array_like
|
||||
Second patch.
|
||||
w : 3-D array_like
|
||||
Array of weights for the different pixels of the patches.
|
||||
s : int
|
||||
Linear size of the patches.
|
||||
|
||||
Returns
|
||||
-------
|
||||
distance : float
|
||||
Gaussian distance between the two patches
|
||||
|
||||
Notes
|
||||
-----
|
||||
The returned distance is given by
|
||||
|
||||
.. math:: \exp( -w (p1 - p2)^2)
|
||||
"""
|
||||
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 > DISTANCE_CUTOFF:
|
||||
return 0.
|
||||
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 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.
|
||||
|
||||
Returns
|
||||
-------
|
||||
result : ndarray
|
||||
Denoised image, of same shape as input image.
|
||||
"""
|
||||
if s % 2 == 0:
|
||||
s += 1 # odd value for symmetric patch
|
||||
cdef int n_row, n_col, n_ch
|
||||
n_row, n_col, n_ch = image.shape
|
||||
cdef int offset = s / 2
|
||||
cdef int row, col, i, j, color
|
||||
cdef int row_start, row_end, col_start, col_end
|
||||
cdef int row_start_i, row_end_i, col_start_j, col_end_j
|
||||
cdef IMGDTYPE [::1] new_values = np.zeros(n_ch).astype(np.float32)
|
||||
cdef IMGDTYPE [:, :, ::1] padded = np.ascontiguousarray(util.pad(image,
|
||||
((offset, offset), (offset, offset), (0, 0)),
|
||||
mode='reflect').astype(np.float32))
|
||||
cdef IMGDTYPE [:, :, ::1] result = padded.copy()
|
||||
cdef float A = ((s - 1.) / 4.)
|
||||
cdef float new_value
|
||||
cdef float weight_sum, weight
|
||||
xg_row, xg_col = np.mgrid[-offset:offset + 1, -offset:offset + 1]
|
||||
cdef IMGDTYPE [:, ::1] w = np.ascontiguousarray(np.exp(
|
||||
-(xg_row ** 2 + xg_col ** 2) / (2 * A ** 2)).
|
||||
astype(np.float32))
|
||||
cdef float distance
|
||||
w = 1. / (n_ch * np.sum(w) * h ** 2) * w
|
||||
|
||||
# Coordinates of central pixel
|
||||
# Iterate over rows, taking padding into account
|
||||
for row in range(offset, n_row + offset):
|
||||
row_start = row - offset
|
||||
row_end = row + offset + 1
|
||||
# Iterate over columns, taking padding into account
|
||||
for col in range(offset, n_col + offset):
|
||||
# Initialize per-channel bins
|
||||
for color in range(n_ch):
|
||||
new_values[color] = 0
|
||||
# Reset weights for each local region
|
||||
weight_sum = 0
|
||||
col_start = col - offset
|
||||
col_end = col + offset + 1
|
||||
|
||||
# Iterate over local 2d patch for each pixel
|
||||
# First rows
|
||||
for i in range(max(-d, offset - row),
|
||||
min(d + 1, n_row + offset - row)):
|
||||
row_start_i = row_start + i
|
||||
row_end_i = row_end + i
|
||||
# Local patch columns
|
||||
for j in range(max(-d, offset - col),
|
||||
min(d + 1, n_col + offset - col)):
|
||||
col_start_j = col_start + j
|
||||
col_end_j = col_end + j
|
||||
# Shortcut for grayscale, else assume RGB
|
||||
if n_ch == 1:
|
||||
weight = patch_distance_2d(
|
||||
padded[row_start:row_end,
|
||||
col_start:col_end, 0],
|
||||
padded[row_start_i:row_end_i,
|
||||
col_start_j:col_end_j, 0],
|
||||
w, s)
|
||||
else:
|
||||
weight = patch_distance_2drgb(
|
||||
padded[row_start:row_end,
|
||||
col_start:col_end, :],
|
||||
padded[row_start_i:row_end_i,
|
||||
col_start_j:col_end_j, :],
|
||||
w, s)
|
||||
|
||||
# Collect results in weight sum
|
||||
weight_sum += weight
|
||||
# Apply to each channel multiplicatively
|
||||
for color in range(n_ch):
|
||||
new_values[color] += weight * padded[row + i,
|
||||
col + j, color]
|
||||
|
||||
# Normalize the result
|
||||
for color in range(n_ch):
|
||||
result[row, col, color] = new_values[color] / weight_sum
|
||||
|
||||
# Return cropped result, undoing padding
|
||||
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).
|
||||
|
||||
Returns
|
||||
-------
|
||||
result : ndarray
|
||||
Denoised image, of same shape as input image.
|
||||
"""
|
||||
if s % 2 == 0:
|
||||
s += 1 # odd value for symmetric patch
|
||||
cdef int n_pln, n_row, n_col
|
||||
n_pln, n_row, n_col = image.shape
|
||||
cdef int offset = s / 2
|
||||
# padd the image so that boundaries are denoised as well
|
||||
cdef IMGDTYPE [:, :, ::1] padded = np.ascontiguousarray(util.pad(
|
||||
image.astype(np.float32),
|
||||
offset, mode='reflect'))
|
||||
cdef IMGDTYPE [:, :, ::1] result = padded.copy()
|
||||
cdef float A = ((s - 1.) / 4.)
|
||||
cdef float new_value
|
||||
cdef float weight_sum, weight
|
||||
xg_pln, xg_row, xg_col = np.mgrid[-offset: offset + 1,
|
||||
-offset: offset + 1,
|
||||
-offset: offset + 1]
|
||||
cdef IMGDTYPE [:, :, ::1] w = np.ascontiguousarray(np.exp(
|
||||
-(xg_pln ** 2 + xg_row ** 2 + xg_col ** 2) /
|
||||
(2 * A ** 2)).astype(np.float32))
|
||||
cdef float distance
|
||||
cdef int pln, row, col, i, j, k
|
||||
cdef int pln_start, pln_end, row_start, row_end, col_start, col_end
|
||||
cdef int pln_start_i, pln_end_i, row_start_j, row_end_j, \
|
||||
col_start_k, col_end_k
|
||||
w = 1. / (np.sum(w) * h ** 2) * w
|
||||
|
||||
# Coordinates of central pixel
|
||||
# Iterate over planes, taking padding into account
|
||||
for pln in range(offset, n_pln + offset):
|
||||
pln_start = pln - offset
|
||||
pln_end = pln + offset + 1
|
||||
# Iterate over rows, taking padding into account
|
||||
for row in range(offset, n_row + offset):
|
||||
row_start = row - offset
|
||||
row_end = row + offset + 1
|
||||
# Iterate over columns, taking padding into account
|
||||
for col in range(offset, n_col + offset):
|
||||
col_start = col - offset
|
||||
col_end = col + offset + 1
|
||||
new_value = 0
|
||||
weight_sum = 0
|
||||
|
||||
# Iterate over local 3d patch for each pixel
|
||||
# First planes
|
||||
for i in range(max(-d, offset - pln),
|
||||
min(d + 1, n_pln + offset - pln)):
|
||||
pln_start_i = pln_start + i
|
||||
pln_end_i = pln_end + i
|
||||
# Rows
|
||||
for j in range(max(-d, offset - row),
|
||||
min(d + 1, n_row + offset - row)):
|
||||
row_start_j = row_start + j
|
||||
row_end_j = row_end + j
|
||||
# Columns
|
||||
for k in range(max(-d, offset - col),
|
||||
min(d + 1, n_col + offset - col)):
|
||||
col_start_k = col_start + k
|
||||
col_end_k = col_end + k
|
||||
weight = patch_distance_3d(
|
||||
padded[pln_start:pln_end,
|
||||
row_start:row_end,
|
||||
col_start:col_end],
|
||||
padded[pln_start_i:pln_end_i,
|
||||
row_start_j:row_end_j,
|
||||
col_start_k:col_end_k],
|
||||
w, s)
|
||||
# Collect results in weight sum
|
||||
weight_sum += weight
|
||||
new_value += weight * padded[pln + i,
|
||||
row + j, col + k]
|
||||
|
||||
# Normalize the result
|
||||
result[pln, row, col] = new_value / weight_sum
|
||||
|
||||
# Return cropped result, undoing padding
|
||||
return result[offset:-offset, offset:-offset, offset:-offset]
|
||||
|
||||
#-------------- Accelerated algorithm of Froment 2015 ------------------
|
||||
|
||||
|
||||
@cython.cdivision(True)
|
||||
@cython.boundscheck(False)
|
||||
cdef inline float _integral_to_distance_2d(IMGDTYPE [:, ::] integral,
|
||||
int row, int col, int offset, float h2s2):
|
||||
"""
|
||||
References
|
||||
----------
|
||||
Jacques Froment. Parameter-Free Fast Pixelwise Non-Local Means
|
||||
Denoising. Image Processing On Line, 2014, vol. 4, p. 300-326.
|
||||
|
||||
Used in _fast_nl_means_denoising_2d
|
||||
"""
|
||||
cdef float distance
|
||||
distance = integral[row + offset, col + offset] + \
|
||||
integral[row - offset, col - offset] - \
|
||||
integral[row - offset, col + offset] - \
|
||||
integral[row + offset, col - offset]
|
||||
distance /= h2s2
|
||||
return distance
|
||||
|
||||
|
||||
@cython.cdivision(True)
|
||||
@cython.boundscheck(False)
|
||||
cdef inline float _integral_to_distance_3d(IMGDTYPE [:, :, ::] integral,
|
||||
int pln, int row, int col, int offset,
|
||||
float s_cube_h_square):
|
||||
"""
|
||||
References
|
||||
----------
|
||||
Jacques Froment. Parameter-Free Fast Pixelwise Non-Local Means
|
||||
Denoising. Image Processing On Line, 2014, vol. 4, p. 300-326.
|
||||
|
||||
Used in _fast_nl_means_denoising_3d
|
||||
"""
|
||||
cdef float distance
|
||||
distance = (integral[pln + offset, row + offset, col + offset] -
|
||||
integral[pln - offset, row - offset, col - offset] +
|
||||
integral[pln - offset, row - offset, col + offset] +
|
||||
integral[pln - offset, row + offset, col - offset] +
|
||||
integral[pln + offset, row - offset, col - offset] -
|
||||
integral[pln - offset, row + offset, col + offset] -
|
||||
integral[pln + offset, row - offset, col + offset] -
|
||||
integral[pln + offset, row + offset, col - offset])
|
||||
distance /= s_cube_h_square
|
||||
return distance
|
||||
|
||||
|
||||
@cython.cdivision(True)
|
||||
@cython.boundscheck(False)
|
||||
cdef inline _integral_image_2d(IMGDTYPE [:, :, ::] padded,
|
||||
IMGDTYPE [:, ::] integral, int t_row,
|
||||
int t_col, int n_row, int n_col, int n_ch):
|
||||
"""
|
||||
Computes the integral of the squared difference between an image ``padded``
|
||||
and the same image shifted by ``(t_row, t_col)``.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
padded : ndarray of shape (n_row, n_col, n_ch)
|
||||
Image of interest.
|
||||
integral : ndarray
|
||||
Output of the function. The array is filled with integral values.
|
||||
``integral`` should have the same shape as ``padded``.
|
||||
t_row : int
|
||||
Shift along the row axis.
|
||||
t_col : int
|
||||
Shift along the column axis.
|
||||
n_row : int
|
||||
n_col : int
|
||||
n_ch : int
|
||||
|
||||
Notes
|
||||
-----
|
||||
|
||||
The integral computation could be performed using
|
||||
``transform.integral_image``, but this helper function saves memory
|
||||
by avoiding copies of ``padded``.
|
||||
"""
|
||||
cdef int row, col
|
||||
cdef float distance
|
||||
for row in range(max(1, -t_row), min(n_row, n_row - t_row)):
|
||||
for col in range(max(1, -t_col), min(n_col, n_col - t_col)):
|
||||
if n_ch == 1:
|
||||
distance = (padded[row, col, 0] -
|
||||
padded[row + t_row, col + t_col, 0])**2
|
||||
else:
|
||||
distance = ((padded[row, col, 0] -
|
||||
padded[row + t_row, col + t_col, 0])**2 +
|
||||
(padded[row, col, 1] -
|
||||
padded[row + t_row, col + t_col, 1])**2 +
|
||||
(padded[row, col, 2] -
|
||||
padded[row + t_row, col + t_col, 2])**2)
|
||||
integral[row, col] = distance + \
|
||||
integral[row - 1, col] + \
|
||||
integral[row, col - 1] - \
|
||||
integral[row - 1, col - 1]
|
||||
|
||||
@cython.cdivision(True)
|
||||
@cython.boundscheck(False)
|
||||
cdef inline _integral_image_3d(IMGDTYPE [:, :, ::] padded,
|
||||
IMGDTYPE [:, :, ::] integral, int t_pln,
|
||||
int t_row, int t_col, int n_pln, int n_row,
|
||||
int n_col):
|
||||
"""
|
||||
Computes the integral of the squared difference between an image ``padded``
|
||||
and the same image shifted by ``(t_pln, t_row, t_col)``.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
padded : ndarray of shape (n_pln, n_row, n_col)
|
||||
Image of interest.
|
||||
integral : ndarray
|
||||
Output of the function. The array is filled with integral values.
|
||||
``integral`` should have the same shape as ``padded``.
|
||||
t_pln : int
|
||||
Shift along the plane axis.
|
||||
t_row : int
|
||||
Shift along the row axis.
|
||||
t_col : int
|
||||
Shift along the column axis.
|
||||
n_pln : int
|
||||
n_row : int
|
||||
n_col : int
|
||||
|
||||
Notes
|
||||
-----
|
||||
|
||||
The integral computation could be performed using
|
||||
``transform.integral_image``, but this helper function saves memory
|
||||
by avoiding copies of ``padded``.
|
||||
"""
|
||||
cdef int pln, row, col
|
||||
cdef float distance
|
||||
for pln in range(max(1, -t_pln), min(n_pln, n_pln - t_pln)):
|
||||
for row in range(max(1, -t_row), min(n_row, n_row - t_row)):
|
||||
for col in range(max(1, -t_col), min(n_col, n_col - t_col)):
|
||||
integral[pln, row, col] = \
|
||||
((padded[pln, row, col] -
|
||||
padded[pln + t_pln, row + t_row, col + t_col])**2 +
|
||||
integral[pln - 1, row, col] +
|
||||
integral[pln, row - 1, col] +
|
||||
integral[pln, row, col - 1] +
|
||||
integral[pln - 1, row - 1, col - 1] -
|
||||
integral[pln - 1, row - 1, col] -
|
||||
integral[pln, row - 1, col - 1] -
|
||||
integral[pln - 1, row, col - 1])
|
||||
|
||||
|
||||
@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
|
||||
2-D input data to be denoised, grayscale or RGB.
|
||||
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.
|
||||
|
||||
Returns
|
||||
-------
|
||||
result : ndarray
|
||||
Denoised image, of same shape as input image.
|
||||
"""
|
||||
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 IMGDTYPE [:, :, ::1] padded = np.ascontiguousarray(util.pad(image,
|
||||
((pad_size, pad_size), (pad_size, pad_size), (0, 0)),
|
||||
mode='reflect').astype(np.float32))
|
||||
cdef IMGDTYPE [:, :, ::1] result = np.zeros_like(padded)
|
||||
cdef IMGDTYPE [:, ::1] weights = np.zeros_like(padded[..., 0], order='C')
|
||||
cdef IMGDTYPE [:, ::1] integral = np.zeros_like(padded[..., 0], order='C')
|
||||
cdef int n_row, n_col, n_ch, t_row, t_col, row, col
|
||||
cdef float weight, distance
|
||||
cdef float alpha
|
||||
cdef float h2 = h ** 2.
|
||||
cdef float s2 = s ** 2.
|
||||
n_row, n_col, n_ch = image.shape
|
||||
cdef float h2s2 = n_ch * h2 * s2
|
||||
n_row += 2 * pad_size
|
||||
n_col += 2 * pad_size
|
||||
|
||||
# Outer loops on patch shifts
|
||||
# With t2 >= 0, reference patch is always on the left of test patch
|
||||
# Iterate over shifts along the row axis
|
||||
for t_row in range(-d, d + 1):
|
||||
# Iterate over shifts along the column axis
|
||||
for t_col in range(0, d + 1):
|
||||
# alpha is to account for patches on the same column
|
||||
# distance is computed twice in this case
|
||||
if t_col == 0 and t_row is not 0:
|
||||
alpha = 0.5
|
||||
else:
|
||||
alpha = 1.
|
||||
# Compute integral image of the squared difference between
|
||||
# padded and the same image shifted by (t_row, t_col)
|
||||
integral = np.zeros_like(padded[..., 0], order='C')
|
||||
_integral_image_2d(padded, integral, t_row, t_col,
|
||||
n_row, n_col, n_ch)
|
||||
|
||||
# Inner loops on pixel coordinates
|
||||
# Iterate over rows, taking offset and shift into account
|
||||
for row in range(max(offset, offset - t_row),
|
||||
min(n_row - offset, n_row - offset - t_row)):
|
||||
# Iterate over columns, taking offset and shift into account
|
||||
for col in range(max(offset, offset - t_col),
|
||||
min(n_col - offset, n_col - offset - t_col)):
|
||||
# Compute squared distance between shifted patches
|
||||
distance = _integral_to_distance_2d(integral, row, col,
|
||||
offset, h2s2)
|
||||
# exp of large negative numbers is close to zero
|
||||
if distance > DISTANCE_CUTOFF:
|
||||
continue
|
||||
weight = alpha * exp(-distance)
|
||||
# Accumulate weights corresponding to different shifts
|
||||
weights[row, col] += weight
|
||||
weights[row + t_row, col + t_col] += weight
|
||||
# Iterate over channels
|
||||
for ch in range(n_ch):
|
||||
result[row, col, ch] += weight * \
|
||||
padded[row + t_row, col + t_col, ch]
|
||||
result[row + t_row, col + t_col, ch] += \
|
||||
weight * padded[row, col, ch]
|
||||
|
||||
# Normalize pixel values using sum of weights of contributing patches
|
||||
for row in range(offset, n_row - offset):
|
||||
for col in range(offset, n_col - offset):
|
||||
for channel in range(n_ch):
|
||||
# No risk of division by zero, since the contribution
|
||||
# of a null shift is strictly positive
|
||||
result[row, col, channel] /= weights[row, col]
|
||||
|
||||
# Return cropped result, undoing padding
|
||||
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.
|
||||
|
||||
Returns
|
||||
-------
|
||||
result : ndarray
|
||||
Denoised image, of same shape as input image.
|
||||
"""
|
||||
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 IMGDTYPE [:, :, ::1] padded = np.ascontiguousarray(util.pad(image,
|
||||
pad_size, mode='reflect').astype(np.float32))
|
||||
cdef IMGDTYPE [:, :, ::1] result = np.zeros_like(padded)
|
||||
cdef IMGDTYPE [:, :, ::1] weights = np.zeros_like(padded)
|
||||
cdef IMGDTYPE [:, :, ::1] integral = np.zeros_like(padded)
|
||||
cdef int n_pln, n_row, n_col, t_pln, t_row, t_col, \
|
||||
pln, row, col
|
||||
cdef int pln_dist_min, pln_dist_max, row_dist_min, row_dist_max, \
|
||||
col_dist_min, col_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_pln, n_row, n_col = image.shape
|
||||
n_pln += 2 * pad_size
|
||||
n_row += 2 * pad_size
|
||||
n_col += 2 * pad_size
|
||||
|
||||
# Outer loops on patch shifts
|
||||
# With t2 >= 0, reference patch is always on the left of test patch
|
||||
# Iterate over shifts along the plane axis
|
||||
for t_pln in range(-d, d + 1):
|
||||
pln_dist_min = max(offset, offset - t_pln)
|
||||
pln_dist_max = min(n_pln - offset, n_pln - offset - t_pln)
|
||||
# Iterate over shifts along the row axis
|
||||
for t_row in range(-d, d + 1):
|
||||
row_dist_min = max(offset, offset - t_row)
|
||||
row_dist_max = min(n_row - offset, n_row - offset - t_row)
|
||||
# Iterate over shifts along the column axis
|
||||
for t_col in range(0, d + 1):
|
||||
col_dist_min = max(offset, offset - t_col)
|
||||
col_dist_max = min(n_col - offset, n_col - offset - t_col)
|
||||
# alpha is to account for patches on the same column
|
||||
# distance is computed twice in this case
|
||||
if t_col == 0 and (t_pln is not 0 or t_row is not 0):
|
||||
alpha = 0.5
|
||||
else:
|
||||
alpha = 1.
|
||||
# Compute integral image of the squared difference between
|
||||
# padded and the same image shifted by (t_pln, t_row, t_col)
|
||||
integral = np.zeros_like(padded)
|
||||
_integral_image_3d(padded, integral, t_pln, t_row, t_col,
|
||||
n_pln, n_row, n_col)
|
||||
|
||||
# Inner loops on pixel coordinates
|
||||
# Iterate over planes, taking offset and shift into account
|
||||
for pln in range(pln_dist_min, pln_dist_max):
|
||||
# Iterate over rows, taking offset and shift into account
|
||||
for row in range(row_dist_min, row_dist_max):
|
||||
# Iterate over columns
|
||||
for col in range(col_dist_min, col_dist_max):
|
||||
# Compute squared distance between shifted patches
|
||||
distance = _integral_to_distance_3d(integral,
|
||||
pln, row, col, offset, s_cube_h_square)
|
||||
# exp of large negative numbers is close to zero
|
||||
if distance > DISTANCE_CUTOFF:
|
||||
continue
|
||||
weight = alpha * exp(-distance)
|
||||
# Accumulate weights for the different shifts
|
||||
weights[pln, row, col] += weight
|
||||
weights[pln + t_pln, row + t_row,
|
||||
col + t_col] += weight
|
||||
result[pln, row, col] += weight * \
|
||||
padded[pln + t_pln, row + t_row,
|
||||
col + t_col]
|
||||
result[pln + t_pln, row + t_row,
|
||||
col + t_col] += weight * \
|
||||
padded[pln, row, col]
|
||||
|
||||
# Normalize pixel values using sum of weights of contributing patches
|
||||
for pln in range(offset, n_pln - offset):
|
||||
for row in range(offset, n_row - offset):
|
||||
for col in range(offset, n_col - offset):
|
||||
# No risk of division by zero, since the contribution
|
||||
# of a null shift is strictly positive
|
||||
result[pln, row, col] /= weights[pln, row, col]
|
||||
|
||||
# Return cropped result, undoing padding
|
||||
return result[pad_size:-pad_size, pad_size:-pad_size, pad_size:-pad_size]
|
||||
@@ -0,0 +1,120 @@
|
||||
import numpy as np
|
||||
from skimage.restoration._nl_means_denoising import _nl_means_denoising_2d, \
|
||||
_nl_means_denoising_3d, \
|
||||
_fast_nl_means_denoising_2d, _fast_nl_means_denoising_3d
|
||||
|
||||
def nl_means_denoising(image, patch_size=7, patch_distance=11, h=0.1,
|
||||
multichannel=True, fast_mode=True):
|
||||
"""
|
||||
Perform non-local means denoising on 2-D or 3-D grayscale images, and
|
||||
2-D RGB images.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
image : 2D or 3D ndarray
|
||||
Input image to be denoised, which can be 2D or 3D, and grayscale
|
||||
or RGB (for 2D images only, see ``multichannel`` parameter).
|
||||
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. 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.
|
||||
multichannel : bool, optional
|
||||
Whether the last axis of the image is to be interpreted as multiple
|
||||
channels or another spatial dimension. Set to ``False`` for 3-D images.
|
||||
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
|
||||
Denoised image, of same shape as `image`.
|
||||
|
||||
See Also
|
||||
--------
|
||||
fast_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.
|
||||
|
||||
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.
|
||||
|
||||
However, the default behavior corresponds to ``fast_mode=True``, for which
|
||||
another version of non-local means [2]_ is used, corresponding to a
|
||||
complexity of
|
||||
|
||||
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. 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] 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
|
||||
--------
|
||||
>>> 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:
|
||||
image = image[..., np.newaxis]
|
||||
multichannel = True
|
||||
if image.ndim != 3:
|
||||
raise NotImplementedError("Non-local means denoising is only \
|
||||
implemented for 2D grayscale and RGB images or 3-D grayscale images.")
|
||||
if multichannel: # 2-D images
|
||||
if fast_mode:
|
||||
return np.squeeze(np.array(_fast_nl_means_denoising_2d(image,
|
||||
patch_size, patch_distance, h)))
|
||||
else:
|
||||
return np.squeeze(np.array(_nl_means_denoising_2d(image,
|
||||
patch_size, patch_distance, h)))
|
||||
else: # 3-D 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))
|
||||
|
||||
@@ -17,6 +17,7 @@ def configuration(parent_package='', top_path=None):
|
||||
cython(['_unwrap_2d.pyx'], working_path=base_path)
|
||||
cython(['_unwrap_3d.pyx'], working_path=base_path)
|
||||
cython(['_denoise_cy.pyx'], working_path=base_path)
|
||||
cython(['_nl_means_denoising.pyx'], working_path=base_path)
|
||||
|
||||
config.add_extension('_unwrap_1d', sources=['_unwrap_1d.c'],
|
||||
include_dirs=[get_numpy_include_dirs()])
|
||||
@@ -28,6 +29,9 @@ def configuration(parent_package='', top_path=None):
|
||||
include_dirs=[get_numpy_include_dirs()])
|
||||
config.add_extension('_denoise_cy', sources=['_denoise_cy.c'],
|
||||
include_dirs=[get_numpy_include_dirs(), '../_shared'])
|
||||
config.add_extension('_nl_means_denoising',
|
||||
sources=['_nl_means_denoising.c'],
|
||||
include_dirs=[get_numpy_include_dirs()])
|
||||
|
||||
return config
|
||||
|
||||
|
||||
@@ -46,7 +46,8 @@ def test_denoise_tv_chambolle_float_result_range():
|
||||
img = astro_gray
|
||||
int_astro = np.multiply(img, 255).astype(np.uint8)
|
||||
assert np.max(int_astro) > 1
|
||||
denoised_int_astro = restoration.denoise_tv_chambolle(int_astro, weight=60.0)
|
||||
denoised_int_astro = restoration.denoise_tv_chambolle(int_astro,
|
||||
weight=60.0)
|
||||
# test if the value range of output float data is within [0.0:1.0]
|
||||
assert denoised_int_astro.dtype == np.float
|
||||
assert np.max(denoised_int_astro) <= 1.0
|
||||
@@ -143,5 +144,76 @@ def test_denoise_bilateral_3d():
|
||||
assert out1[30:45, 5:15].std() > out2[30:45, 5:15].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 = restoration.nl_means_denoising(img, 7, 5, 0.2, fast_mode=True)
|
||||
# make sure noise is reduced
|
||||
assert img.std() > denoised.std()
|
||||
denoised = restoration.nl_means_denoising(img, 7, 5, 0.2, fast_mode=False)
|
||||
# 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])
|
||||
# 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.3, fast_mode=True)
|
||||
# make sure noise is reduced
|
||||
assert img.std() > denoised.std()
|
||||
denoised = restoration.nl_means_denoising(img, 7, 9, 0.3, fast_mode=False)
|
||||
# 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 = restoration.nl_means_denoising(img, 5, 4, 0.2, fast_mode=True,
|
||||
multichannel=False)
|
||||
# make sure noise is reduced
|
||||
assert img.std() > denoised.std()
|
||||
denoised = restoration.nl_means_denoising(img, 5, 4, 0.2, fast_mode=False,
|
||||
multichannel=False)
|
||||
# make sure noise is reduced
|
||||
assert img.std() > denoised.std()
|
||||
|
||||
|
||||
def test_nl_means_denoising_multichannel():
|
||||
img = np.zeros((21, 20, 10))
|
||||
img[10, 9:11, 2:-2] = 1.
|
||||
img += 0.3*np.random.randn(*img.shape)
|
||||
denoised_wrong_multichannel = restoration.nl_means_denoising(img,
|
||||
5, 4, 0.1, fast_mode=True, multichannel=True)
|
||||
denoised_ok_multichannel = restoration.nl_means_denoising(img,
|
||||
5, 4, 0.1, fast_mode=True, multichannel=False)
|
||||
snr_wrong = 10 * np.log10(1. /
|
||||
((denoised_wrong_multichannel - img)**2).mean())
|
||||
snr_ok = 10 * np.log10(1. /
|
||||
((denoised_ok_multichannel - img)**2).mean())
|
||||
assert snr_ok > snr_wrong
|
||||
|
||||
|
||||
def test_nl_means_denoising_wrong_dimension():
|
||||
img = np.zeros((5, 5, 5, 5))
|
||||
assert_raises(NotImplementedError, restoration.nl_means_denoising, img)
|
||||
|
||||
|
||||
def test_no_denoising_for_small_h():
|
||||
img = np.zeros((40, 40))
|
||||
img[10:-10, 10:-10] = 1.
|
||||
img += 0.3*np.random.randn(*img.shape)
|
||||
# very small h should result in no averaging with other patches
|
||||
denoised = restoration.nl_means_denoising(img, 7, 5, 0.01, fast_mode=True)
|
||||
assert np.allclose(denoised, img)
|
||||
denoised = restoration.nl_means_denoising(img, 7, 5, 0.01, fast_mode=False)
|
||||
assert np.allclose(denoised, img)
|
||||
|
||||
|
||||
if __name__ == "__main__":
|
||||
run_module_suite()
|
||||
|
||||
Reference in New Issue
Block a user