Integral and Discrete Transforms in Image Processing
Image Restoration
David Svoboda, Marek Kasfk
email: svoboda@fi.muni.cz Centre for Biomedical Image Analysis Faculty of Informatics, Masaryk University, Brno, CZ
CBIA
November 15, 2022
David Svoboda, Marek Kašik (CBIAOFI)
Image Transforms
autumn 2022
1/42
Outline
Q Motivation
• Blur o Noise
Q Restoration
• Fundamentals
• Unconstrained
• Constrained
• Iterative
David Svoboda, Marek Kašik (CBIAOFI)
Image Transforms
autumn 2022
Motivation • Blur o Noise
jii
• Fundamentals
• Unconstrained
• Constrained
• Iterative
David Svoboda, Marek Kašik (CBIAOFI)
Image Transforms
autumn 2022
Motivation
Images in optical microscopy are affected by blur and by noise. This blur is almost visible in z-axis. Images also suffer from noise due to low light intensities in confocal imaging.
Cell at the start of
optical path Cell after acquisition
David Svoboda, Marek Kašik (CBIAOFI)
Image Transforms
autumn 2022
4/42
Motivation
The task of image restoration = transforming the acquired image to its original form.
Acquired cell Restored cell image
David Svoboda, Marek Kašik (CBIAOFI)
Image Transforms
autumn 2022
5/42
Motivation
Image Blur
The blur can be described by Point Spread Function (PSF). PSF is response of optical setup to an infinitely small point source of light placed to the input. All the points are influenced by this function.
Blur can be caused by different sources: O Move of the camera during acquisition Q Defocus O Physical limits
David Svoboda, Marek Kašik (CBIAOFI)
Image Transforms
autumn 2022
6/42
Motivation
Image Blur
Examples
line - move        disc - defocused camera   Airy disc - physical limit
David Svoboda, Marek Kašik (CBIAOFI)
Image Transforms
autumn 2022
7/42
Motivation
Noise
Noise is present in almost every real image. It can be caused by, for example:
o Environmental conditions during image acquisition (temperature)
• Quality of the sensing elements (hot pixels) 9 Interference during image/data transmission
Noise is "a random change" of pixel values. Our interest is usually focused on the three basic types of noise:
• Additive noise (amplifiers)
• Impulse noise (hot/cold pixels in CCD)
• Poisson noise (photon-shot noise)
David Svoboda, Marek Kašik (CBIAOFI)
Image Transforms
autumn 2022
8/42
Motivation
Noise
Examples
Additive noise Impulse noise
David Svoboda, Marek Kašik (CBIAOFI)
Image Transforms
autumn 2022
9/42
Motivation
Noise
Additive noise
The most common type of noise. Gray values and noise are independent:
g = f + n
where f is the original image, n is the noise, and g is the noisy image
Noise n may have different distributions:
• Gaussian distribution (amplifiers) 9 Rayleigh distribution (radar)
• Exponential distribution (laser imaging)
• Gamma distribution (laser imaging)
David Svoboda, Marek Kašik (CBIAOFI)
Image Transforms
autumn 2022
10/42
Motivation
Noise
Poisson Noise
o Important type of noise in CCD imaging (photon-shot noise, thermal noise)
9 Poisson noise is not additive and depends on the signal.
• The noisy image f of the original one g is given by random Poisson process, which describes photon collection for each pixel position
(U).
• Let photons occur at CCD pixel (/,_/) with the average rate gij (photons/s) and let the original image is observed for t seconds. Then the probability that exactly X photons have occurred at pixel (ij) is given by
Axe
P(X) =
AX.-A
where A = gjj ■ t and X = 0,1,2,...
David Svoboda, Marek Kašik (CBIAOFI)
Image Transforms
autumn 2022
11/42
• Blur
• Noise
Q Restoration
• Fundamentals
• Unconstrained
• Constrained
• Iterative
David Svoboda, Marek Kašik (CBIAOFI)
Image Transforms
autumn 2022
Restoration
Fundamentals
Digital image restoration tries to restore original image from acquired image WITH the knowledge of characteristics of degradations.
Examples of restorations: 9 intensity correction
• chromatic aberration correction
• deconvolution
David Svoboda, Marek Kašik (CBIAOFI)
Image Transforms
autumn 2022
13/42
Deconvolution
Deconvolution is an inverse process to convolution. It tries to remove blur from image. This inverse process has to deal with noise as well.
Classification of deconvolution methods:
=4> According to linearity of image processing:
• Linear - linear filtering is performed
• Non-linear - non-linear filtering is performed
=4> According to knowledge of PSF:
• Blind Deconvolution - PSF is unknown
• Non-blind Deconvolution - PSF is known
David Svoboda, Marek Kašik (CBIAOFI)
Image Transforms
autumn 2022
14/42
Deconvolution
Image degradation:
original image		observed		
	DISTORTION	image	RESTORATION	estimate
