Filters in Image Processing
Image Transforms (II) – Wavelets
David Svoboda
email: svoboda@fi.muni.cz
Centre for Biomedical Image Analysis
Faculty of Informatics, Masaryk University, Brno, CZ
October 11, 2019
David Svoboda (CBIA@FI) Filters in Image Processing autumn 2019 1 / 50
Outline
1 Motivation
2 New basis
3 1D Discrete Wavelet Transform
4 Subband coding
Signal Analysis
From Filter Banks to Wavelets
5 1D Fast Discrete Wavelet Transform
6 2D Discrete Wavelet Transform
7 Wavelet Packets
David Svoboda (CBIA@FI) Filters in Image Processing autumn 2019 2 / 50
1 Motivation
2 New basis
3 1D Discrete Wavelet Transform
4 Subband coding
Signal Analysis
From Filter Banks to Wavelets
5 1D Fast Discrete Wavelet Transform
6 2D Discrete Wavelet Transform
7 Wavelet Packets
David Svoboda (CBIA@FI) Filters in Image Processing autumn 2019 3 / 50
Motivation
Fourier transform and its derivatives have nice properties
reversible
linear operator
easy to construct (data independent)
decorrelation property
Fourier transform has some drawbacks
localization is very expensive
complexity cannot be lower than O(n log n)
David Svoboda (CBIA@FI) Filters in Image Processing autumn 2019 4 / 50
1 Motivation
2 New basis
3 1D Discrete Wavelet Transform
4 Subband coding
Signal Analysis
From Filter Banks to Wavelets
5 1D Fast Discrete Wavelet Transform
6 2D Discrete Wavelet Transform
7 Wavelet Packets
David Svoboda (CBIA@FI) Filters in Image Processing autumn 2019 5 / 50
New basis
Fourier basis originates from sin and cos functions.
Let us introduce so called scaling function ϕ and wavelet function φ.
ϕ and ψ should be strictly localized.
A new basis should originate from ϕ and ψ.
Scaling function ϕ
ϕ(x) =
1 iff 0 ≤ x < 1
0 otherwise
0 1 2 3 4
0
1
2
Wavelet function ψ
ψ(x) =



