Integral and Discrete Transforms in Image Processing
Subband Coding & Fast Wavelet Transform
David Svoboda
email: svoboda@fi.muni.cz
Centre for Biomedical Image Analysis
Faculty of Informatics, Masaryk University, Brno, CZ
October 30, 2022
David Svoboda (CBIA@FI) Image Transforms autumn 2022 1 / 30
Outline
1 Short Revision
2 Troubles when computing DWT
3 Subband coding
Signal Analysis
From Filter Banks to Wavelets
4 1D Fast Discrete Wavelet Transform
David Svoboda (CBIA@FI) Image Transforms autumn 2022 2 / 30
1D Discrete Wavelet Transform (DWT)
Revision
Haar basis generated from φ & ψ:
Scaling function φ
φ(x) =
1 iff 0 ≤ x < 1
0 otherwise
0 1 2 3 4
0
1
2
φj,k(x) = 2j/2
φ 2j x
M
− k
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
ψj,k(x) = 2j/2
ψ 2j x
M
− k
• j . . . scale • k . . . shift • M . . . signal length
David Svoboda (CBIA@FI) Image Transforms autumn 2022 4 / 30
1D Discrete Wavelet Transform (DWT)
Revision
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) Image Transforms autumn 2022 5 / 30
1D Discrete Wavelet Transform (DWT)
Revision
Haar =
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] Haar * fT
David Svoboda (CBIA@FI) Image Transforms autumn 2022 6 / 30
1D Discrete Wavelet Transform (DWT)
Revision
Haar =
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] Haar * fT
David Svoboda (CBIA@FI) Image Transforms autumn 2022 6 / 30
1D Discrete Wavelet Transform (DWT)
Revision
Haar =
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] Haar * fT
David Svoboda (CBIA@FI) Image Transforms autumn 2022 6 / 30
1D Discrete Wavelet Transform (DWT)
Issues
Obstacles when computing DWT:
issue: M (signal length) is not power of 2
possible solution: transform matrix need not be square ⇒ technical
solution but rather complicated
issue: transform matrix is rather sparse ⇒ a lot useless multiplications
possible solution: skipping the zero positions
David Svoboda (CBIA@FI) Image Transforms autumn 2022 8 / 30
1D Discrete Wavelet Transform (DWT)
Issue – Signal Length
Haar transform matrix for the signal of length M = 7 and j ∈ {2}
Haar =
1
√
8











2 2 0 0 0 0 0
0 0 2 2 0 0 0
0 0 0 0 2 2 0
2 0 0 0 0 0 2
2 −2 0 0 0 0 0
0 0 2 −2 0 0 0
0 0 0 0 2 −2 0
−2 0 0 0 0 0 2











⇒
f = [3 5 3 7 0 -1 2] Haar * fT
David Svoboda (CBIA@FI) Image Transforms autumn 2022 9 / 30
1D Discrete Wavelet Transform (DWT)
Issue – Matrix Sparsity
Haar matrix for M = 16 and j0 = 0
David Svoboda (CBIA@FI) Image Transforms autumn 2022 10 / 30
1D Discrete Wavelet Transform (DWT)
Issue – Matrix Sparsity
Questions to be answered
What is the relationship of convolution and matrix multiplication?
Can we subtitute each of them with one another?
Can we express subsampling by a matrix multiplication?
David Svoboda (CBIA@FI) Image Transforms autumn 2022 11 / 30
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) Image Transforms autumn 2022 13 / 30
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) Image Transforms autumn 2022 14 / 30
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) Image Transforms autumn 2022 15 / 30
Subband Coding
Signal Analysis and Synthesis
Filter banks
H . . . high-pass analysis filter
L . . . low-pass analysis filter
H′ . . . high-pass synthesis filter
L′ . . . low-pass synthesis filter
2x
2x
2x
A
D
2x
analysis synthesis
H’H
L’L
f f’
D
A
signaltransmission
David Svoboda (CBIA@FI) Image Transforms autumn 2022 16 / 30
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) Image Transforms autumn 2022 17 / 30
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
Asignaltransmission
David Svoboda (CBIA@FI) Image Transforms autumn 2022 18 / 30
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) Image Transforms autumn 2022 19 / 30
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) Image Transforms autumn 2022 20 / 30
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) Image Transforms autumn 2022 21 / 30
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) Image Transforms autumn 2022 22 / 30
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) Image Transforms autumn 2022 23 / 30
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) Image Transforms autumn 2022 25 / 30
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) Image Transforms autumn 2022 26 / 30
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) Image Transforms autumn 2022 27 / 30
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) Image Transforms autumn 2022 28 / 30
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) Image Transforms autumn 2022 29 / 30
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) Image Transforms autumn 2022 30 / 30