f		g		f
n
Deconvolution tries to invert degradation of an image. Such process is possible only in some cases. Sometimes, the image is irrecoverably damaged and we cannot restore most of details.
Notice: Blur removes some frequencies from the image. In this case, we cannot restore these frequencies.
David Svoboda, Marek Kašik (CBIAOFI)
Image Transforms
autumn 2022
15/42
Deconvolution
Example
An example of blurred image may be an image acquired in confocal microscopy. Picture is blurred with the point-spread function (PSF) of microscope.
1 - & -	0	II * #1	1 1
- * i		~l ■ II	
Cell at the start of      Cell after acquisition        Acquired image optical path (blur + noise) after deconvolution
David Svoboda, Marek Kašik (CBIAOFI)
Image Transforms
autumn 2022
16/42
Deconvolution
Example
Point-spread function used in the previous slide:
1	i w
%	■
Point-spread function of confocal microscope
David Svoboda, Marek Kašik (CBIAOFI)
Image Transforms
autumn 2022
17/42
Convolution x Matrix multiplication
Block circulant matrix
A =
1
3
2 4
Aext —
B =
-1 -2
1
2
1 2 3    4 0 0    0 0
C = A* B =
Bext —
-1
-2 0
-1 -5 -6
1 0
2 0 0 0
1 2
3 8
2 8
Cp = Bb • Ap =
■ -1	0	1	0	0	0	-2	0	2 ■
1	-1	0	0	0	0	2	-2	0
0	1	-1	0	0	0	0	2	-2
-2	0	2	-1	0	1	0	0	0
2	-2	0	1	-1	0	0	0	0
0	2	-2	0	1	-1	0	0	0
0	0	0	-2	0	2	-1	0	1
0	0	0	2	-2	0	1	-1	0
0	0	0	0	2	-2	0	1	-1 .
	- 1 -		- -1 -
	2		-1
	0		2
	3		-5
•	4	—	-3
	0		8
	0		-6
	0		-2
	_ 0 _		8 _
