# Image Approximation with Orthogonal Bases¶

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 uses several orthogonal bases to perform non-linear image approximation.

In [74]:
using PyPlot
using NtToolBox
arequire("NtToolBox")

WARNING: using NtToolBox.compute_wavelet_filter in module Main conflicts with an existing identifier.
WARNING: using NtToolBox.perform_wavortho_transf in module Main conflicts with an existing identifier.
WARNING: using NtToolBox.imageplot in module Main conflicts with an existing identifier.
WARNING: using NtToolBox.div in module Main conflicts with an existing identifier.
WARNING: using NtToolBox.snr in module Main conflicts with an existing identifier.
WARNING: using NtToolBox.rescale in module Main conflicts with an existing identifier.
WARNING: using NtToolBox.load_image in module Main conflicts with an existing identifier.
WARNING: using NtToolBox.plot_wavelet in module Main conflicts with an existing identifier.
WARNING: using NtToolBox.clamP in module Main conflicts with an existing identifier.
WARNING: replacing module NtToolBox


## Best $M$-terms Non-linear Approximation¶

This tours makes use of an orthogonal base $\Bb = \{ \psi_m \}_{m=0}^{N-1}$ of the space $\RR^N$ of the images with $N$ pixels.

The best $M$-term approximation of $f$ is obtained by a non-linear thresholding

$$f_M = \sum_{ \abs{\dotp{f}{\psi_m}}>T } \dotp{f}{\psi_m} \psi_m,$$

where the value of $T>0$ should be carefully selected so that only $M$ coefficients are not thresholded, i.e.

$$\abs{ \enscond{m}{ \abs{\dotp{f}{\psi_m}}>T } } = M.$$

The goal is to use an ortho-basis $\Bb$ so that the error $\norm{f-f_M}$ decays as fast as possible when $M$ increases, for a large class of images.

This tour studies several different orthogonal bases: Fourier, wavelets (which is at the heart of JPEG-2000), cosine, local cosine (which is at the heart of JPEG).

First we load an image of $N = n \times n$ pixels.

In [75]:
n = 512

Out[75]:
512×512 Array{Float32,2}:
0.672897  0.672897  0.672897  0.669782  …  0.705608  0.64486   0.528037
0.672897  0.672897  0.672897  0.669782     0.705608  0.64486   0.528037
0.672897  0.672897  0.672897  0.669782     0.705608  0.64486   0.528037
0.672897  0.672897  0.672897  0.669782     0.705608  0.64486   0.528037
0.672897  0.672897  0.672897  0.669782     0.705608  0.64486   0.528037
0.682243  0.682243  0.658879  0.64486   …  0.514019  0.415888  0.291277
0.668224  0.668224  0.67757   0.658879     0.303738  0.224299  0.196262
0.663551  0.663551  0.64486   0.654206     0.149533  0.17757   0.135514
0.64486   0.64486   0.658879  0.658879     0.126168  0.135514  0.14486
0.64486   0.64486   0.654206  0.658879     0.135514  0.126168  0.116822
0.649533  0.649533  0.649533  0.668224  …  0.126168  0.154206  0.154206
0.649533  0.649533  0.649533  0.663551     0.126168  0.121495  0.140187
0.658879  0.658879  0.654206  0.649533     0.135514  0.140187  0.149533
⋮                                       ⋱            ⋮
0.163551  0.163551  0.182243  0.200935  …  0.17757   0.168224  0.154206
0.149533  0.149533  0.163551  0.182243     0.182243  0.228972  0.224299
0.140187  0.140187  0.182243  0.219626     0.200935  0.257009  0.317757
0.14486   0.14486   0.172897  0.17757      0.247664  0.280374  0.364486
0.140187  0.140187  0.168224  0.14486      0.261682  0.299065  0.336449
0.135514  0.135514  0.17757   0.168224  …  0.317757  0.359813  0.373832
0.140187  0.140187  0.186916  0.163551     0.35514   0.345794  0.378505
0.172897  0.172897  0.140187  0.135514     0.369159  0.373832  0.35514
0.126168  0.126168  0.149533  0.149533     0.411215  0.392523  0.383178
0.11215   0.11215   0.149533  0.135514     0.420561  0.401869  0.392523
0.116822  0.116822  0.17757   0.154206  …  0.420561  0.425234  0.439252
0.116822  0.116822  0.17757   0.238318     0.420561  0.425234  0.439252

