Scikit Image − Non-local means denoising for preserving textures

The non-local means algorithm is a powerful denoising technique that works well for images with specific textures. Its fundamental concept is to average the value of a pixel with values from other pixels in a limited neighborhood, as long as the patches centered on those other pixels are similar to the patch centered on the pixel of interest.

The scikit-image library offers a dedicated function for applying the non-local means denoising algorithm to images. This function provides users with the flexibility to adjust various parameters such as patch_size, patch_distance, and h. Additionally, users can choose between two modes: "fast" and "classic," depending on their computational requirements.

Using the skimage.restoration.denoise_nl_means() function

The restoration.denoise_nl_means() function performs non-local means denoising on 2D-4D grayscale or RGB images.

Syntax

Following is the syntax of this function −

skimage.restoration.denoise_nl_means(image, patch_size=7, patch_distance=11, h=0.1, fast_mode=True, sigma=0.0, *, preserve_range=False, channel_axis=None)

Parameters

Here's an explanation of its parameters −

• image (2D or 3D ndarray): This is the input image to be denoised. It can be 2D or 3D and can be grayscale or RGB. There can be any number of channels (e.g., for RGB images), and you can specify the channel axis using the channel_axis parameter.
• patch_size (int, optional): it specifies the size of patches used for denoising.
• patch_distance (int, optional): The maximal distance in pixels where to search for patches used for denoising.
• h (float, optional): Cut-off distance (in gray levels). A higher value of h results in a smoother image but may blur features. It controls how permissive the algorithm is in accepting patches. For a Gaussian noise choosing the h value will be based on the standard deviation of the noise.
• fast_mode (bool, optional): If True, a faster version of the non-local means algorithm is used. If False, the original version of non-local means is used.
• sigma (float, optional): The standard deviation of the Gaussian noise. If provided, the algorithm computes patch weights that take the expected noise variance into account.
• preserve_range (bool, optional): Whether to keep the original range of values in the image. If True, the image is not rescaled. If False, the input image is converted to a floating-point format ( according to the conventions of img_as_float).
• channel_axis (int or None, optional): If not None, it indicates which axis of the array corresponds to channels (e.g., for RGB images). If None, the image is assumed to be a grayscale (single channel) image. This parameter was added in version 0.19.

Return value

The function returns a denoised image (ndarray) with the same shape as the input image.

Example

The following example applies the non-local means denoising algorithm to images −

import numpy as np
import matplotlib.pyplot as plt
from skimage.restoration import denoise_nl_means

# Create a noisy image
image = np.zeros((40, 40))
image[10:-10, 10:-10] = 1.
rng = np.random.default_rng()
image += 0.3 * rng.standard_normal(image.shape)

# Denoise the image using non-local means
denoised_image = denoise_nl_means(image, patch_size=7, patch_distance=5, h=0.1)

# Create subplots for original and denoised images
fig, axes = plt.subplots(1, 2, figsize=(10, 5))

# Plot the original noisy image
axes[0].imshow(image, cmap='gray')
axes[0].set_title('Noisy Image')
axes[0].axis('off')

# Plot the denoised image
axes[1].imshow(denoised_image, cmap='gray')
axes[1].set_title('Denoised Image')
axes[1].axis('off')

plt.tight_layout()
plt.show()

Output

Example

This example demonstrates the application of non-local means (NLM) denoising to an image to remove noise while preserving image details. It also shows how parameter settings and algorithms (slow vs. fast) can impact the denoising results −

import numpy as np
import matplotlib.pyplot as plt
from skimage import io, img_as_float
from skimage.restoration import denoise_nl_means, estimate_sigma
from skimage.metrics import peak_signal_noise_ratio
from skimage.util import random_noise

# Load the input image and select a region of interest
image = img_as_float(io.imread('Images/Tajmahal.jpg'))

# Add Gaussian noise to the image
sigma = 0.08
noisy = random_noise(image, var=sigma**2)

# Estimate the noise standard deviation from the noisy image
sigma_est = np.mean(estimate_sigma(noisy, channel_axis=-1))
print(f'Estimated noise standard deviation = {sigma_est:.4f}')

# Define parameters for non-local means denoising
patch_kw = dict(patch_size=5, patch_distance=6, channel_axis=-1)

# Apply non-local means denoising with slow algorithm
denoise = denoise_nl_means(noisy, h=1.15 * sigma_est, fast_mode=False, **patch_kw)

# Apply non-local means denoising with slow algorithm and sigma provided
denoise2 = denoise_nl_means(noisy, h=0.8 * sigma_est, sigma=sigma_est, fast_mode=False, **patch_kw)

# Apply non-local means denoising with fast algorithm
denoise_fast = denoise_nl_means(noisy, h=0.8 * sigma_est, fast_mode=True, **patch_kw)

# Apply non-local means denoising with fast algorithm and sigma provided
denoise2_fast = denoise_nl_means(noisy, h=0.6 * sigma_est, sigma=sigma_est, fast_mode=True, **patch_kw)

# Create subplots for visualization
fig, axes = plt.subplots(nrows=2, ncols=3, figsize=(8, 6), sharex=True, sharey=True)
ax0, ax1, ax2, ax3, ax4, ax5 = axes.ravel()

# Plot the noisy image
ax0.imshow(noisy)
ax0.set_title('Noisy')

# Plot denoised images with different configurations
ax1.imshow(denoise)
ax1.set_title('Non-local means\n(Slow)')

ax2.imshow(denoise2)
ax2.set_title('Non-local means\n(Slow, using $\\sigma_{est}$)')

# Plot the original noise-free image
ax3.imshow(image)
ax3.set_title('Original\n(Noise-Free)')

ax4.imshow(denoise_fast)
ax4.set_title('Non-local means\n(Fast)')

ax5.imshow(denoise2_fast)
ax5.set_title('Non-local means\n(Fast, using $\\sigma_{est}$)')

for ax in axes.flat:
    ax.axis('off')
fig.tight_layout()

# Calculate and print PSNR metric for each case
psnr_noisy = peak_signal_noise_ratio(image, noisy)
psnr = peak_signal_noise_ratio(image, denoise)
psnr2 = peak_signal_noise_ratio(image, denoise2)
psnr_fast = peak_signal_noise_ratio(image, denoise_fast)
psnr2_fast = peak_signal_noise_ratio(image, denoise2_fast)

print(f'PSNR (Noisy) = {psnr_noisy:.2f}')
print(f'PSNR (Slow) = {psnr:.2f}')
print(f'PSNR (slow, using sigma) = {psnr2:0.2f}')
print(f'PSNR (Fast) = {psnr_fast:.2f}')
print(f'PSNR (fast, using sigma) = {psnr2_fast:0.2f}')

plt.show()

Output

Estimated noise standard deviation = 0.0770
PSNR (Noisy) = 22.35
PSNR (Slow) = 29.39
PSNR (slow, using sigma) = 29.45
PSNR (Fast) = 28.64
PSNR (fast, using sigma) = 29.09