1 iff 0 ≤ x < 0.5
−1 iff 0.5 ≤ x < 1
0 elsewhere
0 1 2 3 4
1
2
0
−1
David Svoboda (CBIA@FI) Filters in Image Processing autumn 2019 6 / 50
New basis
Scaling function
Let ϕ(x) be a scaling function. We will modify it as follows:
ϕj,k(x) = 2j/2
ϕ 2j x
M
− k
where j, k ∈ Z then
j . . . ϕj,k(x)’s width
controls broadness and height of the function along x and y axis,
respectively.
k . . . shift of ϕj,k(x) along x-axis
M . . . length of the processed signal
Notice: The function is orthogonal to its integer shifts
David Svoboda (CBIA@FI) Filters in Image Processing autumn 2019 7 / 50
New basis
Scaling function (M = 1) – examples
−1
1
1
−1
1
1
−1
1
1
−1
1
1
−1
1
1
−1
1
1
−1
1
1
ϕ0,0(x) = ϕ(x)
ϕ1,1(x) =
√
2ϕ(2x − 1)
ϕ2,0(x) = 2ϕ(4x)
ϕ1,0(x) =
√
2ϕ(2x)
ϕ2,1(x) = 2ϕ(4x − 1) ϕ2,2(x) = 2ϕ(4x − 2) ϕ2,3(x) = 2ϕ(4x − 3)
David Svoboda (CBIA@FI) Filters in Image Processing autumn 2019 8 / 50
New basis
Scaling function (M = 4) – examples
−1
1
4
−1
4
1
−1
4
1
−1
4
1
−1
4
1
−1
4
1
−1
4
1
ϕ0,0(x) = ϕ(x
4)
ϕ1,1(x) =
√
2ϕ(2x
4 − 1)ϕ1,0(x) =
√
2ϕ(2x
4)
ϕ2,2(x) = 2ϕ(4x
4 − 2) ϕ2,3(x) = 2ϕ(4x
4 − 3)ϕ2,0(x) = 2ϕ(4x
4) ϕ2,1(x) = 2ϕ(4x
4 − 1)
David Svoboda (CBIA@FI) Filters in Image Processing autumn 2019 9 / 50
New basis
Wavelet function
Let ψ(x) be a wavelet function. We will modify it as follows:
ψj,k(x) = 2j/2
ψ 2j x
M
− k
where j, k ∈ Z then
j . . . ϕj,k(x)’s width
controls broadness and height of the function along x and y axis,
respectively.
k . . . shift of ϕj,k(x) along x-axis
M . . . length of the processed signal
Notice: The function is orthogonal to its integer shifts
David Svoboda (CBIA@FI) Filters in Image Processing autumn 2019 10 / 50
New basis
Wavelet function (M = 1) – examples
−1
1
1
−1
1
1
−1
1
1
−1
1
1
−1
1
1
−1
1
1
−1
1
1
ψ0,0(x) = ψ(x)
ψ1,1(x) =
√
2ψ(2x − 1)
ψ2,0(x) = 2ψ(4x)
ψ1,0(x) =
√
2ψ(2x)
ψ2,1(x) = 2ψ(4x − 1) ψ2,2(x) = 2ψ(4x − 2) ψ2,3(x) = 2ψ(4x − 3)
David Svoboda (CBIA@FI) Filters in Image Processing autumn 2019 11 / 50
New basis
Wavelet function (M = 4) – examples
−1
4
1
−1
4
1
−1
4
1
−1
4
1
−1
4
1
−1
4
1
−1
4
1
ψ0,0(x) = ψ(x
4)
ψ1,1(x) =
√
2ψ(2x
4 − 1)
ψ2,0(x) = 2ψ(4x
4)
ψ1,0(x) =
√
2ψ(2x
4)
ψ2,1(x) = 2ψ(4x
4 − 1) ψ2,2(x) = 2ψ(4x
4 − 2) ψ2,3(x) = 2ψ(4x
4 − 3)
David Svoboda (CBIA@FI) Filters in Image Processing autumn 2019 12 / 50
New basis
Support of basis functions
Fourier basis
The domain of sin and cos is 0; 1 .
The domain is periodic.
When applying the Fourier basis function ϕm(k) to a transformed
function f of length N, this basis function ϕm(k) (its period) is
stretched to the length N.
New basis
The domain of ϕ(x) and ψ(x) is 0; 1 .
Each function has limited compact support.
When applying the scaling or wavelet function to a transformed
function f of length M, both scaling and wavelet function are
appropriately stretched to the required length.
David Svoboda (CBIA@FI) Filters in Image Processing autumn 2019 13 / 50
New basis
Design of transform matrix
Let us recall that the basis forms the rows of transform matrix!
Sample matrix for signal of length M = 4
HT4 =
1
2




1 1 1 1
1 1 −1 −1√
2 −
√
2 0 0
0 0
√
2 −
√
2




= ϕ0,0
= ψ0,0
= ψ1,0
= ψ1,1
Notice the structure of individual rows.
Explain the value 1
2.
David Svoboda (CBIA@FI) Filters in Image Processing autumn 2019 14 / 50
New basis
Design of transform matrix
Matrices for signal of length M = 8
HT8 =
1
√
8











2 2 0 0 0 0 0 0
0 0 2 2 0 0 0 0
0 0 0 0 2 2 0 0
0 0 0 0 0 0 2 2
2 −2 0 0 0 0 0 0
0 0 2 −2 0 0 0 0
0 0 0 0 2 −2 0 0
0 0 0 0 0 0 2 −2











= ϕ2,0
= ϕ2,1
= ϕ2,2
= ϕ2,3
= ψ2,0
= ψ2,1
= ψ2,2
= ψ2,3
0 1 2 3 4 5 6 7 8
-4
-2
0
2
4
6
8
10
⇒
0 1 2 3 4 5 6 7 8
-4
-2
0
2
4
6
8
10
f = [3 5 3 7 0 -1 2 4] HT8 * f’
David Svoboda (CBIA@FI) Filters in Image Processing autumn 2019 15 / 50
New basis
Design of transform matrix
Matrices for signal of length M = 8
HT8 =
1
√
8












