Wavelet Block Thresholding¶
Process neighboring wavelet coefficients together instead of thresholding each coefficient independently. The experiments show how shared local energy changes the balance between texture preservation and noise removal.
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 presents block thresholding methods, that makes use of the structure of wavelet coefficients of natural images to perform denoising. Theoretical properties of block thresholding were investigated in CaiSilv Cai99 HallKerkPic99
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, plot_wavelet, plt, pylab, snr
import warnings
%matplotlib inline
Generating a Noisy Image¶
Here we use an additive Gaussian noise.
Size of the image of $N=n \times n$ pixels.
n = 256
First we load an image $f_0 \in \RR^N$.
f0 = rescale(load_image("nt_toolbox/data/boat.bmp", n))
Display it.
plt.figure(figsize=(5, 5))
imageplot(f0)
Noise level.
sigma = 0.08
Generate a noisy image $f=f_0+\epsilon$ where $\epsilon \sim \Nn(0,\si^2\text{Id}_N)$.
from numpy import random
f = f0 + sigma * random.standard_normal((n, n))
Display it.
plt.figure(figsize=(5, 5))
imageplot(clamp(f))
Orthogonal Wavelet Thresholding¶
We first consider the traditional wavelet thresholding method.
Parameters for the orthogonal wavelet transform.
Jmin = 4
Shortcuts for the foward and backward wavelet transforms.
from nt_toolbox.perform_wavelet_transf import np, perform_wavelet_transf
wav = lambda f: perform_wavelet_transf(f, Jmin, +1)
iwav = lambda fw: perform_wavelet_transf(fw, Jmin, -1)
Display the original set of noisy coefficients.
plt.figure(figsize=(10, 10))
plot_wavelet(wav(f), Jmin)
plt.show()
Denoting $\Ww$ and $\Ww^*$ the forward and backward wavelet transform, wavelet thresholding $\tilde f$ is defined as
$$ \tilde f = \Ww^* \circ \theta_T \circ \Ww(f) $$
where $T>0$ is the threshold, that should be adapted to the noise level.
The thresholding operator is applied component-wise
$$ \th_T(x)_i = \psi_T(x_i) x_i $$
where $\psi_T$ is an atenuation fonction. In this tour, we use the James Stein (JS) attenuation:
$$ \psi_T(s) = \max\pa{ 0, 1-\frac{T^2}{s^2} } $$
psi = lambda s, T: np.maximum(
1 - T**2 / np.maximum(abs(s) ** 2, 1e-9 * np.ones(np.shape(s))),
np.zeros(np.shape(s)),
)
Display the thresholding function $\th_T$.
s = np.linspace(-3, 3, 1024)
plt.plot(s, s * psi(s, 1))
plt.plot(s, s, "r--")
plt.show()
Thresholding operator.
theta = lambda x, T: psi(x, T) * x
ThreshWav = lambda f, T: iwav(theta(wav(f), T))
Test the thresholding.
T = 1.5 * sigma
plt.figure(figsize=(5, 5))
imageplot(clamp(ThreshWav(f, T)))
Worked example 1
Display the evolution of the denoising SNR when $T$ varies. Store in $f_{Thresh}$ the optimal denoising result.
Tlist = np.linspace(0.5, 2.5, 30) * sigma
snr_stein = [snr(f0, ThreshWav(f, T)) for T in Tlist]
plt.plot([T / sigma for T in Tlist], snr_stein, linewidth=2)
plt.xlabel(r"$T/\sigma$")
plt.ylabel("SNR")
plt.show()
Tmax = np.max(snr_stein)
fThresh = ThreshWav(f, T)
Display the optimal thresolding.
plt.figure(figsize=(5, 5))
imageplot(clamp(fThresh), "SNR = %.1f dB" % snr(f0, fThresh))
Block Thresholding Operator¶
A block thresholding operator of coefficients $x=(x_i)_{i=1}^P \in \RR^P$ is defined using a partition $B$ into a set of blocks $b$
$$ \{1,\ldots,P\} = \bigcup_{b \in B} b. $$ Its definition reads
$$ \forall i \in b, \quad \theta_T(x)_i = \psi_T\left( \norm{x_b}_2 \right) x_i $$ where $ x_b = (x_j)_{j \in B} \in \RR^{\abs{b}} $. One thus thresholds the $\ell^2$ norm (the energy) of each block rather than each coefficient independently.
For image-based thresholding, we use a partition in square blocks of equal size $w \times w$.
The block size $w$.
w = 4
n = 256
Compute indexing of the blocks.
[X, Y, dX, dY] = np.meshgrid(
np.arange(1, n - w + 2, w),
np.arange(1, n - w + 2, w),
np.arange(0, w),
np.arange(0, w),
)
I = (X + dX - 1) + (Y + dY - 1) * n
for k in range(n // w):
for l in range(n // w):
I[k][l] = np.transpose(I[k][l])
Block extraction operator. It returns the set $ \{x_b\}_{b \in B} $ of block-partitioned coefficients.
block = lambda x: np.ravel(x)[I]
Block reconstruction operator.
def assign(M, I, H):
M_temp = M
np.ravel(M_temp)[I] = H
return np.reshape(M_temp, (n, n))
iblock = lambda H: assign(np.zeros([n, n]), I, H)
Check that block extraction / reconstruction gives perfect reconstruction.
from numpy import linalg
print("Should be 0:", linalg.norm(f - iblock(block(f))))
Compute the average energy of each block, and duplicate.
def energy(H):
H_tmp = np.copy(H)
for i in range(n // w):
for j in range(n // w):
H_tmp[i][j] = np.sqrt(np.mean(H_tmp[i][j] ** 2)) # *np.ones([1,1])
return H_tmp
Block thresholding operator.
Thresh = lambda H, T: psi(energy(H), T) * H
ThreshBlock = lambda x, T: iblock(Thresh(block(x), T))
Worked example 2
Test the effect of block thresholding on the image $f_0$ itself, for increasing value of $T$. Of course directly thresholding the image has no interest, this is just to vizualize the effect.
plt.figure(figsize=(10, 10))
tlist = np.linspace(0.3, 0.9, 4)
for i in range(len(tlist)):
T = tlist[i]
imageplot(clamp(ThreshBlock(f, T)), "T = %.1f" % T, [2, 2, i + 1])
Orthogonal Wavelet Block Thresholding¶
Wavelet coefficients of natural images are not independant one from each other. One can thus improve the denoising results by thresholding block of coefficients togethers. Block thresholding is only efficient when used as a soft thresholder. Here we use a Stein soft thresholder.
Display the thresholded coefficients for a threshold value $T$ proportional to the noise level $\si$.
T = 1.25 * sigma
plt.figure(figsize=(10, 10))
plot_wavelet(ThreshBlock(wav(f), T), Jmin)
plt.show()
Define the wavelet block thresholding operator.
ThreshWav = lambda f, T: iwav(ThreshBlock(wav(f), T))
Test the thresholding.
plt.figure(figsize=(5, 5))
imageplot(clamp(ThreshWav(f, T)))
Worked example 3
Display the evolution of the denoising SNR when $T$ varies. Store the optimal denoising result in $f_{Block}$.
tlist = np.linspace(0.5, 2, 30) * sigma
snr_stein = [snr(f0, ThreshWav(f, t)) for t in Tlist]
plt.plot([T / sigma for T in Tlist], snr_stein, linewidth=2)
plt.xlabel(r"$T/\sigma$")
plt.ylabel("SNR")
plt.show()
Tmax = np.max(snr_stein)
fBlock = ThreshWav(f, T)
Display the result.
plt.figure(figsize=(5, 5))
imageplot(clamp(fBlock), "SNR = %.1f dB" % snr(f0, fBlock))
Translation invariant Block Thresholding¶
Block thresholding can also be applied to a translation invariant wavelet transform. It gives state of the art denoising results.
Shortcuts for the foward and backward translation invariant wavelet transforms.
wav = lambda f: perform_wavelet_transf(f, Jmin, +1, ti=1)
iwav = lambda fw: perform_wavelet_transf(fw, Jmin, -1, ti=1)
Foward wavelet transform.
fw = wav(f)
n = 256
# np.shape(fw)[0]
Compute indexing of the blocks.
[X, J, Y, dX, dY] = np.meshgrid(
np.arange(1, n - w + 2, w),
np.arange(1, np.shape(fw)[0] + 1),
np.arange(1, n - w + 2, w),
np.arange(0, w),
np.arange(0, w),
)
I = (X + dX - 1) + (Y + dY - 1) * n + (J - 1) * n**2
for k in range(n // w):
for l in range(n // w):
for m in range(np.shape(fw)[0]):
I[m][k][l] = np.transpose(I[m][k][l])
Forward and backward extraction operators.
block = lambda x: np.ravel(x)[I]
def assign(M, I, H):
M_temp = M
np.ravel(M_temp)[I] = H
return np.reshape(M_temp, (np.shape(fw)[0], n, n))
iblock = lambda H: assign(np.zeros([np.shape(fw)[0], n, n]), I, H)
Compute the average energy of each block, and duplicate.
def energy(H):
H_tmp = np.copy(H)
for i in range(n // w):
for j in range(n // w):
for k in range(np.shape(fw)[0]):
H_tmp[k][i][j] = np.sqrt(np.mean(H_tmp[k][i][j] ** 2))
return H_tmp
Block thresholding operator.
Thresh = lambda H, T: psi(energy(H), T) * H
ThreshBlock = lambda x, T: iblock(Thresh(block(x), T))
Define the wavelet block thresholding operator.
ThreshWav = lambda f, T: iwav(ThreshBlock(wav(f), T))
Test the thresholding.
T = 1.25 * sigma
plt.figure(figsize=(5, 5))
imageplot(clamp(ThreshWav(f, T)))
Worked example 4
Display the evolution of the denoising SNR when $T$ varies. Store the optimal denoising result in $f_{TI}$.
tlist = np.linspace(0.5, 2, 20) * sigma
snr_stein = [snr(f0, ThreshWav(f, t)) for t in Tlist]
plt.plot([T / sigma for T in Tlist], snr_stein, linewidth=2)
plt.xlabel("$T/sigma$")
plt.ylabel("SNR")
plt.show()
Tmax = np.max(snr_stein)
fTI = ThreshWav(f, T)
Display the result.
plt.figure(figsize=(5, 5))
imageplot(clamp(fTI), "SNR = %.1f dB" % snr(f0, fTI))
References and further reading¶
[CaiSil01] T. Cai and B.W. Silverman, [Incorporating information on neighboring coefficients into wavelet estimation][1], Sankhya 63, 127-148, 2001.
[Cai99] T. Cai, [Adaptive wavelet estimation: a block thresholding and oracle inequality approach][2], The Annals of Statistics 27, 898-924, 1999.
[HallKerkPic99] P. Hall, G. Kerkyacharian and D. Picard, On the minimax optimality of block thresholded wavelet estimator, Statistica Sinica 9(1999), 33-49
David L. Donoho. De-noising by Soft-Thresholding. 1995, IEEE Transactions on Information Theory 41(3), 613–627. Why shrinkage of wavelet coefficients suppresses noise.
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.
Stéphane Mallat. A Wavelet Tour of Signal Processing: The Sparse Way. 2009, 3rd ed., Academic Press. Multiresolution analysis, sparse approximation, and wavelet algorithms.
Ingrid Daubechies. Ten Lectures on Wavelets. 1992, SIAM. Compactly supported orthogonal wavelets and their regularity.