Integral and Discrete Transforms in Image Processing
When standard convolution comes short . . .
David Svoboda
email: svoboda@fi.muni.cz
Centre for Biomedical Image Analysis
Faculty of Informatics, Masaryk University, Brno, CZ
November 13, 2022
David Svoboda (CBIA@FI) Image Transforms autumn 2022 1 / 49
Outline
1 Linear Recursive Filters
Motivation
Filter Analysis
Filter Design
Applications
Conclusion
2 Steerable Filters
Motivation
Filter Design
Conclusion
David Svoboda (CBIA@FI) Image Transforms autumn 2022 2 / 49
1 Linear Recursive Filters
Motivation
Filter Analysis
Filter Design
Applications
Conclusion
2 Steerable Filters
Motivation
Filter Design
Conclusion
David Svoboda (CBIA@FI) Image Transforms autumn 2022 3 / 49
Linear Recursive Filters
Motivation
Common properties of linear filters based on convolution
defined via convolution kernel
naive convolution complexity: O(n2)
FFT based convolution complexity: O(n log n)
Idea of an improvement
do not evaluate the convolution process separately for each pixel
include the already convolved neighbouring values into the
convolution at the next pixel
complexity: o(n log n)
David Svoboda (CBIA@FI) Image Transforms autumn 2022 4 / 49
Linear Recursive Filters
An example
Let be given a simple recursive filter:
g : h(n) = αh(n − 1) + (1 − α)f (n)
where α is a real constant, typically α ∈ ⟨0; 1⟩.
filter takes the fraction α from the previously calculated value
filter works in certain direction
left to right – causal filter (this case)
right to left – anti-causal filter
both side – non-causal filter
no convolution kernel is defined, recursion formula is used instead
David Svoboda (CBIA@FI) Image Transforms autumn 2022 5 / 49
Linear Recursive Filters
An example
Impulse response (PSF) for filter g: h(n) = αh(n − 1) + (1 − α)f (n)
0 5 10 15 20
0
0.1
0.2
0.3
0.4
0.5
0.6
α = 0.5
0 10 20 30 40 50
0
0.02
0.04
0.06
0.08
0.1
0.12
0.14
0.16
α = 0.875
Notice: PSF can be generated by passing a brief signal
f (n) = δ(n) = [1, 0, 0, . . . ] through the filter.
Question: What happens if α > 1?
David Svoboda (CBIA@FI) Image Transforms autumn 2022 6 / 49
Linear Recursive Filters
Transfer Function
Impulse response (PSF)
filter output when accepting
a very brief signal
(δ impulse)
usually represented by
convolution kernel
(g(n))
expresses how the input
signal is modified when
passed through the filter
h(n) = f (n) ∗ g(n)
(Optical) Transfer function
Fourier transform of PSF
(G(k))
G(k) = H(k)
F(k)
H(k) = F(k) · G(k)
David Svoboda (CBIA@FI) Image Transforms autumn 2022 7 / 49
Linear Recursive Filters
Definition
Finite Impulse Response (FIR) filters
defined via finite convolution kernel
g =
1
4
1 2 1 ⇒ h(k) =
i
f (k − i)g(i)
Infinite Impulse Response (IIR) filters
defined via recursion formula
h(k) =
m
j=1
bj h(k − j)+
n
i=0
ai f (k − i)
Notice: Any recursive filter can be replaced by a nonrecursive filter (with
a mask of infinite size). Its mask is given by the PSF of the recursive filter.
David Svoboda (CBIA@FI) Image Transforms autumn 2022 8 / 49
Linear Recursive Filters
Common properties
causal – recursion formula uses only previously computed values
h(k) =
m
j=1
bj h(k − j)+
n
i=0
ai f (k − i)
anti-causal – recursion formula goes from right to left
h(k) =
m
j=1
bj h(k + j)+
n
i=0
ai f (k − i)
non-causal – filter “looks” both sides
impulse response is infinite
(we do not have to crop Gaussian hat when smoothing the image)
recursive filters need not be stable in general
(recursion may cumulate small errors)
David Svoboda (CBIA@FI) Image Transforms autumn 2022 9 / 49
Linear Recursive Filters
Tasks to solve
1 How to efficiently design any recursive filter?
2 How to guarantee its stability?
3 It is possible to design a recursive version of standard non-recursive
filters like Gaussian, Sobel, Laplace, . . . ?
Notice: We can find answers for all the questions above, but we need to
be familiar with Z-transform.
David Svoboda (CBIA@FI) Image Transforms autumn 2022 10 / 49
Filter Analysis
Given a general recursive filter:
h(n) = a0f (n) + a1f (n − 1) + a2f (n − 2) + · · · +
b1h(n − 1) + b2h(n − 2) + b3h(n − 3) + . . .
where
f (n) . . . input signal
h(n) . . . output signal
Applying the substitution
a0 = 0.389
a1 = −1.558 b1 = 2.161
a2 = 2.338 b2 = −2.033
a3 = −1.558 b3 = 0.878
a4 = 0.389 b4 = −0.161
we get the following formula:
h(n) = 0.389f (n) − 1.558f (n − 1) + 2.338f (n − 2) − 1.558f (n − 3) + 0.389f (n − 4) +
2.161h(n − 1) − 2.033h(n − 2) + 0.878h(n − 3) − 0.161h(n − 4)
David Svoboda (CBIA@FI) Image Transforms autumn 2022 11 / 49
Filter Analysis
Let us apply Z-transform to the recursion formula
h(n) = a0f (n) + a1f (n − 1) + a2f (n − 2) + · · · +
b1h(n − 1) + b2h(n − 2) + b3h(n − 3) + . . .
/Z{.}/
H(z) = a0F(z) + a1z−1
F(z) + a2z−2
F(z) + · · · +
b1z−1
H(z) + b2z−2
H(z) + b3z−3
H(z) + . . .
If we state H(z) = F(z) · G(z) then:
G(z) =
H(z)
F(z)
=
a0 + a1z−1 + a2z−2 + a3z−3 + . . .
1 − b1z−1 − b2z−2 − b3z−3 − . . .
David Svoboda (CBIA@FI) Image Transforms autumn 2022 12 / 49
Filter Analysis
Let us substitute the particular values:
G(z) =
a0 + a1z−1 + a2z−2 + a3z−3 + . . .
1 − b1z−1 − b2z−2 − b3z−3 − . . .
=
0.389 − 1.558z−1 + 2.338z−2 − 1.558z−3 + 0.389z−4
1 − 2.161z−1 + 2.033z−2 − 0.878z−3 + 0.161z−4
/
z4
z4
/
=
0.389z4 − 1.558z3 + 2.338z2 − 1.588z + 0.389
z4 − 2.161z3 + 2.033z2 − 0.878z + 0.161
/factoring/
=
(z − z1)(z − z2)(z − z3) . . .
(z − p1)(z − p2)(z − p3) . . .
pi . . . poles of transfer function G
zi . . . zeros of transfer function G
David Svoboda (CBIA@FI) Image Transforms autumn 2022 13 / 49
Filter Analysis
Transfer function properties:
poles and zeros are complex numbers
each pole must lie within the unit circle of the z-plane in order to
guarantee filter stability, i.e. |pi | ≤ 1
poles and zeros uniquely define the shape of transfer function G
factoring polynomials of higher degrees is non-trivial task
David Svoboda (CBIA@FI) Image Transforms autumn 2022 14 / 49
Filter Design
From scratch
design a completely new filter with specific conditions
rather complicated
Approximation of an existing filter (e.g. Gauss, Sobel, Laplace, . . . )
analytical approach – direct computation of recursive coefficients
[Jin & Gao, 1997]
numerical approach – search for recursive coefficients by iterative
minimization
[Deriche, 1987], [Young & Vliet, 1995]
Notice: We will focus on Jin’s approach. We will try to a design recursive
version of Gaussian filter.
David Svoboda (CBIA@FI) Image Transforms autumn 2022 15 / 49
Filter Design
Jin’s Approach
The design consists of the following three steps:
1 guarantee of stability
– specify the pole-zero placement in the z-plane;
– the position of poles defines whether the filter converges or diverges.
G(z) =
(z − z1)(z − z2)(z − z3) . . .
(z − p1)(z − p2)(z − p3) . . .
2 guarantee of accuracy
– design the filter transfer function;
– Z-transform of a recursive filter is a rational function. The accuracy
corresponds to the degree of both polynomials.
G(z) =
a0 + a1z−1 + a2z−2 + a3z−3 + . . .
1 − b1z−1 − b2z−2 − b3z−3 − . . .
3 compute the recursion coefficients of the filter
David Svoboda (CBIA@FI) Image Transforms autumn 2022 16 / 49
Filter Design
Jin’s Approach
Task
Let be Gaussian filter
gσ(n) =
1
σ
√
2π
e− n2
2σ2 = k · αn2
where k = 1
σ
√
2π
and α = e− 1
2σ2 fixed terms. The Z-transform pair of
gσ(n) is:
Gσ(z) = k
+∞
n=−∞
αn2
z−n
The task is to design a recursive version of gσ(n) and Gσ(z).
David Svoboda (CBIA@FI) Image Transforms autumn 2022 17 / 49
Filter Design
Jin’s Approach
non-recursive version recursive version
Gσ(z) = k
+∞
n=−∞
αn2
z−n
Jσ(z) = k
+∞
n=−∞
?
Notice: Without lost of generality, let us split bilateral sequence Jσ(z)
into two unilateral sequences (causal & anti-causal):
Jσ(z) = J +
σ (z) + J −
σ (z),
i.e.
J +
σ (z) = k
+∞
n=0
? and J −
σ (z) = k
0
n=−∞
?
David Svoboda (CBIA@FI) Image Transforms autumn 2022 18 / 49
Filter Design
Jin’s Approach
How to design J +
σ (z), which should be a transfer function, i.e. rational
function?
1 Let us use the second and third order polynomial. We get
J +
σ (z) = k
1 + a1z−1 + a2z−2
(1 − pz−1)3
The denominator has a unique pole of order 3 which should guarantee
(|p| ≤ 1) filter stability.
2 Using polynomial division the function J +
σ (z) can be simply
converted into infinite series in power of z:
J +
σ (z) = k
1 + a1z−1 + a2z−2
(1 − pz−1)3
/·
z3
z3
/
= k
z3 + a1z2 + a2z
z3 − 3pz2 + 3p2z + p3
/polynomial division/
= k[z0
+ (3p + a1)z−1
+ (6p2
+ 3a1p + a2)z−2
+ (6a1p2
+ 3a2p + 10p2
)z−3
+ . . . ]
David Svoboda (CBIA@FI) Image Transforms autumn 2022 19 / 49
Filter Design
Jin’s Approach
J +
σ (z) = k[z0
+ (3p + a1)z−1
+ (6p2
+ 3a1p + a2)z−2
+ (6a1p2
+ 3a2p + 10p2
)z−3
+ . . . ]
Gσ(z) = k
+∞
n=−∞
αn2
z−n
= k[· · · + αz−1
+ α4
z−2
+ α9
z−3
+ . . . ]
3 Comparing z-coefficients between J +
σ (z) and Gσ(z)
z−1 : 3p + a1 = α
z−2 : 6p2 + 3a1p + a2 = α4
z−3 : 6a1p2 + 3a2p + 10p2 = α9
we finally get:
p =
α
2
3 − α2
− 9 − 6α2 − 3α4 ,
where α = e− 1
2σ2 .
David Svoboda (CBIA@FI) Image Transforms autumn 2022 20 / 49
Filter Design
Jin’s Approach
Solution – Causal Part
J +
σ (z) = k
1 + a1z−1
+ a2z−2
(1 − pz−1)3
= k
1 + a1z−1
+ a2z−2
1 + b1z−1 + b2z−2 + b3z−3
⇓ /inverted Z-transform/
h+(n) = k{f (n) + a1f (n − 1) + a2f (n − 2)}
−{b1h+(n − 1) + b2h+(n − 2) + b3h+(n − 3)}
where
b1 = −3p
b2 = 3p2
b3 = −p3
a1 = α − 3p
a2 = α4
− 3αp + 3p2
David Svoboda (CBIA@FI) Image Transforms autumn 2022 21 / 49
Filter Design
Jin’s Approach
When dealing with (anti)symmetrical filters, it is unnecessary to apply two
times the design procedure. We can simply mirror the causal part and
eliminate the central point to avoid counting it twice.
Solution – Anti-Causal Part
J −
σ (z) = k
1 + a1z+1
+ a2z+2
1 + b1z+1 + b2z+2 + b3z+3
− 1
= k
1 + a1z+1
+ a2z+2
− (1 + b1z+1
+ b2z+2
+ b3z+3
)
1 + b1z+1 + b2z+2 + b3z+3
= k
(a1 − b1)z+1
+ (a2 − b2)z+2
− b3z+3
1 + b1z+1 + b2z+2 + b3z+3
/a3 = a1 − b1, . . . /
= k
a3z+1
+ a4z+2
+ a5z+3
1 + b1z+1 + b2z+2 + b3z+3
David Svoboda (CBIA@FI) Image Transforms autumn 2022 22 / 49
Filter Design
Jin’s Approach
Solution – Anti-Causal Part
J −
σ (z) = k
a3z + a4z2
+ a5z3
1 + b1z + b2z2 + b3z3
⇓ /inverted Z-transform/
h−(n) = k{a3f (n + 1) + a4f (n + 2) + a5f (n + 3)}
−{b1h−(n + 1) + b2h−(n + 2) + b3h−(n + 3)}
where
b1 = −3p
b2 = 3p2
b3 = −p3
a3 = a1 − b1
a4 = a2 − b2
a5 = −b3
David Svoboda (CBIA@FI) Image Transforms autumn 2022 23 / 49
Filter Design
Jin’s Approach
Final Solution
h+(n) = k{f (n) + a1f (n − 1) + a2f (n − 2)}
−{b1h+(n − 1) + b2h+(n − 2) + b3h+(n − 3)}
h−(n) = k{a3f (n + 1) + a4f (n + 2) + a5f (n + 3)}
−{b1h−(n + 1) + b2h−(n + 2) + b3h−(n + 3)}
h(n) = h+(n) + h−(n)
David Svoboda (CBIA@FI) Image Transforms autumn 2022 24 / 49
Filter Design
Jin’s Approach
−100 −80 −60 −40 −20 0 20 40 60 80 100
0
0.005
0.01
0.015
0.02
0.025
0.03
standard Gaussian
recursive Gaussian
David Svoboda (CBIA@FI) Image Transforms autumn 2022 25 / 49
Recursive Filters
Simple Exercise
A uniform averaging filter
h(k) =
n−1
i=0
f (k − i)
Its computational complexity depends on the width n. The same filter can
be written in the recursive form:
h(k) = h(k − 1) + f (k) − f (k − n)
Exercise: Show (using Z-transform) that it is formally equivalent!
David Svoboda (CBIA@FI) Image Transforms autumn 2022 26 / 49
Applications
Deriche Edge Detector
Deriche used essentially the same reasoning as Canny with one
exception.
While Canny sought an optimal filter of finite width W , Deriche
derived an optimal filter of infinite width using the same optimality
criteria as Canny.
The solution is
g(x) ≈ −cxe−α|x|
−10 −5 0 5 10
−0.8
−0.6
−0.4
−0.2
0
0.2
0.4
0.6
0.8
David Svoboda (CBIA@FI) Image Transforms autumn 2022 27 / 49
Applications
Deriche Edge Detector (1D)
Let h+ and h− denote two 1D arrays. Deriche computes the 1D
gradient along one row using this recursive form:
h+
(m) = f (m − 1) − b1h+
(m − 1) + b2h+
(m − 2)
h−
(m) = f (m + 1) − b1h+
(m + 1) + b2h−
(m + 2)
|∇f (m)| = −ce−α
(h+
(m) + h−
(m))
with
b1 = −2e−α
and b2 = e−2α
The computational load is much smaller than that of the Canny filter.
The computational time is independent of the size of the smoothing
parameter α.
David Svoboda (CBIA@FI) Image Transforms autumn 2022 28 / 49
Applications
Deriche Edge Detector (2D)
Horizontal Edge Map
gv1(x, y) = f (x, y − 1) − b1gv1(x, y − 1) − b2gv1(x, y − 2)
gv2(x, y) = f (x, y + 1) − b1gv2(x, y + 1) − b2gv2(x, y + 2)
ghv (x, y) = a(gv1(x, y) − gv2(x, y))
gh1(x, y) = a0ghv (x, y) + a1ghv (x − 1, y) − b1gh1(x − 1, y)
−b2gh1(x − 2, y)
gh2(x, y) = a2ghv (x + 1, y) + a3ghv (x + 2, y) − b1gh2(x + 1, y)
−b2gh2(x + 2, y)
V (x, y) = gh1(x, y) + gh2(x, y)
David Svoboda (CBIA@FI) Image Transforms autumn 2022 29 / 49
Applications
Deriche Edge Detector (2D)
Vertical Edge Map
gv1(x, y) = f (x − 1, y) − b1gv1(x − 1, y) − b2gv1(x − 2, y)
gv2(x, y) = f (x + 1, y) − b1gv2(x + 1, y) − b2gv2(x + 2, y)
ghv (x, y) = a(gv1(x, y) − gv2(x, y))
gh1(x, y) = a0ghv (x, y) + a1ghv (x, y − 1) − b1gh1(x, y − 1)
−b2gh1(x, y − 2)
gh2(x, y) = a2ghv (x, y + 1) + a3ghv (x, y + 2) − b1gh2(x, y + 1)
−b2gh2(x, y + 2)
H(x, y) = gh1(x, y) + gh2(x, y)
David Svoboda (CBIA@FI) Image Transforms autumn 2022 30 / 49
Applications
Deriche Edge Detector (2D)
Final solution:
|∇f (x, y)| = H(x, y)2 + V (x, y)2
Where the constants in use are:
a = −(1 − e−α
)2
b1 = −2e−α
b2 = e−2α
a0 =
−a
1 − αb1 − b2
a1 = a0(α − 1)e−α
a2 = a1 − a0b1
a3 = −a0b2
and α is the only one parameter.
David Svoboda (CBIA@FI) Image Transforms autumn 2022 31 / 49
Conclusion
When designing recursive filter one meets the following tasks:
replication – given slow (but nice) non-recursive filter, how to design
its recursive counterpart
stability – whether the new filter diverges (poles |pi | > 1) or
converges (poles |pi | ≤ 1)
accuracy
polynomial degree
numerical method error
David Svoboda (CBIA@FI) Image Transforms autumn 2022 32 / 49
Conclusion
Replication
non-recursive filter PSF g(n) with its Z-transform transfer function:
G(z) =
∞
i=0
g(i)z−i
we want to design recursive filter defined using its transfer function
G(z) =
∞
i=0
g(i)z−i
=
n
i=0
ai z−i
1 −
m
j=1
bj z−j
David Svoboda (CBIA@FI) Image Transforms autumn 2022 33 / 49
Conclusion
Stability
Given a simple recursive filter:
h(n) = αh(n − 1) + f (n)
→ Z-transform (applied on both sides of the equation):
H(z) =
1
1 − αz−1
F(z)
→ Z-transform based transfer function:
G(z) =
z
z − α
Let us analyze the problem:
G(z) has one pole at z = α
checking the filter against δ impulse f (n) = [1, 0, 0, 0, . . . ] we get
h(n) = 1, α, α2, α3, α4, . . .
for |α| < 1 the filter is stable (series converges)
for |α| > 1 the filter is unstable (series diverges)
David Svoboda (CBIA@FI) Image Transforms autumn 2022 34 / 49
Conclusion
Accuracy
Given
G(z) =
n
i=0
ai z−i
1 −
m
j=1
bj z−j
we search for ai and bi :
directly – analytical approach (see the example)
iteratively – numerical minimization:
E =
c
G(z) − G(z)
2 dz
2πiz
= /energy theorem/ =
n
|g(n)−g(n)|2
David Svoboda (CBIA@FI) Image Transforms autumn 2022 35 / 49
1 Linear Recursive Filters
Motivation
Filter Analysis
Filter Design
Applications
Conclusion
2 Steerable Filters
Motivation
Filter Design
Conclusion
David Svoboda (CBIA@FI) Image Transforms autumn 2022 36 / 49
Steerable Filters
Motivation
Let us recall template-based edge detection:
The specified filter is rotated and applied n-times
We perform n convolutions
Each subsequent convolution uses kernel rotated by n/360 degrees.
Can we decrease the task complexity?
clockwise rotation clockwise rotationclockwise rotation
A A A A
A
A A
A
A A A
A
A A A A A
A
A
A
A A
A
A
B B B B B
B
B
B
B
B
BB
BB
B
B B
B
B
B
B
B
B
B
David Svoboda (CBIA@FI) Image Transforms autumn 2022 37 / 49
Steerable Filters
Motivation
Gabor filters
Gabor(x, y) = Gaussσ(x, y) · FourierBasisθ
ω(x, y)
where
ω . . . speed of waving
θ . . . orientation of the filter
σ . . . width of Gaussian envelope
The use of Gabor filters
optical flow detection
feature extraction
. . .
How to optimize their computation?
David Svoboda (CBIA@FI) Image Transforms autumn 2022 38 / 49
Steerable Filters
Definition
A steerable filter f θ(x, y) is an orientation-selective convolution kernel
used for image enhancement and feature extraction that can be expressed
via a linear combination of a small set of rotated versions of itself:
f θ
(x, y) =
M
j=1
kj (θ)f θj
(x, y)
where f θj
(x, y) are called basis functions and kj (θ) are interpolation functions.
Notice: We wish the value of M to be the lowest possible.
David Svoboda (CBIA@FI) Image Transforms autumn 2022 39 / 49
Steerable Filters
Example (1st derivative)
Task
We are looking for arbitrary oriented 1st derivative of Gaussian Gθ
1 (x, y).
Consider simple 2D Gaussian function G:
G(x, y) = e−(x2+y2
)
Let us perform the two first-order axis-oriented derivatives:
G0◦
1 (x, y) =
∂
∂x
e−(x2+y2
) = −2xe−(x2+y2
)
G90◦
1 (x, y) =
∂
∂y
e−(x2+y2
) = −2ye−(x2+y2
)
supscript . . . orientation of derivative
subscript . . . derivative order
David Svoboda (CBIA@FI) Image Transforms autumn 2022 40 / 49
Steerable Filters
Example (1st derivative – cont’d)
The first derivative of Gaussian G at any arbitrary orientation θ can be
expressed as:
Gθ
1 (x, y) = cos (θ)G0◦
1 (x, y) + sin (θ)G90◦
1 (x, y)
G0◦
1 and G90◦
1 are called basis functions
Detection of edges in image I at any orientation can be obtained by:
R0◦
1 = G0◦
1 ∗ I
R90◦
1 = G90◦
1 ∗ I
Rθ
1 = cos (θ)R0◦
1 + sin (θ)R90◦
1
Notice: A whole family of filters can be evaluated with very little cost by
first convolving the image with basis functions.
David Svoboda (CBIA@FI) Image Transforms autumn 2022 41 / 49
Steerable Filters
Example (1st derivative – cont’d)
G60◦
1 =
1
2
G0◦
1 (x, y) +
√
3
2
G90◦
1 (x, y)
R60◦
1 =
1
2
R0◦
1 (x, y) +
√
3
2
R90◦
1 (x, y)
David Svoboda (CBIA@FI) Image Transforms autumn 2022 42 / 49
Steerable Filters
Example (2nd derivative)
Task
We are looking for arbitrary oriented 2nd derivative of Gaussian Gθ
2 (x, y).
2nd derivative of Gaussian (≈ Laplacian): G0◦
2 (x, y) = (4x2 − 2)e−(x2+y2
)
G0◦
2 (x, y) G60◦
2 (x, y) G120◦
2 (x, y)
Gθ
2 (x, y) = k1(θ)G0◦
2 (x, y) + k2(θ)G60◦
2 (x, y) + k3(θ)G120◦
2 (x, y)
where
kj (θ) =
1
3
[1 + 2 cos (2 (θ − θj ))]
David Svoboda (CBIA@FI) Image Transforms autumn 2022 43 / 49
Filter Design
Algorithm
Task
Given a function f (x, y) we wish to derive its steerable version when using
the least possible number of basis functions.
1 Assume f (x, y) = W (r)PN(x, y) /r = x2 + y2/
W (r) . . . an arbitrary windowing function (e.g. Gaussian, Hamming)
PN (x, y) . . . Nth
order polynomial in x and y
2 Function f (x, y) rotated to any angle can be synthesized as a linear
combination of 2N + 1 basis functions
PN (x, y) contains only even/odd order terms → N + 1 basis function
are sufficient for synthesis.
David Svoboda (CBIA@FI) Image Transforms autumn 2022 44 / 49
Filter Design
Algorithm
3 The interpolation functions kj (θ) must hold the following:





