INTRODUCTION

---

By image restoration, we seek to recover the original sharp image using a mathematical model of the blurring process. The key issue is that some information on the lost details is indeed present in the blurred image, but this information is hidden and can only be recovered if we know the details of the blurring process. Due to various unavoidable errors in the recorded image, we can not recover the original image exactly. The most important errors are fluctuations in the recording process and approximation errors when representing the image with a limited number of digits (Hansen and James, 2006). We can broadly classify restoration techniques into two classes: the filtering reconstruction techniques and the algebraic techniques. The filtering techniques are rather classical and they make use of the fact that noise signals usually have higher frequencies than image signals. This means that image signals die out faster than noise signals in high frequencies. Examples of the restoration filters are the deconvolution filter, in which the transfer function of the degraded system is inverted to produce a restored image and the Wiener filter that uses the mean-squared error criterion to minimize the error signal between the original and degraded images. The MSE technique treats all errors equally, regardless of their spatial location in the image. A limitation in this filter is that it cannot handle dynamically changing image and noise signals. Khireddine et al. (2007) used this filter in 2D case for digital image restoration. The lucy-richardson algorithm can be used effectively when the point-spread function PSF (blurring operator) is known. The blind de-convolution algorithm maximizes the likelihood that the resulting image, when convolved with the resulting PSF, is an instance of the blurred image, assuming Poisson noise statistics.

The linear algebraic techniques utilize matrix algebra and discrete mathematics for solving the problem of image restoration. Some of the algebraic restoration techniques are: the unconstrained and constrained reconstruction techniques. Regularized de-convolution can be used effectively when constraints are applied on the recovered image (e.g., smoothness) and limited information is known about the additive noise. Optimization methods can be used to solve a large scale constrained linear least-squares optimization problems to recover images. Lamotte and Alt (1994) did a comparative study of four optimization algorithms based on simulated annealing: the Gibbs sampler, the metropolis algorithm, the iterated conditional modes and an original method of random descent. Chen et al. (1999) modeled the restoration problem as an optimization problem and used the genetic algorithm for gray images restoration. Stack filters that are non-linear spatial operators used for noise suppression have been formulated by Undrill and Delibassis (1997) as an optimization problem solved by genetic algorithm for restoring magnetic resonance images corrupted with uncorrelated additive noise. A new method for the restoration of images degraded by noise and spatially invariant blur has been proposed by Landi (2007) in which image restoration problem is replaced by an equality constrained minimization problem. A quasi-Newton method is applied to the first-order optimality conditions. In each quasi-Newton iteration, the hessian of the lagrangian is approximated by a circulant matrix and the fast fourier transform is used to compute the quasi-Newton step. The quasi-Newton iteration is terminated according to the discrepancy principle. A novel method called edge-preserving regularization is presented by Gu and Gao (2008). This method is used to solve an optimization problem whose objective function has a data fidelity term and a regularization term; the two terms are balanced by a parameter λ.

Most of the optimal techniques that have been proposed in literature over the past few decades to solve such problem by iterative optimization procedures are computationally demanding and time consuming. The novel approach (PSO) introduced in this study is to take advantage of swarm intelligence in order to facilitate optimization process in total variation regularized methods.

NEUTRON RADIOGRAPHY

Digital radiological image is a digital image acquired by a certain radiological procedure which can be: X-rays, neutron radiography, gamma camera or a nuclear magnetic resonance image. It is a two-dimensional MXN array of non-negative integers (gray levels). For neutron radiography, Fig. 1, the gray level value represents the relative linear neutron attenuation coefficient of the object. Each of these gray images has 8-bit representations of their intensity levels. Hence, there are 256 gray levels. Degradations in this imaging technique are essentially due to bad situation with respect to randomly distributed neutron flux causing dissimilarity in images taken for the same object, in addition to the presence of gamma radiations causing additive noises (Wong et al., 1995).

Fig. 1:	Neutron imaging system (Kharfi et al., 2005)

Fig. 2:	Simplified model for image degradation/restoration process. The image signal f(x, y) is subjected to a linear degrading function H(x, y) and an arbitrary noise η(x, y) is added to produce the degraded signal g(x, y)

