Image Reconstruction Using Particle Filters and by jasonndcosta


More Info
									IEEE TRANSACTIONS ON IMAGE PROCESSING, VOL. 19, NO. 5, MAY 2010                                                                                    1181

       Image Reconstruction Using Particle Filters and
               Multiple Hypotheses Testing
                          Noura Azzabou, Nikos Paragios, Senior Member, IEEE, and Frédéric Guichard

   Abstract—In this paper, we introduce a reconstruction frame-                     each pixel improves the quality of restoration compared to a
work that explicitly accounts for image geometry when defining the                   fixed size window. A further step in this direction consists in
spatial interaction between pixels in the filtering process. To this                 defining more appropriate means for selecting the spatial scale
end, image structure is captured using local co-occurrence statis-
tics and is incorporated to the enhancement algorithm in a sequen-                  of interaction between pixels. In this paper, we would like to
tial fashion using the particle filtering technique. In this context,                study the geometric structure of the image and if possible define
the reconstruction process is modeled using a dynamical system                      novel mechanisms for pixels selection that explicitly account
with multiple states and its evolution is guided by the prior density               for the observed local image structure. Such a process involves
describing the image structure. Towards optimal exploration of the                  two steps, (i) a learning stage where the image structure is
image geometry, an evaluation process of the state of the system is
performed at each iteration. The resulting framework explores op-
                                                                                    modeled, and (ii) a reconstruction step that consists of a novel
timally spatial dependencies between image content towards vari-                    mathematical approach to encode measurements driven by the
able bandwidth image reconstruction. Promising results using ad-                    geometric image model. Our aim is to introduce a strategy
ditive noise models demonstrate the potentials of such an explicit                  that allows a best possible selection of the pixels contributing
modeling of the geometry.                                                           to the reconstruction process driven by the observed image
  Index Terms—Additive noise, co-occurrence matrices, multiple                      geometry. Such geometry is expected to be anisotropic. More
hypothesis testing, nonparametric densities, particle filtering,                     explicitly, we take into account the image structure to retrieve
speckle noise, statistical models, structure enhancement.                           similar pixels to a reference one (that we want to reconstruct)
                                                                                    without any exhaustive scan of the image domain. To this end,
                        I. INTRODUCTION                                             the reconstruction is done in a sequential fashion using notions
                                                                                    of dynamical systems evolution. The idea is to consider mul-
T      HE core components of an image enhancement approach
       are: (i) the selection of the bandwidth or the scale of
interaction between pixels, and (ii) the way of introducing
                                                                                    tiple random walks/trajectories as potential candidates for the
                                                                                    filtering window and evaluate the pertinence of each trajectory.
                                                                                    Such a process involves two issues (i) the selection of the tra-
interactions between them based on the image content. Prior
                                                                                    jectory, (ii) and the evaluation of the trajectory appropriateness.
art consists of neighborhood filtering methods, PDE-based
                                                                                    Each walk is composed of a number of possible neighboring
methods and energy minimization approaches [3], [6], [22],
                                                                                    sites/pixels in the image which are determined according to
[25], [26]. Neighborhood filters consist in performing a
                                                                                    the observed image structure. Towards optimizing the selec-
weighted sum of noisy pixels over a given neighborhood
                                                                                    tion of candidate pixels within a walk as well as the overall
(filtering window). Such approaches include the sigma filter
                                                                                    performance of the method, image structure at local scale is
[17], the bilateral filter [24] and the NL-mean filter [5] and
                                                                                    considered through a learning stage. It consists in computing a
the fast version of the NL-mean filter [18]. In spite the fact
                                                                                    probability density function that describes the spatial relation
that these methods include the image structure in the weight
                                                                                    between similar image patches in a local scale. Random pertur-
definition, the selection of filtering window is always fixed in
                                                                                    bations according to these densities guide the “trajectories” of
shape and size. The main advantages of these methods is com-
                                                                                    a discrete number of walkers. The weighted integration of the
putational efficiency, but they are constrained by the amount
                                                                                    average values of these walks is the reconstructed image value.
of information present at the considered window. Examples
                                                                                    Image similarities, spatial smoothness as well photometric
of variable size filtering windows are introduced in [2], [10],
                                                                                    consistency for the samples are the criteria used to determine
[14], [21] where the authors propose an iterative algorithm to
                                                                                    the quality of the reconstructed value associated to a given
select the most appropriate bandwidth for each pixel. Their
                                                                                    walk. The overview of our method is presented in Fig. 1.
experimental results show that adapting the bandwidth for
                                                                                       The remainder of the paper is organized as follows. In Sec-
                                                                                    tion II, we present the learning of image structure at a local
   Manuscript received August 24, 2008. First published December 01, 2009;          scale. Random walks and particle filters are presented in Sec-
current version published April 16, 2010. The associate editor coordinating the
review of this manuscript and approving it for publication was Dr. Luminita         tion III. The Section IV is devoted to the application of the
Vese.                                                                               particle filtering to image denoising for the case of additive
   N. Azzabou and N. Paragios are with the Ecole Centrale, Grande Voie des Vi-      Gaussian noise. In Section V, we present experimental results