Display it.

In [76]:
figure(figsize = (5,5))
imageplot(f)


## Fourier Approximation¶

The discrete 2-D Fourier atoms are defined as: $$\psi_m(x) = \frac{1}{\sqrt{N}} e^{ \frac{2i\pi}{n} ( x_1 m_1 + x_2 m_2 ) },$$ where $0 \leq m_1,m_2 < n$ indexes the frequency.

The set of inner products $\{ \dotp{f}{\psi_m} \}_m$ is computed in $O(N \log(N))$ operations with the 2-D Fast Fourier Transform (FFT) algorithm (the pylab function is fft2).

Compute the Fourier transform using the FFT algorithm. Note the normalization by $1/\sqrt{N}$ to ensure orthogonality (energy conservation) of the transform.

In [77]:
fF = (plan_fft(f)*f)/n;

Out[77]:
512×512 Array{Complex{Float32},2}:
256.476+0.0im        -3.95615+20.4749im   …   -3.95615-20.4749im
0.0844855-11.1962im    -15.1921+13.9858im       0.630981+7.94724im
-3.55624-1.76501im    -2.36528-0.595752im       7.28666-0.598839im
4.73903+1.55932im   -0.805839-5.27302im        6.26791-6.04117im
-0.624841-5.30147im     1.55698-1.17529im      -0.988565+0.720974im
-1.61084-4.77705im    0.165628-1.66687im   …   -2.28491+1.02422im
1.22427-0.125188im   -2.10585-2.96785im        1.15793-2.06182im
2.86644-0.98837im   -0.248534+1.02317im       -0.55987-1.09585im
-0.569311-0.075523im   0.416833-0.219624im       1.58716-0.4015im
0.993204-1.84695im    -1.48726-0.244343im      -1.33509+0.742391im
-0.0138412-0.112531im  -0.682093-1.73745im   …    1.08995-0.3994im
0.83187-1.05805im    0.247246+0.45962im       0.700845-0.604206im
-0.225175+0.184748im  -0.398911-0.745511im       1.23564-0.586102im
⋮                                  ⋱
-0.225175-0.184748im    1.23564+0.586102im  …  -0.398911+0.745511im
0.83187+1.05805im    0.700845+0.604206im      0.247246-0.45962im
-0.0138412+0.112531im    1.08995+0.3994im       -0.682093+1.73745im
0.993204+1.84695im    -1.33509-0.742392im      -1.48726+0.244343im
-0.569311+0.075523im    1.58716+0.4015im        0.416833+0.219624im
2.86644+0.98837im    -0.55987+1.09585im   …  -0.248534-1.02317im
1.22427+0.125188im    1.15793+2.06182im       -2.10585+2.96785im
-1.61084+4.77705im    -2.28491-1.02422im       0.165628+1.66687im
-0.624841+5.30147im   -0.988565-0.720973im       1.55698+1.17529im
4.73903-1.55932im     6.26791+6.04117im      -0.805839+5.27302im
-3.55624+1.76501im     7.28667+0.598839im  …   -2.36528+0.595752im
0.0844855+11.1962im    0.630981-7.94724im       -15.1921-13.9858im 

Display its magnitude (in log scale). We use the pylab function fftshift to put the low frequency in the center.

In [78]:
figure(figsize = (5,5))
imageplot(log(1e-5 + abs(fftshift(fF))))