IMAGE DEGRADATION

The degradation process is modeled as a degradation function H(x,y) that, together with an additive noise term η(x,y), operates on an input image f(x,y) to produce a degraded image g(x,y), Fig. 2:

(1)

where, H represents a convolution matrix that models the blurring that many imaging systems introduce. The vectors g, f and η represent the observed, the original and the noise. More specifically, η is a random vector that models the random errors in the observed data.

Given g(x,y), some knowledge about the degradation function H(x,y) and the additive noise ç(x,y), the objective of restoration is to obtain an estimate of the original image f(x, y), as close as possible to the original input image.

Blurring: Blurring images can arise from many sources, such as limitations of the optical system, camera and object motion, astigmatism and environmental effects (Hansen and James, 2006). In some cases the blurring can be described analytically and thus it can be constructed from a function, rather than through experimentation.

Fig. 3:	Some blurring functions: disk, motion and Gaussian. (a) PSF Disk blurring function, (b) PSF motion blurring function and (c) PSF Gaussian blurring

Motion blurs which occur when there is a relative motion between the object and the camera during exposure. For example, the horizontal motion blur, which spreads a point source into a line, can be modeled using the following 1D degradation function:

(2)

where, L is the length of degradation vector.

Atmospheric turbulence blur commonly found in remote sensing and ariel imaging applications and this blur is due to long term exposure through the atmosphere and can be modeled by the Gaussian function:

(3)

where, K is a normalizing constant and σ² is the variance that determines the severity of the blur.

Uniform out-of-focus blur usually found in a variety of imaging systems as a uniform intensity distribution within a circular disk and the same is described by 2D degradation function. Some blurring functions: disk, motion and Gaussian is shown in Fig. 3.

(4)

Noise: In addition to blurring, observed images are usually contaminated with noise. Noise can arise from several sources and can be linear, nonlinear, multiplicative and additive. In our application, we consider a common additive noise model that comes essentially from the following three sources:

•	Photoelectric noise of background photons, from both natural and artificial sources. This kind of noise is typically modeled by a poisson process
•	Noise from electronics used to capture images, modeled usually by a white Gaussian noise, with zero mean and a fixed standard deviation proportional to the amplitude of the noise
•	Film grain noise, from the randomness of silver halid grains in the film used for recording
•	Quantization noise which occurs during image digitization

TOTAL VARIATION DEBLURRING

Tikhonov regularization: The main objective of regularization is to incorporate more information about the desired solution in order to stabilize the problem and find a useful and stable solution. The most common and well-known form of regularization is that of Tikhonov (Stand and Rudnicki, 2007). The Tikhonov regularized minimum norm solution of Eq. 1 is the vector F_δ∈U^N that minimizes the expression:

(5)

where, λ>0 is called a regularization parameter.

We denote:

(6)

Regularization can be understood as a balance between two requirements:

•	f should give a small residual Hf-g
•	f should be small in L2 norm

The regularization parameter λ>0 can be used to tune the balance. Note that in inverse problems there are typically infinitely many f satisfying Eq. 6.

Total variation regularization: Total variation is often used for image filtering and restoration. Total variation based filtering was introduced by Rudin et al. (1992) and computed using many algorithms as a constraint piecewise l1 minimization problem by Li and Santosa (1996). After that, it has been largely used for image restoration and deconvolution. Rudin et al. (1992) introduced the following idea:

Instead of minimizing:

(7)

Let us minimize:

(8)

Recall that:

And:

The idea is that Eq. 8 should allow occasional larger jumps in the reconstruction leading to piecewise smoothness instead of overall smoothness. It turns out that minimizing Eq. 8 is really a powerful method, but numerical minimization is more difficult. Finding the proper value for the parameter λ is an important problem. A large value of λ results in a smoother image and is necessary if the noise variance is high or H is highly ill-conditioned. On the other hand, a small λ blurs out the image details. So, one has to decide between smoothness and detail in the solution. Choosing the regularization parameter for an ill-posed problem is an art based principally on prior knowledge of the noise in the observations.