gnes, 92 295 Chatenay-Malabry, France (e-mail: n.azzabou@institut-myologie.
org;                                                        and comparisons with other filtering methods.
   F. Guichard is with the DxOLabs, 92100 Boulogne Billancourt, France
e-mail:                                                             II. STATISTICAL DESCRIPTION OF IMAGE STRUCTURE
   Color versions of one or more of the figures in this paper are available online
at                                                        During the past three decades, significant efforts were de-
   Digital Object Identifier 10.1109/TIP.2009.2037468                                voted to provide an efficient tool towards characterization of
                                                                 1057-7149/$26.00 © 2010 IEEE
1182                                                                               IEEE TRANSACTIONS ON IMAGE PROCESSING, VOL. 19, NO. 5, MAY 2010

                                               Fig. 1. Overview of “Random Walks” based image enhancement.

natural images content. In the context of computer vision and                 placement vector               . For a given intensity   it is defined
in particular image reconstruction, the performance of the de-                as
signed solutions is highly dependent on the accuracy of such
models. Image models can be either global or local. Global                                            such that                                 (2)
models focus on the statistics over the entire domain, while                  which corresponds to the probability of observing the intensity
local ones aim to capture the local co-dependencies within the                   at the position        knowing that             . Such co-oc-
observed image structure. The most popular work in this field                  currences will likely be sensitive to noise, and, therefore, the
refers to the generalized Laplacian distribution of the wavelet               constraint of exact matching could be relaxed, leading to (3)
coefficients distribution [19] as well as the joint statistics of pairs
or triplets of pixels in the wavelet domain [13], [16]. These sta-                                               such that
tistical models are generic for all natural images but they are
global and they do not provide a clear understanding on how the                                       and                                       (3)
information is spatially distributed in the image. The Gray Level
Co-occurrence Matrix (GLCM) [12] is a similar tool for statis-                with a constant determined according to the observed noise
tical description and characterization of texture. It is defined for           level.
a displacement vector                    by the          matrix                  With this formulation, and for each gray level in the image we
(where is the number of gray levels in the image). The coeffi-                 associate to each spatial transition the probability of leading
cient              is the number of occurrence of the pair of gray            to the similar pixels. It is important to note that this probability
level and after a displacement . In other words, it is pro-                   measure relies only on gray level similarity which is not suffi-
portional to the joint probability of the intensity of two pixels             cient to characterize image structure. For instance, if we con-
that are separated by the displacement vector                                 sider an image that involves texture as well as smooth regions,
                                                                              and assume the same observed intensity in both classes, then
                                                                              a generic model will fail to capture the distinct characteristics
                                                                              of the two classes. Therefore, we will rather consider the local
                                                               (1)            mean instead of the intensity at a pixel level towards dealing
 1 is the observed noisy image, the spatial domain and                        the the noise effect and the local variance towards dealing with
refers to the number of elements in a set. This matrix is of sig-
                                                                              multiple image components structure with similar appearance.
nificant importance since implicitly it provides information on
                                                                              We used the local variance since it is a simple primitive ca-
the geometric structure of the image.
                                                                              pable to describe texture at local scale. Considering both local
   As far as our denoising algorithm is concerned, we aim to find
                                                                              mean and local variance , would result in a more specific pdf
a description of the spatial relationship between similar pixels              describing spatial relation only between patches having close
belonging to the same structure in the image. Incorporating such              mean intensity and variance. The probability measure that de-
information will help the algorithm to cope with different image              scribes the spatial relationship between similar patches knowing
components leading to an adapted filtering window. More ex-                    their local mean and variance is then
plicitly, given the position of an observation, our objective is
to know what is the most likely pixel position in the vicinity                                                 such that
with similar image content. This could help to define the transi-
tions on the image lattice towards finding similar pixels within                                                              and
a window of a fixed radius . Such a move is by definition
anisotropic and is deduced from the observations. For this pur-                                                                                 (4)
pose, we introduce a probability function with respect to the dis-
                                                                              where       (resp       ) corresponds to the function that com-
  1I   = U + n U is the noise free image and n is the additive noise.         putes the local mean (resp local standard deviation) in a given

                                                                           interactions between similar pixels was restricted to a local
                                                                           neighborhood of size . The generalization to larger scale will
                                                                           be done in a sequential fashion and this will be explained in the
                                                                           following section.
                                                                              Given this density, our aim is to determine the most appro-
                                                                           priate set of neighbors to estimate the noise-free intensity of a
                                                                           given pixel without an exhaustive scan of a large image domain.
                                                                           To this end, one needs to define a strategy to generate neigh-
                                                                           borhood candidate windows that takes into account the image
                                                                           content. In the present work, this is done through particle filters
                                                                           technique and this will be the focus of the following section.

                                                                                III. OVERVIEW OF PARTICLE FILTERING TECHNIQUE
Fig. 2. Two pdf distributions p (d) for different values of  and  [top      As stated before, our purpose is to determine the most ap-
( = 39;  = 11:67)], bottom [ = 96;  = 3:55)], and sample genera-
tion according to these pdf (red pixel) for two different positions.       propriate window or neighborhood shape and size to estimate
                                                                           the image intensity in a given position. One attempt to do that
                                                                           was presented in [23] where the authors perform filtering by se-
position.        are constants to be fixed according to the noise           lecting the neighboring pixels in a random fashion but without
level.                                                                     taking image structure into account. The image lattice is repre-
   We must point out that we compute the probabilities of dis-             sented using a graph where the weights of vertices encode the
placement values included in a window of radius which cor-                 gray level difference between pixels. The transition of a walker
responds to the local scale. The selection of such a parameter is          from one position to another is based only on these weights
important and must be adapted to the image content. Choosing a             without considering the previous positions or local image sta-
high value for increases computation complexity while a small              tistics. In our algorithm, we explore the image domain using a
value for the scale parameter would result in a less accurate de-          sequential procedure taking into account the image structure as
scription of the structure.                                                well as the previous observation and positions. To this end, we
   The outcome of this process consists in a probability den-              will consider the trajectory of a particle that starts from the pixel
sity function           that aims to find a spatial representation          that we want to recover and moves to other neighboring ones.
of different structures through the computation of the relative            Each trajectory corresponds to a statistical hypothesis and will
position of similar patches. It is obvious that such a density is          result in a reconstructed value. One should address two aspects
far from being parametric due to the randomness of the observed            in such a context: (i) the trajectories themselves or transitions
image geometry. In practice the estimation of the pdf is based on          between stages, and (ii) the “appropriateness” of each trajec-
a nonparametric kernel density approximation [27] like Parzen              tory given the origin and the associated observations. Such ap-
windows.                                                                   propriateness should depend on the consistency of the observa-
   Let                                              and                    tions, the ability to encode local image structure and the rela-
                . Then, the probability density function estimator         tionship between the origin value and the contributing pixels.