An image is recovered from a set of coefficients $c_m$ using the inverse Fourier Transform (pylab function ifft2)) that implements the formula

$$f_M = \sum_m c_m \psi_m.$$

Perform a thresholding.

In [79]:
T = .3
c = fF.*(abs(fF) .> T);

Out[79]:
512×512 Array{Complex{Float32},2}:
256.476+0.0im        -3.95615+20.4749im   …   -3.95615-20.4749im
0.0844855-11.1962im    -15.1921+13.9858im       0.630981+7.94724im
-3.55624-1.76501im    -2.36528-0.595752im       7.28666-0.598839im
4.73903+1.55932im   -0.805839-5.27302im        6.26791-6.04117im
-0.624841-5.30147im     1.55698-1.17529im      -0.988565+0.720974im
-1.61084-4.77705im    0.165628-1.66687im   …   -2.28491+1.02422im
1.22427-0.125188im   -2.10585-2.96785im        1.15793-2.06182im
2.86644-0.98837im   -0.248534+1.02317im       -0.55987-1.09585im
-0.569311-0.075523im   0.416833-0.219624im       1.58716-0.4015im
0.993204-1.84695im    -1.48726-0.244343im      -1.33509+0.742391im
-0.0-0.0im       -0.682093-1.73745im   …    1.08995-0.3994im
0.83187-1.05805im    0.247246+0.45962im       0.700845-0.604206im
-0.0+0.0im       -0.398911-0.745511im       1.23564-0.586102im
⋮                                  ⋱
-0.0-0.0im         1.23564+0.586102im  …  -0.398911+0.745511im
0.83187+1.05805im    0.700845+0.604206im      0.247246-0.45962im
-0.0+0.0im         1.08995+0.3994im       -0.682093+1.73745im
0.993204+1.84695im    -1.33509-0.742392im      -1.48726+0.244343im
-0.569311+0.075523im    1.58716+0.4015im        0.416833+0.219624im
2.86644+0.98837im    -0.55987+1.09585im   …  -0.248534-1.02317im
1.22427+0.125188im    1.15793+2.06182im       -2.10585+2.96785im
-1.61084+4.77705im    -2.28491-1.02422im       0.165628+1.66687im
-0.624841+5.30147im   -0.988565-0.720973im       1.55698+1.17529im
4.73903-1.55932im     6.26791+6.04117im      -0.805839+5.27302im
-3.55624+1.76501im     7.28667+0.598839im  …   -2.36528+0.595752im
0.0844855+11.1962im    0.630981-7.94724im       -15.1921-13.9858im 

Inverse the Fourier transform.

In [80]:
fM = real((plan_ifft(c)*c)*n);

Out[80]:
512×512 Array{Float32,2}:
0.483716  0.447112  0.418023  0.396734  …  0.585784  0.562059  0.524727
0.51758   0.488234  0.465988  0.451702     0.604017  0.582666  0.551245
0.539193  0.521543  0.510711  0.507189     0.601224  0.584037  0.561614
0.54883   0.547347  0.551568  0.560633     0.575548  0.56518   0.555555
0.549028  0.56703   0.587815  0.608695     0.528842  0.528631  0.535878
0.543316  0.581933  0.618134  0.647772  …  0.466492  0.479865  0.507471
0.534777  0.592552  0.640765  0.67498      0.396211  0.425793  0.475666
0.525114  0.598457  0.654222  0.689204     0.326193  0.372947  0.4447
0.514596  0.598928  0.658206  0.691695     0.263264  0.326025  0.416883
0.50279   0.593895  0.654202  0.685873     0.211683  0.287318  0.392734
0.489592  0.58462   0.645369  0.67632   …  0.172922  0.257102  0.371857
0.475928  0.573684  0.635677  0.667361     0.14634   0.234615  0.35398
0.463739  0.564229  0.628635  0.661802     0.130261  0.219043  0.33957
⋮                                       ⋱            ⋮
0.164341  0.157544  0.158637  0.163874  …  0.177551  0.182946  0.176
0.17064   0.160941  0.159687  0.162059     0.204989  0.201743  0.187224
0.17822   0.165681  0.161262  0.160214     0.230465  0.2189    0.198424
0.187692  0.171487  0.162934  0.158407     0.254148  0.236004  0.211164
0.199676  0.178702  0.165474  0.15798      0.277641  0.25518   0.22694
0.215206  0.188769  0.17102   0.161295  …  0.303283  0.2785    0.247101
0.235894  0.20412   0.182518  0.170985     0.333354  0.307573  0.272892
0.263563  0.227351  0.202647  0.189068     0.369378  0.34327   0.305329
0.299379  0.260039  0.232722  0.216331     0.411582  0.385421  0.34475
0.342859  0.301752  0.272093  0.252257     0.458498  0.432445  0.390143
0.391305  0.349777  0.318321  0.295429  …  0.50677   0.481046  0.43861
0.440074  0.399778  0.368038  0.344128     0.551284  0.526253  0.485399