PARTICLE SWARM OPTIMIZATION (PSO)

The PSO algorithm was first described by Kennedy and Eberhart (1995). The application of PSO can be found in many engineering areas (Tlelo-Cuautle et al., 2010; Rezazadeh et al., 2009; Hassanzadeh and Saleh, 2008; Ying Yang et al., 2009). To enhance its usefulness, some strategies such as differential evolution have been linked to PSO (Otieno et al., 2010; Adeyemo and Otieno, 2009; Fraga and Schutze, 2009; Hernane et al., 2010). Also, fuzzy approaches have been applied Hooshmand (2008) and Flores-Becerra et al. (2009). Although, the PSO-based techniques have evolved greatly making the original version of the algorithm barely recognizable in the current ones, it is a stochastic population-based evolutionary computer algorithm for problem solving. It is a kind of swarm intelligence that is based on social-psychological principles and provides insights into social behavior, as well as contributing to engineering applications.

In a PSO system, a swarm of individuals (called particles) fly through the search space. Each particle represents a candidate solution to the optimization problem. The position of a particle is influenced by the best position visited by itself (i.e., its own experience) and the position of the best particle in its neighborhood (i.e., the experience of neighboring particles). When the neighborhood of a particle is the entire swarm, the best position in the neighborhood is referred to as the global best particle and the resulting algorithm is referred to as a gbest PSO. When smaller neighborhoods are used, the algorithm is generally referred to a lbest PSO. The performance of each particle (i.e., how close the particle is from the global optimum) is measured using a fitness function that varies depending on the optimization problem.

Each particle in the swarm is represented by the following characteristics:

xi	:	The current position of the particle
vi	:	The current velocity of the particle
yi	:	The personal best position of the particle
í	:	The neighborhood best position of the particle

The personal best position of particle i is the best position (i.e., the one resulting in the best fitness value)

visited by particle i so far. Let, F denote the objective function. Then, the personal best of a particle at time step t is updated as:

(9)

For the gbest model, the best particle is determined from the entire swarm by selecting the best personal best position. If the position of the global best particle is denoted by the vector í, then:

(10)

Where:

(11)

where, s denotes the size of the swarm.

The velocity update step is specified for each dimension j:

j∈{1,......,Nd}

Hence, v_i,j represents the jth element of the velocity vector of the ith particle. Thus the velocity of particle i is updated using the Eq. 12:

(12)

Where:

(13)

(14)

where, ω is the inertia weight, C₁ and C₂ are the acceleration constants and r_1,j, r_2,j are random coefficients distributed as:

(15)

The position of particle i, x_i is then updated using the following equation:

(16)

This process is repeated until a specified number of iterations is exceeded, or velocity updates are close to zero. The quality of particles is measured using a fitness function which reflects the optimality of a particular solution. The following steps summarize the basic PSO algorithm (Omran, 2004):

Algorithm:

It is important to clarify that good choice of the initial population can make the PSO faster to the global minimum, for this reason some works used the normal cloud method to find the best initial populations (Wu et al., 2008). In this study, for the presented case, the choice of the initial population is made by a simple instruction:

(17)

where, k is the dimension of the objective function.

APPLICATION TO IMAGE RESTORATION

This study will establish a new approach that can be used to solve a constrained optimization ill-posed problem in order to improve a blurred or noisy image. We will add many types of degradation functions to matrices of different sizes and then try to restore the original. Our starting image is a gray-level image contained in the mxn matrix. Each element in the matrix represents a pixel's gray intensity between black and white (0 and 255).

Assume, we know how fast the blurring function operator is known. The simplest approach is to solve the least squares problem:

(18)

In practice the results obtained with this simple approach tend to be noisy, because this term expresses only the fidelity to the available data g. To compensate for this, a regularization term below is added to improve smoothness of the estimate:

(19)

where, L is the discrete Laplacian, which relates each pixel to those surrounding it. L = del2(X) is a discrete approximation of:

(20)

where, X is the estimated matrix. The matrix L has the same size as X with each element equal to the difference between an element of X and the average of its four neighbors.

