Filters in Image Processing
2nd Generation of Wavelets – Lifting Scheme
David Svoboda
email: svoboda@fi.muni.cz
Centre for Biomedical Image Analysis
Faculty of Informatics, Masaryk University, Brno, CZ
October 21, 2019
David Svoboda (CBIA@FI) Filters in Image Processing autumn 2019 1 / 40
Outline
1 Motivation
2 Introduction to Z-transform
3 Analysis of FWT
4 Lifting scheme (LS)
5 Integer Wavelet transform
6 Applications
David Svoboda (CBIA@FI) Filters in Image Processing autumn 2019 2 / 40
1 Motivation
2 Introduction to Z-transform
3 Analysis of FWT
4 Lifting scheme (LS)
5 Integer Wavelet transform
6 Applications
David Svoboda (CBIA@FI) Filters in Image Processing autumn 2019 3 / 40
Motivation
Practical use of DWT or FWT?
Complexity of:
DWT: O(n2
)
FWT: O(cn)
Can we do it faster?
DWT/FWT are computed in floating point arithmetic.
Can we restrict ourselves to integer number?
DWT/FWT are computed out-of-place.
Can we reduce the memory usage?
Some basis functions in DWT are not orthogonal → dual basis must
be found.
Can we avoid this issue?
David Svoboda (CBIA@FI) Filters in Image Processing autumn 2019 4 / 40
1 Motivation
2 Introduction to Z-transform
3 Analysis of FWT
4 Lifting scheme (LS)
5 Integer Wavelet transform
6 Applications
David Svoboda (CBIA@FI) Filters in Image Processing autumn 2019 5 / 40
Z-Transform
Definition
(bilateral) Z-Transform
F(z) =
∞
n=−∞
f (n)z−n
(unilateral) Z-Transform
F(z) =
∞
n=0
f (n)z−n
Notice: Here z ∈ C, i.e. for z = eiω we get a special case of Z-transform,
which is DFT.
David Svoboda (CBIA@FI) Filters in Image Processing autumn 2019 6 / 40
Z-Transform
Definition
Forward transform
converts discrete series into continuous signal (Z-plane)
F(z) =
∞
n=0
f (n)z−n
Inverse transform
maps continuous signal into discrete series
f (n) = Z−1
(F(z)) =
1
2πi C
F(z)zn−1
dz
where C is a counterclockwise closed path encircling (with radius 1)
the origin
David Svoboda (CBIA@FI) Filters in Image Processing autumn 2019 7 / 40
Z-Transform
Properties
Important notes:
green circle (z = eiϕ ⇒ |z| = 1) reduces Z-transform simply to
discrete Fourier transform
DC (direct current) term is DC term from DFT
Re
Im
r=1
ω
Z−plane
DC term
David Svoboda (CBIA@FI) Filters in Image Processing autumn 2019 8 / 40
Z-Transform
Properties
List of the most important properties
linearity:
Z{af (n) + bg(n)} = aF(z) + bG(z)
delay (shift):
Z{f (n − k)} = F(z)z−k
convolution theorem:
Z{f (n) ∗ g(n)} = F(z) · G(z)
symmetry:
Z{f (−n)} = F(z−1
)
David Svoboda (CBIA@FI) Filters in Image Processing autumn 2019 9 / 40
Z-Transform
Some examples
Z-transform of a constant signal:
f (n) = [1, 1, 1, 1, . . . ]
F(z) = 1 + z−1
+ z−2
+ z−3
+ · · · =
1
1 − z−1
=
z
z − 1
Z-transform of a linear filter:
filter(f (k)) = (−2)f (k − 1) + 3f (k) + 5f (k + 2)
= [−2, 3 , 0, 5] /written as a kernel/
Filter(z) = (−2)z−1
+ 3z0
+ 5z2
Notice: Inverse for linear filters is straightforward.
David Svoboda (CBIA@FI) Filters in Image Processing autumn 2019 10 / 40
1 Motivation
2 Introduction to Z-transform
3 Analysis of FWT
4 Lifting scheme (LS)
5 Integer Wavelet transform
6 Applications
David Svoboda (CBIA@FI) Filters in Image Processing autumn 2019 11 / 40
Analysis of FWT
Let |f | = N = 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 12 / 40
Analysis of FWT
An example (Haar wavelets)
Let f = A2 = [1 4 -3 0] be an input signal:
A1 = f ∗
√
2
2
,
√
2
2
(↓ 2×)
D1 = f ∗
√
2
2
, −
√
2
2
(↓ 2×)
During one transition, we perform only one averaging and difference:
A1(k) = (f (k) + f (k + 1))
√
2
2
(↓ 2×)
D1(k) = (f (k + 1) − f (k))
√
2
2
(↓ 2×)
David Svoboda (CBIA@FI) Filters in Image Processing autumn 2019 13 / 40
Analysis of FWT
An example (Haar wavelets)
>> [LoD, HiD, LoR, HiR] = wfilters(’haar’)
LoD = 0.7071 0.7071
HiD =-0.7071 0.7071
LoR = 0.7071 0.7071
HiR = 0.7071 -0.7071
>> [A1,D1] = dwt([1 4 -3 0], LoD, HiD)
A1 = 3.5355 -2.1213
D1 =-2.1213 -2.1213
>> [A0,D0] = dwt(A1, LoD, HiD)
A0 = 1.0000
D0 = 4.0000
Here c = |LoD| = |HiD| = 2.
David Svoboda (CBIA@FI) Filters in Image Processing autumn 2019 14 / 40
1 Motivation
2 Introduction to Z-transform
3 Analysis of FWT
4 Lifting scheme (LS)
5 Integer Wavelet transform
6 Applications
David Svoboda (CBIA@FI) Filters in Image Processing autumn 2019 15 / 40
Lifting scheme (LS)
Developed in 1996 by Wim Sweldens
Adopted idea of FastDWT
The transition from level j to j − 1 is computed efficiently
Basic idea: split → predict → update → normalize
j−1A
j−1D
predictsplit update
even
odd
−
+
jA
N
NA
D
David Svoboda (CBIA@FI) Filters in Image Processing autumn 2019 16 / 40
Lifting scheme
Split step
Also known as lazy wavelet.
The signal Aj is split into odd and even samples:
(Aj−1, Dj−1) = split(Aj )
j−1A
j−1D
predictsplit update
even
odd
−
+
jA
N
NA
D
David Svoboda (CBIA@FI) Filters in Image Processing autumn 2019 17 / 40
Lifting scheme
Prediction step
j−1A
j−1D
split update
even
odd
−
+
jA
N
NA
D
predict
Dj−1(k) = Dj−1(k) − predict(Aj−1(k))
Also known as dual lifting step
When using Haar wavelets the neighbouring samples are supposed to
be equal, i.e. the predictor is simple: predictHaar (f (k)) = f (k)
Dj−1(k) = Dj−1(k) − Aj−1(k)
David Svoboda (CBIA@FI) Filters in Image Processing autumn 2019 18 / 40
Lifting scheme
Update step
j−1A
j−1D
predictsplit update
even
odd
−
+
jA
N
NA
D
Aj−1(k) = Aj−1(k) + update(Dj−1(k))
Also known as primal lifting step
Update step repairs the wrong estimate of the prediction step.
When using Haar wavelet, we use: updateHaar (f (k)) = 1
2f (k)
Aj−1(k) = Aj−1(k) + 1
2Dj−1(k)
David Svoboda (CBIA@FI) Filters in Image Processing autumn 2019 19 / 40
Lifting scheme
Normalization step
j−1A
j−1D
predictsplit update
even
odd
−
+
jA
N
NA
D
Aj−1(k) = Aj−1(k) · NA
Dj−1(k) = Dj−1(k) · ND
NA · ND = 1
The output signals are normalized to avoid boosting the signal.
David Svoboda (CBIA@FI) Filters in Image Processing autumn 2019 20 / 40
Lifting scheme
Inverse lifting
The forward algorithm can be simply inverted:
Aj−1(k) = Aj−1(k) − update(Dj−1(k))
Dj−1(k) = Dj−1(k) + predict(Aj−1(k))
Aj = merge(Aj−1, Dj−1)
jA
j−1A
j−1D
ND
AN +
−
update predict
even
odd
merge
David Svoboda (CBIA@FI) Filters in Image Processing autumn 2019 21 / 40
Lifting scheme
Inverse lifting (example)
Namely for Haar wavelets we get:
Aj−1(k) = Aj−1(k) · ND
Dj−1(k) = Dj−1(k) · NA
Aj−1(k) = Aj−1(k) − (Dj−1(k)/2)
Dj−1(k) = Dj−1(k) + Aj−1(k)
Aj = merge(Aj−1, Dj−1)
jA
j−1A
j−1D
ND
AN +
−
update predict
even
odd
merge
David Svoboda (CBIA@FI) Filters in Image Processing autumn 2019 22 / 40
Lifting scheme
From Filters to Lifting
There exists algorithm invented by Wim Sweldens (1996):
1 Input: either ’LoD’ and ’HiD’ filters or φ and ϕ functions
2 Convert both filters ’LoD’ and ’HiD’ to Z-domain
3 Create polyphase matrix (2×2)
4 Factorize matrix into simple (lower and upper diagonal) matrices
5 Each simple matrix correspond either to one update or prediction step
6 Convert each matrix from Z-domain to time domain
The algorithm is quite tricky → we utilize Matlab function ’liftwave’
that performs steps 1-5. All we have to do is to enjoy step 6
David Svoboda (CBIA@FI) Filters in Image Processing autumn 2019 23 / 40
Lifting scheme
From Filters to Lifting (Haar)
>> h = liftwave(’haar’)
h = ’d’ [ -1] [0]
’p’ [0.5000] [0]
[1.4142] [0.7071] []
’d’ . . . dual lifting step (predict)
Predict(z) = (−1) · z0 → predict(f (k)) = (−1) · f (k)
’p’ . . . primal lifting step (update)
Update(z) = 0.5 · z0 → update(f (k)) = 0.5 · f (k)
[1.4142] . . . NA
[0.7071] . . . ND
David Svoboda (CBIA@FI) Filters in Image Processing autumn 2019 24 / 40
Lifting scheme
From Filters to Lifting (Daubechies-2)
In general, there may appear two or more predict(.) and update(.) steps:
>> h = liftwave(’db2’)
h = {...
’d’ [ -1.732] [0]
’p’ [ -0.067 0.433] [1]
’d’ [ 1.000] [-1]
[ 1.932] [ 0.518] []
’d’ . . . dual lifting step (predict)
Predict1(z) = (−1.732) · z0 → predict1(f (k)) = (−1.732) · f (k)
’p’ . . . primal lifting step (update)
Update1(z) = (−0.067) · z1 + 0.433 · z0 →
update1(f (k)) = (−0.067) · f (k + 1) + 0.433 · f (k)
’d’ . . . dual lifting step (predict)
Predict2(z) = 1.000 · z−1 → predict2(f (k)) = 1.000 · f (k − 1)
[1.932] . . . NA
[0.518] . . . ND
David Svoboda (CBIA@FI) Filters in Image Processing autumn 2019 25 / 40
Lifting scheme
From Filters to Lifting (Daubechies-2)
NA
j−1A
ND
j−1D
predictsplit update
even
odd
−
+
jA
−
predict
predict1(f (k)) = (−1.732) f (k)
update1(f (k)) = 0.433 f (k) − 0.067 f (k + 1)
predict2(f (k)) = f (k − 1)
David Svoboda (CBIA@FI) Filters in Image Processing autumn 2019 26 / 40
Lifting scheme
Technical/Implementation notes
Lifting ordering (for N = 8)
f(0) f(1) f(2) f(3) f(4) f(5) f(6) f(7)
A2(0) D2(0) A2(1) D2(1) A2(2) D2(2) A2(3) D2(3)
A1(0) D1(0) A1(1) D1(1)
A0(0) D0(0)
Notice: The computation is performed completedly in-place.
From filters to lifting scheme
conversion ’LoD’,’HiD’ → update(·), predict(·) always exists but is
not unique
conversion is performed in frequency domain via Z-transform
(decomposition of polyphase matrices)
David Svoboda (CBIA@FI) Filters in Image Processing autumn 2019 27 / 40
Lifting scheme
FWT and DWT comparison (examples)
Price of one decomposition level using DWT (|f | = N)
family of wavelets multiplications additions
Haar 4N 2N
Daubechies-2 8N 6N
Extra memory usage: one memory buffer of size O(N) needed for convolution.
Price of one decomposition level using LS (|f | = N)
family of wavelets multiplications additions
Haar 2N N
Daubechies-2 3N 2N
Extra memory usage: ∅
David Svoboda (CBIA@FI) Filters in Image Processing autumn 2019 28 / 40
1 Motivation
2 Introduction to Z-transform
3 Analysis of FWT
4 Lifting scheme (LS)
5 Integer Wavelet transform
6 Applications
David Svoboda (CBIA@FI) Filters in Image Processing autumn 2019 29 / 40
Integer Wavelet Transform
Basic idea & Properties
IWT originates from lifting scheme (chain of predictions and updates).
The fixed point arithmetic is guaranteed via ’floor’ function.
The rounding error produced in forward transform is compensated by
mirror ’floor’ in inverse lifting.
The lifting is the same as the standard one but
each multiplication is followed but the truncation
no normalization is present
David Svoboda (CBIA@FI) Filters in Image Processing autumn 2019 30 / 40
Integer Wavelet Transform
An example (Haar)
j−1D
j−1A
split
even
odd
−
+
jA predict update
floor floor
& &
(Aj−1, Dj−1) = split(Aj )
Dj−1(k) = Dj−1(k) − floor (Aj−1(k))
Aj−1(k) = Aj−1(k) + floor (Dj−1(k)/2)
David Svoboda (CBIA@FI) Filters in Image Processing autumn 2019 31 / 40
Integer Wavelet Transform
An example (Haar)
>> h2 = liftwave(’haar’, ’Int2Int’)
h2 = ’d’ [ -1] [0]
’p’ [0.5000] [0]
[1.4142] [0.7071] ’I’
>> [A1,D1] = lwt([1 4 -3 0], h2)
A1 = 2 -2
D1 = 3 3
>> ilwt(A1,D1,h2)
ans = 1 4 -3 0
David Svoboda (CBIA@FI) Filters in Image Processing autumn 2019 32 / 40
1 Motivation
2 Introduction to Z-transform
3 Analysis of FWT
4 Lifting scheme (LS)
5 Integer Wavelet transform
6 Applications
David Svoboda (CBIA@FI) Filters in Image Processing autumn 2019 33 / 40
Applications
Image fusion – in confocal microscopy
⇓
source: B. Z´ıtov´a, UTIA, CAS
David Svoboda (CBIA@FI) Filters in Image Processing autumn 2019 34 / 40
Applications
Image fusion – in photography
⇓
David Svoboda (CBIA@FI) Filters in Image Processing autumn 2019 35 / 40
Applications
Noise removal
source: B. Z´ıtov´a, UTIA, CAS
David Svoboda (CBIA@FI) Filters in Image Processing autumn 2019 36 / 40
Applications
Image compression
original (979 kB) JPEG (6.21 kB) JPEG2000 (1.83 kB)
David Svoboda (CBIA@FI) Filters in Image Processing autumn 2019 37 / 40
Applications
. . . and the others
fusion of images with different resolution
image registration
edge detection
. . .
David Svoboda (CBIA@FI) Filters in Image Processing autumn 2019 38 / 40
Bibliography
Li H., Manjunath B.S., Mitra S.K. Multisensor Image Fusion Using
the Wavelet Transform, Graphical Models and Image Processing,
Volume 57, Issue 3, May 1995, Pages 235-245, ISSN 1077-3169
Jensen A., La Cour-Harbo A. Ripples in mathematics: the discrete
wavelet transform, Springer, Berlin, 2001, ISBN 3-540-41662-5
Sweldens W. The lifting scheme: A custom-design construction of
biorthogonal wavelets. Applied and Computational Harmonic
Analysis. 1996, Vol 3, Issue 2, pp 186200
David Svoboda (CBIA@FI) Filters in Image Processing autumn 2019 39 / 40
You should know the answers . . .
Describe the relationship between Fourier transform and Z-transform.
Can you apply Z-transform to a given signal or a linear filter? Give an
example.
Explain three principal phases of lifting scheme.
Compare the time and spatial complexity of FWT and lifting scheme.
Compute one level of wavelet transform (using Haar’s basis) via lifting
scheme for signal f=[3 5 0 -1 4 2].
What does integer wavelet transform mean?
How does image fusion via DWT work?
How does image denoising via DWT work?
David Svoboda (CBIA@FI) Filters in Image Processing autumn 2019 40 / 40