Nonlocal Means Denoising¶
Find similar image patches and average their center pixels using similarity weights. Explore the effects of patch size and search windows, then examine how a lower-dimensional patch representation changes the estimate.
Run this tour¶
Run the cells in order with a Python 3 kernel. The first cell locates the companion data and toolbox and installs missing dependencies when needed. All worked examples include their implementation directly in this notebook. Random seeds make comparisons reproducible; you can change them to explore other samples.
# Locate the companion toolbox locally, or fetch it for a standalone/Colab copy.
from pathlib import Path
import importlib.util
import os
import subprocess
import sys
working = Path.cwd()
candidates = [working, working / "python", working.parent / "python"]
python_dir = next((p for p in candidates if (p / "nt_toolbox").is_dir()), None)
if python_dir is None:
checkout = working / "numerical-tours-support"
if not checkout.exists():
subprocess.run(
[
"git",
"clone",
"--depth",
"1",
"--branch",
"master",
"https://github.com/gpeyre/numerical-tours.git",
str(checkout),
],
check=True,
)
python_dir = checkout / "python"
os.chdir(python_dir)
if str(python_dir) not in sys.path:
sys.path.insert(0, str(python_dir))
requirements = python_dir / "requirements.txt"
if any(
importlib.util.find_spec(name) is None
for name in [
"numpy",
"scipy",
"matplotlib",
"skimage",
"sklearn",
"pywt",
"ipywidgets",
"cvxpy",
"skfmm",
"autograd",
"progressbar",
"celer",
]
):
subprocess.run(
[sys.executable, "-m", "pip", "install", "-r", str(requirements)], check=True
)
import numpy as np
import matplotlib.pyplot as plt
np.random.seed(0)
plt.rcParams.update(
{
"figure.figsize": (8, 4),
"figure.dpi": 100,
"axes.spines.top": False,
"axes.spines.right": False,
"font.size": 11,
"image.cmap": "gray",
}
)
%matplotlib inline
$\newcommand{\dotp}[2]{\langle #1, #2 \rangle}$ $\newcommand{\enscond}[2]{\lbrace #1, #2 \rbrace}$ $\newcommand{\pd}[2]{ \frac{ \partial #1}{\partial #2} }$ $\newcommand{\umin}[1]{\underset{#1}{\min}\;}$ $\newcommand{\umax}[1]{\underset{#1}{\max}\;}$ $\newcommand{\uargmin}[1]{\underset{#1}{argmin}\;}$ $\newcommand{\norm}[1]{\|#1\|}$ $\newcommand{\abs}[1]{\left|#1\right|}$ $\newcommand{\choice}[1]{ \left\{ \begin{array}{l} #1 \end{array} \right. }$ $\newcommand{\pa}[1]{\left(#1\right)}$ $\newcommand{\diag}[1]{{diag}\left( #1 \right)}$ $\newcommand{\qandq}{\quad\text{and}\quad}$ $\newcommand{\qwhereq}{\quad\text{where}\quad}$ $\newcommand{\qifq}{ \quad \text{if} \quad }$ $\newcommand{\qarrq}{ \quad \Longrightarrow \quad }$ $\newcommand{\ZZ}{\mathbb{Z}}$ $\newcommand{\CC}{\mathbb{C}}$ $\newcommand{\RR}{\mathbb{R}}$ $\newcommand{\EE}{\mathbb{E}}$ $\newcommand{\Zz}{\mathcal{Z}}$ $\newcommand{\Ww}{\mathcal{W}}$ $\newcommand{\Vv}{\mathcal{V}}$ $\newcommand{\Nn}{\mathcal{N}}$ $\newcommand{\NN}{\mathcal{N}}$ $\newcommand{\Hh}{\mathcal{H}}$ $\newcommand{\Bb}{\mathcal{B}}$ $\newcommand{\Ee}{\mathcal{E}}$ $\newcommand{\Cc}{\mathcal{C}}$ $\newcommand{\Gg}{\mathcal{G}}$ $\newcommand{\Ss}{\mathcal{S}}$ $\newcommand{\Pp}{\mathcal{P}}$ $\newcommand{\Ff}{\mathcal{F}}$ $\newcommand{\Xx}{\mathcal{X}}$ $\newcommand{\Mm}{\mathcal{M}}$ $\newcommand{\Ii}{\mathcal{I}}$ $\newcommand{\Dd}{\mathcal{D}}$ $\newcommand{\Ll}{\mathcal{L}}$ $\newcommand{\Tt}{\mathcal{T}}$ $\newcommand{\si}{\sigma}$ $\newcommand{\al}{\alpha}$ $\newcommand{\la}{\lambda}$ $\newcommand{\ga}{\gamma}$ $\newcommand{\Ga}{\Gamma}$ $\newcommand{\La}{\Lambda}$ $\newcommand{\Si}{\Sigma}$ $\newcommand{\be}{\beta}$ $\newcommand{\de}{\delta}$ $\newcommand{\De}{\Delta}$ $\newcommand{\phi}{\varphi}$ $\newcommand{\th}{\theta}$ $\newcommand{\om}{\omega}$ $\newcommand{\Om}{\Omega}$
This numerical tour study image denoising using non-local means. This algorithm has been introduced for denoising purposes in BuaCoMoA05
import numpy as np
import scipy as scp
import pylab as pyl
import matplotlib.pyplot as plt
from nt_toolbox.general import clamp, np, plt, pylab, rescale
from nt_toolbox.signal import imageplot, load_image, np, plt, pylab, snr
%matplotlib inline
Patches in Images¶
This numerical tour is dedicated to the study of the structure of patches in images.
Size $N = n \times n$ of the image.
n = 128
We load a noisy image $f_0\in \RR^N$.
c = [100, 200]
f0 = load_image("nt_toolbox/data/lena.bmp")
f0 = rescale(f0[c[0] - n // 2 : c[0] + n // 2, c[1] - n // 2 : c[1] + n // 2])
Display $f_0$.
plt.figure(figsize=(5, 5))
imageplot(f0)
Noise level $\si$.
sigma = 0.04
Generate a noisy image $f=f_0+\epsilon$ where $\epsilon \times \Nn(0,\si^2\text{Id}_N)$.
from numpy import random
f = f0 + sigma * random.standard_normal((n, n))
Display $f$.
plt.figure(figsize=(5, 5))
imageplot(clamp(f))
We denote $w$ to be the half width of the patches, and $w_1=2w+1$ the full width.
w = 3
w1 = 2 * w + 1
We set up large $(n,n,w_1,w_1)$ matrices to index the the X and Y position of the pixel to extract.
Location of pixels to extract.
[X, Y, dX, dY] = np.meshgrid(
np.arange(1, n + 1), np.arange(1, n + 1), np.arange(-w, w + 1), np.arange(-w, w + 1)
)
X = X + dX
Y = Y + dY
We handle boundary condition by reflexion
X[X < 1] = 2 - X[X < 1]
Y[Y < 1] = 2 - Y[Y < 1]
X[X > n] = 2 * n - X[X > n]
Y[Y > n] = 2 * n - Y[Y > n]
Patch extractor operator.
I = (X - 1) + (Y - 1) * n
for i in range(n // w):
for j in range(n // w):
I[i, j] = np.transpose(I[i, j])
patch = lambda f: np.ravel(f)[I]
Define the patch matrix $P$ of size $(n,n,w_1,w_1)$. Each $P(i,j,:,:)$ represent an $(w_1,w_1)$ patch extracted around pixel $(i,j)$ in the image.
P = patch(f)
Display some example of patches.
from numpy import random
plt.figure(figsize=(5, 5))
for i in range(16):
x = random.randint(n)
y = random.randint(n)
imageplot(P[x, y], "", [4, 4, i + 1])
Dimensionality Reduction with PCA¶
Since NL-means type algorithms require the computation of many distances between patches, it is advantagous to reduce the dimensionality of the patch while keeping as much as possible of information.
Target dimensionality $d$.
d = 25
A linear dimensionality reduction is obtained by Principal Component Analysis (PCA) that projects the data on a small number of leading direction of the covariance matrix of the patches.
Turn the patch matrix into an $(w_1*w_1,n*n)$ array, so that each $P(:,i)$ is a $w_1*w_1$ vector representing a patch.
resh = lambda P: np.transpose((np.reshape(P, (n * n, w1 * w1), order="F")))
Operator to remove the mean of the patches to each patch.
remove_mean = lambda Q: Q - np.tile(np.mean(Q, 0), (w1 * w1, 1))
Compute the mean and the covariance of the points cloud representing the patches.
P1 = remove_mean(resh(P))
C = np.dot(P1, np.transpose(P1))
Extract the eigenvectors, sorted by decreasing amplitude.
from numpy import linalg
D, V = linalg.eigh(C)
I = np.argsort(D)[::-1]
D = D[I]
V = V[:, I]
Display the decaying amplitude of the eigenvalues.
plt.plot(D, linewidth=2)
plt.ylim(0, max(D))
plt.show()
Display the leading eigenvectors - they look like Fourier modes.
plt.figure(figsize=(5, 5))
for i in range(16):
imageplot(abs(np.reshape(V[:, i], (w1, w1))), "", [4, 4, i + 1])
Patch dimensionality reduction operator.
iresh = lambda Q: np.reshape(np.transpose(Q), (n, n, d), order="F")
descriptor = lambda f: iresh(np.dot(np.transpose(V[:, :d]), remove_mean(resh(P))))
Each $H(i,j,:)$ is a $d$-dimensional descriptor of a patch.
H = descriptor(f)
Non-local Filter¶
NL-means applies, an adaptive averaging kernel is computed from patch distances to each pixel location.
We denote $H_{i} \in \RR^d$ the descriptor at pixel $i$. We define the distance matrix $$ D_{i,j} = \frac{1}{w_1^2}\norm{H_i-H_j}^2. $$
Operator to compute the distances $(D_{i,j})_j$ between the patch around $i=(i_1,i_2)$ and all the other ones.
distance = lambda i: (
np.sum((H - np.tile(H[i[0], i[1], :], (n, n, 1))) ** 2, 2) / (w1 * w1)
)
The non-local mean filter computes a denoised image $\tilde f$ as :
$$ \tilde f_i = \sum_j K_{i,j} f_j $$ where the weights $K$ are computed as : $$ K_{i,j} = \frac{ \tilde K_{i,j} }{ \sum_{j'} \tilde K_{i,j'} } \qandq \tilde K_{i,j} = e^{-\frac{D_{i,j}}{2\tau^2}} . $$
The width $\tau$ of the Gaussian is very important and should be adapted to match the noise level.
Compute and normalize the weight.
normalize = lambda K: K / np.sum(K)
kernel = lambda i, tau: normalize(np.exp(-distance(i) / (2 * tau**2)))
Compute a typical example of kernel for some pixel position $(x,y)$.
tau = 0.05
i = [83, 72]
D = distance(i)
K = kernel(i, tau)
Display the squared distance and the kernel.
plt.figure(figsize=(10, 10))
imageplot(D, "D", [1, 2, 1])
imageplot(K, "K", [1, 2, 2])
Localizing the Non-local Means¶
We set a "locality constant" $q$ that set the maximum distance between patches to compare. This allows to speed up computation, and makes NL-means type methods semi-global (to avoid searching in all the image).
q = 14
Using this locality constant, we compute the distance between patches only within a window. Once again, one should be careful about boundary conditions.
selection = lambda i: np.array(
(
clamp(np.arange(i[0] - q, i[0] + q + 1), 0, n - 1),
clamp(np.arange(i[1] - q, i[1] + q + 1), 0, n - 1),
)
)
Compute distance and kernel only within the window.
def distance_0(i, sel):
H1 = H[sel[0], :, :]
H2 = H1[:, sel[1], :]
return np.sum(
(H2 - np.tile(H[i[0], i[1], :], (len(sel[0]), len(sel[1]), 1))) ** 2, 2
) / (w1 * w1)
distance = lambda i: distance_0(i, selection(i))
kernel = lambda i, tau: normalize(np.exp(-distance(i) / (2 * tau**2)))
Compute a typical example of kernel for some pixel position $(x,y)$.
D = distance(i)
K = kernel(i, tau)
Display the squared distance and the kernel.
plt.figure(figsize=(10, 10))
imageplot(D, "D", [1, 2, 1])
imageplot(K, "K", [1, 2, 2])
The NL-filtered value at pixel $(x,y)$ is obtained by averaging the values of $f$ with the weight $K$.
def NLval_0(K, sel):
f_temp = f[sel[0], :]
return np.sum(K * f_temp[:, sel[1]])
NLval = lambda i, tau: NLval_0(kernel(i, tau), selection(i))
We apply the filter to each pixel location to perform the NL-means algorithm.
[Y, X] = np.meshgrid(np.arange(0, n), np.arange(0, n))
def arrayfun(f, X, Y):
n = len(X)
p = len(Y)
R = np.zeros([n, p])
for k in range(n):
for l in range(p):
R[k, l] = f(k, l)
return R
NLmeans = lambda tau: arrayfun(lambda i1, i2: NLval([i1, i2], tau), X, Y)
Display the result for some value of $\tau$.
tau = 0.03
plt.figure(figsize=(5, 5))
imageplot(NLmeans(tau))
Worked example 1
Compute the denoising result for several values of $\tau$ in order to determine the optimal denoising that minimizes $\norm{\tilde f - f_0}$.
plt.figure(figsize=(10, 10))
ntests = 4
tau_list = np.linspace(0.03, 0.25, ntests)
e0 = -float("inf")
for k in range(ntests):
tau = tau_list[k]
f1 = NLmeans(tau)
e = snr(f0, f1)
if e > e0:
fNL = f1
e0 = e
imageplot(clamp(f1), "tau = %.3f, SNR = %.1f dB" % (tau, e), [2, 2, k + 1])
Display the best result.
plt.figure(figsize=(5, 5))
imageplot(clamp(fNL))
Worked example 2
Explore the influence of the $q$ and $w$ parameters.
plt.figure(figsize=(10, 10))
tau = 0.03
q_list = [10, 20]
w_list = [3, 6]
ind_plot = 0
for i_q in range(len(q_list)):
for i_w in range(len(w_list)):
w = w_list[i_w]
q = q_list[i_q]
ind_plot += 1
# patch
w1 = 2 * w + 1
[X, Y, dX, dY] = np.meshgrid(
np.arange(1, n + 1),
np.arange(1, n + 1),
np.arange(-w, w + 1),
np.arange(-w, w + 1),
)
X = X + dX
Y = Y + dY
X[X < 1] = 2 - X[X < 1]
Y[Y < 1] = 2 - Y[Y < 1]
X[X > n] = 2 * n - X[X > n]
Y[Y > n] = 2 * n - Y[Y > n]
I = (X - 1) + (Y - 1) * n
for i in range(n // w):
for j in range(n // w):
I[i, j] = np.transpose(I[i, j])
patch = lambda f: np.ravel(f)[I]
P = patch(f)
# PCA
resh = lambda P0: np.transpose((np.reshape(P0, (n * n, w1 * w1), order="F")))
remove_mean = lambda Q: Q - np.tile(np.mean(Q, 0), (w1 * w1, 1))
P1 = remove_mean(resh(P))
C = np.dot(P1, np.transpose(P1))
D, V = linalg.eigh(C)
I = np.argsort(D)[::-1]
D = D[I]
V = V[:, I]
iresh = lambda Q: np.reshape(np.transpose(Q), (n, n, d), order="F")
descriptor = lambda f: iresh(
np.dot(np.transpose(V[:, :d]), remove_mean(resh(P)))
)
H = descriptor(f)
# NL_means
i = [83, 72]
def distance_0(i, sel):
H1 = H[sel[0], :, :]
H2 = H1[:, sel[1], :]
return np.sum(
(H2 - np.tile(H[i[0], i[1], :], (len(sel[0]), len(sel[1]), 1))) ** 2, 2
) / (w1 * w1)
distance = lambda i: distance_0(i, selection(i))
kernel = lambda i, tau: normalize(np.exp(-distance(i) / (2 * tau**2)))
selection = lambda i: np.array(
(
clamp(np.arange(i[0] - q, i[0] + q + 1), 0, n - 1),
clamp(np.arange(i[1] - q, i[1] + q + 1), 0, n - 1),
)
)
def NLval_0(K, sel):
f_temp = f[sel[0], :]
return np.sum(K * f_temp[:, sel[1]])
NLval = lambda i, tau: NLval_0(kernel(i, tau), selection(i))
NLmeans = lambda tau: arrayfun(lambda i1, i2: NLval([i1, i2], tau), X, Y)
f1 = NLmeans(tau)
imageplot(
clamp(f1),
"q = %i, w = %i, SNR = %.1f dB" % (q, w, snr(f0, f1)),
[2, 2, ind_plot],
)
References and further reading¶
- [BuaCoMoA05] Buades, B. Coll, J.f Morel, [A review of image denoising algorithms, with a new one][1], SIAM Multiscale Modeling and Simulation, Vol 4 (2), pp: 490-530, 2005.
Antoni Buades, Bartomeu Coll, and Jean-Michel Morel. A Non-Local Algorithm for Image Denoising. 2005, CVPR, vol. 2, 60–65. Patch similarity as the basis for nonlocal averaging.
Stéphane Mallat. A Wavelet Tour of Signal Processing: The Sparse Way. 2009, 3rd ed., Academic Press. Multiresolution analysis, sparse approximation, and wavelet algorithms.
Stéfan van der Walt et al.. scikit-image: image processing in Python. 2014, PeerJ 2:e453. Reproducible image processing and image measurements in Python.
Charles M. Stein. Estimation of the Mean of a Multivariate Normal Distribution. 1981, Annals of Statistics 9(6), 1135–1151. Unbiased risk estimation for Gaussian observations.