Image Restoration Using Two-Dimensional Variations

In many applications (consumer and commercial imaging, medical imaging, robotics, space research, and etc.) observed images are often degraded due to atmospheric turbulence, relative motion between a scene and a camera, nonuniform illumination, wrong focus, etc. Image restoration refers to the problem of estimating the ideal image from its observed degraded one. Numerous restoration techniques (linear, nonlinear, deterministic, stochastic, etc.) optimized with respect to different were introduced (Banham & Katsaggelos, 1997 ; Jain, 1989; Sezan & Tekalp, 1990; Bovik, 2005; Gonzalez & Woods 2008). The amount of a priori information about degradation such as the size and shape of blurs, noise level determines how mathematically ill-posed the problem is. A priori information can be used in a variety of ways in modeling and algorithm development. The information about the nature of blur (e.g., linear or nonlinear and space-variant or space-invariant) and noise (additive or multiplicative) is used in modeling the input-output relation of imaging systems. In blur modeling, when the type of blur is known (e.g., out of focus, motion, turbulence), the blurring operator can be parameterized using only a few parameters. In image modeling, the ideal image can be modeled, for instance, on the basis of a priori Markovian assumption. In algorithm development, a priori information is used in defining constraints on the solution and in defining a criterion or a quantitative description of the solution. The blind and non-blind deconvolutions were extensively studied, and many techniques were proposed for their solution (Kundur & Hatzinakos, 1996; Bertero & Boccacci, 1998; Biemond et al., 1990; Sroubek & Flusser, 2003). They usually involve some regularization which assures various statistical properties of the image or constrains on the estimated image and restoration filter according to some assumptions. This regularization is required to guarantee a unique solution and stability against noise and some model discrepancies. One of the most popular fundamental techniques is a linear minimum mean square error method. It finds the linear estimate of the ideal image for which the mean square error between the estimate and the ideal image is minimal. This linear operator acting on the observed image to determine the estimate on the base of a priori second-order statistical information about the image and noise processes. For images with sharp changes of intensity, the appropriate regularization is based on variational integrals (Rudin, et al.,


Introduction
In many applications (consumer and commercial imaging, medical imaging, robotics, space research, and etc.) observed images are often degraded due to atmospheric turbulence, relative motion between a scene and a camera, nonuniform illumination, wrong focus, etc. Image restoration refers to the problem of estimating the ideal image from its observed degraded one. Numerous restoration techniques (linear, nonlinear, deterministic, stochastic, etc.) optimized with respect to different were introduced (Banham & Katsaggelos, 1997 ; Jain, 1989;Sezan & Tekalp, 1990;Bovik, 2005;Gonzalez & Woods 2008). The amount of a priori information about degradation such as the size and shape of blurs, noise level determines how mathematically ill-posed the problem is. A priori information can be used in a variety of ways in modeling and algorithm development. The information about the nature of blur (e.g., linear or nonlinear and space-variant or space-invariant) and noise (additive or multiplicative) is used in modeling the input-output relation of imaging systems. In blur modeling, when the type of blur is known (e.g., out of focus, motion, turbulence), the blurring operator can be parameterized using only a few parameters. In image modeling, the ideal image can be modeled, for instance, on the basis of a priori Markovian assumption. In algorithm development, a priori information is used in defining constraints on the solution and in defining a criterion or a quantitative description of the solution. The blind and non-blind deconvolutions were extensively studied, and many techniques were proposed for their solution (Kundur & Hatzinakos, 1996;Bertero & Boccacci, 1998;Biemond et al., 1990;Sroubek & Flusser, 2003). They usually involve some regularization which assures various statistical properties of the image or constrains on the estimated image and restoration filter according to some assumptions. This regularization is required to guarantee a unique solution and stability against noise and some model discrepancies. One of the most popular fundamental techniques is a linear minimum mean square error method. It finds the linear estimate of the ideal image for which the mean square error between the estimate and the ideal image is minimal. This linear operator acting on the observed image to determine the estimate on the base of a priori second-order statistical information about the image and noise processes. For images with sharp changes of intensity, the appropriate regularization is based on variational integrals (Rudin, et al., 1992;Perona & Malik, 1990;Chan & Wong, 1998). Minimization of the variational integrals preserves edges and fine details in the image. It is obvious that quality of the restored image depends on accuracy of mathematical model of image formation process. In particular, a good estimation of distortion parameters such as speed of camera movement, transparency of atmosphere or water, etc. is very important for restoration (Biemond et al., 1990).
Recently, restoration methods based on image variations were proposed (Milukova et al., 2010a(Milukova et al., , 2010b(Milukova et al., , 2011. In this chapter, image restoration with two-dimensional variations is presented. The restored image minimizes two-dimensional variations defined by Kronrod (Kronrod, 1950). We also consider the identification of distortion operator and estimation of its parameters. It is assumed that a monochrome stationary image distorted by homogeneous integrated transformation. Various physical problems can be modeled by such transformation. The spectral method of identification of distortion parameters uses only degraded image. Computer simulation results illustrate the performance of the proposed method for restoration of blurred images.

Variation concept for image restoration
Image restoration problem is usually formulated as follows. Undistorted (original) image z is recovered from the given equation: where :Z Q A  ( Z,Q are metric spaces) is linear or nonlinear operator, Ζ z  , n is noise, v is observed distorted image. A general approach for image restoration can be formulated using statistical estimation methods and the theory of solving of ill-posed problems (Tikhonov & Arsenin, 1977). The restoration problem is a typical inverse problem of mathematical physics, and, therefore, it can be correctly solved on the base of mathematical methods. The restored image can be obtained by minimization of the following functional: where Q  is a metric in Q . Note that various definitions of a distance Q  between two images may be used. It is easy to show that the solution of the optimization problem in Eq.
(2) is not unique even when the operator A and the distorted image v are exactly known, and no additive noise. A priori information about the original image should be used to obtain a unique and stable solution from the set of solutions. The simplest way to guarantee uniqueness and stability of the solution is to describe the image model with a functional () z  that possesses stabilizing properties. In this case the image restoration problem can be reduced to conditional or unconditional optimization problem, in particular to the Tikhonov's minimization (Tikhonov & Arsenin, 1977), where  is a parameter of regularization. Note that the statistical methods used in image restoration lead to optimization problems, which are similar to that of Eq. (3). For instance, using Bayes' strategy (Kay, 1993) we can obtain the optimal estimation in the following form: It is commonly assumed that the original image is a smooth function with respect to the Sobolev space (Adams, 1975), and the stabilization functional in Eq. (5) is () p q q W zz  . Quadratic forms can be used in order to avoid nonlinear restoration algorithms. Note that a Gaussian image model leads to minimization of a quadratic form. In discrete case it corresponds to the Sobolev norm for 2 p  in Eq. (5). On the other hand, the use of quadratic forms in image restoration often yields undesirable results because of real images are not Gaussian. Now suppose that an image to be restored is a function of bounded variations. Therefore, it may be written as The variation of a 1D function () , [,] fx x ab  is defined as follows: It can be shown, that if the image (,) zxy , (,) xy D  consists of 1D functions of bounded variation along its rows and columns then the image is also a 2D function of bounded variation. Different multidimensional variations were proposed such as variations of Arzela, Vitali, Tonelly, etc. (Vitushkin, 1955). A different approach was suggested by Kronrod, who introduced two functionals in order to describe an image as a function of two variables (Kronrod, 1950). The functionals are given as follows: where t e is t -level set of the function (,) zxy , i.e. a set of points (,) xy with function values equal to t, 0 () t me is the number of components of t e , and 1 () t me is the length of the set t e .
The class of functions of bounded variations given in Eq. (8) is very extensive. These functions possess the following attractive properties: they are differentiable almost everywhere and their Fourier series are convergent almost everywhere. Note that numerous attempts to create a mathematical image model with the help of one functional were unsatisfactory. It can be done on the base of two (independent in a certain way) functionals. It is interesting to point out, that the first variation d 1 in Eq. (8) is a topological characteristic of the image. If the original image is a continuous differentiable function, then the second variation can be represented as If only the second variation is used, the image restoration can be carried out as follows (Perona & Malik,1990): where grad(.) is a gradient operator.
It is of interest to note that this nonlinear method of image restoration minimizes the functional that is identical to the Kronrod's second variation. We propose to minimize the functional in Eq. (10) subject to constraint on the first Kronrod's variation of the image. This approach is referred to as conditional variation approach . Next with the help of computer simulation we illustrate the difference in the performance of two variation methods: minimization of the functional in Eq. (10) and conditional minimization of the same functional. Additionally, the performance of minimum norm image restoration from Eq. (10) without considering variations is also provided.

Restoration of uniformly blurred image with spatial variations
The impulse response of the 1D uniform blur can be expressed as follows: where L determines a blur extension. It is known that point spread functions for such blurs do have zeros in the frequency domain, and they can be uniquely identified by the location of these zero crossings (Cannon, 1976;Gennery, 1973 www.intechopen.com and the linear system of equations in Eq. (1) is rewritten as where   where m e is a basis of kernel A,   m  are real variables, is a particular solution of linear system of equations in Eq. (13). Suppose that 0 L z  , the particular solution can be expressed as The basis m e can be found from Eq. (13)  .
It is of interest to note that the basis of kernel A contains vertical columns of unities; therefore, a general solution of Eq. (13) could contain a periodic structure with the period of blurring (Buades et al., 2006). In our computer simulation we compared three methods: 1) first, substituting Eq. (15) into Eq. (10) and minimizing the functional with respect to   m  ; 2) the second method minimizes the same functional subject to the first Kronrod's variation given in Eq. (8), which is taken close to that of the original image; 3) minimum norm image restoration from Eq. (10) without considering variations. Actually, if the inverse of the blur operator exists, it can be applied to the observed image to obtain an estimate. This is called inverse filtering. The estimate differs from the actual image by the additional error of amplified noise, and depending on the nature of the blur operator and the noise, it may drastically obscure the desired image information. Hence, inverse filtering is extremely noise sensitive. If the inverse operator does not exist, a solution can be found on the basis of a least squares criterion. A least squares solution minimizes the norm of the residual signal Az-v. The least squares solution with minimum norm (energy) is called also the pseudoinverse filtering (Jain, 1989). The first tested method is referred to as Grad method, the second one is called Grad-conditional method, and the last method is named Min-norm method. Fig. 1(a) shows a test input image used in experiments. The size of the image is 256×256 pixels, N=256. The signal range is . The test input scene is homogeneously blurred with L=7. The blurred image is shown in Fig. 1(b).

Figs. 3(a), 3(b), and 3(c) show differences between the original image and that of restored
with the Grad algorithm, the Grad-conditional algorithm, and the Min-norm algorithm, respectively. We see that the second algorithm, which takes into account two Kronrod's variations yields the best recognition performance. A quantitative comparison is given by the peak signal-to-noise ratio (PSNR), The image restored with the Grad method has d1=2050, d2=12015, the image restored with Grad-conditional method has d1=2103, d2=12120, and finally, the image restored by the Min-norm method possesses d1=2218, d2=12343. So, in order to achieve a good restoration, it is important take into account topological characteristics of the image to be restored. These characteristics can be described by the first Kronrod's variation.

Identification of degradation operators and estimation parameters in the Fourier domain
In this section, identification of operator degradation type and estimation of distortion parameters is discussed. Monochrome stationary image distorted by homogeneous linear transformation is considered. Various physical problems can be modeled by such degradations (Bertero & Boccacci, 1998). Identification of the distortion operator is carried out using the Fourier spectrum of the distorted image. Automatic image restoration is performed in three steps, that is, i) identification of distortion operator, ii) estimation of distortion parameters, and iii) image restoration with estimated parameters. Certain types of distortion operators are completely characterized by attributes such as the location of frequency-domain zeros. The techniques (Cannon, 1976;Gennery, 1973) make the following two assumptions: (i) the blurring produces zero crossings in the frequency domain and it can be completely characterized by the location of these zero crossings, and (ii) the location of zero crossings can be determined from the Fourier transform or power cepstrum (the logarithm of the power spectrum) of the observed image. These methods are very simple to use and they can successfully applied in many real-life situations. It is indeed true that the models for motion and focus blurs do have zeros in the frequency domain, and they can be uniquely identified by the location of these zero crossings. On the other hand, blurring models that do not have zero crossings in the frequency domain (e.g., Gaussian modeling atmospheric turbulence) cannot be identified by these techniques. Furthermore, the identification of the zero crossings from the observed image may be quite difficult due to the presence of strong observation noise. Almost all practical implementations of the restoration algorithms assume that the observation noise is a zeromean, white Gaussian process that is uncorrelated to the image signal. In this case, the noise field is completely characterized by its variance, which is commonly estimated by the sample variance computed over a low-contrast local region of the observed image (Yaroslavsky & Eden, 1996).
Let us consider an observed image degraded with a linear spatially invariant system and additive sensor noise, that is, ,s i n c ( ) where L is the motion path, is the Dirac delta function.
Atmospheric turbulence is a common blur in remote sensing and aerial imaging. For long term exposure through the atmosphere Gaussian model is used. So, the impulse response and the frequency response of the linear system of turbulence blur are given, respectively, as where α is the parameter that determines the severity of blur.
Defocusing is another common type of blurring owing to the finite size of the camera aperture. When the defocusing blur is large, the following uniform model is used. The impulse response and the frequency response of the linear system can be expressed, respectively, as where J 1 is the first-order Bessel function.
Image blurring also occurs in image acquisition by scanners in which the image pixels are integrated over the scanning aperture. Example of such degradations can be found in image capturing by radar, beam-forming arrays, and display systems using television raster. The impulse response and the frequency response of the linear system can be written, respectively, as where axb is the rectangular aperture.
It is of interest to compare the spectra of the original and degraded images. From numerous experiments it is known, that the spectral magnitude of a realistic image is not very informative. It contains information about distribution of the signal energy in the frequency domain. For instance, if we exchange the spectral magnitudes of two similar images belonging to the same class and perform the inverse Fourier transform, then the difference in visual appearance of the original and obtained images will be negligible (Yaroslavsky & Eden, 1996). However, the difference between the spectral magnitudes of the original and degraded images may be significant due to the spectral magnitude of the frequency response of the linear system. Figs In order to identify the distortion operator a database containing various images of spectral magnitudes was created. Training elements of database were obtained on the base of mathematical modeling or computer simulation. In practice, the number of degradation operators is not very large. Next, the spectral magnitude of a degraded image is matched to those of the database. This simple recognition system works well to identify the type of degradation operator for common blurs. Figs. 5 illustrate spectral magnitudes for different common blurring operators. One can observe that spectral magnitudes of distorted images contain mainly the information about distortion operators such as zero crossings on the Actually, composite degradations can be considered as a combination of the basic distortion operators. In this case, the Fourier spectrum of a new composite operator is the product of the spectra of used basic operators. Figs. 6(a) and 6(b) show the spectral magnitudes of the test original image degraded with isotropic blur and horizontal motion. www.intechopen.com One can observe that since one of the distortion operators dominates, the spectral magnitudes of composite degradations cannot be synthesized by a simple combination of those of the basic degradation operators. This means that the number of training elements of a matching system should be drastically increased.
In the Fourier representation of images, spectral magnitude and phase tend to play different roles (Oppenheim & Lim, 1981). For instance, in some situations many of the important features of an image are preserved if only the phase is retained. Furthermore, under a variety of conditions, phase information alone is sufficient to completely reconstruct an image to within a scale factor. If noise fluctuation in Eq. (19) is small, the phase of the distortion operator is equal to the difference between the phases of the degraded and original images. If the distortion operator is Gaussian, its phase is zero, and the phases of the distorted and original images coincide If the distortion operator is a finite function, e.g. (,) xy W  , then the phase of the distorted image may differ from the phase of the original image by ±π, and points at which the phase jumps by ±π coincide with the location of zeros of the spectral magnitude of the degraded image. These zeros, for even functions are all located on the real axis. So, the phase of the original image either coincides with that of the distorted image or differs from that of the distorted image by ±π. Fig. 7 shows the differences between the phases of the distorted and original images for rectangular and circular aperture blurs. Therefore, under certain conditions, we can identify the type of the distortion operator and estimate its spectral phase from the observed degraded image.

Conclusion
In this chapter we treated the problem of restoring linearly degraded image using twodimensional image variations. The restored image minimizes the objective functional subject to the Kronrod's variations. In order to achieve a good restoration, it is important take into account topological characteristics of the original image, which are well described by the first Kronrod's variation. The first step in restoring a degraded image is the identification of the type of degradation operator. It can be done by matching of the spectral magnitude of the degraded image with those of created database. Under certain conditions, the phase of the distortion operator may be also estimated from the distorted image. Extensive testing it was shown that the original image can be automatically restored by proper choice of the parameters of the proposed method.

Acknowledgment
This work was supported by the Russian Foundation for Basic Research, project no. 11_07_00361.