is given by the following expression                                       Toward optimizing the trajectory computation, we have to es-
                                                                           timate the probability density function of visiting a position
                                                                    (5)    at time knowing all the past observations that refer to the in-
                                                                           tensities relative to the previously visited positions. This pdf is
                                                                           noted               where is the current position of the particle
with      being a symmetric definite positive kernel. Gaussian              and        are the observation related to the previous and current
kernels are the most common selection of such an approach and              state of the particle that are the different positions and image
they were considered in our case to approximate             . The          properties in these positions.
bandwidth is a diagonal matrix with coefficients equal to                      In order to estimate such a pdf, we use particle filtering tech-
  and           . In terms of the mean and variance quantization           niques that are an efficient sequential Monte-Carlo methods in-
values, we considered a quantification step           . Examples of         troduced in [11]. We will give a brief overview of the use of
these densities are shown in Fig. 2.                                       particle filtering as a way of sequential estimation of posterior
   Once the structural model has been constructed from the                 pdf and focus on the Sequential Importance Resampling tech-
image, we are able starting from a given image position                    nique. For an extended review the reader can refer to [1], [7],
characterized by the couple          , to determine the probable           [8].
transitions towards the pixels that belong to the same structure
as . This property is illustrated in Fig. 2, where several posi-           A. Preliminaries
tions (in green) similar to the origin (in red) are drawn without             Let us introduce first the problem of Bayesian filtering where
scanning its entire neighborhood. The displacements with                   one has to estimate the evolution of a hidden state of a system
important probability are those that guarantee local statistics            given a set of observations. Such a system involves a function
preservation in terms of appearance. We recall that the learning           that relates states with observations and a function that relates
of the probability density function that describes the spatial             the present state with the past ones. To introduce this notion,
1184                                                                       IEEE TRANSACTIONS ON IMAGE PROCESSING, VOL. 19, NO. 5, MAY 2010

let us note     the state sequence and             the sequence of        To summarize, the concept of particle filtering consists of
measurements that verify                                               three main steps:
                                                                          • hypothesis generation according to the transition law
                                                                          • computation of the likelihood of observations generated by
                                                                             this hypothesis             ;
where                                is the transition function,          • weight updating according to                            .
                             is the observation function while         A common problem in this estimation strategy is the samples
          and       are independent and identically distributed        degeneration, where many particles carry on less information
noise sequences.                                                       and the variance in the weight increases. In such a situation,
   Filtering is then equivalent to determine an estimation of the      many particles have their weight close to zero and do not con-
state vector      knowing the set of all available observations        tribute to the approximation of               . To overcome this
       . In the Bayesian framework, one has to compute the pos-        problem, a resampling step is necessary. The basic idea is to
terior pdf             which can be obtained recursively in two        eliminate the particles that have limited capacity to approximate
stages: prediction and update. The main idea is to suppose that        the density and focus on the ones which substantially contribute
the pdf                   is available and calculate          via      to the reconstruction of the density (important weights). Many
the Bayes rule (6)                                                     algorithms were proposed to perform particle selection. In the
                                                                       present work, we considered the systematic resampling scheme
                                                                (6)    [15] which consists in eliminating particle with small weight
                                                                       and making several copies of the other particles according to
                                                                       their weight.
where the prior pdf is computed via the Chapman-Kolmogorov
                                                                          Now towards the application of such an algorithm in the con-
                                                                       text of image denoising we have to define the following ele-
                                                                          • a state vector     at a given time that refers to the state of
                                                                             the reconstruction process. It corresponds to the position
where                                                    ,                   of the new pixel added to the filtering window, as well as
is defined by the observation model and                  by the tran-         the currently reconstructed value;
sition model. One has then to evaluate the expression (7) in order        • an observation which corresponds to local image statistics
to determine the probability of having the state      given the ob-          (patch, local mean, local standard deviation) relative to a
servations        . The recursive computation of the prior and               new candidate pixel to the current state vector;
the posterior pdf leads to the exact computation of the poste-            • a measure that evaluates the contribution of each new can-
rior density. Nevertheless, in practice, it is impossible to com-            didate pixel to the estimation process;
pute explicitly the posterior pdf              , and, therefore, an       • a systematic way of generating new samples given the
approximation method needs to be introduced.                                 existing ones able to account for the expected image
B. Sampling Importance Resampling Filter (SIR)
  The key idea is to approximate the posterior pdf by
                                                                                 IV. APPLICATION TO IMAGE RESTORATION
random state samples                       with     associated
weights                      that reflect the importance of the            Thus, given an origin pixel      the process aims to recover a
associated sample in the pdf. Thus, the posterior pdf can be           number of “random” positions                                 with
approximated by a discrete weighted sum                                similar properties to to reconstruct the corrupted origin value
                                                                              . To this end we have to estimate in a sequential fashion
                                                                       the posterior               using the previously described algo-
                                                                       rithm. One can imagine multiple particles that are moving ran-
                                                                       domly on the image domain. The path (or random walk) of a par-
where    is the Kronecker Delta function defined as              if     ticle corresponds to a candidate filtering window. The weighted
        and           if       .                                       mean of the intensities along a given particle trajectory is used
   The generation of the samples is done through the principle         to restore the intensity of . As far as observations are con-
of Importance Sampling [4], [7]. This approach suggest an itera-       cerned, they correspond to the local patch around the pixel
tive scheme where the particle are generated using the transition      and around the candidate . To perform filtering, our approach
model               and then compute the particle weights itera-       requires the definition of a perturbation model as well as a like-
tively according to                                                    lihood measure that reflects the contribution of a trajectory to
                                                                       the denoising process.
                                                                       A. Transition Model
Hence, particle weights are updated using two main informa-
tion: the observation pdf which reflects the likelihood of having          An important issue for the proposed algorithm consists in
an observation knowing the state        and the transition model       defining an appropriate strategy for samples perturbation. An
which controls the evolution of a particle state.                      intuitive solution is to perform Gaussian perturbations. Such