• Cp, Ap ... linearized versions of Cext and Aext
• Matrix Bb ... block circulant version of Bext
David Svoboda, Marek Kašik (CBIAOFI)
Image Transforms
autumn 2022
18/42
Convolution x Matrix multiplication
Block circulant matrix
Any convolution
g = h * f
can be written as a matrix multiplication. The above equation can be written as
g = Hf
where
• g and f are understood as linearized matrices, i.e. vectors
• H is block circulant version of matrix h
Notice: The math background and complexity is the same as in the case of the convolution. It is only a notation.
David Svoboda, Marek Kašik (CBIAOFI)
Image Transforms
autumn 2022
19/42
Unconstrained Restoration
Unconstrained restoration is a base method for image deconvolution. It is very fast and applicable to data without noise. This method assumes following image formation process:
g = Hf
g H
f
f
acquired image (including degradation) point spread function ideal image we are searching for our estimate of f
We seek to minimize the function
W(0 = WW = U - Hf'f = {g~ Hf)T(g ~ Hf)
where a = v a1 a is the Euclidean norm of vector and e(f) — g — Hf a vector of residuals.
is
David Svoboda, Marek Kašik (CBIAOFI)
Image Transforms
autumn 2022
20/42
Unconstrained Restoration
Setting the derivative of W(f) to zero with respect to f produces:
dW(f) ouT
AX
df
and solving for f yields:
= -2H1 (g - Hf) = 0
f = H-Xg
which can be written (in the case of space-invariant PSF by using of convolution theorem):
f = rr-i I mil)
Disadvantages:
• Problems with zero-amplitude frequencies.
• The approach doesn't deal with noise.
• Inverse matrix may not exists.
David Svoboda, Marek Kašik (CBIAOFI)
Image Transforms
autumn 2022
21/42
Wiener Filtering
Cross-correlation (revision)
Cross-correlation
Cross-correlation function Rfg indicates the relative degree to which two functions f and g agree for various amounts of misalignment. It is given by
oo
Rfg(r)=  f f(t)g(t + r)dt
—OO
Auto-correlation
Auto-correlation is special case of cross-correlation:
oo
Rf(r)= J f(t)f(t + r)dt
—OO
David Svoboda, Marek Kašik (CBIAOFI)
Image Transforms
autumn 2022
22/42
Wiener Filtering
Fourier analysis (revision)
Power spectrum
Power spectrum V of signal f is the FT of autocorrelation of f. FT[Rf{r)] = FT[f(t) * f(-t)] = ?(k)T{-k) = F(k)T*(k) = \F(k)\2 = Vf(k)
Lena Power spectra of Lena
David Svoboda, Marek Kašik (CBIAOFI)
Image Transforms
autumn 2022
23/42
Wiener Filtering
Derivation
We are searching for filter 1/1/ that suppresses noise and keeps the signal quality:
f=|/|/*(f + n)=l/l/*s
• f is noise-free input
• 1/1/ is unknown filter
• n is additive noise
• f is signal s = f + n filtered with W
We try to minimize:
oo
oo
MSE(W)   =    / e2(t)dt= / (f(t)-f(t)) dt =
—OO
—oo
David Svoboda, Marek Kašik (CBIAOFI)
Image Transforms
autumn 2022
24/42
Wiener Filtering
Derivation
oo
MSE( W) =
f(t)-f(t)) dt
—OO
oo
f2(t) - 2f(t)f(t) + f2(t)
dt = Ti + T2 + T3
— OO
where Tl =
oo
J f(t)2dt = Rf(0)
— OO
oo
oo
oo
T2    =    -2      f(t)f(t)dt = -2 / f(t)[(W * s)(t)]dt = -2 / W(t)Rfs(t)dt
■oo
■oo
•oo
oo
oo
oo oo
T3    =     / P(t)dt=  / (W * s)(t)(W * s)(t)dt =
W(t)W(r)Rs(r - t)drdt
•oo
■oo
■oo —oo
David Svoboda, Marek Kašik (CBIAOFI)
Image Transforms
autumn 2022
25/42
Wiener Filtering
Derivation
oo
oo oo
MSE(W)   =   /?f(0)-2 j W(t)Rf5(t)dt+ J J W(ü)W(t)/?s(t - t)drdt
—oo
—oo —oo
We try to minimize MSE(W), i.e.
dMSE( W) dW
0
Solution
Rfs=W*Rs   ^FT^   Vfs = FT(W)-Vs   ->■   FT(W) =
Vfs
David Svoboda, Marek Kašik (CBIAOFI)
Image Transforms
autumn 2022
26/42
Wiener Filtering
Provided noise and signal are NOT correlated the transfer function of filter 'W minimizing MSE is defined:
FT(W) =
1^(01
Vfs =
Vs     Vf + Vn     |FT(0l2 + |FT(n)|2
where:
• Pf is power spectrum of noise-free signal
• Pn is power spectrum of noise
Notice:   'W is called Wiener filter.
J
David Svoboda, Marek Kašik (CBIAOFI)
Image Transforms
autumn 2022
27/42
Wiener Deconvolution
Unconstrained restoration followed by Wiener filtering we get optimal deconvolution respecting the MSE condition:
f = FT
-l
(FT{S) Vf
\FT(h) Vf + Vk
where
vk =
FT(n) FT(h)
is a cross power spectrum of noise and PSF.
There exists more practical version of Wiener deconvolution
f=FT-i (FT(g)Vg-Vn\
David Svoboda, Marek Kašik (CBIAOFI)
Image Transforms
autumn 2022
28/42
Wiener Deconvolution
An example
Lena Deconvolved Lena
David Svoboda, Marek Kašik (CBIAOFI)
Image Transforms
autumn 2022
29/42
Constrained Least Square Restoration
Idea
To suppress noise in the result of deconvolution we should impose some constraints to solution. Such constraint can be the convolution of result with some convolution kernel and the minimization of power of this image. The task is to minimize:
W(f) = \\Qff + X(\\g - Hff - \\n\\2)
g-Hf = n
where
• Q ... convolution kernel, which imposes some constraint on resulting image f
• A ... Lagrange multiplier (magic factor)
Example: Concerning Q, we can use Laplace filter which boosts high frequencies (noise).
David Svoboda, Marek Kašik (CBIAOFI)
Image Transforms
autumn 2022
30/42
Constrained Restoration
Now we need to solve equation:
= 2QTQf - 2XHT(g - Hf) = 0
df
Solving for f then yields:
f = (\HTH+Q' Q)-LXH'g
T     — 1 \ i_iT
David Svoboda, Marek Kašik (CBIAOFI)
Image Transforms
autumn 2022
31/42
Constrained Least Square Deconvolution
An example
Lena Deconvolved Lena
David Svoboda, Marek Kašik (CBIAOFI)
Image Transforms
autumn 2022
32/42
Iterative Methods
In the most cases it is not possible to find solution of problem of restoration by simple linear filtering. Hence there exist different methods for restoration, which are iterative. This methods iteratively improve initial estimate until it is acceptable.
Most common iterative methods:
• EM-MLE (Expectation Maximization - Maximum Likelihood Estimation)
• ICTM (Iterative Constrained Tikhonov-Miller algorithm)
Used criteria:
• number of iterations
• relative change between two iterations
o difference between blurred estimate and input
David Svoboda, Marek Kašik (CBIAOFI)
Image Transforms
autumn 2022
33/42
EM-MLE
EM-MLE is one of the best methods for restoration of images degraded by noise with Poisson statistics. Its principle consists in maximizing certain functional. Such maximization is performed by Expectation-Maximization method, which is iterative numerical method. The functional which is maximized is likelihood functional and is expressed by:
where
• p(x) is a probability density function following Poisson distribution
• Nj is j-th point value of acquired image
• /jl(uj) is blurred version of our estimate (Hf).
Notice:   EM-MLE method is also known as Richardson-Lucy method.
David Svoboda, Marek Kašik (CBIAOFI)
Image Transforms
autumn 2022
34/42
EM-MLE
Illustration
Xo     Xi     X2     X3     X4     X Xo     Xi     X2     X3     X4 X
On the left is acquired image with Poisson distributions on values of g(x). On the right, you can see the same distributions, but used probabilities values (green spots) are from positions of (Hf)(x) (green lines). Product of this probabilities is likelihood value.
David Svoboda, Marek Kašik (CBIAOFI)
Image Transforms
autumn 2022
35/42
EM-MLE
After derivation of this method we get simple expression for one iteration: where
• b is the known background
o Hfk and HT() are matrix-vector multiplications
a all other operations are point-wise
David Svoboda, Marek Kašik (CBIAOFI)
Image Transforms
autumn 2022
36/42
EM-MLE
Example
David Svoboda, Marek Kašik (CBIAOFI)
Image Transforms
autumn 2022
37/42
ICTM
ICTM method considers additive noise with Gaussian distribution. Due to iterative form, we can impose non-linear constraints to result of deconvolution. One of the base constraints is non-negativity of resulting intensities. This cannot be achieved by non-iterative method. In this method we seek to minimize functional of the form:
i
W(f) = -(\\Hf-g\\2+1\\Qf\\2).
Notice: The functional of this form is solvable by algorithm of conjugate gradients (see numerical methods - FI:PV027 Optimization).
David Svoboda, Marek Kašik (CBIAOFI)
Image Transforms
autumn 2022
38/42
ICTM
Algorithm
Computation consists in the following steps:
fk+! = P(fk + akdk),
where dk is direction of k — th step, ak is size of the k — th step, fk is result of the last iteration and P() is projection operator, which clips values to non-negative values. Direction dk+1 is computed in this way:
dk+i =
k 2
■k-1
-dk - r
2 5
where rk = Afk-b,A = HTH + ^QTQ and b = HTg. Termination of algorithm:
threshold >
Wi+1 - 1/1/'' W+1
where 1/1/' is value of l/l/() at iteration /
David Svoboda, Marek Kašik (CBIAOFI)
Image Transforms
autumn 2022
39/42
ICTM
An example
Acquired cell Restored cell image
David Svoboda, Marek Kašik (CBIAOFI)
Image Transforms
autumn 2022
40/42
Bibliography
• Gonzalez, R. C, Woods, R. E., Digital image processing / 2nd ed., Upper Saddle River: Prentice Hall, c2002, pages 793, ISBN 0201180758
• Castleman, R. C, Digital Image Processing, Prentice Hall, 1996
• Verveer, P. J., Computational and Optical Methods for Improving Resolution and Signal Quality in Fluorescence Microscopy, PhD thesis, Technical University Delft, 1998
David Svoboda, Marek Kašik (CBIAOFI)
Image Transforms
autumn 2022
41/42
You should know the answers
• Explain the difference between a constrained and unconstrained image restoration.
Explain the equation g = H * f + n.
9 Explain why we cannot invert the convolution theorem to eliminate the consequences of convolution with known PSF.
• Explain the meaning of the multiplication       in EM-MLE algorithm.
• What is a difference between Wiener filtering and Wiener deconvolution?
9 How do we get Vn in order we could perform the Wiener deconvolution?
9 Are you able to implement Constrained least square restoration in your favourite programming language?
• How do we stop iterative restoration methods?
David Svoboda, Marek Kašik (CBIAOFI)
Image Transforms
autumn 2022
42/42