Transcript Deconvolution, Deblurring, and Restoration

Deconvolution, Deblurring and
Restoration
T-61.182, Biomedical Image Analysis
Seminar Presentation 14.4.2005
Seppo Mattila & Mika Pollari
Overview (1/2)
• Linear space-invariant (LSI) restoration filters
- Inverse filtering
- Power spectrum equalization
- Wiener filter
- Constrained least-squares restoration
- Metz filter
• Blind Deblurring
Overview (2/2)
• Homomorphic Deconvolution
• Space-variant restoration
– Sectioned image restoration
– Adaptive-neighbourhood deblurring
– The Kalman filter
• Applications
- Medical
- Astronomical
Introduction
• Find the best possible estimate of the
original unknown image from the degraded
image.
• One typical degradation process has a form:
g x ,y
h x ,y
f x, y
n x, y
g(x,y) measured image
f(x,y) true (ideal) image
h(x,y) point spread function (PSF)
(impuse response function)
n(x,y) additive random noise
Image Restoration General
• One has to have some a priori knowledge about
the degragation process.
• Usually one needs 1) model for degragation, some
information from 2) original image and 3) noise.
• Note! Eventhough one doesn’t know the original
image some information such as power spectral
density (PSD) and autocorreletion function (ACF)
are easy to model.
Linear-Space Invariant (LSI) Restoration
Filters
• Assume: linear and shift-invariant degrading process
g x ,y
h x ,y
f x, y
n x, y
G u, v
H u, v F u ,v
n u, v
• Random noise statistically indep.
of image-generating process
• Possible to design LSI filters to
restore the image
Inverse Filtering
• Consider degrading process in matrix form: g hf
• Given g and h, estimate f by minimising the squared
error between observed image (g) and g :
g g g hf
where f and g are approximations of f and g
2
T
T
T
T
T
g g f h g g hf
• Set derivative of є2 to zero:
F u,v
G u,v
H u,v
F u,v
G u,v
H u,v
T
T
f h hf
(see Sect. 3.5.3 for details)
(if no noise)
n u,v
H u,v
(if noise present)
Inverse Filtering Examples
f'
•
•
•
•
Works fine if no noise but...
H(u,v) usually low-pass function.
N(u,v) uniform over whole spectrum.
High-freq. Noise amplified!!
G
0.4x
0.2x
Power Spectrum Equalization (PSE)
• Want to find linear transform L such that:
f x,y
L g x ,y
F u,v
i.e.
L u ,v G x, y
• Power spectral density (PSD) = FT(Autocorrelation function)
L u,v
F u,v
PSD( f (u,v)) = PSD(f(u,v))
...
1
H u,v 2
n u ,v
f u,v
H u ,v
2
H u ,v
n u,v
2
1 2
f u,v
1 2
G u,v
H u ,v
The Wiener Filter (1/2)
• Degradation model: g hf  
• Assumtions: Image and noise are secondorder-stationary random processes and they
are statistically independent
• Optimal mean-square error (MSE) criterion
Find Wiener filter (L) which minimize MSE
~
~
2
  E f  f where f  Lg


The Wiener Filter (2/2)
• Minimizing the criterion we end up to
optimal Wiener filter.
T
T
1
Lw   f h (h f h   )
• The Wiener filter depends on the
autocorrelation function (ACF) of the image
and noise (This is no problem).
• In general ACFs are easy to estimate.
Comparison of Inverse Filter,
PSE, and Wiener Filter
Constrained Least-squares Restoration
2
• Minimise: L f with constraint: g h f
where L is a linear filter operator
2
n
2
...
F u,v
H u,v
H u,v
●
●
2
2
L u,v
2
G u,v
H u,v
Similar to Wiener filter but does not require the
PSDs of the image and noise to be known
The mean and variance of the noise needed to set
optimally. If = 0
inverse filter
The Metz Filter
• Modification to inverse filter.