a choice, suggests an isotropic way for particle transitions             where      is the bandwidth which must be carefully se-
which can be justified in the case of homogeneous regions                 lected to get a reliable measure of similarity while being
but not in the presence of some texture or edges. An efficient            computationally efficient.
manner to adapt the particle transition to the image content           • Another alternative to the    distance is the Sobolev dis-
is to rely on the statistical model for image structure intro-           tance between patches defined as
duced in the (Section II). The distribution       determines the
transition model between particle at position         and      .
We recall that            is a probability density function that
indicates for a pixel     (whose local mean and variance are
equal to and ) the probability that the destination position
        according to the transition vector      has close local                                                                     (11)
statistics to   . This pdf is computed in the beginning of
the process for each couple (mean, variance). We recall here               Such a distance accounts also for the image gradient which
its definition when considering local statistics of      . Using            is an important information on image structures and their
                                                            and            orientation. Thus, similar patch have to be similar both in
                         , transition model is                             sense of their intensity and their gradient.
                                                                     In addition to the measure of the quality of a new position in the
                                                              (9)    particle walk, we consider a measure of its consistency within
                                                                     the whole trajectory. This measure is the intravariability of in-
                                                                     tensities in the different positions visited by the particle during
Transitions according to this pdf guarantee that a maximum           its walk. This constraint gives more importance in the filtering
number of particles will explore sites that are similar to the       process to the pixel set composed of element having homoge-
origin pixel . Such a strategy makes possible the use of a re-       neous properties. This is mainly interesting in the case of edges
duced number of particles and could also prevent their degen-        or other linear structures in the image, where walks that follow
eration. The essence of the proposed algorithm is to explore the     exactly the structures would be favored compared to other tra-
image domain in a random fashion and to find the maximum              jectories that fail to take into consideration this structure during
number of similar sites to the origin one. Our approach does         their evolution. This feature can enhance in an efficient manner
perform this task and explore the image domain with a smart          image structures. Hence, in order to account for the intravari-
fashion such that particles trajectories follow image structure      ability of the trajectories, we consider as an observation the walk
without an exhaustive scan of the pixel neighborhood. In this        variance, with respect to the origin value
context, the statistical structural description of the image geom-
etry is a convenient transition model. Once such transitions have
been identified, the next step consists of defining a likelihood                                                                      (12)
measure for each particle based on its trajectory and the inten-
sity observation relative to the particle position.                  Note that such an observation provides a measure of the “uni-
                                                                     formity” of the trajectory and could also be determined within
B. Likelihood Measure                                                a larger neighborhood (not at the pixel level). Finally, the im-
   Within the proposed approach, the random walk that is the set     portance of a new sample given the prior state of the walk and
of the image sites contributing to the reconstruction of a given     the observation that is the patch around   is then defined as an
pixel is determined using a sequential approach. To evaluate the     exponential function of the two metrics
relevance of a new particle state, a likelihood measure must be
introduced. This measure evaluates the significance of a new
candidate pixel position. The first criterion that is considered is
                                                                     where     and      are constants that determine the bandwidth of
the similarity between the image content observed in the origin
                                                                     the weight computation function. The choice of these values de-
pixel     and the current path position . We adopt the image
                                                                     pends on the noise level of the image. High values of     and
patches comparison, motivated by the fact that measuring simi-
                                                                     would relax the constraints on trajectories evolution and would,
larities between image patches is an efficient and popular metric
                                                                     therefore, result in a loss of image details.
used in many computer vision tasks such as texture synthesis,
restoration etc [5], [9]. In such a context we define the observa-    C. Intensity Reconstruction
tion model as follows.
   • The      distance between the local patch around the origin        After having defined the transition model as well as the
      pixel    and the one corresponding to a current particle       observation model, we focus on the intensity reconstruction
      position . Thus, if we call         this similarity measure,   process. To this end, for each pixel          of the image, we
      we have                                                        generate      particles by applying        perturbations starting
                                                                     from the initial position      and according to the transition
                                                                     law                                    . Iterating the transition
                                                                     process time for each particle yields a walk that is the set
                                                                     of successive positions of the particle                     . We
1186                                                                              IEEE TRANSACTIONS ON IMAGE PROCESSING, VOL. 19, NO. 5, MAY 2010

                                                                              noise model, but this framework can be used with different noise
                                                                              models. The extension of the method to cope with different
                                                                              noise models should involve both the statistical modeling of
                                                                              co-dependencies between pixels as well as the likelihood used
                                                                              to determine the importance of a new sample.
Fig. 3. Example of three particles walks (red pixel) starting from the same
origin position (green pixel), the walks follow image structure.              D. Comparison With NL-Mean
                                                                                 In the recent years the NL-mean [5] algorithm has become
                               TABLE I                                        the state of the art in image denoising. In this section we will
                                                                              provide the theoretical differences between the two algorithms.
                                                                                 The NL-mean algorithm assumes fixed size with respect to
                                                                              the local filtering window (          window centered at the origin
                                                                              pixel). Then based on the similarity between the center patch
                                                                              and the candidate patches it performs filtering. This is not the
                                                                              case of the proposed approach where the filtering window is
                                                                              adapted to the pixel content. More explicitly, the set of candidate
                                                                              pixels in our approach is not fixed and changes per pixel location
                                                                              according to local pixel properties. In other words the filtering
                                                                              window is reconstructed on line to ensure a selection of greater
                                                                              number of similar local patches to be used in the reconstruction
                                                                              with minimum computation effort. In our filtering process we
                                                                              consider local neighborhood then we enlarge it iteratively fol-
define then the walk value          as the average weighted value              lowing the structure present in the image. This is done through
along the walk. It corresponds to the value used to reconstruct               the transition step of the particle filter where we select in each
the original pixel according to the walk of the particle                      position the most appropriate transition to find similar patches.
                                                                              The way the candidate pixels are selected is made according
                                                                              to the probability of co-occurrence of similar pixels at a given
                                                                      (14)    scale and given orientation. Towards recovering these probabil-
                                                                              ities, we perform a nonparametric density approximation per
