Explicit Algorithms for a New Time Dependent Model Based on Level Set Motion for Nonlinear Deblurring and Noise Removal.
In this paper we formulate a time dependent model to approximate the solution to the nonlinear total variation optimization problem for deblurring and noise removal introduced by Rudin and Osher [ Total variation based image restoration with free local constraints, in Proceedings IEEE Internat. Conf. Imag. Proc., IEEE Press, Piscataway, NJ, (1994), pp. 31--35] and Rudin, Osher, and Fatemi [ Phys. D, 60 (1992), pp. 259--268], respectively. Our model is based on level set motion whose steady state is quickly reached by means of an explicit procedure based on Roe's scheme [ J. Comput. Phys.,  43 (1981), pp. 357--372], used in fluid dynamics. We show numerical evidence of the speed of resolution and stability of this simple explicit procedure in some representative 1D and 2D numerical examples.
Introduction
. The classical algorithms for image deblurring and/or denoising
have been mainly based on least squares, Fourier series and other L 2 -norm approxi-
mations, and, consequently, their outputs may be contaminated by Gibbs' phenomena
and do not approximate well images containing edges. Their computational advantage
comes from the fact that they are linear, thus fast solvers are widely available. How-
ever, the effect of the restoration is not local in spatial scale. Other bases of orthogonal
functions have been introduced in order to get rid of those problems, e.g., compactly
supported wavelets. However, Gibbs' phenomenon, (ringing), is still present for these
norms.
The Total Variation (TV) deblurring and denoising models are based on a variational
problem with constraints using the total variation norm as a nonlinear non-differentiable
functional. The formulation of these models was first given by Rudin,
Osher and Fatemi in ([19]) for the denoising model and Rudin and Osher in ([18]) for
the denoising and deblurring case. The main advantage is that their solutions preserve
edges very well, but there are computational difficulties. Indeed, in spite of the fact
that the variational problem is convex, the Euler-Lagrange equations are nonlinear
and ill-conditioned. Linear semi-implicit fixed-point procedures devised by Vogel and
Oman, (see [26]), and interior-point primal-dual implicit quadratic methods by Chan,
Golub and Mulet, (see [6]), were introduced to solve the models. Those methods give
good results when treating pure denoising problems, but the methods become highly
ill-conditioned for the deblurring and denoising case where the computational cost is
very high and parameter dependent. Furthermore, those methods also suffer from the
undesirable staircase effect, namely the transformation of smooth regions (ramps) into
piecewise constant regions (stairs).
In this paper we present a very simple time dependent model constructed by evolving
the Euler-Lagrange equation of the Rudin-Osher optimization problem, multiplied
by the magnitude of the gradient of the solution. The two main analytic features of
y Department of Mathematics, University of California, Los Angeles, 405 Hilgard Av-
enue, Los Angeles, CA 90095-1555 and Departament de Matem'atica Aplicada, Universitat de
Dr. Moliner, 50, 46100 Burjassot, Spain. E-mail addresses: marquina@uv.es, URL:
http://gata.uv.es/~marquina. Supported by NSF Grant INT9602089 and DGICYT Grant PB97-
1402.
Department of Mathematics, University of California, Los Angeles, 405 Hilgard Avenue, Los
Angeles, CA 90095-1555. E-mail address: sjo@math.ucla.edu. Supported by NSF Grant DMS
this formulation are the following: 1) the level contours of the image move quickly
to the steady solution and 2) the presence of the gradient numerically regularizes
the mean curvature term in a way that preserves and enhances edges and kills noise
through the nonlinear diffusion acting on small scales. We use the entropy-violating
Roe scheme, ([16]) for the convective term and central differencing for the regularized
mean curvature diffusion term. This makes a very simple, stable, explicit procedure,
computationally competitive compared with other semi-implicit or implicit procedures.
We show numerical evidence of the power of resolution and stability of this explicit
procedure in some representative 1D and 2D numerical examples, consisting of noisy
and blurred signals and images, (we use Gaussian white noise and Gausssian blur).
We have observed in our experiments that our algorithm shows a substantially reduced
staircase effect.
2. Deblurring and Denoising. A recording device or a camera would record
a signal or image so that 1) the recorded intensity of a small region is related to the
true intensities of a neighborhood of the pixel, through a degradation process usually
called blurring and 2) the recorded intensities are contaminated by random noise.
To fix our ideas we restrict the discussion to R 2 . An image can be interpreted as
either a real function defined on \Omega\Gamma a bounded and open domain of R 2 , (for simplicity
we will
assume\Omega to be the unit square henceforth) or as a suitable discretization of
this continuous image. Our interest is to restore an image which is contaminated with
noise and blur in such a way that the process should recover the edges of the image.
Let us denote by u 0 the observed image and u the real image. A model of blurring
comes from the degradation of u through some kind of averaging. Indeed, u may be
blurred through the application of a kernel: k(x; s; r) by means of
Z
\Omega
u(s; r) k(x; s; ds dr (2.1)
and, we denote this operation by v u. The model of degradation we assume is
where n is Gaussian white noise, i.e., the values n i of n at the pixels i are independent
random variables, each with a Gaussian distribution of zero mean and variance oe 2 .
If the kernel k is translation invariant, i.e., there is a function j(x; y), (also called
a kernel), such that k(x; s; and the blurring is defined as a
'superposition' of j 0 s:
Z
\Omega
u(s; r) ds dr (2.3)
and this isotropic blurring is called convolution. Otherwise, if the kernel k is not
translation-invariant we call this blurring anisotropic. For the sake of simplicity, we
suppose that the blurring is coming from a convolution, through a kernel function j
such that j   u is a selfadjoint compact integral operator. Typically, j has the following
properties, goes to 1 and R
For any ff ? 0 the so-called heat kernel, defined as
is an important example that we will use in our numerical experiments.
The main advantage of the convolution is that if we take the Fourier transform of
(2.3) we get
then, to solve the model (2.2) with we take Fourier transform and we arrive at
To recover u(x; y), we need to deconvolve, i.e., this means that we have to divide the
equation (2.6) by - j(k; l) and to apply the inverse Fourier transform. This procedure
is generally very ill-posed. Indeed, j is usually smooth and j(x; y) ! 0 rapidly
as goes to 1, thus large frequencies in u 0 get amplified considerably.
The function u 0 is generally piecewise smooth with jumps in the function values and
derivatives; thus the Fourier method approximation gives global error estimates of
order O(h), (see ([11])) and suffers from Gibbs' phenomenon. Discrete direct methods
dealing with the linear integral equation (2.6) have been designed by different authors,
(see ([13] and references therein).
One way to make life easier is to consider a variational formulation of the model
that regularizes the problem. Our objective is to estimate u from statistics of the noise,
blur and some a priori knowledge of the image (smoothness, existence of edges). This
knowledge is incorporated into the formulation by using a functional R that measures
the quality of the image u, in the sense that smaller values of R(u) correspond to better
images. The process, in other words, consists in the choice of the best quality image
among those matching the constraints imposed by the statistics of the noise together
with the blur induced by j.
The usual approach consists in solving the following constrained optimization
problem:
min
subject to jjj
since E(
denotes the expectation of the
random variable X) imply that jjj
R\Omega (j
Examples of regularization functionals that can be found in the literature are,
r is the gradient and \Delta is the Laplacian, see Refs. [22,
8]. The main disadvantage of using these functionals is that they do not allow discontinuities
in the solution, therefore the edges can not be satisfactorily recovered.
In [19], the Total Variation norm or TV-norm is proposed as a regularization
functional for the image restoration problem:
Z
\Omega
Z
\Omega
y dx: (2.8)
The norm does not penalize discontinuities in u, and thus allows us to recover
the edges of the original image. There are other functionals with similar properties
introduced in the literature for different purposes, (see for instance, [7, 5, 25, 2]). The
restoration problem can be thus written as
min
Z
\Omega
subject to 1i Z
\Omega
(j
Its Lagrangian is
Z
\Omega
\Omega
(j
and its Euler-Lagrange equations, with homogeneous Neumann boundary conditions
for u, are:
\Omega
(j
There are known techniques, (see [3]), for solving the constrained optimization
problem (2.9) by exploiting solvers for the corresponding unconstrained problem,
whose Euler-Lagrange equations are (2.11) for - fixed. Therefore, for the sake of
clarity, we will assume the Lagrange multiplier - to be known throughout the exposi-
tion. For
, we can then write the equivalent unconstrained problem as
min
Z
\Omega
(j
and its Euler-Lagrange equation in the more usual form:
We call (2.14) the nonlinear deconvolution model. The linear deconvolution model
would be
that comes from the Euler-Lagrange equation of the corresponding unconstrained problem
with the norm
Since the equation (2.14) is not well defined at points where due to the
presence of the term 1=jruj, it is common to slightly perturb the Total Variation
functional to become:
Z
\Omega
where fi is a small positive parameter, or,
Z
\Omega
with the notation
3. The time dependent model. Vogel and Oman and Chan, Golub and Mulet
devised direct methods to approximate the solution to the Euler-Lagrange equation
with an a priori estimate of the Lagrange multiplier and homogeneous Neumann
boundary conditions. Those methods work well for denoising problems but the
removal of blur becomes very ill-conditioned with user-dependent choice of parame-
ters. However, stable explicit schemes are preferable when the steady state is quickly
reached because the choice of parameters is almost user-independent. Moreover, the
programming for our algorithm is quite simple compared to the implicit inversions
needed in the above mentioned methods.
Usually, time dependent approximations to the ill-conditioned Euler-Lagrange
equation are inefficient because the steady state is reached with a very small time
step, when an explicit scheme is used. This is the case with the following formulation
due to Rudin, Osher and Fatemi (see [19]) and Rudin and Osher (see [18]):
with u(x; given as initial data, (we have used as initial guess the original blurry
and noisy image u 0 ) and homogeneous Neumann boundary conditions, i.e., @u
@n
the boundary of the domain. As t increases, we approach to a restored version of our
image, and the effect of the evolution should be edge detection and enhancement and
smoothing at small scales to remove the noise. This solution procedure is a parabolic
equation with time as an evolution parameter and resembles the gradient-projection
method of Rosen (see [17]). In this formulation we assume an a priori estimate of the
Lagrange multiplier, in contrast with the dynamic change of - supposed in the Rosen
method, (see section 6 for details). The equation (3.1) moves each level curve of u
normal to itself with normal velocity equal to the curvature of the level surface divided
by the magnitude of the gradient of u, (see ([23]), ([15]) and ([20])). The constraints
are included in the -term and they are needed to prevent distortion and to obtain a
nontrivial steady state.
However, this evolution procedure is slow to reach steady state and is also stiff
since the parabolic term is quite singular for small gradients. In fact, an ad hoc rule
of thumb would indicate that the timestep \Deltat and the space stepsize \Deltax need to be
related by
\Deltat
for fixed c ? 0, for stability. This CFL restriction is what we shall relax. These issues
were seen in numerous experiments. In order to avoid these difficulties, we propose
a new time dependent model that accelerates the movement of level curves of u and
regularizes the parabolic term in a nonlinear way. In order to regularize the parabolic
term we multiply the whole Euler-Lagrange equation (2.14) by the magnitude of the
gradient and our time evolution model reads as follows:
We use as initial guess the original blurry and noisy image u 0 and homogeneous
Neumann boundary conditions as above, with an a priori estimate of the Lagrange
multiplier. From the analytical point of view this solution procedure approaches the
same steady state as the solution of whenever u has nonzero gradient. The effect
of this reformulation, (i.e. preconditioning) is positive in various aspects:
1. The effect of the regularizing term means that the movement of level curves
of u is pure mean curvature motion, (see [15]).
2. The total movement of level curves goes in the direction of the zeros of j   u\Gammau 0
regularized by the anisotropic diffusion introduced by the curvature term.
3. The problem for the denoising case is well-posed in the sense that there exists
a maximum principle that determines the solution, (see ([15])).
4. There are simple explicit schemes, such as Roe's scheme, that behave stably
with a reasonable CFL restriction for this evolution equation. Let us remark
that explicit schemes could also be applied for the 'anisotropic blurring' case.
5. This procedure is more morphological, (see [1]), in the pure denoising case,
i.e., it operates mainly on the level sets of u and u 0 . This is easily seen if
we replace u by h(u) and u 0 by h(u 0 ) with equation (3.3) is
invariant, except that replaced by (h(u) \Gamma h(u 0 ))=h 0 (u).
The anisotropic diffusion introduced in this model is a nonlinear way to discriminate
scales of computation. This never occurs with a linear model, (e.g. the linear
deconvolution model), because in this case we would have the linear heat equation
with constant diffusion. Thus, our model (3.3) can be seen as a convection-diffusion
equation with morphological convection and anisotropic diffusion.
4. Explicit numerical schemes for the 1D model. The 2D model described
before is more regular than the corresponding 1D model, because the 1D original
optimization problem is barely convex. For the sake of understanding the numerical
behavior of our schemes, we also discuss the 1D model. The Euler-Lagrange equation
in the 1D case reads as follows:
x
This equation can be written either as
x
using the small regularizing parameter fi ? 0 introduced at the end of the previous
section or
using the ffi -function.
The Rudin-Osher-Fatemi model, (ROF model), in terms of the ffi -function will read
as follows
Our model in 1D will be
x
regularizing parameter. The parameter fi ? 0 plays a more
relevant role in this case than in the 2D model. We can also state our model in terms
of the ffi function as
where a convolution of the ffi function must be used in practice. The intensity of this
kind of convolution decides which scale acts on the diffusion term. In this paper, we
always approximate ffi by
A radical way to make the coefficient of u xx nonsingular is to solve the evolution
model:
This model works in such a manner that away from extrema we have a large multiplier
of \Gammaj   (j   and at extrema it is just the heat equation.
These evolution models are initialized with the blurry and noisy signal u 0 and
homogeneous Neumann boundary conditions, and with a prescribed Lagrange multi-
plier. We estimated - ? 0 near the maximum value such that the explicit scheme is
stable under appropriate CFL restrictions, (see below).
In order to convince the reader about the speed and programming simplicity of
our model, we shall give the details of the first order scheme for the 1D pure denoising
model, i.e.,
x
Let u n
j be the approximation to the value u(x
Then, the scheme for the problem (4.9) will be
where
and ug j is the upwind gradient, i.e.,
\Deltax
if
\Deltax
if
Our general explicit scheme has the following features:
1. We use central differencing for u xx ,
2. The convolution operator j is computed by evolving the heat equation u
with the explicit Euler method in time and central differencing in space with
corresponding to a oe of the 1D heat kernel:
-oe
3. We use upwind Roe differencing, (see [16], [10]), checking the direction of
propagation by computing the sign of the derivative of the coefficient of j
(j   respect to u x times the sign of this term. Indeed, for our
evolution model (4.5) it is enough to check the sign of u x
For the model (4.8) we get the same direction of propagation as before. We
note that there is no notion of "entropy condition satisfying" discontinuities
in image processing; thus we omit the usual "entropy-fix" applied to the Roe
solver in this work.
4. The CFL condition depends on - and fi.
Indeed, the parabolic term in our model (4.5) gives a CFL restriction
\Deltat
x
and the convection term gives
\Deltat
s
x
for fixed c. These restrictions are reasonable at local extrema and near edges, compared
with the parabolic CFL restriction that corresponds to the reaction-diffusion
ROF model, (4.4):
\Deltat
which is too stiff along flat regions or at local extrema. The CFL restriction coming
from the convection term in the radical model (4.8) is better but also unfortunate
\Deltat
Thus, our model is more convenient from this point of view.
5. Explicit numerical schemes for the 2D model. We can express our 2D
model in terms of explicit partial derivatives as:
x
y
using u 0 as initial guess and homogeneous Neumann boundary conditions, (i.e., absorbing
boundary).
The denominator, u 2
y , appearing in the diffusion term may vanish or be small
along flat regions or at local extrema, when it is computed. Then, we can use either
the regularizing parameter fi ? 0, (small enough to perform floating point division),
or make the diffusion term equal to zero when gradient is smaller than a tolerance,
(we can also use parameter fi small as tolerance cut-off). Our choice in this paper
was the cut-off option, following a suggestion by Barry Merriman. Thus, concerning
stability and resolution the role of parameter fi is almost irrelevant in 2D calculations.
Let u n
ik be the approximation to the value u(x
and \Deltay and \Deltat are the spatial stepsizes and the time stepsize,
respectively. We denote by v
ik ). We point out that
we used for j, the convolution with the 2D heat kernel, (2.4), in our experiments,
aproximated by evolving the 2D heat equation u by means of the explicit
Euler method in time and central differencing in space. Then our first order scheme
reads as follows:
ik
ug x
ik
ik- (w n
where the second order term is defined by
if g x
ik
ik! fi and
ik
ik
otherwise, where
x
y
ik
2\Deltay
xx
ik
ik
yy
\Deltay 2
xy
ik
2\Deltax\Deltay
ug x
ik
is the upwind gradient in the x-direction, i.e.,
ug x
\Deltax (5.10)
if g x
ik
(w n
ug x
ik
\Deltax (5.11)
if g x
ik
(w n
ik
is the upwind gradient in the y-direction, i.e.,
ug y
\Deltay (5.12)
if g y
ik (w n
ug y
ik
\Deltay
if g y
ik (w n
A very simple way to extend this scheme to get high order accuracy is to follow
Shu-Osher prescription, (see [21]). Thus, we consider a method of lines, using an
explicit high order Runge-Kutta method in time and using a method of spatial ENO
reconstruction, (see [24], [9], [21] and [12]), of the same order, for the convection term,
applied on every time substep.
We have tested the Van Leer second order MUSCL spatial reconstruction using the
minmod function as slope-limiter together with classical second order Runge-Kutta
method and the third order PHM spatial reconstruction as in [12], using as slope-
limiter the harmod function, consisting of the harmonic mean of the lateral slopes
when they have the same sign and zero when they have different sign, together with the
third order Shu-Osher Runge-Kutta method of [21]. We have found that these explicit
methods are stable and give high accuracy under the same CFL restrictions as the first
order scheme.
As a sample we shall describe the second order MUSCL method. Since the Runge-Kutta
methods used here are linear combination of first order explicit Euler timesteps,
it is enough to formulate one Euler step, (in fact, in this case it is Heun's method
which is the arithmetic mean of two Euler timesteps). Following the notation used
above we have:
ik
ik
ik- (w n
where the reconstructed upwind gradients rug x
ik
and rug y
ik
are computed in the following
way. We reconstruct the left x-gradient in from the linear function:
\Deltax (5.15)
where
computed in x i , i.e.
gl x
where the minmod function is defined as
being sgn the sign function. Analogously, we have the reconstructed right x-gradient,
gr x
gr x
where
\Deltax (5.20)
where
Then the reconstructed upwind gradient in the x-direction is defined from the mean
value
as
if gm x
ik
if gm x
The procedure in the y-direction is similar.
-50501500 50 100 150 200 250 300
-5050150Fig. 6.1. Left, original vs. noisy 1D image; right original vs. recovered 1D image
6. Numerical Experiments. In this section, we perform some numerical experiments
in 1D and 2D.
We have used 1D signals with values in the range [0; 255]. The signal of (6.1, left)
represents the original signal versus the noisy signal with SNR - 5. The signal of
(6.1, right) represents the original signal versus the recovered signal after 80 iterations
with first order scheme with CFL 0:25. The estimated computed as
the maximum value allowed for stability, using the explicit Euler method in time. We
have used this experiment in order to achieve the appropiate amount of
difusion at small scales. In pure denoising 1D problems the choice of the value of fi
in our model depends on the SNR. Let us observe the very reduced staircase effect,
compared with the usual one obtained with either fixed-point iterative methods or
nonlinear primal-dual methods, (see [4]).
Now, we present a pure deblurring problem in 1D. The signal of (6.2, left) represents
the original signal versus the blurred signal with (as in 4.11. The
signal of (6.2, right) represents the original signal versus the recovered signal after 40
iterations with first order scheme with CFL 0:1. The estimated computed
as the maximum value allowed for stability, using the explicit Euler method in time.
We use 0:01 in this experiment.
The signal of (6.3, left) represents the original signal versus the blurred and noisy
signal with (as in 4.11), and SNR - 5. The signal of (6.2, right) represents the
original signal versus the recovered signal after 80 iterations with first order scheme
Fig. 6.2. Left, original vs. blur 1D image; right original vs. recovered 1D image
-50501500 50 100 150 200 250 300
-5050150Fig. 6.3. Left,original vs. noisy and blurred 1D signal ; right, original vs. recovered 1D signal
with CFL 0:25. The estimated computed as the maximum value allowed
for stability, using explicit Euler method in time. The - used for the current denoising
and deblurring problem is smaller than the one used in the above pure deblurring
problem, as we expected. We use this experiment to get the correct degree
of difusion at small scales. This shows that the 1D problem is quite sensitive to the
choice of fi, in contrast with the 2D case where the size of this parameter becomes
irrelevant. Let us also observe a very reduced staircase effect. We performed many
other experiments with 1D signals, obtaining similar results.
All our 2D numerical experiments were performed on the original image (Fig 6.4,
left) with 256 \Theta 256 pixels and dynamic range in [0; 255].
The third order scheme we used in our 2D experiments was based on the third
order Runge-Kutta introduced by Shu and Osher, (see [21]), to evolve in time with
a third order spatial approximation based on the PHM reconstruction introduced in
([12]).
Our first 2D experiment was made on the noisy image, (6.4, right), with a SNR
which is approximately 3. Details of the approximate solutions using the Chan-Golub-
Mulet primal-dual method and our time dependent model using the third order Roe's
scheme, (described above), are shown in Fig. 6.5. We used - 0:0713 and we perform
50 iterations with CFL number 0:1. We used the same estimated - as the one used for
the primal-dual method, and we observed that this value correponds to the largest we
50 100 150 200 250100200Fig. 6.4. Left: original image, right: noisy image, SNR- 3.
Chan-Golub-Mulet Primal-Dual
Resolution 256x256, SNR \approx 3, Estimated  =0.0713
50 100 150 200 250100200ROE-ORDER3-RK3
50 100 150 200 250100200Fig. 6.5. Left: image obtained by the Chan-Golub-Mulet primal-dual method, right: image
obtained by our time evolution model,with 50 timesteps and CFL-0.1
allowed for stability with this CFL restriction. We also remark that the third order
Runge-Kutta method used enhances the diffusion at small scales. The contour plots
are shown in Fig 6.6. We can infer from these contours that the edges obtained by
the new model are sharper than the ones obtained by the primal-dual method. This
might seem surprising, since the steady state satisfies the same equation (2.14) on the
analytic level. Numerically they are quite different because the approximation of the
convection term involves hyperbolic upwind ideas.
Our second 2D experiment is a pure deblurring problem. Fig (6.7, left), corresponds
to the original image blurred with Gaussian blur where as in (2.4). We
remark that we computed the convolution operator j by evolving the 2D heat equation
with explicit Euler method in time and central differencing in space with a CFL
number of 0.125, in order to test our model in practical conditions. In Fig (6.7, right),
we represent the approximation using our third order Roe's scheme where we perform
50 iterations with CFL number 0:1. We have used (the maximum value that
allows stability for the above CFL restriction), and We observe that the
scheme is not sensitive to the choice of fi provided the value be small enough, (smaller
than 0:1). This behavior is justified from the fact that the 2D problem is more regular.
Fig. 6.6. Left: isointensity contours of part of the image obtained by the primal-dual method,
right: isointensity contours of part of the image obtained by our time evolution model.
50 100 150 200 25010020050 100 150 200 250100200Fig. 6.7. Left: image blurred with Gaussian blur with image restored with our
model, using third order Roe's scheme with 50 timesteps and CFL-0.1.
The isointensity contours showed in (6.8) make clear the edge enhancement obtained
through our algorithm.
Our 2D critical experiment was performed on the blurry and noisy image represented
in Fig (6.9, left), with Gaussian blur where as in (2.4) and SNR - 5.
We have used the 0:01. We performed 50 iterations with a CFL
number of 0:1, using our third order Roe's scheme, obtaining the approximation represented
in figure (6.9, right). Let us observe the denoising and deblurring effect in
the isointensity contours picture represented in figure (6.10).
Finally, we shall include the convergence history of the two 1D experiments corresponding
to the pure denoising problem and a denoising and deblurring problem
presented above. In Figs 6.11 and 6.12 we represent the semilog plot of the L 2 -norm
of the differences between consecutive iterates versus the number of iterations and the
plot of the evolution of the total variation of the solution, respectively. We observe
'superlinear' convergence along the first third part of the evolution and linear convergence
along the remainder. We pointed out that all our experiments were performed
with a constant timestep and thus, the computational cost is very low compared with
the semi-implicit methods. These usually require one third of the number of iterations
we needed, but every step of the semi-implicit method requires about five iterations
Fig. 6.8. Left: isointensity contours of part of the blurred image, right: isointensity contours
of part of the image restored by using our time evolution model.
50 100 150 200 25010020050 100 150 200 250100200Fig. 6.9. Left: image blurred with Gaussian blur with noisy with SNR - 10, right:
image restored with our model, using third order Roe's scheme with 50 timesteps and CFL-0.1.
of the preconditioned conjugate gradient method to invert.
7. Concluding remarks. We have presented a new time dependent model to
solve the nonlinear TV model for noise removal and deblurring together with a very
simple explicit algorithm based on Roe's scheme of fluid dymamics. The numerical
algorithm is stable with a reasonable CFL restriction, it is easy to program and it
converges quickly to the steady state solution, even for deblurring and denoising prob-
lems. The algorithm is fast and efficient since no inversions are needed for deblurring
problems with noise. Our time dependent model is based on level set motion that
makes the procedure morphological and appears to satisfy a maximum principle in the
pure denoising case, using as initial guess the noisy image. We also have numerical
evidence, (through our numerical tests), of this stability in the deblurring case, using
the noisy and blurred image as initial guess.



--R


A variational method in image recovery
Modular solvers for constrained image restoration problems
Extensions to total variation denoising
Image recovery via total variation minimization and related problems
A nonlinear primal-dual method for total variation-based image restoration
Constrained restoration and the recovery of discontinuities
The theory of Tikhonov regularization for Fredholm integral equations of the first kind
Uniformly high order accurate essentially non-oscillatory schemes III
Numerical methods for conservation laws
The Fourier method for nonsmooth data
reconstructions for nonlinear scalar conservation laws
Restoring images degraded by spatially variant blur

Fronts propagating with curvature dependent speed: algorithms based on a Hamilton-Jacobi formulation

The gradient-projection method for nonlinear programming: Part II
Total variation based image restoration with free local constraints
Nonlinear total variation based noise removal algorithms
Cambridge University Press
Efficient implementation of essentially non-oscillatory shock capturing schemes II
Solutions of ill-posed problems
On the numerical solution of Fredholm integral equations of the first kind by the inversion of the linear system produced by quadrature
Towards the ultimate conservative difference scheme V.
Variational problems and PDE's for image analysis and curve evolution
Iterative methods for total variation denoising
--TR

--CTR
John Steinhoff , Meng Fan , Lesong Wang , William Dietz, Convection of Concentrated Vortices and Passive Scalars as Solitary Waves, Journal of Scientific Computing, v.19 n.1-3, p.457-478, December
Youngjoon Cha , Seongjai Kim, Edge-Forming Methods for Image Zooming, Journal of Mathematical Imaging and Vision, v.25 n.3, p.353-364, October   2006
Ronald P. Fedkiw , Guillermo Sapiro , Chi-Wang Shu, Shock capturing, level sets, and PDE based methods in computer vision and image processing: a review of Osher's contributions, Journal of Computational Physics, v.185 n.2, p.309-341, March
upwind schemes;total variation norm;image restoration;nonlinear diffusion;level set motion