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")
WARNING: Method definition ndgrid(AbstractArray{T<:Any, 1}) in module NtToolBox at /Users/gpeyre/.julia/v0.5/NtToolBox/src/ndgrid.jl:3 overwritten at /Users/gpeyre/.julia/v0.5/NtToolBox/src/ndgrid.jl:3. WARNING: Method definition ndgrid(AbstractArray{#T<:Any, 1}, AbstractArray{#T<:Any, 1}) in module NtToolBox at /Users/gpeyre/.julia/v0.5/NtToolBox/src/ndgrid.jl:6 overwritten at /Users/gpeyre/.julia/v0.5/NtToolBox/src/ndgrid.jl:6. WARNING: Method definition ndgrid_fill(Any, Any, Any, Any) in module NtToolBox at /Users/gpeyre/.julia/v0.5/NtToolBox/src/ndgrid.jl:13 overwritten at /Users/gpeyre/.julia/v0.5/NtToolBox/src/ndgrid.jl:13. WARNING: Method definition ndgrid(AbstractArray{#T<:Any, 1}...) in module NtToolBox at /Users/gpeyre/.julia/v0.5/NtToolBox/src/ndgrid.jl:19 overwritten at /Users/gpeyre/.julia/v0.5/NtToolBox/src/ndgrid.jl:19. WARNING: Method definition meshgrid(AbstractArray{T<:Any, 1}) in module NtToolBox at /Users/gpeyre/.julia/v0.5/NtToolBox/src/ndgrid.jl:33 overwritten at /Users/gpeyre/.julia/v0.5/NtToolBox/src/ndgrid.jl:33. WARNING: Method definition meshgrid(AbstractArray{#T<:Any, 1}, AbstractArray{#T<:Any, 1}) in module NtToolBox at /Users/gpeyre/.julia/v0.5/NtToolBox/src/ndgrid.jl:36 overwritten at /Users/gpeyre/.julia/v0.5/NtToolBox/src/ndgrid.jl:36. WARNING: Method definition meshgrid(AbstractArray{#T<:Any, 1}, AbstractArray{#T<:Any, 1}, AbstractArray{#T<:Any, 1}) in module NtToolBox at /Users/gpeyre/.julia/v0.5/NtToolBox/src/ndgrid.jl:44 overwritten at /Users/gpeyre/.julia/v0.5/NtToolBox/src/ndgrid.jl:44.
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")
PyObject <matplotlib.text.Text object at 0x32561aa10>
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 0x325746250>
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 0x3264d2490>
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
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 0x326df8050>
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 0x32aa4d050>
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$.
#run -i nt_solutions/introduction_4_fourier_wavelets/exo1
include("NtSolutions\\introduction_4_fourier_wavelets\\exo1.jl")
could not open file /Users/gpeyre/Dropbox/github/numerical-tours/julia/NtSolutions\introduction_4_fourier_wavelets\exo1.jl in include_from_node1(::String) at ./loading.jl:488 in include_from_node1(::String) at /Applications/Julia-0.5.app/Contents/Resources/julia/lib/julia/sys.dylib:?
## 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()
UndefVarError: fM not defined
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 0x32afa2150>
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 0x3257dac50>
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")
could not open file /Users/gpeyre/Dropbox/github/numerical-tours/julia/NtSolutions\introduction_4_fourier_wavelets\exo2.jl in include_from_node1(::String) at ./loading.jl:488 in include_from_node1(::String) at /Applications/Julia-0.5.app/Contents/Resources/julia/lib/julia/sys.dylib:?
## 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")
Display the transformed coefficients.
figure(figsize = (10, 10))
NtToolBox.plot_wavelet(fw)
title("Wavelet coefficients")
PyObject <matplotlib.text.Text object at 0x32af8a810>
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")
could not open file /Users/gpeyre/Dropbox/github/numerical-tours/julia/NtSolutions\introduction_4_fourier_wavelets\exo3.jl in include_from_node1(::String) at ./loading.jl:488 in include_from_node1(::String) at /Applications/Julia-0.5.app/Contents/Resources/julia/lib/julia/sys.dylib:?
## 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);
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 0x33055f150>
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")
could not open file /Users/gpeyre/Dropbox/github/numerical-tours/julia/NtSolutions\introduction_4_fourier_wavelets\exo4.jl in include_from_node1(::String) at ./loading.jl:488 in include_from_node1(::String) at /Applications/Julia-0.5.app/Contents/Resources/julia/lib/julia/sys.dylib:?
## 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()