√
2
√
2
√
2
√
2 0 0 0 0
0 0 0 0
√
2
√
2
√
2
√
2√
2
√
2 −
√
2 −
√
2 0 0 0 0
0 0 0 0
√
2
√
2 −
√
2 −
√
2
2 −2 0 0 0 0 0 0
0 0 2 −2 0 0 0 0
0 0 0 0 2 −2 0 0
0 0 0 0 0 0 2 −2












= ϕ1,0
= ϕ1,1
= ψ1,0
= ψ1,1
= ψ2,0
= ψ2,1
= ψ2,2
= ψ2,3
0 1 2 3 4 5 6 7 8
-4
-2
0
2
4
6
8
10
⇒
0 1 2 3 4 5 6 7 8
-4
-2
0
2
4
6
8
10
f = [3 5 3 7 0 -1 2 4] HT8 * f’
David Svoboda (CBIA@FI) Filters in Image Processing autumn 2019 15 / 50
New basis
Design of transform matrix
Matrices for signal of length M = 8
HT8 =
1
√
8












1 1 1 1 1 1 1 1
1 1 1 1 −1 −1 −1 −1√
2
√
2 −
√
2 −
√
2 0 0 0 0
0 0 0 0
√
2
√
2 −
√
2 −
√
2
2 −2 0 0 0 0 0 0
0 0 2 −2 0 0 0 0
0 0 0 0 2 −2 0 0
0 0 0 0 0 0 2 −2