Display the approximation.

In [81]:
figure(figsize = (5,5))
imageplot(clamP(fM))


Exercise 1

Compute a best $M$-term approximation in the Fourier basis of $f$, for $M \in \{N/100, N/20\}$. Compute the approximation using a well chosen hard threshold value $T$.

In [82]:
include("NtSolutions/coding_1_approximation/exo1.jl")

In [10]:
## Insert your code here.


The best $M$-term approximation error is computed using the conservation of energy as

$$\epsilon[M]^2 = \norm{f-f_M}^2 = \sum_{ \abs{\dotp{f}{\psi_m}} \leq T } \abs{\dotp{f}{\psi_m}}^2.$$

If one denotes by $\{ c_R[k] \}_{k=0}^{N-1}$ the set of coefficients magnitudes $\abs{\dotp{f}{\psi_m}}$ ordered by decaying magnitudes, then this error is easily computed as $$\epsilon[M]^2 = \sum_{k=M}^{N-1} c_R[k]^2 = \norm{f}^2 - \sum_{k=0}^{M-1} c_R[k]^2.$$ This means that $\epsilon^2$ is equal to $\norm{f}^2$ minus the discrete primitive of $c_R^2$.

Exercise 2

Compute and display in log scales the ordered coefficients $c_R$. Hint: a discrete primitive can be computed using the numpy function cumsum.

In [83]:
include("NtSolutions/coding_1_approximation/exo2.jl")

In [12]:
## Insert your code here.


Exercise 3

Compute and display in log-scale the non-linear approximation error $\epsilon[M]^2$. Store the values of $\epsilon[M]^2$ in a vector $err\_fft$.

In [84]:
include("NtSolutions/coding_1_approximation/exo3.jl")

In [14]:
## Insert your code here.


## Wavelet Approximation¶

The Wavelet basis of continuous 2-D functions is defined by by scaling and translating three mother atoms $\{\psi^H,\psi^V,\psi^D\}$: $$\psi_{j,n}^k(x) = \frac{1}{2^j}\psi^k\pa{\frac{x-2^j n}{2^j}}$$

Non-linear wavelet approximation is a the heart of the JPEG-2000 compression standard.

The set of inner products $\{ \dotp{f}{\psi_m} \}_m$ is computed in $O(N)$ operations with the 2-D Fast Wavelet Transform algorithm.

Perform a wavelet transform. Here we use a daubechies wavelet transform.

In [85]:
Jmin = 1

h = compute_wavelet_filter("Daubechies",10)

fW = perform_wavortho_transf(f, Jmin, + 1, h);

