Scale-Space Methods and Regularization for Denoising and Inverse Problems

ADVANCES IN IMAGING AND ELECTRON PHYSICS, VOL. 128

Scale-Space Methods and Regularization for Denoising and Inverse Problems OTMAR SCHERZER Department of Computer Science, Universit6t Innsbruck, TechnikerstraJ3e 25, A-6020 Innsbruck, Austria

I. I n t r o d u c t i o n . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . II. I m a g e S m o o t h i n g a n d R e s t o r a t i o n via Diffusion Filtering . . . . . . . . . . . . A. Level Set M o d e l i n g . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . B. M o r p h o l o g i c a l Diffusion Filtering . . . . . . . . . . . . . . . . . . . . . . . . C. A p p l i c a t i o n s of Diffusion Filtering . . . . . . . . . . . . . . . . . . . . . . . D. Scale-Space T h e o r y . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . III. R e g u l a r i z a t i o n o f Inverse P r o b l e m s . . . . . . . . . . . . . . . . . . . . . . . . . A. T i k h o n o v - T y p e R e g u l a r i z a t i o n M e t h o d s . . . . . . . . . . . . . . . . . . . . B. R e g u l a r i z a t i o n M o d e l s for D e n o i s i n g . . . . . . . . . . . . . . . . . . . . . . C. R e l a t i o n s b e t w e e n R e g u l a r i z a t i o n a n d P e r o n a - M a l i k Diffusion F i l t e r i n g . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . D. N u m e r i c a l E x p e r i m e n t s . . . . . . . . . . . . . . . . . . . . . . . . . . . . . IV. M u m f o r d - S h a h Filtering . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . V. R e g u l a r i z a t i o n a n d Spline A p p r o x i m a t i o n . . . . . . . . . . . . . . . . . . . . . VI. Scale-Space M e t h o d s for Inverse P r o b l e m s . . . . . . . . . . . . . . . . . . . . . A. D e b l u r r i n g w i t h a Scale-Space M e t h o d . . . . . . . . . . . . . . . . . . . . . B. N u m e r i c a l S i m u l a t i o n s . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . VII. N o n c o n v e x R e g u l a r i z a t i o n M o d e l s . . . . . . . . . . . . . . . . . . . . . . . . . A. P e r o n a - M a l i k R e g u l a r i z a t i o n . . . . . . . . . . . . . . . . . . . . . . . . . . B. Relative E r r o r R e g u l a r i z a t i o n . . . . . . . . . . . . . . . . . . . . . . . . . . VIII. Discrete BV R e g u l a r i z a t i o n a n d T u b e M e t h o d s . . . . . . . . . . . . . . . . . . A. Discrete BV R e g u l a r i z a t i o n ( S a m p l i n g ) . . . . . . . . . . . . . . . . . . . . . B. F i n i t e V o l u m e BV R e g u l a r i z a t i o n . . . . . . . . . . . . . . . . . . . . . . . . C. T h e T a u t String A l g o r i t h m . . . . . . . . . . . . . . . . . . . . . . . . . . . D. M u l t i d i m e n s i o n a l Discrete BV R e g u l a r i z a t i o n . . . . . . . . . . . . . . . . . E. N u m e r i c a l Test E x a m p l e s . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1. O n e - D i m e n s i o n a l Test E x a m p l e . . . . . . . . . . . . . . . . . . . . . . . 2. T w o - D i m e n s i o n a l B e n c h - M a r k P r o b l e m . . . . . . . . . . . . . . . . . . IX. W a v e l e t S h r i n k a g e . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . A. D a u b e c h i e s ' W a v e l e t s . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . B. D e n o i s i n g b y W a v e l e t S h r i n k a g e . . . . . . . . . . . . . . . . . . . . . . . . 1. R e l a t i o n to Diffusion F i l t e r i n g . . . . . . . . . . . . . . . . . . . . . . . . X. R e g u l a r i z a t i o n a n d Statistics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . XI. C o n c l u s i o n s . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Acknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

446 447 453 455 458 458 460 464 465 466 468 472 474 478 481 484 493 493 494 500 502 504 505 506 508 509 510 510 511 514 514 517 522 523 523

445 Copyright 9 2003 Elsevier Inc. All rights reserved. 1076-5670/2003 $35.00

446

OTMAR SCHERZER I. INTRODUCTION

Inverse problems and imaging are two of the fastest growing areas in applied mathematics. Such problems appear in a variety of applications such as medical imaging and nondestructive evaluation. A typical example is computerized tomography (CT) where the density of a body is determined from X-ray measurements at the boundary. Inverse problems can be vaguely characterized as the problems of estimating the cause for an observed effect; in CT the cause is the density of the body and the observed effect is the X-ray data at the boundary of the object. With inverse problems one typically associates ill-posedness, that is, that there may not exist a solution, the solution is nonunique, or the solution does not depend continuously on the input data. In order to overcome these difficulties Tikhonov suggested approximating the ill-posed problem by a scale of well-posed variational problems. This initiated the work on regularization methods for the solution of ill-posed problems. Partial differential equations (PDEs) have proved to be efficient methods in image processing and computer vision. They are mainly used for smoothing and restoration, in particular noise removal. Their success is partly due to the fact that the approximation is independent of the underlying numerical method. The success of PDE methods in image processing has stimulated the development of new efficient numerical algorithms for the solution of inverse problems by constructing variational methods based on the energy formulations of PDEs. Nowadays the interaction between PDE models and variational formulations is subtle and has led to a fruitful interaction of inverse problems and image processing with splines, wavelets, morphology, and statistics. A goal of this survey is to review these interactions. The second goal of this survey is to compare various reconstruction algorithms. The outline of this work is as follows. In Section II we review image smoothing and restoration with PDEs. We compare several noise removal (denoising) techniques and show the effect of filtering as a prerequisite step of image analysis, such as segmentation. Moreover, we use the analogy of fluid flow to motivate PDEs for diffusion filtering. In Section III we review regularization methods for the solution of inverse problems and establish the connection between PDEs and variational methods for denoising. Section IV is devoted to the Mumford-Shah filtering method, which is a combined method for image smoothing and segmentation. In Section V we review the interaction between approximate spline filtering and variational methods. Section VI establishes a diffusion framework for the solution of inverse problems, which is linked to Section VII where nonconvex

DENOISING AND INVERSE PROBLEMS

447

variational problems are considered. In Section VIII we introduce a discrete framework for regularization and in Section IX we highlight the relation of variational methods and diffusion filtering with wavelets. In Section X we review the interactions of regularization and statistics.

II. IMAGE SMOOTHING AND RESTORATION VIA DIFFUSION FILTERING PDE-based models have proved to be efficient in a variety of image processing and computer vision areas such as restoration, denoising, segmentation, shape from shading, histogram modification, optical flow, and

stereo vision. To demonstrate the efficiency of diffusion filtering we recall a few models. To this end let u s be an image defined on the open domain : - (0, 1) • (0, 1). (1) The simplest and best investigated diffusion filtering technique for image smoothing is the linear heat equation 0u

--=Au, Ot

(1)

associated with homogeneous Neumann boundary data 0u

--=0;

Ov

here and in the following Ou/Ov denotes the derivative of u in normal direction v to the boundary Of2 of f2. As initial data we use the input image

u(O, x) - uS(x) for x ~ if2.

(2)

It is well known that the heat equation blurs the initial data and spurious noise is filtered. Figure 1 shows heat equation filtering of ultrasound data at specified times. (2) The heat equation is equally efficient in removing noise and destroying image details such as edges and corners. The total variation flow equation is able to preserve edges and denoise the image simultaneously. Here the differential equation

448

OTMAR SCHERZER

FIGURE 1. Solution of the heat equation at time t - O, O. 1, 1, 1O, 100, 1000.

together with homogeneous Neumann data is applied to the initial data u8. Here and in the following I. ] denotes the Euclidean norm. A detailed rigorous mathematical analysis of this partial differential equation (see [14-17,23]) has been given. The mathematical analysis impressively supports the remarkable properties of this filtering technique (cf. Figure 2). (3) The Binghamfluidflow equation is able to preserve flat regions and denoise the image simultaneously. This filtering technique requires one to solve the differential equation

D E N O I S I N G A N D INVERSE P R O B L E M S

449

FIGURE 2. Solution of the total variation flow equation at time t = 0, 0.01,0.05, 0.1,0.5, 1.

together with homogeneous Neumann conditions and initial data us. The parameters r and • are strictly positive. Bingham fluid flow is a widely investigated model in fluid mechanics (see, e.g., [80]), in which the parameters r and ~; have the physical meaning of yield stress and plastic viscosity. The particular properties of Bingham fluids make them extremely useful for image denoising (see [77]). Figure 3 shows the solution of Equation (4) at specified times. (4) Linear anisotropic diffusion filtering is based on matrix-valued diffusivity. Let u~ be a smooth approximation of us, then a linear anisotropic diffusion equation is Obl -

-

=

Ot

v . (D(Vu.)Vu),

(5)

450

O T M A R SCHERZER

FIGURE 3. Bingham filtering technique

with K = 0 . 0 5

and r = 0 . 0 5

at t = 0 , 0 . 2 ,

2, 10, 100, 1000.

where -k-/x2I , IVvl 2 + 2/x 2

-(Ov/Ox)

(6)

-~y'

with /x > 0. In image processing the differential Equation (5) is considered with homogenous Neumann boundary conditions (D(u.)Vu)

. v =

0

and initial data US. 1 Uo- is considered an approximation of the filtered data which can be obtained, for instance, by solving the heat 1The product of two vectors x

Zin--_l xiYi.

=

(x 1 .....

Xn) and y = (Yl . . . . . y.) is defined by x . y =

DENOISING AND INVERSE PROBLEMS

451

equation with initial data ua up to a certain time. In anisotropic diffusion models the matrix D(Vu~) is designed in such a way that the eigenvectors ~1 and ~2 a r e parallel, respectively orthogonal to Vuo. These methods prefer diffusion along edges to diffusion perpendicular to them. (5) Nonlinear anisotropic diffusion utilizes a matrix-valued diffusivity which itself depends on the solution. A typical example is 0u

at

(7)

= V . (D(Vu)Vu)

where D(. ) is as defined in (6). In Figure 4 we have evolved an image according to the nonlinear anisotropic diffusion equation. (6) The classical Perona-Malik filter [134,135] is the oldest nonlinear diffusion filter. It is based on the equation

ou

(

Ot =

1

)

