G.M.P. van Kempen and L.J. van Vliet, Improving the restoration of textured objects with prefiltering, in: H.E. Bal, H.
Improving the restoration of textured objects
with prefiltering
Geert M.P. van Kempen, Lucas J. van Vliet
Pattern Recognition Group, Faculty of Applied Physics,
Delft University of Technology,
Lorentzweg 1, 2628 CJ Delft, The Netherlands
{geert, L.J.vanVliet}@ph.tn.tudelft.nl
keywords: Image restoration, Gaussian filtering, median filtering, prefiltering
Abstract
The restoration of images acquired by a fluorescence microscope is in the presence of noise known to be a difficult problem. This presence of noise can result in restored images that contain severe noise artifacts, or that are smoothed by a strongly imposed regularization. In both cases, small textures present in the image, are lost in the restoration procedure.
In this paper we propose a prefiltering method to reduce the noise in the image without hampering the object. We have applied a Gaussian filter and a median filter prior to the restoration, and show how we can compensate for the extra blurring of the Gaussian in the restoration procedure. We show how this prefiltering improves the performance of the restoration algorithms on textured objects. The experiments were performed on textured disks convolved with the point spread function of a fluorescence microscope and distorted with Poisson noise.
1. Introduction
The restoration of images acquired by a fluorescence microscope is, in the presence of noise, known to be a difficult problem. In previous work, we have investigated the performance of the Richardson-Lucy and iterative constrained Tikhonov-Miller (ICTM) algorithm. We found that the performance of the two algorithms is strongly dependent on the signal-to-noise ratio in the image. Furthermore, we found that the results of the Richardson-Lucy algorithm are often hampered by noise artifacts, whereas the imposed regularization of the ICTM algorithm produces smoothened results. Therefore, textures present in the original image could be lost during restoration.
We have therefore investigated filtering techniques that decrease the amount of noise by suppressing those parts of the image spectrum that do not contain any signal information (or where the noise contribution is far larger
than the signal contribution). These frequencies will prevent signal recovering and only amplify noise in the final result.
We propose two methods to suppress these (high frequency) parts of the spectrum: by convoluting the recorded image with a Gaussian, and by filtering it with a median filter. Although, both the median and Gaussian filter will mainly suppress high frequencies, lower (object) frequencies are also effected. However, being a linear filter, we can compensate for the extra blurring imposed by the Gaussian filter on the image by convolving the PSF with the same Gaussian.
We have tested these pre-filtering techniques for both the Richardson-Lucy algorithm and ICTM algorithm in a simulation experiment. In this experiment we have convolved a textured disk with a computed point spread function (PSF) and tested the performance of these algorithms as function of the signal-to-noise ratio. The texture has been made by modulating the object with a sine in all dimensions. The performances of the restoration algorithms has been measured with the mean-square-error and I-Divergence measures.
2. Image restoration
2.1 Image acquisition
The incoherent nature of emitted fluorescence light allows us to model the image formation of the fluorescence microscope as a convolution of the object function f(x) with the point spread function (PSF) of the microscope h(x), with x being a two-dimensional coordinate in the space X. The image g(x) formed by an ideal noise-free fluorescence microscope can thus be written as
g x =h x ⊗f x (1)
with ⊗ denoting the convolution operator. However, noise induced by photon-counting (Poisson noise), by
the readout of the camera (Gaussian), and by the analog-to-digital conversion (uniform), disturbs the image g(x). Furthermore, it is common in fluorescence microscopy to measure a non-zero background level, arising from auto-fluorescence, inadequate removal of fluorescent staining material, dark current, spurious light sources or electronic sources. We model the noise distortion and background here in a general form,
m x =N g x +b x (2)
with m(x) being the recorded image, b(x) the background and N(.) the noise distortion function.
2.2 The Richardson-Lucy and ICTM algorithms Generally, restoration methods yield an estimate of the original image ( )f x given an imaging model, a noise model, and additional criteria. These criteria depend on the imposed regularization, and constraints implied on the solution found by the restoration algorithm.Although the methods we investigate in this paper share the imaging model, they differ significantly due to the different modeling of noise distortion on the image, imposed constraints and regularization.
The performance of state-of-the-art scientific cameras is limited by photon counting noise. The image formation of the microscope is than described as a translated Poisson process1. The log likelihood function of such a Poisson process, is given by
L f g x dx g x b x m x dx
X X
= − + ln + (3)
where we have dropped all terms that are not dependent on f(x). The maximum likelihood estimator (MLE) ( )f x for restoring f(x) given h(x) and m(x) can be found using the Expectation-Maximization (EM) algorithm, which results in the following iterative formula,
f x f x h x h x f x b x m x k k k + = − ⊗ + ⊗ 1 (4)
This iterative algorithm, often referred to as EM-MLE, (4) is identical to the Richardson-Lucy algorithm2 and was first adopted to fluorescence microscopy by Holmes3. The EM algorithm ensures a non-negative solution when a non-negative initial estimate
( )
f0 x is used.
The convergence of the Richardson-Lucy algorithm has been observed to be slow, we have therefore used an acceleration method based on the idea introduced to microscopy by Holmes3. Rewriting equation (4) as
fkacc+1 = fk+α∆fk+1 with f∆k+1= fk+1−fk (5) an α≥1 has to be searched for that maximizes (3) under the condition that fkacc+1 remains non-negative. This condition restricts the maximum value of α to be searched for to αmax min α = = − < ∀x + + k k k x f x f x f x ∆ 1 ∆ 1 0 (6)
The optimum value for α in the interval [1, αmax] can be found with a small number of Newton iterations. We have used only one iteration to find α, in which case the Newton iteration is equal to a Taylor expansion.
The iterative constrained Tikhonov-Miller (ICTM) algorithm is based on the assumption that the noise distortion function N(.) can be modeled as an additive noise function,
′ = − = ⊗ +
m x m x b x h x f x n (7)
with n the additive noise distortion with zero mean. Finding an estimate ( )f x from (7) is known as an ill-posed problem4. To solve this ill-posedness, Tikhonov
defined the regularized solution ( )f x of (7) as the one that minimizes the well-known Tikhonov functional4
Φ f = m x′ −h x ⊗f x 2+λ f x 2 (8) with λ the regularization parameter and
f x f x dx
X
2 2
=
the integrated square modules of f(x). Since the intensity of an imaged object represents signal energy, it is positive. To insure that the solution f( )x is non-negative, the Tikhonov functional has to be minimized in an iterative manner5. The ICTM algorithm finds the minimum of (8) by using conjugate gradients. The so-called conjugate gradient direction of (8) is given by
p x r x p x r x r x k k k k k k k = + − = − γ 1 γ 2 1 2 , (9)
with rk(x) denoting the steepest descent direction,
r xk f f x h x m x g x f x = − ∇1 = ⊗ ′ − − 2 Φ λ
with g x defined as ( ) f x convolved with the PSF h(x). A new conjugate gradient estimate is now found with
fk+1 x = f x +βkp xk (10)
After each iteration the ICTM algorithm imposes its non negativity constraint by setting all negative intensities to zero. Therefore, the optimal ßk can not be calculated analytically but must be searched for iteratively. In our implementation, we have used three Newton iterations to find the step size ßk, which is a faster method than the golden section rule line-search we used previously6. 2.3 Prefiltering
In previous experiments6, we observed that the performance of the Richardson-Lucy and the ICTM algorithms are strongly dependent on the signal-to-noise ratio (SNR) of the acquired image. Furthermore, these
the Richardson-Lucy algorithm is very sensitive to noise, whereas the ICTM algorithm produced a relative smooth solution, induced by a strong regularization. The obvious way to improve the result of image restoration is therefore to increase the signal-to-noise ratio of the image. Given that the noise induced by photon counting is the dominant source of noise, this can be achieved by collection more photons per pixel. However, this is not always possible, limits on the total acquisition time, saturation of the fluorescence molecules, bleaching, and the capacity of a CCD camera’s potential well can limit the SNR of the acquired image.
PSF and noise Gaussian
PSF convolved with Gaussian
Object Measured Image Filtered Image Image Acquisition Pre-filtering
Image Resoration
Figure 1. The incorporation of the proposed Gaussian pre-filtering in an image restoration procedure
We have therefore investigated three methods that reduce this noisy sensitivity after the image has been acquired by suppressing those parts of the image spectrum that do not contain any signal information (or where the noise contribution is far larger than the signal contribution). These frequencies will prevent signal recovering and only amplify noise in the final result. We have suppressed these (high frequency) parts of the spectrum by convolving the acquired image with a Gaussian, by filtering it with a median filter, and by a combination of these two methods. Although, both the median and Gaussian filter will mainly suppress high frequencies, lower (object) frequencies are also effected. However, being a linear filter, we can compensate for the extra blurring of the Gaussian filter by convolving the PSF with the same Gaussian, as shown in Figure 1.
3. Experiments & results
We have performed a simulation experiment which tests the performance of the Richardson-Lucy and ICTM algorithm on restoring a textured disk convolved with a PSF and distorted with Poisson noise. The performance of the algorithms is measured with the mean-square-error and I-Divergence as a function of the signal-to-noise ratio.
3.1 Object and noise generation
We generated the disks using an analytical description of their Fourier transform, as given byvan Vliet7,
S J r r disk ω x y ω ω ω ω ω = 2 1 = 2+ 2 (11)
with r the radius of the disk. The Fourier transform of the disk is multiplied by a Gaussian transfer function (sigma of 1 pixel in the spatial domain) to ensure bandlimitation. Generated in this way, the disks are free from aliasing effects that arise from sampling non-bandlimited analytical objects below the Nyquist rate. We have added texture to the disk by modulating the disk in both dimensions with a sine. We have used sine of with a period equal to 1/7th of the image size, and a modulation depth of 20% of the image intensity.
The point spread function of a fluorescence microscope can be generated using an analytical formula for the optical transfer function (OTF), the Fourier transform of the PSF, OTF c c c ( )ω arccos π ω ω ω ω ω ω = 2 − 1− 2 (12) with ω π λ c NA = 4
and NA the numerical aperture and λ the emission wave length. For our simulations, we have selected the microscope parameters corresponding to typical working conditions: a numerical aperture of 1.3 and an emission wavelength of 530.
The performance of the restoration algorithms is measured as a function of the signal-to-noise ratio which we define as
SNR= E
ε (13)
with ε the total power of the noise and E the total power of the convolved object. The simulated images are distorted by Poisson noise, which is generated by using the intensity of the convolved disk as an average of a spatially variant Poisson process. We have varied the signal-to-noise ratio of the simulated images by changing the photon-conversion efficiency.
3.2 Performance measures
We have used the mean-square-error (MSE) and the I-Divergence as performance measures. The MSE measures the difference in energy between the two compared signals and is defined as
MSE f f, f x f x
= − 2 (14)
The I-Divergence, introduced by Csiszár8,
I a b a a
b a b
, = ln − −
measures the distance of a function b to a function a. He has postulated a set of axioms of regularity (consistency, distinctness, and continuity) and locality that a distance measure should posses. He concluded that for functions that are required to be non-negative, the I-Divergence is the only consistent distance measure. Snyder9 has shown
that maximizing the mean of the log likelihood of (3) is equal to minimizing Csiszár's I-Divergence,
I f f L f E L f g x g x g x g x g x dx X , ln = − = − + (15)
3.3 Iterative optimization: where to start and when to stop
The two investigated restoration algorithms need a first guess to start their iterations. We have used the measured image m(x) as a first estimate f x to start the( ) ICTM algorithm, and a image with a constant intensity equal to the average intensity of m(x) for the Richardson-Lucy algorithm. The performance of the algorithms is highly dependent on the number of iterations. We have allowed the ICTM to converge, by stopping the algorithm when the change of the functional (fk+1−fk) fk was smaller than a preset threshold value. We have used a threshold value of 10-6. We have used a similar approach for the Richardson-Lucy algorithm in previous research6. However, we found that method for stopping the Ricardson-Lucy algorithm did not always produce solutions that minimized the I-0.001 0.01 0.1 I-D iv e rg e n c e 100 1000 signal-to-noise ratio R-L median Gauss R-L median R-L Gauss Richardson-Lucy 0.001 0.01 0.1 I-D iv e rg e n c e 100 1000 signal-to-noise ratio
ICTM median Gauss ICTM median ICTM Gauss ICTM
Figure 2. The I-Divergence performance of the Richardson-Lucy algorithm (left) and the ICTM algorithm (right) with and without Gaussian, median or combined median-Gaussian prefiltering.
10 100 m e a n s q u a re e rr o r 100 1000 signal-to-noise ratio R-L median Gauss R-L median R-L Gauss Richardson-Lucy 10 100 m e a n s q u a re e rr o r 100 1000 signal-to-noise ratio
ICTM median Gauss ICTM median ICTM Gauss ICTM
Figure 3. The mean-square-error performance of the Richardson-Lucy algorithm (left) and the ICTM algorithm (right) with and without Gaussian, median or combined median-Gaussian prefiltering.
Divergence between f(x) and f x . (The Richardson-( ) Lucy algorithm minimizes the I-Divergence between m(x) and f x convolved with h(x). Under noisy( ) circumstances, adaptation to noise can decrease the latter measure, whereas it will increase the former). To allow for a fair comparison, we have therefore stopped the Richardson-Lucy algorithm when the minimum of the I-Divergence between f(x) and f x was reached. This is( ) method of stopping the iteration is of course only possible in a simulation experiment, not in an experiment with real data. Therefore the question of how to stop the Richardson-Lucy algorithm needs to be investigated further. We the used the method of constrained least squares, as described by Galatsanos10,
to determine the value for the regularization parameter λ.
3.4 Restoration of textured disks
This experiment tests the influence of the proposed pre-filtering methods on the performance of the Richardson-Lucy and ICTM algorithm to restore textured disks convolved with a PSF and distorted by Poisson noise. The performance of the algorithms is measured with the MSE and I-Divergence measure as a function of the SNR. We have generated disks with a radius of 1.0 µm, an object intensity is 200.0 and a background of 20 (Figure 5). The images 256 x256 pixels in size, the SNR ranges from 16 to 1024 (12 dB - 30 dB). We have restored the textured disks with using the two restoration methods with and without using three methods of pre-filtering: Gaussian pre-filtering (with a sigma of 2.0 pixels), median prefiltering (with a filter size of 5
pixels), and a combined median-Gaussian filtering (with a sigma of 1.4 pixels, and a filter size of 3 pixels). Figure 2 shows the I-Divergence performance of the Richardson-Lucy and ICTM algorithm as function of the signal-to-noise ratio. The mean-square-error performance of both algorithms are shown in Figure 3. The graphs clearly show the strong dependence of the performance of these algorithms on the SNR. Furthermore it shows the improvement of the three pre-filtering methods on both algorithms. The restoration results, for a image with a SNR of 16.0, are shown in Figure 4.
Figure 5. The textured disk (left) and its image (right) Cross-sections of the center line of the images shown in Figure 4 together with that of the textured disk are shown in Figure 6.
4. Conclusions
We have tested the performance of the Richardson-Lucy and ICTM restoration algorithms on simulated fluorescence microscope images. Both algorithms
Richardson-Lucy R-L median R-L Gauss R-L median Gauss
ICTM median Gauss ICTM Gauss
ICTM median ICTM
Figure 4 The results of Richardson-Lucy (top) and the ICTM algorithm (bottom) are shown without prefiltering and with median, Gaussian, and combined median-Gaussian prefiltering.
reduce the distortions of these images. However we found a strong dependence of their performance on the signal-to-noise ratio of the images. We propose a prefiltering method to improve the restoration performance of images which SNR can not be improved during acquisition. This method reduces the noise by suppressing those parts of the image spectrum that do not contain any signal information or which are dominated by noise. We implemented this in three ways: by using a median filter, a Gaussian filter and a combination of the two. The extra blurring of the Gaussian filter can be compensated for by filtering the PSF with the same Gaussian. We have shown that performance, measured by both the mean-square-error as the I-Divergence measure, is greatly improved for both algorithms by using pre-filtering.
5. Acknowledgment
This work was partially supported by the Royal Dutch Academy of Sciences (L. v V.). We would like to thank Peter Verveer for many discussions on image restoration and for contributing a considerable part of the software used in this research.
6. References
[1] D. L. Snyder and M. I. Miller, Random Point Processes in Time and Space, Springer Verlag, Berlin, 1991.
[2] W. H. Richardson, Bayesian-based iterative method of image restoration, J. Opt. Soc. of Am., 62, pp. 55-59, 1972.
[3] T. J. Holmes and Y. H. Liu, Acceleration of maximum-likelihood image restoration for fluorescence microscopy and other noncoherent imagery, J. Opt. Soc. of Am. A, 8, pp. 893-907, 1991.
[4] A. N. Tikhonov and V. Y. Arsenin, Solutions of Ill-Posed Problems, Wiley, New York, 1977. [5] R. L. Lagendijk and J. Biemond, Iterative
identification and restoration of images, Kluwer Academic Publishers, Boston/Dordrecht/London, 1991.
[6] G. M. P. v. Kempen, L. J. v. Vliet, P. J. Verveer, and H. T. M. v. d. Voort, A Quantitative Comparison of Image Restoration methods for confocal microscopy, J. Micros., 185, 1997. [7] L. J. v. Vliet, Grey-scale Measurements in
Multi-Dimensional Digitized Images, Delft University Press, Delft, 1993.
[8] I. Csisár, Why Least Squares and maximum entropy? An axiomatic approach to inference for linear inverse problems, The Annals of Statistics, 19, pp. 2032-2066, 1991.
[9] D. L. Snyder, T. J. Schutz, and J. A. O'Sullivan, Deblurring subject to Nonnegative Constraints, IEEE Trans. Signal Processing, 40, pp. 1143-1150, 1992.
[10] N. P. Galatsanos and A. K. Katsaggelos, Methods for choosing the regularization parameter and estimating the noise variance in image restoration and their relation, IEEE Trans. Image Processing, 1, pp. 322-336, 1992. 0 50 100 150 200 250 0 50 100 150 200 250 300 0 50 100 150 200 250 0 50 100 150 200 250 300 0 50 100 150 200 250 0 50 100 150 200 250 300 0 50 100 150 200 250 0 50 100 150 200 250 0 50 100 150 200 250 0 50 100 150 200 250 300 0 50 100 150 200 250 0 50 100 150 200 250 300 0 50 100 150 200 250 0 50 100 150 200 250 0 50 100 150 200 250 0 50 100 150 200 250
Richardson-Lucy R-L Gauss R-L median R-L median Gauss
ICTM ICTM Gauss ICTM median ICTM median Gauss
in te n s it y in te n s it y
Figure 6. Cross sections of the center line of the textured disk together with the result of the Richardson-Lucy algorithm (left) and the ICTM algorithm (right) with and without Gaussian, median or combined median-Gaussian prefiltering.