As we discussed above, the pertinence of a walk as well as the                each couple (intensity/variance). Hence, for each couple (inten-
associated reconstructed value is computed using the likelihood               sity/variance), we define the probability of stepping into a sim-
measure         defined in expression (13). Once the weights of                ilar pixel given a specific transition. This aspect of our algorithm
all the particles are determined, we normalize them so that their             is the fundamental difference with respect to the NL-mean ap-
sum is equal to one. Next, we perform resampling in order to                  proach as well as its fast version proposed in [18]. Actually, the
replace the particles with insignificant weights. These particles              technique proposed in [18] assumes a fixed size window, where
are not consistent with the reference pixel      and using them               local image characteristics like variance and gradient orientation
will introduce a bias in the reconstruction process. Finally, the             are only used to make pixels preselection to avoid unnecessary
reconstructed value is a weighted average of the mean intensity               computation of weight values. The number of pixels considered
of each walk and defined as                                                    in the denoising is simply defined by the initial size of the fil-
                                                                              tering window. In our case the local variance is used to describe
                                                                              the image structure and the spatial relationship between similar
                                                                      (15)    pixels. This information will be used to reconstruct iteratively
                                                                              the filtering window. To illustrate this concept, we considered a
The whole process (transition, weight computation, and resam-                 textured image and one pixel in this image (marked in red cross)
pling) is repeated (T) times. An example of particle walks is                 and we tried to compute the set of candidate using NL-mean
shown in Fig. 3, where the pixel in green refers to the origin pixel          with a fixed window size centered on the pixel (19 19) and
we want to denoise and each trajectory is a candidate filtering                our algorithm with similar complexity using (60 particle and
window. In practice we use                particles, with                     six iterations). The results are presented in Fig. 4. We can see
pixels contributing to each walk. To illustrate the random walks              that one iteration of our algorithm using 60 particles detects an
filtering an overview of the whole process is presented in Fig. 1              equivalent number of similar pixels to the central one [Fig. 4(a)
and the denoising algorithm is given in Table I. We should point              and (b)]. After six iterations, the particle filter is able to recover
out here, that the process is performed individually per pixel;               more candidate pixels at larger scale. This is the main strength
therefore, the reconstruction is always done using the observed               of our algorithm, where the candidate pixels are sampled from
image and there is no interaction between reconstructed values.               a distribution that takes into account the local pixels statistics.
   It is important to note that the model described above and                 Regarding the NL-mean, the number of candidates is related to
mainly the likelihood measure assumes an additive zero mean                   the research window size and the image content is not taken
AZZABOU et al.: IMAGE RECONSTRUCTION USING PARTICLE FILTERS AND MULTIPLE HYPOTHESES TESTING                                                         1187

                                                                                     algorithm to suppress noise while preserving image details. In
                                                                                     addition to that, we considered a more objective indicator which
                                                                                     is the Peak Signal to Noise Ration (PSNR) defined as

                                                                                     In the first part, we will present comparison with other state
                                                                                     of the art techniques. Next, we will focus on the performance
                                                                                     difference between the NL-mean algorithm and the proposed

                                                                                     A. Comparison With State of the Art
                                                                                        For performance comparison, we considered competitive ap-
                                                                                     proaches based on filters that encode implicitly the image struc-
                                                                                     ture through the definition of the weights such as the bilateral
Fig. 4. Example of textured image where the pixel in red refers to the cen-          filter [24], the NL-mean [5], and a neighborhood filtering tech-
tral pixel that must be denoised. (a) The set of candidate pixels provided by        nique using variable bandwidth size [14]. Furthermore, we have
the NL-mean. (b) The set of candidate pixels provided by one transition of the       considered PDE’s based approaches that are the total variation
particle filter. (d) The set of candidate pixels provided by six transitions of the
particle filter.                                                                      [22] and the anisotropic filtering [20] using an edge stopping
                                                                                     function of the type                                    . Our sim-
                                                                                     ulation consists in adding a synthetic white Gaussian noise of
into account for the selection of the size and shape of the fil-                      known standard deviation                 to the set of images fre-
tering window. This is also the case of the fast denoising algo-                     quently adopted by the image restoration community for val-
rithm based on NL-mean [18]. We have to emphasize that, if we                        idation. Their parameters were tuned manually to achieve op-
consider the very simple case where the co-occurrence densi-                         timal performance in terms of PSNR and visual quality. Re-
ties are set to uniform with respect to the displacement vector,                     garding our approach, each walk consists of six steps and for
then we can assume to a certain extent that our approach is                          each pixel we have considered 60 particles                        .
equivalent in terms of bandwidth and filtering window with the                        Given this choice we are able to explore 6 60 positions in the
NL-mean. However, in the presence of structure and textured                          image for each pixel. The radius of the explored region is of the
component, our technique adapts the search window to these                           order          where is the maximum distance that a particle
structures. As far as weight computation is concerned, like the                      can achieve in a transition (in our case       ). The window size
NL-mean case we can also encode any imaginable patch-based                           used for likelihood computation is 7 7                 the parame-
similarity (intensity-driven or probability-driven) in the process                   ters of the Gaussian used for likelihood computation (13) are set
through a slight modification of the cost related with the ap-                        to                  . The same set of parameters was considered
propriateness of the considered walk. Furthermore, opposite to                       when the Sobolev distance between patches was used. In such a
                                                                                     case, the image gradient is determined using a smoothed version
simple patch-based similarities weights being considered in the
                                                                                     of the image by a Gaussian kernel. For the bilateral filter, the pa-
case of NL-mean, we adopt a metric that combines this mea-
                                                                                     rameters are           and          . The NL-mean filter was used
sure with a consistency of the walk. Such consistency is de-
                                                                                     with            , the patch size used for comparison is 7 7 and
fined according to the homogeneity of the samples being vis-
                                                                                     the weighted average was computed over a local neighborhood