1  1  H (u, v)
LM 
H (u, v)
• Supress the high frequency noise instead of
amplyfying it.
• Select factor  so that mean-square error
(MSE) between ideal and filtered image is
minimized.
2
Motion Deblurring – Simple
Model
• Assume simple in plane movement during
the exposure
T
g ( x, y)   f ( x  a(t ), y  b(t ))dt
• Either PSF or MTF is needed for restoration
0
T
G (u, v)  F (u, v)  exp( j 2 (ua(t )  vb(t ))dt
0
T
 H (u, v)  exp( j 2 (ua(t )  vb(t ))dt
0
Blind Deblurring
• Definition of deblurring.
g ( x, y) PSF  f ( x, y)   ( x, y)
• Blind deblurring: models of PSF and noise
are not known – cannot be estimated
separately.
• Degragated image (in spectral domain)
consist some information of PSF and noise
but in combined form.
Method 1 – Extension to PSE
• Broke image to M x M size segment where
M is larger than dimensions of PSF then
• Average of PSD of these segments tend
toward the true signal and noise PSD
• This is combined information of blur
function and noise which is needed in PSE
• Finaly, only PSD of image is needed
Extension to PSE Cont...
g l (m, n) psf (m, n)  f l (m, n)  l (m, n) wherel  1,...,Q 2
 gl (u , v) | H (u , v) |2  fl (u , v)  l (u , v)
Q2
1
~
2 ~


)
v
,
u
(

|
)
v
,
u
(
H
|
 (u , v )
f
2 
Q l 1
is denominat or of LPSE :
Average:


 f (u , v)
LPSE (u , v)  

2
)
v
,
u
(


)
v
,
u
(

|
)
v
,
u
(
H
|



f
1/ 2
Method 2 – Iterative Blind
Deblurring
• Assumptation: MTF of PSF has zero phase.
• Idea: blur function affects in PSD but phase
information preserves original information
from edges.
g ( x, y) psf ( x, y)  f ( x, y)   ( x, y)
set noise to zero and look themag./phaseof FT
M g (u, v)  M PSF (u, v)Mf (u, v) and  g (u, v)   f (u, v)
Iterative Blind Deblurring Cont...
• Fourier transform of restored image is
~
~
F(u, v)  Mf (u, v)exp(jg (u, v))
• Note that smoothing operator S[] has small
effect to smooth functions (PSF). This leads
to iterative update rule
~ l1
Mf  M g
~l
S[ M f ]
S[ M g ]
~0
M f  M g  1 or
where
Mg
~0
Mf  M g 
S[ M g ]
Examples of Iterative Blind
Deblurring
Homomorphic deconvolution
G u, v H u, v F u ,v
• Start from:
• Convert convolution operation to addition:
log G u , v
●
●
log H u , v
log F u , v
Complex cepstrum: g x , y FT 1 log G u , v
Complex cepstra related: g x , y h x , y f x , y
Practical application, however, not simple...
Steps involved in deconvolution using
complex cepstrum
Space-variant Image Restoration
• So far we have assumed that images are
spatially (and temporaly) stationary
• This is (generally) not true – at the best
images are locally stationary
• Techniques to overcome this problem:
– Sectioned image restoration
– Adaptive neighbourhood deblurring
– The Kalman filter (the most elegant approach)
Sectioned Image Restoration
• Divide image into small [P x P] rectangular,
presumably stationary segments.
g l (m, n)  h(m, n)  f l (m, n)  l (m, n)
• Centre each segment in a region, and pad
the surrounding with the mean value.
• For each segment apply separately image
restoration (e.g. PSE or wiener).
Adaptive-neighborhood deblurring (AND)
• Grow adaptive neighborhood regions:
g m,n p , q
h p,q
Centered on (m,n)
f m,n p , q
n m,n p , q
Pixel locations
within the region
• Apply 2D Hamming window to each region:
g m,n p , q w H p , q
h p ,q
f m,n p , q w H p , q
n m,n p ,q w H p ,q
• Estimate the noise spectrum:
n m,n u , v
Am ,n u , v Gm,n u , v
A is a freq. domain scale factor that depends on the
spectral characterisics of the region grown
etc.
AND segmentation
Adaptive-neighborhood deblurring
(AND) Cont…
• Frequency-domain estimate of the uncorrupted
adaptive-neighborhood region:
Fm ,n u , v
1
n
gm , n
u,v
u,v
1 2
Gm , n u , v
H u,v
• Obtain estimate for deblurred adaptive neigborhood
region f m,n(p,q) by FT-1
• Run for every pixel in the input image g(x,y)
Deblurred image
Comparison of Sectioned and
AND-technique
Kalman Filter
• Kalman filter is a set of mathematical
equations.
• Filter provides recursive way to estimate the
state of the process (in non-stationary
environment), so that mean of squared
errors is minimized (MMSE).
• Kalman filter enables prediction, filtering,
and smoothing.
Kalman Filter State-Space
• Process Eq.
f(n  1)  a(n  1, n)f(n) d
• Observation Eq.
g(n)  h(n)f(n) O
• Innovation process:
 (n)  g (n)  g~(n | G )
n1
Kalman Filter in a Nutshell (1/2)
• Data observations are available
• System parameters are known
– a(n+1,n), h(n), and the ACF of driving and
observation noise
• Initial conditions
~
f (1| G 0 )  E  f (1)  0
p (1,0)  D0 Diagonalmatrixelementsof order102
• Recursion
Kalman Filter in a Nutshell (2/2)
1)
Compute the Kalman gain K(n)
2)
Obtain the innovation process
3)
Update
4)
Compute the ACF of filtered state error
5)
Compute the ACF of predicted state error
K(n)  a(n 1,n)p (n, n 1)hT (n) [h(n)p (n, n 1)hT (n)   0 (n)]1
~
 (n)  g(n)  h(n) f (n | Gn1 )
~
~
f (n 1| Gn )  a(n 1, n) f (n | Gn1 )  K (n) (n)
p (n)  p (n, n 1)  a(n, n 1)K (n)h(n)p (n, n 1)
p (n 1, n)  a(n 1, n)p (n)aT (n 1, n)  nd (n)
Wiener Filter Restoration of
Digital Radiography
Astronomical applications
• Images blurred by atmospheric turbulence
• Observing above the atmosphere very expensive (HST)
• Improve the ground-based resolution by
– Suitable sites for the observatory (@ 4 km height)
– Real time Adaptive optics correction
– Deconvolution
Point Spread Function (PSF) in
Astronomy
Iobserved
= Ireal ⊗PSF
Easy to measure and model from several stars usually present in astroimages
Determines the spatial resolution of an image
Commonly used for image matching and deconvolution
Ideal PSF if no atmosphere
FWHM ~ 1.22x/D
< 0.1" (8m telescope)
Atmospheric turbulence broadens the PSF
Gaussian PSF with FWHM ~ 1"
Richardson-Lucy deconvolution
• Used in both fields: astronomy & medical imaging
• Start from Bayes's theorem, end up with:
fi
1
x
g x
fi x h x
h
x fi x
• Takes into account statistical fluctuations in the
signal, therefore can reconstruct noisy images!
• In astronomy the PSF is known accurately
• From an initial guess f0(x) iterate until converge
Astro-examples
Observed
PSF
Inverse filter
Richardson-Lucy