Image Sequence Analysis via Partial Differential Equations.
This article deals with the problem of restoring and motion segmenting noisy image sequences with a static background. Usually, motion segmentation and image restoration are considered separately in image sequence restoration. Moreover, motion segmentation is often noise sensitive. In this article, the motion segmentation and the image restoration parts are performed in a coupled way, allowing the motion segmentation part to positively influence the restoration part and vice-versa. This is the key of our approach that allows to deal simultaneously with the problem of restoration and motion segmentation. To this end, we propose a theoretically justified optimization problem that permits to take into account both requirements. The model is theoretically justified. Existence and unicity are proved in the space of bounded variations. A suitable numerical scheme based on half quadratic minimization is then proposed and its convergence and stability demonstrated. Experimental results obtained on noisy synthetic data and real images will illustrate the capabilities of this original and promising approach.
Introduction
Automatic image sequence restoration is clearly a
very important problem. Applications areas include
image surveillance, forensic image process-
ing, image compression, digital video broadcast-
ing, digital -lm restoration, medical image pro-
cessing, remote sensing example, the
recent work done within the European projects,
fully or in part, involved with this important
problem (Automated Restoration of
Film and Video Archives), NOBLESSE 2 (Nonlin-
Model-Based Analysis and Description of Images
for Multimedia Application), IMPROOFS 3
(IMage PROcessing Operations for Forensic Sup-
sequence restoration is tightly
coupled to motion segmentation. It requires to
extract moving objects in order to separately restore
the background and each moving region
along its particular motion trajectory. Most of
the work done mainly involves motion compensated
temporal -ltering techniques with appropriate
2D or 3D Wiener -lter for noise suppres-
sion, 2D/3D median -ltering or more appropriate
?morphological operators for removing impulsive
noise [16, 38, 39, 31, 27, 52, 19, 17]. However,
and due to the fact that image sequence restoration
is an emerging domain compared to 2D image
restoration, the literature is not so abundant
than the one related to the problem of restoring
just a single image. For example, numerous PDE
based algorithms have been recently proposed for
noise removal, 2D image enhancement and 2D image
restoration in real images with a particular
emphasis on preserving the grey level discontinuities
during the enhancement/restoration process.
These methods, which have been proved to be very
eOEcient, are based on evolving nonlinear partial
dioeerential equations (PDE's) (See the work of
Alvarez et al [4], Aubert et al. [8], Chambolle
Lions [21], Chan [14, 67] Cohen [23], Cottet &
Germain [24], Kornprobst & Deriche [44, 43, 42],
Morel [3, 51], Nordstr#m [54], Osher & Rudin [60],
Perona & Malik [58], Proesman et al. [59], Sapiro
et al. [20, 61, 62, 12, 63], Weickert [71, 72], You et
al. This methodology provides several
advantages. Firstly, we can justify on a theoretical
point of view the model, using the theory of viscosity
solutions or the calculus of variations. Sec-
ondly, it provides some suitable numerical schemes
for which convergence may be proved. Finally, it
permits to obtain results of high quality.
It is the aim of this article to consider the important
problem of image sequence restoration by
applying this PDE based methodology, which has
been proved to be very successful in anisotropically
restoring images. To our knowledge, few
litterature exists on the analysis of sequences of
images using Partial Dioeerential Equations. How-
ever, we mention that this methodology has been
previously used in the context of multiscale analysis
of movies (see the works of Guichard [35] and
Moisan [50]). In all this work, we will assume that
the background is static. We recall that the background
will be de-ned as the most often observed
part over the sequence. Our goal will be to obtain
the motion segmentation and the restored background

Therefore, considering the case of an image sequence
with some moving objects, we have to
consider both motion segmentation and image
restoration problems. Usually, these two problems
are treated separately in image sequence analysis.
However, it is clear that these two problems should
be treated simultaneously in order to achieve better
results. This is the key of our approach that
allows to deal simultaneously with the problem of
restoration and motion segmentation.
The organization of the article is as follows.
In Sect. 2, we make some precise recalls about
one of our previous approach for denoising a single
image [26, 8, 43]. The formalism and the methods
introduced will be very useful in the sequel.
Section 3 is then devoted to the presentation of
our new approach to deal with the case of noisy
images sequences. We formulate the problem into
an optimization problem.
The model is theoretically justi-ed in Sect. 4
we prove the existence and the unicity of the
solution to our problem in the space of bounded
variation functions.
A suitable algorithm is then proposed in Sect.
5 to approximate numerically the solution. We
prove its convergence and its stability.
We propose in Sect. 6 some experimental results
obtained on noisy synthetic and real data that will
illustrate the capabilities of this new approach.
We conclude in Sect. 7 by recalling the speci-
-cities of that work and giving the future developments

2. Restoring a single image
In Sect. 2.1, we recall a classical method in image
restoration formulated as a minimization problem
[26, 11, 8]. Section 2.2 presents a suitable
algorithm called the half quadratic minimization
which will also be used in the sequel.
2.1. A Classical Approach for Image Restoration
be a given noisy image de-ned for
2\Omega ae R 2 which corresponds to the
domain of the image. r: is the gradient operator.
We search for the restored image as the
solution of the following minimization problem :
I
Z\Omega
Z\Omega
OE(jrI j)dx
is the usual euclidian norm, ff r is a constant
and OE is a function still to be de-ned. Notice
that if we recognize the Tikhonov-
Arsenin regularization term [68]. How can we interpret
this minimization with this choice? In fact,
we search for the function I which will be simultaneously
close to the initial image N and smooth
(since we want the gradient as small as possible).
However, this method is well known to smooth
the image isotropically without preserving discontinuities
in intensity. The reason is that with
the quadratic function, gradients are too much
penalized. One solution to prevent the destruction
of discontinuities but allows for isotropically
smoothing uniform areas, is to change the above
quadratic term. This point have been widely discussed
[13, 53, 64, 66, 11, 8]. We refer to [26]
for a review. The key idea is that for low gra-
dients, isotropic smoothing is performed, and for
high gradient, smoothing is only applied in the
direction of the isophote and not across it. This
condition can be mathematically formalized if we
look at the Euler-Lagrange Equation (2), associated
to energy (1) :
Notice that Neumann conditions are imposed on
the boundaries. Let us concentrate on the divergence
term associated to the term 2 of (1). If
we note
the normal vector to j
1), we can show that [26] :
div
cj
I jj (3)
where I jj (respectively I - ) denotes the second order
derivate in the direction j (respectively -). It
is interesting to notice that most dioeusions operators
used for image restoration may also be decomposed
as the weighted sum of the second directional
derivatives I - and I jj . We refer to [41]
for more details. As for operator (3), if we want
a good restoration as described before, we would
like to have the following properties
lim
lim
But it is clear that the two conditions in (5) are
incompatible. So, we will only impose for high
gradients [26, 11, 8] :
Many functions OE have been proposed in the literature
that comply with the conditions (4) and
(see [26]). From now on, OE will be a convex
function with linear growth at in-nity which veri-
-es conditions (4) and (6). For instance, a possible
choice could be the hypersurface minimal function
In that case, existence and unicity of problem
(1) has recently been shown in the Sobolev space
2.2. The Half Quadratic Minimization
Solving directly the minimization problem (1) by
solving directly its Euler Lagrange equation (2),
is something hard because this equation is highly
non linear.
To overcome the diOEculty, the key idea is to
introduce a new functional which, although de-
-ned over an extended domain, has the same
minimum in I and can be manipulated with linear
algebraic methods. The method is based on
the half quadratic minimization theorem, inspired
from Geman and Reynolds [30]. The general idea
is that under some hypotheses on OE (mainly OE(
concave), we can write it as an in-mum :
d
where d will be called the dual variable associated
to x, and where /(\Delta) is a convex and decreasing
function. We refer to the Appendix A for more
details. We can verify that the function proposed
in (7) can be written as in (8). Consequently, the
problem (1) is now to -nd I and its dual variable
d I minimizing the functional F(I ; d I ) de-ned by :
Z
Z
\Omega (d I jrI
It is easy to check that for a -xed I , the functional
F is convex in d I and for a -xed d I , it is convex in
I . These properties are used to perform the algorithm
which consists in minimizing alternatively
in I and d I :
I
I
I
d n+1
d I
To perform each minimization, we simply solve the
Euler-Lagrange equations, which can be written as
I
I
d n+1
with discretized Neumann conditions at the
boundaries. Notice that (13) gives explicitly d n+1
I
while for (12), for a -xed d n
I , I n+1 is the solution
of a linear equation. After discretizing in space,
we have that (I n+1
(i;j)2\Omega is solution of a linear
system which is solved iteratively by the Gauss-Seidel
method for example. We refer to the Appendix
B for more details about the discretization
of the divergence operator. We also mention that
the convergence of the algorithm has been proved
[69].
3. Dealing with Noisy Images Sequences
denotes the noisy images sequence
for which the background is assumed to be static.
A simple moving object detector can be obtained
using a thresholding technique over the
inter-frame dioeerence between a so-called reference
image and the image being observed. Decisions
can be taken independently point by point
[73]. More complex approaches can also be used
[55, 57, 56, 1, 36, 45, 16, 38, 39, 31, 27, 52]. How-
ever, in our application, we are not just dealing
with a motion segmentation problem neither just
a restoration problem. In our case, the so-called
reference image is built at the same time while
observing the image sequence. Also, the motion
segmentation and the restoration are done in a
coupled way, allowing the motion segmentation
part to positively inAEuence the restoration part
and vice-versa. This is the key of our approach
that allows to deal simultaneously with the problem
of restoration and motion segmentation.
We -rst consider that the data is continuous in
time. This permits us to present the optimization
problem we want to study (Sect. 3.1). In Sect.
3.2, we rewrite the problem when the sequence is
given only by a -nite set of images. This leads to
the Problem 2.
3.1. An Optimization Problem
denotes the noisy images sequence
for which the background is assumed to be static.
Let us describe the unknown functions and what
we would like them ideally to be :
restored background,
the sequence which will indicate
the moving regions. Typically, we would like that
if the pixel belongs to a
moving object at time t, and 1 otherwise.
Our aim is to -nd a functional depending on
so that the minimizers
verify previous statements. We propose to solve
the following problem :
Problem 1. Let N(x 1 the given noisy image
sequence. We search for the restored background
and the motion segmented sequence
as the solution of the following
minimization problem :
Z\Omega
Z
Z\Omega
Z\Omega
c
Z
Z\Omega
where OE 1 and OE 2 are convex functions that comply
conditions (4) and (6) , and ff c ; ff r
c are positive
constants. We will specify later the spaces over
which the minimization runs.
Getting the minimum of the functional means that
we want each term to be small, having in mind the
phenomena of the compensations.
The term 3 is a regularization term. Notice
that the functions OE 1 ,OE 2 have been chosen as in
Sect. 2 so that discontinuities may be kept.
If we consider the term 2, this means that
we want the function C(x 1 to be close to
one. In our interpretation, this means that we
give a preference to the background. This is
physically correct since the background is visible
most of the time. However, if the data
too far from the supposed background
t, then the dioeerence
will be high, and to
compensate this value, the minimization process
will force to be zero. Therefore, the
function can be interpretated as a motion
detection function. Moreover, when searching
for we will not take into account
1). This
exactly means that B(x 1 will be the restored
image of the static background.
What about regularizing in time
the functions? As we can notice, the term 3
is a spatial smoothing term and we may suggest
to add some temporal smoothing for the sequence
C. However, there are two diOEculties to keep in
ffl the sequence has to be well sampled in time.
This temporal regularization term will have no
real interpretation in cases of images taken at
very large times (as it can be the case in video-
surveillance).
ffl In the same spirit as before, the discretization
of the regularization operator (in time) will
be hard because it will depend strongly on the
kind of movement in the sequence. In fact this
kind of regularization is, in some way, equivalent
to -nd the optical AEow, which we wanted to avoid.
We can think that this term could be usefull and
well discretized in a multiscale approach.j
3.2. The Temporal Discretized Problem
In fact, we have only a -nite set of images. Con-
sequently, we are going to rewrite the Problem 1,
taking into account that the sequence N(x 1
is represented during a -nite time by T images
noted
Problem 2. Let N be the noisy se-
quence. We search for B and C as the
solution of the following minimization problem :
Z\Omega
Z\Omega
Z\Omega
c
Z\Omega
Before going further, one may be interested in
the link between this method and the variational
method developed for image restoration in section
2. To this end, let us consider a sequence of the
same noisy image. More generally, we can consider
a sequence of the same static image corrupted with
dioeerent noises. If we admit the interpretation of
the functions C h , we will have C h j 1. After few
computations, (15) may be re-written :
Z\Omega
ff r
Z\Omega
Consequently, if we observe the energy (1) proposed
for the image restoration problem, we can
consider B as the restored version of the mean in
time of the sequence. Notice that if the sequence
is simply T times the same image, both methods
correspond exactly. Therefore, this model devoted
to sequences of images can be considered as a natural
extension of the previous one for single image
restoration.
Now that we have justi-ed the proposed model,
let us prove that it is mathematically well posed.
It is the purpose of the next section.
4. A Rigorously Justified Approach in
The Space of Bounded Variations
Section 4.1 presents the mathematical background
of our problem : the space of bounded variations
which is suitable to most problems in vision
[60, 22]. Roughly speaking, the idea is to generalize
the classical Sobolev space W
1;1(\Omega\Gamma so that
discontinuities along hypersurfaces may be consid-
ered. After having precisely speci-ed the problem
in Sect. 4.2, we -rst prove the existence of a solution
in a constrained space (See Sect. 4.3). Using
this result, we -nally prove the existence and the
unicity of a solution over the space in bounded
variations in Sect. 4.4.
4.1. The Space
In this section we only recall main notations and
de-nitions. We refer to [2, 28, 33, 29, 75] for the
complete theory.
Let\Omega be a bounded open set in R N , with
Lipschitz-regular boundary
@\Omega . We denote by L N
or dx the N-Lebesgue dimensional measure in R N
and by H ff the ff\Gammadimensional Hausdoroe measure.
We also set jEj = L N (E), the Lebesgue measure
of a measurable set E ae R N .
B(\Omega\Gamma denotes the
family of the Borel subsets
of\Omega . We will respectively
denote the strong, the weak and weak? convergences
in a space
V(\Omega\Gamma by \Gamma\Gamma\Gamma!
, \Gamma\Gamma\Gamma*
\Gamma\Gamma\Gamma*
Spaces of vector valued functions will be noted by
bold characters.
Working with images requires that the functions
that we consider can be discontinuous along
curves. This is impossible with classical Sobolev
spaces such as W 1;1
(\Omega\Gamma . This is why we need
to use the space of bounded variations (noted
Z\Omega
where C 1(\Omega\Gamma is the set of dioeerentiable functions
with compact support
in\Omega . We will note :
sup
aeZ\Omega
oe
If
and Du is the gradient in the sense
of distributions, then Du is a vector valued Radon
measure and
jDuj(\Omega\Gamma is the total variation of Du
on\Omega . The set of Radon measure is noted
The product topology of the strong topology of
for u and of the weak? topology of measures
for Du will be called the weak? topology of BV ,
and will be denoted by BV \Gamma w?.
\Gamma\Gamma\Gamma*
We recall that every bounded sequence in
admits a subsequence converging in BV \Gamma w?.
We de-ne the approximate upper limit u
and the approximate lower limit
lim
lim
is the ball of center x and radius
ae. We denote by S u the jump set, that is to say
the complement of the set of Lebesgue points, i.e.
the set of points x where dioeerent
namely
After choosing a normal n u pointing
toward the largest value of u, we recall the
following decompositions ([5] for more details):
where ru is the density of the absolutely continuous
part of Du with respect to the Lebesgue mea-
sure,
jSu is the Hausdoroe measure of dimension
restricted to the set S u and C u is the Cantor
part. We then recall the de-nition of a convex
function of measures. We refer to the works
of Gooeman-Serrin [34] and Demengel-Temam [25]
for more details. Let OE be convex and -nite on
R with linear growth at in-nity. Let OE 1 be the
asymptote (or recession) function :
then for u 2
using classical notations, we
de-ne Z
Z
Z
\Omega nSu
jC
Z
We -nally mention that this function is lower
semi-continuous for the BV \Gamma w?-topology.
4.2. Setting the problem
Let us recall the problem. Notice that the derivatives
will be now considered as distributional
derivatives. Consequently, the problem is to minimize
over
BV(\Omega\Gamma T+1 the functional E de-ned by
Z
Z\Omega
Z\Omega
c
Z\Omega
We recall that the regularization terms are interpretated
as convex functions of measures (see
(18)). The precise hypotheses on the functions
is an even and strictly convex
function, nondecreasing on R + and there
exist constants c ? 0 and b - 0 such that
As for the data (N h ) h=1::T , we will assume that :
and we will denote mN and MN the constants de-
-ned by :
where ess \Gamma inf (resp. ess \Gamma sup) is the essential
in-mum (resp. supremum).
4.3. Existence of a solution in a constrained
space
Let us consider the problem :
BV(\Omega\Gamma T+1 be a minimizing
sequence of E. Thanks to the property (20), one
may bound the derivatives of B and C h , and the
second term of E (see (19)) permits us to obtain
a bound for C h . However, nothing can be said
about the norm of B because of the product in
the -rst term of E (functions C h may be zero).
To overcome this diOEculty, let us introduce the
restricted space
E(\Omega\Gamma de-ned by :
Then, we have the following theorem :
Theorem 1. Given a sequence of images N h
verifying (22)-(23), the minimization problem :
where OE j verify (20)-(21), admits a solution in the
set
E(\Omega\Gamma .
Proof: The proof of this theorem is based on
classical arguments. As mentionned at the begining
of this section, the idea is to bound unifor-
maly a minimizing sequence, extract a converging
subsequence and pass to the limit. Notice that
working on this restricted space permits to obtain
a uniform bound for B. We refer to [10] for the
complete proof.
4.4. Existence and unicity of a solution over
The previous theorem establishes the existence
of a solution on a restricted space. However,
this result is not satisfying because working in a
constrained space is not easy to handle because
the optimality conditions are inequations and not
equations. In fact, even if these constraints are
natural (with regard to the interpretation of the
variables), we would like to avoid them. This is
the aim of Theorem 2 but we -rst need a preliminary
Lemma 1. Let u 2
verifying
hypotheses (20)-(21), and ' ff;fi the cut-ooe
function de-ned by :
ff if x - ff
Then we have :
Z
Z
\Omega OE(Du)
This Lemma is very intuitive, however we have
to deal with distributional derivatives and functions
of bounded variation. Consequently, we have
to deal with jump sets and Cantor parts. We refer
to the Appendix C where the complete proof
is sketched.
Using this Lemma, we can state the following
Theorem 2. Under hypotheses (20)-(21)and
(22)-(23), the minimization problem :
admits a solution in
where the constants mN ; MN are de-ned by (23),
then the solution is unique.
Proof: Existence is proved showing that the
minimization problem (26) over
E(\Omega\Gamma is equivalent
to the same problem posed over
that is to say without any constraint (this is a direct
consequence of Lemma 1). This remark will
permit us to prove the existence of a solution.
As for unicity, the diOEculty comes from the apparent
non convexity of the function :
with respect to all variables (Notice that it is convex
with respect to each variable). However, if ff C
is large enough, we prove that this functional is
convex over E which permits to conclude.
We refer to [10] for the complete proof.
This theorem is important since it permits
to consider the minimization problem over all
any constraint. On a numerical
point of view, this remark will be also important
since we will not have to handle with Lagrange
multipliers. We can also remark that the
condition (29) is in fact natural : it means that
the background must be suOEciently taken into account

5. The Minimization Algorithm
In the preceding section, we saw that there was a
unique solution in
BV(\Omega\Gamma T+1 of the minimization
problem (28). The aim of this section is to propose
a suitable algorithm to approximate numerically
this solution.
Before begining, we would like to insist on
the fact that working numerically with
BV(\Omega\Gamma is
something hard. Firstly, we cannot write Euler-Lagrange
equations. Anzellotti [7] proposes an extension
of Euler-Lagrange equation but they are
variational inequalities. In an image restoration
background, Vese [69] gives a characterisation of
the solution using a dual formulation. However,
both of them cannot be used, for the time being,
numerically.
Secondly, discretizing directly functions in
is still an opened question. For theses reasons, we
propose an algorithm with two steps :
- Section 5.1 : we de-ne a functional E ffl on
a more regular space. We show that the associated
minimization problem admits a unique solution
in W 1;2
that the functional E ffl \Gamma\Gammaconverges to E for the
refer to [32, 47] for more
details about the notion of \Gamma-convergence). Con-
sequently, converge for the
-strong topology to the unique solution of the
initial problem.
- Section 5.2 : For a -xed ffl, we are going to
construct a sequence (B
converging
to for the L 2 \Gammastrong topology.
It will be found as a minimizing sequence of an
extended functional. This part usually referenced
as the half quadratic minimization.
Consequently, we are able to construct a sequence
converging to the unique
minimum of the functional E for the L 2 \Gammastrong
topology. We will end this section by presenting in
section 5.3 the precise discretized algorithm. Its
stability will be proved using the -xed point theorem

5.1. A Quadratic Approximation
We -rst extend an idea developed in [21]. For a
function f having hypotheses (20)-(21), we de-ne
the odd function f ffl by :
We observe that for ffl ? 0, f ffl - f and for all t, we
f(t). Using this de-nition,
let us denote by OE 1;ffl and OE 2;ffl the two functions
associated to OE 1 and OE 2 . We then de-ne the function
Z\Omega
Z\Omega
Z\Omega
c
Z\Omega
Then, using same ideas than for Theorem 2,
we can prove that there exists a solution in
1;2(\Omega\Gamma T+1 of the problem :
where the constants mN ; MN are de-ned by (23),
then the solution is unique. We will denote by
the unique minimizer. We have
the following proposition :
Proposition 1. The sequence of functionals
\Gamma\Gammaconverges to the functional E for the
\Gammastrong topology as ffl goes to zero. The
sequence of the unique minimum of E ffl , noted
\Gammastrong to
the unique minimum of E.
Proof: By construction, the sequence E ffl is a
decreasing sequence converging pointwisely to the
functional e
de-ned by :
e
e
Thanks to [47] (proposition 5.7), we can deduce
that \Gamma-converges to the lower semi continuous
envelope of e
(for the L 2 T+1 \Gammastrong topology)
noted R(E). We then show that in fact
using some compacity results developed for instance
in [25, 15].
5.2. An extension using dual variables
be the unique minumum of
the functional E ffl over W
. For a -xed
ffl, our aim is to approximate it. To this end,
we need the result recalled in the Appendix A
and already used for the image restoration problem
(see Sect. 2.2) : let us apply Theorem 3 to
the functions OE 1;ffl and OE 2;ffl which ful-l desired hypotheses
(Typically OE i;ffl (
are concave). We will
denote by \Psi 1;ffl and \Psi 2;ffl the associated functions
\Psi. We then de-ne the functional E d
ffl de-ned over
(\Omega\Gamma3 \Theta W 1;2
Z
\Omega \Theta C 2
dx
Z
\Omega \Theta dBjrBj 2
dx
c
Z
\Omega \Theta
dx
where we have introduced the variables dB , dC1 ,.,
dCT associated to B, C 1 ,., CT respectively. To
minimize the functional E d
ffl , the idea is to minimize
successively with respect to each variable :
given the initial conditions (B
Ch ), we
iteratively solve the following system :
ffl (B; d n
d n+1
d n+1
Equalities (37)-(38) are written for
that the order of the minimization procedure
is not important for all the results presented be-
low. The way to obtain each variable like described
in (35) to (38) consists in solving the associated
Euler-Lagrange equations. As we will
see in section 5.3, the dual variables d n+1
(d n+1
are given explicitly, while B n+1 and
are solutions of linear systems. Any-
way, before going further, we need to know more
about the convergence of this algorithm : does it
converges and does (B
This is the purpose of the following
proposition
Proposition 2. Let (B
Ch ) be the
initial condition in W 1;2
. Then the sequence
de-ned by the system (35)-(36)-(37)-
(38) is convergent in L 2
\Gammastrong. More-
over, the sequence (B
\Gammastrong to the unique minimum of E ffl
in W
that is to say (B
Proof: The basis of the proof is to write
the variational optimality conditions associated to
each step and to pass to the limit into them. To
this end we needed some results about non-linear
elliptic equations [48, 49] and we used a trick of
Minty (see for instance [18, 21]). For more details,
we refer to [21, 9, 40] where such kind of ideas have
been developed.
5.3. The discretized algorithm
Let us write explicitely the equations that the
system (35)-(36)-(37)-(38) implies. Starting from
an initial estimate (B 0 ; d 0
Ch ), the equations
that will be solved are the following
b div(d n
d n+1
\Theta ff c
c div(d n
d n+1
As we said in the previous section, (40) and (42)
give explicitely the values of d n+1
B and d n+1
Ch while
are solutions of a linear system.
Once discretized using -nite dioeerences, the linear
system can be solved by a Gauss-Seidel method for
instance.
We next prove that the discretized algorithm
described by (39) to (42) is unconditionally stable.
Proposition 3.
Let\Omega d correspond to the dis-
cretization
of\Omega . Let E
be the space of discrete
functions (B;
where
c
Then, for a given (B
there exists a unique
(\Omega\Gamma such that (39)-
are satis-ed.
Remark that the bounds (43) and (44) can be
justi-ed if we consider the continuous case (see the
proof of the Theorem 2). As for condition (45), it
is also very natural if we admit the interpretations
of the variables C h : if this condition is false, this
would mean that the background is never seen at
some points which we refuse.
Proof: Let us sketch the proof. The -rst step is
to express the discretized equations (39) and (41).
Using equations (40) and (42) and the Appendix
B for the divergence terms (see the de-nition of
coeOEcients (p i+k;j+l
k;l C n+1
ff r
ff r
ff C
We then show that we have a contractive application
and conclude by applying the -xed point
theorem. We refer to [10] for the complete proof
which is mainly technical.
During this proof we needed to write explicitely
the discretized equations to be solved. We give
below a sum-up of the precise algorithm. Notice
that it is not necessary to compute explicitely the
dual variables because they are directly replaced
into the divergence operator.
1. /* Initializations (may be changed) */
2.
3. /* General loop */
4. for(It=0;It-ItNumber;It++) f
5. /* Minimizing in B */
corresponding to the divergence discretization for

Appendix

7. - Solve the linear system (47) by an iterative
method (Gauss-Seidel) to -nd B n+1
8. /* Minimizing in C h */
9.
corresponding to the divergence discretization for
(see

Appendix

11. - Solve the linear system (48) by an iterative
method (Gauss-Seidel) to -nd C n+1
12. g /* Loop on h */
13. g /* Loop on It */
To conclude this section, we will notice that if
ff r
0, the functions (C n+1
are in fact obtained
explicitly by :
As we can imagine, this case permits important
reduction of the computational cost since T linear
systems are replaced by T explicit expressions. We
will discuss in Sect. 6 if it is worth regularizing or
not the functions C h .
6. The Numerical Study
This section aims at showing quantitative and
qualitative results about this method. Synthetic
noisy sequences will be used to estimate rigorously
the capabilities of our approach. In all experi-
ments, we will -x the weights of dioeerents terms
to
and we will discuss
about the opportunity to choose a non zero coef-
-cient ff r
c . The purpose of Sect. 6.1 is the quality
of the restoration. The Sect. 6.2 is devoted to the
motion detection and its sensibility with respect
to noise. We will conclude in Sect. 6.3 by real
sequences.
6.1. About the Restoration
To estimate the quality of the restoration, we used
the noisy synthetic sequence presented in Fig. 4
(a)(b).

Figure

4 (c) is a representation of the
noisy background without the moving objects. We
mentioned the value of the Signal to Noise Ratio
(SNR) usually used in image restoration to quantify
the results quality. We refer to [43] for more
details. We recall that the higher the SNR is, the
best the quality is. Classically used to extract the
foreground from the background, the median (see
Fig. 4 (d)) appears to be ineOEcient. The average
in time of the sequence (see Fig 4 (e)), although
it permits a noise reduction, keeps the trace of the
moving objects. The Fig. 4 (f) is the result that
we obtained.
To conclude that section, let us mention that
we also tried the case ff r
that is to say we do
not regularized the functions C h . The resulting
SNR was 14, to be compared with 14.4 (ff r
c 6= 0).
This kind of results has been observed in all experiments
regularizing the functions C h does not
seem to inAEuence the quality of the restored back-
ground. Naturally if we are just interested to the
movement detection, this regularization may be
important. However, this point has to be better
investigated and more experimental results have
to be considered before to conclude.
6.2. The Sensitivity of Motion Detection With
Respect to Noise
In this section, we aim at showing the robustness
of our method with respect to noise. To this end,
we choose a synthetic sequence (see Fig. 5) where
a grey circle is translating from left to right in
front of a textured background.
To estimate the sensitivity of the algorithm, we
corrupted the sequence by gaussian noise of different
variance (from 5 to 50). We give in Fig. 1
the value of the SNR of the corrupted sequences
for each variance.

Figure

7 presents -ve typical results obtained
for dioeerent values of oe (oe=5,15,25,35,45). The
second one gives qualitative informations concerning
the quality of the restoration and the motion
detection. The criterion used to decide whether
a pixel belongs to the background or not is : if
then the pixel (i; j) of the image
number h belongs to the background. Other-
wise, it belongs to a moving object. The threshold
has been -xed to 0.25 in all experiments.
We can observe that when the SNR of the data
is more than 8 (corresponding to oe - 25), results
are particularly precise : The SNR of the background
is more than 20 (see Fig. 2) and the error
detections are less than 5 percent (see Fig. 3).
When the SNR of the data is less than 8, the motion
detection errors grow rapidly but the quality
of the restored background still remains correct.
See for instance the last row of Fig. 7 obtained for
the triangles on both sides are well recovered
(observe the strong noise in the sequence).
Finally, notice that same parameters (ff r
ff r
c ) have been used for all experiments. Generally
speaking, we remarked that the algorithm
performs well on a wide variety of sequences
with the same set of parameters (ff r
Fig. 1. Signal to Noise Ratio of the data as a function of the variance.
5 153025Fig. 2. SNR of the background as a function of the SNR of the data.
Fig. 3. dotted (resp. plain) line : percentage of bad detections for the moving regions (resp. static background) as a
(a) (b) (c) SNR=9.5
(d) SNR=5.7 (e) SNR=9.8 (f) SNR=14.4
Fig. 4. Results on a synthetic sequence (5 images) (a) Description of the sequence (-rst image) (b) Last image of the
sequence (c) The noisy background without any objects (d) Mediane (e) Average (f) Restored background (ff r
c 6=
Fig. 5. Three images of the initial synthetic sequence (35 images are available)
Fig. of the background as a function of the SNR of the data. Right : dotted (resp. plain) line : percentage
Fig. 7. Left : One image of the noisy sequence. Middle : The motion detection based on variable C h at the same
time. Right : The restored background B. From top to bottom : Results for dioeerent variances of the gaussian noise
(5,15,25,35,
6.3. Results on Real Sequences
Numerous real sequences have been tested using
this methodology. We will present some results
where the background of the scene is seen
most of the time. To be more precise, we mention
that some experiments have been done where
some people were hiding the background more
than sixty percent of the time. In that case, the
background found does not correspond to the real
static regions and takes into account some people
for the reconstruction. One possible way to avoid
this could be to add some a priori information of
the mouvement.
The -rst real sequence is presented in Fig. 8
(a)-(b). A small noise is introduced by the
camera and certainly by the hard weather condi-
tions. Notice the reAEections on the ground which
is frozen. We show in Fig. 8 (c) the average in
time of the sequence. The restored background is
shown in Fig. 8 (d). As we can see, it has been
very well found and enhanced. Figure 8 (e) is a
representation of the function C h (using a threshold
of 0.5) and we show in Fig 8 (f) the associated
dual variable dCh .
The second sequence is more noisy than the -rst
one. Its description is given in Fig. 9 (a). To
evaluate the quality of the restoration, we show a
close-up of the same region for one original image
(see Fig. 9 (b)), the average in time (see Fig. 9
(c)) and the restored background B (see Fig. 9
(d)). The detection of moving regions is displayed
in Fig. 9 (e). Notice that some sparse motion
have been detected at the right bottom and at
the left side of the two persons. They correspond
to the motion of a bush and the shadow of a tree
due to the wind.
The last sequence is taken from an highway (see
Fig. 10). We give two images (Fig. 10 (a) and
(b)) and the corresponding motion detection below
Finally, we show in
Fig. 10 (e) the restored background. Notice that
there is a black zone at the top of the road which
comes from the fact that there are always cars in
that region.
Notice that corresponding animations are available
in the Pierre Kornprobst's home page 4 .
7. Conclusion
We have presented in this article an original
coupled method for the problem of image sequence
restoration and motion segmentation. A theoretical
study in the space of bounded variations
showed us that the problem was well-posed. We
then proposed a convergent stable algorithm to
approximate the unique solution of the initial minimization
problem.
This original way to restore image sequence
has been proved to give very promising result.
A straightforward extension to color image sequences
has recently been developed. To complete
this work, several ideas are considered : use
the motion segmentation part to restore also the
moving regions, think about possible extensions
for non-static cameras. This is the object of our
current work.


Appendix

A
The Half Quadratic Minimization Theorem
This theorem has been inspired by Geman and
Reynolds [30] and proposed by Aubert [8].
Theorem 3. Let ' : [0; +1[! [0; +1[ be such
'(
is concave on ]0; +1[. (1)
Let L and M be de-ned as:
and
. Then, there exists a convex
and decreasing function
such that
where:
and
Moreover, for every -xed t - 0
the value d t for which the minimum is reached is
unique and given by:
In addition, we can give the expression of the
function \Psi with respect to OE. If we note
Fig. 8. Sweeden Sequence : (a) and (b) Description of the sequence (55 images available). Two people are walking from
top to bottom. This sequence is available from the web site http://www.ien.it/is/is.html. (c) The average over the
time. (d) The restored background B. (e) Function C h associated to the image (a) (a threshold of 0.5 has been used). (f)
The dual variable d C h
associated to the image (a).
Fig. 9. INRIA Sequence : (a) Description of the sequence (12 images available). (b) Zoom on a upper right part of
the original sequence (without objects). (c) Zoom on the mean image. (d) Zoom on the restored background B. (e) The
function C h thresholded. (f) The dual variable d C h .
Fig. 10. Highway Sequence : (a) and (b) Two images from the sequence (90 images available). (c) and (d) Corresponding
C h functions. (e) The restored background.
However, notice that this expression will never be
used explicitely.


Appendix

On Discretizing the Divergence Operator
Let d and A given at nodes (i; j). The problem
is to get an approximation of div(drA) at the
node (i; j). We denote by ffi x1 and ffi x2 the -nite
dioeerence operators de-ned by :
1Using that notation, Perona and Malik [58] proposed
the following approximation :
@
d
@A
@
d
@A
-@
(1)
where the symbol ? denotes the convolution and
S P is the sum of the four weights in the principal
directions. Notice that we need to estimate the
function d at intermediate nodes. Our aim is to
extend this approximation so that we could take
into account the values of A at the diagonal nodes
A ? A i;j (2)
where ff P and ff D are two weights to be discussed,
and S D is the sum of the four weights in the diagonal
directions. Approximation (2) is consistent if
and only if :
there remains one degree of freedom. Two
possibilities have been considered :
functions of d (See Fig. 10) (5)
a
Fig. B.1. ff
is the direction of the gradient of d. Notice that ff D can be
deduced from the consistency condition is then computed
thanks to the consistency condition.
Before going further, remark that any kind of
discretization leads to :
div
rA
and
To compare these dioeerent discretizations, we
made numerical experiments with the image
restoration problem where such kind of operator
have to be discretized. We recall that for a given
I , we need to -nd I n+1 such that :
I
I
We refer to section 2 for more details. The value
of d n
I
at intermediate nodes is computed
by interpolation (see [58]).
We tested these dioeerent discretizations on a
noisy test image using quantitative measures. We
checked that (2) permits to restore identically
edges in principal or diagonal directions. More-
over, we observed that choosing ff P adaptatively
gave more precise results than (4). We used
this approximation (5) in our experiments.

Appendix

Proof of Lemma 1
Proof: Let us -rst recall the Lebesgue decomposition
of the measure OE(Du) :
Z
Z
\Omega OE(jruj)dx
Z
Z
\Omega =Su
We are going to show that cutting the fonction
u using the fonction ' ff;fi permits to reduce each
term. To simplify notations, we will sometimes
use the notation -
u for the troncated function
fig
=\Omega =\Omega c . Thanks to [37], we have
Z
Z
OE(jruj)dx. Consequently :
Z
Z
Z\Omega
c
dx
Z\Omega
OE(jruj)dx (1)
2: using results proved in [5], we know
that :
Thanks to these results, and since ' ff;fi is Lipschitz
continuous with a constant equals to 1, we have :
Z
Z
Z
3: we need to understand how is the Cantor
part of the distributional derivative of the composed
function ' ff;fi (u). Vol'pert [70] -rst proposed
a chain rule formula for functions
for
BV(\Omega\Gamma and when ' is continuously dif-
ferentiable. Ambrosio and Dal Maso [6] gave extended
results for functions ' uniformely Lipschitz
continuous. Since u is scalar, it is demonstrated
in [6] that we can
on\Omega =S u
where ~
u is the approximate limit of u de-ned by :
lim
r \GammaN
Z
where B(x; r) is the closed ball with center x and
radius r. Moreover, we have :
Z\Omega
Z\Omega
Z
Notice that the second integral equals to zero because
the Hausdoroe dimension of the set S u =S - u
is at most N \Gamma 1 and we know that for any
BV(\Omega\Gamma and any set S of Hausdoroe dimension
at most N \Gamma 1, we have C v
using the chain rule formula (3), we have :
Z\Omega
Z\Omega
jC
Z\Omega
Finally, using results (1), (2), (5) permits to write
Z\Omega
Z\Omega
This concludes the proof.
Notes
1. http://www.ina.fr/INA/Recherche/Aurora/index.en.html
2. http://www.spd.eee.strath.ac.uk/users/harve/noblesse.html
3. http://www.esat.kuleuven.ac.be/ konijn/improofs.html
4. http://www.inria.fr/robotvis/personnel/pkornp/pkornp-eng.html



--R

Bayesian algorithms for adaptive change detection in image sequences using markov random
Analysis of bounded variation penalty methods for ill-posed p roblems
Image selective smoothing and edge detection by nonlinear dioeusion (ii).
Signal and image restoration using shock
A compactness theorem for a new class of functions of bounded variation.
A general chain rule for distributional derivatives.
The euler equation for functionals with linear growth.
Deterministic edge-preserving regularization in computed imaging
A mathematical study of the regularized optical AEow problem in the space BV(
A variational method and its mathematical study in image sequence analysis.
A variational method in image recovery.
Robust anisotropic dioeusion.
Visual Recon- struction
Color tv: Total variation methods for restoration of vector-valued images

Noise reduction of image sequences using adaptative motion compensated frame averaging.
Simultaneous recursive displacement estimation and restoration of noisy-blurred image sequences

Deterioration detection for digital

Image recovery via total variation minimization and related problems.
Nonlinear variational method for optical AEow computation.
Auxiliary variables and two-step iterative algorithms in computer vision problems
Image processing through reaction combined with nonlinear dioeusion.
Convex functions of a measure and applications.
"Traitement du Signal"
Noise reduction in image sequences using motion-compensated temporal -ltering
Measure Theory and Fine Properties of Functions.
Geometric Measure Theory.
Constrained restoration and the recovery of discontinuities.
McClure, and Donald Ge- man
un tipo di conver- genza variazionale
Minimal Surfaces and Functions of Bounded Variation.
Sublinear functions of measures and variational integrals.
Axiomatisation des analyses multi- #chelles d'images et de -lms
Moving object segmentation based on adaptive reference images.
An Introduction to Variational Inequalities and Their Applications.
Reconstruction of severely degraded image sequences.
A system for re-construction of missing data in image sequences using sampled 3d ar models and mrf motion priors

Image restoration via PDE's.
Image cou- pling
Nonlinear operators in image restoration.

Motion detection in spatio-temporal space
Image processing: Flows under min/max curvature and mean curvature.
An introduction to

Some results on regularity for solutions of non-linear elliptic systems and quasi-regular functions
Traitement num
Segmentation of images by variational methods: A constructive approach.
Image Sequence Restoration using Gibbs Distributions.
Optimal approximations by piecewise smooth functions and associated variational problems.


IEEE Computer Society Press.
Detection and localization of moving objects in image sequences.
Detecting multiple moving targets using deformable con- tours

Coupled Geometry-Driven Dioeusion Equations for Low-Level Vision
Total variation based image restoration with free local constraints.
Contrast enhancement via image evolution AEows.
Anisotropic dioeusion of multivalued images with applications to color
Experiments on geometric image enhancement.
Unique reconstruction of piecewise-smooth images by minimizing strictly convex non-quadratic functionals
A common framework for curve evolu- tion
Discontinuity preserving regularization of inverse visual problems.
Spatially and scale adaptive total variation based regularization and anisotropic dioeusion in image processing.
Solutions of Ill-posed Problems

The spaces BV and quasilinear equa- tions
Anisotropic Dioeusion in Image Process- ing
Anisotropic Dioeusion in Image Process- ing
Motion detection from image informa- tion
Analysis and Design of Anisotropic Dioeusion for Image Processing.
Weakly Dioeerentiable Functions.
--TR
Visual reconstruction
Motion detection in spatio-temporal space
Weakly differentiable functions
Scale-Space and Edge Detection Using Anisotropic Diffusion
Biased anisotropic diffusion
Image selective smoothing and edge detection by nonlinear diffusion. II
Constrained Restoration and the Recovery of Discontinuities
A nonlinear filter for film restoration and other problems in image processing
Signal and image restoration using shock filters and anisotropic diffusion
Image processing
A Variational Method in Image Recovery
Contrast enhancement via image evolution flows
A System for Reconstruction of Missing Data in Image Sequences Using Sampled 3D AR Models and MRF Motion Priors
Reconstruction of Severely Degraded Image Sequences
Deterioration detection for digital film restoration
Non-linear operators in image restoration
A PDE-Based Level-Set Approach for Detection and Tracking of Moving Objects

--CTR
Jong Bae Kim , Hang Joon Kim, GA-based image restoration by isophote constraint optimization, EURASIP Journal on Applied Signal Processing, v.2003 n.1, p.238-243, January
Jong Bae Kim , Hang Joon Kim, Region removal and restoration using a genetic algorithm with isophote constraint, Pattern Recognition Letters, v.24 n.9-10, p.1303-1316, 01 June
A. Ben Hamza , Hamid Krim , Josiane Zerubia, A nonlinear entropic variational model for image filtering, EURASIP Journal on Applied Signal Processing, v.2004 n.1, p.2408-2422, 1 January 2004
Etienne Mmin , Patrick Prez, Hierarchical Estimation and Segmentation of Dense Motion Fields, International Journal of Computer Vision, v.46 n.2, p.129-155, February 2002
Thomas Corpetti , tienne Mmin , Patrick Prez, Dense Estimation of Fluid Flows, IEEE Transactions on Pattern Analysis and Machine Intelligence, v.24 n.3, p.365-380, March 2002
Stphanie Jehan-Besson , Michel Barlaud , Gilles Aubert, A 3-step algorithm using region-based active contours for video objects detection, EURASIP Journal on Applied Signal Processing, v.2002 n.1, p.572-581, January 2002
Daniel Cremers , Stefano Soatto, Motion Competition: A Variational Approach to Piecewise Parametric Motion Segmentation, International Journal of Computer Vision, v.62 n.3, p.249-265, May 2005
discontinuity preserving regularization;space of functions of bounded variation;variational approaches;sequence image restoration;motion segmentation