ited during the walk. Such a choice introduces a more robust
                                                                                     of size 11 11. As far as the variable bandwidth method is con-
behavior with respect to isolated and unstructured data where
                                                                                     cerned, the results were provided by the authors and one can find
the notion of geometric consistency across candidate pixels for                      more details in [14]. In Table II, we provided the PSNR values
filtering is not considered. Last, but not least the method is mod-                   relative to each method. From this table, we can conclude that
ular, gradient-free and can be used to deal with various models                      in general, computing similarity on image patches is better than
of noise. Finally, to compare performance of the two algorithms,                     the pixel-wise comparison. The Bilateral filter, the total varia-
we will present experimental results in the next section.                            tion minimization technique (TV), and the anisotropic diffusion
                                                                                     (AD) have small PSNR values because they are not based on
                      V. EXPERIMENTAL RESULTS                                        patch comparison which provides a best measure of similarity
   In order to evaluate the performance of the proposed frame-                       than the pixel wise comparison or the gradient information. The
work, the case of additive Gaussian noise was considered and                         approach we presented and the NL-mean algorithm have about
compared with some of the state of the art methods with assump-                      the same performance, while the method presented in [14] out-
tions/aims similar to the ones of our approach. We used natural                      performs all the others which is natural since the underlying idea
images corrupted by a additive Gaussian noise with different                         of this approach is to minimize a bound of the quadratic error
standard deviations. Performance evaluation involves a qualita-                      mean. Moreover, our algorithm selects the patches for compar-
tive criterion (visual assessment) that reflect the ability of the                    ison randomly while being guided at the same time by the image
1188                                                                                IEEE TRANSACTIONS ON IMAGE PROCESSING, VOL. 19, NO. 5, MAY 2010

                                                                 TABLE II

Fig. 5. (a) Noise free image. (b)Image corrupted with Gaussian noise  = 20.   Fig. 6. (a) Restoration result obtained with NL-mean algorithm. (b) NL-mean
(c) Bilateral filter restoration. (d) Residual of the bilateral filter.          residual. (c) Restoration result obtained with random walks algorithm (L dis-
                                                                               tance). (d) Residual of the random walk based approach.

structure. The number of pixels contributing to the walks has a
tremendous impact on the performance of the method. When
only a few samples have an important weight in the filtering
process, it produces under smoothing in homogeneous areas and
the texture is not well reconstructed. On the other hand, the vari-
able bandwidth approach [14] involves all the pixels within a
neighborhood in the filtering process and uses more candidates
than our method. We recall that for our method, in homogeneous
zones, particles evolve in an isotropic fashion which is equiva-
lent with considering a window of a fixed bandwidth. In tex-
tured regions, hypothesis generation relies on image structure,
and, therefore, capturing local image characteristics is favored.
In terms of photometric distance comparison, one can observe
that the Sobolev norm has a slightly better performance than the
    distance because it involves other features than the image
   The qualitative results presented in Figs. 5–7 show that the
restoration quality of our algorithm when compared to the one                  Fig. 7. (a) Restoration result obtained with random walks algorithm (Sobolev
using an adaptive window size as well as the NL-mean algo-                     distance). (b) Residual of the random walk based approach. (c) Restoration re-
rithm is satisfactory. In Fig. 8, restoration results emphasizing              sult obtained with KB06 [14]. (d) KB06 [14] residual.
an area with a very fine texture scale are presented. This figure
[as well as Fig. 7(c)] shows the presence of some artifacts in the
denoised image using [14] due to over smoothing. Such a result                 of important noise level the efficiency of the learning step of the
is due to the large window size selected in this region that re-               image structures is compromised. In such a situation the transi-
duced excessively the variance of the estimator. This is not the               tion model is not too robust with respect to noise. Considering
case for our algorithm, as it yields results with a more natural               such a model for sampling new particle positions could result
aspect. In terms of qualitative results, our approach does equally             in small numbers of relevant pixels which may not be sufficient
well if compared to the NL-mean algorithm. Note that in case                   for the restoration of some structures.
AZZABOU et al.: IMAGE RECONSTRUCTION USING PARTICLE FILTERS AND MULTIPLE HYPOTHESES TESTING                                                        1189

                                                             TABLE III

                                                                                   small. We can notice also that for equivalent complexity, our
                                                                                   method outperforms globally the NL-mean algorithm. However,
                                                                                   in the case of images that are composed of irregular texture and a
                                                                                   dominant local smooth component (like lena and house), an ex-
                                                                                   haustive search of local neighborhood (as done in the NL-mean)
                                                                                   outperforms our approach. In our case the sampling is not ex-
                                                                                   haustive because we use a small number of particles. Besides,
                                                                                   for stochastic texture our model reaches its limit because spa-
                                                                                   tial relationship between similar patches cannot be extended to
                                                                                   larger scale. For the case of images that present a significant
                                                                                   amount of texture (barbara, fingerprint, baboon), our method
                                                                                   performs better since we consider larger scale then the NL-mean
                                                                                   for texture restoration.
Fig. 8. Zoom on the Baboon image. (a) Original image. (b) Random walk
restoration. (c) NL-Mean restoration. (d) Adaptive window size restoration [14].                           VI. CONCLUSION
                                                                                      In this paper, we introduced a new denoising technique that
                                                                                   explicitly accounts for geometric relationships of image con-
B. Comparison With the NL-Mean
                                                                                   tent. This was achieved using a particle filtering technique with
   In order to better understand the potential of the method, we                   a perturbation model determined using the observed statistical
