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 overviews the use of Fourier and wavelets for image approximation.
using PyPlot
using NtToolBox
# using Autoreload
# arequire("NtToolBox")
Note: to measure the error of an image $f$ with its approximation $f_M$, we use the SNR measure, defined as
$$ \text{SNR}(f,f_M) = -20\log_{10} \pa{ \frac{ \norm{f-f_M} }{ \norm{f} } }, $$which is a quantity expressed in decibels (dB). The higer the SNR, the better the quality.
First we load an image $ f \in \RR^N $ of $ N = N_0 \times N_0 $ pixels.
n0 = 512
f = rescale(load_image("NtToolBox/src/data/lena.png", n0));
Display the original image.
figure(figsize = (5,5))
imageplot(f, "Image f");
Display a zoom in the middle.
figure(figsize = (5,5))
imageplot(f[Int(n0/2 - 32) : Int(n0/2 + 32), Int(n0/2 - 32) : Int(n0/2 + 32)], "Zoom")
PyObject <matplotlib.text.Text object at 0x000000001E7DA240>
An image is a 2D array, it can be modified as a matrix.
figure(figsize = (8,8))
imageplot(-f, "-f", [1, 2, 1])
imageplot(f[end:-1:1, :], "Flipped", [1, 2, 2])
PyObject <matplotlib.text.Text object at 0x000000001EF47898>
Blurring is achieved by computing a convolution $f \star h$ with a kernel $h$.
Compute the low pass kernel.
k = 9; #size of the kernel
h = ones(k, k)
h = h/sum(h); #normalize
9×9 Array{Float64,2}: 0.0123457 0.0123457 0.0123457 … 0.0123457 0.0123457 0.0123457 0.0123457 0.0123457 0.0123457 0.0123457 0.0123457 0.0123457 0.0123457 0.0123457 0.0123457 0.0123457 0.0123457 0.0123457 0.0123457 0.0123457 0.0123457 0.0123457 0.0123457 0.0123457 0.0123457 0.0123457 0.0123457 0.0123457 0.0123457 0.0123457 0.0123457 0.0123457 0.0123457 … 0.0123457 0.0123457 0.0123457 0.0123457 0.0123457 0.0123457 0.0123457 0.0123457 0.0123457 0.0123457 0.0123457 0.0123457 0.0123457 0.0123457 0.0123457 0.0123457 0.0123457 0.0123457 0.0123457 0.0123457 0.0123457
Compute the convolution $f \star h$.
fh = conv2(Array{Float64, 2}(f), h);
Display.
figure(figsize = (5,5))
imageplot(fh, "Blurred image")
PyObject <matplotlib.text.Text object at 0x000000001F218C50>
The Fourier orthonormal basis is defined as $$ \psi_m(k) = \frac{1}{\sqrt{N}}e^{\frac{2i\pi}{N_0} \dotp{m}{k} } $$ where $0 \leq k_1,k_2 < N_0$ are position indexes, and $0 \leq m_1,m_2 < N_0$ are frequency indexes.
The Fourier transform $\hat f$ is the projection of the image on this Fourier basis
$$ \hat f(m) = \dotp{f}{\psi_m}. $$The Fourier transform is computed in $ O(N \log(N)) $ operation using the FFT algorithm (Fast Fourier Transform). Note the normalization by $\sqrt{N}=N_0$ to make the transform orthonormal.
F = plan_fft(f)
F = (F*f)/n0;
We check this conservation of the energy.
println(@sprintf("Energy of Image: %f", norm(f)))
println(@sprintf("Energy of Fourier: %f", norm(F)))
Energy of Image: 262.554108 Energy of Fourier: 262.554138
Compute the logarithm of the Fourier magnitude $ \log\left(\abs{\hat f(m)} + \epsilon\right) $, for some small $\epsilon$.
L = fftshift(log(abs(F) + 1e-1));
Display. Note that we use the function fftshift to put the 0 low frequency in the middle.
figure(figsize = (5,5))
imageplot(L, "Log(Fourier transform)")
PyObject <matplotlib.text.Text object at 0x00000000244FA978>
An approximation is obtained by retaining a certain set of index $I_M$
$$ f_M = \sum_{ m \in I_M } \dotp{f}{\psi_m} \psi_m. $$Linear approximation is obtained by retaining a fixed set $I_M$ of $M = \abs{I_M}$ coefficients. The important point is that $I_M$ does not depend on the image $f$ to be approximated.
For the Fourier transform, a low pass linear approximation is obtained by keeping only the frequencies within a square.
$$ I_M = \enscond{m=(m_1,m_2)}{ -q/2 \leq m_1,m_2 < q/2 } $$where $ q = \sqrt{M} $.
This can be achieved by computing the Fourier transform, setting to zero the $N-M$ coefficients outside the square $I_M$ and then inverting the Fourier transform.
Number $M$ of kept coefficients.
M = Base.div(n0^2, 64);
Exercise 1
Perform the linear Fourier approximation with $M$ coefficients. Store the result in the variable $f_M$.
include("NtSolutions/introduction_4_fourier_wavelets/exo1.jl")
PyObject <matplotlib.text.Text object at 0x0000000024A79908>
## Insert your code here.
Compare two 1D profile (lines of the image). This shows the strong ringing artifact of the linea approximation.
figure(figsize = (7, 6))
subplot(2, 1, 1)
plot(f[: , Base.div(n0, 2)])
xlim(0, n0)
title("f")
subplot(2, 1, 2)
plot(fM[: , Base.div(n0, 2)])
xlim(0, n0)
title("f_M")
show()
Non-linear approximation is obtained by keeping the $M$ largest coefficients. This is equivalently computed using a thresholding of the coefficients $$ I_M = \enscond{m}{ \abs{\dotp{f}{\psi_m}}>T }. $$
Set a threshold $T>0$.
T = .2;
Compute the Fourier transform.
F = plan_fft(f)
F = (F*f)/n0;
Do the hard thresholding.
FT = F .* (abs(F) .> T);
Display. Note that we use the function fftshift to put the 0 low frequency in the middle.
L = fftshift(log(abs(FT) + 1e-1))
figure(figsize = (5,5))
imageplot(L, "thresholded Log(Fourier transform)")
PyObject <matplotlib.text.Text object at 0x0000000024EC62E8>
Inverse Fourier transform to obtain $f_M$.
fM = plan_ifft(FT)
fM = real(n0*(fM*FT));
Display.
figure(figsize = (5,5))
imageplot(clamP(fM), @sprintf("Linear, Fourier, SNR = %.1f dB", snr(f, fM)))
PyObject <matplotlib.text.Text object at 0x0000000029BFF5C0>
Given a $T$, the number of coefficients is obtained by counting the non-thresholded coefficients $ \abs{I_M} $.
m = sum(FT .!= 0)
print(@sprintf("M/N = 1/%d" ,(n0^2)/m))
M/N = 1/31
Exercise 2
Compute the value of the threshold $T$ so that the number of coefficients is $M$. Display the corresponding approximation $f_M$.
include("NtSolutions/introduction_4_fourier_wavelets/exo2.jl")
PyObject <matplotlib.text.Text object at 0x0000000024D681D0>
## Insert your code here.
A wavelet basis $ \Bb = \{ \psi_m \}_m $ is obtained over the continuous domain by translating and dilating three mother wavelet functions $ \{\psi^V,\psi^H,\psi^D\} $.
Each wavelet atom is defined as $$ \psi_m(x) = \psi_{j,n}^k(x) = \frac{1}{2^j}\psi^k\pa{ \frac{x-2^j n}{2^j} } $$
The scale (size of the support) is $2^j$ and the position is $2^j(n_1,n_2)$. The index is $ m=(k,j,n) $ for $\{ j \leq 0 \}$.
The wavelet transform computes all the inner products $ \{ \dotp{f}{\psi_{j,n}^k} \}_{k,j,n} $.
Set the minimum scale for the transform to be 0.
Jmin = 0;
Perform the wavelet transform, $f_w$ stores all the wavelet coefficients.
fw = NtToolBox.perform_wavelet_transf(f, Jmin, + 1);
# using NPZ
# test = npzread("pgsh.npy")
512×512 Array{Float64,2}: 233.032 19.1135 -13.0311 … 0.0418207 -0.131423 -23.1363 1.65645 -14.2758 0.0413636 -0.128621 -4.93937 15.4845 1.19021 0.0394703 -0.133461 -5.88592 -3.01506 -3.32998 -0.00103448 -0.0385822 0.922625 -2.16059 0.014336 -0.0159323 -0.00210343 -0.335541 5.64581 -9.66819 … -0.00869429 0.000518759 -0.518646 -0.502312 5.0826 0.00102269 0.0117643 4.98911 -5.50287 9.23286 -0.00281264 -0.00380574 0.964176 0.19277 0.00966503 -0.0331573 0.012915 -3.20033 0.128032 -0.512992 0.00189917 0.0302912 0.319147 0.377833 0.160483 … 0.0115657 0.0148096 -0.273935 -0.155208 0.166632 0.00603326 0.0116551 -0.603172 0.766374 -1.38124 -0.00246024 0.0250189 ⋮ ⋱ ⋮ -0.00091898 -0.0039071 -0.0059023 … -0.0110338 -0.00267938 0.0181912 -0.00505674 0.00412992 -0.0279641 -0.0170076 -0.0186527 -0.0144188 -0.0191199 0.00409298 -0.00667346 -0.0275172 0.00556199 -0.0102931 -0.0184985 0.000831558 0.0513395 -0.00900993 0.00858931 0.018305 0.0116154 0.0127468 -0.0135157 0.027699 … -0.00627344 0.0268372 0.00047219 -0.0162774 -0.0247788 -0.00592464 -0.0204156 0.00738821 -4.32492e-5 -0.00901067 0.011207 0.0168448 -0.0102837 0.00511243 -0.00709129 -0.00368487 -0.0161069 0.0455343 -0.0165399 -0.000882116 0.00323547 -0.0198884 -0.010534 -0.0135025 -0.00207893 … 0.00426511 -0.00434588 2.50898e-5 0.0201377 0.0201784 -3.24896e-5 -0.000651598
Display the transformed coefficients.
figure(figsize = (10, 10))
NtToolBox.plot_wavelet(fw)
title("Wavelet coefficients")
PyObject <matplotlib.text.Text object at 0x0000000024C7FE80>
Linear wavelet approximation with $M=2^{-j_0}$ coefficients is obtained by keeping only the coarse scale (large support) wavelets:
$$ I_M = \enscond{(k,j,n)}{ j \geq j_0 }. $$It corresponds to setting to zero all the coefficients excepted those that are on the upper left corner of $f_w$.
Exercise 3
Perform linear approximation with $M$ wavelet coefficients.
include("NtSolutions/introduction_4_fourier_wavelets/exo3.jl")
PyObject <matplotlib.text.Text object at 0x000000002AB4D8D0>
## Insert your code here.
A non-linear approximation is obtained by keeping the $M$ largest wavelet coefficients.
As already said, this is equivalently computed by a non-linear hard thresholding.
Select a threshold.
T = .15
0.15
Perform hard thresholding.
fwT = fw.*(abs(fw) .> T);
512×512 Array{Float32,2}: 260.699 17.7005 -11.932 … -0.0 -0.0 0.0 -0.0 0.0 -0.0 -24.4527 1.57447 -14.2721 -0.0 -0.0 0.0 -0.0 0.0 -0.0 -4.41654 14.5837 1.60971 -0.0 -0.0 0.0 -0.0 0.0 -0.0 -6.88045 -2.49392 -3.37559 -0.0 -0.0 -0.0 0.0 0.0 -0.0 0.850178 -2.0469 -0.165711 0.0 0.0 -0.0 -0.0 -0.0 -0.0 -0.334239 5.67381 -9.23927 … 0.0 0.0 -0.0 0.0 -0.0 0.0 -0.45224 -0.560716 5.35115 -0.0 -0.0 -0.0 -0.0 0.0 0.0 4.92688 -5.87634 9.30895 -0.0 0.0 -0.0 -0.0 -0.0 -0.0 0.93541 0.201458 0.0 -0.0 0.0 0.0 0.0 -0.0 0.0 -3.12512 0.0 -0.473662 -0.0 -0.0 -0.0 -0.0 0.0 0.0 0.295174 0.361619 0.168442 … 0.0 -0.0 -0.0 0.0 0.0 0.0 -0.26721 -0.163369 0.348265 0.0 0.0 -0.0 0.0 0.0 0.0 -0.637094 0.747106 -1.44393 0.0 0.0 0.0 -0.0 -0.0 0.0 ⋮ ⋱ ⋮ -0.0 -0.0 -0.0 … 0.0 0.0 0.0 0.0 -0.0 -0.0 0.0 -0.0 0.0 0.0 -0.0 -0.0 0.0 -0.0 -0.0 -0.0 -0.0 -0.0 0.0 -0.0 -0.0 -0.0 0.0 -0.0 -0.0 0.0 -0.0 0.0 0.0 -0.0 0.0 -0.0 0.0 0.0 -0.0 0.0 -0.0 0.0 -0.0 -0.0 0.0 0.0 0.0 -0.0 0.0 … -0.0 0.0 -0.0 -0.0 -0.0 0.0 -0.0 -0.0 -0.0 0.0 -0.0 -0.0 0.0 -0.0 -0.0 0.0 0.0 -0.0 0.0 -0.0 -0.0 0.0 0.0 0.0 -0.0 0.0 -0.0 0.0 0.0 -0.0 0.0 -0.0 -0.0 0.0 -0.0 -0.0 -0.0 0.0 0.0 0.0 0.0 -0.0 -0.0 -0.0 -0.0 … -0.0 -0.0 -0.0 -0.0 0.0 -0.0 -0.0 0.0 0.0 0.0 0.0 0.0 -0.0 -0.0 -0.0
Display the thresholded coefficients.
figure(figsize = (15,15))
subplot(1, 2, 1)
plot_wavelet(fw)
title("Original coefficients")
subplot(1, 2, 2)
plot_wavelet(fwT)
title("Thresholded coefficients")
show()
Perform reconstruction.
fM = NtToolBox.perform_wavelet_transf(fwT, Jmin, -1)
512×512 Array{Float32,2}: 0.686229 0.685407 0.683721 0.681588 … 0.662526 0.615727 0.606579 0.685079 0.684285 0.682657 0.6806 0.673313 0.618527 0.607662 0.682756 0.682019 0.68051 0.678606 0.678797 0.612439 0.599023 0.679865 0.6792 0.677839 0.676126 0.672799 0.596129 0.58039 0.676523 0.675941 0.674754 0.673262 0.633161 0.552921 0.536134 0.673179 0.672681 0.671667 0.670399 … 0.510544 0.446753 0.432775 0.669654 0.669246 0.668416 0.667383 0.36318 0.322226 0.312339 0.665888 0.665576 0.664943 0.664164 0.215849 0.197916 0.192186 0.66205 0.661835 0.661404 0.660886 0.0984763 0.0994317 0.0972883 0.658284 0.658166 0.657935 0.657673 0.0765254 0.0812711 0.0803829 0.654587 0.654564 0.654531 0.654523 … 0.0905613 0.0935602 0.0930404 0.651038 0.651108 0.651266 0.651503 0.111606 0.112075 0.111789 0.64759 0.64775 0.648096 0.648574 0.132954 0.130535 0.130229 ⋮ ⋱ ⋮ 0.147979 0.148415 0.150817 0.155315 … 0.222486 0.223004 0.223461 0.143876 0.145378 0.149306 0.155092 0.231788 0.232565 0.233201 0.141104 0.143624 0.148991 0.15598 0.240186 0.241207 0.242012 0.138028 0.141741 0.148859 0.15739 0.249573 0.250832 0.251801 0.135549 0.140233 0.148773 0.158573 0.264317 0.265806 0.266919 0.134055 0.138998 0.147884 0.157995 … 0.294538 0.296256 0.297475 0.133131 0.138055 0.146874 0.156933 0.328499 0.330425 0.331727 0.132573 0.137414 0.14608 0.156022 0.360417 0.362523 0.36389 0.132268 0.136994 0.145464 0.15525 0.386299 0.388543 0.389955 0.131977 0.136721 0.145208 0.155011 0.395695 0.397999 0.399436 0.131802 0.136584 0.145123 0.154972 … 0.398126 0.400452 0.4019 0.131744 0.136531 0.145076 0.15493 0.399811 0.402148 0.403601
Display approximation.
figure(figsize = (5, 5))
imageplot(clamP(fM), @sprintf("Approximation, SNR, = %.1f dB", snr(f, fM)))
PyObject <matplotlib.text.Text object at 0x0000000031911080>
Exercise 4
Perform non-linear approximation with $M$ wavelet coefficients by chosing the correct value for $T$. Store the result in the variable $f_M$.
include("NtSolutions/introduction_4_fourier_wavelets/exo4.jl")
PyObject <matplotlib.text.Text object at 0x00000000322BCDA0>
## Insert your code here.
Compare two 1D profile (lines of the image). Note how the ringing artifacts are reduced compared to the Fourier approximation.
figure(figsize = (7, 6))
subplot(2, 1, 1)
plot(f[:, Base.div(n0, 2)])
title("f")
subplot(2, 1, 2)
plot(fM[:, Base.div(n0, 2)])
title("f_M")
show()