IMAGE DENOISING AND SEGMENTATION VIA NONLINEAR
DIFFUSION*
YUNMEI CHENt, BABA C. VEMURI 1, AND LI WANG
Abstract.
Image dnoising and segmentation are fundamental problems in the field of image processing and
computer vision with numerous i.,.;., .... In this paper we present a novel nonlinear diffusion
model augmented with reactive terms that yields quality denoising and segmentation results on a
variety of images. We present a proof for the existence, uniqueness and stability of the viscosity
solution of this PDEbased model. To achieve a faster implementation, we embed the the model in
a scale space and the solution is achieved via a dynamic system governed by a coupled system of
first order differential equations. The dynamic system finds the solution at a coarse scale and tracks
it continuously to a desired fine scale. We implement this scalespace tracking using a .... .il1 ;.
technique and demonstrate the smoothing and segmentation results on several images.
Key words. Nonlinear Diffusion, Image Processing, Segmentation, PDEs, Scalespace
1. Introduction. Image denoising and segmentation are fundamental problems
in the field of image processing and computer vision. Image denoising (or noise re
moval) is a technique that enhances images by reducing ;,I degradations that 11i ,
be present. The most common degradation source is the noise from the image acqui
sition  i. in and can be modeled as Gaussian random noise in most cases. Another
source of degradation is the so called saltandpepper noise that can occur due to a
random bit error in a communication channel during transmission. In the context
of image segmentation, the most important step is detection of region boundaries or
edges. An edge in an image !1 be defined as a location in the image at which a
significant change occurs in image i,, I i Segmented images or edge maps contain
very useful information and are used very often to  the essential content of an
image. Such image representations are useful in object recognition, low bitrate image
coding  1 in and various other applications [20].
The problem of image denoising and segmentation can be posed in either a deter
ministic or a stochastic framework. 4, .., I I,, methods are quite effective in achieving
segmentation but are limited by their intense computational requirements [12, 11].
We will therefore limit ourselves to the deterministic formulations, specifically, par
tial li!. t. i I equation (PDE) based methods that lend themselves to fast numerical
implementations.
Image denoising and segmentation can be formulated using variational principles
which in turn require solutions to PDEs. Recently, there has been a flurry of activity
in the PDEbased segmentation schemes. In [29], Perona and Malik developed an
anisotropic 1' T.,, .... scheme for image smoothing and segmentation. The basic idea
of this nonlinear smoothing scheme was to smooth the image while preserving the
edges in it. This was done by using the following equation It = div(c(VI)VI), where
I is the image to be smoothed and It describes its evolution over time, and c(VI)
is a decreasing function of VI. Segmentation was achieved by finding edges in this
*This work was supported by the Whitaker Foundation, the NIH under grant RO1LM05944 and
the National Science Foundation ERC under the grant EEC9402989 and its industrial partners.
tDepartment of Mathematics, University of Florida ( .i r.:' At r ie .
+Department of Computer and Information Science and Engineering, University of Florida
(vemuri~cise.uf edu).
Department of Computer and Information Science and Engineering, University of Florida
(liwangQcise.ufl.edu).
Y. CHEN, B. C. VEMURI AND L. WANG
smoothed image. Catte et al., [6], Nitzberg and Mumford [26] and Alvarez et al. [2]
recognized the illposedness of the PeronaMalik diffusion and proposed modifications
to overcome the same. Since then, several nonlinear diffusion methods have been
developed and a good account of these can be found in [27, 30, 33]. In [17], Kimia et
al. proposed an elegant reactiondiffusion based theory which describes the shapes of
objects in an entropy scalespace This theory was later used by Tek and Kimia [31]
for image segmentation applications.
Image segmentation can also be achieved by approaches based on curve evolution.
Malladi et al. [23, 24] and Caselles et al. [4] used curve evolution for recovering
shapes from 2D and 3D images. The curve evolution equation was implemented by
embedding the initial curve as a level curve in a surface and allowing all the level curves
of the surface to evolve simultaneously. This levelset method has the advantage of
being able to elegantly represent topological changes during the evolution of the curve
and thereby allows recovery of shapes without a priori knowledge of their I. 1. .1. 1 .
The evolution equation used in [4, 23, 24] was 08/8t = g(VI I '. I c + r), where 5
is the embedding surface for the curve evolution, g(VI) = 1/(1+ I(G, I ) is
a stopping term applied on the curve evolution, = div( ) is the curvature of
the level curves of 5, and c is a constant speed evolution term. This method was
generalized more recently in Caselles et al. [5] and Kichenassamy et al. [16] who also
established the link between the curve evolution based methods and the very popular
S1 ,I i i based snakes (active contour models) [15, 32, 25] used for segmentation in
computer vision and image processing literature. In [22], Malladi and Sethian propose
a unified approach to noise removal and image segmentation using the concept of min
max curvature flow. Based on the image data, a min/max switch was designed to
select min(r,, 0.0) or maz(r,, 0.0) so that the curvature based curve evolution smoothes
out small oscillations, but maintains the essential properties of the shape. Results of
implementation were shown on a i Ii of images yielding !,L ,ll noise removal
and image segmentation. Anisotropic diffusion till. that use a tensorial .lt!! IL, I1
parameter was introduced in Weickert [33]. These till, i can be tailored to enhance
image structures (edges, parallel lines, curves etc.) that occur in preferred directions.
More recently, Kimmel et. al., [18] presented a very general flow called the Beltrami
flow as a general framework for image smoothing and show that most flowbased
smoothing schemes 11 i be viewed as special cases in their framework.
In [30], ,! I! developed a common framework for curve evolution, image denoising
and segmentation, and anisotropic diffusion. In this work, a new segmentation func
tional was developed which lead to a coupled I 1 i i of PDEs, one of them performed
nonlinear smoothing of the input image and the other smoothed an ". l,. i !i.I I !i
function. I !, iI [30] demonstrated that all the existing curve evolution and anisotropic
diffusion schemes reported in literature can be viewed as special cases of his method.
Each of the methods we discussed above is in a variational form that minimizes
an i functional which in some cases is nonconvex. This nonconvex minimization
is hard and most computationally feasible methods lead to suboptimal solutions. In
[35, 34], a coarsetofine scale space tracking technique was proposed as a means
to tin , i 11i achieve a significant optimum for the nonconvex optimization. In this
approach, the desired solution is found by first finding a solution at a significantly
coarser scale and then tracking it down through finer scales (see 1 !i,.ii 1). It was
demonstrated that this technique can find significant minima of practical interest that
exist over a large range of scale [35, 34]. In our implementation, to achieve better
computational performance, we use a similar approach to l ;]
IMAGE DENOISING AND SEGMENTATION
E(%, )
FIG. 1. Scale space tracking: the desired solution is obtained by finding a solution at a coarser
scale and then tracking it down through finer scales.
In this paper, we present a new PDEbased image denoising and segmentation ap
proach which is based on a nonlinear [' li ti. ,I equation with additive reactive terms.
We prove the existence, uniqueness and l .1 ilil of the ;... i . ,Ii. i.i for this equa
tion and present implementation details along with several experiments demonstrating
the effectiveness of the proposed denoising and segmentation scheme.
The rest of the paper is organized as follows: in the next section, we present our
new nonlinear [d!Iti! i model which is used for selective smoothing of images. In
section 3, we prove the existence, uniqueness and l I1 1il of the ;... i solution
of the model. Section 4 contains a description of the numerical methods used to
implement the model equation. In section 5, examples of results obtained using our
new technique on a i;. of image data are presented. We conclude the paper in
section 6.
2. Nonlinear diffusion model. Numerous models of linear and nonlinear dif
fusion have been proposed in literature for achieving image smoothing and segmenta
tion. A survey of various nonlinear methods is discussed in Weickert [33]. The linear
models involved the standard heat equation dtu = Au with u(x,0) = f(x) where
f : R2  R is a scalar valued image. This I i" of linear [!It!..! blurs important
image features such as edges. In addition, it displaces the edges i.e., when moving
from finer to coarser scales, it dislocates the edges. In general, relating structures
across scales is nontrivial due to the presence of bifurcations.
Nonlinear anisotropic diffusion has been proposed by ii 1, ~ researchers [29, 6, 2, 7].
All of these are nonlinear models and [',ti, i in the diffusivity .... Itin Ii and/or the
diffusion term. Some of them are also supplemented with a reactive term. In the
following we present a nonlinear [' tli t !. I equation supplemented with reactive terms
for achieving edge preserving smoothing and segmentation.
Our model takes the following form:
Ou Vu
(1) = g(VG, u) Vuldiv( ) + Vg(VG, u) Vu
at Vu
(2) IVu(u I)
Ou
(3) IR = 0, o= I
On
Y. CHEN, B. C. VEMURI AND L. WANG
where I(x, y) is the !ii,. I 1i image to be processed, u(x, y, t) is it's smoothed version,
[3 is a weighting parameter, G,(x,y) = 'exp{ (x2 + y2) /(2 7r2)} is a Gaussian
smoothing kernel with a prespecified a, and f .) = 1/(1 I 12/K) is a nonincreasing
real valued function (for s > 0) which tends to zero as s ooc with a constant K.
This equation can be posed as the minimization of the following i. i functional
for fixed v, where v = VG, u. Note that the variation is taken with respect to u for
a fixed v.
(4) F (u) = f {g(v) u + 3(u )2 }ddy
From this functional, it is easy to see that:
1. the w.... i, i. of the first term namely, g(VG, u) serves the purpose of
selecting the locations in the image for smoothing. For instance, at image
locations having large values of gradient, this .... ti. I 1 takes on a small
value thereby reducing the smoothing performed at these locations since ,/ )
is a nonincreasing function of s.
2. The term Vu regularizes the solution u and is responsible for the nonlinear
term IVu div ( U) in equation 2. It diffuses u in a direction away from that
of Vu.
3. The term (u I)2 forces u to be a close approximation to the data I.
Note that our model is .[l!!. i! I! from the one reported in Alvarez et. al., [2]
in that it yields two additional terms in the EulerLagrange PDE 2 a.k.a. the
gradient descent equation for the variation principle 4. These two reactive terms are
responsible for forcing, the evolutionbased image smoothing to stop at the edges and
the resulting solution to be close approximation to the original image data respectively.
Our model is also .[l!!. i! I! from the one reported in Alvarez and Esclarin [1] where
in the reactive term was used to accommodate the image quantization application.
Our model is in spirit similar the smoothing equation in [30] however, the tili 1 11i
*.. tih w ii is very tl!, i, W I in our case.
For a fast implementation and to obtain ,!ti Ill smoothing and segmentation
results, the scale parameter a in the functional 4 can be varied from large (coarse) to
small (fine) values using a scale space tracking scheme. The scale space tracking can be
posed as the steady state solution of a differential equation obtained by t I i I ii i,
the equation representing the equilibrium condition of the functional 4 with respect to
a. The resultant differential equation describing the scale space trajectory is however
very complicated. A simpler coupled  , i of equations for tracking the solution
across scale (from coarsetofine a) for this nonlinear, L[I !,. !! model can be obtained
using the approach described in ;,] The coupled 1 11, of scale space equations are
given by
Ou
(5) o= Vr ,
VFt c_ cIVE.(u,a)
at
where cl and c2 are prespecified constants. For a set of initial no and go which is
far away from the solution, VF (u, a) is large while cle~C27E(u a) is small, thus
equation (5) is solved for u at a nearly constant scale until a solution at this scale
is found. Once a solution is obtained, it satisfies 0 = 0 and = c. Equation
(5) can then be used to track the solution down to finer scales. This technique is
IMAGE DENOISING AND SEGMENTATION
computationally tin i. i and yields !i'iL li' smoothing and segmentation results as is
demonstrated in our experiments. Note that the proof of existence, uniqueness and
.l ,.1;1 of a  ... i solution of our model presented in the following section does
not hold for the above described coupled  1. i of scalespace tracking equations.
3. Existence, Uniqueness and Stability. Since our model (3) is highly non
linear and degenerate, we need the notion of socalled i ... i solution (see [9]). In
this section we will prove the existence, uniqueness and l l.l;! for the ... ;
solution to the equation (3).
The wellposedness of the ..;i solution of the mean curvature flow ut 
Vu div( ) =) 0, in R" x R+ and for the generalized mean curvature flow ut 
Vu (div( ') + 7) = 0, in R" x R+, 7 E R, were studied by Evans and Spruck
[10] and by ('I!. !GigaGoto [8], respectively. In [2], the existence, uniqueness and
I 1.I;l; for the i ; i; solution of the following highly nonlinear'I w!! Li .i equation:
ut g(VG, u) Vu (div( ) = 0, in R" x R+, where g(i. = 1/1 + '!'." was
established. Vi i ; 1 for the level set form of the active contour model
ut g(VG, I) Vu (div( 7) +) + Vg. Vu 0, in R2 x R+, 7 e R, where I(x, y)
is the initial surface embedding the initial contour, was briefly discussed in [5] and
[16].
Our model has a similar structure, but more nonlinear terms than in the models
mentioned above. The proof of well posedness for our model is inspired by the idea
of [2] and the definition of the  ...i; solutions for our model follows the notion
established in [9]. We first prove the existence of a solution. Our solution is obtained
by ;, 1.I' i, limits to the approximate solutions of the penalized uniformly parabolic
equations. We then prove the uniqueness and the l I1.l; 1 of our solution. The proof
follows the ideas in [2]. However, since our model has more nonlinear factors or terms
than the models mentioned above, more careful and delicate estimates are required, es
pecially in getting the uniform L"norm estimate for the gradient of the approximate
solutions and in establishing the estimate Li' .. r] u v < C i.' (c=0r Iu v ,
where u and v are two i ... ; solutions of (3).
Our model equation (2) is in two dimensions, mathematically we can study this
problem for ndimensional cases. a, / and K are constants in (2), and they do not
effect the proof of wellposedness. To simplify the presentation, we shall consider a =
p = K = 1 and work with periodic boundary conditions similar to the presentation
in [2]. Then, by periodic extension we consider the following Cauchy problem
F= g(VG u)aij(V ,;,
(6) + (VG *u)[(VG,, *u) Vu] IVu(u I), x E R", t R+
u(x, 0) = I(), x e R"
where G exp{(x2+y2)/4}, g(,; = aij (~,' = 6ij and the summation
convention is used.
Sii , let us recall the definition of ...i ; subsolution of (6). A function u E
C([0, T] x R") for some T > 0 is said to be a i... solution of (6), if for all
0 E C2(R2 x R), the following condition holds at ;~ point (xo, to) E R" x (0, T], at
which (u 0) attains a local maximum.
(7) (Xo, to) ((VG u)(xo, to))aij (V(xo, to)).'. (xo, to)
+ a((VG u)(zo, to))[(VG', u)(xo,to) V(zxo, to)]
Y. CHEN, B. C. VEMURI AND L. WANG
V(xo, to) (u I) (zo, to) < 0, if VO(xo, to) 0
(8) (xo,to) g((VG *u)(xo, to))limsupaij (,I.'. (xo,to) < 0,
at p+o
if VO(xo,to)= 0
A ; ...i supersolution is similarly defined by substituting "I ... 11 I. i.!11i iin
for "local !iii!i,!woi!,! "< 0" for "> 0", and "I1!i.' for "liminf" in equations (7)
and (8) respectively. A i .. i solution is a continuous function which is both a
subsolution and a supersolution. We now state and prove the main theorem of this
paper.
THEOREM 3.1. 1!. Cauchy problem (6) has a unique viscosity solution u E
C(R" x [0, T]) n L"(0, T; WIV'(R")) for ....., T [0, oo), and infR I < u(x, t) <
supR, I, provided that I is Lipschitz continuous and continuous on R".
Moreover, if v E C(R" x R+) is a viscosity solution of (6) with I replaced by a
Lipschitz continuous function Ii, then for all T E [0,+oo), there exists a constant
C > 0, depending (..l., on I, I1, and T, such that
sup .,i t)v(x,li .(R^) I < [ Il IL(R" )
O
j ...... In this proof of the theorem, we use the technique of ...i solutions
theory discussed in [1, 2, 8, 9, 10]. The proof is in several stages.
Step 1. We first show that if u is a I... solution of (6) on R" x R+, then,
(9) inf I < u < sup I, on R" x [0, oo)
R" R^
Let = supR.n I + 5t (where 5 > 0) in (7) and assume that u 0 attains a
local maximum at (xo, to) with ,, > 0, then Vo(xo,to) = 0 and from equation (7),
(xo, to) < 0. This contradicts = 5 > 0 on R" x [0, oo). Therefore, u must
attain its maximum at to = 0. So,
u < sup(I sup I)
R"
u ^{
R"
Similarly we have (from the definition of supersolution)
u > inf I St.
R"
Leting 6 * 0 proves (9).
Step 2. Next we prove the gradient estimate for the approximate solution. Con
sider the following Cauchy problem:
a = gF(VG uF)a(. (Vu))u,
(10) + (VG *u)[(VG', *u'). Vu']
b"(VuF)(u" IF), x e R", t R+
u (x,0) = I(x), x e R"
IMAGE DENOISING AND SEGMENTATION
where
0
g (s) = ) +
aF(i,. = (E + 1)6ij Pi P32
b (, = I 12
IF E Co (Rn) (periodic) such that I  I uniformly and
I I .1 (R" ) < 'I 11 (R"),> I I'(R)< l .III (R" ).
From the theory of quasilinear uniformly parabolic equations ([19] 6, thm. 4.4),
the problem (10) admits a smooth solution uF E C"(R' x R+). Since ;:i smooth
solution is a ; ... i solution, by an argument similar to that in step 1, we know
that,
(11) u1 < M for (x,t) E R' x [0, oo)
where M > 0 is a constant depending only on I. Now we shall show a uniform
estimate for IVu6L(R ).
Differentiating (10) with respect to Xk yields,
at 01
(1 = g G VG ,., u?(VG u) (Vu"x <
+ 1(VG u)(G, u" .G (Vu",
02 1
(VG *ug ) ( *
F (VG u')(Gxn,* u')(^Gx, u' Vu')
S81 (VG u)[(VGox u) Vu] + 8(VG u')[(VG,, *u') Vu]
ob" (Vu')
am Xmes (
b (Vu')(u 1 k), x E R", t > 0.
, ill.li ;!d by 2u' on both sides of (12) and taking summation w.r.t. k, we get
qVu 2 82 1 V Vu I2
at axiaxj
g(VG u') (Vu,.:
01 ax,
,,V .b"( (Vu) F OVu" 2
(VG* u)(G *u u) ((VV" / )
01 Om ax,,
=2 (VG *u)(Gx, *u", (Vu' ,2 '
82 qF
+ 2 (VG u')(Gxmk u')(VGx, *u" Vu'.,
+2 2 (VG uF)[(VG,,,, uF) VU],,
2g(VG u',.,, (Vu '., k, 'U
2b (Vu')(Uk J X
Y. CHEN, B. C. VEMURI AND L. WANG
;From the definitions of a1., b6, g", and G, we have
a (Vu )u) 2 < 2a, (Vu')u' ui,
aI(VG u)2 < 2g6(VG u)
and for ;. i multiindex a with a < 2,
supR R+ JVaG u < C
supR xR+ IVg (VG u) < C,
where C > 0 is a constant depending only on M in (11).
Inserting these estimates into (13) and using Cauchy's ;!1,,l l i we have
(14) RHS of (13) < C(Vu 2 + 1) in R x R+
where C > 0 is a constant depending only on M in (11), hence C depends only on I.
A1! 1;, the maximum principle [3] to (14) yields for all t E [0, T] (for .!! T < oo),
17U (),1I Io,(J(l) < I171 I.(wR")
(15) <. I' I,.(R ) < CT,
where CT > 0 depends only on T, and I. This implies that
u'(x, t) u'(y,t)l < CT x y for yx, y E R" and Vt E [0, T].
By the same argument used in [9], we have
u (x,s) u'(x,t) < CTIt s 1, for Vx E R' and Vs,t E [0, T].
Then, by the AscoliArzela theorem, there exists a subsequence ui" of u", and a
function u E C(R1" x [0,1T]) n L (O, ; WI'IX(R")) such that as Ek 0,
(16) uS" u locally uniformly in R" x R+.
Step 3. Existence of .... i solution (16).
We assert now that u obtained in (16) is a i ..if solution of (6) in the sense
of (7) and (8).
Let EC C2 (R x R) and assume u 0 has a strict local maximum at a point
(xo, to) E R" x R+. As u" > u uniformly near (Xo, to), uS has a local maximum
at a point (xk, tk) with
(17) (Xk,tk) > (xo,to) as k > oc
and at (k, tk)
(18) Vu" = V Uki = 4i, aj (Vul',, < a1 (V '.V'
Therefore, (10) implies that at (Xk, tk),
(19) g" (VG 'u G( '..  (VG' u )[VGr *u V ]
+T b" (VO) (u" Ir) < 0
IMAGE DENOISING AND SEGMENTATION
If Vo(xo,to) # 0, from to (17), for ,Lm. i. ,II large k, VO(xk,tk) # 0. One .,
apply limits in (19) to obtain (recalling the definitions of a g", b6, and (16) (17)
.s/ ag
(,g g(VG u)aij(V '(VG u)[(VG, u) V] + b(V)(u I) 0
at (xo,to)
which is the same as (7).
If Vo(o, to) = 0, let
hk = V(Xk,tk)
I V(ZkA,tk) 2 2
then (19) reduces to
(21) 00 g(VG U )((Ek + 1)6i h.i
at Ij
S(VG u")[(VG,, u ). V ]
+ b" (V )(u Il6) < 0 at (xk,tk)
Since VO(xk,tk) > 0, Ek > 0 as k + o0, hence, bk (VO(Xk, tk)) * 0. Moreover,
becuase ,.' < 1, there is a subsequence of hk, also denoted by hk, such that as
k + oo, h'  h in R" x R with lhl < 1. A,1I, ;i_ limits to (21), we get
(22) 0 g(VG u)(Siy hih,' < 0 at (xo, to)
This is the same as (8). If u 0 has a local maximum, but not necessarily a strict
local maximum at (xo, to), we just need to repeat the argument above with O(x, t)
replaced by ((x, t) = O(x, t) + x xo 4 + (t to)4. Therefore, u is a subsolution of
(6). Similarly, we can show that u is a supersolution. Hence, u is a i ; ... i solution
of (6).
Step 4. Uniqueness.
We will now prove the uniqueness and I l.1 li for the i... i solution of (6).
This proof is based on Theorem 8.3 in [9]. Let u be a i ... ;i solution of (6) with
Lipschitz continuous initial data I and v be a i i solution of (6) with I replaced
by a Lipschitz continuous function I1. Let
w (x, y, t) = u(x, t) v (y, t) (465)1x j At, t E [0, T], x,y E R"
where 6 > 0 and A > 0 are constants to be determined later.
Claim: w(x, y, t) attains maximum at t = 0 for an arbitrary positive constant A.
Indeed, if w(x, y, t) attains its maximum at some point (xo, Yo, to) with to > 0, by
theorem 8.3 in [9], for each p > 0, there exist X and Y, (n x i'i  !i! r i;, matrices,
and a, 3 E R, such that
(23) a 3 = A,
(24) ( 0 )
Y. CHEN, B. C. VEMURI AND L. WANG
and
(25) a g((VG* u)(xo,to))aijS(1 zo I o 2(zo yo))X
((VG u)(xo,to))[(VG,, u)(xo,to) 51'o o 2o o)]
+ 6l o yo 3(u(xo,to) I(xo)) < 0
(26) 3 g((VG v)(yo, to))aij (51 o o 2(o jo))
((VG v)(yo, to))[(VG,, v)(yo,to) 51 Io o 2(xo o)]
+ 6 1o .,, (v(yo, to) I1 (o)) > 0
where
(27) A = (Aij),x and
Aij = 651 Xo Yo 26j + 261 (xo yo)i(xo yo)j,
here (xo yo)i stands for the i th component of xo yo.
We observe that xo 7 yo. Indeed, if xo = yo, then from (27), A = 0, hence X < 0
and Y > 0 from (24). Thus, (25) (26) leads to a < 0 and f > 0, which contradicts
a = > 0. We now choose
(28) p = 6izo yo 2
,From (27) and (24) after some algebra we have,
(X ( B B
(29) ( Y B B '
where
Bij = Xo  Y 26ij + 5(xo yo)i(xo yo)j, 1 < i,j < n.
Let
U= (VG*u)(xo,to), V= (VG*v)(yo,to),
D = (aij(S~1 zo *,,,"' zo Yo))l
and
Q g(U)D fg(U)g(V)D
( vg(U)g(V)D g(V)D
Noting that Q is a nonnegative symmetric matrix, from (29) we have
Q( X 5 ( B BB
Q O  26Q B B
Taking the trace we get
(30) g(U)DijXj g(V)DijYij < 265 l( U) g))2trace(DB)
< 461 (f ( ',,, yo2.
IMAGE DENOISING AND SEGMENTATION
Then, from (23) and (25) ('t i)
I = a A <
where
I = g(U)DiXij 
< 46( g(
g(V)Dij^j
 1r),, yo 2
II= [a(U)(VG, u)(xo, to)
( V)(VG, v)(yo,to)] 6 xo o Yo o)
III = 1 o .,. u(xo, to) v(yo, to)) + (I(zo) (yo))].
We now estimate (32)(34). 1 i! I
(.;) UV < (VG u)(xo,to) (VG *v)(xo,to)
+ (VG v)(xo, to) (VG v)(yo, to)
< C( sup u v + zo o)
R"x [0,T]
where constant C > 0 depends only on I and 11 (as evident from (9)) and the Lipschitz
constants for u and v. Then, by the mean value theorem and the fact that u and v
are Lipschitz continuous, we have
q(U) (V < CU V
R" x 0.T1
v + Xzo Yol).
Similarly,
< (U)(VG,, *u)(xo,to) (V)(VGx, *v)(yo,to)
01 01
a+ I(U)(VG, ))(o,to) (1) v)(yo,to)
+ I (V)(VGx, v)(zo, to) 5 (V)(VGx, v)(yo,to)
< C( sup u v + zo yol),
R x [O,T]
here we used (. ;,) and the facts that < 2 and v are Lipschitz continuous.
It is easy to see that
(38) u(xo, to) v(yo, to) + II(xo) I (yo)
(39) < u(xo, to) v(xo, to)\ + 1v(xo, to) v(yo, to)
(40) + II(xo) Ii(xo)I + Ii1(xo) II(Yo)
R"x [0,T]
The constants C > 0 in (.;,) (37) depend only on I, I1 and the Lipschitz constants
for u and v.
Y. CHEN, B. C. VEMURI AND L. WANG
Inserting ( ;.)(38) into (31) (34) yields
(41) A< C51[( sup u vl ,, yo+ Xo yo4
R"x [0,T]
+( sup u v)zxo yo3].
R x [0,T]
On the other hand, since (xo, Yo, to) is the maximum point of w(x, y, t),
u(Xo, to) v(xo, to) (46) ' o  yo 4  to > u(yo, to) v(yo, to) Ato.
This leads to
(42) xo Yo < (46L)3
where L is a Lipschitz constant for u in R" x [0,T]. Combining this with (41) yields
(43) A
R"x [0,T] R" x[0,T]
Rx [0,T] Rx [0,T]
where Co > 0 depends only on I, I1 and the Lipschitz constants of u and v. We now
set
(44) 5=L4( sup juv v)
R"x [0,T]
and from (43) we obtain
(45) A
R"x [0,T]
Let
(46) A = Co(L4 + L + 2) sup u v
This leads to a contradiction with (45). Therefore for the choice (46) of A, w(x, y, t)
attains its maximum at t = 0, which is our claim. Hence,
(47) u(x,t) v(y,t) (46) x At
< sup (u(x,0) v(y,0) (4)1x y4)
x,yER"
< sup (I(y) I (y) + I(x) (x) (45)1 )
x,yER"
< supI i + sup (I(x)I(y)(45)I x )
R" \xyL>o
R" \xyL>o
Noticing that sup,>o Lr (46)lr4) is achieved at r = (6L)1/3, and letting x = y in
(47), from (44) and (46) we get
3
(48) sup u v < sup I i + sup u v
Rx [0,T] R" 4 R" x[,T]
+Co(L +L + 2)T sup u v
R"x [O,T
IMAGE DENOISING AND SEGMENTATION
Therefore, there exists To > 0, ,tl!, !,' small (To < 4 ) such that
8C'o(L3 +L 3 +2)
from (48) we have
(49) sup Iu v < 8 sup I I I1 .
R [O,To] R"
For large t, by iteration, we esily obtain
sup u v < C(T) sup lI I
R" x [0,T] R"
This proves the uniqueness and I.1 ,l1 ; for u. D
4. Numerical implementation. The numerical implementation of the nonlin
ear diffusion equation (2) is based on the upwind finite 1'T ...... scheme developed
by Osher and Sethian [28] for curve evolution via level sets. Implementing the time
derivative and the diffusive term g(VG, I) Vu div( ~ ) presents no" '!!!il! i 
and is straight forward. u is approximated by forward 1i!!. i. i. 
u11+1[i,j] ,, [,j]
At
and the diffusive term is approximated using the usual central differences, with
VU U2U 2UUyxy + Uyy
Vuldiv( )V 2uu u
w Vu u2 +u2
What does require special care is the implementation of the data term (in fact, an infla
tion term i.e., constant speed of expansion) /3Vu ) and the ".1ii.., I term VgVu
[5]. The inflation/ballooning term permits the development of firstorder shocks, i.e.,
discontinuities in orientation of the boundary of a shape, where the derivative is not
defined. Thus, we have to approximate the spatial derivative using the upwind finite
difference scheme [28]. The [.1., 'I I I term permits the development of discontinuities
which indicate the presence of object boundaries. In this case, we can not use central
finite differences but have to use forward or backward finite 1!tl! i. i adaptively so
that their directions are always away from the discontinuities. Let
D+, [,j]= [ +l,j], [,j]
D , [,j]= ,* [,j] , [ 1,j]
DJ, [,]=, [,J] ,] [j1]
D , [,j]= , [,j 1], [, jl]
D,. [ j]= (, [ +,j] [ 1,])/2
D,., [j=(, [J + 1] [,j1])/2
then,
(Vg Vu)[i,j] = max(D., j], )D1u [i,j]
+ min(D [ j],0)D+., [ ,]
+ max(D ., j],0)D u"[i, j]
+ min(D., [ j],O)D+., [ j]
and
Vu = {(max(D1 u[i,j],0))2 + (mrin(D. [ ], 0))2
+ (max(D ., [ j],0))2 + (min(Df+u"[i,j],0))2}
Y. CHEN, B. C. VEMURI AND L. WANG
For a detailed discussion on this scheme, we refer the reader to Osher and Sethian
[28], and Malladi et al. [21].
The discrete equations for the scalespace tracking are realized using an explicit
Euler step in which time derivatives are approximated using forward ,itl, i, I . for
mulas in time with a step size At. This gives us the following discrete equations:
(50) at+1 = t (At)cl exp c . '. (u, a)
(51) ut+1 = u_ (At)VEt (uta).
These equations were implemented using a multigrid technique [35, 13, 34]. The
multigrid implementation speeds up the solution and leads to ii li; smoothing and
segmentation results.
5. Experimental results. In this section we present the application of our
model for smoothing and segmentation of several image data sets. To test the ef
fectiveness of this model, we first choose a !l!. Gr, image. I i,ii.. 2(a) shows the
noiseless image which contains a triangle, a rectangle, a circle and a thin ellipse. Fig
ure 2(b) is obtained by adding Gaussian noise to (a), and the signaltonoise ratio
(SNR: the ratio between the variance of the noisefree image and the variance of noise
[20]) is 1:2. 1 ,i ,i. 2(c)(e) shows the smoothing and segmentation results of the i1. .i
synthetic image using our nonlinear h!, LI!.t model, where (c) is the smoothed image
u, (d) shows the magnitude of g(Vu), and (e) depicts the edge map which is the local
minima in g(Vu) obtained using the nonminima suppression method [14]. The vari
ous parameter settings used to obtain these results were, / = 0.01, K = 200.0 (note
that t[K is the contrast parameter), number of iterations is 250, and the total com
puting time on an Ultra Sparc1 is 180 seconds. Note that our smoothing model has
preserved the edges with minimal rounding of the sharp corners. The local maxima
in g(Vu) yields a very good il ;i segmentation considering the amount of noise in
the input data. In this and all the following examples, the stopping criterion for our
nonlinear smoothing algorithm was a user specified tolerance (103) on the norm of
the h!lti. i. ii.. between two consecutive iterates of u.
The second example is the smoothing and segmentation results of the popular
Lenna image wherein we artificially added Gaussian noise (SNR = 3 : 1) to the
image. The original Lenna image is shown in figure 3(a), (b) is the i. . version of
(a) obtained by adding Gaussian noise, (c) depicts the smoothed image u, (d) and
(e) shows the magnitude of g(Vu) and the edge map. The parameter settings used
to obtain these results were, / = 0.01, K = 200.0, number of iterations is 120, and
the total computing time on an Ultra Sparc1 is 85 seconds. Once again, the results
of smoothing and segmentation are of high '.lL ,li and the later i!i be used as a
caricature of the original image leading to considerable image data compression.
I ,L.' i 4 presents the smoothing and segmentation of a CT chest scan. In this
figure, (a) is the CT chest scan, (b) depicts the smoothed image u, (c) shows the mag
nitude of g(Vu) and (d) is the edge map. The parameter settings used to obtain these
results were, 3 = 0.005, K = 200.0, number of iterations is 10, and the total compute
time on an Ultra Sparc1 is 15 seconds. We can see that some of the "!ip".i Ii ,"1
details which are unclear in the original image are enhanced in the smoothing and
segmentation results.
Each of the following examples in figures 5 7 presents four images: (a) is the
original image, (b) shows the smoothing result, and (c) and (d) are the segmentation
results from our nonlinear [!liti.,1 model. These images are infrared images for
IMAGE DENOISING AND SEGMENTATION
automatic target tracking applications and some of them are the Texas Instruments
(TI) high value target acquisition images.
6. Conclusion. We have proposed a new PDEbased image denoising and seg
mentation algorithm based on nonlinear 'tIl t,.!!I augmented by reactive terms. We
prove the existence, uniqueness and I 1i.1i of the i...i solution of this model.
For a fast implementation, we embed this model in a scale space and achieve scale
space tracking via a dynamical I. !1, of coupled t1! i ii. I1 1! equations. The scale
space tracking is implemented using a multigrid scheme. We present experiments
depicting the performance of our model on several image data sets. Our future efforts
will be focused on a 3D implementation of our model and testing on volume images
e.g., magnetic resonance (MR) and CT images.
Acknowledgments. Authors would like to thank Dr. P. Mergo for providing
the CT images, and Dr. J. N. \\ !. !! for providing the infrared images.
REFERENCES
[1] L. Alvarez and J. Esclarin. Image quantization using reactiondiffusion equations. SIAM J.
Appl. Math., 57(1):153175, 1997.
[2] L. Alvarez, PL. Lions, and JM. Morel. Image selective smoothing and edge detection by
nonlinear diffusion. ii. SIAM J. Numer. Anal., 29(3):845866, June 1992.
[3] H. Brezis. Analyse Fonctionnelle, Theorie et Applicatons. Masson, Paris, 1987.
[4] V. Caselles, F. Catte, T. Coll, and F. Dibos. A geometric method for active contours. Nu
merische Mathematik, 66, 1993.
[5] V. Caselles, R. Kimmel, and G. Sapiro. Geodesic active contours. In Fifth International
Conference on Computer Vision, 1995.
[6] F. Catte, P. L. Lions, J. M. Morel, and T. Coll. Image selective smoothing and edge detection
by nonlinear diffusion. SIAM Journal of Numerical Analysis, 29:182193, 1992.
[7] P. Charbonnier, L. BlancFeraud, G. Aubert, and M. Barlaud. Two deterministic : .1i i _ ..i i .
regularization algorithms for computed imaging,. In in Proc. of the IEEE Intl. Conf. on
Image Processing (Ih !! volume 2, pages 168172. IEEE Computer Society Press, 1994.
[8] YG. Chen, Y. Giga, and S. Goto. Uniqueness and existence of viscosity solutions of . .". .11; .
mean curvature flow equations. Journal of L' .' Geometry, 33:749786, 1991.
[9] M. G. Crandall, H. Ishii, and P. L. Lions. User's guide to viscosity solution of second order
partial differential equation. Bulletin of the American Mathematical Society, 27(1):167,
1992.
[10] L. C. Evans and J. Spruck. Motion of level sets by mean curvature. i. J. D : .' Geometry,
33:635681, 1991.
[11] D. Geiger and A. Yuille. A common framework for image segmentation,. International Journal
of Computer Vision,, 6:227243, 3 1991.
[12] S. Geman and D. Geman. Stochastic relaxation, gibbs distributions and bayesian restoration
of images,. IEEE Trans. on Pattern Analysis and Machine Intelligence,, 6:721741, 1984.
[13] W. Hackbusch. MultiGrid Methods and Applications,. SpringerVerlag,, Berlin, 1985,.
[14] R. Jain, R. Kasturi, and B. G. Schunck. Machine Vision. McGrawHill, Inc., 1995.
[15] M. Kass, A. Witkin, and D. I.! ..... i.. Snakes: Active contour models. Int. Journal of
Computer Vision, 1(4):321332, 1988.
[16] S. Kichenassamy, A. Kumar, P. Olver, A. Tannenbaum, and A. Yezzi. Gradient flows and
geometric active contour models. In Fifth International Conference on Computer Vision,
1995.
[17] B. B. Kimia, A. R. Tannenbaum, and S. W. Zucker. 1i .I" shocks, and deformations i: The
components of shape and the reactiondiffusion space. Technical Report LEMS105, LEMS,
Division of Engineering, Brown University, 1992.
[18] R. Kimmel, N. Sochen, and R. Malladi. Images as ... 1.1i;. maps and minimal sur
faces:movies, color and volumetric medical images,. In in Proc. of the IEEE Conf. on
Computer Vision and Pattern I. pages 350355, June 1997.
[19] O. A. Ladyzhenskaya, V. A. Solonnikov, and N. N. Uraltseva. Linear and Quasilinear Equations
of Parabolic American Mathematical Society, Providence, RI, 1968.
Y. CHEN, B. C. VEMURI AND L. WANG
[20] J. S. Lim. TwoDimensional .' and Image Processing. Prentice Hall, Englewood Cliffs,
NJ, 1990.
[21] R. Malladi, J. Sethian, and B. C. Vemuri. A fast level set based algorithm for topology
independent shape : ..i.. .i; J. Mathematical Imaging and Vision, 6:269289, 1996.
[22] R. Malladi and J. A. Sethian. A unified approach to noise removal, image enhancement and
shape recovery. IEEE Trans. on Image Processing, 5(11):15541568, 1996.
[23] R. Malladi, J. A. Sethian, and B. C. Vemuri. Evolutionary fronts for topology independent
shape : ..i.. .i;, and recovering. In Proc. ECCV, 1995.
[24] R. Malladi, J. A. Sethian, and B. C. Vemuri. II ... : 1...1. Ii with front propagation: A level
set approach. IEEE Trans. PAMI, 17(2):158175, 1995.
[25] T. Mclnerney and D. !. ..I...,i.. Deformable models in medical image analysis: A survey.
Medical Image Analysis, 1(2), 1996.
[26] M. ,1 i.. and T. Shiota. Nonlinear image iii. .... ii. edge and corner enhancement,. IEEE
Transactions on Pattern Analysis and Machine Intelligence, 14(8):826832, 1992.
[27] P. J. Olver, G. Sapiro, and A. Tannenbaum. Invariant geometric evolutions of surfaces and
volumetric smoothing. Preprint, 1995.
[28] S. J. Osher and J. A. Sethian. Fronts propagation with curvature dependent speed: Algorithms
based on hamiltonjacobi formulations. Journal of Computational 79:1249, 1988.
[29] P. Perona and J. Malik. Scalespace and edge detection using anisotropic diffusion. IEEE
Trans. PAMI, 12(7):629639, 1990.
[30] J. Shah. A common framework for curve evolution, segmentation and anisotropic diffusion. In
IEEE Conf. on Computer Vision and Pattern I. 1996.
[31] H. Tek and B. B. Kimia. Image segmentation by reactiondiffusion bubbles. In Fifth Interna
tional Conference on Computer Vision, 1995.
[32] D. !. i .....i.. A. Witkin, and M. Kass. Constraints on deformable models: Recovering 3d
shape and nonrigid motion. Artificial Intelligence, 36:91123, 1988.
[33] J. Weickert. A review of nonlinear diffusion 1,11i ..... In (Eds.) B. ter Haar Romney, L. Florack,
J. Koenderink, and M. i .. i editors, Scalespace theory in computer vision,, volume
1252, of Lecture Notes in Computer Science,, pages 328. SpringerVerlag, 1997.
[34] G. Whitten. Scale space tracking and deformable sheet models for computational vision. IEEE
Trans. Patt. Anal. Machine Intell., 15(7):697706, July 1993.
[35] A. Witkin, D. i .....i,,.. and M. Kass. .... 1 matching through scale space. Int. J. Comput.
Vision, 1:133144, 1987.
IMAGE DENOISING AND SEGMENTATION
(a) (b)
> (d
(d)
(e)
FIG. 2. Smoothing and segmentation results of a noisy synthetic image (256 x 256). (a) Noise
free image; (b) noisy image; (c) smoothed image u; (d) of g(Vu); (e) edge map.
Y. CHEN, B. C. VEMURI AND L. WANG
(a) (b)
(c) (d)
/
s .1 /2
xY
FIG. 3. Smoothing and segmentation results of the noisy Lenna image (256 x 256) using single
PDEbased model. (a) Original Lenna image; (b) Lenna image with Gaussian noise; (c) smoothed
image u; (d) of g(Vu); (e) edge map.
IMAGE DENOISING AND SEGMENTATION
(b) (c)
13 '
FIG. 4. CT chest scan smoothing and segmentation. (a) Original image (512 x 512); (b)
smoothed image u; (c) f .. ofg(Vu); (d) edge map.
Y. CHEN, B. C. VEMURI AND L. WANG
(b) (c)
(d)
FIG. 5. Smoothing and segmentation results of tanks image. (a) Original image (256 x 256);
(b) smoothed image u; (c) . of g(Vu); (d) edge map. Parameters: /f = 0.05, K = 200.0,
number of iterations is 20, and the total computing time on an Ultra Sparc 1 is 7 seconds.
IMAGE DENOISING AND SEGMENTATION
(b)
(c)
,(d)
(d)
FIG 6. Smoothing and segmentation results of one of the TI InfraRed HighValue Target
i. (IRHVTA) Data. (a) Original image (256 256); (b) smoothed image u; (c)
of g(Vu); (d) edge map. Parameters: /[ = 0.05, K = 100.0, number of iterations is 20, and the
total computing time on an Ultra Sparc1 is 7 seconds.
... .. ......~
~ Tr
rr
iinr, ~~
Y. CHEN, B. C. VEMURI AND L. WANG
"I""  ""
0,
i
r~hLZ~
,,
~
 L
(C)
 V ......
.(d) .
(d)
FIG. 7. Smoothing and segmentation results of one of TI InfraRed HighValue Target Acqui
sition (IRHVTA) Data. (a) Original image (256 x 256); (b) smoothed image u; (c) .. of
g(Vu); (d) edge map. Parameters: /3 = 0.05, K = 200.0, number of iterations is 10, and the total
computing time on an Ultra Sparc1 is 4 seconds.
}
