EXPLICIT ALGORITHMS FOR A NEW TIME DEPENDENT MODEL BASED ON LEVEL SET MOTION FOR NONLINEAR DEBLURRING AND NOISE REMOVAL
ANTONIO MARQUINAy AND STANLEY OSHER
Dedicated to the memory of Emad Fatemi
Abstract. 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, ( 18]), and Rudin, Osher and Fatemi, ( 19]), 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, ( 16]), used in uid 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.
1. Introduction. The classical algorithms for image deblurring and/or denoising have been mainly based on least squares, Fourier series and other L2 -norm approximations, 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. However, the e ect 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 nondi erentiable functional. The formulation of these models was rst 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 di culties. Indeed, in spite of the fact that the variational problem is convex, the Euler-Lagrange equations are nonlinear and ill-conditioned. Linear semi-implicit xed-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 su er from the undesirable staircase e ect, 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 Avenue, Los Angeles, CA 90095-1555 and Departament de Matematica Aplicada, Universitat de Valencia, 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 PB971402. 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 9706827. 1
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 di usion acting on small scales. We use the entropy-violating Roe scheme, ( 16]) for the convective term and central di erencing for the regularized mean curvature di usion term. This makes a very simple, stable, explicit procedure, computationally competitive compared with other semi-implicit or implicitprocedures. 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 e ect. 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 x our ideas we restrict the discussion to R2. An image can be interpreted as either a real function de ned on , a bounded and open domain of R2, (for simplicity we will assume 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 u0 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; y; r) by means of v0 (x; y) =
Z
u(s; r) k(x; s; y; r) ds dr
(2.1)
and, we denote this operation by v0 = k u. The model of degradation we assume is k u + n = u0 ; (2.2) where n is Gaussian white noise, i.e., the values ni of n at the pixels i are independent random variables, each with a Gaussian distribution of zero mean and variance 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; y; r) = j(x ? s; y ? r) and the blurring is de ned as a 'superposition' of j 0 s: v0 (x; y) = (j u)(x; y) =
Z
u(s; r) j(x ? s; y ? 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 R properties, j(x; y) 0, j(x; y) ! 0 as (x2 +y2 )1=2 goes to 1 and 2 j(x; y) dx dy = 1. For any > 0 the so-called heat kernel, de ned as j(x; y) = 4 1 e?(x2 +y2 )=4 (2.4) is an important example that we will use in our numerical experiments.
R
2
The main advantage of the convolution is that if we take the Fourier transform of (2.3) we get v0 (k; l) = ^(k; l)^(k; l) ^ j u (2.5) then, to solve the model (2.2) with k = j we take Fourier transform and we arrive at ^(k; l)^(k; l) + n(k; l) = u0(k; l) j u ^ ^ (2.6) To recover u(x; y), we need to deconvolve, i.e., this means that we have to divide the equation (2.6) by ^(k; l) and to apply the inverse Fourier transform. This procedure j is generally very ill-posed. Indeed, j is usually smooth and j(x; y) ! 0 rapidly as (x2 + y2 )1=2 goes to 1, thus large frequencies in u0 get ampli ed considerably. The function u0 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 su ers from Gibbs' phenomenon. Discrete direct methods dealing with the linear integral equation (2.6) have been designed by di erent 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 R(u) u (2.7) subject to jjj u ? u0 jj2 2 = j j 2; L R since n = u0 ? j u and E( n2 dx) = j j 2 (E(X) denotes the expectation of the R random variable X) imply that jjj u ? u0 jj2 2 = (j u ? u0 )2 dx j j 2. L Examples of regularization functionals that can be found in the literature are, R(u) = k ukL2 ; krukL2 , where r is the gradient and 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: TV (u) =
Z
jruj dx =
Z q
u2 + u2 dx: x y
(2.8)
The TV 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 di erent purposes, (see for instance, 7, 5, 25, 2]). The restoration problem can be thus written as subject to 1 2 min u
Z
jruj dx;
Z
(j u ? u0)2 dx ? j j 2 = 0:
3
(2.9)
Its Lagrangian is
Z
jruj dx + 2
Z
(j u ? u0 )2 dx ? j j 2
(2.10)
and its Euler-Lagrange equations, with homogeneous Neumann boundary conditions for u, are:
ru 0 = ?r jruj + (j (j u ? u0 )) 1?Z 0 = 2 (j u ? u0)2 dx ? j j 2 :
(2.11) (2.12)
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 xed. Therefore, for the sake of clarity, we will assume the Lagrange multiplier to be known throughout the exposition. For = 1 , we can then write the equivalent unconstrained problem as
Z 1 min ( jruj + 2 (j u ? u0)2 ) dx u
(2.13)
and its Euler-Lagrange equation in the more usual form: 0 = ?r
ru jruj + j (j u ? u0 ):
(2.14)
We call (2.14) the nonlinear deconvolution model. The linear deconvolution model would be 0 = ? u + j (j u ? u0): (2.15) that comes from the Euler-Lagrange equation of the corresponding unconstrained problem with the norm R(u) = krukL2 , Since the equation (2.14) is not well de ned at points where ru = 0, due to the presence of the term 1=jruj, it is common to slightly perturb the Total Variation functional to become:
Z p
jruj2 + dx;
(2.16)
where is a small positive parameter, or,
Z
jruj dx;
p
(2.17)
with the notation (x 2 R, v 2 R2)
jxj = x2 + ; jvj = jvj2 + :
4
p
(2.18)
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 (2.14) 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 parameters. 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 (2.14) are ine cient 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]): r ut = ? j (j u ? u0 ) + r jruj : (3.1) u with u(x; y; 0) given as initial data, (we have used as initial guess the original blurry @u and noisy image u0) and homogeneous Neumann boundary conditions, i.e., @n = 0 on the boundary of the domain. As t increases, we approach to a restored version of our image, and the e ect 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 sti since the parabolic term is quite singular for small gradients. In fact, an ad hoc rule of thumb would indicate that the timestep t and the space stepsize x need to be related by t c j uj; (3.2) x2 for xed c > 0, for stability. This CFL restriction is what we shall relax. These issues were seen in numerous experiments. In order to avoid these di culties, 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: ru (3.3) ut = ?jruj j (j u ? u0 ) + jruj r jruj : We use as initial guess the original blurry and noisy image u0 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 (2.14) whenever u has nonzero gradient. The e ect of this reformulation, (i.e. preconditioning) is positive in various aspects:
5
1. The e ect 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?u0 regularized by the anisotropic di usion 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 u0 . This is easily seen if we replace u by h(u) and u0 by h(u0 ) with h0 > 0. Then, equation (3.3) is invariant, except that u ? u0 gets replaced by (h(u) ? h(u0 ))=h0 (u). The anisotropic di usion 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 di usion. Thus, our model (3.3) can be seen as a convection-di usion equation with morphological convection and anisotropic di usion. 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: u 0 = ? juxj + j (j u ? u0 ): (4.1) x x This equation can be written either as u (4.2) 0 = ? ju xj + j (j u ? u0 ): x x using the small regularizing parameter > 0 introduced at the end of the previous section or 0 = ? (ux )uxx + j (j u ? u0): (4.3) using the -function. The Rudin-Osher-Fatemi model, (ROF model), in terms of the -function will read as follows ut = ? j (j u ? u0) + (ux ) uxx (4.4) Our model in 1D will be ut = ?jux j j (j u ? u0 ) + + u2 uxx (4.5) x where > 0 is the small regularizing parameter. The parameter > 0 plays a more relevant role in this case than in the 2D model. We can also state our model in terms of the function as ut = ?juxj j (j u ? u0) + juxj (ux ) uxx (4.6)
6
where a convolution of the function must be used in practice. The intensity of this kind of convolution decides which scale acts on the di usion term. In this paper, we always approximate by (z) (z 2 + )?3=2 (4.7) A radical way to make the coe cient of uxx nonsingular is to solve the evolution model: 1 (4.8) ut = ? (u ) j (j u ? u0) + uxx: x This model works in such a manner that away from extrema we have a large multiplier of ?j (j u ? u0) and at extrema it is just the heat equation. These evolution models are initialized with the blurry and noisy signal u0 and homogeneous Neumann boundary conditions, and with a prescribed Lagrange multiplier. 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 rst order scheme for the 1D pure denoising model, i.e., ut = ?jux j (u ? u0) + + u2 uxx (4.9) x Let un be the approximation to the value u(xj ; tn), where xj = j x and tn = n t. j Then, the scheme for the problem (4.9) will be un+1 ? un un ? 2un + un j j = ?jugj j (un ? u0(xj )) + + g2 j +1 xj2 j ?1 j t j (4.10) un ? un gj = j +1 x j ?1 2 and ugj is the upwind gradient, i.e., un ? un ugj = j xj ?1 if gj (un ? u0(xj )) > 0 and j un ? un ugj = j +1 x j if gj (un ? u0(xj )) < 0 j Our general explicit scheme has the following features: 1. We use central di erencing for uxx , 2. The convolution operator j is computed by evolving the heat equation ut = uxx with the explicit Euler method in time and central di erencing in space with CFL = 0:25 corresponding to a of the 1D heat kernel: 2 1 (4.11) j(x) = 2p e?x =4
7
where
3. We use upwind Roe di erencing, (see 16], 10]), checking the direction of propagation by computing the sign of the derivative of the coe cient of j (j u ? u0 ) with respect to ux times the sign of this term. Indeed, for our evolution model (4.5) it is enough to check the sign of ux j (j u ? u0 ). 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- x" applied to the Roe solver in this work. 4. The CFL condition depends on and . Indeed, the parabolic term in our model (4.5) gives a CFL restriction t + u2 ; x (4.12) 2 x2 and the convection term gives t x c
s
1 + u2 ;
x
(4.13)
for xed c. These restrictions are reasonable at local extrema and near edges, compared with the parabolic CFL restriction that corresponds to the reaction-di usion ROF model, (4.4): 1 ; t (4.14) x2 2 (ux) which is too sti along at regions or at local extrema. The CFL restriction coming from the convection term in the radical model (4.8) is better but also unfortunate t (4.15) x 3juxj( + u2 )1=2 ; x 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: q u u2 ? 2u u u + u u2 ut = ? u2 + u2 j (j u ? u0) + xx y u2xy+ x 2y yy x : (5.1) x y x uy using u0 as initial guess and homogeneous Neumann boundary conditions, (i.e., absorbing boundary). The denominator, u2 +u2 , appearing in the di usion term may vanish or be small x y along at regions or at local extrema, when it is computed. Then, we can use either the regularizing parameter > 0, (small enough to perform oating point division), or make the di usion term equal to zero when gradient is smaller than a tolerance, (we can also use parameter small as tolerance cut-o ). Our choice in this paper was the cut-o option, following a suggestion by Barry Merriman. Thus, concerning stability and resolution the role of parameter is almost irrelevant in 2D calculations. Let un be the approximation to the value u(xi; yk ; tn), where xi = i x, yk = k y ik and tn = n t, where x, y and t are the spatial stepsizes and the time stepsize, n respectively. We denote by v0 = j u0 and wik = j j (un ). We point out that ik
8
we used for j, the convolution with the 2D heat kernel, (2.4), in our experiments, aproximated by evolving the 2D heat equation ut = uxx +uyy by means of the explicit Euler method in time and central di erencing in space. Then our rst order scheme reads as follows: un+1 ? un = ?qugx 2 + ugy 2 (wn ? v (x ; y )) + sn ik ik (5.2) 0 i k ik ik ik ik t where the second order term is de ned by sn := 0; (5.3) ik
y x if gik 2 + gik 2 < and xy x y yy x 2 xx y 2 sn := gik gik ? 2gik gik gik + gik gik ; ik y x gik 2 + gik 2
(5.4) (5.5) (5.6) (5.7) (5.8) (5.9)
otherwise, where
un ? un x gik = i+1;k2 x i?1;k ;
y gik =
un +1 ? un ?1 i;k i;k ; 2 y
n ? 2un + un ik xx u gik = i+1;k x2 i?1;k ; yy gik =
un +1 ? 2un + un ?1 i;k ik i;k ; y2
x ugik is the upwind gradient in the x-direction, i.e., un ? un x (5.10) ugik = ik xi?1;k x n if gik (wik ? v0 (xi ; yk )) > 0 and un ? un x (5.11) ugik = i+1;k x ik y x n if gik (wik ? v0 (xi ; yk )) < 0, and ugik is the upwind gradient in the y-direction, i.e., un ? un y ugik = ik yi;k?1 (5.12) y n if gik (wik ? v0 (xi ; yk )) > 0 and un ? un y (5.13) ugik = i;k+1 y ik
n ? un +1 ? un + un i xy u gik = i+1;k+1 i?1;k2 x y+1;k?1 i?1;k?1 ;
9
y n if gik (wik ? v0 (xi; yk )) < 0. 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 slopelimiter the harmod function, consisting of the harmonic mean of the lateral slopes when they have the same sign and zero when they have di erent 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 rst order scheme. As a sample we shall describe the second order MUSCL method. Since the RungeKutta methods used here are linear combination of rst 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: un+1 ? un = ?qrugx 2 + rugy 2 (wn ? v (x ; y )) + sn ik ik (5.14) 0 i k ik ik ik ik t where the reconstructed upwind gradients rugx and rugy are computed in the folik ik lowing way. We reconstruct the left x-gradient in (xi ; yk ) from the linear function: un ? u n (5.15) pl(x) := mi?1 (x ? xi?1=2) + ik xi?1;k where xx mi?1 = minmod(gixx1;k; gi;k) (5.16) ? computed in xi , i.e. glix := pl(xi); (5.17) where the minmod function is de ned as 1 minmod(r; s) := 2 min(jrj; jsj) (sgn(r) + sgn(s)); (5.18) being sgn the sign function. Analogously, we have the reconstructed right x-gradient, grix , as grix := pr(xi); (5.19) where un ? u n (5.20) pr(x) := mi (x ? xi+1=2) + i+1;k x i;k where xx mi = minmod(gi;k ; gixx ;k) (5.21) +1
10
Then the reconstructed upwind gradient in the x-direction is de ned from the mean value x x gmx := gli + gri (5.22) i 2 as
x rugik = glix ; n if gmx (wik ? v0 (xi; yk )) > 0 and i x rugik = grix ; n if gmx (wik ? v0 (xi; yk )) < 0. The procedure in the y-direction is similar. i
200 200
(5.23) (5.24)
150
150
100
100
50
50
0
0
?50
0
50
100
150
200
250
300
?50
0
50
100
150
200
250
300
Fig. 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 rst order scheme with CFL 0:25. The estimated = 0:05 was computed as the maximum value allowed for stability, using the explicit Euler method in time. We have used = 15 in 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 in our model depends on the SNR. Let us observe the very reduced staircase e ect, compared with the usual one obtained with either xed-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 = 10, (as in 4.11. The signal of (6.2, right) represents the original signal versus the recovered signal after 40 iterations with rst order scheme with CFL 0:1. The estimated = 1:5 was 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 = 5, (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 rst order scheme
11
200
200
180
180
160
160
140
140
120
120
100
100
80
80
60
60
40
40
20
20
0
0
50
100
150
200
250
300
0
0
50
100
150
200
250
300
Fig. 6.2. Left, original vs. blur 1D image; right original vs. recovered 1D image
200 200
150
150
100
100
50
50
0
0
?50
0
50
100
150
200
250
300
?50
0
50
100
150
200
250
300
Fig. 6.3. Left,original vs. noisy and blurred 1D signal ; right, original vs. recovered 1D signal
with CFL 0:25. The estimated = 0:25 was 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 = 10 in 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 , in contrast with the 2D case where the size of this parameter becomes irrelevant. Let us also observe a very reduced staircase e ect. 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 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 rst 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-GolubMulet 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
12
SNR=3
50
50
100
100
150
150
200
200
250 50 100 150 200 250
250 50 100 150 200 250
Fig. 6.4. Left: original image, right: noisy image, SNR 3.
Resolution 256x256, SNR \approx 3, Estimated \lambda=0.0713
ROE?ORDER3?RK3
50
50
100
100
150
150
200
200
250 50 100 150 Chan?Golub?Mulet Primal?Dual 200 250
250 50 100 150 200 250
Fig. 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 di usion 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 satis es the same equation (2.14) on the analytic level. Numerically they are quite di erent 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 = 5 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 di erencing 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 = 1:5, (the maximum value that allows stability for the above CFL restriction), and = 0:01. We observe that the scheme is not sensitive to the choice of provided the value be small enough, (smaller than 0:1). This behavior is justi ed from the fact that the 2D problem is more regular.
13
PRIMAL?DUAL
ROE?ORDER3?RK3
70
70
60
60
50
50
40
40
30
30
20
20
10
10
10
20
30
40
50
60
70
10
20
30
40
50
60
70
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
50
100
100
150
150
200
200
250 50 100 150 200 250
250 50 100 150 200 250
Fig. 6.7. Left: image blurred with Gaussian blur with = 5, right: 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 = 5 as in (2.4) and SNR 5. We have used the = 1:5 and = 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 gure (6.9, right). Let us observe the denoising and deblurring e ect in the isointensity contours picture represented in gure (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 L2-norm of the di erences 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 rst 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 ve iterations
14
70
70
60
60
50
50
40
40
30
30
20
20
10
10
10
20
30
40
50
60
70
10
20
30
40
50
60
70
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
50
100
100
150
150
200
200
250 50 100 150 200 250
250 50 100 150 200 250
Fig. 6.9. Left: image blurred with Gaussian blur with = 5 and 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 uid 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 problems. The algorithm is fast and e cient 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.
REFERENCES 1] L. Alvarez, F. Guichard, P.L. Lions and J. M. Morel Axioms and fundamental equations of image processing Arch. Rational Mechanics and Anal., v. 16, IX, (1993), pp. 199-257. 15
70
70
60
60
50
50
40
40
30
30
20
20
10
10
10
20
30
40
50
60
70
10
20
30
40
50
60
70
Fig. 6.10. Left: isointensity contours of part of the blurred and noisy image, right: isointensity contours of part of the image restored by using our time evolution model.
10
2
3500
3000
10
1
2500
2000
10
0
1500
1000
10
?1
0
10
20
30
40
50
60
70
80
500
0
10
20
30
40
50
60
70
80
Fig. 6.11. Left, semilog di erences of iterates vs. number of iterations for the pure denoising problem; right, total variation vs. number of iterations
2] G. Aubert and L. Vese, A variational method in image recovery, SIAM J. Numer. Anal., 34, (1997), pp 1948{1979. 3] P. Blomgren and T.F. Chan, Modular solvers for constrained image restoration problems, UCLA CAM report 97-52, (1997), 4] P. Blomgren and T. F. Chan and P. Mulet, Extensions to total variation denoising, Proc. SPIE 97, San Diego, (1997). 5] A. Chambolle and P.-L. Lions, Image recovery via total variation minimization and related problems, Numerische Mathematik, 76 (1997), pp. 167{188. 6] T. Chan, G. Golub, and P. Mulet, A nonlinear primal-dual method for total variation-based image restoration, SISC, (1998). To appear. 7] D. Geman and G. Reynolds, Constrained restoration and the recovery of discontinuities, IEEE Trans. on Pat. An. and Mach. Intel., 14 (1992), pp. 367{383. 8] C. W. Groetsch, The theory of Tikhonov regularization for Fredholm integral equations of the rst kind, Pitman, Boston, 1984. 9] A. Harten, B. Engquist, S. Osher and S. Chakravarthy, Uniformly high order accurate essentially non-oscillatory schemes III, J. Comput. Phys., v. 71 No. 2, (1987), pp. 231-303. 10] R.J. Leveque, Numerical methods for conservation laws, Birkhauser Verlag, Zuerich, (1990). 11] A. Majda, J. McDonough and S. Osher, The Fourier method for nonsmooth data, Math. Comp., 22, (1978), pp 1041-1081. 12] A. Marquina, Local piecewise hyperbolic reconstructions for nonlinear scalar conservation laws,, SIAM J. Scienti c Comp., v. 15, (1994) pp. 892-915. 13] James G. Nagy and Dianne P. O'leary, Restoring images degraded by spatially variant blur, SIAM J. Sci. Comput., 19, (1998), pp 1063-1082. 16
10
2
2600
2400
2200
10
1
2000
1800
1600
10
0
1400
1200
1000
10
?1
0
10
20
30
40
50
60
70
80
800
0
10
20
30
40
50
60
70
80
Fig. 6.12. Left, semilog di erences of iterates vs. number of iterationsfor the denoising and deblurring problem; right, total variation vs. number of iterations
14] S. J. Osher and L. I. Rudin, Feature-oriented image enhancement using shock lters, SIAM J. Numer. Anal., 27, (1990), pp 919{940. 15] S. J. Osher and J. A. Sethian, Fronts propagating with curvature dependent speed: algorithms based on a Hamilton-Jacobi formulation, J. Comput. Phys., 79, (1988), pp 12{49. 16] P.L. Roe, Approximate Riemann solvers, parameter vectors, and difference schemes , J. Comput. Phys., v. 43, (1981) pp. 357-372. 17] J.B. Rosen, The gradient-projection method for nonlinear programming: Part II. Nonlinear constraints, J. Soc. Indust. Appl. Math., v. 9, (1961) pp. 514-532. 18] L. Rudin and S. Osher, Total variation based image restoration with free local constraints, Proc. IEEE Internat. Conf. Imag. Proc., (1994), pp. 31{35. 19] L. Rudin, S. Osher, and E. Fatemi, Nonlinear total variation based noise removal algorithms, Physica D, 60 (1992), pp. 259{268. 20] J. A. Sethian, Level set methods, Cambridge University Press, (1996). 21] C. W. Shu and S. J. Osher, Ef cient implementation of essentially non-oscillatory shock capturing schemes II, J. Comput. Phys., v. 83, (1989) pp. 32-78. 22] A. N. Tikhonov and V. Y. Arsenin, Solutions of ill-posed problems, John Wiley, New York, 1977. 23] S. Twomey, On the numerical solution of Fredholm integral equations of the rst kind by the inversion of the linear system produced by quadrature, J. ACM, Vol. 10, (1963), pp. 97{107. 24] B. Van Leer, Towards the ultimate conservative di erence scheme V. A second order sequel to Godunov's method, J. Comput. Phys.,32,(1979), p.101-136. 25] L. Vese, Variational problems and PDE's for image analysis and curve evolution, PhD thesis, University of Nice, 1996. 26] C. R. Vogel and M. E. Oman, Iterative methods for total variation denoising, SIAM J. Sci. Statist. Comput., 17 (1996), pp. 227{238.
17