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. 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:
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:
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:
where, K is a normalizing constant and σ2 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.
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δ∈UN that minimizes the
expression:
where, λ>0 is called a regularization parameter.
We denote:
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:
Let us minimize:
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:
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:
Where:
where, s denotes the size of the swarm.
The velocity update step is specified for each dimension j:
Hence, vi,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:
Where:
where, ω is the inertia weight, C1 and C2 are the acceleration constants and r1,j, r2,j are random coefficients distributed as:
The position of particle i, xi is then updated using the following equation:
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):
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:
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:
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:
where, L is the discrete Laplacian, which relates each pixel to those surrounding it. L = del2(X) is a discrete approximation of:
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:
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.