have considered for a fixed number of candidate pixels different                    behavior of the image. The structural model was learned in a
geometric exploitation strategies. This is achieved by varying                     first step in order to guide the particle evolution toward image
the number of particles       and the number of transitions to                     structure such as edges. Based on such a model, the approach
build the trajectories. In terms of computation complexity this                    acts like an isotropic filter when considering smooth areas, and
is equivalent to considering a window containing             pixels                anisotropic when dealing with edges. Besides, when we con-
for the NL-mean algorithm. The use of a small number of par-                       sider the case of regular texture with a repetitive pattern, the
ticles with a large number of iterations (or number of pixels per                  presented technique is able to recover similar patches in a large
walk) would have the ability to perform a sparse and elongated                     neighborhood using a small number of particles while visiting a
sampling and to seek similar structure far from the origin pixel.                  limited number of positions using a learning step at local scale
On the other hand, while considering small number of pertur-                       (Section II). However, in the case of random texture, gener-
bation, we perform local dense sampling to seek similar con-                       alizing at larger scale the spatial relationship between similar
tent within limited-bandwidth anisotropic window. In this set                      pixels learned at local scale is not straightforwrd. In such a case,
of experiments we compared the performance of the NL-mean                          the evolution of particles is not optimal and a higher number
using             , a patch size equal to 7 7 for weight com-                      of particles is needed. However, both the theoretical framework
putation while the weighted average was computed over a local                      and the experimental validation suggest that the use of the image
neighborhood of size 11 11. For our algorithm the patch size                       structure to determine the perturbation model and the considera-
was set to 7 7 and                             . The tested set of                 tion of multiple hypotheses testing can yield interesting results.
couple              are selected in order to have the same com-                    We have also shown that our method is flexible with respect to
plexity with NL-mean. In Table III, we provide SNR value of re-                    the choice of the similarity measure between image patches.
stored images. It is clear that the choice made between sparsity                      Improving the learning stage of the image structural model
versus dense sampling affects significantly the results. For in-                    and guiding the particles to the most appropriate directions
stance, considering a small number of particles and an important                   could be a step toward increasing the efficiency of particle
number of iterations will result in a large scale research window.                 transitions. On the other hand, an automatic technique to
In that case (          and           ), the SNR values show that                  capture the scale of texture should be considered. It enables
non local sampling is not always a good choice which is a nat-                     to determine the patch size considered when computing the
ural outcome since similar patches are more likely to be found at                  likelihood of a new particle state on one hand and the size
local scale in natural images. This statement is also valid for the                of the maximum radius of transitions on the other hand. The
NL-mean algorithm since using a large size neighborhood is ef-                     likelihood measure could be also modified to be more specific
ficient only with an adequate choice of         which must be very                  to each noise distribution and more robust. Last, but not least,
1190                                                                                       IEEE TRANSACTIONS ON IMAGE PROCESSING, VOL. 19, NO. 5, MAY 2010

the ability to perform the process in parallel for all pixels                            [26] L. Vese and S. Osher, “Modeling textures with total variation minimiza-
and benefit from the reconstructed values might improve the                                    tion and oscillating patterns in image processing,” J. Sci. Comput., vol.
                                                                                              19, pp. 553–572, 2003.
