Important: Please read the installation page for details about how to install the toolboxes. $\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{\umin}[1]{\underset{#1}{\min}\;}$ $\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{\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
using PyPlot
using NtToolBox
# using Autoreload
# arequire("NtToolBox")
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("NtToolBox/src/data/boat.png", n));
Display it.
figure(figsize = (5,5))
imageplot(f0)
Noise level.
sigma = .08;
Generate a noisy image $f=f_0+\epsilon$ where $\epsilon \sim \Nn(0,\si^2\text{Id}_N)$.
using Distributions
f = f0 + sigma.*rand(Normal(), n, n);
Display it.
figure(figsize = (5,5))
imageplot(clamP(f))
We first consider the traditional wavelet thresholding method.
Parameters for the orthogonal wavelet transform.
Jmin = 4;
Shortcuts for the foward and backward wavelet transforms.
wav = f -> NtToolBox.perform_wavelet_transf(f, Jmin, +1)
iwav = fw -> NtToolBox.perform_wavelet_transf(fw, Jmin, -1);
Display the original set of noisy coefficients.
figure(figsize = (10, 10))
plot_wavelet(wav(f), Jmin)
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 = (s, T) -> max(1 - T^2./max(abs(s).^2, 1e-9.*ones(size(s))), zeros(size(s)));
Display the thresholding function $\th_T$.
s = linspace(-3, 3, 1024)
plot(s, s.*psi(s, 1))
plot(s, s, "r--")
show()
Thresholding operator.
theta = (x, T) -> psi(x, T).*x
ThreshWav = (f, T) -> iwav(theta(wav(f), T));
Test the thresholding.
T = 1.5*sigma
figure(figsize = (5, 5))
imageplot(clamP(ThreshWav(f, T)))
Exercise 1
Display the evolution of the denoising SNR when $T$ varies. Store in $f_{Thresh}$ the optimal denoising result.
include("NtSolutions/denoisingwav_4_block/exo1.jl")
256×256 Array{Float64,2}: 0.625968 0.694511 0.664028 0.685825 … 0.587188 0.58559 0.556904 0.666833 0.713415 0.698709 0.719822 0.588342 0.582681 0.582041 0.687508 0.723447 0.757062 0.723292 0.584483 0.58857 0.618832 0.852152 0.727595 0.629546 0.656954 0.629153 0.629052 0.646233 0.6322 0.617701 0.726676 0.7196 0.673415 0.663018 0.661216 0.568712 0.733649 0.651538 0.684211 … 0.657218 0.644178 0.645722 0.705647 0.644722 0.693751 0.672634 0.630126 0.610606 0.614323 0.682008 0.659418 0.667178 0.680056 0.565555 0.605219 0.578848 0.663885 0.678039 0.67714 0.650876 0.583889 0.570246 0.527736 0.688729 0.649929 0.651663 0.739553 0.592867 0.606361 0.58689 0.715614 0.633788 0.634491 0.801354 … 0.580648 0.553483 0.55032 0.696123 0.666799 0.706528 0.635616 0.535158 0.571715 0.581837 0.685458 0.682204 0.665456 0.725485 0.598165 0.580181 0.576446 ⋮ ⋱ ⋮ 0.636211 0.488664 0.280411 0.418032 0.528137 0.556917 0.515754 0.569564 0.605694 0.741009 0.75532 … 0.582789 0.527077 0.473523 0.514487 0.604099 0.664932 0.647077 0.512832 0.506804 0.456944 0.631128 0.666075 0.698192 0.704423 0.356548 0.510958 0.480715 0.765495 0.714023 0.667206 0.721657 0.495711 0.436455 0.445046 0.718283 0.718393 0.717137 0.689142 0.561636 0.498873 0.468732 0.757006 0.747529 0.738718 0.753412 … 0.631288 0.560233 0.498495 0.710838 0.714638 0.715273 0.722023 0.598957 0.586513 0.437358 0.672412 0.678607 0.68473 0.693114 0.528807 0.546961 0.507658 0.672463 0.673304 0.689778 0.692533 0.545252 0.493636 0.645247 0.686829 0.671221 0.702466 0.699128 0.529137 0.63748 0.435184 0.643881 0.713313 0.683531 0.696336 … 0.489809 0.560741 0.183504
## Insert your code here.
Display the optimal thresolding.
figure(figsize = (5, 5))
imageplot(clamP(fThresh), @sprintf("SNR = %.1f dB", snr(f0, fThresh)))
PyObject <matplotlib.text.Text object at 0x00000000292193C8>
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.
include("NtToolBox/src/ndgrid.jl")
(X, Y, dX, dY) = ndgrid(1:w:n-w+1, 1:w:n-w+1, 0:w-1, 0:w-1)
I = (X + dX) + (Y + dY - 1).*n
for k in 1:Base.div(n, w)
for l in 1:Base.div(n, w)
I[k, l, :, :] = transpose(I[k, l, :, :])
end
end
Block extraction operator. It returns the set $ \{x_b\}_{b \in B} $ of block-partitioned coefficients.
block = x -> x[I];
Block reconstruction operator.
function assign(M,I,H)
M_temp = M
M_temp[I] = H
return reshape(M_temp, n, n)
end
iblock = H -> assign(zeros(n, n), I, H);
Check that block extraction / reconstruction gives perfect reconstruction.
print(string("Should be 0:", string(norm(f - iblock(block(f))))))
Should be 0:0.0
Compute the average energy of each block, and duplicate.
function energy(H)
H_tmp = copy(H)
for i in 1:Base.div(n, w)
for j in 1:Base.div(n, w)
H_tmp[i, j, :, :] = sqrt(mean(H_tmp[i, j, :, :].^2))#*np.ones([1,1])
end
end
return H_tmp
end;
Block thresholding operator.
Thresh = (H, T) -> psi(energy(H), T).*H
ThreshBlock = (x, T) -> iblock(Thresh(block(x), T));
Exercise 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.
include("NtSolutions/denoisingwav_4_block/exo2.jl")
## Insert your code here.
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
figure(figsize = (10, 10))
plot_wavelet(ThreshBlock(wav(f), T), Jmin)
show()
Define the wavelet block thresholding operator.
ThreshWav = (f, T) -> iwav(ThreshBlock(wav(f), T));
Test the thresholding.
figure(figsize = (5, 5))
imageplot(clamP(ThreshWav(f, T)))
Exercise 3
Display the evolution of the denoising SNR when $T$ varies. Store the optimal denoising result in $f_{Block}$.
include("NtSolutions/denoisingwav_4_block/exo3.jl")
256×256 Array{Float64,2}: 0.708276 0.716079 0.721681 0.705112 … 0.594622 0.594432 0.594191 0.689898 0.685507 0.685854 0.707782 0.594116 0.593946 0.593714 0.676101 0.698417 0.722323 0.707287 0.593171 0.593031 0.59281 0.735642 0.699185 0.659755 0.681072 0.592193 0.592076 0.591864 0.658546 0.679994 0.705833 0.702673 0.591 0.590899 0.59069 0.685041 0.684093 0.68207 0.677891 … 0.589333 0.589235 0.589019 0.671367 0.681083 0.691338 0.682948 0.58778 0.587669 0.587437 0.671851 0.658516 0.649705 0.681245 0.586528 0.586348 0.586087 0.668965 0.675646 0.681652 0.671425 0.586267 0.586041 0.585747 0.670751 0.672871 0.675187 0.673316 0.58829 0.588204 0.587895 0.673557 0.673727 0.674567 0.676771 … 0.591353 0.591355 0.591016 0.677324 0.677102 0.677436 0.679677 0.594976 0.595035 0.594666 0.681583 0.681799 0.682248 0.682823 0.598643 0.598275 0.59778 ⋮ ⋱ ⋮ 0.591496 0.514052 0.434554 0.44353 0.498052 0.496636 0.481191 0.611982 0.627389 0.660556 0.643108 … 0.506399 0.506142 0.507296 0.616264 0.609159 0.619226 0.602805 0.529565 0.519757 0.523831 0.651678 0.677453 0.720804 0.724378 0.477549 0.499257 0.491822 0.73234 0.710128 0.687241 0.707617 0.485855 0.449755 0.470535 0.729241 0.725592 0.720386 0.725833 0.539397 0.508234 0.503808 0.710149 0.721543 0.732099 0.721406 … 0.584279 0.548427 0.523942 0.68637 0.695571 0.7063 0.707024 0.533141 0.541157 0.494518 0.662679 0.666339 0.67409 0.689402 0.482314 0.494347 0.497923 0.651602 0.65748 0.668284 0.684492 0.499455 0.470907 0.499921 0.646073 0.654807 0.668746 0.683093 0.497069 0.523742 0.473121 0.643925 0.653195 0.667823 0.682205 … 0.4853 0.519448 0.423335
## Insert your code here.
Display the result.
figure(figsize = (5, 5))
imageplot(clamP(fBlock), @sprintf("SNR = %.1f dB", snr(f0, fBlock)))
PyObject <matplotlib.text.Text object at 0x000000002BDAAEF0>
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 = f -> NtToolBox.perform_wavelet_transf(f, Jmin, + 1, "9-7", 0, 1)
iwav = fw -> NtToolBox.perform_wavelet_transf(fw, Jmin, -1,"9-7", 0, 1);
Foward wavelet transform.
fw = wav(f)
n = 256;
256
Compute indexing of the blocks.
(X, Y, J, dX, dY) = ndgrid(1 : w : n - w + 1, 1 : w : n - w + 1, 1 : size(fw)[3], 0 : w - 1, 0 : w - 1)
I = (X + dX) + (Y + dY - 1).*n + (J - 1).*n^2
for k in 1 : Base.div(n, w)
for l in 1 : Base.div(n, w)
for m in 1 : size(fw)[3]
I[k, l, m, :, :] = transpose(I[k, l, m, :, :])
end
end
end
Forward and backward extraction operators.
block = x -> x[I]
function assign(M, I, H)
M_temp = M
M_temp[I] = H
return reshape(M_temp,n, n, size(fw)[3])
end
iblock = H -> assign(zeros(n, n, size(fw)[3]), I, H);
WARNING: Method definition assign(Any, Any, Any) in module Main at In[20]:2 overwritten at In[34]:5.
Compute the average energy of each block, and duplicate.
function energy(H)
H_tmp = copy(H)
for i in 1 : Base.div(n, w)
for j in 1 : Base.div(n, w)
for k in 1 : size(fw)[3]
H_tmp[i, j, k, :, :] = sqrt(mean(H_tmp[i, j, k, :, :].^2))
end
end
end
return H_tmp
end;
WARNING: Method definition energy(Any) in module Main at In[22]:2 overwritten at In[35]:2.
Block thresholding operator.
Thresh = (H, T) -> psi(energy(H), T).*H
ThreshBlock = (x, T) -> iblock(Thresh(block(x), T));
Define the wavelet block thresholding operator.
ThreshWav = (f, T) -> iwav(ThreshBlock(wav(f), T));
Test the thresholding.
T = 1.25*sigma
figure(figsize = (5, 5))
imageplot(clamP(ThreshWav(f, T)))
Exercise 4
Display the evolution of the denoising SNR when $T$ varies. Store the optimal denoising result in $f_{TI}$.
include("NtSolutions/denoisingwav_4_block/exo4.jl")
# tlist = linspace(.5, 2, 20).*sigma
# snr_stein = collect(snr(f0, ThreshWav(f,t)) for t in Tlist)
# plot(collect(T/sigma for T in Tlist), snr_stein, linewidth = 2)
# xlabel(L"$T/sigma$")
# ylabel("SNR")
# show()
# Tmax = maximum(snr_stein)
# fTI = ThreshWav(f, T)
256×256 Array{Float64,2}: 0.705959 0.708559 0.707965 0.702877 … 0.603356 0.58059 0.565643 0.694658 0.695238 0.69265 0.701108 0.596645 0.588432 0.582934 0.694262 0.705138 0.713384 0.706671 0.586699 0.602027 0.613183 0.708009 0.700663 0.685735 0.693351 0.606078 0.62087 0.630828 0.702548 0.681536 0.70551 0.694023 0.634547 0.633418 0.632037 0.673386 0.707069 0.683062 0.694029 … 0.62847 0.622934 0.619181 0.693209 0.680631 0.696625 0.686612 0.595181 0.599137 0.602777 0.687071 0.684781 0.679226 0.687507 0.579531 0.590985 0.599508 0.675471 0.682667 0.680491 0.679325 0.591727 0.582115 0.576856 0.682105 0.674623 0.673461 0.687863 0.595415 0.611112 0.619196 0.693863 0.662993 0.675944 0.696263 … 0.588999 0.568881 0.560253 0.695392 0.663971 0.682741 0.693498 0.578174 0.590491 0.595126 0.6893 0.678009 0.686724 0.686026 0.582681 0.57764 0.575819 ⋮ ⋱ ⋮ 0.624488 0.433774 0.358613 0.41849 0.511859 0.483121 0.452036 0.701078 0.575388 0.723412 0.718359 … 0.515765 0.486101 0.465175 0.625001 0.569137 0.662932 0.585952 0.493854 0.492315 0.489494 0.651638 0.699196 0.761278 0.728811 0.454787 0.490745 0.504276 0.708001 0.721193 0.694245 0.689323 0.508013 0.453997 0.503434 0.732676 0.723376 0.701544 0.695993 0.56831 0.457493 0.471463 0.724078 0.718564 0.711483 0.703552 … 0.567619 0.508575 0.436394 0.703213 0.704356 0.706554 0.70648 0.525425 0.531478 0.441264 0.680183 0.681951 0.686566 0.691357 0.488311 0.502686 0.528685 0.660508 0.662877 0.66909 0.676029 0.52151 0.421022 0.477251 0.646613 0.649947 0.65867 0.669151 0.480349 0.569425 0.404329 0.641423 0.645237 0.655289 0.667905 … 0.413961 0.534786 0.180515
## Insert your code here.
Display the result.
figure(figsize = (5,5))
imageplot(clamP(fTI), @sprintf("SNR = %.1f dB", snr(f0, fTI)));
PyObject <matplotlib.text.Text object at 0x000000002BA152B0>