Transcript Gaussian derivatives

Mathematische Grundlagen in Vision &
Grafik (710.100)
more appropriate title:
Scale Space and PDE methods in image
analysis and processing
Arjan Kuijper
07.04.2008
Johann Radon Institute for Computational and Applied
Mathematics (RICAM)
Austrian Academy of Sciences
Altenbergerstraße 56
A-4040 Linz, Austria
[email protected]
Johann Radon Institute for Computational and Applied Mathematics: www.ricam.oeaw.ac.at
1/48
Summary of the previous time
• Observations are necessarily done through a finite aperture.
• Observed noise is part of the observation.
• The aperture cannot take any form.
• We have specific physical constraints for the early vision
front-end kernel.
• We are able to set up a 'first principle' framework from which
the exact sensitivity function of the measurement aperture
can be derived.
• There exist many such derivations for an uncommitted
kernel, all leading to the same unique result: the Gaussian
kernel.
• Differentiation of discrete data is done by the convolution
with the derivative of the observation kernel.
Johann Radon Institute for Computational and Applied Mathematics: www.ricam.oeaw.ac.at
2/48
Summary of the previous time
•
The Gaussian kernel...
•
Many functions can not be differentiated.
•
A well know variational form of regularization is given by the socalled Tikhonov regularization.
–
–
–
–
–
–
–
–
is normalized
can be cascaded
can be made dimensionless using 'natural coordinates'.
is the 'blurred version' of the Dirac Delta function.
is the results of the central limit theorem
is anisotropic if the scales are different for the different dimensions.
acts as a low-pass filter in Fourier space.
is described by the diffusion equation
– The solution, due to Schwartz, is to regularize the data by convolving them
with a smooth test function.
– Taking the derivative of this 'observed' function is then equivalent to
convolving with the derivative of the test function.
– A functional is minimized in sense with the constraint of well behaving
derivatives.
– Tikhonov regularization with inclusion of the proper behavior of all
derivatives is essentially equivalent to Gaussian blurring.
Johann Radon Institute for Computational and Applied Mathematics: www.ricam.oeaw.ac.at
3/48
Today – part 1
• Gaussian derivatives
• Natural limits on observations
• Deblurring Gaussian blur
• Multi-scale derivatives: implementations
Johann Radon Institute for Computational and Applied Mathematics: www.ricam.oeaw.ac.at
4/48
Gaussian derivatives
• Shape and algebraic structure
• Gaussian derivatives in the Fourier domain
• Zero crossings of Gaussian derivative functions
• The correlation between Gaussian derivatives
• Discrete Gaussian kernels
• Other families of kernels
Taken from B. M. ter Haar Romeny,
Front-End Vision and Multi-scale Image Analysis,
Dordrecht, Kluwer Academic Publishers, 2003.
Chapter 4
Johann Radon Institute for Computational and Applied Mathematics: www.ricam.oeaw.ac.at
5/48
Shape and algebraic structure
• When we take derivatives to x (spatial derivatives) of
the Gaussian function repetitively, we see a pattern
emerging of a polynomial of increasing order, multiplied
with the original (normalized) Gaussian function again.
Order 0
0.4
0.3
0.2
0.1
4
2
Order 1
0.2
0.1
4
2
4
Order 4
0.5
2
0.5
2
4
Order 5
2
1
1
4
20.1
0.2
4
2
4
4
2 1
2
Order 2
Order 3
0.1
0.4
0.2
20.1
0.2
0.3
0.4
2
4
4
Order 6
4
2
2
4
4
22
4
6
2
20.2
0.4
2
4
2
4
Order 7
10
5
4
4
Johann Radon Institute for Computational and Applied Mathematics: www.ricam.oeaw.ac.at
2 5
10
6/48
Hermite polynomials
• The Gaussian function itself is a common element of all
higher order derivatives.
• These polynomials are closely related to the Hermite
polynomials:
Johann Radon Institute for Computational and Applied Mathematics: www.ricam.oeaw.ac.at
7/48
Gaussian envelope
• The amplitude of the Hermite polynomials explodes for
large x, but the Gaussian envelope suppresses any
polynomial function.
30
1000
20
500
4
2
10
2
500
4
4
2
2
4
10
20
1000
30
Johann Radon Institute for Computational and Applied Mathematics: www.ricam.oeaw.ac.at
8/48
Gaussian envelope
• Due to the limiting extent of the Gaussian window
function, the amplitude of the Gaussian derivative
function can be negligible at the location of the larger
zeros.
2 108
1 108
4
2
2
4
1 108
2 108
NB: The Gaussian derivatives are not normalized.
Johann Radon Institute for Computational and Applied Mathematics: www.ricam.oeaw.ac.at
9/48
Orthogonal functions
• The Hermite polynomials belong to the family of
orthogonal functions on R w.r.t. their weight function
For n, m 0…3:
• This does not hold for the Gaussian derivatives:
Johann Radon Institute for Computational and Applied Mathematics: www.ricam.oeaw.ac.at
10/48
Fourier domain
• The Fourier transform of the derivative of a function is
(-) times the Fourier transform of the function.
Johann Radon Institute for Computational and Applied Mathematics: www.ricam.oeaw.ac.at
11/48
Zero crossings of Gaussian derivative
functions
• We can define the 'width' of a Gaussian derivative
function as the distance to the outermost zerocrossing.
Zeros of
HermiteH
Width of Gaussian
derivative in
7.5
5
10
2.5
5
Order
5
10
15
20
2.5
Order
10
20
30
40
50
5
5
7.5
• Only an estimation of the exact analytic solution can be
given.
Johann Radon Institute for Computational and Applied Mathematics: www.ricam.oeaw.ac.at
12/48
The correlation between Gaussian
derivatives
• Higher order Gaussian derivative kernels tend to
become more and more similar.
Order 20
Order 24
100
3000
2000
50
1000
6
4
2
2
4
6
6
4
2
2
4
6
1000
50
100
2000
3000
• the correlation coefficient r between two Gaussian
derivatives of order n and m:
Johann Radon Institute for Computational and Applied Mathematics: www.ricam.oeaw.ac.at
13/48
Correlation
• How does it look like:
15
n
10
5
1
0.75
Abs 0.5
r n,m 0.25
0
15
10
m
5
Johann Radon Institute for Computational and Applied Mathematics: www.ricam.oeaw.ac.at
14/48
Correlation
• For n,m 0…4 in Fourier space:
• The correlation is unity when n=m, as expected, is
negative when n-m=2, and is positive when n-m=4,
and is complex otherwise.
4
Order 20
Order 21
Order 22
600
100
200
400
50
100
200
2
2
50
4
4
2
200
100
3000
2000
1000
2
100
Order 24
4
4
2
2
200
4
4
2
1000
400
2000
600
3000
2
4
• Recall:
Johann Radon Institute for Computational and Applied Mathematics: www.ricam.oeaw.ac.at
15/48
Correlation
• The correlation coefficient between a Gaussian
derivative function and its even neighbour up quite
quickly tends to unity for high differential order:
Correlation
coefficient
Order
5
10
15
20
0.95
0.9
0.85
0.8
(taking the absolute value)
• Gaussian derivatives are not suitable as a basis.
Johann Radon Institute for Computational and Applied Mathematics: www.ricam.oeaw.ac.at
16/48
Discrete Gaussian kernels
• The optimal kernel for the discretized Gaussian kernel
is " normalized modified Bessel function of the first
kind".
• This function is almost equal to the Gaussian kernel for
>1.
2
0.2
Gauss
0.15
Bessel
0.1
0.05
2
4
6
8
Johann Radon Institute for Computational and Applied Mathematics: www.ricam.oeaw.ac.at
17/48
Other families of kernels: Gabor
• Add constraints: specific frequency:
• This gives the Gabor family of receptive fields
2
2
2 2
1 k
2 2
2
2 2
2
• In real space:
x2
2 2
x
2
2
30
20
0.04
0.04
0.02
0.02
10
10
20
30
30
20
10
10
0.02
0.02
0.04
0.04
Johann Radon Institute for Computational and Applied Mathematics: www.ricam.oeaw.ac.at
20
18/48
30
Other constraints:  scale space
• Relax the seperability constraint (the power 2 in Fourier
space)
filter
PDE:
F ¡1 (e¡k!k2s )
F ¡1 (e¡k!ks )
F ¡1 (e¡k!k2® s )
Taken from the PhD thesis of R. Duits
Johann Radon Institute for Computational and Applied Mathematics: www.ricam.oeaw.ac.at
19/48
Poisson scale space
• A comparison for
and
 = ½ (Poission kernel)
 = 1 (Gaussian kernel):