performance of the method.                                                               [27] M. Wand and M. Jones, Kernel Smoothing. New York: Chapman &
                                                                                              Hall, 1995.
    [1] M. Arulampalam, S. Maskell, N. Gordon, and T. Clapp, “A tutorial on                                     Noura Azzabou received the engineering diploma
        particle filters for on-line non-linear/non-gaussian Bayesian tracking,”                                 from the Ecole Polytechnique of Tunisia in 2001, the
        IEEE Trans. Signal Process., vol. 50, pp. 174–188, 2002.                                                M.S. degree in applied mathematics and computer vi-
    [2] N. Azzabou, N. Paragios, and F. Guichard, “Random walks, con-                                           sion from the Ecole Normale Supérieure de Cachan,
        strained multiple hypothesis testing and image enhancement,” in Proc.                                   France, in 2004, and the Ph.D. degree from the Ecole
        Eur. Conf. Computer Vision, 2006, pp. 379–390.                                                          Nationale de Ponts et Chaussees, France, in 2008.
    [3] N. Azzabou, N. Paragios, F. Guichard, and F. Cao, “Variable bandwidth                                   Her Ph.D. work was carried in the computer vision
        image denoising using image-based noise models,” in Proc. IEEE Int.                                     group in the Applied Mathematic Department, Ecole
        Conf. Computer Vision and Pattern Recognition, 2007, pp. 1–7.                                           Centrale de Paris and DxOLabs.
    [4] N. Bergman, “Recursive Bayesian Estimation: Navigation and                                                 Her research activities were focused on variable
        Tracking Applications,” Ph.D. dissertation, Linkoping Univ.,                                            bandwidth image models for texture-preserving en-
        Linkoping, Sweden, 1999.                                                      hancement of natural images. She is currently a research scientist in the RMN
    [5] A. Buades, B. Coll, and J.-M. Morel, “A non-local algorithm for image         Laboratory in the Myology Institute (France), working on medical image pro-
        denoising,” in Proc. IEEE Int. Conf. Computer Vision and Pattern              cessing and enhancement. Her research interests include registration, segmen-
        Recognition, 2005, pp. 60–65.                                                 tation, inhomogeneity correction, and denoising. She is the author of several
    [6] V. Caselles, J. M. Morel, G. Sapiro, and A. Tannenbaum, “Introduc-            paper published in intrenatioal conference and journals.
        tion to the special issue on partial-differential equations and geom-
        etry-driven diffusion in image-processing and analysis,” IEEE Trans.
        Image Process., vol. 7, no. 3, pp. 269–273, Mar. 1998.
    [7] A. Doucet, On Sequential Simulation-Based Methods for Bayesian Fil-
        tering, Dept. Eng., Cambridge Univ., Cambridge, U.K., 1998, Tech.                                        Nikos Paragios (SM’03) received the B.Sc. and
        Rep. CUED/F-INFENG/TR. 310.                                                                              M.Sc. degrees from the University of Crete, Greece,
    [8] A. Doucet, J. de Freitas, and N. Gordon, Sequential Monte Carlo                                          in 1994 and 1996, respectively, the Ph.D. degree
        Methods in Practice. New York: Springer-Verlag, 2001.                                                    from I.N.R.I.A./University of Nice/Sophia Antipolis,
    [9] A. Efros and T. Leung, “Texture synthesis by non-parametric sam-                                         France, in 2000, and the HDR (Habilitation a Diriger
        pling,” in Proc. Int. Conf. Computer Vision, 1999, pp. 1033–1038.                                        de Recherches) from the University of Nice/Sophia
   [10] K. Egiazarian, V. Katkovnik, and J. Astola, “Adaptive window size                                        Antipolis in 2005.
        image denoising based on ICI rule,” in Proc. IEEE Int. Conf. Acoustic,                                      He is currently a Professor at the Ecole Centrale
        Speech and Signal Processing, 2001, pp. 1869–1872.                                                       de Paris—one of most exclusive engineering schools
   [11] N. Gordon, D. Salmond, and A. Smith, “Novel approach to non-                                             “Grange Ecoles”—leading the Medical Imaging and
        linear/non-gaussian bayesian state estimation,” in Proc. IEEE Radar                                      Computer Vision Group at the Applied Mathematics
        and Signal Processing, 1993, vol. 140, pp. 107–113.                           Department. He is also affiliated with INRIA Saclay Ile-de-France, the French
   [12] R. M. Haralick, K. Shanmugam, and I. Dinstein, “Textural features for         Research Institute in Informatics and Control heading the GALEN group. Prior
        image classification,” IEEE Trans. Syst., Man, Cybern., vol. SMC-6,            to that, he was Professor/Research Scientist (2004–2005) at the Ecole Nationale
        pp. 610–621, 1973.                                                            de Ponts et Chaussees-France affiliated with Siemens Corporate Research
   [13] J. Huang and D. Mumford, “Statistics of natural images and models,” in        (Princeton, NJ) as a project manager (smart cameras and machine vision
        Proc. IEEE Int. Conf. Computer Vision and Pattern Recognition, 1999,          group) (2002–2004), senior research scientist (2004), and research scientist
        pp. 541–547.                                                                  (1999–2003). In 2002, he was an adjunct Professor at Rutgers University and
   [14] C. Kervrann and J. Boulanger, “Unsupervised patch-based image reg-            in 2004 at New York University and a visiting Professor at Yale University,
        ularization and representation,” in Proc. Eur. Conf. Computer Vision,         New Haven, CT, in 2007. He has co-edited four books, published more than
        2006, pp. 555–567.                                                            100 papers in the most prestigious journals and peer reviewed conferences of
   [15] G. Kitagawa, “Monte carlo filter and smoother for non-gaussian non-            computer vision and medical imaging, and holds nine U.S. issued patents and
        linear state space models,” J. Comput. Graph. Statist., vol. 5, pp. 1–25,     more than 20 pending. His research interests include computer vision, medical
        1996.                                                                         image analysis, remote sensing, and human computer interaction.
   [16] A. Lee, K. Pedersen, and D. Mumford, “The nonlinear statistics of                Prof. Paragios serves as a member of the Editorial Board of the International
        high-contrast patches in natural images,” Int. J. Comput. Vis., pp.           Journal of Computer Vision (IJCV), the Computer Vision and Image Under-
        83–103, 2003.                                                                 standing Journal (CVIU), and Medical Image Analysis (MedIA). In 2008, he
   [17] S. Lee, “Digital image smoothing and the sigma filter,” CVGIP, vol.            was the laureate of one of Greece’s highest honor for young academics and sci-
        24, no. 2, pp. 255–269, Nov. 1983.                                            entists, the Bodossaki Foundation Prize in the field of applied sciences, while in
   [18] M. Mahmoudi and G. Sapiro, “Fast image and video denoising via non-           2006, he received the TR35 prize from the MIT Technology Review journal.
        local means of similar neighborhoods,” IEEE Signal Process. Lett., vol.
        12, pp. 839–842, 2005.
   [19] S. Mallat, “A theory for multiresolution signal decomposition: The
        waveletrepresentation,” IEEE Trans. Pattern Anal. Mach. Intell., vol.                                    Frédéric Guichard graduated in mathematics from
        11, pp. 674–693, 1989.                                                                                   the Ecole Normale Supérieure, Paris, France, and is
   [20] P. Perona and J. Malik, “Scale space and edge detection using                                            also a member of the “Corps des Ponts et Chaussées”.
        anisotropic diffusion,” IEEE Trans. Pattern Anal. Mach. Intell., vol.                                    He received the Ph.D. degree in applied mathematics
        12, pp. 629–639, 1990.                                                                                   and image processing from the University Paris IX in
   [21] J. Polzehl and V. Spokoiny, “Adaptive weights smoothing with applica-                                    1994.
        tions to image restoration,” J. Roy. Statist. Soc. B, vol. 62, pp. 335–354,                                 He is currently Chief Scientist and CTO at DxO
        2000.                                                                                                    Labs, which he cofounded in 2003. He is responsible
   [22] L. Rudin, S. Osher, and E. Fatemi, “Nonlinear total variation based                                      for the image processing research and development
        noise removal,” Phys. D, vol. 60, pp. 259–268, 1992.                                                     for all products of the company. Prior to that, he oc-
   [23] B. Smolka and K. Wojciechowski, “Random walk approach to image                                           cupied several academic and industrial research po-
        enhancement,” Signal Process., vol. 81, pp. 465–482, 2001.                    sitions. His research interests include digital photography, image processing,
   [24] C. Tomasi and R. Manduchi, “Bilateral filtering for gray and color im-         computer vision, and applied math. He has published over 100 papers in inter-
        ages,” in Proc. Int. Conf. Computer Vision, 1998, pp. 839–846.                national journals, conferences, or patents.
   [25] D. Tschumperlé, “PDE’s Based Regularization of Multivalued Images                Dr. Guichard has received several awards, such as the “Prix Sciences et
        and Applications,” Ph.D. dissertation, Univ. Nice Sophia Antipolis,           Défense” in 1996 and a “grand prix” from the French Academy of Sciences in
        France, 2002.                                                                 2006.

To top