Out[85]:
512×512 Array{Float32,2}:
130.682        142.046         23.0967       …  -0.0406469     0.0246702
96.9038       143.32          -2.83519         -0.016554      0.0195549
7.4268        -5.12091      -15.7467           0.0234772     0.00472445
4.11817       -1.29105        1.42537          0.0293722     0.0384212
0.0566384     -1.73085        0.428896         0.0393595    -0.00673362
-6.64421       11.2755         3.53238      …   0.0199681     0.0232348
5.74178        9.06315       -0.550003        -0.0177949     0.0120984
6.68185       -5.23018       -1.82195          0.000925797  -0.010417
0.820846      -1.09493       -1.861            0.0400613     0.0276619
1.91916        0.913518       2.53322          0.018259     -0.00182923
1.62441       -0.901944      -3.52287      …   0.00866102    0.00382185
-2.74687        3.32436        3.65457         -0.0347214     0.00603406
2.21268       -3.21147        0.576862         0.0249128     0.0243874
⋮                                          ⋱   ⋮
0.0244523     -0.00824016    -0.00759765   …   0.0161912     0.0022017
0.00646857    -0.000968072    0.00924794       0.00881646   -0.0119001
0.00273045     0.0153158      0.00559765       0.0113431     0.00890303
-0.00496361     0.00441209    -0.00164017      -0.0269199     0.0237037
0.0266946     -0.0142672     -0.0303202       -0.00764149   -0.0236405
0.000783565    0.0103738     -0.0105912    …   0.000633202   0.00325321
0.0131437     -0.0114727      0.00803314      -0.00283124   -0.018369
0.0128851     -0.0333292      0.0088044        0.00507466   -0.0270831
0.015645       0.00695539    -0.000900542      0.000450602   0.00291813
0.0229567      0.0167443      0.0104381       -0.0149523     0.0132251
-0.0172143      0.0177822      0.00987626   …  -0.000377828  -0.00135188
0.0279684      0.00113695     0.00962025       0.00178862    0.0188452 
In [23]:
# using NPZ

Out[23]:
512×512 Array{Float64,2}:
142.692        98.3625        13.3814       …   0.237184    -0.113969
142.945       128.952         -3.23211          0.0465118   -0.0287645
-1.41844       3.26086      -11.5858          -0.151489     0.051181
-7.05074      10.0609         1.12759         -0.138798     0.049865
2.58432      -3.79339        7.71082         -0.287621     0.101513
3.95298      11.1895         4.43363      …  -0.322094     0.103109
-3.73472       2.47893       -3.37804         -0.317403     0.100171
4.85028      -4.14531        5.21644         -0.311572     0.103892
-0.62423      -2.07634        3.26276         -0.327977     0.0993851
1.33123       3.2964        -3.53485         -0.300084     0.0979851
1.46328       1.87355        2.31762      …  -0.292651     0.0915352
-3.27029      -2.82297       -0.00836166      -0.295894     0.0948614
-2.82005       1.85118       -2.18305         -0.309516     0.106956
⋮                                         ⋱   ⋮
-0.0172059    -0.00283857    -0.0074868    …  -0.00768957   0.00561585
0.0189631    -0.0176611     -0.00816675       0.00516067  -0.0203371
0.0131052     0.0251774      0.00513254       0.0270514   -0.014152
-0.00317476    0.0492549      0.00600047      -0.0123553    0.0360208
-0.0178114     0.0179188     -0.0151002        0.0335839   -0.0124882
0.0133063    -0.00707918    -0.00425823   …  -0.0205506    0.0122348
0.0148492     0.0252244      0.00765644       0.0226339   -0.00848701
0.0139278     0.000284363   -0.000902768     -0.0104011    0.00998178
-0.0298214     0.0190962      0.0054781       -0.021531     0.0266814
0.00232894   -0.0207997     -0.0166564       -0.00244105  -0.000962089
-0.12244      -0.200709      -0.317618     …   0.112873    -0.0401913
0.0381981     0.0641532      0.101473        -0.0375089    0.0140014  

Display the coefficients.

In [86]:
figure(figsize = (8,8))

plot_wavelet(fW, Jmin)
title("Wavelet coefficients")

show()