1
eiθ
...
eiNθ





=





1 1 . . . 1
eiθ1 eiθ2 . . . eiθM
...
...
...
...
eiNθ1 eiNθ2 . . . eiNθM










k1(θ)
k2(θ)
...
kM(θ)





Use only the lines corresponding to the degree of non-zero coefficients
from PN(x, y)
4 Solve the above system. For reasons of symmetry and robustness
against noise, the angles are equally sampled in the range 0 to π.
5 f θ(x, y) = M
j=1 kj (θ)f θj (x, y), where θ = {0, π
M , 2π
M , . . . , (M−1)π
M }
David Svoboda (CBIA@FI) Image Transforms autumn 2022 45 / 49
Filter Design
Example
Task
Assume we want to make the 1st order derivative of 2D Gaussian steerable:
1 G0◦
1 (x, y) = −2xe−(x2+y2
)
W (r) = e−(x2
+y2
) . . . windowing function
PN (x, y) = −2x . . . first order odd polynomial
2 N = 1 → we need 2(= N + 1) basis functions
3 Use only the complex exponential constraints corresponding to the
degree of non-zero coefficients from PN(x, y)
eiθ = eiθ1 eiθ2
k1(θ)
k2(θ)
4 Solving the system pro particular values θ1 = 0◦ and θ2 = 90◦ we
obtain:
k1(θ) = cos (θ)
k2(θ) = sin (θ)
5 Gθ
1 (x, y) = cos (θ)G0◦
1 (x, y) + sin (θ)G90◦
1 (x, y)
David Svoboda (CBIA@FI) Image Transforms autumn 2022 46 / 49
Conclusion
All functions that are bandlimited in angular frequency, are steerable,
given enough basis functions.
The most useful functions require small number of basis functions.
David Svoboda (CBIA@FI) Image Transforms autumn 2022 47 / 49
Bibliography
Smith S. W. ”The Scientist and Engineer’s Guide to Digital Signal
Processing, 1998, www.DSPguide.com
Jin, J. S., Gao Y. Recursive Implementation of LoG Filtering,
Real-Time Imaging 3, 59-65 (1997)
Deriche R. Separable recursive filtering for efficient multi-scale edge
detection, in Proc. Int. Workshop Machine Vision and Machine
Intelligence, Tokyo, Japan, 1987, pp. 18-23
Young I. T., Vliet L. J. van. Recursive implementation of the
Gaussian filter. Signal Process. 44, 2 (Jun. 1995), 139-151
Freeman W. T., Adelson E. H. ”The Design and Use of Steerable
Filters,” IEEE Transactions on Pattern Analysis and Machine
Intelligence, vol. 13, no. 9, pp. 891-906, September, 1991
David Svoboda (CBIA@FI) Image Transforms autumn 2022 48 / 49
You should know the answers . . .
Check, how the filter g : h(n) = αh(n − 1) + (1 − α)f (n) behaves for
α ∈ {0, 0.5, 1, 1.5}.
Describe the difference between the transfer function of FIR and IIR
filters.
What is the direction of computation of causal filters?
How do we check the stability of an existing filter?
Prove that h(k) =
n−1
i=0
f (k − i) is equal to
h(k) = h(k − 1) + f (k) − f (k − n)
What is the time-complexity of recursive filters (compared to standard
FIR filters)?
How do the steerable filters speed up the computation?
Show, how to make the steerable version of the first derivative of 2D
Gaussian.
David Svoboda (CBIA@FI) Image Transforms autumn 2022 49 / 49