TALLINN UNIVERSITY OF TECHNOLOGY Institute of Cybernetics Laboratory of Systems Biology The Analysis of Richardson-Lucy Deconvolution Algorithm with Application to Microscope Images Richardson-Lucy dekonvolutisooni algoritmi anal¨ u¨ us ja rakendus mikroskoopias by Martin Laasmaa Thesis Submitted to Tallinn University of Technology for the degree of Master of Science in Natural Sciences Supervisors: Marko Vendelin, PhD Pearu Peterson, PhD Laboratory of Systems Biology Institute of Cybernetics at TUT Tallinn 2009
Transcript
TALLINN UNIVERSITY OF TECHNOLOGYInstitute of Cybernetics
Laboratory of Systems Biology
The Analysis of Richardson-Lucy Deconvolution Algorithm
with Application to Microscope Images
Richardson-Lucy dekonvolutisooni algoritmi analuus ja
rakendus mikroskoopias
by
Martin Laasmaa
Thesis
Submitted to Tallinn University of Technology
for the degree of
Master of Science in Natural Sciences
Supervisors:Marko Vendelin, PhDPearu Peterson, PhDLaboratory of Systems BiologyInstitute of Cybernetics at TUT
Tallinn
2009
Declaration:
Hereby I declare that this master’s thesis, my original investigation and achievement,
submitted for the master’s degree at Tallinn University of Technology has not been sub-
mitted for any degree or examination.
Deklareerin, et kaesolev magistritoo, mis on minu iseseisva too tulemus, on esitatud
Tallinna Tehnikaulikooli magistrikraadi taotlemiseks ja selle alusel ei ole varem taotletud
akadeemilist kraadi.
Martin Laasmaa:
Supervisor Marko Vendelin, PhD:
Co-supervisor Pearu Peterson, PhD:
Date: 15th June, 2009
Kokkuvote
Dekonvolutisoon on tohus vahend nii fluorestsents- kui konfokaalmikroskoopia digitaalsete
piltide parendamiseks. Kuigi konfokaalmikroskoobis on punkti hajuvuse funktsioon kul-
laltki vaike vorreldes fluorestsentsmikroskoobi omaga ja seetottu on pildid teravamad,
saab dekonvolutsiooni protsessi abil veelgi parandada nende kontrastsust ning vahendada
mura.
Praeguseks on valja pakutud 3D mikroskoopia jaoks mitmeid dekonvolutsiooni algo-
ritme. Antud toos kasutatakse Richardson-Lucy (RL) iteratiivset algoritmi, mis on tule-
tatud eeldades Poisson’i mura. Kuna konfokaalmikroskoobi piltidel salvestatud mura
vastab Poisson’i murale. RL-i iteratiivne protsess ei koondu dekonvoleeritavas infos
leiduva mura tottu alati soovitud tulemusele. Koonduvust saab parandada, kui RL
kombineerida taieliku variatsiooni (TV) regulariseerimisega, kasutades selle optimaalset
regulariseerimise parameetrit.
Kaesoleva too eesmargiks on leida regulariseerimise parameetri optimaalne vahemik.
Selleks loodi tehislik pilt, holmates endas mitmeid erinevaid geomeetrilisi kujundeid
erinevate intensiivsustega, mida konvoleeriti ja millele rakendati Poisson’i mura. An-
tud piltide dekonvoleerimise tulemusena, kasutades erinevaid regulatsiooni parameetri λ
vaartusi, leiti, et antud algoritm annab antud rakenduse jaoks rahuldavaid tulemusi, kui
parameeter λ on vahemikus 0.008 kuni 0.03. Leitud vahemikus vaiksemate λ vaartuste
korral on objektid ”silutud”, intensiivsuse vead jarskudel intensiivsuste uleminekutel on
aga suhteliselt vaikesed. Vastukaaluks, suuremate λ vaartuste korral antud vahemikust
saavutati tulemusi, mis viisid kull objektide aarte parema tuvastamiseni, kuid ka suure-
mate intensiivsuste koikumisteni pildi uhtlastel aladel.
Lisaks dekonvolutsiooni algoritmi rakendamisele sunteetilistel piltidel, testiti algo-
ritmi ka konfokaalmikroskoopia piltidel. Konfokaalmikroskoopia pildid olid tehtud roti
ja vikerforelli kardiomuotsuutidest, milles erinevad rakustruktuurid olid margistatud
erinevate fluorestsentsmarkeritega. Nende dekonvoleerimiseks kasutati antud algoritmi
parameetri λ = 0.01 korral, mis viis piltide kvaliteedi olulise paranemiseni. Lisaks nai-
dati, et dekonvoleerimise protsess tuleb lopetada parast teatud arv iteratsioone, vastasel
juhul voib see hakata voimendama mura ja tekitada ebareaalseid efekte piltidel.
Lahitulevikus plaanitakse dekonvoleerimise protsessi kiirendamise eesmargil teha tark-
varale taiendavaid modifikatsioone, et arvutusi labi viia NVIDIA graafikakaartidel ja/voi
arvutiklastris. Lisaks sellele antakse valja tarkvarapakett, mis baseerub antud tool.
Abstract
Deconvolution is an efficient tool for enhancing both fluorescence and confocal microscopy
images. Although in confocal microscopy the point spread function (PSF ) is rather small
and images are much sharper compared to fluorescence microscopy, deconvolution can
improve image contrast and reduce noise considerably in confocal microscopy.
Several deconvolution methods have been proposed for 3D microscopy. In this work,
we used Richardson-Lucy (RL) iterative algorithm assuming Poisson noise (because the
noise on confocal microscope images corresponds to Poisson noise). However, RL does
not always converge to a suitable solution. Convergence is improved, when combining
RL with total variation regularization, using optimal regularization parameter.
The aim of this work was to find the range of optimal values of the regularization
parameter. For that, a synthetic image representing different shapes and intensities
were generated, convolved and added Poisson noise. We deconvolved degraded images
with different λ values. As result, from iteration history, we found that the range of
TV regularization parameter λ = 0.008 to 0.03 gives satisfying results with respect to
sharpness and the homogeneity of details. At smaller values of λ, all objects are smoothed
out, but the oscillations near the sharp transitions of intensity are relatively small. In
contrast, larger values of λ lead to a better edge detection but with the larger oscillations
in the homogeneous regions.
After the test on synthetic images, we applied the deconvolution algorithm to experi-
mental data. Rat and trout cardiomyocytes were labeled with different fluorescence dyes,
recorded with a confocal microscope and deconvolved. The deconvolution algorithm with
the regularization parameter of λ = 0.01 led to significant improvement in image quality.
We showed that the deconvolution process should be stopped after a certain number of
iterations. Prolonged iterations can lead to misleading results, start amplifying noise or
producing artifacts.
In the near future, we plan to make a software bridge, which uses NVIDIA graphics
cards (CUDA support) and/or cluster computing to enhance the computational speed.
The software package composed on the basis of this work will be released as an open
source project.
Contents
Abbreviations used 6
1 Introduction 7
2 Methods 9
2.1 Description of the Iterative Process . . . . . . . . . . . . . . . . . . . . . 9
2.2 Determination of Point Spread Function . . . . . . . . . . . . . . . . . . 10
2.3 Analysis of the Deconvolution Process . . . . . . . . . . . . . . . . . . . 10
A Derivation of the Richardson-Lucy Algorithm with Total-Variation Reg-
ularization 32
B List of Available Deconvolution Software 35
Acknowledgments 37
5
Abbreviations used
WF widefield fluorescence microscopyCLSM confocal laser scanning microscopeRL Richardson-Lucy deconvolution algorithmTV total variation regularizationPSF the point spread functiono observable objecth the PSF of the optical systemi degraded object o by PSF and noiseP(i) Poisson noise with average i
λ regularization parameter for TVMSE mean square errorI I-divergencev voxel coordinates in tree dimensional image v = (i, j, k)o(n) the nth estimation of o
τ (n+1) relative change between two successive estimates o(n+1)
and o(n), see (2.6)t threshold value for the stopping criterion
f the Fourier transform of f
f the inverse Fourier transform of f
f the complex conjugate of f
f the adjoint of f , f =ˇf
6
1 Introduction
In biosciences, microscopy is an extremely useful and important method for studying
living organisms. Usually microscopy is divided to three well-known branches: opti-
cal, electron and scanning probe microscopy. In optical microscopes, the sample can
be observed or recorded, when passing light is reflected, transmitted through object, or
emitted due to fluorescence from object and guided to eye-piece(s) or an electronic cap-
ture device such as a CCD camera. In electron microscopes, the specimen is illuminated
by an electron beam, which causes specimen to emit X-rays at a characteristic frequency
and emitted rays can be detected by the electron microprobe. In scanning probe micro-
scopes, the image is formed by recording the physical probes and surface interaction as
a function of position.
Using fluorescent materials in samples, the specimen can be excited in two ways:
the entire specimen at once or scanning over the specimen with a fine beam. In the
first approach the emitted light is also recorded at once and that type of microscopy
is called widefield fluorescence microscopy (WF). In the second approach the emitted
light is collected over the specimen point by point scanning and exiting with laser beam
focusing to such points. This technique is called confocal laser scanning microscope
(confocal microscope or CLSM)
Confocal microscopy has several advantages over traditional widefield fluorescence
microscopy. The main advantage is the ability to produce in-focus images of thick spec-
imens via elimination or reduction of background information outside of the focal plane
and ability to control the depth of field (within the accuracy of Airy disk size). Despite
the advantages over WF, confocal images are still far from being perfect. There still exist
aberrations that arise from imperfections in the optical pathway, residual out-of-focus
light, and noise from electronics as well as from quantum effects.
In this work we focus on image enhancement of confocal microscope images. Each
microscope changes the apparence of specimens in a specific way. Image formation can be
described by the mathematical operation of convolution, where “true” image is convolved
with distortion effects from microscope. We use deconvolution, which is the reverse
operation of convolution, to undo the aberrations caused by the optical train and remove
contributions from of out-of-focus objects. Deconvolution takes into account microscope
parameters and the nature of noise. Therefore, it is a method that can efficiently enhance
both WF and confocal microscopy images. It can considerably improve image contrast
and reduce noise in microscopy images.
Several deconvolution algorithms have been proposed by others for 3D microscopy.
7
Deconvolution algorithms can be divided into two categories: non-iterative and iterative.
Nearest-neighbor algorithms and Wiener filtering are one of the simplest non-iterative al-
gorithms for restoring microscope images. Usually these methods do not provide optimal
image quality. This is due to noise, which is always present and non-iterative algorithms
amplify noise and alter signal amplitudes [Cannell et al., 2006]. For this reason non-
iterative approaches are replaced by nonlinear iterative deconvolution algorithms. For
example, Fourier-wavelet regularized deconvolution [Neelamani et al., 2004], Jansson-van
Cittert algorithm [Cannell et al., 2006, Abdelhak and Sedki, 1992], iterative constrained
Tikhonov-Miller [van Kempen et al., 1997], Carrington algorithm [Carrington et al.,
1995], and the Richardson-Lucy algorithm [Richardson, 1972, Lucy, 1974].
In this work, we use of the Richardson-Lucy (RL) iterative algorithm assuming Pois-
son noise [Cannell et al., 2006]. The assumption of Poisson noise comes from the fact
that in confocal microscopes a photodetection device (a photomultiplier tube or avalanche
photodiode) counts the number of photons emitted from the specimen over the period
of time when the laser beam is focused in certain point. Due to the quantum nature
of light the number of photons observed behaves as a Poisson process whose variance is
equal to the mean of counted photons.
Richardson-Lucy algorithm belongs to a class of maximum-likelihood algorithms. It is
a common deconvolution algorithm that is used in astrophysical and microscopy imaging
[Dey et al., 2006]. One problem with the RL algorithm is that, in presence of noise,
it will converge to solution which is dominated by the noise [Dey et al., 2004]. As an
option, it can be used on images that are pre-filtered [Cannell et al., 2006]. Another
option is to apply regularization terms to the RL algorithm like Tikhonov-Miller and
maximum entropy regularization [de Monvel et al., 2003]. Algorithms, which are based
on Tikhonov-Miller regularization are often used for deconvolving 3D images. Such
algorithms avoid noise amplification by causing smear of the object edges. Alternatively
to obtain a result where object borders are sharp and homogeneous areas are smooth, is
to apply total variation regularization to RL [Dey et al., 2006].
The aim of this work was to find the range of optimal values of the total variation
regularization parameter. For that, we developed a free software package that is usable
by a wider scientific community.
8
2 Methods
In microscopy, a general model of image formation consists of tree parts. First, the object
itself. Second, optical system that transfers the information between object and image
plane. Last, noise which is always present mainly due to electronical recording devices.
Interplay between optical path and fluorescent light from object can be characterized
by the point spread function (PSF ). Generally, the PSF is a system’s impulse response,
the reaction of any dynamic system in response to some external change. In microscopy,
the PSF describes the response of an imaging system to a point source or point object.
It can be computed from the physical and optical properties of the imaging system or
measured by imaging objects, which are smaller than the wavelength of light.
2.1 Description of the Iterative Process
Assuming that the noise is of Poisson nature in confocal microscope, an image formation
can be mathematically presented as follows [Dey et al., 2006]:
i = P(o ⊗ h) (2.1)
where i is the recorded image stack represented as three dimensional array where each
value corresponds to the intensity of a measured voxel, o is the object, h is the PSF ,
⊗ denotes convolution, P is Poisson noise. Our aim is to find o, knowing h and i. To
find o, iterative process can be used. One option is to use Richardson-Lucy algorithm.
From equation (2.1) a multiplicative gradient type RL algorithm for one iteration can be
derived [Richardson, 1972, Lucy, 1974, Dey et al., 2004]. In this algorithm one iteration
step is given by:
o(n+1) =
(
i
o(n) ⊗ h⊗ h
)
· o(n), (2.2)
where o(n+1) is new estimate from previous estimate o(n) whereby o(0) = i and h is the
adjoint of h. For derivation, see Appendix A.
The RL algorithm does not always converge to a suitable solution in the presence
of noise. In fact, when n → +∞ the result will only comprise of noise [Dey et al.,
2006]. In order to get a better convergence, RL is combined with total variation (TV)
regularization [Dey et al., 2004]:
o(n+1) =
(
i
o(n) ⊗ h⊗ h
)
· o(n)
1 − λ div(
∇o(n)
|∇o(n)|
) (2.3)
9
where λ is regularization parameter, div stands for divergence with respect to voxels,
∇o(n) is gradient of o(n) and |∇o(n)| is the length of ∇o(n). See Appendix A for derivation.
2.2 Determination of Point Spread Function
For deconvolution one needs the PSF of the microscope. Microscope used in this studies
was confocal laser scanning microscope Zeiss LSM510 META. The PSF for this micro-
scope was found by recording 3D images of fluorescent microspheres with a diameter of
with a voxel size of 0.027 × 0.027 × 0.191µm. After averaging the intensity profiles of
different microspheres the PSF for the system was obtained (Figure 1) [Vendelin and
Birkedal, 2008].
Figure 1: Point spread function by averaged intensity profiles of microspheres that describesa confocal Zeiss LSM510 META. Figure taken form [Vendelin and Birkedal, 2008].
2.3 Analysis of the Deconvolution Process
To study the effects of deconvolution, synthetic image with different shapes were gener-
ated, convolved with PSF (Figure 1) and degraded with Poisson noise. To quantify the
quality of the deconvolution on synthetic image, mean square error (MSE) and the I-
divergence criteria (I) between original object deconvolved images were used [Dey et al.,
2006],[van Kempen et al., 1997]:
MSEA,B =∑
v
(Av − Bv)2 (2.4)
and
IA,B =∑
v
[
Av · ln(
Av
Bv
)
− Av + Bv
]
, (2.5)
10
where A,B are images to be compared and v = (i, j, k) is the position coordinate of voxel
in three dimensional image.
For measured myocytes, the actual object is unknown and therefore MSE nor I
cannot be used as quality criteria. In addition, it is difficult to decide when is the proper
time to stop iterative deconvolution process. For these reasons, we use a relative change
between two successive estimates as a stopping criterion. The criterion is defined as
follows [Cannell et al., 2006]:
τ (n+1) =
∑
v
∣
∣
∣o(n+1)v − o
(n)v
∣
∣
∣
∑
v o(n)v
< t, (2.6)
where t is a threshold parameter (an appropriate threshold is found in the next section).
To speed up computation, the convolution operations were computed in Fourier space.
Then we can replace in equation (2.3) the convolution (⊗) with multiplication [Press
et al., 2007]:
o(n+1) =
i(
ˆo(n) · h)∨
∧
· ¯h
∨
· o(n)
1 − λ div(
∇o(n)
|∇o(n)|
) . (2.7)
Where f denotes the Fourier transform of f , f is the inverse Fourier transform and f
denotes the complex conjugate of f . The discrete Fourier transform and its inverse of
the array f in three dimensional space are given as:
fαβγ =
N1,N2,N3∑
ijk=1
exp
(
−2π√−1
(
αi
N1
+ βj
N2
+ γk
N3
))
fijk (2.8)
and
fijk =1
N1N2N3
N1,N2,N3∑
αβγ=1
exp
(
2π√−1
(
iα
N1
+ jβ
N2
+ kγ
N3
))
fαβγ, (2.9)
where i, j, k denotes to element coordinates in three dimensional index space and α, β, γ
denotes to element coordinates in three dimensional Fourier index space, where i, j, k
and α, β, γ can have values from 1 to N1, N2, N3, respectively.
2.4 Numerical Methods
Because of the digital nature of images, we used the discrete Fourier transform (DFT ).
For computing DFT the FFTW library (http://www.fftw.org/) were used. This al-
11
lowed us to parallelize the deconvolution process.
For finding div(
∇f
|∇f |
)
, a stable numerical scheme from [Dey et al., 2004] was used.
The discrete form of this scheme is:
div
( ∇f
|∇f |
)
ijk
=1
hx
∆x−
∆x+fijk
√
(∆x+fijk)2 + m(∆y
+fijk, ∆y−fijk)2 + m(∆z
+fijk, ∆z−fijk)2
+
+1
hy
∆y−
∆y+fijk
√
(∆y+fijk)2 + m(∆x
+fijk, ∆x−fijk)2 + m(∆z
+fijk, ∆z−fijk)2
+
+1
hz
∆z−
∆z+fijk
√
(∆z+fijk)2 + m(∆x
+fijk, ∆x−fijk)2 + m(∆y
+fijk, ∆y−fijk)2
,
(2.10)
where
∆x+fijk = h−1
x (f(i+1)jk − fijk), ∆x−fijk = h−1
x (fijk − f(i−1)jk),
∆y+fijk = h−1
y (fi(j+1)k − fijk), ∆y−fijk = h−1
y (fijk − fi(j−1)k),
∆z+fijk = h−1
z (fij(k+1) − fijk), ∆z−fijk = h−1
z (fijk − fij(k−1)),
(2.11)
m(a, b) =sign a + sign b
2min(|a|, |b|), (2.12)
and hx, hy, hz are voxel sizes in three dimensional space.
In boundary points the following relations are used:
f0jk = f1jk, f(Nx+1)jk = fNxjk,
fi0k = fi1k, fi(Ny+1)k = fiNyk,
fij0 = fij1, fij(Nz+1) = fijNz,
(2.13)
where Nx, Ny, Nz denote to the dimensions of matrix f .
The deconvolution algorithm was implemented in Python using NumPy and PyFFTW
package and computations were performed on a Linux/ AMD64 computer with 4 CPUs
and 16GB of RAM.
12
3 Results
In this section we report:
1. how to find the optimal range of total variation regularization parameter for our
application purposes,
2. how to find the threshold value for stopping criterion,
3. the results of deconvolved confocal microscope images.
3.1 Analysis on Synthetic Images
To test the algorithm (2.2) and (2.3), we generated artificial images representing different
geometrical shapes and intensities. The following shapes were used: a filled cylinder, a
cylinder surrounded with higher intensity cylinder, two different bars and a polygon (see
Figure 2(a)). Those shapes cover several possible intracellular structures induced by
mitochondria, nucleus, sarcoplasmic reticulum and etc. Artificial images was convolved
with a PSF that was measured on Zeiss LSM510 META [Vendelin and Birkedal, 2008]
and degraded with Poisson noise (Figure 2(b)). To quantify deconvolution algorithms,
we compared original images with deconvolved images visually and through calculation
of several norms (MSE and I).
In the first stage of the study, we applied the RL without TV regularization (λ =
0, in (2.3)) to degraded image Figure 2(b). As a result, the algorithm brought out
edges of shapes and reduced noise (Figure 2(c) and 3, blue lines). As an artifact,
the algorithm produced relatively high intensity fluctuations in regions where intensity
should be constant. For example in Figure 2(c), the filled cylinder has wavy structure
inside and intensities on edges have high values. The same effect is visible on Figure
3. The reason for this effect is that the RL algorithm does not converge to a satisfying
solution, because the inversion problem is ill-posed and maximum likelihood estimator
is not regularized [Dey et al., 2006]. As can be seen on Figure 4 (black lines: λ = 0),
I-divergence and mean square error decreases to point where the number of iteration is
around 25 to 40 and after that the norms start to increase due to noise amplification.
In the second stage of this project, we tested the RL algorithm with TV regularization
with different regularization parameter λ values. As it turned out, TV regularization
preserves objects edges and overcomes the oscillation issues that we had previously. This
result can be seen on Figure 2(d) and 3(b). For different regularization parameter values,
13
the convergence and the speed of convergence are different. As shown on Figure 4(a)
and 4(b), the process continues converging after 200 iterations for some λ values (for
example, λ = 0.01). From that, where the process continues converging, we found
that the optimal range of regularization parameter values is 0.008 < λ < 0.03 for our
application. If λ < 0.001 the algorithm started to behave like standard RL, since the
influence of TV regularization became too small. As a result, the deconvolved image
is further away from original image. In contrast, if λ is larger than 0.08, the dominant
part in the algorithm is TV regularization. Then the denominator in equation (2.3)
can become zero or negative. In the case of approaching to zero, the denominator will
be small and therefore it creates high intensity points, that will be amplified after every
iteration.
When deconvolving images of heart muscle cells, for instance, rat and trout cardiomy-
ocytes, there is no original image with which we can compare our results. For this reason
we have to find criteria to stop the deconvolution process. For example, it is possible
to use relative change (τk from equation (2.6)) between two successive estimates during
deconvolution. We found that on synthetic images sensible values for threshold t are
around 0.0001 (see Figure 5). For higher parameter λ values the threshold should be
also higher. If the threshold is set too low then the program cannot reach a point where
it should stop. Also we noticed that the τk may start to oscillate between two iterations.
Thus, from analysis of synthetic images, we found the optimal range of TV regular-
ization parameter λ = 0.008 to 0.03 for our application and a suitable threshold value
(around τk = 0.0001) for a stopping criterion. In this range, for smaller values of λ, all
objects are smoothed out, but the oscillations near the sharp transitions of intensity are
relatively small. In contrast, larger values of λ lead to a better edge detection but with
the larger oscillations in the homogeneous regions.
Next we apply the RL algorithm with TV regularization to images which were
recorded by confocal microscope. For λ, we used value in the optimal parameter range
found from this section.
14
(a) (b)
(c) (d)
Figure 2: Tests of deconvolution algorithm on a synthetic image. Section (a) shows generatedsynthetic image, (b) corresponds to the synthetic image that was convolved and degraded withPoisson noise. Image (c) is deconvolved from (b) using RL algorithm with no TV regularization.(d) is for comparison to the image (c) that is deconvolved using RL with TV regularization(regularization parameter λ = 0.01). It can be seen from (c) and (d) that deconvolution withTV regularization is giving better results – intensity in homogeneous areas is more even andsmoother on (d) than (c). (Red lines indicate cross-sections on XZ plane, see Figure 3.)
15
0 50 100 150 200 2500
50
100
150
200
250Synthetic image
Degraded image
Deconvolved image with RL
Deconvolved image with RL TV
(a)
0 50 100 150 200 2500
50
100
150
200
250Synthetic image
Degraded image
Deconvolved image with RL
Deconvolved image with RL TV
(b)
Figure 3: Intensity profiles of Figure 2. On subplot (a) are intensities that corresponds toupper red line on Figure 2 and (b) corresponds lower red line. It is easy to see the effect thatTV regularization has on RL algorithm if we compare blue (RL) and red (RL with TV, whenλ = 0.01) line.
Figure 4: Mean square error and I-divergence dependence on iteration. Both (a) and (b)show that the nearest results to original (non-degraded) image is achieved in regularizationparameter range 0.008 < λ < 0.03. For this particular case the best result was obtained whenλ = 0.01.
Figure 5: Relative change between two successive estimates during deconvolution process, τk.Note that reasonable values for threshold t are around 0.0001. Furthermore, if we compareτk values in the range 0.008 < λ < 0.03, we can see that for higher parameter λ values thethreshold t should be also higher.
18
3.2 Analysis on Confocal Microscope Images
We applied RL algorithm with TV regularization to confocal microscope images of rat
and trout cardiomyocytes which were recorded on a Zeiss LSM510 META. Both, rat
and trout myocytes were stained with Mitotracker Red CMXRos (561nm excitation,
>575nm emission, Invitrogen) and Di-8-ANEPPS (489-nm excitation, >550nm emission,
Invitrogen) fluorescent dyes.
Figure 6(a) shows an image of a rat cardiomyocyte. Here, red color marks mitochon-
dria and green sarcolemma. The image is deconvolved with regularization parameter
value λ = 0.01. For 200 iterations we computed the relative change between two con-
secutive iterations and found that τk was minimal at the 80th iteration (below threshold
value which was found previously). After that τk stared slowly to increase (Figure 8).
From that we considered 80th iteration as the best solution which was also affirmed by
visual observation. Furthermore, it is clear that the noise level is decreased as shown
in the histograms on Figure 9. The first part of the histograms describes the image
background and the nature of its noise. Firstly, as the background intensity is reduced,
the peak is shifted towards zero. Secondly, counts for the peak is gained and the width of
noise distribution is shortened. The result of deconvolution can be seen on Figure 7(a).
Note that deconvolution brought out mitochondria and sarcolemma. Moreover, fuzzy
XZ cross-section became rather clear. For particular case further iterations did not
change the result noticeably (see Figure 7(b)).
The effect of deconvolution on trout cardiomyocyte is shown on Figure 10. We
carried out 200 iterations from where it turned out that the best result was achieved
in 33 iterations (Figure 10(a)). At the 33rd iteration value of τk was minimal and
visual assessments of nearby iteration results was conducted (Figure 11). The recorded
image of trout myocyte Figure 6(b) has quite low fluorescence signal, with significantly
high background intensity and noise level. These properties are common in practice and
thus this image is a good test of the effectiveness of the deconvolution process. From
histograms of images before and after deconvolution we see how the level of background
intensity and noise level are decreased remarkably (Figure 12). By inspection of trout
images after deconvolution, we see the clear borders of the cell, the core of mitochondria,
and sarcolemma which are all resolved.
Let us inspect what will happen if deconvolution process is not interrupted after τk
is reached its minimum. Figure 7 shows results after 80th and 200th iterations and Fig-
ure 10 after 33rd and 200th iterations. As it clear from the figures, significant aberrations
were induced by prolonged deconvolution process. In Figure 10(b) the mitochondrial
19
core seems to be smeared, also the inner part starts to vanish. If we look at sarcolemma
(green) we can see new membrane-like but unrealistic structures starting to appear. This
is in contrast to deconvolved rat myocyte, where continued iterations did not change the
result as significantly. The recorded image of trout myocyte had a high background
intensity, noise level and low fluorescence signal, all of which might cause decreased of
image quality as deconvolution process continues past the optimal solution.
(a) Rat cardiomyocyte (b) Trout cardiomyocyte
Figure 6: Rat and trout cardiomyocytes before deconvolution algorithm was applied. On thisfigure upper images show cross-sections on XY plane and lower images correspond to cross-sections on XZ plane. Both displayed cross-sections are taken from the middle of image stacks.Subfigure (a) shows rat and (b) trout myocyte, where red color marks mitochondria and greensarcolemma. Note the size and cellular structure differences between rat and trout myocytes.
20
(a) 80th iteration (b) 200th iteration
Figure 7: Effect of deconvolution on rat myocyte. (a) shows the 80nd iteration (τk is minimal)that was considered as the best result of recorded image of rat cardiomyocyte, (b) is the 200thiteration. Cross-sections XY and XZ are displayed.
21
0 50 100 150 200Iteratsion
0.00004
0.00006
0.00008
0.00010
0.00012
0.00014
0.00016
0.00018
�k
Figure 8: The relative change between two successive iterations during deconvolution processof the rat cardiomyocyte. The τk reaches its minimum at 80th iteration, that corresponds tothe best result for image presented Figure 7(a).
22
Figure 9: Histograms of the rat myocyte, where blue corresponds to recorded image, redto 80th iteration and green to 200th iteration (images are shown on Figures 6(a), 7(a) and7(b), respectively). Background intensity and noise are reduced after deconvolution, the peakis shifted towards zero and has larger number of counts.
23
(a) 33rd iteration (b) 200th iteration
Figure 10: Effect of deconvolution on trout cardiac cell. (a) shows the 33rd iteration (τk isminimal) that was considered as the best of recorded image of trout cardiomyocyte, (b) is the200th iteration.
24
Figure 11: The relative change between two successive iterations during deconvolution processof the trout cardiomyocyte. The τk reach its minimum at 33rd iteration, that corresponds tothe best result for image presented in Figure 10(a).
0 50 100 150 200 250Intensity
100
101
102
103
104
105
106
107
Counts
recorded33th iteration200th iteration
Figure 12: Histograms of the trout myocyte from the deconvolution process. Blue correspondsto recorded image, shown on Figure 6(b), red corresponds to 33th iteration on Figure 10(a),and green corresponds to on Figure 10(b) before and after deconvolution, respectively.
25
4 Discussion
In this work, we tested both the standard Richardson-Lucy (2.2) and the Richardson-
Lucy with total variation regularization (2.3) algorithms on synthetic images (Figure 2)
as well as on experimental data (Figure 6). From the analysis of synthetic images we
found optimal range of parameter λ for TV regularization and, suitable threshold value
for stopping criterion. These optimal parameters were used for deconvolving confocal
microscope images.
The main finding of this work is that the RL algorithm with TV regularization with
carefully selected λ value gives satisfying result on recorded images. In addition, a new
software package has been developed that implements deconvolution algorithms used in
present work. The package will be released as an open source project that would allow
others to benefit from this work.
In the following, in this section we give an overview of earlier works, discuss on the
results that we got on the analysis of synthetic and confocal microscope images. We
give a brief summary of advantages and disadvantages of available software and give
arguments for the necessity of the open source software package that we plan to release
as a part of this work.
4.1 Comparison with Earlier Works
Dey et al., (2006) have shown that the result of deconvolution using the RL algorithm
with TV regularization depends on the value of the regularization parameter λ [Dey et al.,
2006]. Values of λ used in their work were similar to optimal values found in this work.
We used a measured PSF in contrast to [Dey et al., 2006] where it was computed using
microscope objective NA. In addition, we found optimal λ from a systematic analysis of
the deconvolution process on synthetic images.
4.2 Analysis of the Richardson-Lucy Algorithm with Total Vari-
ation Regularization
In this work we presented an analysis of the dependence of TV regularization parameter λ
on the deconvolution process. To estimate this parameter value we ran the deconvolution
algorithm several times with different values of λ. We compared each iteration with
original image through calculation of several norms (MSE and I) shown on Figure 4(a)
and 4(b). For each λ value, there exists an optimum number of iterations. For small
26
number of iteration, the results are similar for deconvolution with different λ values.
The noticeable difference between the results starts to occur around 10th iteration. If
deconvolved images are compared with several norms to the original image, we can
see that in some cases, the difference dependences on λ and iteration step goes closer
to original image in contrast with other cases. The minimum of the difference can
be reached at optimal λ. In this case, optimal λ is around 0.01, which leads to the
decrease of difference even after 200 iterations (see Figure 4). As result, an optimal
range 0.008 < λ < 0.03 for parameter λ was found from analysis on iteration history.
Recommended threshold for stopping criterion is around 10−5 to 10−4 [Dey et al.,
2006, Cannell et al., 2006, van Kempen et al., 1997], which is close to the values in
our study. But keep in mind that for larger images in size the threshold should also be
larger. When deconvolving images from confocal microscope, the relative change between
two estimates may be larger than given threshold, even when the deconvolution process
should be stopped. For this reason, the results of a deconvolution process on confocal
microscope images should be also verified visually.
4.3 Deconvolved Myocyte Images
We applied the deconvolution algorithm to rat and trout myocytes confocal images.
Deconvolution on these images was performed with λ = 0.01, as found form analysis
of synthetic images. From the results of deconvolution, we see remarkable improvement
of contrast and of noise reduction. All cell structures, labeled with fluorescent dyes,
became clearer. In agreement with previous studies, mitochondria are highly organized
in the rat cardiomyocyte (Figure 7) if compared to trout (Figure 10) [Birkedal et al.,
2006]. There is a noticeable difference in organization of t-tubules between rat and trout
cardiomyocytes. T-tubules (transverse tubule) are deep invaginations of the sarcolemma
allowing depolarization of the membrane to quickly penetrate to the interior of the cell.
In trout, t-tubules are absent, in contrast to rat. This is clear from the inspection of the
corresponding deconvolved images (Figure 7(a) and 10(a)).
As mentioned before, if the iterative deconvolution process is not stopped at the right
time, it can lead to creation of artifacts and noise dominant solutions. For example, in
Figure 10(b) we can see that new membrane-like structures start to appear and some
structures begin to disappear. Most likely this is caused by the low quality of recorded
image, where the noise and background level were quite high compared to the signal, as
in the case for trout images.
27
4.4 Available Software
There are both commercial and free deconvolution computer programs available for
microscopy image enhancement that use different methods and algorithms (some are
brought out in Appendix B).
Commercial image restoration software solutions are giving good results in image
enhancement and are easy to use, but, as a drawback, they are expensive and as a
rule, they do not support testing alternative deconvolution algorithms due to the closed
source development. There is a lack of open source software, that can compete with
commercial ones. One commonly used free image enhancement software package is BIG
- 3D Deconvolution. However, it is slow when deconvolving big data series. The Jansson-
van-Cittert algorithm used in the BIG - 3D Deconvolution plugin is one of simplest
iterative algorithms, where noise can lead to uncertainty as to whether the algorithm
actually achieves the optimal solution [Cannell et al., 2006]. We could not find the
details regarding about mentioned regularized deconvolution algorithm. So it is not clear,
which algorithm was used and what terms were applied for regularization. Therefore we
cannot compare it with our work. However, from our experience BIG - 3D Deconvolution
is not a viable alternative to commercial deconvolution programs. Because deconvolution
process with BIG - 3D Deconvolution takes lot of time and implemented algorithms to
not provide quality as commercial ones. Our open source software package allows one
to use Richardson-Lucy with total variation regularization algorithm that should lead to
well deconvolved images in practice.
Parallel to our work, open source Clarity Deconvolution Library was released in De-
cember 2008. Clarity Deconvolution has three different algorithms for image restoration
(see Appendix B). Clarity uses maximum likelihood iterative deconvolution, but in con-
trast to algorithm presented in this work, it is not regularized. In presence of NVIDIA
graphics card in computer, that supports CUDA (Compute Unified Device Architec-
ture), Clarity takes advantage of it and speeds up computations. In addition, Clarity is
included into freeware software called ImageSurfer (http://imagesurfer.org/) which
has reasonable graphical user interface and makes image restoration easy to use. All in
all, it can be said that Clarity Deconvolution Library seems to be one of the promising
choices open source for deconvolution.
The software package which will be released as a part of this work will be compli-
mentary to Clarity Deconvolution Library. In contrast to Clarity our package will have
commandline interface for batch jobs. That would simplify its use in clusters loading to
streamlining analysis of multiple images. We implemented the software in Python pro-
28
gramming language, that will allow to use a strength of Python language. Since Python
language has proven itself as an excellent fast-development platform and is simple to use,
the deconvolution software package should be a valuable platform for testing different
deconvolution algorithms.
4.5 Future Work
We plan to apply our software package to wide-field microscope images. This will require
measurement of PSF , analysis of deconvolution process on synthetic WF images, similar
to the analysis performed here on synthetic confocal images. After finding optimal λ value
and threshold for stopping criterion τ , we will deconvolve WF images of rat and trout
cardiomyocytes.
After that, we will prepare the software package for open source release. In open
source format, we plan to implement CUDA support and a cluster version.
29
5 Conclusion
As result of this work a software package for deconvolution using Richardson-Lucy algo-
rithm with total variation regularization was developed. The optimal range of regular-
ization parameter λ for our application were found. The performance of the algorithm
was tested on synthetic and confocal microscope images. Finally, it was shown that
prolonged iterations can lead to misleading results.
30
6 Personal Contribution
The author of this thesis carried out following tasks:
• developed a software that uses the RL algorithm with TV regularization
• analyzed the RL algorithm with TV regularization on synthetic images and found
the optimal parameter λ range and the threshold value for the stopping criterion
• recorded images of rat cardiomyocytes
• deconvolved and analyzed images of rat and trout myocytes
• wrote up the thesis.
31
A Derivation of the Richardson-Lucy Algorithm with
Total-Variation Regularization
In this appendix we derive the RL algorithm with TV regularization under the assump-
tion of a Poisson noise. As noted before, image formation can be mathematically de-
scribed as follows:
i(v) = P((o ⊗ h)(v)), (A.1)
where h = h(v) is the PSF that convolves acquired object o = o(v) to degraded image
i = i(v) with Poisson P , and v denotes voxel coordinates in three dimensional space.
In the following derivation we assume that voxel coordinates vary continuously, and the
result (2.3) will be obtained by the discretization of (A.11).
To derive the algorithm, we use Bayesian rule:
P (o|i) = P (o)P (i|o)P (i)
, (A.2)
where P (o) and P (i) are the prior probabilities of o and i, respectively; P (o|i) is the
likelihood probability of o, when i is given. In the case of image formation, the object
o can be statistically described by its prior probability function P (o). The likelihood
probability function P (i|o) is introduced to model blurred image i given the object o.
The RL algorithm is based on maximizing the probability function P (i|o) and estimating
the maximum value of the o given the known image i.
When noise follows the Poisson distribution, the likelihood probability is equal to the
product of individual probabilities because each pixel in a recorded image is statistically
independent from the others. Since at the detector the object is given as o⊗ h then the
equation for the likelihood probability can be written as follows:
P (i|o) =∏
v
(o ⊗ h)(v)i(v) e−(o⊗h)(v)
i(v)!(A.3)
Maximizing the equation (A.3) is equivalent to minimizing its negative logarithm. There-
fore, we take negative logarithm from both sides:
L(o) = − ln P (i|o ⊗ h) = − ln
(
∏
v
(o ⊗ h)(v)i(v) e−(o⊗h)(v)
i(v)!
)
=
= −ˆ
v
ln(o ⊗ h)i e−o⊗h
i!dv =
ˆ
v
(o ⊗ h − i ln(o ⊗ h) + ln i!) dv,
(A.4)
32
where L(o) is the likelihood estimation function of the object o.
The functional of total variation LTV (o) for any function o is defined as:
LTV (o) = λ
ˆ
v
|∇o|dv. (A.5)
where λ is regularization parameter, ∇o is gradient of o and |∇o| is the length of ∇o.
For deriving the Richardson-Lucy algorithm with total variation regularization we
have to minimize the L(o) + LTV (o).
First, let us find necessary conditions for L(o) minimum. We consider a small pertur-
bation ǫs of o, where ǫ ≪ 1, and s is test function. Notify that ln i! is constant relative
to o, then equation (A.4) reads (up to constant term):
L(o + ǫs) =
ˆ
v
(o + ǫs) ⊗ h − i ln [(o + ǫs) ⊗ h] dv =
=
ˆ
v
o ⊗ h + ǫ(s ⊗ h) − i ln [o ⊗ h + ǫ(s ⊗ h)] dv =
=
ˆ
v
o ⊗ h + ǫ(s ⊗ h) − i ln
[
o ⊗ h
(
1 + ǫs ⊗ h
o ⊗ h
)]
dv =
=
ˆ
v
o ⊗ h + ǫ(s ⊗ h) − i ln (o ⊗ h) − i ln
(
1 + ǫs ⊗ h
o ⊗ h
)
dv =
= L(o) +
ˆ
v
ǫ(s ⊗ h) − i ln
(
1 + ǫs ⊗ h
o ⊗ h
)
dv =
= L(o) +
ˆ
v
ǫ(s ⊗ h) − iǫs ⊗ h
o ⊗ hdv + O(ǫ2) =
= L(o) + ǫ
ˆ
v
(
1 − i
o ⊗ h
)
(s ⊗ h) dv + O(ǫ2) =
=
from the definition of convolution integral:´
vf s ⊗ h dv =
´
vs f ⊗ h dv,
where h is the adjoint of h
=
= L(o) + ǫ
ˆ
v
s
(
1 − i
o ⊗ h
)
⊗ h dv + O(ǫ2)
(A.6)
Second, we find necessary conditions for LTV minimum. Thus the equation (A.5) is
given as:
L(o + ǫs) = λ
ˆ
v
|∇(o + ǫs)| dv =
= λ
ˆ
v
√
|∇o|2 + 2ǫ∇o · ∇s + ǫ2|∇o|2 dv =
33
= λ
ˆ
v
|∇o|√
1 + 2ǫ∇o · ∇s
|∇o|2 dv + O(ǫ2) =
= λ
ˆ
v
|∇o| dv + ǫ
ˆ
v
∇o
|∇o| · ∇s dv + O(ǫ2) =
=
[
from the Gauss-Ostrogradsky theorem:´
v∇f · φ dv = −
´
vf div φ dv
]
=
= LTV − λǫ
ˆ
v
s div
( ∇o
|∇o|
)
dv + O(ǫ2)
(A.7)
In summary, the necessary conditions for L(o) + LTV (o) minimum is:
ˆ
v
s
[(
1 − i
o ⊗ h
)
⊗ h − λdiv
( ∇o
|∇o|
)]
dv = 0, (A.8)
that must hold for all test functions s. Therefore we have
1 − i
o ⊗ h⊗ h − λ div
( ∇o
|∇o|
)
= 0, (A.9)
where we have used the PSF property 1 ⊗ h = 1.
To obtain an iterative scheme for computing o from i and h, we define o(n) as the
result of nth iteration. Using the convergence condition o(n+1)
o(n) = 1 we write (A.9) as:
o(n+1)
o(n)
[
1 − λ div
( ∇o(n)
|∇o(n)|
)]
=i
o(n) ⊗ h⊗ h, (A.10)
from where we obtain the iterative Richardson-Lucy algorithm with total variation reg-
ularization assuming Poisson noise:
o(n+1) =
(
i
o(n) ⊗ h⊗ h
)
· o(n)
1 − λ div(
∇o(n)
|∇o(n)|
) . (A.11)
For λ = 0, the scheme (A.11) results to the well-known iterative Richardson-Lucy algo-
rithm:
o(n+1) =
(
i
o(n) ⊗ h⊗ h
)
o(n) (A.12)
34
B List of Available Deconvolution Software
• Huygens Deconvolution Software (http://www.svi.nl/) is one of the high-
quality commercial restoration and analysis programs for microscopy. This software
uses the following deconvolution algorithms:
– Classic Maximum Likelihood Estimation
– Iterative Maximum Likelihood Estimation
– Quick Maximum Likelihood Estimation
– Iterative Constrained Tikhonov-Miller
– Quick Tikhonov-Miller
– Costs 4 500 – 17 000 EUR; commercial license
• AutoDeblur (http://www.mediacy.com/) is Media Cybernetic’s deconvolution
software product which uses the following deconvolution algorithms:
– The 3D Inverse Filter
– The No/Nearest Neighbor
– DIC Restoration
– Non-Blind Deconvolution
– 3D – Adaptive-Blind Deconvolution
– 2D – Adaptive Blind Deconvolution
– Cost around 10 000 EUR; commercial license
• Volocity (http://www.improvision.com/products/volocity/) is Improvision
deconvolution software, which allows to use
– Fast Restoration
– Iterative Restoration: based on Maximum Entropy techniques