= ϕ0,0
= ψ0,0
= ψ1,0
= ψ1,1
= ψ2,0
= ψ2,1
= ψ2,2
= ψ2,3
0 1 2 3 4 5 6 7 8
-4
-2
0
2
4
6
8
10
⇒
0 1 2 3 4 5 6 7 8
-4
-2
0
2
4
6
8
10
f = [3 5 3 7 0 -1 2 4] HT8 * f’
David Svoboda (CBIA@FI) Filters in Image Processing autumn 2019 15 / 50
1 Motivation
2 New basis
3 1D Discrete Wavelet Transform
4 Subband coding
Signal Analysis
From Filter Banks to Wavelets
5 1D Fast Discrete Wavelet Transform
6 2D Discrete Wavelet Transform
7 Wavelet Packets
David Svoboda (CBIA@FI) Filters in Image Processing autumn 2019 16 / 50
1D Discrete Wavelet Transform (DWT)
Forward
Aj0 (k) =
1
√
M
M−1
m=0
f (m)ϕj0,k(m)
Dj (k) =
1
√
M
M−1
m=0
f (m)ψj,k(m)
Inverse
f (m) =
1
√
M
2j0 −1
k=0
Aj0 (k)ϕj0,k(m) +
1
√
M
J−1
j=j0
2j −1
k=0
Dj (k)ψj,k(m)
ϕ, ψ . . . orthogonal scaling and wavelet
function, respectively
Aj0
(k) . . . scaling coefficients
(approximations)
Dj (k) . . . wavelet coefficients (details)
M = 2J . . . number of samples in
function f
j ∈ {j0, ..., J − 1} . . . level of detail,
where j0 ≥ 0
k ∈ {0, 1, . . . , 2j − 1}
David Svoboda (CBIA@FI) Filters in Image Processing autumn 2019 17 / 50
1D Discrete Wavelet Transform (DWT)
An example
Input:
f = [1, 4, −3, 0]
|f | = M = 4 ⇒ J = 2
decomposition level: j0 = 1
)
PPPPPPPPq
A1 D1
f = A2
A1(0) =
1
2
3
m=0
f (m)ϕ1,0(m) =
1
2
[1 ·
√
2 + 4 ·
√
2 + (−3) · 0 + 0 · 0] =
5
√
2
A1(1) =
1
2
3
m=0
f (m)ϕ1,1(m) =
1
2
[1 · 0 + 4 · 0 + (−3) ·
√
2 + 0 ·
√
2] =
−3
√
2
D1(0) =
1
2
3
m=0
f (m)ψ1,0(m) =
1
2
[1 ·
√
2 + 4 · (−
√
2) − 3 · 0 + 0 · 0] = −
3
√
2
D1(1) =
1
2
3
m=0
f (m)ψ1,1(m) =
1
2
[1 · 0 + 4 · 0 − 3 ·
√
2 + 0 · (−
√
2)] = −
3
√
2
David Svoboda (CBIA@FI) Filters in Image Processing autumn 2019 18 / 50
1D Discrete Wavelet Transform (DWT)
An example
Input:
f = [1, 4, −3, 0]
|f | = M = 4 ⇒ J = 2
decomposition level: j0 = 0
)
PPPPPPq
PPPPPPPPq
9
D1
f = A2
D0A0
A0(0) =
1
2
3
m=0
f (m)ϕ0,0(m) =
1
2
[1 · 1 + 4 · 1 − 3 · 1 + 0 · 1] = 1
D0(0) =
1
2
3
m=0
f (m)ψ0,0(m) =
1
2
[1 · 1 + 4 · 1 − 3 · (−1) + 0 · (−1)] = 4
D1(0) =
1
2
3
m=0
f (m)ψ1,0(m) =
1
2
[1 ·
√
2 + 4 · (−
√
2) − 3 · 0 + 0 · 0] = −
3
2
√
2
D1(1) =
1
2
3
m=0
f (m)ψ1,1(m) =
1
2
[1 · 0 + 4 · 0 − 3 ·
√
2 + 0 · (−
√
2)] = −
3
2
√
2
David Svoboda (CBIA@FI) Filters in Image Processing autumn 2019 19 / 50
1D Discrete Wavelet Transform (DWT)
An example
DWT(f ) = [1, 4, −3
2
√
2, −3
2
√
2], i.e.
f (m) = IDWT([1, 4, −
3
2
√
2, −
3
2
√
2])
=
1
2
A0(0) · ϕ0,0(m) +
1
2
(D0(0) · ψ0,0(m) + D1(0) · ψ1,0(m) + D1(1) · ψ1,1(m))
=
1
2
· 1 · ϕ0,0(m) +
1
2
4 · ψ0,0(m) −
3
2
√
2 · ψ1,0(m) −
3
2
√
2 · ψ1,1(m)
Utilization of the same basis functions in forward and inverse transforms is
conditioned to orthonormality of selected functions.
David Svoboda (CBIA@FI) Filters in Image Processing autumn 2019 20 / 50
1D Discrete Wavelet Transform (DWT)
An example
After submitting signal
f = f(0) f(1) f(2) f(3)
to 1D-DWT, we obtain separately approximations and details of the signal:
for j0 = 2:
no decomposition
for j0 = 1:
DWT(f) = A1(0) A1(1) D1(0) D1(1)
for j0 = 0:
DWT(f) = A0(0) D0(0) D1(0) D1(1)
Notice: The output signal is always of the same length as the input signal.
David Svoboda (CBIA@FI) Filters in Image Processing autumn 2019 21 / 50
1D Discrete Wavelet Transform (DWT)
Variability & Issues
The common families of scaling (father) and wavelet (mother) functions
Haar (already introduced)
Daubechies: db1, db2, db3, db4, . . .
Meyer
Coiflets: coif1, coif2, coif3, . . .
Symlets: sym2, sym3, sym4, . . .
Biorthogonal: bior1, bior2, bior3, . . .
Complexity of 1D-DWT
matrix multiplication – O(n2)
the whole transform matrix typically built only for Haar wavelets
other wavelets computed iteratively
(one matrix per one level of decomposition) ⇒ iterations × O(n2)
can we speed it up?
David Svoboda (CBIA@FI) Filters in Image Processing autumn 2019 22 / 50
1 Motivation
2 New basis
3 1D Discrete Wavelet Transform
4 Subband coding
Signal Analysis
From Filter Banks to Wavelets
5 1D Fast Discrete Wavelet Transform
6 2D Discrete Wavelet Transform
7 Wavelet Packets
David Svoboda (CBIA@FI) Filters in Image Processing autumn 2019 23 / 50
Subband Coding
Signal Analysis
Any signal f can be decomposed into two parts:
approximation (A) . . . obtained by low-pass filtering of the original
signal
detail (D) . . . obtained by high-pass filtering of the original signal
f
DA
Filters
low−pass high−pass
David Svoboda (CBIA@FI) Filters in Image Processing autumn 2019 24 / 50
Subband Coding
Signal Analysis
The filtered signal must be downsampled (↓2×) to avoid data redundancy.
D − high frequencies
A − low frequencies
256 data points
128 data points
128 data points
2x
2x
f
David Svoboda (CBIA@FI) Filters in Image Processing autumn 2019 25 / 50
Subband Coding
Signal Analysis and Synthesis
The decomposed signal may be reconstructed:
detail (D) is upsampled (↑2×) and then high-pass filtered
approximation (A) is upsampled (↑2×) and then low-pass filtered
results are added → f
2x
2x
2x
A
D
2x
analysis synthesis
H’H
L’L
f f’
D
A
signaltransmission
Notice: We would like to have f = f
David Svoboda (CBIA@FI) Filters in Image Processing autumn 2019 26 / 50
Subband Coding
Signal Analysis and Synthesis
Filter banks
H . . . high-pass analysis filter (FIR)
L . . . low-pass analysis filter (FIR)
H . . . high-pass synthesis filter (FIR)
L . . . low-pass synthesis filter (FIR)
2x
2x
2x
A
D
2x
analysis synthesis
H’H
L’L
f f’
D
A
signaltransmission
David Svoboda (CBIA@FI) Filters in Image Processing autumn 2019 27 / 50
Subband Coding
Signal Analysis and Synthesis
Filter banks
If f = f then the filters L, L , H, H are called perfect reconstruction filters
and they must fulfill one of the following conditions:
H (n) = (−1)n
L(n)
L (n) = (−1)n+1
H(n)
H (n) = (−1)n+1
L(n)
L (n) = (−1)n
H(n)
2x
2x
2x
A
D
2x
analysis synthesis
H’H
L’L
f f’
D
A
signaltransmission
David Svoboda (CBIA@FI) Filters in Image Processing autumn 2019 28 / 50
Subband Coding
Signal Analysis and Synthesis
Filter banks
H and L are mutually cross-modulated
H and L are mutually cross-modulated
H, H , L, L are called quadrature mirror filters (QMF)
2x
2x
2x
A
D
2x
analysis synthesis
H’H
L’L
f f’
D
A
signaltransmission
David Svoboda (CBIA@FI) Filters in Image Processing autumn 2019 29 / 50
Subband Coding
Signal Analysis and Synthesis
Filter banks
Biorthogonal filters
We need to define two filters
H and L. The remaining H
and L are derived by
cross-modulation.
2x
2x
2x
A
D
2x
analysis synthesis
H’H
L’L
f f’
D
A
signaltransmission
Orthogonal filters
We define only one filter H . The
remaining filters fulfill:
L (n) = (−1)n
H (length − 1 − n)
H(n) = H (length − 1 − n)
L(n) = L (length − 1 − n)
where
length = size(H ) &
is even(length) = true
Notice: We will focus namely on the orthogonal filters.
David Svoboda (CBIA@FI) Filters in Image Processing autumn 2019 30 / 50
Subband Coding
Recursive Signal Analysis
Once the input signal is decomposed into two parts (A and D), its
approximation (A) can be further decomposed. In the reverse order, the
same is valid for reconstruction.
2x
2x
2x
2x
f
2x
2x
f’
2x
2x
H
L
H
L
......
256
64
128
64
H’
L’
H’
L’
analysis synthesis
signaltransmission
Notice: Let us assume we have already employed (bi)orthogonal filters.
David Svoboda (CBIA@FI) Filters in Image Processing autumn 2019 31 / 50
Subband Coding
The most common orthogonal filters
. . . and their scaling and wavelet functions
Notice: Useful web-pages: http://wavelets.pybytes.com/
David Svoboda (CBIA@FI) Filters in Image Processing autumn 2019 32 / 50
Subband Coding
The most common orthogonal filters
. . . and their scaling and wavelet functions
Notice: Useful web-pages: http://wavelets.pybytes.com/
David Svoboda (CBIA@FI) Filters in Image Processing autumn 2019 32 / 50
Subband Coding
The most common orthogonal filters
. . . and their scaling and wavelet functions
Notice: Useful web-pages: http://wavelets.pybytes.com/
David Svoboda (CBIA@FI) Filters in Image Processing autumn 2019 32 / 50
Subband Coding
The most common orthogonal filters
. . . and their scaling and wavelet functions
Notice: Useful web-pages: http://wavelets.pybytes.com/
David Svoboda (CBIA@FI) Filters in Image Processing autumn 2019 32 / 50
From Filter Banks to Wavelets
Cascade algorithm for ϕ function (numerical solution)
Algorithm
1: L ← fetch low-pass synthesis filter from the selected filter bank
2: hϕ = fliplr(L )
3: ϕ ← Dirac delta impulse
4: while (ϕ is converging) do
5: ϕ ← conv(ϕ, hϕ)
6: ϕ ← upsample(ϕ, 2×)
7: end while
8: OUTPUT ← ϕ
0 0.5 1 1.5 2 2.5 3
−0.2
0
0.2
0.4
0.6
0.8
1
1.2
⇒
0 0.5 1 1.5 2 2.5 3
−0.4
−0.2
0
0.2
0.4
0.6
0.8
1
1.2
1.4
⇒
0 0.5 1 1.5 2 2.5 3
−0.4
−0.2
0
0.2
0.4
0.6
0.8
1
1.2
1.4
⇒
0 0.5 1 1.5 2 2.5 3
−0.4
−0.2
0
0.2
0.4
0.6
0.8
1
1.2
1.4
David Svoboda (CBIA@FI) Filters in Image Processing autumn 2019 33 / 50
From Filter Banks to Wavelets
Cascade algorithm for ψ function (numerical solution)
Algorithm
1: ϕ ← call Cascade algorithm to get ϕ function
2: H ← fetch high-pass synthesis filter from the selected filter bank
3: hψ = fliplr(H )
4: ψ ← conv(ϕ, hψ)
5: ψ ← upsample(ψ, 2×)
6: OUTPUT ← ψ
0 0.5 1 1.5 2 2.5 3
−0.2
0
0.2
0.4
0.6
0.8
1
1.2
⇒
0 0.5 1 1.5 2 2.5 3
−0.4
−0.2
0
0.2
0.4
0.6
0.8
1
1.2
1.4
⇒
0 0.5 1 1.5 2 2.5 3
−0.4
−0.2
0
0.2
0.4
0.6
0.8
1
1.2
1.4
⇒
0 0.5 1 1.5 2 2.5 3
−1.5
−1
−0.5
0
0.5
1
1.5
2
David Svoboda (CBIA@FI) Filters in Image Processing autumn 2019 34 / 50
1 Motivation
2 New basis
3 1D Discrete Wavelet Transform
4 Subband coding
Signal Analysis
From Filter Banks to Wavelets
5 1D Fast Discrete Wavelet Transform
6 2D Discrete Wavelet Transform
7 Wavelet Packets
David Svoboda (CBIA@FI) Filters in Image Processing autumn 2019 35 / 50
1D Fast Discrete Wavelet Transform
Definition
Dj (k) =
|Aj+1|−1
r=0
H (2k + 1 − r)Aj+1(r)
Aj (k) =
|Aj+1|−1
r=0
L (2k + 1 − r)Aj+1(r)
AJ(k) = f (k)
Each step in FWT corresponds to convolution with high-pass and low-pass
analysis filter followed by down-sampling (↓2×).
1D-DWT ≡ Subband coding
David Svoboda (CBIA@FI) Filters in Image Processing autumn 2019 36 / 50
1D Fast Discrete Wavelet Transform
Basic scheme
Let |f | = M = 8 = 23 = 2J and j0 = 0
f (0) f (1) f (2) f (3) f (4) f (5) f (6) f (7)
A3(0) A3(1) A3(2) A3(3) A3(4) A3(5) A3(6) A3(7)
⇓
A2(0) A2(1) A2(2) A2(3) D2(0) D2(1) D2(2) D2(3)
⇓
A1(0) A1(1) D1(0) D1(1)
⇓
A0(0) D0(0)
David Svoboda (CBIA@FI) Filters in Image Processing autumn 2019 37 / 50
Fast Wavelet Transform
An example
Given f (k) = [1, 4, −3, 0] and Haar scaling and wavelet coefficients
L (k) =
1/
√
2 k = 0, 1
0 otherwise
/and/ H (k) =