1 + [Vul2/~. Vu

(8)

with a positive parameter X. In comparison with total variation flow the diffusivity is smaller near edges. So far no completely successful mathematical analysis for this model has been obtained. A few results concerning the existence of a solution have been given in [102]. (7) Mean curvature motion is a widely inspected model in applied mathematics describing phenomena such as crystal growth and polymer processing. The mean curvature equation

ou 0-7 = IV u l V .

Vu

(9)

is a paradigm of morphological differential equations. The image evolution according to the mean curvature motion is shown in Figure 5. It is common to distinguish between two classes of diffusion filtering techniques: 9 Giving tribute to Perona and Malik [134] for initiating the use of nonlinear diffusion filtering we call any of the differential equations Obl - - = V . (D(Vu)Vu) Ot

(10)

452

OTMAR SCHERZER

FIGURE 4. Solution of nonlinear anisotropic diffusion at time t = 0 , 10,500, 1000, 5000, 10,000; /x = 10-4.

Perona-Malik diffusion filtering. Prototypes are the heat equation, the total variation, the Bingham model, as well as anisotropic diffusion. Typically one differentiates also between isotropic Perona-Malik filtering models, where D is a one dimensional function and anisotropic Perona-Malik, where D is a matrix-valued function with nonzero entries in the off-diagonals. 9 Morphological partial differential equations are invariant under image transformation such as gray level modification and data set deformation. A paradigm of a morphological differential equation is the mean curvature flow equation. The use of Perona-Malik diffusion filtering can be motivated by showing its analogy to fluidflow.

DENOISING

453

AND INVERSE PROBLEMS

FIGURE 5. Solution o f the m e a n c u r v a t u r e e q u a t i o n at time t = 0, 0.1, 1, 10, 100, 1000.

A. Level Set Modeling We consider the movement of a fluid in ~2 = (0, 1) x (0, 1) over time. The conservation of mass principle (see, e.g., [48]) is expressed by the differential equation

Op

(t, x) + v . (prO(t, x) - o,

(11)

where p is the density of the fluid and ~ is the velocity field of the fluid. For a fixed time t we consider a level curve of the density. A level curve x(r) = (Xl (r), x2(r)), parameterized by r, in ~2 satisfies

p(t, x(r)) = constant for all r.

454

OTMARSCHERZER

Differentiating this equation with respect to r we find that

8p (t, x(O) Oxl Op Ox2 -g/- (0 + b-x2 (t, x(0)-g-r (0 - 0.

(12)

The direction of the tangent at the level curve at x(r) is

--, ( (oqX1/oq'g(T)) ) t "-

(ax2/Or) (r) "

If 7-r 0 and

Op Op ) Vp(t, x(r)) -- O-XXl'8-~2 (x(0) r 0, then from (12) it follows that the vectors Vp and 7T are orthogonal; here and in the following 9r denotes the transpose of a vector or matrix. Fick's law states that the velocity ~ is orthogonal to the level curves, i.e.,

9- - C V p with - C < O; the negative sign indicates that the direction of the flow is from regions of high density to regions of low density. Different models can be imagined by adequately choosing C: (1) If C = 1/p, then we get the diffusion filtering technique (1). This choice of C represents the fact that an object in a fluid of higher density moves slower than an object in a fluid of lower density. (2) If C = 1/(p[Vp[), then we get the diffusion filtering technique (3). The choice of C represents the fact that an object in a fluid moves faster at smooth portions of level curves. (3) If C - (1/p)D(Vp~) the flux is biased both in the tangential direction and in normal direction to the level curve. To derive the analogy of image diffusion filtering and fluid flow we identify the gray value image data with the density of a fluid. Since the fluid flow equations have been derived from the conservation of mass principle, we have mean gray value invariance for diffusion filtering in image processing, that is

fa p(t, x) dx - constant over time.

DENOISING AND INVERSE PROBLEMS

455

Without external forces a fluid flow is subsequently simplified: for instance, we might expect that the entropy of the density increases over time and that finally the density approaches a constant function. The analogy between fluid flow and image diffusion filtering suggests a subsequent simplification of the gray value data leading to a scale space of images. In image processing a quantization of these phenomena via Lyapunov functionals has been given by Weickert [172].

B. Morphological Diffusion Filtering

Morphological diffusion filtering techniques, such as (9), are closely related to shape and curve evolutions. To illustrate this connection we recall the definition of curvature. Let c:[0, 2rt) ~

[~2

Z-~__._(> xl (z')) x2(r) be a closed parameterized curve in ~2, then the standard definition (see, e.g., [33]) of curvature is

(OXl /O'C) (02X2/O'g 2) -- (02Xl /OV22)(OX2/O'C) ((OXl/O'g) 2 + (0X2/0"C)2) 3/2 where 0 - / 0 r denotes the derivative with respect to the curve parameter r. Let C : [0,cx~) x [0,2rt) ~ ~2 (xl(t, r ) ) (t, r) ~

x2(t, r)

be a temporally varying oriented closed curve. We consider the curvaturebased evolution process OC ~ ( t , r) =/3(/C)(t, r)v(t, r), Ot

(13)

where v denotes the normal vector to the curve C, and/~ is an appropriate scalar-valued function.

456

OTMAR SCHERZER

Let f 6 C2([0, OO) • ~ ) . We assume that the zero

level set

/~(t) "-- {X -- (Xl, X 2 ) ' f ( t , X) -- O}

can be parameterized by a curve C(t, 9) which evolves according to (13) and that f is locally invertible in a neighborhood N" of the zero-level set, i.e.,

Vf-

Oxi'0-~2

#OinN'.

Then f(t, C(t, r ) ) - 0 for all r e [0,2rt) and t 6 [0, oc). Consequently, by differentiation with respect to r and t we get

OC (t, r) +57 Of (t, C(t, r)) _ O, vf(t, c(t, r)).-~ OC Vf (t, C(t, r)) .-~r (t, r) - O,

(14)

for all r 6 [0,2r0, t e [0, oe). The latter equation shows that Vf(t, C(t, r)) and the tangential vector on the level curve (OC/Or)(t,r)- (OXl/Or, Ox2/Or)r(t,r) are orthogonal, which implies that Vf(t,C(t,r)) is proportional to the normal vector (Ox2/Or,-OXl/Or)r(t,r) on the level curve, that is

gOfl '

vf(t, c(t, ~)) -

Of )(t, C(t, r))

(ox2 OXl)

= gr(t, r) \ Or '

Or

Since we assumed that Ig/I ~ 0 we find that differentiation of (15) with respect to r we get

(15)

(t, r). 7t(t,r)#0.

Then

by

a~x2

~ ~f af ozf Of) O~ Of (t, C(t, r)) ax~ ax2 aXlaX2 ~ (t, c(t, ~))-~(t, ~)5-~x~

a2Xl

~ ~f af

~2(t, ~) ~ ~2 ( t, ~) ~

(t, r)( t, ~) -

ax~ ax

_

a2/ _

af (t, c(t, ~)) +

aXlaX~~

a~ (t, r) af (t, c(t, r)) ~ ~ "

DENOISING AND INVERSEPROBLEMS

457

Consequently, we have

/(7=

1 { ( ~ f ) 2~2f [V-fl3 ~x2 ~ -

02f

~fOf 2 O~ gx 2 0X l ~

~f(~f)2} {- -~x~ ~X l

-1 "

(16)

Let 7t < 0, then

--

V IK Or J

I,l x2 1

nt- ~, 072 ,1

OXl ( t, T)

-

-gg~

--

vf IVfl"

Note that iff(t,. ) is monotonically increasing into the interior of a domain (2(0 with boundary E(t), then - V f / l V f l points in outside direction of (2(0, which implies that 7t < 0. Using the abbreviations H i for the Hessian o f f it follows that

curv(f) "- V.

I-~

]Vfl2Af- gfr gf vf

IVfl 3 =/C.

(17)

This, together with (14), shows that the level set formulation of (13) is

~f

- - ( t , x) -- fi(curv(f)(t, x))lVf (t, x)]. Ot

(18)

Examples of curvature-based morphological processes are summarized in the following: 9 fl(t) 91 9 -1 9t 9 t 1/3

morphological process dilation erosion mean-curvature flow affine invariant mean-curvature flow.

Morphological diffusion filtering methods have been derived axiomatically in [2-4,37].

458

OTMAR SCHERZER

C. Applications of Diffusion Filtering The design and mathematical analysis (including existence, uniqueness, and stability results for solutions) of diffusion filtering techniques are active research areas. Appropriate filtering is important as a prerequisite step for image segmentation and edge detection, to mention but a few applications. In the following we apply Canny's edge detection algorithm from the software package MATLAB [112] to the filtered ultrasound examples in Figures 1, 2, 4, and 5, respectively. There are various parameters to be tuned in the implementation of Canny's edge detection algorithm, which might of course have considerable effect on the detected edges. For reasons of comparison we have used the standard MATLAB setting, which does not require input parameters. Canny's edge detector is a sophisticated algorithm to extract edges in image data. For more background on this method we refer the reader to [36,108]. As can be realized from Figures 6-10, appropriate filtering is an important prerequisite for edge detection.

D. Scale-Space Theory Images contain structures at a variety of scales. Any feature can optimally be recognized at a particular scale. This has already been observed in the edge detection example above. If the optimal scale is not available a priori, it is desirable to have an image representation at multiple scales. A scale space is an image representation at a continuum of scales, embedding the image u s into a family

{Tt(uS)" t >_ 0} of gradually simplified versions satisfying:

(1) Fiaetity: To(u s) -- u s. (2) Causality."

Tt+s(u s) - Tt(Ts(uS)) for all s, t _> 0. (3) Regularity." lim Tt(u ~) - u 8.

t---~O+

DENOISING AND INVERSE PROBLEMS

459

FIGURE 6. Canny's edge detector applied to the solution of the heat equation at time t - 0,0.1, 1, 10, 100, 1000. At t = 1000 the image is so blurred that no edges could be detected.

The differential equations introduced above satisfy these properties with Tt(u ~) -- u(t,.) .

In mathematics, a family of operators Tt satisfying fidelity, causality, and regularity is called s e m i - g r o u p . For more background on semi-group theory we refer the reader to Pazy [133] (in the linear case) and Br6zis [32] (in the nonlinear case).

460

OTMAR SCHERZER

FIGURE7. Canny's edge detector applied to the solution of the total variation flow equation at time t = 0, 0.01,.0.05,0.1,0.5, 1.

III. REGULARIZATION OF INVERSE PROBLEMS A vague characterization of inverse problems is that they are concerned with determining causes for a desired or an observed effect. Such problems appear in a variety of applications like

(1) Medical imaging such as CT (see, e.g., [25,94,123,170]). A mathematical framework for CT has been analyzed by R a d o n [140]. The theory has been applied in other areas including radioastronomy (e.g., [28]) and electron microscopy (e.g., [79]). (2) Signal and Image processing, such as the extrapolation of band-limited functions (see, e.g., [35]).

DENOISING AND INVERSE PROBLEMS

461

FIGURE 8. Canny's edge detector applied to the solution of the Bingham filtering technique at time t = 0, 0.2, 2, 10, 100, 1000.

It is well known that many inverse problems violate Hadamard's principle of well-posedness, that is, at least one of the following postulates is violated: (1) There exists a solution. (2) The solution is unique. (3) The solution depends continuously on the input data. If one of these properties is violated the problem is said to be ill-posed or improperly posed. Regularization methods are numerical algorithms for solving ill-posed problems in a stable way. In the linear setting Torre and Poggio [165] emphasize that differentiation is ill-posed, and that applying suitable regularization strategies approximates linear diffusion filtering or--equivalently--Gaussian convolution. Much of the linear-scale-space

462

OTMAR SCHERZER

FIGURE 9. Canny's edge detector applied to the solution of nonlinear anisotropic diffusion at time t = 0, 10, 500, 1000, 5000, 10,000; tz -- 10 -4.

literature is based on the regularization properties of convolutions with Gaussians. In particular, differential geometric image analysis is performed by replacing derivatives by Gaussian-smoothed derivatives; see, e.g., [76,106,126,156] and references therein. In order to present a general framework of regularization methods it is convenient to consider an inverse problem as the problem of solving an illposed operator equation f ( u ) = y0. Here 9c: D(9c) _c X ~ D(U) of a space X.

(19)

Y is an operator defined on an appropriate subset

DENOISING AND INVERSE PROBLEMS

463

FIGURE 10. C a n n y ' s edge detector applied to the solution of the mean curvature motion at time t = 0, 0.1, 1, 10, 100, 1000.

We use the following terminology: 9 Linear inverse problems: if 3c is a linear operator. (1) If r the identity operator, then the linear inverse problem is called denoising. (2) If ~ is a convolution, that is

fu(x) - f~

k(lx - yl2)u(y) @,

(20)

with k being a smooth function, and l" I being the Euclidean distance, then the problem is referred to as deblurring.

464

OTMAR SCHERZER

(3) If f ' i s the Radon transform (see, e.g., [123]), then the problem of solving (19) is the problem of computerized tomography. 9 Nonlinear inverse problems" if 9r is a nonlinear operator.

Regularization methods were first considered by Tikhonov in 1930. Since that time regularization theory has developed systematically. (1) During the 1980s there was success in a rigorous analysis of linear ill-posed problems. We mention the books of Louis [107], Groetsch [86], Tikhonov and Arsenin [164], Morozov [116], Nashed [122], Engl and Groetsch [66], Natterer [123,124], Bertero and Boccacci [25], Kirsch [103], and Colton and Kress [52,53,105]. See also Groetsch [83,84] for some elementary introduction in the topic of inverse problems. (2) Since 1989, starting with three fundamental papers of Seidman and Vogel [154] and Engl and co-workers [68,125], regularization theory for nonlinear inverse problems developed systematically. Some expository books on this topic are [21,67,99,116,117], to name but a few. (3) Acar and Vogel [1] and Geman and Yang [78] proposed a novel framework of nondifferentiable regularization of Tikhonov type. This work stimulated the development of regularization methods for efficiently recovering discontinuous solutions in inverse problems.

A. Tikhonov-Type Regularization Methods Tikhonov proposed approximating the solution of the operator Equation (19) by the minimizer of the functional (Tikhonov regularization)

f(u) "- II~'(u) - Yall~ + otllu - u, [IX" 2

(21)

Here, u, E X is some initial (a priori selected) guess on the desired solution and y8 is an approximation of the right-hand-side data in (19). The classical theory of regularization methods assumes a Hilbert space setting, that is X and Y are Hilbert spaces and that 9v" D(.F') c X ~ Y is

(1) continuous and (2) weakly (sequentially) closed, that is for any sequence {Un}nE~J~ D ( f ) , Xn~x x and .T(Xn)~ r Y imply x E D(9r) and ~ - ( x ) - y The Hilbert space X and Y are associated with inner products and norms, (',')x

and

(',')r,

]].]lx and

I].]]r,

DENOISING AND INVERSE PROBLEMS

465

respectively. In almost every application considered in practice Y--L2(~) with inner product

~C,g)_ j~ fg is used; typically for X Sobolev spaces of weakly differentiable functions are used. There exists a variety of results showing that Tikhonov's approach in fact yields a regularization method, that is (1) there exists minimizer of (21) and (2) for fixed a > 0 the minimizers are stable with respect to perturbations in ya. (3) Even more, it has been proved that for an appropriate choice of a the minimizer is an approximation of the solution of (19).

B. Regularization Models for Denoising For denoising we have ~ = I and ya = ua. Tikhonov-type regularization methods for denoising consist in minimizing a functional

f(u)'-- f (u-u~)2+ot f g([Vu[2).

(22)

(1) g(t)--t is refered to as HI-semi-norm regularization. (2) A popular specific energy functional arises from unconstrained total variation denoising [1,41,43,46]. Here g(t)= ,47. This method is called BV-semi-norm regularization. (3) The combination of H 1- and BV-semi-norm regularization gives g(t) = ~ct+ ~x/t. This method exhibits similar filtering properties as the Bingham fluid flow. (4) The regularization counterpart to linear anisotropic diffusion filtering consists in minimizing the functional

f(

u

--

US)2 -[- 0t f~ [A 1/2 VVu[2,

where D(Vu~) =

V TA V

(23)

OTMAR SCHERZER

466

is the singular value decomposition of D(Vuo) with

:xl/2 0 1 A1/2 _ [ "'I 0 11/2 ""2 with hi and ~2 the singular values of D(Vu~). (5) Minimizing the functional in (23) where A and V are dependent on u results in the regularization counterpart to nonlinear anisotropic

diffusion. C. Relations between Regularization and Perona-Malik Diffusion Filtering Let us assume that the functional (22) is defined on the Sobolev space HI(F2), that is the space of weakly differentiable functions. Moreover, we assume that there exists a minimizer, which is denoted by u~. Then for any h ~ Hl(f2) and any real number t the definition of u~ implies that

f (u~ + t h ) - f (u~) > O, which is equivalent to

fo {('-+ ,~-uS)2_(ua

_

4- ot ff2 {g([Vua[2 4- 2 t V u ~ .

uS)2} V h 4-

t2[Vh[2) -g([Vu~[2)} > 0.

If g is twice differentiable, then by making a Taylor series expansion, we find

g (IVUa]2 4- 2t(Vu~.

Vh) + t2 ]Vh]2)

= g(lVual 2) 4- (2tVu~. Vh + t2lVhl2)g'(lVu=l2) 4- O(t 2) = g(IVu~l 2) + 2t(Vu=. Vh)g'(IVu=l2) + O(t2). Therefore, for t > 0 we have

0 < if2(2h(uot-u')4-th2)4- 2.f.(Vu..Vh)g'(IVu.l 2)+ O(t). Taking the limit t ~ 04- shows

0 < f h(ur u~) 4-~ f (Vu~. Vh)g'(lVu~12). d ~2

d ~2

467

D E N O I S I N G AND INVERSE PROBLEMS

Simulating the above calculation with - t instead of t gives

O> ffa h(u~ - ua) + ~ fa

(Vu,~.

Vh)g'(lVua]2).

(Vu~.

Vh)g'(lVuo~12).

Thus, in total, we have

0 -- fa h(u,~- u~) +or fe

Then, by using Green's formula, we get

O- fe h(u,~-ua)-ot f v. (g'(lVuotl2)Vu~)h-+-otfoa g'(lVu~176

=

fo

hgt([Vuot]2)(gua.

v)h ~),

(24)

where v denotes the unit normal vector on O~. If g' > 0, then since (24) holds for all h 6 H 1(~), we find Uor - -

U~

= V. (g'([Vu~[2)Vu~) on OUot

0 ----

Ov =

(Vu~. v) on aS"2.

fa, (25)

In particular, setting a - A t , u(O)--u a, u(At)--u~ shows that Tikhonov regularization with small regularization parameter At provides an approximation of the solution of the diffusion filtering method oqb/

at 0U

Ov

-- V "

=0on

(gt([VU[2)VU)o n

~2,

0~,

u(0) - u~ on ~2, at time At. In other words, the regularization parameter and the diffusion time can be identified if regularization is regarded as time-discrete diffusion filtering with a single implicit time step [115,138,139,145,148,158]. Moreover,

468

OTMAR SCHERZER

iterated regularization with small regularization parameters, consisting in subsequently minimizing the functionals

f(k)(u) "-- ~ (U -- u(k-1))2+ ot s g(]VU]2),

k = 1,2...

(26)

and denoting the minimizers by u(k) approximates a diffusion process. The basic connection between regularization and diffusion filtering methods is the basis of both practical considerations and fundamental mathematical theory, such as (nonlinear) semi-group theory (see [31,133]).

D. Numerical Experiments The numerical experiments presented below have been considered in [138,139,148] and illustrate the behavior of different regularization strategies. For more details on the numerical implementation we refer to these papers. Figure 11 shows three common test images and a noisy variant of each of them: an outdoor scene with a camera, a magnetic resonance (MR) image of a human head, and an indoor scene. Gaussian noise with zero mean has been added. Its variance was chosen to be a quarter, equal to, and four times the image variance, respectively. We applied linear and total variation regularization to the three noisy test images, used 1, 4, and t6 regularization steps, and varied the regularization parameter until the optimal restoration was found. Discretizing stabilized total variation regularization with g(x)

-

4'/32+ x

leads to a nonlinear system of equations. The system of nonlinear equations was solved numerically for /3--0.1 by combining convergent fixed point iterations as outer iterations [62] with inner iterations using the GaussSeidel algorithm for solving the linear system of equations. The results are shown in Figures 12 and 13. Figure 14 shows BV-denoised and rendered 3D ultrasound data. This gives rise to the following conclusions (Figure 15): 9 In all cases total variation (BV) regularization performed better than Tikhonov regularization. As expected, total variation regularization leads to visually sharper edges. The BV-restored images consist of piecewise almost constant patches.

DENOISING AND INVERSE PROBLEMS

469

FIGURE 11. Test images. Top left: Camera scene. Top right: Gaussian noise added. Middle left: Magnetic resonance image. Middle right: Gaussian noise added. Bottom left: Office scene. Bottom right: Gaussian noise added.

9 In the linear case, iterated T i k h o n o v regularization p r o d u c e d better restorations t h a n noniterated. Visually, n o n i t e r a t e d regularization results in images with m o r e high-frequency fluctuations. I m p r o v e m e n t s caused by iterating the regularization were mainly seen between 1 a n d 4

470

OTMAR SCHERZER

FIGURE 12. Optimal restoration results for Hl-regularization (cf. (22)). Top left: Camera, 1 iteration. Top right: Camera, 16 iterations. Middle left: MR image, 1 iteration. Middle right: MR image, 16 iterations. Bottom left: Office, 1 iteration. Bottom right: Office, 16 iterations.

iterations. Increasing the iteration number to 16 hardly leads to further improvements. 9 It appears that the theoretical and experimental results in the linear setting do not necessarily carry over to the nonlinear case with total variation regularization. For the slightly degraded camera

DENOISING AND INVERSE PROBLEMS

471

FIGURE 13. Optimal restoration results for total variation regularization. Top left: Camera, 1 iteration. Top right: Camera, 16 iterations. Middle left: MR image, 1 iteration. Middle right: MR image, 16 iterations. Bottom left: Office, 1 iteration. Bottom right: Office, 16 iterations.

image, iterated regularization performed worse than noniterated regularization. For the MR image, the differences are negligible, and the highly degraded office scene allows better restoration results with iterated regularization.

472

OTMAR SCHERZER

FIGURE 14. Bounded variation seminorm denoising with g(t) = v/t + f12 and fl = 0.001 of three-dimensional ultrasound data (top). The left column shows the renderings for noniterated, the right column for iterated regularization. The regularization parameter for iterated regularization was ot = 2.

IV. MUMFORD--SHAH FILTERING The M u m f o r d - S h a h filtering technique has been proposed in [118] for simultaneous filtering and edge detection of noisy piecewise continuous data. Since then the M u m f o r d - S h a h technique has received considerable interest, theoretically, due to the challenging m a t h e m a t i c s involved (see, e.g.,

DENOISING AND INVERSE PROBLEMS

473

FIGURE 15. Results for the MR image from Figure 1 l(a) with noniterated and iterated regularization (fl = 0.001). The left column shows the results for noniterated, the middle column for iterated regularization. The images in the right column depict the modulus of the differences between the results for the iterated and noniterated method.

[6-10,115] and references therein), for segmentation applications, and its numerical implementation (see, e.g., [29,30,38,39,59,60]). Formally the Mumford-Shah segmentation model looks like a regularization functional, and consists in minimizing the functional

f(u,K)'-f(u-ua)2-k-(otlf~

]Vu[2 -+-c~27"/1(K)).

(27)

\K

Here, IYl > 0, Ot2 > 0, and K is the discontinuity set (edges and corners) of u, which is assumed to be of finite one-dimensional Hausdorff measure, i.e., 7-[l(K) < cx); f2 \ K denotes the set fa excluded by the discontinuity set K. For instance for a rectifiable curve, the one-dimensional Hausdorff measure is the length of the curve. The minimizer u :-- u~,,~2 is the filtered data with discontinuity set K := K~,,~2.

474

OTMAR SCHERZER

The functional 0/lf~\g [Vu[2§ serves as functional and is designed to simultaneously penalize for

a

penalization

(1) high oscillations of the filtered data u outside the discontinuity set and (2) complex (long) discontinuity sets K. For the numerical minimization of f(u,K) a common tool is to use non-local approximations, like the Ambrosio-Tortorelli approximation (see [9-11]), where the minimizer off(u, K) is approximated by the minimizer of the functional f ( u , W) "--

fofofo(') (U -- US)2 § 0/1

wZ]Vu] 2 § 0/2

fllVw]2 + ~(1 - w) 2 . (28)

This functional is minimized with respect to (u, w) E H1(~2) • Hi(g2) and no longer involves tedious minimization over a family of discontinuity sets K. The functional (1/fl)f~ ( 1 - w)2 in (28) penalizes for w r 1. Eventually, for fl ~ 0+ this term becomes dominant and the set where w ~ 1 becomes one-dimensional, e.g., a curve with finite length. For fl --+ 0§ the set {w ~ 1} eventually becomes the discontinuity set K of the minimizer of the Mumford-Shah functional (27). The minimizer (u, w) of the functional (28) satisfies the optimality condition, which is a system of coupled partial differential equations (U -- U8) -- 0/1V" (w2Vu) = 0, w]Vul2

-

0/2fl AW 0/1

0/2

0 / l f i ( 1 - w) -

0,

(29)

together with homogeneous Neumann boundary data for both u and w. Figure 16 shows some numerical simulations for Mumford-Shah segmentation and filtering by solving the system of coupled differential equations (29).

V. REGULARIZATION AND SPLINE APPROXIMATION

So far, the regularization models have been presented in an infinite dimensional setting. In this section we review a relation between regularization and cubic spline approximation, by using a semi-infinite dimensional setting.

DENOISING AND INVERSE PROBLEMS

475

FIGURE 16. Top: Test data u8. Bottom left: u solving (29). Bottom right: w approximating the discontinuity set.

S u p p o s e u ~ is a s m o o t h function on 0 _< x <_ 1 a n d noisy samples u i of the values u~ are k n o w n at the points o f a u n i f o r m grid A={0--x0
X i

<

...


be the mesh size o f the grid a n d suppose

lu~ -- uO(xi)l ~ 3,

(30)

where 3 is a k n o w n level o f noise in the data. F o r the sake o f simplicity of p r e s e n t a t i o n we assume that the b o u n d a r y d a t a are k n o w n exactly:

U~o--u~

and

u 8n _ u0(1).

476

OTMAR SCHERZER

We are interested in finding a smooth approximation Ou/Ox of Ou~ in (0, 1), from the given data u 8i. To make the computations concrete we have to quantize the terminology "smooth," which we characterize by the size of the second derivative, i.e., we consider a function to be smooth if the second derivative is small. Then, this approximation problem can be formulated as a constraint optimization problem. Problem 5.1 Minimize fo (OZu/Ox2) 2 among all smooth functions u satisfying u(O) = u ~ u(1) = u ~(1), and 1

n-1

Z n--li=

(U~ -- U(Xi))2 _~< 32.

(31)

1

Then, take the derivative Ou,/Ox of the minimizing element u, as an approximation of Ou~ In fact, given the uncertainty in the data, all functions u satisfying (31) can be considered as solution candidates. The minimizer of Problem 5.1 is the particular candidate that is "smoothest." If the minimizing element u, of Problem 5.1 satisfies the constraint (31) with strict inequality (i.e., the constraint (31) is inactive) then u. (x) = u ~

+ x(u ~

(32)

- u~

i.e., it is the straight line interpolating the two boundary values. This case occurs if and only if u, satisfies the constraint (31). Excluding this trivial case, the minimizer u, satisfies (31) with equality and hence can be calculated using the method of Lagrange. If 1/or denotes the corresponding Lagrange multiplier for constraint (31), the equivalent formulation of Problem 5.1 is: Problem 5.2 Minimize

9__

f (u)

1 n-1 s n - 1 Z (u~ - u(xi)) 2 + ot i=1

(02U~ 2 \OxZ ]

among all smooth functions u satisfying u(0) = u~ u(1) = u~ such that the minimizing element u~ of (33) satisfies 1

(33) where ot is

n--1

n--li= 1

C-

u (xi)) -

(34)

DENOISING

AND INVERSE PROBLEMS

477

The derivative Ou~/Ox of the minimizing function u~ is then an approximation of Ou~ The model essentially differs from the regularization models considered in Section III, since here discrete sample data are available. Also the scope is to find a smooth approximation of the derivative of us and not just a denoised approximation of u s as in the regularization methods considered in Section III.B. Problem 5.2 is a special instance of Tikhonov regularization. The way of choosing the regularization parameter ~ in Problem 5.2 is called the discrepancy principle [87]. Except for the interpolatory constraints at the boundary of the interval, (33) has been investigated and solved by Schoenberg [151] and Reinsch [142], who showed that the solution of Problem 5.2 is a natural cubic spline over the grid A. Reinsch also gives a constructive algorithm for calculating this spline. A more comprehensive level on the interaction between cubic spline approximation and numerical differentiation can be found in [92], see also Hanke [90]. The interaction between regularization and spline approximation is not limited to cubic splines. In fact it can be shown that the optimal solution u~ of the functional with m = 1 , 2 , . . . i67/

is a combination of B-splines of order n = 2m - 1,2 i.e., Uot(X) -- ~

blifl 2m-1 (X -- k ) ,

k~Z

where fl2m-l(.) denotes the B-spline of order 2 m - 1, which is defined as follows fl2m-1 __ ~ O , f l O ,

... ,/~0, .r 2m times

where 9 denotes convolution and

/~o_ 2Z denotes

1,

--12 ~ x ~ 1,

89

I x l - ~ ,1

O,

otherwise.

the set of integer n u m b e r s , - . . . . - 1, O, 1 . . . .

478

OTMAR SCHERZER

A survey on the interaction of splines and regularization can be found in Unser [168]. The topic of numerical differentiation has been studied extensively before (see, e.g., Torre and Poggio [165], Murio [119], Groetsch [82-84]). For more background on the topic of spline approximation we refer to Schoenberg [151], Reinsch [142], Schultz [152], Strang and Fix [157], de Boor [58], Schumaker [153], Rice and Rosenblatt [143], and Wahba [169] to mention a few. VI. SCALE-SPACE METHODS FOR INVERSE PROBLEMS As has been shown in [147] the concept of diffusion filtering cannot be used directly for the solution of ill-posed operator equations. The argumentation is outlined below. For the moment we restrict our attention to Tikhonov functionals defined on the Sobolev space of differentiable functions Hi(g2) where the stabilization term [lu- u, ll2 in (21) is replaced by

f g(Ivul). Then, arguing as in Section III.C the minimizer u~ of (22) satisfies for any h 6 Hi(f2) and any real number t f (u~ + t h ) - f

(u~) > 0,

which is equivalent to

+ ot f~2 g(lVua]2 -+"2t(Vu~ 9Vh) + t2]Vh] 2) - g(lVuot] 2) > 0. A Taylor series expansion of (f(u~ + th) - y~)2 gives

(~'(Ua -4- th) - y8)2 = (.~L"(Uot)- y8)2 -4- 2 t ( f ' ( u ~ ) h ) ( f (u~) - yS) + O(t2). Then by similar arguments as in Section III.C we find

479

D E N O I S I N G A N D INVERSE P R O B L E M S

Using Green's formula, we get

0 - f~ h.Y"(u~)*2(.y'(uo~)-ya)_ ot f~ V. (g'([Vu,~]2)Vu,~)h Jr-Olfofa hgt ([Vu~176

" 1J).

Here f"(u~) .2 denotes the L2-adjoint of f"(u~), i.e.,

f f'(u~)*2(v)w- fu vT(u~)(w)

for all v E L2(~), w E Hi(f2).

This shows that the optimality criterion for the minimizer u~ of (21) is .T'(u~)*2(.T(u~) - ya) - otV.

OUot= 0 0 n Ov

(g'(lVu~lZ)Vu~)

on f2,

0f2.

In the case of noise-flee attainable data, that is for 9 - y O _ we have T ( u ~ ) , 2 (f'(u,~)

~-(ut),

-ot "Y'(u*)) - V . (g'([Vu~[2)Vu,~)

and there exists an associated diffusion-type methodology

.T'(u~) .2 H(u)--O'" _ V. (g'([Vul2)Vu) on Ot Ou = 0 on (0, oo) x Of2, Ov u(0) - u t on f2.

(0, ~ ) x f2,

(35)

Due to the ill-posedness of the operator Equation (19) there will generally not exist a solution of (19) when y0 is replaced by ya ~ y0. The ill-posedness thus prohibits an a priori estimation of an approximation of u t. Thus method (35) is inappropriate for calculating a scale space of an inverse problem. The relation to diffusion filtering becomes apparent if we use 9 ~'-I, 9 Y - L2(~), the space of square integrable functions, 9 X - Hi(f2), the Sobolev space of weakly differentiable functions, and 9 the Hl(f2)-seminorm for regularization, that is g(x) - x.

480

OTMAR SCHERZER

In this setting the minimizer u~ of the Tikhonov functional satisfies ld a ~

U8

= Au~ on ~2, Ol Ouu

~=Oon

Ov

0ft.

Thus for ot --+ 0 + the diffusion filtering Equation (1) is approximated. The iterative Tikhonov-Morozov method is a variant of Tikhonov regularization for solving inverse problems. This method consists in iteratively minimizing the sequence of functionals

f(k)(u) "-- II~(u) -- Yall2 -+-OtkllU-- u
k -- 1 , 2 , . . .

(36)

and denoting the minimizer by/d (k). If the functionals f(k) are convex, then the minimizers U(k) satisfy

.T'(u)*XY(.T(u (k)) - ya) + Otk(U(k) -- u (k-l)) -- 0, k - 1 , 2 , . . .

(37)

Here .T"(u) *xY denotes the adjoint of .T'(u) with respect to the spaces X and Y, that is

(.T'(u)*XY(v), W)x= (v,.T'(u)(w)) r for all v ~ Y, w e X. Typically in the T i k h o n o v - M o r o z o v method one sets u (~ - 0. But any other choice is suitable as well. For example, a priori information on the solution m a y be incorporated in the initial approximation u {~ Taking Olk - - 1 / ( t k - tk-1) shows that U(k) and U(k-l) c a n be considered as approximations of the solution u of the asymptotic Tikhonov-Morozov filtering technique Obl

at

= -Y"(u) *X~(Y'(u) - ya) in (0, e~) x ~,

(38)

u(0, .) - u {~ - 0 in ~2.

F o r . T " - I, the embedding operator from H I ( ~ ) into L2(~), the iterative T i k h o n o v - M o r o z o v method, where we use the H l - s e m i n o r m for regularization instead of the full norm, generates minimizers u (k) of the functionals

f(k)(u) "-- s (U-- Ua)2 +Otk s

IVu-- Vu(k-1)[ 2, k - 1, 2 , . . .

(39)

481

DENOISING AND INVERSE PROBLEMS

Accordingly, the asymptotic Tikhonov-Morozov solving the differential equation of third order u-

0U

u s - A-3- 7 in (0, ~ ) 0U

av

method

consists

in

x f2,

= 0 on (0, oe) x Of2,

u(0, .) - 0 on S2.

(40)

Figure 17 shows the evolution of the solution of the differential equation (40). It starts with a completely diffused image; at t = e~ the input data is restored. In analogy to scale-space theory (cf. Section II.D) we call this method the inverse scale-space method, since it generates a data representation at a continuum of scales, embedding the input data u s into a family of gradually simplified versions initialized with a totally blurred imaged. In Section VI.A we discuss the asymptotic Tikhonov-Morozov method for deblurring images. In this case, )r is a linear integral operator. For this particular model problem we can motivate preferences of different numerical methods in inverse problems and image processing.

A. Deblurring with a Scale-Space Method We consider a problem of deblurring data to recover a function u t on ~2 - (0, 1)2 given (blurred) data t"

yS _ f u t + noise "-- Ja k(]. - y l ) u t (y)dy + noise on ~2. To formulate the Tikhonov-Morozov method we have to specify a similarity measure for the data and an appropriate function space containing u t. In this section we restrict our attention to those u* in one of the following three spaces: (1) The Sobolev space Hi(f2), that is the Hilbert space of weakly differentiable functions u that satisfy ]]UI[H, "--

(fo

]VU] 2 -+- colU] 2


with an appropriate positive weighting parameter co > 0.

482

OTMAR SCHERZER

FIGURE 17. The result of nonstationary regularization (39) and regularization parameters cq--100,000, o~2-- 50,000 , c~3=25,000, ot4---12,500, ot5=6250, c~6=3125, oe7=1500, ot8=750, 0(9 ~- 3 0 0 , O~10 ~-- 150, Otll -- 75, and oq2 = 30 (the last image is not visually different to the input data).

(2) The more general Banach space Wl,P(f2), with p > 1, of functions u satisfying Ilullw,,p " -

(fo)lip IVulp +

colul p

< c~,

with an appropriate positive weighting parameter oJ > 0. (3) The space B V(f2) of functions of bounded variation. That is the class of functions u satisfying [lullBv(a) := f~ (IVul + oJlul) < e~.

483

DENOISING AND INVERSE PROBLEMS

For a function u E BV(S2) the term fa measure (see [73]).

IVul has

to be understood as a

An appropriate choice for the similarity measure is the L2(S2)-norm. Depending on a priori information on u t it is instructive to study the T i k h o n o v - M o r o z o v method in a variety of settings. 9 If u t 6 Hi(f2), it is appropriate to consider r as an operator from Hl(~2) into L2(S2). Accordingly, the iterated T i k h o n o v - M o r o z o v method consists in minimizing

f ~ ) ( u ) "--117u -- y~]]L2(a) 2 + Otkllu - u (k- 1) II2~(a)"

(41)

Instead of the H I - n o r m the H l - s e m i n o r m can be used if ~ does not annihilate constant functions. In particular, for denoising images, that is if ~" = I the Ha-seminorm is suitable. For ill-posed problems, such as deconvolution problems, this seminorm may lead to some numerical difficulties. 9 For u t 6 WI,p(f2), p > 1, the corresponding T i k h o n o v - M o r o z o v method consists in minimizing the functional (k)

fwl,p(U ) "--][.f'U

--

y~

1122(~'2)+ OtkllU -- U(k-l) l]pW,,p(~ ).

(42)

9 For u t 6 B V(f2) the T i k h o n o v - M o r o z o v method consists in minimizing

(k) y~ 2 1) f Bv(U) "--IIUU Ilta(a) + Otk IlU- u
(43)

Since the operator 9c is self-adjoint on L2(f2), that is , ~ " ' 2 -_.- , ~ " the asymptotic T i k h o n o v - M o r o z o v method in the Ha-setting reads as follows Oil

(3CTu)(t, x) - (S'ya)(x) - (A - coI)-z-: (t, x) for (t, x) ~ (0, c~) x f2, ot

Olg

--(t,x) Ov

- 0 for (t,x) ~ (0, cxz) x of 2,

u(0, x) = 0 for x 6 f2.

(44)

The minimizer u (k) ~"J "~ ,~(k) Wl,, has to satisfy

.~(3cu(k) _ ya) _ otkP

v.

_ u

- - Otk P -~Co Iu(k) -- U ( k -

1)) 1 -2V(u k

_

l) IP-Z(u -- u(k-1)).

(45)

484

OTMAR SCHERZER

Introducing the relation 2 1 otk -- P (tk -- tk-1)p-1

(46)

between the regularization parameters and the time discretization we derive the asymptotic Tikhonov-Morozov method on WI,P(f2): p-2

7(7u-

yS) _ v .

(47)

v

For p = 1 the relation (46) degenerates, indicating that there is no asymptotic integro-differential equation for the Tikhonov-Morozov method on B V(f2). One of the most significant differences between diffusion filtering and iterative Tikhonov-Morozov regularization is that a small timestep size in the diffusion filtering method results in very large regularization parameters. This is not inconsistent with standard regularization theory since we consider an iterative regularization technique which uses the information of the previous iteration cycle. In our numerical simulations an exponentially decreasing sequence Otk for the iterative regularization algorithms (41)-(43) leads to a visually attractive image sequence. This, in turn, implies that the time steps tk of the diffusion filtering method (47) are exponentially increasing. This compensates for the fact that in the beginning the diffusion process is rather strong and a small step size is required. As the diffusion progresses the image starts to stagnate and a large timestep size becomes appropriate. B. Numerical Simulations

The following test cases have been considered in [147]. We discuss the numerical implementation of the asymptotic Tikhonov-Morozov method and present some numerical simulations for deblurring images. In the numerical simulations presented below we have used the kernel function (t 2 -- 82) 4

k(t) -

e8

for t ~ [-e, e]

and

k(t) - 0 otherwise.

For the numerical solution of the integro-differential Equation (44) we discretize in time and use a finite element ansatz of products of linear splines on f2. Let N V(tk, X l , X2) -

~ c i j ( t k ) v O ( X l , X2) id=O

DENOISING

AND INVERSE PROBLEMS

485

be the approximation of the solution of (44) where

Uij(Xl, X2)

=

Ui(Xl)Uj(X2)

and Ui is a spline of order 1, that is v i ( j / n ) = 60- for i = 0 , . . . , N and vi is piecewise linear on [0, 1]. For the approximation of the time derivative of v we use a backward difference operator, that is

v(tk, X) -- v(tk-1, X) OV ~, --(tk, X). t k - tk-1 at Using otk = 1/(tk -- tk-1) the discretized system for an approximation of (44) at time tk requires solving the following linear equation for the coefficients co.(tk+l) from given coefficients cij(tk)

~ciy(tk+1)(Fo"kl + ~k;Z~'kl) -- L yam(vkvt) + ~k~Ciy(tk);ZO.~,kt O O"

(48)

for all l,k 9 {0,..., N}. Here

and

'

'

O,kl

The solution of the unregularized Equation (48) (that is with Ofk --" 0 ) is illconditioned. This becomes clear when the singular values of the matrix F are plotted (cf. Figure 18); most of the singular values are comparatively small. Errors in components of the data corresponding to singular functions with singular value near zero are then exceedingly amplified. Thus, it is prohibitive to calculate the solution of the unregularized equation. Example VI.1 In the first example we aim to reconstruct the pattern (top image in the first row of Figure 19) from the blurred and additionally noisy data (cf. Figure 19). Figures 20-22 show the inverse scale-space method for reconstructing the pattern from blurred data. When the blurred data is additionally distorted with Gaussian noise the ill-posedness of the problems

486

OTMAR SCHERZER

becomes apparent. Only for a relatively short period of time is the reconstruction visually attractive. For t--+ cx~ the reconstruction becomes useless. This effect is more significant the more error we have in the data as a comparison of Figures 20-22 shows. One of the major concerns in regularization theory is the estimation of appropriate regularization parameters needed to stop the iteration process before the image becomes hopelessly distorted by noise. For some references on appropriate stopping rules for the Tikhonov-Morozov method we refer the reader to [85,88,91,145].

Example VI.2 Here we aim to compare the Tikhonov-Morozov method on Hi(f2) and B V(f2). We have chosen a piecewise constant function on a rectangle as a paradigm of a function that is in B V(f2) but not in Hi(f2) (cf. Figure 23).This has the effect that the reconstruction with the (asymptotic) Tikhonov-Morozov method on Hi(f2) always has a blurry character (cf. Figure 24). Figure 25 shows the reconstruction with the Tikhonov-Morozov method on B V(f2). This method performs worse than the asymptotic Tikhonov-Morozov method on Hi(f2). This numerically supports the fact that there is no inverse scales space method on B V(f2). This section has been devoted to highlighting the controversial behavior of scale-space methods for the solution of inverse problems and image smoothing and restoration. One of the significant differences in inverse scale space theory for inverse problems is the choice of an adequate stopping

x 104 1.5

1

0.5

0

0

100

200

300

400

500

600

700

800

FIGURE18. The singular values of the matrix F.

900

DENOISING AND INVERSE PROBLEMS

487

FIGURE 19. Top: The test pattern. This pattern is aimed to be recovered from the blurred data (middle left), the blurred data which is additionally distorted with medium noise (middle right), and distorted with high noise (bottom).

488

OTMAR SCHERZER

FIGURE 20. Reconstruction from blurred data without noise by the inverse scale-space method (44). The images show the solution u of (44) at specified time with exponentially decreasing time-steps. At a certain time the test pattern can be completely recovered. The inverse scale-space methods stagnated at the test pattern.

DENOISING AND INVERSE PROBLEMS

489

FIGURE 21. Reconstruction from blurred data with medium noise using the inverse scalespace method (44). The images show the solution u of (44) at specified time. Top left shows the optimal time for recovery with medium noise. After that time the reconstruction gets worse, showing the importance of determining an optimal stopping time for the inverse scale-space method.

490

OTMAR SCHERZER

FIGURE 22. Reconstruction from blurred data with high noise using the inverse scale-space method (44). Middle right shows the optimal time for recovery with high noise. After that time the reconstruction algorithm diverges extremely fast (cf. the scales of the images).

491

DENOISING A N D INVERSE PROBLEMS o o3s

o o25

so

FIGURE 23. Test-data for comparing the Tikhonov-Morozov method on Hi(f2) and

BV(S2). Left: Image to be reconstructed. Right: The available blurred data, from which we intend to recover the left image.

x 10 4

3L,

30

30

30

FIGURE 24. Reconstruction with the asymptotic Tikhonov-Morozov method on H I ( ~ ) at specified time.

492

OTMAR SCHERZER

x 10 -~

x 10 "4

~ ~ i!~ii,!~! ......... iiiii!il ~ iiiili~i i ~i~ ~ii~~,~i ,~i , ~'~!i!'! ii

i~i,~,~ ~

~'"~ .... ~"~'~!!i!ii!!iii~'~'!!!!~,~ ~,~!~i ~'~I~i~ ~!'~!'!i!~!ii!!!ii!? ~~'

FIGURE 25. Reconstruction with the asymptotic Tikhonov-Morozov method on specified time.

BV(f2) at

DENOISING AND INVERSE PROBLEMS

493

time; after a certain time the noise is considerably amplified in the reconstruction. This is not an issue in image smoothing, where the effect of noise is weakened over time. We also remark that for image smoothing and restoration the total variation flow filtering in almost all documented cases performed significantly better than the heat equation. This is not always true for inverse problems. VII. NONCONVEX REGULARIZATION MODELS In Section III.C we considered regularization functionals of the general form

f (u) "-- fe(u -- ua)2+ fe g(Vu). The existence of a minimizer is relatively easy to establish under the essential assumptions that g is convex with respect to the gradient variable 7u and the functional is coercive (see, e.g., Dacorogna [54,55] and Aubert and Kornprobst [19]). The analysis of regularization functionals becomes considerably more involved if the functional f is nonconvex. Such models are outlined below.

A. Perona-Malik Regularization In the classical Perona-Malik filter [134,135] we have D(Vu) =

1 1 + ]Vu] 2 "

The corresponding variational technique consists in minimizing the functional

f~ (u-uS)2-']-Ot lift ln(1 +

'Vu'2)-fo

The function ~ is nonconvex with respect to the variable Vu. In this case it is well known from the calculus of variations (see, e.g., [55]) that the optimization problem is not well-posed in the sense that there need not exist a minimizer. Therefore, additional regularization concepts are involved in the functional, such as

fR_pM(U)'--fa(u--uS)2+, fa

ln(1 + IVL•

where L• is a linear convolution operator with a smooth kernel.

(49)

494

OTMAR SCHERZER

The minimizer of the regularized Perona-Malik functional satisfies

(50)

i uIV The corresponding nonlinear regularization technique is

diffusion

Ou_ LyV. at -

(

process

VL•

1 + IVL•

associated

)

2 "

with

this

(51)

Regularized Perona-Malik filters have been considered in the literature [22,37,130,171,172]. Catt6 et al. [37], for instance, investigated the nonlinear diffusion process

Ou_ a / - - V. 1 + i-~•

2 .

(52)

This technique (as well as other previous regularizations) does not have a corresponding formulation as an optimization problem. In an experiment we juxtapose the regularizations (51) and (52) of the Perona-Malik filter. Both processes have been implemented using an explicit finite difference scheme. The results using the M R image from Figure 11 are shown in Figure 26, where different values for V, the standard deviation of the Gaussian, have been used. For small values of V, both filters produce rather similar results, while larger values lead to a completely different behavior. For (51), the regularization smoothes the diffusive flux, so that it becomes close to zero everywhere, and the image remains unaltered. The regularization in (52), however, creates a diffusivity which gets closer to one for all image locations, so that the filter creates blurry results resembling linear diffusion filtering.

B. Relative Error Regularization The noise in data detected with common measurement devices frequently correlates with the exact data. Here, relevant situations are when the noise locally correlates with the amplitude or the variation of the data. Assuming correlation between the data and noise we are led to fit-to-data terms of the form

l f (u- uS)2

l f lu-u~l p

lul-----7--, p

IVulP- ~

p-

1,2,...

DENOISING AND INVERSE PROBLEMS

495

FIGURE 26. C o m p a r i s o n of two regularizations of the P e r o n a - M a l i k filter (t = 250). T o p left: Filter (51), y - 0.5. T o p right: Filter (52), y - 0.5. Middle left: Filter (51), y - 2. Middle right: Filter (52), y -- 2. B o t t o m left: Filter (51), y -- 8. B o t t o m right: Filter (52), y - 8.

496

OTMAR SCHERZER 1.5

s!

.

.

.

.

.

.

.

.

FIGURE 27. Correlated noise in 1D signals. Top left: Noise-free data. Top right: Uncorrelated noise. Middle left: f.[(u-u~)2/lu[]--t~ 2. Middle right: f~[(u-us)E/Iu[ 2] = ~2. Bottom left: f.[(u - uS)2/lVu[] -- t~2. Bottom right: f.[(u- uS)4/lVu[ 3] -- ~2.

In Figures 27 and 28 we have plotted noisy data revealing the difference between uncorrelated and correlated noise. We concentrate on Tikhonov type regularization models with BV-seminorm stabilizing functional. This leads to regularization models of the form (relative error regularization)"

lfo,uu,,

lul---------7---+ a

fo

IVul;

p-

fouu, p fo ~iVulp-1 + a

IVul

In order to put this work into context with diffusion filtering techniques it is convenient to consider iterative relative error regularization.

497

D E N O I S I N G AND INVERSE PROBLEMS

FIGURE 28. Correlated noise in 2D signals. Top left: Noise-free data. Top right: Uncorrelated noise. Middle left: fa[(u-uS)2/lul] = 32; Middle right: fa[(u-uS)Z/[ul 2] = 32; Bottom left: f~[(u- u3)Z/[vu[] = 32; Bottom right: fa[(u u6)4/IVu[ 3] = t~2. -

In particular, we consider the models of iteratively minimizing the functionals:

-~

lulP

+ c~

fo

IVul,

P

[VulP- 1

+c~

fo

IVul

(53)

and denoting the minimizers (presuming they exist) by u(k); moreover, we again use the convention u (~ := u 8. Since the functionals in (53) are nonconvex and thus quite delicate to handle analytically and numerically, it is convenient to consider semiimplicit variants such as the models of minimization

lfo,u_u l, fo

~

]u(k_l)]p

+or

IVul;

l/o,u_u l , fo

p-

]Vu(k-1)] p-1 +or

[Vul

(54)

The functionals in (54) are convex and straightforward to analyze (see [146]).

OTMAR SCHERZER

498

Minimization of the second functional in (54) with p = 2 can be considered as a semi-implicit time step with step size At = 0/ for the mean curvature flow equation (9); for p = 4 it is a semi-implicit method for solving the affine invariant mean curvature flow equation

0U

((V]_~~])) 1/3

a t = IVul v .

,

(55)

The first functional in (54) corresponds to a semi-implicit time step for solving

a t = lulPV

"

(56)

The Euler equation for the minimizer of 1 L lu - u(k-1)12 L IVul + o~ IVul

IVul

-

v.

a

2

[Vu[ 2

(57)

,] ~

(58)

"

Note that (58) is only formal since the regularization functional fa IVul is not differentiable. Division of the equation (58) by 0/gives

U-- u(k-1) ~--IVulV. 0/

((1

1 (U-- u(k-1))2 0/ ) VI-~UU]) 2

0/2

IVul 2

Taking the formal limit 0/-+ 0+ and considering again U(k) ~" U(tk) and 0/= t k - tk-1 gives again

"

U(k-l) ~

U(tk-1),

Ot the mean curvature flow equation. Since (58) can be considered to be a Perona-Malik model with positive and negative diffusion, the solution is illposed. The ill-posedness in the optimality condition reflects the fact that the underlying energy functional (57) is nonconvex with respect to the gradient

DENOISING AND INVERSE PROBLEMS

499

FIGURE29. Original image (top) and filter images: mean curvature flow (middle left); affine mean curvature flow (middle right); implicit regularization (bottom left); B V regularization (bottom right). variable. By employing generalized solution concepts such as convexification or F-limits the ill-posedness (see [146]) disappears. We present two numerical experiments for relative error denoising: (1) We use the artificially generated data set at the top left of Figure 28. The several reconstructions in Figure 29 have been created with bounded variation regularization, mean curvature filtering, affine mean curvature filtering, and implicit error regularization. The stopping time in the diffusion filtering method and the regularization parameters are selected such that all reconstructions have about the same amplitudes. (2) The second example is concerned with denoising of ultrasound data sets.

500

OTMAR SCHERZER

FIGURE 30. Original image (top) and filter images: mean curvature flow (middle left); affine mean curvature flow (middle right); implicit regularization (bottom left); B V regularization (bottom right).

From the numerical reconstructions one finds that mean curvature flow and implicit error regularization produce very similar results if the regularization parameter and the diffusion time are identified (Figure 30).

V I I I . DISCRETE BV REGULARIZATION AND TUBE METHODS

So far we have presented regularization models in infinite dimensional settings (cf. Section III) and in semi-infinite dimensional setting (cf. Section V). In this section we concentrate on completely discrete settings for bounded variation regularization.

501

DENOISING AND INVERSE PROBLEMS

The derivation of discrete variants is not straightforward. Several numerical realizations of discrete bounded variation regularization can be derived, some of which are outlined below. For a piecewise constant function g in [0, 1] of the form

g-

gi in f2i "- ( i -

' i),

i-1,...,n

(59)

we define Tg

"-

((Tg)i)i=o

..... n '

where (Zg)i

._

--

gin t- gi+l for i -- 1 . n -- 1 2 " '

(Tg)o := gl, (Tg)n :-- gn. We call Tg the traces of g. A piecewise constant function and its traces are plotted in Figure 31. Using these ingredients we are able to formulate two discrete variants of BV regularization. We restrict our attention to minimization of functionals over the set of piecewise constant functions n

S "-- {U" U ( X ) -

~CiXf2i

with Icil < ~ } ,

(60)

i=1

where Xa/denotes the characteristic function of the interval

,

o

,

o.,

~

I

o.~

i

o.,

o'.~

o:,

i

o.,

o:,

o;

~"~i.

!

,

FIGURE 31. A piecewise constant function with values gi, i = 1 . . . . . n and the traces (Tg)g, i = 0 . . . . . n; symbolized by *

502

OTMAR SCHERZER

Discrete bounded variation regularization functionals differ by the way of interpreting the available discrete data. Two possibilities are considered. (1) The data can be interpreted as the measurement data of traces (Tf)(i/n) of a BV-function f. This is a sampling problem. In typical sampling problems one interprets the data as the function value f(i/n). Since in our setting f may be discontinuous at i/n, and point evaluation is not possible, we are forced to use trace evaluation. This leads us to consider minimization of the functional

TBV-d(U) -- 2(n 1+ 1) Zn [(Tu)i - (Tf)il2 -I- ot f01 luxl

(61)

i=0

over S, where (Tf ) i is the given data at i/n. (2) Alternatively to assuming available sampled data, one can interpret them as values of a piecewise constant function n

f

-

Zfix~.2i , i=1

Given measurement data f., i - - 1 , . . . , n, this suggests minimization of the functional

ini~l(Ci.=

TBV_d2(b/) .__ ~

_j,})2 _it_o/

fol luxl

n

----|Z(Ci--fi) 2n/=1

2

-Jf-Ol~--~~lr 1 -- Ci] i-1

(62)

over S. These two possibilities will be utilized below.

A. Discrete B V Regularization (Sampling) The functional TBV-d is well-posed, i.e., there exists a unique minimizer in S (see [95]). To further analyze properties of the minimizer u~ - ~in=l uixi of TBV-d it is instructive to study the optimality condition for the coefficients

503

DENOISING AND INVERSE PROBLEMS

Ui, i = 1 , . . . , n. Setting fl = 4na, the coefficients of u~ satisfy the set-valued equations

(u/ui+,) Ui

(Ui--1 -'[- Ui) "-Jl-(Ui "Jl'-U/q-l) --[- 1~ 0r

Ui-1

_.[..

[ui - - U i - l [

[ui - -

Ui+l [

+ft') Jr- (/~ -'["fi+l)

for i - - 2 , . . . , n 5Ul + u: + ~

1,

Ul -- U2

~ ~ 5fl + f z, lul - u21

5Un -q- Un-1 nt- fl

Un ~ bin- 1

[Un - Un-l [

(63)

~ 5fn + L - l ,

where we use the abbreviation {1} {-1}

if if

e<0,

{[-1,11}

if

C --

e = le]

which is the subgradient of lel. We observe from (63) that for j -

e>0,

1,... , n -

0,

1

j-1

Uj -- b/j+ 1

(Uj -t- Uj+I) + 2 Z

(Ui + Ui+I) + 4Ul nt- fl [Uj -- Uj+I[

i=1 j-1

(J~ +j}+l) + 4fl.

(64)

n-1 (Ui-[-" b/i+l) -[- 4Ul -- 4fn + 2 ~ (f" +ft'+l) -'[- 4fl. i=1 i=1

(65)

(fj" -[-fj+l)

"nt-

2Z i=1

and n-1

4u, + 2 Z Let

j-1

~'u(J'/n) "-- (uj + Uj+l) + 2 ~ ( u i +

Ui+I) -t-

4Ul,

j -- 1 , . . . ,n - 1,

i=1 n-1

/~u(1) "-- 4Un + 2E(ui-1"- Ui+l) + 4Ul. i=1

(66)

504

OTMAR SCHERZER

Moreover, let /~f+, /~f- be linear splines with respect to the nodes { j / n : j = 0 , . . . , n} interpolating

? f (o) - o, F f ( j / n ) - F f ( j / n ) 4- fl j - 1 , . . . , n - 1,

E l ( l ) -/~f(1).

(67)

The region between the two linear splines is referred to as the tube T./~fand/~f+ mark the lower and upper bounds of T. Since uj ~uj-1 <1 luj - Uj-ll we find from (64) and (65) that /~u 6 7-.

(68)

In other words, the antiderivative 3 /~u of the minimizer u of the discrete bounded variation regularization formulation is in the tube 7-. We therefore refer to (61) as a tube method.

B. Finite Volume B V Regularization

Minimization of TBv-a2 as defined in (62) over S is a standard method of formulating discrete bounded variation regularization (cf. Mallat [108]). Note that in this case the data3'} are interpreted as coefficients of a piecewise constant function, while in Section VIII.A the data (Tf)i are interpreted as sampling data. Since the first term in TBv-a2 is strictly convex it is immediate that the functional TBv-a2 has a unique minimizer. To specify the optimality criteria for a minimizer of the functional TBv-a2 let /3 = not. Then the minimizer u~ of TBV-d2 can be represented as

3We refer to fo f(s)ds as the anti-derivative of the functionf.

DENOISING AND INVERSE PROBLEMS blot "-- Y ~ i n l

Ui,)( i with ui

505

satisfying N

( Ui - Ui--1

Ui+I )

Ui-

ft"

ui + fl ]ui - ui-l [ + lui - ui+l l

for i = 2 , . . . , n -

1,

i

Ul -- u2

Ul + ~ ~ g

[Ul -- UZ] Un -- Un-1

Un+~

lu. - Un-l l

fl, (69)

~ fn.

Let J Fu(O) -- O, Ff(j/n) -- ~ f k

forj - 1,..., n

k=l

F f (O) -- O, F f (j/n) -- Ff(j/n) q- fl

f o r j -- 1 , . . . , n - 1, (70)

F?(1) -- Ff(1).

W i t h f w e associate the tube T bounded by the linear splines F f connecting the values F f ( j / n ) . Thus the minimizer u~ has the property that its antiderivative Fu satisfies

Fu ~ T , and thus the finite volume BV regularization is a tube method as well. Again, the antiderivative Fu of the minimizer u is in the tube T, i.e., it is a tube method.

C. The Taut String Algorithm In this section we recall the taut string algorithm (see [57,109]) for denoising discrete one-dimensional data. We choose a description which allows generalization to higher-dimensional data. T denotes the tube from the previous section. The taut string "algorithm" is actually the solution to a minimization problem, which we specify next. Algorithm VIII.1 Let solution to: -

v--ZiL1 l)ix~i

with antiderivative

-

1 + Iv(r)l 2 dr--+ min

V denote the

(71)

n i=1

over all continuous and piecew&e linear functions V on [0,1] with function values in T.

506

OTMAR SCHERZER

Physically speaking V is a string of minimal length contained in the tube T,, connecting (0,0) and (1,F/(1)), i.e., it is taut. In particular in regions between two contact points of V with the boundary of T a taut string is affine linear and v is a piecewise constant function. The values vi approximate the input data fi and v constitutes a denoised approximation to f. In [57] an algorithm for computing the solution to (71) was presented which proceeded iteratively from one nodal value of the tube to the next. The solution method that we shall propose will be completely different. The taut string, as determined from Algorithm VIII.l, and the finite volume BV-regularized solution with c~ = ~/n are both contained in the same 77. This in particular shows that the graph of the finite BV-regularized solution is at least as long as that of the taut string solution. Finite volume BV regularization and the taut string algorithm share the property that they preserve homogeneous regions of the original data. This is easily seen for the taut string algorithm, since in a flat region of the original data f the function FU is linear and consequently the taut string is linear in this region too, showing that the flat regions of the filtered data (i.e., the derivative of the taut string) either correspond with the input data or are enlarged. For finite volume BV regularization (as well as other methods) this statement was addressed with rigor in [128].

D. Multidimensional Discrete B V Regularization In this section we present multidimensional analogs of sampling and finite volume BV regularization. Moreover, we propose a multidimensional analog of the taut string algorithm. Let g2 = (0, 1)• (0, 1) and let f be piecewise constant with respect ~"2/j = ~"2i X ~"~j.

To introduce the sampling BV regularization in ~2 we proceed as in Section VIII.A and model the fit-to-data term as in (61) as the sum of both components in the xl and x2 directions, separately. F o r f - ~ina.=l c0x~,j the BV sampling method consists in minimization of the functional 1

n-1 n-1

Ci+lj -~- Cij

1 n-1 n-1 9 _ _

J~+lj - [ - ~ j

2

.=

9

2

cij -~- Cid+l

2

2

fij + fij+~ + ~ f~ IVul 2

(72)

D E N O I S I N G A N D INVERSE PROBLEMS

507

over

A multidimensional analog of the finite volume BV regularization in a higher dimension consists in minimizing the functional TBV-d2f(U) "-- ~

1 s163 i=1 j = l

jl2+~flVul

(73)

over S. To propose an extension of the taut string algorithm to the case of twodimensional data it will be useful to reconsider the taut string algorithm in the following form: (1) Integration of j~, i - 1 , . . . , n gives a linear spline FU. (2) Determination of the taut string Fu in the tube FU -ot and Ff + or. (3) Differentiation of Fu to obtain the reconstruction for f. Generalization to higher dimensions is impeded by the fact that there is no obvious analog for integration step 1 above. To overcome this difficulty, we proceed by introducing an appropriate potential equation and consider the one-dimensional case first. Given f, we define 9 as a solution to

~xx--f ~x-0

on f2, on{0,1},

and set F f ~ (~X"

This replaces step 1 above. Step 2 can then be realized by solving the contact

problem: (1) Find a function w 6 BV(0, 1) satisfying w(0) - Fu(O) and w(1) - FU(1) (in the sense of traces) that minimizes

fol V/l + w2 subject to the constraint

Ff - ~ <_w <_Fi +~.

(74)

508

OTMAR SCHERZER

(2) Set

u-

Wx.

This approach can be generalized to higher dimensions in a straightforward way: (1) Solve

A~ --f

in f2

~=0

on 0~2.

Ov

Define FS = V * .

(2) In

~2 find two functions

Wi E

BV(f2), i -

1,2, minimizing

/1 + [Vwil2

(75)

subject to the constraints (Ff ) i - o l ~ w i ~ (Ff ) i --[-13l, i - - 1,2,

(76)

and w = Ff on Og2.

This is a contact problem for finding a "minimal surface" in the layer bounded by (Ff ) i - o t and (Ff ) i + or. (3) Set u = V. w. The choice of proper boundary conditions for q~ and w is not obvious. We tested alternatives to our choice and found that they have no significant influence on the numerical reconstruction. In the one-dimensional case we could have chosen the boundary conditions in such a way that they are consistent with Algorithm VIII.1. This choice, however, has no clear multidimensional analog. E. Numerical Test E x a m p l e s

The practical realization of the taut string algorithm in ~2 requires the solution of the bilateral obstacle problem (cf. Algorithm VIII.l), which can be solved efficiently using active set strategies. The particular implementation has been considered in [95], where also general references on active set

DENOISING AND INVERSE PROBLEMS

509

strategies and its numerical implementation can be found. In the following we present some numerical simulations with this method.

1. One-Dimensional Test Example We consider the function

f(x)'-

1

in [0, 1/4)

2 4-4x

in (1/4, 1/2) in ( 1 / 2 , 3 / 4 )

1

in (3/4, 1]

In Figure 32 we display the results of several test runs for the onedimensional example with three different absolute noise levels, i.e., 31 = 0 , ~2~-0.1, and 3 3 = 0 . 5 from left to right. The respective u values are or1 = 1.0 • 10 -5, u2 = 5.0 • 10 -3, and u 3 - 2.5 • 10 -2. In our test, these values produce the best reconstructions. In Figure 32 column i (i = 1,2, 3) corresponds to a test run with (3i, t~i). The first row in Figure 32 shows the .

9 L.~,

.; .

Data

.

.

.~

d,

.~,

Reconstrucllon

.

.

.

.

.

.

.;

.

.~,

~,

Data

.

.:

,

9

.~,

.:,

for a - 1.0E--5

.

.

.~

.',

.',

R ,econst, n ~ t ~ n

u~'

Data

.'.

for , G - 5 . 0 E - 3 .

"

.i,

,~,

~

,

n,.o!se 0.1

i

.

" .

"

i,

"

Reconstruction 1or o . - 2 . 5 E - 2 , , , . . . .

"

.

noise -0.5 , .

!! w

w

w

o,

/

/

",,

", .u,

,, ",\ ,,'

FIGURE 32. One-dimensional tests. Column i (i = 1,2, 3) corresponds to

(t~i, 0//).

,

.,

510

OTMAR SCHERZER Reconstruction for o - 1 . 0 e - 2 ,

Reconstruction for 0-1.0E-2 .

.

.

.

.

.

.

.

.

.

,

,

noise .

= 0.1

Reconstruction for o,=7.5E-2, noise = 0 . 5

I

FIGURE 33. One-dimensional tests. Column i (i= 1,2, 3) corresponds

\

to

(~i,Ofl).

data, the second presents the reconstructions. In the third row we plot the string w (solid) together with its bounds (dashed). In Figure 33 we study the effect of or. We have selected ot values which are larger than the values in Figure 32. In fact, we have ot~- 1.0 • l0 -2, ot~- 1.0 • 10 -2, and ot~ - 7 . 5 • 10 -2. Since we loosen the barriers of the tube the string becomes more flat. We next summarize some of the features observed for denoising of onedimensional images with the taut string algorithm. They are similar to those obtained by nonlinear BV-regularized reconstructions. (1) The mean value of the registered image intensity is preserved by the filtering method. This important feature of diffusion filtering methods (cf. [172]) and nonlinear regularization models (cf. [148]) does not hold for instance for discrete morphological filters, such as the median filter. (2) Spurious noise is removed. (3) Edges are preserved. (4) The taut string algorithm produces images which are damped in height. The magnitude of damping is comparable to that observed for BV-regularized solutions. 2. Two-Dimensional Bench-Mark Problem

Here the solution bilateral contact problem is shown for the bench-mark image in Figure 34 (upper left).

IX. WAVELET SHRINKAGE

In this section we review the interactions of wavelet filtering, diffusion filtering, and variational methods. For this purpose it is convenient to briefly review orthonormal wavelets.

D E N O I S I N G A N D INVERSE PROBLEMS

511

FIGURE 34. Exact data and noisy data in the first row, reconstruction in the second row.

A. Daubechies' Wavelets We review Daubechies' construction of orthonormal wavelets (see [56]). The construction is based on the existence of a scalingfunction ok, such :---- 2-m/ZdP(2-mx --k), k ~ Z, are that for m 6 7/ the functions r orthonormal with respect to the norm on L2([~). Moreover, 4~ is chosen in such a way that for m ~ 7/

Vm 1=

span{qbm,k : k E Z}

:={ZakCm'k'ak~:Of~176 ' k e ~ form a multiresolution analysis on L2(~), that is

Vm C Vm-1,

mET/,

with ~'~ Vm- {0} mET/

and

U mEZ

Vm - L2(IR).

512

OTMAR SCHERZER

The wavelet spaces that is

W m are

the orthogonal complements of Vm in

Vm-1,

Wm "-- V~m n Vm_l. The mother wavelet qt is chosen such that the functions l/tm, k "--

2-m/Z~(2-mx- k),

k ~ Z,

form an orthonormal basis of W m. Since ~b- 4~0,0 ~ V0 c V-l, the scaling function 4~ must satisfy the dilation equation

c~(x) - Z hkc~(2x -- k),

(77)

kEZ

where the sequence {hk} is known as the filter sequence of the wavelet

7r(x) = Z ( - 1)khl_k~(2X -- k).

(78)

k~Z

The filter coefficients have to satisfy certain conditions in order to guarantee that the scaling functions q~ and qt fulfill certain properties. In orthogonal wavelet theory due to Daubechies the desired properties on the scaling functions and wavelets are: (1) For fixed integer N > 1 the scaling function 4~ has support in the interval [ 1 - N , N ] . This, in particular, holds when the filter coefficients satisfy

hk = O,

for k < 1 - N and for k > N.

(79)

(2) The existence of a scaling function q~ satisfying (77) requires that hk -- 2.

(80)

k~Z

(3) In order to impose orthonormality of the integer translates of the scaling function q~, that is f~ c~(x-l)c~(x)dx-3o,t, the filter coefficients {hk} have to satisfy

Z hkhk-2t -- 230,t, l -- 0 , . . . , N - 1.

(81)

k6Z

(4) The wavelet qt is postulated to have N vanishing moments, that is

f

x tqt(x) d x - O ,

1-0,...,N-1

(82)

DENOISING AND INVERSE PROBLEMS

513

which require the filter sequence to satisfy

E(--1)khl-kkt -- 0, l - - 0 , . . . , N -

1.

(83)

k6Z

The oldest wavelet is the Haar wavelet where h0 = hi = 1. In this case the scaling function is

r

1

for x e [0, 1]

0

otherwise

The wavelet function ~ is given accordingly by

ap(x)-

1

for x 6 [0, 1/2),

-1

for x e [1/2,1],

0

otherwise.

l

Since the functions ~[tm,k form an orthonormal basis of L2(N), any function f 6 L2(N) can be expanded in terms of this basis:

f (x) - E f j',k7tj,k(X). j,kE7/

Orthogonal wavelets on L2(N) form the basis to construct wavelets on compact intervals [50,113] (for a summary of this topic we also refer the reader to [49]). A family of orthonormal scaling functions and wavelets on multidimensional domains can be constructed from products of one-dimensional scaling and wavelet functions:

OZ,(x l, x2) - r

(xl)r

(x2),

~),~,(xl,x2) = ~jl,kl (xl )g,j~,k~(x2), l~2]~(xl, x2) -- ~gjl,kl(Xl)r 1/r3,/~(Xl, X2) -- ~j,,k, (Xl)l~rj2,k2(X2), where we use the convention that j = (jl,j2), k - (kl,k2). The functions Cj,~ are called multidimensional scaling functions and the functions ~p~s i = 1, 2, 3 are called multidimensional wavelets.

514

OTMAR SCHERZER B. Denoising by Wavelet Shrinkage

Donoho and Johnstone [63] introduced a wavelet-based denoising algorithm, the so-called wavelet shrinkage algorithm. This algorithm consists in calculating the wavelet expansion (see, e.g., [56]) 3

s163 r

U~f ,, i - - i

x2)

i= 1 (f,/~)EZ2 x 7/2

of the input data and manipulating its coefficients, to be precise the 3'i ] with coefficients u~'/are approximated by 0 ~(U).s j,k

t-or O~(t) "-

0 t+~

t>a Itl _< o~. t <-~

Figure 35 shows the wavelet denoising algorithm with the Daubechies-2 wavelet (with four coefficients h_l, h0, hi, h2) and the Canny edge detector applied to the filtered data (cf. Figure 36). The method of wavelet shrinkage has been paid considerable attention in the literature (see, e.g., [40,42,63,64,114]) and has been applied for the solution of many practically important problems. 1. Relation to Diffusion Filtering

In [40,42] the relation between wavelet shr&kage and regularization methods on the Besov spaces has been established. For the sake of simplicity of presentation we assume that the image data u8 is available on ~2. This is not quite consistent with the overall presentation where we assumed image data on the bounded domain f2 = (0, 1) x (0, 1). In principal one can proceed as outlined below, if instead of wavelets, periodic wavelets are used. The Besov space Bi (L 1(~2)) can be characterized as follows (note that this is not the standard definition): f e BI(LI(~2)) if and only if 3

oo w i t h f } s i=1 (Xs

iR ~fO~s

x7/2

Here ~r/are smooth, orthonormal wavelet functions, and thus f / - denote ,/z j,k the wavelet coefficients of the function f. j~

9

.

DENOISING AND INVERSE PROBLEMS

5 15

FIGURE 35. Wavelet shrinkage with ot = 0, 1, 5, 10, 50, 100.

Formally, the derivative of 4~(f) can be calculated as follows. Let h ~ BI(LI([~2)) f-) L2(N2), then the derivative of 4~(f) in direction h is given by O~(f)(h) - lim ~ ( f 4- t h ) - 4)(f) t--+O t

3

f./.

:

s

hlll'~~7'

which in turn implies that 3 i= 1 (/,/~)~ze • 772

J'--3-k ~ ,

516

OTMAR SCHERZER

FIGURE 36. Canny's edge detector applied to the wavelet shrinkage data with or=0, 1, 5, 10, 50, 100. o

I o

I

where, of course, the meaning off..';/lfj; I.,, is set-valued, as explained in Section VIII.A. Thus the wavelet coefficients u/ . of the minimizer u~ of the regularization functional j,k ,

.,,

,

(84)

12fR2 (u - uS)2 + c~4~(u) satisfy the optimality condition 3

E

i--1 (Z/~7/2 x7/2

(uif, fc -- u~'i] i j,k! ~

~

3 --Ol~ /=1 (/.k~)~7/2 x ~2

u i .

ui_._.

j,k

517

DENOISING AND INVERSE PROBLEMS

Since the functions apt; are orthonormal we find that u/ . j,k j,k 3 --c~ ~ for all (j, k) ~ Z 2 x 7/2, u/ .

i = 1,2,3.

(85)

j,k

This shows that u5 / > a

if u/- > 0,

ua.' / < - a

ifu/-<0,

uj,k ~'(- 0

if lug-.] < or. j,k

j,k j,k

j,k

j,k

Consequently from (85) it follows

d..j,k

u5 ( _ j,k

if u5 / > oe

0

if ua.'/ < or.

uS~ +

if ua.' / < -or

j,k

j,k

j,k

--

(86)

j,k

This shows that Besov space regularization is D o n o h o ' s wavelet shrinkage algorithm. Proceeding as in Section III with the regularization technique (84) the diffusion filtering Ou - - + &k(u) ~ O f o r t > 0 Ot u(0)

=

u~

is associated.

X. REGULARIZATION AND STATISTICS

There has been considerable interest in incorporating statistical a priori information in regularization techniques. In this section we outline the basic principle. There are several publications in the literature devoted to this topic; an extremely useful overview article is [89], where also adequate references can be found. For an elementary introduction to statistics we refer the reader to [98]. Let blS(Xi) = bl(Xi) .Jr- n(xi) be the measured image intensity at the pixel Xi, which is degraded data u with noise n. In the stochastic framework the

518

OTMAR SCHERZER

intensities ua(xi) and n(xi) a r e considered registered intensities of random variables U and N. We denote by P({U < u}) and P ( { N < n}) the probabilities that the random variables U and N are less than u, n, respectively. The probability density functions are accordingly denoted by P({ u ~ [u, u + du)}) du

P({U=u})=

lim du-+O+

P ( { N - n}) -

lim P ( { N ~ [n, n + dn)}) dn--+O+ dn

The notation P({ U -= u}), P ( { N = n}) is typically used in the case of discrete random variables. We find it instructive to use this notation for continuous r a n d o m variables too. The goal is to recover the image intensity u such that the conditional probability P({ U - u} n {N - u 8 - u}) - P({ U - u})P({N - u a - u})

(87)

is maximized with respect to u. Note that the last identity requires that the random variables U and N are independent. If the noise is normally distributed with mean value zero and variance a, then P({N

-- u a -

u}) _ ~ 1

e-[(u'~-u)2/(2~2)]

(88)

~/2rtcr 2 It is convenient to set P({ U = u}) := e "7(u),

(89)

where f" is a nonnegative function. In this case maximization of (87) is equivalent to minimizing

1

"~(U)(Xi) -'~ ~ff2 (blg(Xi) --

U(Xi))2"

(90)

F o r image processing applications it is necessary to take into account neighborhood relations of image intensities between pixels. This can for instance be achieved by using a nonnegative probabilistic model ~ which is dependent on gradient approximations of the image intensity.

519

DENOISING AND INVERSE PROBLEMS

In order to minimize (90) for all pixel values functional

(

9~ ' ( U ) ( X i )

+ ~

,

Xi, i 6 Z,

(U 8 ( x i ) -- U ( X i ) ) 2

)

we minimize the

9

ieZ

To realize the interaction between Tikhonov-type regularization models it is convenient to note that the sum is a quadrature rule approximation of the integral

s

20"2 f ' ( u ) ( x ) -k- (ua(x)

-- u(x)) 2 d x .

Using .T(u) = fllVu[ 2, with fl > 0, the stochastic approach is equivalent to Tikhonov regularization with regularization parameter ot -- 2fl0"2; ~'(u) =/3lVul is bounded variation regularization; f(u)=/3lVulloglVul is entropy regularization. For Tikhonov regularization, .T(u)= [Vu[2, the associated probability density function P({ U - u}) -- e -f(u)

- - e -[Vu[2

is large in regions where u is almost constant, and small in regions of high oscillations. Or, in other words, the image intensity u is considered to be reliable if the gradient is low. In establishing the link between stochastic models and Tikhonov-type regularization we assumed Gaussian white noise (88) and (89). Following the derivation above, minimization principles get much more complicated if we skip the assumption of Gaussian white noise. To highlight the arising complications we consider exemplarily R a y l e i g h distributed noise, that is

P({N-

u ~ - u}) -

lU3 -- U[

0"2

e

-[(uS-u)2/(2a2)]

(91)

Proceeding as above we find that maximization of the conditional probability results in minimization of the functional

/ ( vu2 +

,og u u,)

92,

OTMAR SCHERZER

520

There is no scale space associated with this regularization functional. However, an inverse scale-space method can be constructed, by considering iterative minimization of the functionals

L (fllV(U- u(k-1))12 -.1_(U82o.2--U)2 -- log(lu 8 -

ul)),

k - 1, 2, . . . ,

(93)

and denoting a minimizer by u (k), which then satisfies the optimality condition

f l A ( u (k) --

1

u(k-1)) __ "2~2 _

1

)

21ua - u(k)l 2'

(u(k)

_ uS).

Setting tk = kfl, k = 1 , . . . , and U (k) = u(kAt), we get by taking the limit fl --+ O+ the inverse scale-space method: Ou A -~(t,x)--

(1

2~r2

1 21Ua(X)_u(t,X)12

u(0, x ) = 0 for x e r

)(u(t,x)

u~ --

(X))

for (t,x)e(O, (x)) x g2, (94)

We present filtering of data degraded with Rayleigh noise, with a = 0.2 (cf. Figure 37). Figure 38 shows the solution of (94) at specified time. Finally we compare the stochastic regularization method with well-established diffusion filtering methods: the quality of the stochastic regularization is completely different from diffusion-type filtering. We have selected test data that are extremely distorted by Rayleigh noise. The filtered images in Figure 39 were obtained with a large as possible stopping time (to reduce the effect of noise), so that still a number of details, like the tripod of the cameraman, could be recovered. The selected stopping time is too small for the mean curvature motion filtering and the anisotropic diffusion filtering to completely smear out the noise. We find that the number of preserved details in the filtered images is optimal for the total variation flow and the stochastic method. The good performance of these methods is due to the fact that 9 the stochastic method uses a priori information on the noise and 9 the total variation filtering optimally incorporates information on the image data, which is a blocky image.

DENOISING AND INVERSE PROBLEMS

521

FIGURE 38. Denoising of the image represented in Figure 37 with the inverse scale-space method (94) at time t = 0.04, 0.12, 0.32, 0.56, 1.2, 2.4.

522

OTMAR SCHERZER

FIGURE 39. Denoising. Top left: Heat equation. Top right: Total variation flow. Middle left: Perona-Malik diffusion. Middle right: Anisotropic diffusion. Bottom left: Mean curvature motion. Bottom right: Stochastic regularization.

XI. CONCLUSIONS In this chapter we reviewed interactions between variational methods, diffusion filtering for denoising and image smoothing, and reviewed links to splines, wavelets, and statistical methods. We presented an introduction to inverse problems, such as deblurring and deconvolution, and highlighted numerical methods for their solution. Various other important image processing applications which are solved by variational methods and partial differential equations have not been touched on: 9 O p t i c a l f l o w models [5,18,27,51,75,97,100,104,120,121,149,150,173,174].

DENOISING AND INVERSE PROBLEMS

9 Computational

anatomy

523

and image registration [12,13,65,75,81,96,

162,163]. 9 Inpainting [20,24,44,110,111,141]. 9 Diffusion and regularization o f vector-valued data, such as color images and tensor-valued medical images data [26,136,167]. 9 Blind deconvolution [34,45,47,144]. 9 Level set methods [69-72,93,132,137,155]. 9 Surface smoothing [61,144]. 9 Active contours [161]. 9 Other variational techniques [129,130,159,160].

We did not attempt to give a complete list of references on these topics, since they are not within the main goal of this chapter. We apologize for any reference that has been omitted. For the reader interested in these topics it should be possible to get a complete account from the references listed in these papers.

ACKNOWLEDGMENTS

This work has been supported by the Austrian Fonds zur F6rderung der Wissenschaftlichen Forschung (FWF), grant Y-123 INF. The author thanks H. Grossauer, M. Haltmeier, W. Hinterberger, R. Kowar, J. Ktinstle, S. Leimgruber, and G. Regensburger. Moreover, the author is grateful for the agreement of Ch. Groetsch, M. Hintermtiller, K. Kunisch, M. Oehsen, E. Radmoser, and J. Weickert to use some data of previous joint publications.

REFERENCES 1. Acar, R., and Vogel, C. R. (1994). Analysis of bounded variation penalty methods for illposed problems. Inverse Probl. 10, 1217-1229. 2. Alvarez, L., Guichard, F., Lions, P.-L., and Morel, J.-M. (1993). Axioms and fundamental equations of image processing. Arch. Ration. Mech. Anal. 123, 199-257. 3. Alvarez, L., Lions, P.-L., and Morel, J.-M. (1992). Image selective smoothing and edge detection by nonlinear diffusion. II. S I A M J. Numer. Anal. 29, 845-866. 4. Alvarez, L., and Morel, J.-M. (1994). Formalization and computational aspects of image analysis. Acta Numerica, 1-59. 5. Alvarez, L., Weickert, J., and Sanchez, J. (1999) A scale-space approach to nonlocal optical flow calculations, in [127], pp. 235-246. 6. Ambrosio, L. (1989). A compactness theorem for a new class of functions of bounded variation. Boll. Un. Mat. Ital. B 3, 857-881.

524

OTMAR SCHERZER

7. Ambrosio, L. (1989). Variational problems in SBV and image segmentation. Acta Appl. Math. 17, 1-40. 8. Ambrosio, L. (1990). Existence theory for a new class of variational problems. Arch. Rational Mech. Anal. 111, 291-322. 9. Ambrosio, L., Fusco, N., and Pallara, D. (2000). Functions of Bounded Fariation and Free Discontinuity Problems. New York: Oxford University Press. 10. Ambrosio, L., and Tortorelli, V. M. (1990). Approximation of functionals depending on jumps by elliptic functionals via F-convergence. Commun. Pure Appl. Math. 43, 999-1036. 11. Ambrosio, L., and Tortorelli, V. M. (1992). On the approximation of free discontinuity problems. Boll. Un. Mat. Ital. B (7) 6, 105-123. 12. Amit, Y. (1994). A nonlinear variational problem for image matching. S I A M J. Sci. Comput. 15, 207-224. 13. Amit, Y., Grenander, U., and Piccioni, M. (1991). Structural image restoration through deformable templates. J. Amer. Statist. Assoc. 86, 376-387. 14. Andreu, F., Ballester, C., Caselles, V., and Maz6n, J. M. (2000). Minimizing total variation flow. C.R. Acad. Sci. Paris S~r. I Math. 331, 867-872. 15. Andreu, F., Ballester, C., Caselles, V., and Maz6n, J. M. (2001). The Dirichlet problem for the total variation flow. J. Funct. Anal. 180, 347-403. 16. Andreu, F., Ballester, C., Caselles, V., and Maz6n, J. M. (2001). Minimizing total variation flow. Differential Integral Equations 14, 321-360. 17. Andreu, F., Caselles, V., Diaz, J. I., and Maz6n, J. M. (2002). Some qualitative properties for the total variation flow. J. Funct. Anal. 188, 516-547. 18. Aubert, G., Deriche, R., and Kornprobst, P. (1999). Computing optical flow via variational techniques. S I A M J. Appl. Math. 60, 156-182. 19. Aubert, G., and Kornprobst, P. (2002). Mathematical Problems in Image Processing. New York: Springer-Verlag. 20. Ballester, C., Caselles, V., Verdera, J., Bertalmio, M., and Sapiro, G. (2001). A variational model for filling-in gray level and color images, in Proceedings of the Eighth International Conference On Computer Vision (ICCV-01), Los Alamitos, CA: IEEE Computer Society, pp. 10-16. 21. Banks, H. T., and Kunisch, K. (1989). Estimation Techniques for Distributed Parameter Systems. Basel: Birkh~iuser. 22. Barenblatt, G. I., Bertsch, M., Dal Passo, R., and Ughi, M. (1993). A degenerate pseudoparabolic regularization of a nonlinear forward-backward heat equation arising in the theory of heat and mass exchange in stably stratified turbulent shear flow. S l A M J. Math. Anal. 24, 1414-1439. 23. Bellettini, G., Caselles, V., and Novaga, M. (2001). The total variation flow in R N. Technical report, Sezione de analisi matematica e probabilit~t, Universita di Pisa. 24. Bertalmio, M., Sapiro, G., Caselles, V., and Ballester, C. (2000). Image inpainting, in Proceedings of the Computer Graphics Conference 2000 (SIGGRAPH-O0). New York: ACM Press, pp. 417-424. 25. Bertero, M., and Boccacci, P. (1998). Introduction to Inverse Problems in Imaging. London: lOP Publishing. 26. Black, M., Sapiro, G., Marimont, D., and Heeger, H. (1997). Robust anisotropic diffusion and sharpening of scalar and vector images, in [166], pp. 263-266. 27. Black, M. J., and Anandan, P. (1991). Robust dynamic motion estimation over time, in Proc. IEEE Comp. Soc. Conf. on Computer Vision and Pattern Recognition (CVPR '91). Los Alamitos, CA: IEEE Computer Society Press, pp. 292-302. 28. Bracewell, R. N., and Riddle, A. C. (1967). Inversion of fan-beam scans in radio astronomy. Astrophys. J. 150, 427-434.

DENOISING AND INVERSE PROBLEMS

525

29. Braides, A. (2000). Free discontinuity problems and their non-local approximation, in Calculus of Variations and Partial Differential Equations. Berlin: Springer, pp. 171-180. 30. Braides, A., and Dal Maso, G. (1997). Non-local approximation of the Mumford-Shah functional. Calc. Var. Partial Differential Equations 5, 293-322. 31. Brezis, H. (1973). Operateurs Maximaux Monotones et semi-groupes de contractions dans les espaces de Hilbert. Amsterdam: North-Holland. 32. H. Brezis (1983). Analyse fonctionnelle. Theorie et applications. Collection Mathematiques Appliquees pour la Maitrise. Paris: Masson. 33. Bronstein, I. N., Semendjajew, K. A., Musiol, G., and Muehlig, H. (1997). Taschenbuch der Mathematik (Handbook of Mathematics), 3rd edition. Frankfurt am Main: Deutsch. 34. Burger, M., and Scherzer, O. (2001). Regularization methods for blind deconvolution and blind source separation problems. Math. Control Signals Systems 14, 358-383. 35. Cadzow, J. A. (1979). An extrapolation procedure for band-limited signals. IEEE Trans. Acoust. Speech Signal Process. 27, 4-12. 36. Canny, J. F. (1986). A computational approach to edge detection. IEEE Trans. Pattern Anal. Machine Intell. PAMI-8(6), 679-697. 37. Catt6, F., Lions, P.-L., Morel, J.-M., and Coll, T. (1992). Image selective smoothing and edge detection by nonlinear diffusion. S I A M J. Numer. Anal. 29, 182-193. 38. Chambolle, A. (1999). Finite-differences discretizations of the Mumford-Shah functional. M 2 A N Math. Model. Numer. Anal. 33, 261-288. 39. Chambolle, A., and Dal Maso, G. (1999). Discrete approximation of the Mumford-Shah functional in dimension two. M 2 A N Math. Model. Numer. Anal. 33, 651-672. 40. Chambolle, A., DeVore, R. A., Lee, N., and Lucier, B. J. (1998). Nonlinear wavelet image processing: variational problems, compression, and noise removal through wavelet shrinkage. IEEE Trans. Image Process. 7, 319-335. 41. Chambolle, A., and Lions, P. L. (1997). Image recovery via total variation minimization and related problems. Numer. Math. 76, 167-188. 42. Chambolle, A., and Lucier, B. J. (2001). Interpreting translation invariant wavelet shrinkage as a new image smoothing scale space. IEEE Trans. Image Process. 10, 993-1000. 43. Chan, T. F., Golub, G. H., and Mulet, P. (1999). A nonlinear primal-dual method for total variation-based image restoration. S I A M J. Sci. Comput. 20, 1964-1977 (electronic). 44. Chan, T. F., and Shen, J. (2002). Mathematical models for local nontexture inpaintings. S I A M J. Appl. Math. 62, 1019-1043 (electronic). 45. Chan, T. F., and Wong, C. K. (2000). Convergence of the alternating minimization algorithm for blind deconvolution. Linear Algebra Appl. 316, 259-285. 46. Chan, T. F., Golub, G. H., and Mulet, P. (1996). A nonlinear primal-dual method for total variation-based image restoration, in Proceedings ICAOS "96, Berlin: Springer, pp. 241252. 47. Chan, T. F., and Wong, C. K. (1998). Total variation blind deconvolution. IEEE Trans. Image Process 7, 370-375. 48. Chorin, A. J., and Marsden, J. E. (1993). A Mathematical Introduction to Fluid Mechanics, 3rd edition. New York: Springer-Verlag. 49. Chyzak, F., Paule, P., Scherzer, O., Schoisswohl, A., and Zimmermann, B. (2001). The constrution of orthonormal wavelets using symbolic methods and a matrix analytical approach for wavelets on the interval. Exp. Math. 10, 67-86. 50. Cohen, A., Daubechies, I., and Vial, P. (1993). Wavelets on the interval and fast wavelet transforms. Appl. Comput. Harmon. Anal. 1, 54-81. 51. Cohen, I. (1993). Nonlinear variational method for optical flow computation, in Proe. Eighth Scandinavian Conf. on Image Analysis, Vol. 1, Tromso, pp. 523-530.

526

OTMAR SCHERZER

52. Colton, D., and Kress, R. (1983). Integral Equation Methods in Scattering Theory. New York: Wiley. 53. Colton, D., and Kress, R. (1992). Inverse Acoustic and Electromagnetic Scattering Theory. New York: Springer-Verlag. 54. Dacorogna, B. (1982). Weak Continuity and Weak Lower Semi-Continuity of Non-Linear Functionals. Berlin: Springer-Verlag. 55. Dacorogna, B. (1989). Direct Methods in the Calculus of Variations. Berlin: Springer-Verlag. 56. Daubechies, I. (1992). Ten Lectures on Wavelets. Philadelphia, PA: SIAM. 57. Davies, P. L., and Kovac, A. (2001). Local extremes, runs, strings and multiresolution. Ann. Statist. 29, 1-65 (with discussion and rejoinder by the authors). 58. De Boor, C. (1978). A Practical Guide to Splines. New York: Springer. 59. De Giorgi, E., Carriero, M., and Leaci, A. (1989). Existence theorem for a minimum problem with free discontinuity set. Arch. Ration. Mech. Anal. 108, 195-218. 60. Dibos, F., and $6r6, E. (1997). An approximation result for the minimizers of the Mumford-Shah functional. Boll. Un. Mat. Ital. A (7) 11, 149-162. 61. Diewald, U., Preusser, T., Rumpf, M., and Strzodka, R. (2000). Diffusion models and their accelerated solution in image and surface processing. Acta Math. Univ. Comenian (NS) 70, 15-31. 62. Dobson, D. C., and Vogel, C. R. (1997). Convergence of an iterative method for total variation denoising. S I A M J. Num. Anal. 34, 1779-1791. 63. Donoho, D. L., and Johnstone, I. (1995). Adapting to unknown smoothness via wavelet shrinkage. J. Amer. Statist. Assoc. 90, 1200-1224. 64. Donoho, D. L., Johnstone, I. M., Kerkyacharian, G. B., and Picard, D. (1996). Density estimation by wavelet thresholding. Ann. Statist. 24, 508-539. 65. Dupuis, P., Grenander, U., and Miller, M. I. (1998). Variational problems on flows of diffeomorphisms for image matching. Quart. Appl. Math. 56, 587-600. 66. Engl, H. W., and Groetsch, C. W., eds (1987). Inverse and Ill-Posed Problems. Boston: Academic Press. 67. Engl, H. W., Hanke, M., and Neubauer, A. (1996). Regularization of Inverse Problems. Dordrecht: Kluwer Academic. 68. Engl, H. W., Kunisch, K., and Neubauer, A. (1989). Convergence rates for Tikhonov regularization of nonlinear ill-posed problems. Inverse Probl. 5, 523-540. 69. Evans, L. C., and Spruck, J. (1991). Motion of level sets by mean curvature. I. J. Differ. Geom. 33, 635-681. 70. Evans, L. C., and Spruck, J. (1992). Motion of level sets by mean curvature. II. Trans. Am. Math. Soc. 330, 321-332. 71. Evans, L. C., and Spruck, J. (1992). Motion of level sets by mean curvature. III. J. Geom. Anal. 2, 121-150. 72. Evans, L. C., and Spruck, J. (1995). Motion of level sets by mean curvature. IV. J. Geom. Anal. 5, 79-116. 73. Evans, L. C., and Gariepy, R. F. (1992). Measure Theory and Fine Properties of Functions. Boca Raton: CRC Press. 74. Fasano, A., and Primicerio, M., eds. (1994). Proc. Seventh European Conf. on Mathematics in Industry. Stuttgart: Teubner. 75. Fischer, B., and Modersitzki, J. (1999). Fast inversion of matrices in image processing. Numer. Algor. 22, 1-11. 76. Florack, L. (1997). Image Structure. Dordrecht: Kluwer. 77. Frigaard, I. A., Ngwa, G., and Scherzer, O. (2002). On effective stopping time selection for visco-plastic nonlinear BV diffusion filters. S I A M J. Appl. Math. accepted for publication.

DENOISING AND INVERSE PROBLEMS

527

78. Geman, D., and Yang, C. (1995). Nonlinear image recovery with half-quadratic regularization. IEEE Trans. Image Process. 4, 932-945. 79. Gilbert, P. (1972). Iterative methods for the three-dimensional reconstruction of an object from projections. J. Theor. Biol. 36, 105-117. 80. Glowinski, R. (1984). Numerical Methods for Nonlinear Variational Problems. Berlin: Springer. 81. Grenander, U., and Miller, M. I. (1998). Computational anatomy: an emerging discipline. Quart. Appl. Math. 56(4), 617-694. 82. Groetsch, C. W. (1991). Differentiation of approximately specified functions. Amer. Math. Monthly 98, 847-850. 83. Groetsch, C. W. (1993). Inverse Problems in the Mathematical Sciences. Vieweg Mathematics for Scientists and Engineers. Braunschweig: Friedr. Vieweg. 84. Groetsch, C. W. (1999). Inverse Problems. Washington, DC: Mathematical Association of America. 85. Groetsch, C. W., and Scherzer, O. (2000). Nonstationary iterated Tikhonov-Morozov method and third order differential equations for the evaluation of unbounded operators. Math. Meth. Appl. Sci. 23, 1287-1300. 86. Groetsch, C. W. (1983). Comments on Morozov's discrepancy principle, in Improperly Posed Problems and Their Numerical Treatment, edited by G. H/immerlin and K. H. Hoffmann, Basel: Birkhfiuser, pp. 97-104. 87. Groetsch, C. W. (1984). The Theory of Tikhonov Regularization for Fredholm Equations of the First Kind. Boston: Pitman. 88. Groetsch, C. W., and Scherzer, O. (1993). Optimal order of convergence for stable evaluation of differential operators. Electronic J. Differential Equations, 4, 1-10 (http:// ejde.math.unt.edu). 89. Hamza, A. B., Krim, H., and Unal, G. B. (2002). Unifying probabilistic and variational estimation. IEEE Signal Processing Magazine 19, 37-47. 90. Hanke, M. (2002). Grundlagen der Numerischen Mathematik und des Wissenschaftlichen Rechnens. Stuttgart: Teubner. 91. Hanke, M., and Groetsch, C. W. (1998). Nonstationary iterated Tikhonov regularization. J. Optim. Theory Appl. 98, 37-53. 92. Hanke, M., and Scherzer, O. (2001). Inverse problems light: numerical differentiation. Amer. Math. Monthly 108, 512-521. 93. Harabetian, E., and Osher, S. (1998). Regularization of ill:posed problems via the level set approach. SIAM J. Appl. Math. 58, 1689-1706. 94. Herman, G. T. (1980). Image Reconstruction from Projections." The Fundamentals of Computed Tomography. New York: Academic Press. 95. Hinterberger, W., Hintermfiller, M., Kunisch, K., von Oehsen, M., and Scherzer, O. (2002). Tube methods for BV regularization, J. Math. Imag. Vision, accepted for publication. 96. Hinterberger, W., and Scherzer, O. (2001). Models for image interpolation based on the optical flow. Computing 66, 231-247. 97. Hinterberger, W., Scherzer, O., Schn6rr, Ch., and Weickert, J. (2002). Analysis of optical flow models in the framework of calculus of variations. Num. Funct. Anal. Opt. 23, 69-90. 98. Hoel, P. G. (1960). Elementary Statistics. New York: John Wiley. 99. Hofmann, B. (1999). Mathematik inverser Probleme (Mathematics of Inverse Problems). Stuttgart: Teubner. 100. Horn, B., and Schunck, B. (1981). Determining optical flow. Artif. Intell. 17, 185-203. 101. Kerckhove, M., ed. (2001). Scale-Space and Morphology in Computer Vision, Notes in Computer Science, LNCS 2106. New York: Springer Verlag.

328

OTMAR SCHERZER

102. Kichenassamy, S. (1997). The Perona-Malik paradox. SIAM J. Appl. Math. 57, 13281342. 103. Kirsch, A. (1996). An Introduction to the Mathematical Theory of Inverse Problems. New York: Springer-Verlag. 104. Kornprobst, P., Deriche, R., and Aubert, G. (1999). Image sequence analysis via partial differential equations. J. Math. Imaging Vision 11, 5-26. 105. Kress, R. (1989). Linear Integral Equations. Berlin: Springer-Verlag. 106. Lindeberg, T. (1994). Scale-Space Theory in Computer Vision. Boston: Kluwer. 107. Louis, A. K. (1989). Inverse und Schlecht Gestellte Probleme. Stuttgart: Teubner. 108. Mallat, S. (1999). A Wavelet Tour of Signal Processing, 2nd edition. San Diego, CA: Academic Press. 109. Mammen, E., and van de Geer, S. (1997). Locally adaptive regression splines. Ann. Statist. 25, 387-413. 110. Masnou, S. (2002). Disocclusion: a variational approach using level lines. IEEE Trans. Image Process. 11, 68-76. 111. Masnou, S., and Morel, J. M. (1998). Level lines based disocclusion, in [175], pp. 259-263. 112. The MathWorks, http://www.mathworks.com/. MATLAB. 113. Meyer, Y. (1991). Ondeletts sur l'intervalle. Rev. Mat. Iberoam. 7(2), 115-133. 114. Meyer, Y. (2001). Oscillating Patterns in Image Processing and Nonlinear Evolution Equations, Vol. 22 of University Lecture Series. Providence, RI: American Mathematical Society. 115. Morel, J. M., and Solimini, S. (1995). Variational Methods in Image Segmentation. Boston: Birkh~iuser. 116. Morozov, V. A. (1984). Methods for Solving Incorrectly Posed Problems. New York: Springer Verlag. 117. Morozov, V. A. (1993). Regularization Methods for Ill-Posed Problems. Boca Raton: CRC Press. 118. Mumford, D., and Shah, J. (1989). Optimal approximations by piecewise smooth functions and associated variational problems. Commun. Pure Appl. Math. 42, 577-685. 119. Murio, D. A. (1993). The Mollification Method and the Numerical Solution of Ill-Posed Problems. New York: John Wiley. 120. Nagel, N. H. (1987). On the estimation of optical flow: relations between new approaches and some new results. Artif. Intell. 33, 299-324. 121. Nagel, N. H., and Enkelmann, W. (1986). An investigation of smoothness constraints for the estimation of displacement vector fields from image sequences. IEEE Trans. Pattern Anal. Mach. Intell. 8, 565-593. 122. Nashed, M. Z., ed. (1976). Generalized Inverses and Applications. New York: Academic Press. 123. Natterer, F. (1986). The Mathematics of Computerized Tomography. Stuttgart: Teubner. 124. Natterer, F., and Wiibbeling, F. (2001). Mathematical Methods in Image Reconstruction. Philadelphia, PA: Society for Industrial and Applied Mathematics (SIAM). 125. Neubauer, A. (1989). Tikhonov regularization for non-linear ill-posed problems: optimal convergence rates and finite-dimensional approximation. Inverse Probl. 5, 541-557. 126. Nielsen, M., Florack, L., and Deriche, R. (1997). Regularization, scale-space and edge detection filters. J. Math. Imag. Vision 7, 291-307. 127. Nielsen, M., Johansen, P., Olsen, O. F., and Weickert, J., eds. (1999). Scale-Space Theories in Computer Vision. Lecture Notes in Computer Science, Vol. 1683. Springer Verlag. 128. Nikolova, M. (2000). Local strong homogeneity of a regularized estimator. S l A M J. Appl. Math. 61, 633-658.

DENOISING AND INVERSE PROBLEMS

529

129. Nitzberg, M., Mumford, D., and Shiota, T. (1993). Filtering, Segmentation and Depth, Lecture Notes in Computer Science. New York: Springer. I30. Nitzberg, M., and Shiota, T. (1992). Nonlinear image filtering with edge and corner enhancement. IEEE Trans. Pattern Anal. Mach. Intell. 14, 826-833. 131. Orphanoudakis, S., Trahamnias, P., Crowley, J., and Katveas, N., eds. (1998). Proc. Computer Vision and Mobile Robotics Workshop, CMVR'98, Santorini. 132. Paragios, M., ed. (2001). Variational and Level set methods in Computer Vision. Los Alamitos, CA: IEEE Computer Society. 133. Pazy, A. (1983). Semigroups of Linear Operators and Applications to Partial Differential Equations. New York: Springer-Verlag. 134. Perona, P., and Malik, J. (1987). Scale space and edge detection using anisotropic diffusion, in Workshop on Computer Vision. Washington, DC: IEEE Computer Society Press, pp. 16-22. 135. Perona, P., and Malik, J. (1990). Scale space and edge detection using anisotropic diffusion. IEEE Trans. Pattern Anal. Mach. Intell. 12, 629-639. 136. Pollak, I., Krim, H., and Willsky, A. (1998). Stabilized inverse diffusion equations and segmentation of vector-valued images, in [175], pp. 246-248. 137. Preusser, T., and Rumpf, M. (2002). A level set method for anisotropic geometric diffusion in 3D image processing. S I A M J. Appl. Math. 62, 1772-1793 (electronic). 138. Radmoser, E., Scherzer, O., and Weickert, J. (1999). Scale-space properties of regularization methods, in [127], pp. 211-222. 139. Radmoser, E., Scherzer, O., and Weickert, J. (2000). Scale-space properties of nonstationary iterative regularization methods. J. Visual Commun. Image Representation 11, 96-114. 140. Radon, J. (1917). lJber die Bestimmung von Funktionen durch ihre Integralwerte l~ngs gewisser Mannigfaltigkeiten. Ber. Verh. Sachs. Akad. Wiss. Leipzig Math. Phys. Kl. 69. 141. Ramasubramanian, M., Pattanaik, S. N., and Greenberg, D. P. (1999). A perceptually based physical error metric for realistic image synthesis, in Alyn Rockwood, editor, SIGGRAPH 99 Conference Proceedings, Annual Conference Series, Addison Wesley, pp. 73-82. 142. Reinsch, Ch. (1967). Smoothing by spline functions. Numer. Math. 10, 177-183. 143. Rice, J., and Rosenblatt, M. (1983). Smoothing splines: regression, derivatives and deconvolution. Ann. Statist. 11, 141-156. 144. Sapiro, G. (2001). Geometric Partial Differential Equations and Image Analysis. Cambridge: Cambridge University Press. 145. Scherzer, O. (1997). Stable evaluation of differential operators and linear and nonlinear multi-scale filtering. Electronic J. Differential Equations 15, 1-12 (http://ejde.math.unt.edu). 146. Scherzer, O. (2002). Explicit versus implicite relative error regularization on the space of functions of bounded variation, in "Inverse Problems, Image Analysis, and Medical Imaging,"Contemp. Math. 313, 177-198. Providence, RI: American Mathematics Society. 147. Scherzer, O., and Groetsch, C. W. (2001). Inverse scale space theory for inverse problems. In [101], pp. 317-325. 148. Scherzer, O., and Weickert, J. (2000). Relations between regularization and diffusion filtering. J. Math. Imag. Vision 12, 43-63. 149. Schn6rr, C. (1994). Segmentation of visual motion by minimizing convex non-quadratic functionals, in Proc. 12th Int. Conf. on Pattern Recognition, Vol. A, pp. 661-663. 150. Schn6rr, Ch. (1991). Funktionalanalytische Methoden zur Gewinnung yon Bewegungsinformation aus TV-Bildfolgen. PhD thesis, Fakult~it ffir Informatik, University of Karlsruhe. 151. Schoenberg, I. J. (1964). Spline functions and the problem of graduation. Proc. Natl. Acad. Sci. USA 52, 947-950.

530

OTMAR SCHERZER

152. Schultz, M. H. (1973). Spline Analysis. Englewood Cliffs, NJ: Prentice-Hall. 153. Schumaker, L. L. (1981). Spline Functions. Basic Theory. New York: Wiley. 154. Seidman, T. I., and Vogel, C. R. (1989). Well posedness and convergence of some regularisation methods for non-linear ill posed problems. Inverse Probl. 5, 227-238. 155. Sethian, J. A. (1999). Level Set Methods and Fast Marching Methods, 2nd edition. Cambridge: Cambridge University Press. 156. Sporring, J., Nielsen, M., Florack, L., and Johansen, P., eds. (1997). Gaussian Scale-Space Theory. Dordrecht: Kluwer. 157. Strang, G., and Fix, G. J. (1973). An Analysis of the Finite Element Method. Englewood Cliffs, NJ: Prentice-Hall. 158. Strong, D., and Chan, T. F. (1996). Exact solutions to the total variation regularization problem. Technical report, University of California, Los Angeles, CAM 96-41. 159. Terzopoulos, D. (1983). Multilevel computational processes for visual surface reconstructions. Computer Vision, Graphics and Image Processing 24, 52-96. 160. Terzopoulos, D. (1988). The computation of visible-surface representations. IEEE Trans. Pattern Anal. Mach. Intell. 10, 417-438. 161. Terzopoulos, D., Witkin, A., and Kass, M. (1988). Constraints on deformable models: recovering 3D shape and nonrigid motion. Artif. Intell. 36, 91-123. 162. Thirion, J.-P. (1995). Fast non-rigid matching of 3D medical images. Technical Report RR-2547, Inria, Institut National de Recherche en Informatique et en Automatique. 163. Thirion, J.-P. (1996). Non-rigid matching using demons. Preprint, INRIA, France. 164. Tikhonov, A. N., and Arsenin, V. Y. (1977). Solutions of Ill-Posed Problems. Washington, DC: Wiley (translation editor: Fritz John). 165. Torre, V., and Poggio, T. A. (1986). On edge detection. IEEE Trans. Pattern Anal. Mach. Intell. 8, 48-163. 166. Torwick, I., ed. (1997). Proceedings of the 1997 IEEE International Conference on Image Processing (ICIP-97). Los Alamitos, CA: IEEE Computer Society. 167. Tschumperl6, D., and Deriche, R. (2002). Diffusion pdes on vector-valued images. IEEE Signal Processing Magazine 19, 16-25. 168. Unser, M. (1999). Splines--a perfect fit for signal and image processing. IEEE Signal Processing Magazine 16, 22-38. 169. Wahba, G. (1990). Spline Models for Observational Data, Vol. 59. Philadelphia, PA: Society for Industrial and Applied Mathematics (SIAM). 170. Webb, S., ed. (1988). The Physics of Medical Imaging. Bristol: Institute of Physics Publishing. 171. Weickert, J. (1994). Anisotropic diffusion filters for image processing based quality control, in [74]. 172. Weickert, J. (1998). Anistropic Diffusion in Image Processing. Stuttgart: Teubner. 173. Weickert, J. (1999). On discontinuity-preserving optic flow, in [131], pp. 115-252. 174. Weickert, J., and Schn6rr, Ch. (.'2,001). Variational optic flow computation with a spatiotemporal smoothness constraint. J. Math. Imag. Vision, 14, 245-255. 175. Werner, B., ed. (1998). Proceedings of the 1998 IEEE International Conference on Image Processing (ICIP-98). Los Alamitos, CA: IEEE Computer Society.