Since, we know we are looking for a gray intensity, we also impose the constraint that the elements of X must fall between 0 and 255.

To obtain the deblurred image, we want to solve for X:

(21)

We can implement our objective function using this expression; the number of variables in this objective function to be minimized will be mxn which is the size of the original matrix representing the original image.

EXPERIMENTAL RESULTS

We have carried out computer simulations to validate the applicability of our algorithm for image restoration. Some simple images of sizes: 16x16, 24x24 and 40x40 created by the MATLAB function checkerboard are used with different PSO parameters as follows: C1 = 1.5; C2 = 4-C1; minInertia = 0.3; maxInertia = 0.95; Swarm Size = 10, 20, 50, 120; Maximum Iterations = 20, 50, 100, 200. We run the algorithm using Intel Pentium 4 PC with 1.80 GHZ CPU and memory size of 1Go. The average processing time is dependent upon computation machine, image size and choice of PSO algorithm parameters (varies from few seconds to few minutes).

To evaluate restoration performance of our approach quantitatively, we record evolution of the Root Mean Squared Error (RMSE) and the Peak Signal to Noise Ration (PSNR) on Table 1, 2 and Fig. 4:

and

In order to evaluate the performance of the proposed approach when applied to neutron radiography images, we make use of corrupted images due to degradation (neutron flux perturbation) and noise (gamma radiations) taken by neutron radiography in a hostile site. Computer simulations prove that PSO algorithm yields excellent results and presents good efficiency in neutron images restoration (Fig. 5, 6).

Using the (a) checkerboard test image, (b) blurred and noisy and (c, d, e, f) restored using different swarm sizes and process iterations, with and without Laplacian regularization constraint in Fig. 7 and 8. We can observe the effect of regularization both on image quality and cost function evolution, Fig. 9a, b.

Table 1:	Evolution of the Mean Squared Error (MSE) and the Peak Signal to Noise Ratio (PSNR) with swarm size and number of iterations (without regularization)

Table 2:	Evolution of the (MSE) and the Peak Signal to Noise Ratio (PSNR) with swarm size and number of iterations (with regularization)

Fig. 4:	RMSE and PSNR evolution with swarm size and number of iterations

Fig. 5:	Neutron radiography image restoration with regularization and swarm size and iterations: (120, 200), (a) original, (b) degraded and (c) restored

Fig. 6:	Neutron radiography image restoration (with regularization) of mixing light water (H2O) with heavy water (D2O), (a) original images, (b) blurred images and (c) restored images

Fig. 7:	Restoration of blurred and noisy images without regularization constraint: (a) original, (b) noisy and (c-f) restored with swarm size and iterations: (10,20), (20,50), (50,100) and (120,200)

Fig. 8:	Restoration of blurred and noisy images with regularization constraint: (a) Original, (b) noisy and (c-f) restored with swarm size and iterations: (10,20), (20,50), (50,100) and (120,200)

Fig. 9:	Evolution of the cost functions with Iterations: (a) Without regularization and (b) with regularization

In Fig. 5 we make use of corrupted images due to neutron flux perturbation and noise from gamma radiations taken by neutron radiography in a hostile site. In Fig. 6, we extend our test to another important neutron radiography application, imaging a mixture of light and heavy water and see the difference at atomic level. Computer simulations illustrate that PSO algorithm yields excellent results and good efficiency in noisy images restoration. These results indicate that the use of Particle swarm optimization algorithm enables us to perform neutron radiography images restoration properly.

CONCLUSION

In this study, we have introduced the Particle Swarm Optimization (PSO) algorithm in image restoration to solve the ill-posed minimization problem based on total variation. The Laplacian constraint has been used for regularization to smooth deblurred images in the presence of noise. In our experiment, restoration of neutron radiography gray images, we can judge that the method always converges to acceptable results, from quality point of view, with computation time proportional to matrix (image) size. Different types of blurring and noises are tested with optimal regularization parameter λ chosen based on many trials and PSNR progress. To achieve further improvements in visual quality, intensive research is needed that could significantly improve performance in practical applications.