Transcript Digital Image Processing
Digital Image Processing
Image Restoration and Reconstruction (Linear Restoration Methods) Christophoros Nikou [email protected]
University of Ioannina - Department of Computer Science
2 Contents In this lecture we will look at linear image restoration techniques – Differentiation of matrices and vectors – Linear space invariant degradation – Restoration in absence of noise • Inverse filter • Pseudo-inverse filter – Restoration in presence of noise • Inverse filter • Wiener filter • Constrained least squares filter C. Nikou – Digital Image Processing (E12)
3 Differentiation of Matrices and Vectors Notation:
A
is a
M
x
N
matrix with elements
a ij
.
x
is a
N
x1 vector with elements
x i
.
f(
x
) is a scalar function of vector
x
.
g
(
x
) is a
M
x1 vector valued function of vector
x
.
C. Nikou – Digital Image Processing (E12)
4 Differentiation of Matrices and Vectors (cont...) Scalar derivative of a matrix.
A
is a
M
x
N
matrix with elements
a ij
.
A
t
a t
11
a M
t
1
a
1
t N
a MN
t
C. Nikou – Digital Image Processing (E12)
5 Differentiation of Matrices and Vectors (cont...) Vector derivative of a function (gradient).
x
is a
N
x1 vector with elements
x i
.
f(
x
) is a scalar function of vector
x
.
f
x
x
f
f
x
1
f
x N
T
C. Nikou – Digital Image Processing (E12)
6 Differentiation of Matrices and Vectors (cont...) Vector derivative of a vector (Jacobian):
x
is a
N
x1 vector with elements
x i
.
g
(
x
) is a
M
x1 vector valued function of vector
x
.
g
x
g x
g
x
1
M
1 1
x g N
g
x M N
1 C. Nikou – Digital Image Processing (E12)
7 Differentiation of Matrices and Vectors (cont...) Some useful derivatives.
x
is a
N
x1 vector with elements
x i
.
b
is a
N
x1 vector with elements
b i
.
x
b
It is the derivative of the scalar valued function
b
T
x
with respect to vector
x
.
C. Nikou – Digital Image Processing (E12)
8 Differentiation of Matrices and Vectors (cont...) Some useful derivatives.
x
is a
N
x1 vector with elements
x i
.
A
is a
M
x
N
matrix with elements
a ij
.
x
T
x Ax A
A
T
x
If A is symmetric:
x
T
x Ax
2
Ax
C. Nikou – Digital Image Processing (E12)
9 Differentiation of Matrices and Vectors (cont...) Some useful derivatives.
x
is a
N
x1 vector with elements
x i
.
b
is a
N
x1 vector with elements
b i
.
A
is a
M
x
N
matrix with elements
a ij
.
x
2 2
A
T
It may be proved using the previous properties.
C. Nikou – Digital Image Processing (E12)
10 Linear, Position-Invariant Degradation We now consider a degraded image to be modelled by: where
h(x, y)
is the impulse response of the degradation function ( i.e.
point spread function
blurring the image). The convolution implies that the degradation mechanism is linear and position invariant (it depends only on image values and not on location).
C. Nikou – Digital Image Processing (E12)
11 Linear, Position-Invariant Degradation (cont...) In the Fourier domain: where multiplication is element-wise.
In matrix-vector form:
g = Hf + η f
where
H
,
g
, and
η
is a doubly block circulant matrix and are vectors (lexicographic ordering).
C. Nikou – Digital Image Processing (E12)
12 Linear, Position-Invariant Degradation (cont...)
g = Hf + η
If the degradation function is unknown the problem of simultaneously recovering
f(x,y)
and
h(x,y)
is called
blind deconvolution
.
C. Nikou – Digital Image Processing (E12)
13 Linear Restoration Using the imaging system
g = Hf + η
we want to estimate the true image from the degraded observation with known degradation
H.
A linear method applies an operator (a matrix)
P
to the observation
g
to estimate the unobserved noise-free image
f
: C. Nikou – Digital Image Processing (E12)
14 Restoration in Absence of Noise The Inverse Filter When there is no noise:
g = Hf
an obvious solution would be to use the inverse filter:
P = H -1
yielding
-1 -1
C. Nikou – Digital Image Processing (E12)
15 Restoration in Absence of Noise The Inverse Filter (cont...)
-1
For a
N
x
N
image,
H
is a
N
2 x
N
2 matrix! To tackle the problem we transform it to the Fourier domain.
H
is doubly block circulant and therefore it may be diagonalized by the 2D DFT matrix
W
:
-1 H = W ΛW
C. Nikou – Digital Image Processing (E12)
16 Restoration in Absence of Noise The Inverse Filter (cont...) where
Λ
Therefore:
-1 H = W ΛW -1
H
-1 -1
-1 -1 -1 -1
)}
-1
C. Nikou – Digital Image Processing (E12)
17 Restoration in Absence of Noise The Inverse Filter (cont...) Which is the vectorized form of the DFT of the image:
-1
=
Take the inverse DFT and obtain
f(m,n)
.
Problem: what happens if
H(k,l)
has zero values?
Cannot perform inverse filtering!
C. Nikou – Digital Image Processing (E12)
18 Restoration in Absence of Noise The Pseudo-inverse Filter A solution is to set: ( , )
=
0 , , 0 0 which is a type of pseudo-inversion.
Notice that the signal cannot be restored at locations where
H(k,l)
=0.
C. Nikou – Digital Image Processing (E12)
19 Restoration in Absence of Noise The Pseudo-inverse Filter (cont...) A pseudo-inverse filter also arises by the unconstrained least squares approach.
Find the image
f
, that, when it is blurred by
H
, it will provide an observation as close as possible to
g
, i.e. It minimizes the distance between
Hf
and
g
. C. Nikou – Digital Image Processing (E12)
20 Restoration in Absence of Noise The Pseudo-inverse Filter (cont...) This distance is expressed by the norm:
J
Hf - g
2 min
f
J
J
f
0
f
Hf - g
2 0 2
H
T
Hf - g
0 2
T
H Hf
2
T
H g
T
H H
1
T
H g
C. Nikou – Digital Image Processing (E12)
21 Restoration in Presence of Noise The Inverse Filter Recall the imaging model with spatially invariant degradation and noise
g = Hf + η
C. Nikou – Digital Image Processing (E12)
22 Restoration in Presence of Noise The Inverse Filter (cont...) Applying the inverse filter in the Fourier domain: Even if we know
H(k,l)
we cannot recover
F(k,l)
due to the second term.
If
H(k,l)
has small values the second term dominates (it goes to infinity if
H(k,l)
=0!).
C. Nikou – Digital Image Processing (E12)
23 Restoration in Presence of Noise The Inverse Filter (cont...) One approach to get around the problem is to limit the ratio
G
(
k,l
) /
H
(
k,l
) to frequencies near the origin that have lower probability of being zero.
We know that
H(0,0)
is usually the highest value of the DFT.
Thus, by limiting the analysis to frequencies near the origin we reduce the probability of encountering zero values.
C. Nikou – Digital Image Processing (E12)
24 Restoration in Presence of Noise The Inverse Filter (cont...) Blurring degradation C. Nikou – Digital Image Processing (E12)
25 Restoration in Presence of Noise The Inverse Filter (cont...) Inverse filter with cut-off C. Nikou – Digital Image Processing (E12)
26 Restoration in Presence of Noise Wiener Filter So far we assumed nothing about the statistical properties of the image and noise.
We now consider image and noise as random variables and the objective is to find an estimate of the uncorrupted image
f
such that the mean square error between the estimate and the image is minimized: min
f
ˆ
E
(
f
f
ˆ ) 2 where
E
[
x
] is the expected value of vector
x
.
C. Nikou – Digital Image Processing (E12)
27 Restoration in Presence of Noise Wiener Filter (cont...) Recall also the definition of the correlation matrix between two vectors
x
and
y
:
R xy
E
xy
T
N
1 1 2
N
2
N
1
N N
We assume that the image and the noise are uncorrelated:
R ηf
R fη
0
C. Nikou – Digital Image Processing (E12)
28 Restoration in Presence of Noise Wiener Filter (cont...) We are looking for the best estimate: min
f
ˆ
E
(
f
f
ˆ ) 2 Let’s confine our estimate to be obtainable by a linear operator on the observation:
f
ˆ
Pg
and the goal is to find the best matrix
P
.
C. Nikou – Digital Image Processing (E12)
29 Restoration in Presence of Noise Wiener Filter (cont...)
J
f
E
(
f
f
ˆ ) 2
E
f
f
ˆ 2
E
Denoting by
p
T n
the
n
-th row of
P
:
f
Pg
2
J
f
E
n
f
n
T
p g
n
2
n E
f
n
T
p g
n
2 C. Nikou – Digital Image Processing (E12)
30 Restoration in Presence of Noise Wiener Filter (cont...)
J
f
n E
f
n
T
p g f
n
n
T
p g
n
T
n E
f f
n n T
T
p gf
n n T
T
f g p
n n
T n T
p gg p
n
n E
f f
n n T
p
T n E
gf
n T
E
f g
n T
p
n
p
T n E
gg
T
n
R
2
T
p R
n
gf
n
T
p R p
n
gg
n
p
n
C. Nikou – Digital Image Processing (E12)
31 Restoration in Presence of Noise Wiener Filter (cont...) We can now minimize the sum with respect to each term:
p
n
R
2
T
p R
n
gf
n
T
p R p
n
gg
n
0 2
R gf
n
2
R p gg
n
0
p
n
1
R R gg gf
n
p
T n
R R
1
gg
C. Nikou – Digital Image Processing (E12)
32 Restoration in Presence of Noise Wiener Filter (cont...) Assembling the rows together:
P
R R fg
1
gg
We have to compute the two matrices:
R gg
E
gg
T
E
(
Hf
)
T
E T
Hff H
T
Hfη
T
T
HR H ff
T
HR fη
R H
T
R ηη
C. Nikou – Digital Image Processing (E12)
33 Restoration in Presence of Noise Wiener Filter (cont...) Assumming noise is uncorrelated with image:
R gg
HR H ff
T
R ηη
Also,
R fg
E
fg
T
E
η
)
T
R H ff
T
Finally the matrix we are looking for is
P
1
R R fg gg
R H ff
T
HR H ff
T
R ηη
1 C. Nikou – Digital Image Processing (E12)
34 Restoration in Presence of Noise Wiener Filter (cont...) The estimated uncorrupted image is
f
ˆ
Pg f
ˆ
R H ff
T
HR H ff
T
R ηη
1
g
which may be also expressed as
f
ˆ
T
1
ηη
1
ff
1
T
1
H R g
This result is known as the Wiener filter or the minimum mean square error (MMSE) filter.
C. Nikou – Digital Image Processing (E12)
35 Restoration in Presence of Noise Wiener Filter (cont...) Special cases: No blur (
H=I, g=f+η
):
f
ˆ
R ff
R ff
R ηη
1
g
No noise (
R ηη
=
0
,
g=Ηf
): This is the inverse filter.
f
ˆ No blur, no noise (
H=I, R ηη
=
0
):
f
ˆ
g
Do nothing on the observation.
C. Nikou – Digital Image Processing (E12)
36 Restoration in Presence of Noise Wiener Filter (cont...) The size of the matrix to be inverted poses difficulties and Wiener filter is implemented in the Fourier domain.
This occurs when
H
is doubly block circulant (represents convolution) and the image
f
noise
η
are wide-sense stationary (w.s.s).
and Definition of a w.s.s. signal: 1)
E
[
f
(
m,n
)]
=μ
, independent of
m,n
.
2)
E
[
f
(
m,n
)
f(k,l
)]
=r
(
m-k,n-l
) , independent of location.
C. Nikou – Digital Image Processing (E12)
37 Restoration in Presence of Noise Wiener Filter (cont...) Reminder: the inverse DFT complex exponential matrix diagonalizes any circulant matrix:
H
1
W Λ W H
The columns of
W
-1 circulant matrix
H
.
are the eigenvectors of any The corresponding eigenvalues are the DFT values of the signal producing the circulant matrix.
Remember also that and that
W
T T
W W
1 C. Nikou – Digital Image Processing (E12)
38 Restoration in Presence of Noise Wiener Filter (cont...) We will employ the following relations:
H
T
H
1
W Λ W H
1
W Λ W H
T
WΛ W H
1 If
H H
T
is real:
H
T
(
N
W
1 )
Λ
*
H WΛ W H
1
* 1
N
*
H
*
H
* C. Nikou – Digital Image Processing (E12)
39 Restoration in Presence of Noise Wiener Filter (cont...) The Wiener solution is now transformed to
f
ˆ the Fourier domain:
R H ff
T
HR H ff
T
R ηη
1
g
( 1
ff
)( 1 *
H
1
H
)( 1
ff
)( 1
H W Λ W
) 1
g
1
ff
*
H
1 (
Λ Λ Λ H ff
*
H
Λ ηη
) 1
g
Wf
ˆ
ff H
(
H ff
*
H
Λ ηη
) 1
Wg
Notice that the matrices are diagonal.
C. Nikou – Digital Image Processing (E12)
40
Wf
ˆ
ff
Restoration in Presence of Noise Wiener Filter (cont...)
H
(
H ff
*
H
Λ ηη
) 1
Wg
S ff S ff
( , ) * ( , ) ( , ) 2
S
S ff
(
k,l
)=DFT (
R ff
(
m,n
) ) is the power spectrum of the image
f
(
m,n
)
.
S ηη
(
k,l
)=DFT (
R ηη
of the noise (
m,n η
(
m,n
)
.
) ) is the power spectrum C. Nikou – Digital Image Processing (E12)
41 Restoration in Presence of Noise Wiener Filter (cont...) If
S ff
(
k,l
) is not zero we may define the Signal to Noise Ratio in the frequency domain:
S ff S
and the Wiener filter becomes: * ( , ) 2 -1
k l
C. Nikou – Digital Image Processing (E12)
42 Restoration in Presence of Noise Wiener Filter (cont...) A well known estimate of
S ff
(
k,l
) is the periodogram (the ML estimate of
R ff
when
f
assumed Gaussian): ia
S ff
1
MN
2 In practice, as
F
(
k,l
) is unknown, we use
S
ˆ ( , )
ff k l
1
MN
2 C. Nikou – Digital Image Processing (E12)
43 Restoration in Presence of Noise Wiener Filter (cont...) C. Nikou – Digital Image Processing (E12)
44 Restoration in Presence of Noise Wiener Filter (cont...) Degraded image Inverse Wiener Noise variance one order of magnitude less.
Noise variance ten orders of magnitude less.
C. Nikou – Digital Image Processing (E12)
45 Restoration in Presence of Noise Constrained Least Squares Filter When we do not have information on the power spectra the Wiener filter is not optimal.
Another idea is to introduce a smoothness term in our criterion.
2 We define smoothness by the quantity
Qf
where
Q
is a high pass filter operator, e.g. the Laplacian, (
Q
is a doubly block circulant matrix representing convolution by the Laplacian).
2 We look for smooth solutions minimizing
Qf
C. Nikou – Digital Image Processing (E12)
46 Restoration in Presence of Noise Constrained Least Squares Filter (cont...) We have the following constrained least squares (CLS) optimization problem: Minimize
Qf
2 subject to
Hf
g
yielding the Lagrange multiplier minimization of the function:
J
Hf
g
2
Qf
2 Data fidelity term Smoothness term C. Nikou – Digital Image Processing (E12)
47 Restoration in Presence of Noise Constrained Least Squares Filter (cont...)
f
J
0
T
2 (
f
T
H H
T
Q Q
1
T
H g g
T
Q Qf
0 Parameter
λ
controls the degree of smothness: 0,
f
T
H H
1
T
H g
pseudo-inverse, ultra rough solution ,
f
0 ultra smooth solution C. Nikou – Digital Image Processing (E12)
48 Restoration in Presence of Noise Constrained Least Squares Filter (cont...) In the Fourier domain, the constrained least squares filter becomes: * ( , ) 2 2 Keep always in mind to zero-pad the images properly.
C. Nikou – Digital Image Processing (E12)
49 Restoration in Presence of Noise Constrained Least Squares Filter (cont...) Low noise: Wiener and CLS generate equal results.
High noise: CLS outperforms Wiener if
λ
is properly selected.
It is easier to select the scalar value for
λ
than to approximate the SNR which is seldom constant.
C. Nikou – Digital Image Processing (E12)
50 Restoration Performance Measures Original image
f
and restored image
f
ˆ Mean square error (MSE): MSE 1
MN
1
M m
0
N n
0 1 ( , ) 2 C. Nikou – Digital Image Processing (E12)
51 Restoration Performance Measures (cont...) The Signal to Noise Ratio (SNR) considers the difference between the two images as noise: SNR
m
1
M
N
0
n
1 0 1
M m
0
N n
0 1 2 ( , ) 2 C. Nikou – Digital Image Processing (E12)