−1/
√
2 k = 0
1/
√
2 k = 1
0 otherwise
we can evaluate the following:
level 2: A2(k) = f (k) = [1, 4, −3, 0]
level 1: A1(k) =
3
r=0
L (2k + 1 − r)A2(r) = [5/
√
2, −3/
√
2]
D1(k) =
3
r=0
H (2k + 1 − r)A2(r) = [−3/
√
2, −3/
√
2]
level 0: A0(k) =
1
r=0
L (2k + 1 − r)A1(r) = [1]
D0(k) =
1
r=0
H (2k + 1 − r)A1(r) = [4]
David Svoboda (CBIA@FI) Filters in Image Processing autumn 2019 38 / 50
Comparison of FWT and FFT
Fast Fourier Transform
complexity: O(n log n)
existence: at any time
time versus frequency
domain
Fast Wavelet Transform
complexity O(cn)
c . . . support of L filter (typically small)
existence: depends upon the availability
of scaling function and the
orthogonality of the scaling function
time & frequency changes
simultaneously
David Svoboda (CBIA@FI) Filters in Image Processing autumn 2019 39 / 50
1 Motivation
2 New basis
3 1D Discrete Wavelet Transform
4 Subband coding
Signal Analysis
From Filter Banks to Wavelets
5 1D Fast Discrete Wavelet Transform
6 2D Discrete Wavelet Transform
7 Wavelet Packets
David Svoboda (CBIA@FI) Filters in Image Processing autumn 2019 40 / 50
2D Discrete Wavelet Transform
An extension of scaling function and wavelets to 2D in straightforward:
ϕ(x) → ϕ(x, y)
ψ(x) → ψH
(x, y), ψV
(x, y), ψD
(x, y)
where all the 2D functions are separable in the following manner:
ϕ(x, y) = ϕ(x)ϕ(y)
ψH
(x, y) = ψ(x)ϕ(y)
ψV
(x, y) = ϕ(x)ψ(y)
ψD
(x, y) = ψ(x)ψ(y)
David Svoboda (CBIA@FI) Filters in Image Processing autumn 2019 41 / 50
2D Discrete Wavelet Transform
What is the meaning of new wavelets?
ψH(x, y) . . . intensity variations for image columns
ψV (x, y) . . . intensity variations along rows
ψD(x, y) . . . intensity variations along diagonals
Corollary:
ϕj,m,n(x, y) = 2j/2
ϕ 2j x
M
− m, 2j y
N
− n
ψi
j,m,n(x, y) = 2j/2
ψi
2j x
M
− m, 2j y
N
− n , i = {H, V , D}
David Svoboda (CBIA@FI) Filters in Image Processing autumn 2019 42 / 50
2D Discrete Wavelet Transform
Definition
Forward
Aj0 (m, n) =
1
√
MN
M−1
k=0
N−1
l=0
f (k, l)ϕj0,m,n(k, l)
Di
j (m, n) =
1
√
MN
M−1
k=0
N−1
l=0
f (k, l)ψi
j,m,n(k, l)
Inverse
f (k, l) =
1
√
MN m n
Aj0 (m, n)ϕj0,m,n(k, l)
+
1
√
MN i=H,V ,D
J−1
j=j0 m n
Di
j (m, n)ψi
j,m,n(k, l)
where i = {H, V , D}
David Svoboda (CBIA@FI) Filters in Image Processing autumn 2019 43 / 50
2D Discrete Wavelet Transform
Practical implementation
1-level decomposition as a 2-step process
AJ(m, n)
DH
J−1(m, n)AJ−1(m, n)
DD
J−1(m, n)DV
J−1(m, n)
A J−1(m, n)
D V
J−1(m, n)
n-th level decomposition as a iterative process
AJ(m, n)
DH
J−1(m, n)
DD
J−1(m, n)DV
J−1(m, n)
AJ−1(m, n)
David Svoboda (CBIA@FI) Filters in Image Processing autumn 2019 44 / 50
2D Discrete Wavelet Transform
An example – DWT using Haar wavelets
Level of detail: j0 = J
David Svoboda (CBIA@FI) Filters in Image Processing autumn 2019 45 / 50
2D Discrete Wavelet Transform
An example – DWT using Haar wavelets
Level of detail: j0 = J − 1
David Svoboda (CBIA@FI) Filters in Image Processing autumn 2019 45 / 50
2D Discrete Wavelet Transform
An example – DWT using Haar wavelets
Level of detail: j0 = J − 2
David Svoboda (CBIA@FI) Filters in Image Processing autumn 2019 45 / 50
2D Discrete Wavelet Transform
An example – DWT using Haar wavelets
Level of detail: j0 = J − 3
David Svoboda (CBIA@FI) Filters in Image Processing autumn 2019 45 / 50
2D Discrete Wavelet Transform
An example – DWT using Haar wavelets
Level of detail: j0 = J − 4
David Svoboda (CBIA@FI) Filters in Image Processing autumn 2019 45 / 50
2D Discrete Wavelet Transform
An example
DWT → modification → IDWT
David Svoboda (CBIA@FI) Filters in Image Processing autumn 2019 46 / 50
1 Motivation
2 New basis
3 1D Discrete Wavelet Transform
4 Subband coding
Signal Analysis
From Filter Banks to Wavelets
5 1D Fast Discrete Wavelet Transform
6 2D Discrete Wavelet Transform
7 Wavelet Packets
David Svoboda (CBIA@FI) Filters in Image Processing autumn 2019 47 / 50
Wavelet Packets
Problem to solve:
Traditional wavelet transform decomposes the (image) data always in the
same manner.
Solution:
Decompose those parts of the data which need it.
An example:
The lowest entropy lead to better compression. Let us split those parts of
the image (not only Aj (m, n)) which need it → which division causes
entropy reduction.
David Svoboda (CBIA@FI) Filters in Image Processing autumn 2019 48 / 50
Bibliography
Burt P. J., Adelson E. H. The Laplacian Pyramid as a Compact Image
Code, IEEE Trans. on Communications, pp. 532–540, April 1983
Gonzalez, R. C., Woods, R. E. Digital image processing / 2nd ed.,
Upper Saddle River: Prentice Hall, 2002, pages 793, ISBN
0201180758
Klette R., Zamperoni P. Handbook of Image Processing Operators,
Wiley, 1996, ISBN-0471956422
Strang G., Nguyen T. Wavelets and Filter Banks,
Wellesley-Cambridge Press, 1997, ISBN 0-9614088-7-1
David Svoboda (CBIA@FI) Filters in Image Processing autumn 2019 49 / 50
You should know the answers . . .
Explain the difference between Fourier basis functions and scaling and
wavelet functions.
Given a signal of fixed length and given a particular scaling a wavelet
function we can perform discrete wavelet transform. The result is
however not unique. Which parameter controls the behaviour of
DWT? Demonstrate on some sample data.
Explain the meaning of A and D coefficients.
Derive the complexity for DWT and separately for FWT.
What would happen if the quadrature mirror filters are not perfect
reconstruction filters.
Describe the Cascade algorithm.
Design an algorithm for computing 2D-FWT.
David Svoboda (CBIA@FI) Filters in Image Processing autumn 2019 50 / 50