Taken from the PhD thesis of R. Duits
Johann Radon Institute for Computational and Applied Mathematics: www.ricam.oeaw.ac.at
20/48
 scale space
• A comparison for
 = ½, ¾, 1:
Taken from the PhD thesis of R. Duits
Johann Radon Institute for Computational and Applied Mathematics: www.ricam.oeaw.ac.at
21/48
Summary
• The Gaussian derivatives are characterized by the
product of a polynomial function, the Hermite
polynomial, and a Gaussian kernel.
• The order of the Hermite polynomial is the same as the
differential order of the Gaussian derivative.
• Gaussian derivatives are not orthogonal kernels.
• The Gaussian kernel is a special case of the Gabor
kernels.
• The Gaussian scale space is a special case of the
 scale spaces.
Johann Radon Institute for Computational and Applied Mathematics: www.ricam.oeaw.ac.at
22/48
Gaussian derivatives,
Natural limits on observations, deblurring,
implementations
Taken from B. M. ter Haar Romeny,
Front-End Vision and Multi-scale Image Analysis,
Dordrecht, Kluwer Academic Publishers, 2003.
Chapter 7
Johann Radon Institute for Computational and Applied Mathematics: www.ricam.oeaw.ac.at
23/48
Scale limit
• For a given order of differentiation there is a limiting
scale-size below which the results are no longer exact.
xL
1.3
1.25
1.2
1.15
1.1
1.05
0.4
0.6
0.8
1.2
• The value of the derivative starts to deviate for scales
smaller than  = 0.6.
Johann Radon Institute for Computational and Applied Mathematics: www.ricam.oeaw.ac.at
24/48
Scale
• In the Fourier domain leakage (aliasing) occurs
2
2
1.5
1.5
1
1
0.5
0.5
0.4
0.3
0.75
fft0.2
gauss 0.5
0.1
0.25
0
0
2
2
0
2
0
x
Johann Radon Institute for Computational and Applied Mathematics: www.ricam.oeaw.ac.at
2
25/48
Leakage
• Leakage occurs for small scales:
,
.5
,
0.4
0.4
0.3
0.3
0.2
0.2
0.1
0.1
1
• Error measure:
error n,
100
I
2 n Fourier gauss
,
2
0 I
2 n Fourier gauss
,
2
Johann Radon Institute for Computational and Applied Mathematics: www.ricam.oeaw.ac.at
26/48
Accuracy & order
• Relation between scale , order of differentiation n, and
accepted error (left: 5%, right 1, 5, and 10%):
n
2
10
7.5
5
1.75
2.5
0
Scale in pixels
1.5
75
error%
50
1.25
1
0.75
25
0.5
0
0.25
0.5
1
1.5
2
error n,
100
0
2
4
6
Order of differentiation
I
2 n Fourier gauss
,
2
0 I
2 n Fourier gauss
,
2
Johann Radon Institute for Computational and Applied Mathematics: www.ricam.oeaw.ac.at
8
27/48
10
Limits
• Data:
• There is a limit to the order of differentiation for a
given scale of operator and required accuracy.
• The limit is due to the no longer 'fitting' of the Gaussian
derivative kernel in its Gaussian envelop, known as
aliasing.
Johann Radon Institute for Computational and Applied Mathematics: www.ricam.oeaw.ac.at
28/48
Gaussian derivatives,
limits, deblurring, implementations
Johann Radon Institute for Computational and Applied Mathematics: www.ricam.oeaw.ac.at
29/48
Deblurring Gaussian blur
• Modeling
• Wiener Filter
• Deblurring with a scale-space approach
Taken from
A. Kuijper,
Image Restoration in Forensic Research using Minimal Total Variation and Maximum
Entropy, M.Sc. Thesis, 1995.
B. M. ter Haar Romeny,
Front-End Vision and Multi-scale Image Analysis,
Dordrecht, Kluwer Academic Publishers, 2003.
Chapter 16
Johann Radon Institute for Computational and Applied Mathematics: www.ricam.oeaw.ac.at
30/48
Model
• Suppose we know our obtained image is blurred by
some convolution:
g = a(f)
• Reconstruction
is easy:
F = A ¡1 G
(especially in Fourier Space) isn’t it?
• No 
Johann Radon Institute for Computational and Applied Mathematics: www.ricam.oeaw.ac.at
31/48
Noise
• Problem 1: The blurring kernel can become
infinitesimally small, so we get division by zero
F = 1G
A
• Solution 1: Use the complex conjugacy:
F = A¤¤ G = jA¤j G
A A
A2
Johann Radon Institute for Computational and Applied Mathematics: www.ricam.oeaw.ac.at
32/48
Wiener Filter
• Problem 2: Actually, there is always noise:
g = a(f ) + n
• Solution 2: regularize the filter:
F = j jA¤ G
A 2 +R2
0.7
0.6
• This is the Wiener filter,
R = signal-to-noise ratio.
Max at “A=R”
0.5
0.4
0.3
0.2
• Demo
0.1
2
4
6
Johann Radon Institute for Computational and Applied Mathematics: www.ricam.oeaw.ac.at
8
33/48
10
Why regularizing?
• Recall last week:
Johann Radon Institute for Computational and Applied Mathematics: www.ricam.oeaw.ac.at
34/48
Deblurring with a scale-space approach
• We create a Taylor series expansion of the scale-space
L(x,y,t) to third order around the point t=0:
L x, y, 0
L 0,0,1
x, y, 0 t
1
2
L 0,0,2
• Use the diffusion equation
1
x, y, 0 t2
L
t
6
L 0,0,3
2L
2L
x2
y2
L x, y, 0
t L 0,2,0 x, y, 0
L 2,0,0 x, y, 0
1 2
t
L 0,4,0 x, y, 0
2 L 2,2,0 x, y, 0
L 4,0,0 x, y, 0
2
1 3
t
L 0,6,0 x, y, 0
3 L 2,4,0 x, y, 0
3 L 4,2,0 x, y, 0
6
L x, y, 0
1
6
t Lxx
t3 Lxxxxxx
Lyy
3 Lxxxxyy
1
2
t2 Lxxxx
3 Lxxyyyy
2 Lxxyy
x, y, 0 t3
L 6,0,0
O t 4
x, y, 0
Lyyyy
Lyyyyyy
Johann Radon Institute for Computational and Applied Mathematics: www.ricam.oeaw.ac.at
35/48
Choosing the right scale
•
•
•
•
For deblurring, a negative time is taken.
The estimated blur is est.
The kernel size is operator.
The total deblurring distance is tdeblur
L x, y, 0
1
6
t Lxx
t3 Lxxxxxx
Lyy
3 Lxxxxyy
1
2
t2 Lxxxx
3 Lxxyyyy
2 Lxxyy
2
est
2
operator
2
Lyyyy
Lyyyyyy
• Note: It is a well known fact in image processing that
subtraction of the Laplacian (times some constant
depending on the blur) sharpens the image.
This is the first order Taylor approximation!
Johann Radon Institute for Computational and Applied Mathematics: www.ricam.oeaw.ac.at
36/48
Result
order 4
order 8
order 16
order 32
Johann Radon Institute for Computational and Applied Mathematics: www.ricam.oeaw.ac.at
37/48
• With noise:
noisyblur
order 4
order 8
order 16
order 4
order 8
order 16
Johann Radon Institute for Computational and Applied Mathematics: www.ricam.oeaw.ac.at
38/48
Summary
• With a model of the blur (and noise), one can try to
deblur images.
• Deblurring is instable, and can only be carried out
analytically when no data is lost, for example through
finite intensity representation (8 bit), noise of other
pixel errors.
• The Wiener Filter traditionally does a good job, but
requires estimation of a parameter.
• Deblurring can also be done by expanding the scale
space of a blurred image into the negative scale
direction by means of a Taylor expansion.
It avoids Fourier transforms.
Johann Radon Institute for Computational and Applied Mathematics: www.ricam.oeaw.ac.at
39/48
Gaussian derivatives,
limits, deblurring, implementations
Johann Radon Institute for Computational and Applied Mathematics: www.ricam.oeaw.ac.at
40/48
Implementations
• Implementation in the spatial domain
• Separable implementation
• Implementation in the Fourier domain
• Boundaries
Taken from
B. M. ter Haar Romeny,
Front-End Vision and Multi-scale Image Analysis,
Dordrecht, Kluwer Academic Publishers, 2003.
Chapter 5
Johann Radon Institute for Computational and Applied Mathematics: www.ricam.oeaw.ac.at
41/48
Implementation in the spatial domain
• In the spatial domain, the Gaussian is to be sampled.
• A sample range from +/- 4  suffices.
The larger scale is chosen, the larger this kernel
becomes.
• Truncation order is O(10-5).
The higher the order of differentiation, the higher the
error becomes.
• Beware of the convolution rules your programming
language uses!
Johann Radon Institute for Computational and Applied Mathematics: www.ricam.oeaw.ac.at
42/48
Separable implementation
• The fastest implementation exploits the separability of
the Gaussian kernel.
–
–
–
–
Apply convolution an the matrix representing the image
Transpose the matrix
Apply convolution again
Transpose again.
Johann Radon Institute for Computational and Applied Mathematics: www.ricam.oeaw.ac.at
43/48
Orders of differentiation
x
0,
y
0
x
1,
y
0
x
0,
y
1
x
2,
y
0
x
1,
y
1
x
0,
y
2
x
3,
y
0
x
2,
y
1
x
1,
y
2
x
0,
y
3
x
4,
y
0
x
3,
y
1
x
2,
y
2
x
1,
y
3
x
0,
y
4
x
5,
y
0
x
4,
y
1
x
3,
y
2
x
2,
y
3
x
1,
y
4
• 0
• 1
• 2
• 3
• 4
x
0,
y
5
• 5
Johann Radon Institute for Computational and Applied Mathematics: www.ricam.oeaw.ac.at
44/48
Implementation in the Fourier domain
• The spatial convolutions are not exact. The Gaussian
kernel is truncated.
• A convolution of two functions in the spatial domain is a
multiplication of the Fourier transforms of the functions
in the Fourier domain, and take the inverse Fourier
transform to come back to the spatial domain.
f x
h x
F
H
g x
. G
• The filter is computed with the same size as the image.
• Beware of the rules of your programming language for
built-in Fourier transforms!
Johann Radon Institute for Computational and Applied Mathematics: www.ricam.oeaw.ac.at
45/48
Boundaries
• What is a boundary?
Johann Radon Institute for Computational and Applied Mathematics: www.ricam.oeaw.ac.at
46/48
• Boundary effect due to periodicity in the Fourier space:
(original, x-derivative, y-derivative)
Johann Radon Institute for Computational and Applied Mathematics: www.ricam.oeaw.ac.at
47/48
What to choose?
• Tilling or mirrored tilling?
• Padding with zero’s?
• Rule of thumb: tilling & ignore things that are at a
distance  of the boundary.
Johann Radon Institute for Computational and Applied Mathematics: www.ricam.oeaw.ac.at
48/48