Anisotropic Diffusion Partial Differential Equations for Multichannel Image Regularization: Framework and Applications

ADVANCES IN IMAGING AND ELECTRON PHYSICS, VOL. 145

Anisotropic Diffusion Partial Differential Equations for Multichannel Image Regularization: Framework and Applications DAVID TSCHUMPERLÉ1 AND RACHID DERICHE2 1 Image Team, GREYC/ENSICAEN–UMR CNRS 6072, 14050 Caen Cedex, France 2 Odyssée Project Team, INRIA/ENPC/ENS–INRIA, 06902 Sophia Antipolis, France

Preliminary Notations . . . . . . . . . . . I. Introduction . . . . . . . . . . . . . . A. Defining a Local Geometry for Multichannel Images . 1. Local Geometric Features . . . . . . . . . 2. Geometry From a Scalar Feature . . . . . . . 3. Di Zenzo Multivalued Geometry . . . . . . . II. PDE-Based Smoothing of Multivalued Images: A Review . A. Variational Methods . . . . . . . . . . . B. Divergence-Based Diffusion PDEs . . . . . . . C. Oriented Heat Flows . . . . . . . . . . . D. Trace-Based Diffusion PDEs . . . . . . . . E. Links Between Existing Regularization Methods . . . III. Curvature-Preserving PDEs . . . . . . . . . . A. The Single-Direction Case . . . . . . . . . B. Curvature-Preserving PDEs and Line Integral Convolutions C. Between Traces and Divergences . . . . . . . D. Extension to Multidirectional Smoothing . . . . . IV. Implementation Considerations . . . . . . . . . V. Applications . . . . . . . . . . . . . . A. Color Image Denoising and Artifact Removal . . . B. Color Image Inpainting . . . . . . . . . . C. Color Image Interpolation . . . . . . . . . D. Flow Visualization . . . . . . . . . . . VI. Conclusion . . . . . . . . . . . . . . . Appendix A . . . . . . . . . . . . . . Appendix B . . . . . . . . . . . . . . Appendix C . . . . . . . . . . . . . . References . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . . . . . .

150 151 154 154 155 156 160 160 163 165 167 172 174 174 177 178 180 181 183 184 185 189 191 193 195 197 198 203

149 ISSN 1076-5670 DOI: 10.1016/S1076-5670(06)45004-7

Copyright 2007, Elsevier Inc. All rights reserved.

150

TSCHUMPERLÉ AND DERICHE

P RELIMINARY N OTATIONS Throughout this chapter, we represent a multichannel or multivalued image by a continuous function I :  → Rn , where  ⊂ R2 is the definition domain of the image (basically a two-dimensional (2D) rectangle W × H ) and n ∈ N+ is the dimension of each vector-valued image pixel I(X) located at X = ( x y )T ∈ . The notation Ii stands for the ith channel of the image I. Note that Ii can be considered itself as a scalar-valued image Ii :  → R. Thus, we have ∀X = (x, y) ∈ ,

I(X) = ( I1(X)

I2(X)

· · · In(X) )T .

For the common case of color images, we naturally get n = 3, that is, three vector components (R, G, B) per pixel, retrieved respectively from the red (I1 ), green (I2 ), and blue (I3 ) channels of a color image I. We also intensely use second-order diffusion tensors in equations. A diffusion tensor D is assimilated to a 2 × 2 symmetric and positive-definite matrix, having then two positive eigenvalues λ1 , λ2 and two associated orthonormal eigenvectors u1 ⊥u2 . The shape of a tensor D may be seen as an ellipse, oriented by the vector basis u1 ⊥u2 and elongated by λ1 and λ2 , as illustrated below.  D=

a b

b c

 = λ1 u1 uT1 + λ2 u2 uT2

When λ2  λ1 (lengthened ellipse), the tensor D is said to be anisotropic and has u2 as its principal orientation. When λ1 = λ2 = β, the tensor D is isotropic and thus equal to a weighted version of the 2 × 2 identity matrix Id   β 0 λ1 = λ2 = β ⇒ D = βId = . 0 β An isotropic tensor D has no privileged orientations, all vectors of R2 being possible eigenvectors of D. Finally, we denote Gσ , a normalized 2D Gaussian function with a standard deviation of σ :  2  1 x + y2 exp − . Gσ (x, y) = 2π σ 2 2σ 2

ANISOTROPIC DIFFUSION PARTIAL DIFFERENTIAL EQUATIONS

151

I. I NTRODUCTION Obtaining regularized versions of noisy or corrupted image data has always been a desirable goal in the fields of computer vision and image processing. Removing noise or scratches from degraded images is indeed a fundamental preprocessing step that can possibly ease the further analysis of the image data by higher-level algorithms such as detectors of important image features (edges, corners, objects, motion). The ability to create simplified versions of the image data also is very interesting, when considering the analysis of the images at multiple scales. In a more general manner, image regularization is one of the key stages of most computer vision algorithms since it plays a fundamental role for solving ill-posed computer vision problems (Hadamard, 1923), including restoration, segmentation, registration, surface reconstruction, and so on. This explains why many image regularization formalisms have already been proposed and studied in the literature. Perona and Malik (1990) in their pioneering work in the early 1990s were the first to imagine image regularization in terms of anisotropic diffusion partial differential equations (PDEs). Their method, applied on scalar-valued images (one value by pixel), has particularly raised a strong interest for PDE-based formulations, since it succeeded in smoothing image data in a nonlinear way, removing the noise quite well while allowing the preservation of significant image features, such as contours and corners (discontinuities of the signal), despite an initial formulation that later was proved to be unstable (Kichenassamy, 1997; Weickert and Benhamouda, 1997a, 1997b). First created to describe physical laws and natural motions of mechanic objects and fluids (strings, water, wind (Wesseling, 2000)), diffusion PDEs had already been widely studied, and interesting theoretical results from the fields of physics and mathematics have found interesting implications for the purpose of data regularization. Actually, PDEs are local formulations and thus, they are well adapted to deal with degraded images where sources of data corruption are local or semilocal. This is not restrictive: Gaussian noise, scratches, or compression artifacts are, for instance, local degradations usually encountered in digital (original or digitized) images. Following the way opened by Perona and Malik (1990), many authors have proposed variants of diffusion PDEs for image regularization since then, primary for the restoration of scalar-valued datasets. Important theoretical contributions in this field concern the way the classical isotropic diffusion equation (heat flow) has been extended to deal with anisotropic smoothing (Black et al., 1998; Krissian, 2000; Perona and Malik, 1990; Sapiro, 2001; ter Haar Romeny, 1994; Weickert, 1998; Yezzi, 1998), how diffusion PDEs may be seen as gradient descents of various energy functionals (Alvarez and Mazorra, 1994; Aubert and Kornprobst, 2002; Blanc-Feraud et al.,

152

TSCHUMPERLÉ AND DERICHE

1995; Charbonnier et al., 1994; Chambolle and Lions, 1997; Charbonnier et al., 1997; Gilboa et al., 2002; Kimmel et al., 2000; Rudin et al., 1992; Teboul et al., 1998), and the link between regularization PDEs and the concept of nonlinear scale spaces (Alvarez et al., 1993; Lindeberg, 1994; Nielsen et al., 1997). Extensions of these techniques to deal with color images, and more generally multichannel datasets, have been more recently attempted by (Blomgren and Chan, 1998; Chan et al., 2000; Kimmel et al., 2000; Pardo and Sapiro, 2000; Sapiro, 2001; Sapiro and Ringach, 1996; Tschumperlé, 2002; Tschumperlé and Deriche, 2003; Tschumperlé and Deriche, 2005; Weickert, 1998; Weickert, 1999; Weickert and Brox, 2002) (among others), leading to more elaborated expressions: A coupling term between image channels generally appears in the equations. Diffusion equations dealing with constrained multidimensional datasets have been also proposed, allowing regularization of images of unit vectors (Coulon et al., 2001; Kimmel and Sochen, 2002; Perona, 1998; Tang et al., 1998; Vese and Osher, 2001), orthonormal matrices (Chefd’hotel et al., 2004; Chu, 1990a, 1990b; Tschumperlé and Deriche, 2001c, 2002b), positive-definite matrices (Chefd’hotel et al., 2004; Tschumperlé and Deriche, 2001b), or image data defined on implicit surfaces (Bertalmio et al., 2001; Chan and Shen, 2000a; Tang et al., 2000). Usually these types of constrained PDEs simply add an extra constraint term to the corresponding unconstrained equation and are not discussed here. Despite the wide range of existing constrained and unconstrained PDE formalisms for scalar and multichannel images, all proposed methods have something in common—a nonlinear regularization PDE such as ∂I ∂t = R locally smooths the image I along one or several directions of the plane that are different at each image point, depending on the local image configuration. Typically, the principal smoothing direction is always chosen to be parallel to the image contours, resulting in an anisotropic regularization that does not destroy the edges. This has an interesting interpretation in terms of scalespace: As the image data are gently regularized step by step, a continuous sequence of smoother images I(t) is generated while the evolution time t of the PDE goes by. Obviously, anisotropic regularization algorithms must let the less significant data features disappear first (preferably noise), while the interesting image details (edges) are preserved as long as they become unimportant themselves within the image (Alvarez et al., 1993; Lindeberg, 1994; Nielsen et al., 1997; Perona and Malik, 1990; Witkin, 1983). Roughly speaking, regularization PDEs may be seen as iterative and nonlinear filters that simplify the image little by little and then minimize the image variations (Figure 1). Note that such equations generally do not converge toward a very interesting solution. Basically, the image obtained at convergence (t → ∞)

ANISOTROPIC DIFFUSION PARTIAL DIFFERENTIAL EQUATIONS

(a)

(b)

(c)

(d)

153

(e)

F IGURE 1. Nonlinear regularization PDEs and the notion of anisotropic scale-space. (a) Initial image I(t=0) , (b) t = 50, (c) t = 250, (d) t = 1000, (e) t = 3000.

is constant everywhere, corresponding to an image without any variations. This is indeed the most simplified image obtainable. To avoid this undesired oversimplification, regularization algorithms are usually based on a modified PDE velocity R = R + α(Inoisy − I ), including a so-called data fidelity term weighted by a user-defined parameter α ∈ R+ . It prevents the expected solution (regularized image) at convergence from being too different from the original noisy image (not constant, by the way). Another classical restoration technique is done by stopping the pure regularization flow ∂I ∂t = R after a finite number of iterations (which thus becomes a parameter of the method). Here, the main interest is in the regularization term R itself rather than the one containing the fidelity term R . For a broader mathematical study of linear or nonlinear fidelity terms, please refer to (Meyer, 2001; Nikolova, 2001; Nikolova and Ng, 2001). As it is clear that local and oriented image smoothing is one of the key ideas used by most PDE-based regularization methods, this leads to the problem of defining a coherent geometry from a multichannel image, which must be the first aim of a good regularization algorithm. Following this simple and general principle, recent contributions (Tschumperlé and Deriche, 2003; Tschumperlé and Deriche, 2005; Weickert, 1998) proposed two different and generic PDE-based frameworks able to design regularization processes from any given underlying local smoothing geometry. These methods have two main interests: On the one hand, they unify many previously proposed equations into generic diffusion PDEs and provide a local geometric interpretation of the corresponding regularizations. On the other hand, they clearly separate the design of the smoothing geometry from the smoothing process itself. In a first step, the geometry of the structures inside the image is retrieved (generally by the computation of the so-called structure tensor field). Then a local geometry of the desired smoothing is defined by the mean of a second field of diffusion tensors, depending on the first one. Finally, one step of the smoothing process itself is performed through one or several iterations of a specific diffusion PDE. This procedure is repeated until the image is sufficiently regularized.

154

TSCHUMPERLÉ AND DERICHE

This chapter first discusses the definition of a local geometry for multichannel images, by reviewing and comparing proposed solutions in the literature (Blomgren, 1998; Blomgren and Chan, 1998; Di Zenzo, 1986; Sapiro, 1996; Weickert, 1998) (Section I.A). This is followed by a review of important works already proposed for scalar and multichannel image regularization within a diffusion PDE framework. These methods may be classified into three different approaches: (1) variational formulations, (2) divergence expressions, and (3) oriented Laplacians. The main focus is on the interpretation of the algorithms in terms of local smoothing (Section II). We particularly note the advantages and drawbacks of each equation in real cases. Then we focus on a very recent alternative, formulated as a tensor-driven diffusion that regularizes multichannel images while taking specific curvature constraints into account (Section III). This formulation is mathematically positioned between previous existing equations in such a way that it solves most issues encountered with classical regularization methods. Moreover, we show that a theoretical interpretation of the curvature-constrained formalism exists in terms of line integral convolutions (LICs), which is a simple filtering technique originally proposed by Cabral and Leedom (1993). This direct analogy allows the design of an explicit numerical scheme that implements the regularization PDE by successive integrations of pixel values along integral lines (Section IV). This iterative scheme has two main advantages compared with classical PDE implementations. On one hand, it preserves thin image structures remarkably well, since it naturally works at a subpixel accuracy, thanks to the use of a fourth-order Runge–Kutta integration. On the other hand, the algorithm is able to run up to three times faster than classical explicit schemes since it is unconditionally stable, even for large PDE time steps. The described method makes diffusion PDEs a generic and very efficient approach to solving image-processing problems needing multichannel image regularization. Finally, we illustrate this effectiveness, in terms of computational speed and visual quality, with results on color image restoration, color image inpainting, and nonlinear resizing, among all possible applications in the area of image regularization (Section V). A. Defining a Local Geometry for Multichannel Images 1. Local Geometric Features As stated in the introduction, image regularization may be considered a filter that reduces local pixel variations. More precisely, one wants to smooth a multichannel image I :  → Rn while preserving its edges (discontinuities in the image intensities), that is, performing a local smoothing mostly along

ANISOTROPIC DIFFUSION PARTIAL DIFFERENTIAL EQUATIONS

155

directions of the edges, avoiding a smoothing orthogonal to these edges. At first glance, a naive approach would be to apply a scalar-valued regularization filter on each channel Ii of the multichannel image I, doing so independently for each i = 1, . . . , n. However, the correlation between image channels would be ignored in this case, and it might cause important disparities in the smoothing behavior, because local smoothing directions and amplitudes could be very different from each channel to another. Such decoupled regularization methods generally lead to undesirable oversmoothing effects, destroying significant edge structures in the image. Multichannel image regularization is based instead on a coherent image smoothing that locally uses the same smoothing directions and amplitudes for all image channels Ii . Naturally, this means that the local geometry of a multichannel image I must be measured first. Such a geometry consists in the definition of these important features at each image point X = (x, y) ∈  of I: + − • Two orthogonal directions θ(X) , θ(X) ∈ S1 (unit vectors of R2 ) directed respectively across and along the edges (generally the maximum and minimum variations of the image intensities at X). The direction θ − generally corresponds to the edge direction, when there is one, while θ + naturally extends the notion of gradient direction for multichannel images. • A corresponding variation norm N (X) measuring the local strength of an edge. This is the extension of the vector gradient norm for multichannel images.

In order to construct such a vector geometry, different approaches have been considered and are detailed in the following text. 2. Geometry From a Scalar Feature One simple method consists of first computing a scalar image f (I), using a vector to scalar function f : Rn → R that would ideally model the human perception of vector-valued edges. It is particularly conceivable for color images: One may choose, for instance, the lightness function (perceptual response to the luminance) coming from the CIELAB color base (Poynton, 1995): f = L∗ = 116g(Y ) − 16

with Y = 0.2125R + 0.7154G + 0.0721B

where g : R → R is defined by  √ g(s) = 3 s g(s) = 7.787s +

16 116

if s > 0.008856, else.

156

TSCHUMPERLÉ AND DERICHE

(a)

(b)

(c)

(d)

(e)

F IGURE 2. Using lightness L∗ to detect geometry of a color image fails for isolightness contours. (a) Red channel R, (b) green channel G, (c) blue channel B, (d) color image (R, G, B), (e) lightness (scalar) image L∗ .

Thus, we may define a vector-valued local vector geometry {N , θ+ , θ− } of I by choosing  ∇f (I)   θ+ = ∇f (I) , and N = ∇f (I). θ− ⊥θ+ , However, this method has two major drawbacks. First, it is not always possible to easily define a significant function f for multichannel images (particularly when the number of channel is n > 3). Second, there are mathematically no functions f that can detect all possible vector-valued variations. For instance, the lightness function defined above cannot detect isolightness vector contours in a color image. This is the case for the image shown on Figure 2: The contours inside the colored yin-yang symbol will not be detected by N = ∇f (I), since f (I) is constant therein. As a consequence, the smoothing performed will be either isotropic or oriented in an incorrect direction: the existing color edges inside the yin-yang symbol will probably be blurred. 3. Di Zenzo Multivalued Geometry A very elegant solution has been proposed by Di Zenzo (1986) to overcome this limitation. He considers a multichannel image I :  → Rn as a vector field, and looks for the local variations of the vector norm dI2 , mainly given by a variation matrix G = (gi,j ). This yields dI = Ix dx + Iy dy,

where Ix =

∂I ∂I and Iy = ∂x ∂y

(∈ Rn ),

then dI2 = dIT dI = Ix 2 dx 2 + 2ITx Iy dx dy + Iy 2 dy 2 , that is, dI = dX G dX, 2

T

where G =

n  i=1

 ∇Ii ∇IiT

and dX =

dx dy

 ,

ANISOTROPIC DIFFUSION PARTIAL DIFFERENTIAL EQUATIONS

157

G is denoted as the structure tensor. It sums variation contributions from each image channel Ii . It is easy to see that G is a 2 × 2 symmetric and semipositive-definite matrix. Its coefficients (gi,j ) are as follows: ⎧  ⎪ g11 = ni=1 Ii2x , ⎪ ⎨  g12 = g21 = ni=1 Iix Iiy , ⎪ ⎪ ⎩ g22 = n I 2 . i=1 iy

In the common case of color images I = (R, G, B), G is defined as   Rx2 + G2x + Bx2 Rx Ry + Gx Gy + Bx By G= . Rx Ry + Gx Gy + Bx By Ry2 + G2y + By2

(1)

The interesting point about G is that its positive eigenvalues λ+/− give the maximum and the minimum values of dI2 , while the orthogonal eigenvectors θ+ and θ− are the corresponding orientations of these extrema, and are formally given by √   g11 + g22 ±

2g12 √ and θ+/−  λ+/− = , (2) g22 − g11 ±

2 2 . The vectors θ are normalized to the unit where = (g11 − g22 )2 + 4g12 ± vector afterward. With this simple and efficient approach, Di Zenzo opened a natural way to deal with the local vector geometry of multichannel images, through the use of the oriented orthogonal basis (θ+ , θ− ) and the variations measures (λ+ , λ− ). A slight variant has been proposed by Weickert (1998). He proposed instead to study the eigenvalues and eigenvectors of a Gaussian-smoothed version Gσ of the structure tensor G:

Gσ =

n 

  ∇Iiα ∇IiTα ∗ Gσ ,

where ∇Iiα = ∇(Ii ∗ Gα ),

(3)

i=1

where Gα and Gσ are 2D Gaussian kernels with variances respectively equal to α and σ . User-defined parameters α and σ have an influence on the smoothness of the obtained structure tensor field, and by extension, on the regularity of the retrieved vector-valued image geometry. It is noteworthy that eigenvalues of Gσ are well adapted to discriminate different geometric cases: • When λ+  λ−  0, there are very few vector variations around the current point X = (x, y); the region is almost flat and does not contain any edges or corners (as is the case for the inside of the strips in Figure 3a). For this configuration, the variation norm N that must be defined should be low.

158

TSCHUMPERLÉ AND DERICHE

(a)

(b)

(c)

(d)

F IGURE 3. Comparing possible vector variation norms N , N− , and N+ for a synthetic √ √ λ+ , (c) N− = λ+ − λ− , color image. (a) Color checkerboard (size 40 × 40), (b) N = √ (d) N+ = λ+ + λ− .

• When λ+  λ− , there are many vector variations. The current point may be located on a vector edge (as for the edges of the strips in Figure 3a). For this configuration, the variation norm N should be high. • When λ+  λ−  0, the variation are located on a saddle point of the vector surface, which could be a corner structure in the image (for instance, the intersections of the strips in Figure 3a). In this case, N should be even higher than for the previous configuration. Regularization algorithms indeed have a tendency to smooth corners quickly. A very high variation measure estimated on corner points would attenuate the smoothing there, which is often a desired effect. Many proposed regularization algorithms acting on multichannel images have implicitly or explicitly based their smoothing behavior from these Di Zenzo attributes. In particular, three different choices of vector gradient norms N have been proposed so far in the literature to measure vector-valued variations: √ • N = λ+ , as a natural extension of the scalar gradient norm viewed as the value of maximum variations (Blomgren, 1998; Sapiro, 1996; Sapiro, 1997) (Figures 3b and 4b). This norm will not give particular importance to corners√compared with straight edges. λ+ − λ− , also called the coherence norm, has been chosen • N− = by (Sapiro and Ringach, 1996; Weickert, 1996a; Weickert, 1997a). Note that this norm fails to detect discontinuities that are saddle points of the vector-valued surface. This is illustrated on the intersections of the strips (Figure 3c), as well as in the center and the left and right parts of the child’s eye (Figure 4c). This will mainly perturb any regularization process that uses this norm since some colored sharp corners, considered as homogeneous regions, will be probably oversmoothed. √ • N+ = λ+ + λ− , is also denoted by ∇I, is often chosen (Bertalmio et al., 2000a; Blomgren and Chan, 1998; Pardo and Sapiro, 2000; Tang et al.,

ANISOTROPIC DIFFUSION PARTIAL DIFFERENTIAL EQUATIONS

(a)

(b)

(c)

159

(d)

F IGURE 4. Comparing possible vector variation norms N , N− , and N+ for a real color √ √ image. (a) Color photograph (small portion 60 × 40), (b) N = λ+ , (c) N− = λ+ − λ− , √ (d) N+ = λ+ + λ− .

1998; Tschumperlé and Deriche, 2001b; Tschumperlé and Deriche, 2001c) because it detects edges and corners well and is easy to compute. It does not require an eigenvalue decomposition of G as the other norms did, because   n   N+ = ∇I = trace(G) =  ∇Ii 2 . (4) i=1

Moreover, the norm N+ has the interesting property of giving preferences to certain corners (Figure 3d). This is very valuable for image restoration purposes, since the smoothing can be attenuated on high-curvature structures that are classically difficult to preserve. Note that for the scalar case (n = 1), the structure tensor calculus reduces to: dI 2 = dX G1 dX,   2 Ix Iy Ix where G1 = ∇I ∇I T = . Ix Iy Iy2

when n = 1,

1 In this case, the eigenvectors θ+/− and the eigenvalues λ1+/− of G1 are   1 ∇I ⊥ θ−1 = ξ = ∇I λ− = 0, , associated to ∇I λ1+ = ∇I 2 . θ+1 = η = ∇I , 

Basically, the three above-defined norms N+ , N− , and N all reduce to ∇I  in the case of scalar-valued images, which is a desired property. Once a local vector geometry is defined, it can be used as a measure in many image analysis processes involving multichannel images (not only for regularization algorithms). For instance, color edge detection may be performed by finding threshold local maxima of the N+ norm (Figure 5 and Koschan, 1995;

160

TSCHUMPERLÉ AND DERICHE

(a)

(b)

F IGURE 5. Using a vector variation norm for color edge detection. (a) Color image, (b) detect√ ing color edges with the norm N+ = λ+ + λ− .

Tschumperlé and Deriche, 2001a; Tschumperlé and Deriche, 2002a). This vector geometry computation has also been integrated as a measure of contours in some multichannel image segmentation methods (Sapiro, 1996; Sapiro, 1997). √ For all the reasons listed above, the norm N+ = λ+ + λ− associated with the Di Zenzo geometry is probably one of the best measures for detecting local variations in multichannel images and is considered in the next sections of this chapter.

II. PDE-BASED S MOOTHING OF M ULTIVALUED I MAGES : A R EVIEW We review and propose a classification of classical smoothing methods based on diffusion PDEs into three different approaches, related to different interpretation levels of the regularization processes, from the most global to the most local ones. Each section starts with a description of the original idea for scalar-valued images and is then extended for multichannel datasets. A. Variational Methods Contrary to the formulation of the original Perona–Malik equation, several methods have been proposed to resolve the problem of image regularization as a global minimization procedure, within a variational framework. Formalisms described by (Aubert and Kornprobst, 2002; Chambolle and Lions, 1997; Charbonnier et al., 1997; Sapiro, 2001; Weickert, 1998), among numerous references, contributed to the definition of generic energy functionals measuring global image variations. The goal is to minimize adapted variation

ANISOTROPIC DIFFUSION PARTIAL DIFFERENTIAL EQUATIONS

Function name

φ(s)

Reference

Tikhonov Perona–Malik Minimal surfaces Geman–McClure Total variation Green

s2

Tikhonov, 1963 Perona and Malik, 1990 Charbonnier et al., 1994 Geman and McClure, 1985 Rudin et al., 1992 Green, 1990

F IGURE 6.

1 − exp(−s 2 /K 2 ) 2 1 + s2 − 2 s 2 /(1 + s 2 ) s 2 log(cosh(s))

161

List of different φ-functions and corresponding references.

functionals to flatten low image variations (then gradually remove the noise), while preserving the high ones (avoiding the smoothing of image contours). The formulation of the φ-functionals gathers some of these approaches in a general framework and yields a very unifying way to proceed. A noisy scalar image Inoisy can be regularized by minimizing the following φ-functional: 

 (5) min E(I ) = φ ∇I  d, I : →R



where φ : R → R is an increasing function, directing the regularization behavior and penalizing high gradient norms. The minimization is performed via the corresponding diffusion PDE evolution, coming from the Euler– Lagrange equations of E(I ):  I(t=0) = Inoisy ,

φ  (∇I )  (6) ∂I ∂t = div ∇I  ∇I . Different choices of functions φ lead to different proposed regularization methods, especially the simple isotropic smoothing (equivalent to a Gaussian convolution), as introduced by Tikhonov (Tikhonov, 1963), as well as the well-known Perona–Malik (Perona and Malik, 1990) and total variation (TV) anisotropic flows (Rudin et al., 1992). Many regularization methods acting on scalar-valued images have been unified by the φ-function formalism (Figure 6). A natural extension of the φ-functionals for the regularization of multichannel images I could consist in minimizing the following cost functional E(I) measuring a global multichannel image variation: 

 (7) min n E(I) = φ N (I) d, I : →R



where N (I) is one of the three local variation norms defined in Section I.A.

162

TSCHUMPERLÉ AND DERICHE

But more generally, as vector-valued images possess two distinct variation measures λ+ and λ− (eigenvalues of the structure tensor G) contrary to a single measure ∇I  for scalar images, it seems natural to minimize a functional defined by a function ψ : R2 → R of two variables instead of a single one. The ψ-functional below is thus a more complete extension of the φ-function formulation for multichannel images.  min n E(I) = ψ(λ+ , λ− ) d. (8) I : →R



The Euler–Lagrange equations of Eq. (8) can be derived and reduce to a simple form of divergence-based PDE (see Appendix A for details about this Euler–Lagrange derivation):    ∂Ii ∂ψ ∂ψ T T = div (i = 1, . . . , n). (9) θ+ θ + + θ− θ− ∇Ii ∂t ∂λ+ ∂λ− The choice of specific cases of ψ-functions leads to previous vector-valued regularization approaches defined as variational methods, such as the whole range of vector-valued φ-functionals (Blomgren and Chan, 1998; Pardo and Sapiro, 2000; Tang et al., 2000):  ψ(λ+ , λ− ) = φ( λ+ + λ− ) or the Beltrami flow framework (Kimmel et al., 1998; Kimmel et al., 2000; Kimmel and Sochen, 1999; Sochen et al., 2001; Sochen et al., 1998; Sochen, 2001):  ψ(λ+ , λ− ) = (1 + λ+ )(1 + λ− ). Note that this last approach is also equivalent to defining the minimizing functional E(I) as a Polyakov action, which is a physical measure of the area of the image I seen as a 2D surface embedded in a (n + 2)D space. This geometric interpretation helps in understanding how functional minimization can play a role in smoothing images by forcing them to be more regular, thereby finding a minimal surface (Figure 7). Despite the interesting global geometric interpretation of variational formulations, such methods clearly lack flexibility. Indeed, they are formulated as global minimization processes, despite the local geometric smoothing properties that are intrinsically desired for regularization purposes. Such PDEs are obtained by the Euler–Lagrange derivation of a functional and cannot thus be finely tuned to adapt themselves to local geometric cases (e.g., contours, corners). Unfortunately, this adaptability is primordial in many situations, especially when the noise level is high.

ANISOTROPIC DIFFUSION PARTIAL DIFFERENTIAL EQUATIONS

(a1)

(a2)

(b1)

(b2)

163

F IGURE 7. Example of image denoising by surface area minimization. (a1) Noisy image, (a2) corresponding surface. (b1) Restored image, (b2) corresponding surface.

B. Divergence-Based Diffusion PDEs One level of flexibility for designing regularization PDEs has been reached with the introduction of more generic divergence expressions (Aubert and Kornprobst, 2002; Alvarez et al., 1992; Kornprobst et al., 1997; Sapiro, 2001; Weickert, 1998). Basically, the idea was to replace the function φ  (∇I )/∇I  in the divergence of the scalar-valued PDE (6) by expressions depending on more appropriate image features. On one hand, this gives more freedom to design regularization PDEs that better fit local constraints. On the other hand, the global interpretation of the regularization process is often lost; generally, equations so designed no longer correspond to a functional minimization. Historically, authors (Alvarez et al., 1992) first proposed to use a diffusivity function g(∇I ∗Gσ ) depending on the convolved gradient norm ∇I ∗Gσ , rather than simply considering ∇I  as a measure of image variations, for the

164

TSCHUMPERLÉ AND DERICHE

regularization of scalar-valued images:

  ∂I = div g ∇I ∗ Gσ  ∇I , ∂t where Gσ is a 2D normalized Gaussian function. This has initially been done to ensure that the regularization formulations are well posed. However, it also appeared to respect a more coherent local diffusion geometry by involving a larger neighborhood in the computation of the local image variations that influence the smoothing process. A major generalization of divergence-based equations for scalar and multichannel images has been more recently proposed by Weickert (Weickert, 1996b; Weickert, 1997a; Weickert, 1997b; Weickert, 1998). Basically, the idea consists of considering image pixels as chemical concentrations or temperatures that diffuse with respect to some physical laws (Fick’s law and continuity equations). He proposed this very generic divergence-based equation, parameterized by a field D :  → P (2) of 2 × 2 diffusion tensors: ∂Ii = div(D∇Ii ) (i = 1, . . . , n), (10) ∂t The tensor field D defines a gradient flux and controls the local diffusion behavior of the smoothing process [Eq. (10)]. Note that the ψ-functional formalism described in Section II.A is a particular instance of the PDE (10) with D defined as follows: ∂ψ ∂ψ D= θ+ θ+T + θ− θ−T . ∂λ+ ∂λ− More specifically, Weickert proposed to design the diffusion tensor D for each image point X = (x, y), by selecting its two eigenvectors u, v and eigenvalues λ1 , λ2 as functions of the spectral elements of the smoothed structure tensor Gσ [Eq. (3)] such that ⎧ λ = β, ⎪  ⎨ 1  u = θ + , and β (if λ+ = λ− ), (11)

−C  ⎪ v = θ −, ⎩ λ2 = β + (1 − β) exp else 2 (λ −λ ) +

−

(C > 0 and β ∈ [0, 1] are user-fixed parameters of the method). The tensor D is computed at each image point as: D = λ1 uuT + λ2 vvT . Note that the tensor field D is the same for all image channels Ii , ensuring that all Ii are smoothed by a common multichannel geometry, which takes the correlation between image channels into account (since D depends on Gσ ), contrary to an uncorrelated channel-by-channel approach. Weickert assumed that the tensor shape at each point X = (x, y) of the field D gives the preferred smoothing geometry at this point. The idea behind the choice in Eq. (11) was then:

ANISOTROPIC DIFFUSION PARTIAL DIFFERENTIAL EQUATIONS

165

• On almost constant regions, we have λ+  λ−  0 and then we derive λ1  λ2  β, that is, D  αId . The tensor D is then defined as isotropic in flat regions. • Along image contours, we have λ+  λ−  0

and as a result, λ2 > λ1 > 0.

Here, the diffusion tensor D is anisotropic, directed mainly by the smoothed direction θ − of the image contours. However, it is important that the amplitudes and directions of the local smoothing performed by the divergence-based PDE (10) are not precisely defined by the eigen characteristics (shapes) of the diffusion tensor D at X. This may lead to a smoothing behavior that is not expected, as illustrated by the following simple example. Suppose we want to anisotropically smooth a ∇I scalar image I :  → R everywhere along the gradient direction ∇I  with a constant strength of 1. (This is for illustration purposes, since all image discontinuities would be destroyed by choosing such a smoothing geometry.) Intuitively, we would define D at each point X ∈  as    ∇I T ∇I ∀X ∈ , D(X) = ∇I  ∇I  leading to the simplification of Eq. (10) as   1 ∂I T = div ∇I ∇I ∇I = div(∇I ) = I, ∂t ∇I 2 2

2

∂ I ∂ I where I = ∂x 2 + ∂y 2 stands for the Laplacian of I . As noticed in (Koenderink, 1984), the evolution of this well-known heat flow equation is similar to the convolution √ of the image I by a normalized Gaussian kernel Gσ with a variance σ = 2t. So, this particular choice of anisotropic tensors D leads to an isotropic smoothing behavior, without preferred smoothing orientations. Note that choosing D = Id (identity matrix) would yield the same result: Different tensors fields D with very different shapes (isotropic or anisotropic) may define the same regularization behavior. Actually, the divergence is a differential operator, so Eq. (10) implicitly depends on the spatial variations of the tensor field D. Clearly, the divergence Eq. (10) hampers the design of a significant pointwise smoothing behavior.

C. Oriented Heat Flows Oriented heat flows, also named oriented Laplacian formulations, consider that a local smoothing process can be decomposed into two orthogonal

166

TSCHUMPERLÉ AND DERICHE

F IGURE 8. directions.

Principle of oriented Laplacians. Two 1D smoothings are done along adapted

and unidimensional heat flows respectively oriented along two directions u1 and u2 (these vectors forming an orthonormal basis) associated with two smoothing amplitudes c1 and c2 . The smoothing amplitudes and orientations are naturally different for each image point, since they adapt themselves to the local configuration of the image (Figure 8). The resultant equation is written as the sum of these two heat flows: ∂I = c1 Iu1 u1 + c2 Iu2 u2 , (12) ∂t where u1 and u2 are unit orthogonal vectors and c1 , c2  0. Iu1 u1 and Iu2 u2 denote the second derivatives of I in the directions u1 and u2 and their vector components are formally defined as follows: ∀i = 1, . . . , n,

Iiu1 u1 = uT1 Hi u1

and Iiu2 u2 = uT2 Hi u2 ,

where Hi is the Hessian of Ii , defined on each point X ∈  by ⎛ 2 ⎞ ∂ Ii ∂ 2 Ii   2 ∂x∂y Iixx Iixy ∂x ⎠. =⎝ 2 Hi = ∂ 2 Ii Iixy Iiyy ∂ Ii ∂x∂y

(13)

∂y 2

Here, the diffusion behavior is defined entirely by the knowledge of the smoothing directions u1 , u2 and the corresponding weights c1 and c2 . This formulation was first proposed by (Kornprobst, 1998; Kornprobst et al., 1996) for the regularization of scalar-valued images I , with the following choices for c1 , c2 and u1 , u2 :   ∇I ⊥ u1 = ξ = ∇I c1 = 1, , and ∇I c2 = g(∇I ∗ Gσ ), u2 = η = ∇I  , where g : R → R is a function decreasing to 0 (the pixel diffusion must vanish on high gradients). It allows a permanent anisotropic smoothing along

ANISOTROPIC DIFFUSION PARTIAL DIFFERENTIAL EQUATIONS

167

the edges ξ , even on very high gradients, since c1 = 1 everywhere in the image. The general formalism of oriented Laplacians [Eq. (12)] also allows other well-known equations to be calculated, such as the mean curvature flow ∂I ∂t = Iξ ξ , obtained with c1 = 1, c2 = 0, u1 = ξ and v2 = η (Deriche and Faugeras, 1997). Note that this was not possible with divergence-based expressions [Eq. (10)]. Other works for scalar image regularization have used similar variants of this technique (Carmona and Zhong, 1998; Krissian, 2000). Sapiro and Ringach (1996) proposed an extension of the mean curvature flow It = Iξ ξ for multichannel images, using an oriented Laplacian formulation [Eq. (12)]. They naturally used the Di Zenzo attributes to incorporate informations on the multichannel geometry in their proposed equation: ∂I = g(λ+ − λ− )Iθ− θ− , (14) ∂t where g : R → R is a positive decreasing function, avoiding the smoothing of high gradients regions. It was one of the first attempts to construct an oriented Laplacian PDE directly from a local vector-valued geometry. Indeed, all channels Ii are smoothed along a common vector edge direction with a common intensity. However, some drawbacks remain: √ • The coherence norm N− = λ+ − λ− is used here as a measure of vectorvalued variations in order to reduce diffusion on image contours. This may not be a good choice since some corner structures do not respond to the norm N− , as explained in Section I.A, and will then be oversmoothed. • In flat regions (N− → 0), the diffusion is made along a single direction θ− , which is directed mainly by the noise since no coherent structures exist in these regions. Undesired texture effects result from this monodirectional smoothing. This is particularly true here, since contrary to decoupled regularizations, vector components are not blended with this method (the diffusions in all image channels Ii follow a common direction). Isotropic smoothing would be more adapted to remove noise in such flat regions. D. Trace-Based Diffusion PDEs A simple generalization of oriented Laplacians has been proposed (Tschumperlé, 2002; Tschumperlé and Deriche, 2003; Tschumperlé and Deriche, 2005). The concept relies on the use of a generic diffusion tensor field T :  → P (2) to describe the diffusion geometry of Eq. (12), instead of separately describing local directions θ+ , θ− and amplitudes c1 , c2 of smoothing. Actually, the proposed equation was simply a rewrite of the previous PDE (12), using a trace operator: ∂Ii = trace(THi ), ∀i = 1, . . . , n, (15) ∂t

168

TSCHUMPERLÉ AND DERICHE

where Hi is the Hessian of Ii [Eq. (13)] and the tensor field T is computed as: ∀X ∈ ,

T(X) = c1 u1 uT1 + c2 v2 vT2 .

Note that in this case, each channel Ii of I is also smoothed with a common tensor field T. Equations (12) and (15) are strictly equivalent, but Eq. (15) clearly separates the smoothing geometry (defined by the tensor field T) from the smoothing itself. This is close to Weickert’s method that led to the divergence PDE (10); the regularization problem now simplifies to the design of a tensor field T adapted to the considered application. But in the case of trace-based PDEs, the tensor field that defines the local smoothing behavior has the interesting property of unicity: Two different tensor fields will necessarily lead to different smoothing behaviors. Indeed, Eq. (15) has a simple geometric interpretation in terms of local filtering with oriented Gaussian kernels. Indeed, let us consider first that T is a constant tensor field. It can then be demonstrated that the formal solution of the PDE (15) is (see Appendix B for details): Ii(t) = Ii(t=0) ∗ G(T,t)

(i = 1, . . . , n),

(16)

where ∗ designates for the convolution operator and G(T,t) is an oriented Gaussian kernel, defined by   XT T−1 X 1 (T,t) exp − with X = ( x y )T . (X) = (17) G 4π t 4t This is in fact a generalization of Koenderink’s idea (Koenderink, 1984), who proved this property in the field of computer vision for the simpler case of the isotropic diffusion tensor T = Id , resulting in the well-known heat flow i equation: ∂I ∂t = Ii . Figure 9 illustrates three Gaussian kernels G(T,t) (x, y) respectively obtained with isotropic and anisotropic tensors T and the corresponding evolutions of the diffusion PDE (15) on a color image. The Gaussian kernels G(T,t) produce the classical representations of the tensors T with ellipses. Conversely, it is clear that the tensors T represent the exact geometry of the smoothing performed by the PDE (15). When T is not constant (which is generally the case), that is, it represents a field  → P(2) of variable diffusion tensors, the PDE (15) becomes nonlinear but can be viewed as the application of temporally and spatially varying local masks GT,t (X) over the image I. Figure 10 illustrates three examples of spatially varying tensor fields T, represented with fields of ellipsoids, and the corresponding evolutions of Eq. (15) on a color image. As before, the shape of each tensor T gives the exact geometry of the local smoothing process

ANISOTROPIC DIFFUSION PARTIAL DIFFERENTIAL EQUATIONS

F IGURE 9. kernels.

169

Trace-based PDEs [Eq. (15)] viewed as convolutions by oriented 2D Gaussian

F IGURE 10.

Trace-based PDEs [Eq. (15)] with nonconstant diffusion tensor fields T.

performed by the trace-based PDE (15) point by point. Because the trace is not a differential operator, the local interpretation of the smoothing process as a convolution with an oriented Gaussian mask is valid here. In the same way that structure tensors code for each image pixel X the main directions of the edges θ− as well as the edge strength λ+ + λ− , the diffusion tensor field T similarly will code the preferred local smoothing directions, as well as the desired smoothing amplitudes along these directions,

170

TSCHUMPERLÉ AND DERICHE

for each image pixel X. Of course, T(X) depends on the local geometry of I and is thus defined from the spectral elements λ− , λ+ and θ − , θ + of the smoothed structure tensor Gσ . For image denoising purposes, the choice proposed by (Tschumperlé, 2002; Tschumperlé and Deriche, 2003; Tschumperlé and Deriche, 2005) is 1 (1 + λ+ + λ− )p1 1 = (1 + λ+ + λ− )p2

− c1 = f(λ = + ,λ− )

and

+ c2 = f(λ + ,λ− )

with p1 < p2 ,

(18)

where p1 , p2 ∈ R are parameters of the proposed method. At this point, the desired smoothing behavior is intended to be: • If a pixel X is located on an image contour (λ+ (X) is high), the smoothing on − X would be performed primarily along the contour direction θ(X) (since + − f(.,.)  f(.,.) ), with a smoothing strength inversely proportional to the contour strength. • If a pixel X is located on a homogeneous region (λ+ (X) is low), the smoothing on X would be performed in all possible directions (isotropic smoothing), + − since f(.,.)  f(.,.) and then T  Id (identity matrix). This is one possible choice for f − , f + in order to satisfy basic image denoising requirements. Actually, it is quite natural to design a smoothing behavior from the image structure before applying the regularization process itself. The trace-based Eq. (15) is a good attempt to separate the smoothing geometry from the smoothing process itself, while providing a geometric interpretation on how the smoothing is performed. It proved some natural links between PDEs and other local filtering techniques, such as bilateral filtering (Barash, 2000; Tomasi and Manduchi, 1998). Another similar approach based on non-Gaussian convolution kernels has been also proposed for the specific case of the Beltrami flow (Sochen et al., 2001). Nevertheless, the fact that the trace equation [Eq. (15)] behaves locally as an oriented Gaussian smoothing whose strength and orientation is directly related to the tensor T(X) has a major drawback. Indeed, on curved structures (like corners), this Gaussian behavior is not desirable: When the local variation of the edge orientation θ − is high, a Gaussian filter tends to round corners, even by conducting the smoothing only along θ − . This is due to the fact that an oriented Gaussian mask is not curved itself. This classical behavior is best known as the “mean curvature flow” effect, characterized by ∂2I the PDE ∂I . This problem is illustrated in Figures 11b and 12b where −2 ∂t = ∂θ

ANISOTROPIC DIFFUSION PARTIAL DIFFERENTIAL EQUATIONS

(a)

(b)

(c)

171

(d)

F IGURE 11. Problems encountered using trace-based PDEs [Eq. (15)] on curved image structures. (a) Noisy synthetic color image; (b) applying trace-based PDE (15), with p1 = 0.5, p2 = 1.2; (c) applying trace-based PDE (15), with p1 = 0.9, p2 = 1.2; (d) applying our constrained PDE (29), with p1 = 0.5, p2 = 1.2.

Eq. (15) has been applied on synthetic and real color images and T has been defined as in Eq. (18) (then f − = 0). It is obvious that image structures are subject to the mean curvature flow effect, resulting in rounding the corners of the square in Figure 11b, or in blending parallel thin curved structures in Figure 12b. To avoid this oversmoothing effect, one may try to stop the action of the diffusion PDE on corners (by vanishing tensors T(X) there, that is, f − = f + = 0). But this implies the detection of curved structures on noisy or corrupted images, which is generally imprecise in the presence of noise, even when using the Di Zenzo geometry. Conversely, image undersmoothing on edges may occur when limiting the diffusion too much in regions with high-intensity variations (Figure 11c). There is a difficult trade-off between complete noise removal and preservation of curved structures when using trace-based PDEs [Eq. (15)]. These types of regularization processes are not concerned about the curvature of the smoothing directions, and by extension, of the curvature of the image contours. Taking this curvature into account is very desirable and has motivated the work presented later. In Section III, we propose a class of trace-based regularization PDEs that smooth an image I along a tensor field T, while implicitly taking curvatures of specific integral curves of T into account. Generally, the method locally filters the image with curved Gaussian kernels when necessary in order to better

172

TSCHUMPERLÉ AND DERICHE

(a)

(b)

(c)

F IGURE 12. Comparisons between trace-based PDEs [Eq. (15)] and our new curvature-preserving PDEs [Eq. (29)] on a real image. (a) Image of a fingerprint; (b) applying trace-based PDE (15), with p1 = 0.5, p2 = 1.2; (c) applying our constrained PDE (29), with p1 = 0.5, p2 = 1.2.

preserve image structures. For comparison purposes, results of the curvaturepreserving equation are shown on Figures 11d and 12c. E. Links Between Existing Regularization Methods The link between these three formulations is generally not trivial, especially for vector-valued images. It is well known for the classical case of φfunctional regularization of scalar images (n = 1). One can start from a regularizing functional minimization (A) and find the corresponding divergencebased (B) and oriented Laplacian (C) formulations as follows: 

 (A): min φ ∇I  d I :→R



⇒

(B):

⇒

(C):

   φ (∇I ) ∂I = div ∇I ∂t ∇I 

 ∂I φ  (∇I ) = Iξ ξ + φ  ∇I  Iηη , ∂t ∇I 

(19)

where η = ∇I /∇I  and ξ ⊥η. Note that this regularization generally leads to anisotropic smoothing (in the sense that it is performed in privileged spatial directions with different weights), despite the isotropic shape of the

ANISOTROPIC DIFFUSION PARTIAL DIFFERENTIAL EQUATIONS

173



(∇I ) underlying tensor D = φ ∇I  Id in the divergence expression. Also note that this global-to-local path (from variational to trace-based equations) can rarely be followed in the inverse order. For multichannel images, this link also can be found:  (A): min n ψ(λ+ , λ− ) d I : →R



⇒

(B):

⇒

(C):

∂Ii ∂ψ ∂ψ = div(D∇Ii ), where D = θ+ θ+T + θ− θ−T ∂t ∂λ+ ∂λ− n   

∂Ii (20) = trace δij D + Qij Hj , ∂t j =1

where the δij is the Kronecker symbol (δij = 1 when i = j , and 0 elsewhere), Qij designates a family of n2 tensors (i, j = 1, . . . , n), defined as the T symmetric parts of the following matrices Pij (i.e., Qij = (Pij + Pij )/2):   ∂α ∂α ij T T T P = α∇Ii ∇Ij Id + 2 θ+ θ + + θ− θ− ∇Ij ∇IiT G ∂λ+ ∂λ−      ∂β ∂β θ+ θ+T + α + θ− θ−T ∇Ij ∇IiT +2 α+ ∂λ+ ∂λ− with α=

f1 (λ+ , λ− ) − f2 (λ+ , λ− ) λ+ − λ −

and

β=

λ+ f2 (λ+ , λ− ) − λ− f1 (λ+ , λ− ) λ+ − λ −

and f1 =

∂ψ ∂λ+

and f2 =

∂ψ . ∂λ−

The development (A) ⇒ (B) from the functional to the divergence formulation is detailed in Appendix A. The development (B) ⇒ (C) from the divergence to the trace-based equation is detailed in Appendix C. This last development initially proposed by (Tschumperlé and Deriche, 2003; Tschumperlé and Deriche, 2005) unifies an entire range of previously proposed vector-valued regularization algorithms (variational and divergence based PDEs) into an extended trace-based equation, composed of several channel-diffusion contributions that have direct geometric interpretations in terms of local filtering by Gaussian kernels. Although the geometric interpretation of the overall sum of trace equations is not direct, it is interesting to see that additional diffusion tensors Qij clearly appear in the trace expressions, and contribute to modify the inner tensor D, which is not representative

174

TSCHUMPERLÉ AND DERICHE

of the smoothing behavior. More generally, tensors appearing in traces and divergences lead to different smoothing behaviors.

III. C URVATURE -P RESERVING PDE S The framework of curvature-preserving PDEs, first introduced by (Tschumperlé, 2006), defines a specific variant of multichannel diffusion PDEs. Its goal is to provide a generic tensor-driven regularization method as the divergence-based PDE (10) and trace-based PDE (15), but also focuses on the preservation of thin curved structures. We review this very efficient formalism and show how it can be understood from a local smoothing geometry viewpoint. A. The Single-Direction Case To illustrate the general idea of curvature-preserving PDEs, we first focus on image regularization along a vector field w :  → R2 instead of a tensor field T. We then consider a local smoothing everywhere along a single w direction w , with a smoothing strength w. The two spatial components of w are denoted by w(X) = ( u(X) v(X) )T . The curvature-preserving regularization PDE that smoothes I along w is defined by:

 ∂Ii = trace wwT Hi + ∇IiT Jww , ∂t where Jw is the Jacobian of w, and Hi is the Hessian of Ii . ⎛ 2 ⎞ 2   ∀i = 1, . . . , n,

Jw =

∂u ∂x ∂v ∂x

∂u ∂y ∂v ∂y

and

Hi = ⎝

∂ Ii ∂x 2 ∂ 2 Ii ∂x∂y

∂ Ii ∂x∂y ∂ 2 Ii ∂y 2

(21)

⎠.

The PDE (21) adds a term ∇IiT Jww to the trace-based Eq. (15) that smoothes I along w with locally oriented Gaussian kernels (see Section II.D). This extra term naturally depends on the variation of the vector field w. The following shows how Eq. (21) is related to w. X be the curve defining the integral curve of w, starting from X and Let C(a) parametrized by a ∈ R: ⎧ X = X, ⎨ C(0) (22) X ⎩ ∂ C(a) = w(C X ). (a) ∂a

ANISOTROPIC DIFFUSION PARTIAL DIFFERENTIAL EQUATIONS

(a)

175

(b)

F IGURE 13. Integral curve C X of a vector field w :  → R2 . (a) Integral curve of a general field w. (b) Example of integral curves when w is the lowest eigenvector of the structure tensor G of a color image I (one block is one color pixel).

X is tracked forward, and backward When a → +∞, the integral curve C(a) when a → −∞ (Figure 13). F denotes the family of integral curves of w. X around a = 0 is: A second-order Taylor development of C(a) X X C(h) = C(0) +h

X ∂ C(a)

∂a

= X + hw(X) +

|a=0

+

2 X

 h2 ∂ C(a) + O h3 2 2 ∂a |a=0

 h2 Jw(X) w(X) + O h3 , 2

with h → 0, and O(hn ) = hn n . Then, we can compute a second-order Taylor X ) around a = 0, which corresponds to the variations of development of Ii (C(a) the image intensity near X when following the integral curve C X :  



X  h2 Ii C(h) = Ii X + hw(X) + Jw(X) w(X) + O h3 2   h T = Ii (X) + h∇Ii (X) w(X) + Jw(X) w(X) 2 +



 h2 trace w(X) wT(X) Hi (X) + O h3 . 2

The term trace(w(X) wT(X) Hi (X) ) = derivative of Ii along w.

∂ 2 Ii ∂w2

corresponds to the second directional

176

TSCHUMPERLÉ AND DERICHE

X ) at a = 0 is then: The second derivative of the function a → Ii (C(a) X ) ∂ 2 Ii (C(a)   ∂a 2 

= lim a=0

h→0

X 

X  1 X  Ii C(h) + Ii C(−h) − 2Ii C(0) h2

1 2 T h ∇Ii Jw(X) w(X) h2



 + h2 trace w(X) wT(X) Hi (X) + O h3

 = trace w(X) wT(X) Hi (X) + ∇IiT Jw(X) w(X) . = lim

h→0

(23)

Note that this is precisely the right term in the proposed curvature-preserving PDE (21). Eq. (21) can be seen individually for all integral curves of F instead of each point X ∈ ; consider another point Y ∈ C X . Then there exist  ∈ R such X . Indeed, C X and C Y describe the same curve [Eq. (22)] with that Y = C() Y = CX different parameterization: ∀a ∈ R, C(a) (+a) . As Eq. (21) is verified on ∂Ii (C X )

∂ 2 Ii (C X )

Y, then ∂t(a) |a= = ∂a 2(a) |a= . This is obviously true for  ∈ R since Eq. (21) is verified for all points Y lying on the integral curve C X . Then, the PDE (21) may be also written as ∀C ∈ F , ∀a ∈ R,

∂ 2 Ii (C(a) ) ∂Ii (C(a) ) = . ∂t ∂a 2

(24)

We may recognize in Eq. (24) a one-dimensional heat flow constrained on C . This is very different from a heat-flow oriented by w, as in the ∂ 2 Ii i formulation ∂I ∂t = ∂w2 since the curvatures of integral curves of w are now implicitly taken into account. In particular, the constrained Eq. (21) has the interesting property of vanishing when image intensities are perfectly constant on the integral curve C , regardless of the curvature of C . In this context, defining a field w that is tangent everywhere to the image structures will allow the preservation of these structures, even if they are curved (such as corners). This is not the case with divergence or trace-based PDEs [Eqs. (10) and (15)]. This curvature-preserving property of Eq. (21) is illustrated in Figures 11d and 12b). The constrained Eq. (21) is an elliptic PDE since the matrix wwT is positive definite. The existence and unicity of the solutions of Eq. (21) are not directly approached here. Section III.B, shows that its solution can be approximated by the technique of LICs, which is a well-posed analytical approach.

ANISOTROPIC DIFFUSION PARTIAL DIFFERENTIAL EQUATIONS

177

B. Curvature-Preserving PDEs and Line Integral Convolutions LICs were first introduced (Cabral and Leedom, 1993) as a technique to render a textured image ILIC that represents a vector field w :  → R2 . The idea, originally expressed under a discrete formula, consists in smoothing an image Inoise (containing only noise) by averaging its pixel values along the integral curves of w. A continuous formulation of an LIC is then: ∀X ∈ ,

ILIC (X)

1 = N

+∞

X  f (p)Inoise C(p) dp

(25)

−∞

where f : R → R is an even function (strictly decreasing to 0 on R+ ) and C X is defined as the integral curve [Eq. (22)] of w through X. The normalization factor N allows the preservation of the average pixel value along C X and is  +∞ equal to N = −∞ f (p) dp. As noted in Section III.A, the curvature-preserving PDE (21) can be seen as the one-dimensional heat flow [Eq. (24)] constrained on the integral curve X ), Eq. (24) can be C X ∈ F . Using the variable substitution L(a) = I(C(a)  [t] also written as ∂L ∂t (a) = L(a) . The solution L at time t is known to be the convolution of L[t=0] by a normalized Gaussian kernel Gt (see Deriche and Faugeras, 1997; Koenderink, 1984): L[t] (a)

+∞ = L[t=0] (p) Gt (a−p) dp

with Gt (p)

−∞

 p2 . (26) =√ exp − 4t 4π t 1



X = X and Substituting L in Eq. (26) with a = 0, and remembering that C(0) Gt (−p) = Gt (p) :

∀X ∈ ,

I[t] (X)

+∞

X  = I[t=0] C(p) Gt (p) dp.

(27)

−∞

Eq. (27) is a particular form of the continuous LIC-based formulation [Eq. (25)] with a Gaussian weighting function f = Gt . Here, the nor+∞ malization factor is N = −∞ Gt (p) dp = 1. Intuitively, the evolution of the curvature-preserving PDE (21) may be seen as the application of local convolutions by normalized one-dimensional Gaussian kernels along integral curves C of w. Thus this type of anisotropic image smoothing considers a curved filtering, instead of just an oriented one. Applying this setting on a multichannel image I, with w being the lowest eigenvector of the structure tensor field G (i.e., the contour direction) allows

178

TSCHUMPERLÉ AND DERICHE

the anisotropic smoothing of I with edge preservation, even if these edges are curved. This is shown in Figure 13b, where few integral lines C X are computed around a typical T-junction structure. Note how the streamlines rotate when arriving at the junction, with a subpixel precision. The streamlines have been computed with a fourth-order Runge–Kutta scheme. Eq. (27) is an analytical solution of Eq. (21) when w does not evolve over time. This property is generally not verified when dealing with general nonlinear regularization PDEs, where the smoothing geometry is reevaluated at each time step (thus defining a temporal nonlinearity). To eliminate this nonlinearity, we perform several successive iterations of the LIC scheme [Eq. (27)], where the vector field w is updated at each iteration. This is a good method of approximating Eq. (21). Classical explicit schemes usually consider the smoothing geometry w as constant between two successive PDE iterations I[t] and I[t+dt] . Thus, the curvature-preserving Eq. (21) is efficiently discretized by several iterations of the LIC formulation [Eq. (27)]. This is detailed in Section IV. C. Between Traces and Divergences The curvature-preserving PDE (21) may be compared to trace and divergence expressions [Eqs. (10) and (15)], for the case of single-direction smoothing T = wwT . In this case, the divergence PDE (10) may be developed as follows:   2 ∂Ii i u ∂x + uv ∂I

T  ∂y div ww ∇Ii = div 2 ∂Ii i uv ∂I ∂x + v ∂y   ∂ 2 Ii ∂ 2 Ii ∂ 2 Ii + v2 2 = u2 2 + 2uv ∂x∂y ∂x ∂y  ∂u ∂v ∂u  2u ∂x + u ∂y + v ∂y + ∇IiT ∂v ∂u 2v ∂v ∂y + u ∂x + v ∂x  ∂u   ∂u ! u ∂x + v ∂u u ∂x + u ∂v

T  ∂y ∂y T = trace ww Hi + ∇Ii + ∂v ∂u u ∂x + v ∂v v ∂x + v ∂v ∂y ∂y

T  = trace ww Hi + ∇IiT Jww + div(w)∇IiT w. This equation is a sum of three different terms that have the following meaning: • The first term corresponds to the trace PDE (15), which smoothes locally I along w, using oriented Gaussian kernels.

ANISOTROPIC DIFFUSION PARTIAL DIFFERENTIAL EQUATIONS

179

• The two first terms correspond to the curvature-constrained regularization PDE (21), which smooths locally I along w while taking the curvature of integral curves C of w into account. • The three terms together correspond to the classical divergence PDE (10), which performs local diffusions of I along w. This last term, div(w)∇IiT w, is mainly responsible for the perturbations of the effective smoothing direction, as described in Section II.B. It is not desirable for image regularization purposes. It is interesting to observe that the curvature-constrained PDE (21) is then “mathematically” positioned between the trace [Eq. (15)] and divergence formulations [Eq. (10)], and allows at the same time to consistently smooth along the predefined smoothing directions w, while preserving curved images structures. The curvature-preserving PDE (21) can also be written as a divergencebased PDE minus a constraint term:  

trace wwT Hi + ∇IiT Jww = div wwT ∇Ii − div(w)∇IiT w. Two particular cases of directions w merit study, in the case of scalar-valued images (n = 1): ⊥

∇I T • When w = ∇I  (isophote direction), then ∇I Jww = −Iww , eliminating the velocity of the curvature-preserving evolution Eq. (21), by counterbalancing the trace-based term (which is nothing more than the mean curvature motion in this case). No smoothing is performed. This is natural since pixels along the isophotes have constant values, so averaging those values should not modify the image. Note by comparison that the velocity of the corresponding divergence-based expression div(wwT ∇Ii ) also is eliminated here. ∇I T • When w = ∇I  (gradient direction), then ∇I Jww = 0, and the velocity of the curvature-preserving PDE (21) becomes simply Iww , which corresponds to a smoothing of the image along the gradient direction [the same as the unconstrained trace-based PDE (15)]. Note by comparison that the velocity of the corresponding divergence-based expression is I in this case, which corresponds to an isotropic smoothing of the image, instead of an anisotropic one.

These two particular cases allow better understanding of the difference of regularization behaviors among the trace, divergence, and curvaturepreserving formulations. Note that when w is a divergence-free field (i.e., div(w) = 0), the divergence-based PDE (10) and the curvature-preserving formulation Eq. (21) are strictly equivalent. This is very rarely the case.

180

TSCHUMPERLÉ AND DERICHE

D. Extension to Multidirectional Smoothing In (Tschumperlé, 2006), the single-direction smoothing PDE (21) has been extended so that it can be applied to a tensor-valued geometry T :  → P(2), instead of a vector-valued geometry w. This is important since a diffusion tensor describes much more complex smoothing behaviors than single directions. In particular, it may represent both anisotropic or isotropic regularization behaviors. The extension of the curvature-preserving PDE (21) is not straightforward; the notions of curvature and integral curves of tensorvalued fields T are not as natural as with direction fields w. To resolve this problem, we proposed to locally decompose a tensordriven smoothing process into several vector-driven smoothing processes along different orientations. We first notice that π aα aαT dα =

π Id 2

 where aα =

cos α sin α

 .

α=0

Then, any 2 × 2 tensor T may be written as follows:  π  √ 2√ T aα aαT dα T, T= π α=0

  √ T is the square root of T = f + uuT + where T = f + uuT + f − vv√ √ √ − T f vv . It can easily be verify that ( T )2 = T and ( T )T = T. Thus, the tensor T may be decomposed as 2 T= π

π √ α=0

Taα aαT

√

2 T dα = π T

π √ √ ( Taα )( Taα )T dα.

(28)

α=0

√ √ The tensor T has been split into a sum of atomic tensors ( Taα )( Taα )T ; each vector √ is purely anisotropic and directed only along the direction of the vector Taα ∈ R2 . Eq. (28) suggests a means to decompose any tensordriven regularization PDE into a sum of single-direction smoothing processes, each of them respecting the overall geometry T. For instance: • If √ T = Id (identity matrix), the tensor is isotropic and: ∀α ∈ [0, π ], Taα = aα . The resulting smoothing will be then performed in all directions aα of the plane with the same strength. • If T =√uuT (where u ∈ S1 ), the tensor is purely anisotropic and: ∀α ∈ [0, π ], Taα = (uT aα )u. The resulting smoothing will be then performed only along the direction u of the tensor T.

ANISOTROPIC DIFFUSION PARTIAL DIFFERENTIAL EQUATIONS

181

Then, using Eq. (28) and considering that each single-direction smoothing must be done with a curvature-preserving approach [Eq. (21)], the following curvature-constrained regularization PDE results, acting on a multichannel image I :  → Rn and driven by a tensor-valued smoothing geometry T: ∀i = 1, . . . , n, π √ √

√  2 ∂Ii = trace ( Taα )( Taα )T Hi + ∇IiT J√Taα Taα dα, ∂t π α=0

which can be simplified as ∀i = 1, . . . , n,

2 ∂Ii = trace(THi ) + ∇IiT ∂t π

π

√ J√Taα Taα dα, (29)

α=0

J√Taα

stands for the Jacobian of the where aα = ( cos α sin α ) , and √ vector field  → Taα . This type of smoothing decomposition along all orientations of the plane can be also found in (Weickert, 1994). As in the single-direction smoothing case, Eq. (29) may be seen as a trace-based equation [Eq. (15)], where an extra term has been added to respect the curvature of all integral lines passing through the tensor-valued geometry T. T

IV. I MPLEMENTATION C ONSIDERATIONS To implement the regularization method in Eq. (29), the LIC-based interpretation of curvature-preserving PDEs presented in Section III.B can be used. Indeed, we can explicitly discretize Eq. (29) by the following Euler scheme:  N −1   √ 2dt R( Taα ) I[t+dt] = I[t] + N k=0

where α = kπ/N (in the interval [0, π ]), dt is the usual temporal discretization step, and R(w) represents a discretization of the unidirectional smoothing PDE velocity [Eq. (21)] that preserve curvatures along a vector √ field w. If this N −1 [t] expression is written as I[t+dt] = N1 ( k=0 I + 2dt R( Taα )), it may be expressed √ as the averaging of different Gaussian-pondered LICs along vector fields Taα :   N−1  [t] 1 √ , I[t+dt] = I LIC( Taα ) N k=0

182

TSCHUMPERLÉ AND DERICHE

where each Gaussian variance has a standard deviation dt. Basically, the difficulty lies is the LIC computation, which requires the tracking of integral curves of a vector field. Here we used a very simple method based on the classical Runge–Kutta (Press et al., 1992) integration scheme. Faster LIC implementations have been proposed by (Stalling and Hege, 1995) but do not deal with Gaussian-pondering functions, as needed here. This simple observation leads then to the following fast algorithm for the implementation of one iteration of the curvature-preserving PDE (29): • Compute the smoothed structure tensor field Gσ from I[t] : ⎛ ∂I [t] 2

∂Ii[t]  ∂Ii[t]  ⎞ i n  ∂x ∂x ∂y ⎝ ⎠ Gσ = Gσ ∗

∂Ii[t]  ∂Ii[t] 

∂Ii[t] 2 i=1

∂x

∂y

∂y

σ will depend on the noise scale. We used relatively low values (between 0 and 1.5) for our experiments in Section V. • Compute the eigenvalues λ+ , λ− and eigenvectors θ + , θ − of Gσ . • Compute the smoothing geometry tensor field T from Gσ : T=

1 1 T T θ −θ − + θ +θ + . (1 + λ+ + λ− )p1 (1 + λ+ + λ− )p2

• For all α in [0, π ] (discretized with √ a user-fixed step dα ): – Compute the vector field w = Taα . – Perform an LIC of I[t] along C X in the forward and backward directions. • Average all LICs computed in previous steps. The main parameters of the algorithm are p1 , p2 , σ , and dt and the number of PDE iterations nb that are applied. The characteristics of this scheme, compared to the classical finite-difference model are as follows: • It allows the preservation of thin image structures from a numerical point of view; the smoothing is performed along integral curves of w, with a subpixel accuracy. Precise fourth-order Runge–Kutta interpolation (Press et al., 1992) is used to track the integral curves C in the image. • It allows the choice of very large time steps dt, since the scheme we proposed is unconditionally stable. Indeed, dt is simply proportional to the overall smoothing variance of the Gaussian-pondering convolutions done along C ∈ F . • As a result, the regularization algorithm performs very fast. Very few iterations are necessary to obtain the result, even if each iteration is more time consuming. For our applications (presented in Section V), we were able to choose only nb = 1 iteration with very large time steps dt. In fact,

ANISOTROPIC DIFFUSION PARTIAL DIFFERENTIAL EQUATIONS

(a)

(b)

183

(c)

F IGURE 14. Comparisons between classical explicit PDE schemes and LIC-based implementation of the PDE (29). (a) Noisy color image, (b) regularization using a finite-difference scheme (stopped at t = 100), (c) regularization using the LIC-based scheme (stopped at t = 100).

this leads to a rough approximation of Eq. (29) since we lost the temporal nonlinearity property of the PDE. But for images with little noise, this gave surprisingly good results. The spatial nonlinearity seems to play a more important role than the temporal nonlinearity in the PDE evolution. The smoothing is done as an averaging of multiple LICs in different directions α. The choice of the discretization step dα is important in this context. In regions where the smoothing needs to be primarily anisotropic, only a few values of α are necessary since in all cases, the smoothing is done along the same single direction. But in homogeneous regions needing isotropic smoothing, a smaller dα gives much better results. Practically speaking, we chose dα = 45◦ , which is sufficient for good precision for isotropic smoothing. Figure 14 shows the efficiency of the scheme compared with the classical finite-difference one. A synthetic noisy image is anisotropically smoothed with the PDE (29), with p1 = 0.01 and p2 = 100 (smoothing mostly along isophotes θ− , with a strength of 1). The LIC-based scheme (Figure 14c) clearly better preserves the structure along time t. This is due to the important role played by the subpixel accuracy property of the underlying LIC computation.

V. A PPLICATIONS A wide variety of problems can be resolved by the regularization PDEs proposed throughout this chapter. Particularly, the curvature-preserving regularization approach described in Section III is used in the following examples dealing with color images. We show results of color image denoising, inpainting, and resizing by nonlinear interpolation. Given processing times have been obtained on a 2.8 GHz i686 Intel processor.

184

TSCHUMPERLÉ AND DERICHE

A. Color Image Denoising and Artifact Removal Image denoising is a direct application of regularization methods. Sensor inaccuracies, digital quantifications, or compression artifacts are indeed some of the various noise sources that can affect a digital image, and suppressing them is a desirable goal. Such artifacts generally lead to small random variations that affect pixels of the image and that must be removed. Figures 15 through 19 show how the curvature-preserving PDE framework [Eq. (29)] can be successfully applied to remove such artifacts while preserving the essential structures of the processed images. • Figure 15 shows an application of the regularization method on the 512 × 512 Baboon color image, artificially degraded by adding uncorrelated Gaussian noise on (R, G, B). This color image has been then denoised with Eq. (29). Thanks to the proposed LIC-based numerical implementation, only one PDE iteration has been necessary to denoise the image, with parameters p1 = 0.5, p2 = 0.7, σ = 1.5, and dt = 50. Processing time is 19.3 seconds for the entire image. • Figure 16 illustrates a real case where a color photograph has been digitized from a grainy paper, leading to the apparition of watered effects on the digital picture. Using the PDE-based regularization method [Eq. (29)] allows removal of the grains while preserving fine structures (palm tree leaves). The image shown is a 152 × 133 portion of the original one. Only one PDE iteration was necessary, with p1 = 0.5, p2 = 0.7, σ = 1, and dt = 10. Processing time is 11 seconds for the entire image (size 586 × 367). • Figure 17 shows the restoration of a digital photograph shot under low-light conditions by a cellular telephone. Such devices usually have poor-quality cameras, leading to significant digital noise (more precisely, Poisson noise) on the acquired color images. The processed color image has a size of 262 × 280 and was restored in 4 seconds (one PDE iteration), with parameters p1 = 0.2, p2 = 0.5, σ = 2, and dt = 120. Note that the curvature-preserving PDE (29) adapts itself locally to the multichannel image geometry in order to preserve thin structures while removing the noise quite well. • Image regularization can also be useful when dealing with other types of noise. The enhancement of JPEG-compressed images is an example of interest. Figure 18 illustrates the suppression of compression artifacts in color images. A JPEG-compressed color image with a size 283 × 249, where the JPEG quality ratio has been set to 10%, is processed by the multichannel image regularization algorithm. Usual block effects inherent to the Discrete Cosine Transform (DCT) compression are visible on the

ANISOTROPIC DIFFUSION PARTIAL DIFFERENTIAL EQUATIONS

185

F IGURE 15. Denoising of the color image Baboon corrupted with artificial Gaussian noise. Noisy color image (left), denoised image (middle), zoom on the eye (right).

compressed image. One PDE iteration is applied then, with p1 = 0.5, p2 = 0.9, σ = 2, and dt = 200, to obtain the regularized result (right). Processing time is 5.4 seconds for the entire image. • Figure 19 illustrates another regularization with a different type of noise. In this case, the regularization method is used to improve a digital true-color image quantified in 256 colors by the Floyd–Steinberg algorithm (size = 355 × 287), which introduces some dithering effects in the image. One PDE iteration has been applied here, with p1 = 0.5, p2 = 0.8, σ = 1, and dt = 30. A 136 × 118 portion of the image is shown. Processing time is 12.8 seconds for the entire image. This kind of reconstruction may be interesting for image-compression algorithms, allowing them to retrieve true color images, even by storing them with a quantified palette. B. Color Image Inpainting Image inpainting is a very new and challenging application that consists of filling in missing image regions (defined by the user) by guessing pixel values such that the reconstructed image still looks natural. Basically, the user provides one color image I :  → R3 , and one mask image M :  → {0, 1}. The inpainting algorithm must fill in the regions where M(X) = 1, by means of some intelligent interpolations. For example, inpainting algorithms can be used to remove various structures in images (scratches, logos, or real objects) that are usually bigger than other image artifacts. Image inpainting was first proposed as a method based on a variational formulation by Masnou and Morel (Masnou and Morel, 1998), followed by many solutions based on diffusion or transport PDEs (Bertalmio et al., 2000b; Bertalmio et al., 2003; Chan and Shen, 2000b; Chan and Shen, 2001; Chan et al., 2002; Tschumperlé

186

TSCHUMPERLÉ AND DERICHE

F IGURE 16. Denoising of the color image Tunisia desert, containing watered effects. Noisy color image (left), denoised image (right); details are shown on the bottom row.

and Deriche, 2003; Tschumperlé and Deriche, 2005). Other papers related to inpainting without use of PDEs (Ashikhmin, 2001; Criminisi et al., 2003; Jia and Tang, 2003; Wei and Levoy, 2000) among others. In this chapter, we view the inpainting process as a direct application of the proposed curvature-preserving PDE (29). We apply the diffusion equation only on the regions to inpaint, allowing the neighbor pixels of these regions to diffuse inside; a nonlinear completion of the image data along isophote directions is thus naturally done, reconstructing the missing parts of the image, since the performed smoothing tries to follow a coherent multivalued image geometry computed from the known parts of the image. The concept of inpainting with PDEs is shown in Figures 20 through 24. • Figure 20 illustrates a simple case of removal of text from a color image, by guessing the pixel colors behind the text. The mask used here is easily detected from the degraded image by considering only green pixels. The result of the image inpainting is shown on the right. Note how structures have been naturally completed in a coherent way. This example also shows the limitation of the reconstruction with regard to reconstruction of textured regions. For instance, the algorithm has not been able to reconstruct one eye of the woman on the left. This is very difficult task, and it is not

ANISOTROPIC DIFFUSION PARTIAL DIFFERENTIAL EQUATIONS

187

F IGURE 17. Denoising of the color image Lolotte, real color photograph. Noisy color image (left), denoised image (right); details are shown on the bottom row.

surprising that diffusion equations that perform only local smoothing fail in such complex cases. • Figure 21 shows how the PDE-based inpainting technique can be also used to remove real objects from digital photographs. A 500 × 500 color image (left) is inpainted with a user-defined mask (middle), corresponding to the region overlapping the initial man’s glasses. The inpainted image (right) is obtained in 4 minutes 11 seconds, after 200 iterations of the PDE (29) with parameters p1 = 0.001, p2 = 100, σ = 4, and dt = 150. Note that p1  p2 encourages smoothing only along the isophote directions with a strength of 1 everywhere. Another similar example of real object removal can be seen in Figure 24, where a cage automatically disappears in a color photograph of a parrot.

188

TSCHUMPERLÉ AND DERICHE

F IGURE 18. Removing blocs artifacts from the JPEG-compressed color image Baby. JPEG color image (left), regularized image (right); details are shown on the bottom row.

F IGURE 19. Removing quantization noise in the quantized color image Penguin. Quantized color image (left), regularized image (right) (details).

• Figures 22 and 23 show the reconstruction capabilities of our inpainting technique. Here, half of the pixels of two color images have been suppressed by masking them with a checkerboard-shaped mask with squares

ANISOTROPIC DIFFUSION PARTIAL DIFFERENTIAL EQUATIONS

189

F IGURE 20. Removing undesirable text in color image by an inpainting method. Corrupted color image (left), inpainted image (right). Details are shown on the bottom row.

of respective sizes 16 × 16 and 32 × 32. The images are then reconstructed using several iterations of the PDE (29) applied only inside the inpainting masks, with parameters p1 = 0.001, p2 = 100, σ = 4, and dt = 50. This application is very interesting for image transmission within a network. It can be used to generate coherent reconstructed images even if corrupted network packets have been received. For each inpainting result shown, the initialization of the pixel values inside the inpainting masks at t = 0 has been done by white noise. The inpainting algorithm does not depend much on of the initialization step; Eq. (29) diffuses neighborhood values inside the inpainting mask until convergence, and there is then a strong border condition. Different types of initialization (noise, zerofilling, or linear interpolation) did not produce much difference in results. C. Color Image Interpolation These same techniques can easily be used to perform image magnification by edge-preserving interpolation. Starting from a linear or bicubic interpolation

190

TSCHUMPERLÉ AND DERICHE

F IGURE 21. Removing a real object in a color image by an inpainting method. Original color image (left), inpainting mask (middle), inpainted image (right).

F IGURE 22. Reconstruction of 50% of a color image by an inpainting method. Original color image (left), masked color image (middle), reconstruction (right).

of a small image and applying the PDE (29) on the image (except on the original known pixels that form a sparse inpainting mask), allows computation of magnified images that have been regularized while taking their local geometry into account. It allows removal of usual block or jagging effects inherent to classical linear interpolation techniques. This technique is very similar to image inpainting with a very sparse grid for the inpainting mask. • Figure 25 shows one example of image resizing. An original 195 × 173 color image is resized by a factor of 4 with classical nearest-neighbor, linear, and bicubic interpolations, and then by the PDE-based technique (29). The nonlinear regularization filter, driven by the image geometry, clearly allows removal of the aliasing effects usually encountered with simple

ANISOTROPIC DIFFUSION PARTIAL DIFFERENTIAL EQUATIONS

191

F IGURE 23. Reconstruction of 50% of a color image by an inpainting method. Masked color image (left), reconstruction (right).

F IGURE 24. Removing a real object in a color image by an inpainting method. Original color image (left), inpainting mask (middle), inpainted image (right).

interpolation methods, while correctly preserving the small structures of the image. Note that the original known points of the images are always preserved during the regularizing PDE flow, ensuring that resizing the image back to its original dimension (subsampling) always results in the original input data. D. Flow Visualization A final application of regularization PDEs is presented here for visualization purposes. A 2D vector field F :  → R2 can be visualize in several ways. We can first use vector graphics (Figure 26) (left), but we need to subsample the field since this type of representation is not adapted to represent very dense flows. A better solution is as follows. We smooth a completely noisy (color) image I, with a regularizing flow equivalent to Eq. (29) but where T is directed by the directions of F , instead of the local geometry of I. The tensor field D

192

TSCHUMPERLÉ AND DERICHE

(a)

(b)

(c)

(d)

(e) F IGURE 25. Comparisons of image resizing methods. Original thumbnail image (a), nearest-neighbor (b), linear (c), bicubic (d) and PDE-based (e) interpolations.

ANISOTROPIC DIFFUSION PARTIAL DIFFERENTIAL EQUATIONS

193

F IGURE 26. Visualization of a 2D vector field using arrows (left), after 5 PDE iteration (middle), after 15 PDE iteration (right).

in Eq. (29) must be defined as ∀X ∈ ,

D(X) =

1 FF T . F 

(30)

This is a field of fully anisotropic tensors, each time oriented along the flow F . This technique is in fact equivalent to using a single LIC filter (Cabral and Leedom, 1993), but it is interesting to note that when the PDE evolution time t goes by, a visualization scale-space of F is explicitly constructed (Figure 27). Here, the regularization Eq. (30) used ensures that the smoothing of the pixels is done exactly in the direction of the flow F . This is not the case in (Becker et al., 2000; Buerkle et al., 2001; Diewald et al., 2000; Preusser and Rumpf, 1999), where the authors proposed a similar idea using a divergencebased expression. Using similar divergence-based techniques raises a risk of smoothing the image in false directions, as noted in Section II.B.

VI. C ONCLUSION Multichannel image regularization is a fundamental process in many image processing and computer vision algorithms. It is then necessary to gain full control of this process, as well as to understand how the equations function from a geometric point of view. This chapter has described the most classical PDE-based methods proposed for the regularization of multichannel images and introduced a very efficient curvature-preserving framework that generally outperforms its competitors. This is due not only to the particular aim of preserving fine and curved structures, but also thanks to the proposed numerical scheme that is especially efficient since it works at a subpixel level. Clearly, this kind of multichannel image regularization technique can play a role in many image-processing applications. The processing time, which was

194

TSCHUMPERLÉ AND DERICHE

F IGURE 27.

Visualization scale-space generated with regularization PDEs.

one of the main drawbacks of PDE-based methods, is no longer problematic. All these reasons make the framework of multivalued diffusion PDEs a very good choice for image regularization purposes. This chapter has shown the results on color image denoising, inpainting, and resizing, but many other applications may benefit from the proposed curvature-preserving framework. Other application results of the curvature-preserving algorithm can be found at the following web page: http://www.greyc.ensicaen.fr/~dtschump/ greycstoration/. The binaries of the algorithm can be also downloaded and tested on different architectures, as well as the source code (C++), which are available as a part of the open source image-processing library: The CImg Library.1 1 Tschumperlé, D., The CImg Library. Available at: http://cimg.sourceforge.net. The C++ Template Image Processing Library.

ANISOTROPIC DIFFUSION PARTIAL DIFFERENTIAL EQUATIONS

195

A PPENDIX A In this Appendix, we demonstrate how the minimization of the functional  (A.1) min n E(I) = ψ(λ+ , λ− ) d I : →R



can be performed by a gradient descent and its corresponding PDE flow. The Euler–Lagrange equations corresponding to the functional Eq. (A.1) are:  ∂ψ  ∂Iix ∂Ii (i = 1, . . . , n). (A.2) = div ∂ψ ∂t ∂Iiy

∂ψ ( ∂I ix

∂ψ T ∂Iiy )

Actually, the vector , can be written in a more comprehensive form. From the chain-rule property of the derivation, we have:  ∂ψ   ∂λ+ ∂λ−   ∂ψ  ∂Iix ∂ψ ∂Iiy

=

∂Iix ∂λ+ ∂Iiy

∂Iix ∂λ− ∂Iiy

∂λ+ ∂ψ ∂λ−

(A.3)

.

∂ψ We know formally the expressions ∂λ since the function ψ is directly defined ± from the λ± . ∂λ± ± Finding the ∂λ ∂Iix and ∂Iiy is more difficult. Here is a simple way to proceed. As the λ± are the eigenvalues of the structure tensor G = (gkl ), we may decompose its derivatives (with respect to Iix and Iiy ), in terms of derivatives with respect to the gkl : ∂λ±  ∂λ± ∂gkl ∂λ±  ∂λ± ∂gkl = and = . (A.4) ∂Iix ∂gkl ∂Iix ∂Iiy ∂gkl ∂Iiy k,l

k,l

∂gkl The expressions ∂I and ix ⎧ 11 ⎨ ∂g ∂Ii = 2Iix , x

⎩ ∂g11 = 0, ∂Ii y

∂gkl ∂Iiy

are particularly simple: ⎧ ⎧ 12 22 ⎨ ∂g ⎨ ∂g ∂Iix = Iiy , ∂Iix = 0, and and ⎩ ∂g12 = Iix , ⎩ ∂g22 = 2Iiy , ∂Ii ∂Ii y

that is, Eq. (A.4) can be written as:  ∂λ±   ∂Iix ∂λ± ∂Iiy

=

∂λ± 2 ∂g 11 ∂λ± ∂g12

y

∂λ± ∂g12 ∂λ± 2 ∂g 22

 ∇Ii .

(A.5)

± Thus, one last obstacle remains: finding the formal expressions of ∂λ ∂gkl . Remember that λ± and θ± are the eigenvalues and eigenvectors of the

196

TSCHUMPERLÉ AND DERICHE

structure tensor G: G = λ+ θ+ θ+T + λ− θ− θ−T . The derivation of this tensor, with respect to one of its coefficient gkl is: ∂θ+ T ∂θ− T ∂G ∂λ+ ∂λ− = θ+ θ+T + θ− θ−T + λ+ θ + + λ− θ ∂gkl ∂gkl ∂gkl ∂gkl ∂gkl − + λ+ θ+

∂θ+T ∂θ T + λ− θ− − . ∂gkl ∂gkl

(A.6)

Moreover, as the θ± are unitary and orthogonal eigenvectors, we have 

θ+T θ+ = θ−T θ− = 1, θ+T θ− = θ−T θ+ = 0,

and

⎧ T ⎨ ∂θ+ θ+ = θ+T ∂θ+ = 0, ∂gkl ∂gkl ⎩ ∂θ−T

T ∂θ− ∂gkl θ− = θ− ∂gkl = 0.

(A.7)

We first multiply Eq. (A.6) by θ±T at the left, by θ± at the right, and then use the properties of Eq. (A.7). It allows high simplifications and leads to these two relations: ∂G ∂λ+ = θ+T θ+ ∂gkl ∂gkl

and

∂G ∂λ− = θ−T θ− . ∂gkl ∂gkl

(A.8)

The relations in Eq. (A.8) formally show how eigenvalues of a diffusion tensor G vary with respect to a particular coefficient gkl of G. This interesting property can be proved for any symmetric matrix. For instance, authors of (Papadopoulo and Lourakis, 2000) proposed a similar demonstration in a purely matrix form, leading to the same result. They used it to deal with general covariance matrices. ∂G are very simple to write: Moreover, in our case the matrices ∂g kl ∂G = ∂g11



1 0

0 0

 ,

∂G = ∂g12



0 1

1 0

 ,

and

∂G = ∂g22



0 0 0 1

 .

With all these elements, we can express (A.5) as  ∂λ+  ∂Iix ∂λ+ ∂Iiy

 ∂λ−  = 2θ+ θ+T ∇Ii

and

∂Iix ∂λ− ∂Iiy

= 2θ− θ−T ∇Ii .

(A.9)

Finally, replacing Eq. (A.9) in the Euler–Lagrange Eqs. (A.3) and (A.2) gives the vector-valued gradient descent of the functional [Eq. (8)]:

ANISOTROPIC DIFFUSION PARTIAL DIFFERENTIAL EQUATIONS

197

 min

ψ(λ+ , λ− ) d

I : →Rn 

   ∂ψ ∂Ii ∂ψ T T = 2 div ⇒ θ+ θ + + θ− θ− ∇Ii ∂t ∂λ+ ∂λ− (for i = 1, . . . , n).

(A.10)

Note that Eq. (A.10) is a divergence-based equation such that ∂Ii = div(D∇Ii ), ∂t

where D = 2

∂ψ ∂ψ θ+ θ+T + 2 θ− θ−T ∂λ+ ∂λ−

D ∈ P(2) is then a 2 × 2 diffusion tensor, whose eigenvalues are λ1 = 2

∂ψ ∂λ+

and

λ2 = 2

∂ψ ∂λ−

associated with these corresponding orthonormal eigenvectors: u1 = θ +

and u2 = θ− .

Computing this gradient descent is done exactly in the same way when dealing with image domains  defined in higher-dimensional spaces ( ⊂ Rp where p > 2). More particularly, the case of 3D volume regularization (p = 3) can be written as  min n ψ(λ1 , λ2 , λ3 ) d I : →R



⇒

   ∂ψ ∂Ii ∂ψ ∂ψ T T T = 2 div θ1 θ + θ2 θ + θ3 θ ∇Ii . ∂t ∂λ1 1 ∂λ2 2 ∂λ3 3

In this case, the λ1,2,3 are the three eigenvalues of the 3 × 3 structure tensor G, and θ1,2,3 are the corresponding orthonormal eigenvectors.

A PPENDIX B This appendix demonstrates that the solution of the generic trace-based PDE ∀i = 1, . . . , n,

∂Ii = trace(THi ) ∂t

is the convolution of the image I Ii(t) = Ii(t=0) ∗ G(T,t)

(i = 1, . . . , n)

198

TSCHUMPERLÉ AND DERICHE

by an oriented Gaussian kernel G(T,t) defined as   1 XT T−1 X G(T,t) (X) = exp − with X = ( x 4π t 4t

y )T .

To demonstrate this, we simply derive the kernel expression G(T,t) in time and in space    XT T−1 X 1 XT T−1 X ∂G(T,t) exp − =− 1 − ∂t 4t 4t 4π t 2 and

⎧

 ⎨ ∇G(T,t) = − 1 2 exp − XT T−1 X T−1 X, 4t 8πt ⎩ H (T,t) = − 1 exp − XT T−1 X T−1 I − d G 4t 8πt 2

XXT T−1 2t



,

where ∇G(T,t) and HG(T,t) are, respectively, the gradient and the Hessian of G(T,t) . This means that     1 XT T−1 X XXT T−1 trace Id − trace(THG(T,t) ) = − exp − 4t 2t 8π t 2    T −1 T −1 1 X T X X T X =− 2− exp − 4t 2t 8π t 2 ∂G(T,t) . ∂t As the convolution is a linear operation, we have =

∂(Ii0 ∗ G(T,t) ) ∂G(T,t) = Ii0 ∗ = Ii0 ∗ trace(THG(T,t) ) ∂t ∂t = trace(THIi ∗G(T,t) ) 0

as well as

 lim Ii(t) ∗ G(T,t) = Ii0 ,

t→0

which tells us that the initial condition at t = 0 is coherent both for the PDE and the convolution process, since the Gaussian function G(T,t) is normalized. This statement is thus true for each instant t of the PDE flow.

A PPENDIX C This appendix develops tensor-driven divergence PDEs into their tracebased counterpart. Most divergence-based regularization PDEs acting on

ANISOTROPIC DIFFUSION PARTIAL DIFFERENTIAL EQUATIONS

199

multivalued images have the following form: ∂Ii = div(D∇Ii ) (i = 1, . . . , n), (C.1) ∂t where D is a diffusion tensor based only onfirst-order operators. D is often computed from the structure tensor G = nj=1 ∇Ij ∇IjT and depends mainly on the spatial derivatives Iix and Iiy . Intuitively, as the divergence ∂ ∂ + ∂y is itself a first-order derivative operator, we should be able div( ) = ∂x to write Eq. (C.1) only with first and second spatial derivatives Iix , Iiy , Iixx , Iixy , and Iiyy . Thus, it could be expressed with oriented Laplacians in each image channel Ii as well, that is an expression based on the trace operator ∂Ii ∂t = trace(DHi ). We want to make the link between the two different diffusion tensors D and T in the divergence-based and trace-based regularization PDEs, in the case when D is not constant: ∂Ii ∂Ii = div(D∇Ii ) and = trace(THi ). ∂t ∂t As noted in the previous section, these two formulations are almost equivalent, up to an additional term depending on the variation of the tensor field D: −−→

−−→

div(D∇Ii ) = trace(DHIi ) + ∇IiT div(D),

(C.2)

where div( ) is the matrix divergence. −−→ A natural progression is to then decompose the additional term ∇IiT div(D) into oriented Laplacians, expressed with additional diffusion tensors Q in the trace operator. For this purpose, we consider that the divergence tensor D is defined at each point X ∈  by D = f1 (λ+ , λ− )θ+ θ+T + f2 (λ+ , λ− )θ− θ−T

with f1/2 : R2 → R. (C.3)

It means that D is only expressed from the eigenvalues λ± and the eigenvectors θ± of the structure tensor G: G = λ+ θ+ θ+T + λ− θ− θ−T . This is indeed a very generic hypothesis that is verified by the majority of the proposed vector-valued regularization methods, for instance, the one proposed in Appendix A: ⎧ ⎨ f1 (λ+ , λ− ) = 2 ∂ψ , ∂Ii ∂λ+ = div(D∇Ii ) with Eq. (C.3) and ⎩ f2 (λ+ , λ− ) = 2 ∂ψ . ∂t ∂λ−

200

TSCHUMPERLÉ AND DERICHE −−→

In order to develop the additional diffusion term ∇IiT div(D) in Eq. (C.2), we propose to write D as a linear combination of G and Id : D = α(λ+ , λ− )G + β(λ+ , λ− )Id

(C.4)

that is, we separate the isotropic and anisotropic parts of D, with f1 (λ+ , λ− ) − f2 (λ+ , λ− ) and λ+ − λ − λ+ f2 (λ+ , λ− ) − λ− f1 (λ+ , λ− ) β= . λ+ − λ −

α=

(C.5)

Thus  f1 − f2 λ+ θ+ θ+T + λ− θ− θ−T λ+ − λ −  λ+ f2 − λ− f1 + θ+ θ+T + θ− θ−T λ+ − λ −  1 = θ+ θ+T (λ+ f1 − λ− f1 ) + θ− θ−T (λ+ f2 − λ− f2 ) λ+ − λ −

αG + βId =

= f1 θ+ θ+T + f2 θ− θ−T = D. Here we assumed that λ+ = λ− (i.e., the structure tensor G) is anisotropic. If G is isotropic, an isotropic diffusion tensor D is usually chosen too, in the divergence operator of Eq. (C.2), i.e., f1 (λ+ , λ− ) = f2 (λ+ , λ− ). In this case, we choose α = 0 and β = f1 (λ+ , λ− ). −−→ This decomposition is useful to rewrite the matrix divergence div(D) as: −−→

−−→

div(D) = α div(G) + G∇α + ∇β and the additional term of Eq. (C.2) would be computed as follows: 

−−→ −−→ ∇I T div(D) = trace div(D)∇IiT

−−→  = α trace div(G)∇IiT 

+ trace G∇α∇IiT

 + trace ∇β∇IiT .

(C.6)

(C.7) (C.8) (C.9)

In the following, we propose to find formal expressions of Eqs. (C.7)–(C.9). • First, remember that the structure tensor G is defined as: G=

n  j =1

∇Ij ∇IjT .

ANISOTROPIC DIFFUSION PARTIAL DIFFERENTIAL EQUATIONS

201

We have then: −−→

div(G) =

 2 n  Ijx −−→ div Ijx Ijy j =1

=

Ijx Ijy Ij2y



 n   2Ijx Ijxx + Ijx Ijyy + Ijy Ijxy Ijx Ijxy + Ijy Ijxx + 2Ijy Ijyy j =1

=

   n   Ijx (Ijxx + Ijyy ) I I + Ijy Ijxy + jx jxx Ijx Ijxy + Ijy Ijyy Ijy (Ijxx + Ijyy ) j =1

=

n 

Ij ∇Ij + Hj ∇Ij ,

j =1

where Ij and Hj are, respectively, the Laplacian and the Hessian of the image component Ij . Then, we can write Eq. (C.7) as n



−−→    α trace Hj ∇IiT ∇Ij Id + ∇Ij ∇IiT . α trace div(G)∇IiT = j =1

(C.10)

• We finally need to compute ∇α and ∇β in Eqs. (C.8) and (C.9). This can be done by the decomposition ∇α =

∂α ∂α ∇λ+ + ∇λ− ∂λ+ ∂λ−

and

∇β =

∂β ∂β ∇λ+ + ∇λ− ∂λ+ ∂λ− (C.11)

and as the λ± , eigenvalues of the structure tensor G, depend on the Ijx and Ijy :  λ±x λ± y ⎛ ∂λ ⎞ ∂λ± ± n  ∂Ijx Ijxx + ∂Ijy Ijxy ⎝ ⎠ = ∂λ± ∂λ± I + I j =1 ∂Ijx jxy ∂Ijy jyy   ∂λ ± n  ∂Ix HIj ∂λ±j . = 

∇λ± =

j =1

∂Iyj

202

TSCHUMPERLÉ AND DERICHE

In Appendix A, we derived eigenvalues of a structure tensor G, with respect to the spatial image derivatives. We derived the following relation:  ∂λ±  ∂Ixj ∂λ± ∂Iyj

= 2θ± θ±T ∇Ij .

Then, ∇λ± =

n  j =1

2Hj θ± θ±T ∇Ij .

(C.12)

We can replace Eq. (C.12) into the expressions of Eq. (C.11) to find the spatial gradients of α and β:

∂α    ∂α ∇α = nj=1 2Hj ∂λ θ+ θ+T + ∂λ θ+ θ+T ∇Ij , + − (C.13)

∂β   T + ∂β θ θ T ∇I . ∇β = nj=1 2Hj ∂λ θ θ + + j + + ∂λ− + Using Eq. (C.13), we finally compute−−→ the two missing parts of Eqs. (C.8) and (C.9) of the additional term ∇IiT div(D):

∂α     ∂α θ+ θ+T + ∂λ θ− θ−T ∇Ij ∇IiT , trace(G∇α∇IiT ) = nj=1 trace 2GHj ∂λ + −

∂β    ∂β θ+ θ+T + ∂λ θ− θ−T ∇Ij ∇IiT . trace(∇β∇IiT ) = nj=1 trace 2Hj ∂λ + − (C.14) • The final step consists in combining Eqs. (C.10) and (C.14) to express the −−→ additional term ∇IiT div(D) in the PDE (C.2). −−→

∇IiT div(D) =

n 



trace Hj Pij ,

(C.15)

j =1

where the Pij are the following 2 × 2 matrices: Pij = α∇IiT ∇Ij Id   ∂α ∂α T T +2 θ+ θ + + θ− θ− ∇Ij ∇IiT G ∂λ+ ∂λ−      ∂β ∂β T T θ+ θ + + α + θ− θ− ∇Ij ∇IiT . +2 α+ ∂λ+ ∂λ− (C.16) Note that the indices i, j in the notation Pij do not designate the coefficients of a matrix P, but the parameters of the family consisting of n2 matrices Pij (each of them is a 2 × 2 matrix).

ANISOTROPIC DIFFUSION PARTIAL DIFFERENTIAL EQUATIONS

203

The matrices Pii are symmetric, but generally not the Pij (where i = j ), since the gradients ∇Ii and ∇Ij are not usually aligned. However, we want to express Eq. (C.15) only with symmetric matrices, in order to interpret it as a sum of local smoothing processes oriented by diffusion tensors. Fortunately, the trace operator has this simple property:   A + AT trace(AH) = trace H , 2 where (A + AT )/2 is a 2 × 2 symmetric matrix (the symmetric part of A). Thus, we define the symmetric matrices Qij , corresponding to the symmetric parts of the Pij : Pij + Pij Q = 2

T

ij

(C.17)

and derive −−→

∇IiT div(D) =

n 

 trace Hj Qij .

j =1

Finally, the divergence-based PDE (C.2) can be written as: div(D∇Ii ) =

n 

 

trace δij D + Qij Hj ,

(C.18)

j =1

where δij is the Kronecker’s symbol:  0 if i =  j, δij = 1 if i = j. This completes the link between divergence PDEs and sums of atomic trace-based PDEs. A direct geometric interpretation of Eq. (C.18) is not direct.

R EFERENCES Alvarez, L., Guichard, F., Lions, P.L., Morel, J.M. (1993). Axioms and fundamental equations of image processing. Arch. Ration. Mech. Anal. 123 (3), 199–257. Alvarez, L., Lions, P.L., Morel, J.M. (1992). Image selective smoothing and edge detection by nonlinear diffusion (II). SIAM J. Numer. Anal. 29, 845– 866.

204

TSCHUMPERLÉ AND DERICHE

Alvarez, L., Mazorra, L. (1994). Signal and image restoration using shock filters and anisotropic diffusion. SIAM J. Numer. Anal. 31 (2), 590–605. Aubert, G., Kornprobst, P. (2002). Mathematical Problems in Image Processing: Partial Differential Equations and the Calculus of Variations. In: Applied Mathematical Sciences, vol. 147. Springer-Verlag. Ashikhmin, M. (2001). Synthesizing natural textures. In: ACM Symposium on Interactive 3D Graphics. Research Triangle Park, NC, pp. 217–226. Barash, D. (2000). Bilateral filtering and anisotropic diffusion: Towards a unified viewpoint. Technical Report, HP Laboratories Israel. Becker, J., Preusser, T., Rumpf, M. (2000). PDE methods in flow simulation post processing. Comput. Vis. Sci. 3 (3), 159–167. Bertalmio, M., Cheng, L.T., Osher, S., Sapiro, G. (2000a). Variational problems and partial differential equations on implicit surfaces: The framework and examples in image processing and pattern formation. UCLA Research Report. Bertalmio, M., Sapiro, G., Caselles, V., Ballester, C. (2000b). Image inpainting. In: Proceedings of the SIGGRAPH. ACM Press, Addison–Wesley– Longman, pp. 417–424. Bertalmio, M., Cheng, L.T., Osher, S., Sapiro, G. (2001). Variational problems and partial differential equations on implicit surfaces. Comput. Vis. Sci. 174 (2), 759–780. Bertalmio, M., Vese, L., Sapiro, G., Osher, S. (2003). Simultaneous structure and texture image inpainting. IEEE Trans. Image Process. 12 (8), 882–889. Black, M.J., Sapiro, G., Marimont, D.H., Heeger, D. (1998). Robust anisotropic diffusion. IEEE Trans. Image Process. 7 (3), 421–432. Blanc-Feraud, L., Charbonnier, P., Aubert, G., Barlaud, M. (1995). Nonlinear image processing: Modeling and fast algorithm for regularization with edge detection. In: Proceedings of the International Conference on Image Processing (ICIP), Washington, DC, 1995, pp. 474–477. Blomgren, P. (1998). Total Variation Methods for Restoration of Vector Valued Images. PhD dissertation, Department of Mathematics, University of California, Los Angeles. Blomgren, P., Chan, T.F. (1998). Color TV: Total variation methods for restoration of vector-valued images. IEEE Trans. Image Process. 7 (3), 304–309. Buerkle, D., Preusser, T., Rumpf, M. (2001). Transport and diffusion in timedependent flow visualization. In: Proceedings IEEE Visualization. Cabral, B., Leedom, L.C. (1993). Imaging vector fields using line integral convolution SIGGRAPH’93. Computer Graphics 27, 263–272. Carmona, R., Zhong, S. (1998). Adaptive smoothing respecting feature directions. IEEE Trans. Image Process. 7 (3), 353–358. Chambolle, A., Lions, P.L. (1997). Image recovery via total variation minimization and related problems. Numer. Math. 76 (2), 167–188.

ANISOTROPIC DIFFUSION PARTIAL DIFFERENTIAL EQUATIONS

205

Chan, T., Kang, S.H., Shen, J. (2002). Euler’s elastica and curvature based inpainting. SIAM J. Appl. Math. Chan, T.F., Kang, S.H., Shen, J. (2000). Total variation denoising and enhancement color images based on the CB and HSV color models. J. Vis. Comm. Image Represent. 12 (4). Chan, T.F., Shen, J. (2001). Non-texture inpainting by curvature-driven diffusions (CDD). J. Vis. Comm. Image Represent. 12 (4), 436–449. Chan, T., Shen, J. (2000a). Variational restoration of non-flat image features: Models and algorithms. SIAM J. Appl. Math. 61 (4), 1338–1361. Chan, T., Shen, J. (2000b). Mathematical models for local deterministic inpaintings. Technical Report 00-11. Department of Mathematics, UCLA, Los Angeles. Charbonnier, P., Aubert, G., Blanc-Féraud, M., Barlaud, M. (1994). Two deterministic half-quadratic regularization algorithms for computed imaging. In: Proceedings of the International Conference on Image Processing (ICIP), vol. II, pp. 168–172. Charbonnier, P., Blanc-Féraud, L., Aubert, G., Barlaud, M. (1997). Deterministic edge-preserving regularization in computed imaging. IEEE Trans. Image Process. 6 (2), 298–311. Chefd’hotel, C., Tschumperlé, D., Deriche, R., Faugeras, O. (2004). Regularizing flows for constrained matrix-valued images. J. Math. Imaging Vision 20 (2), 147–162. Chu, M.T. (1990a). A list of matrix flows with applications. Technical Report. Department of Mathematics, North Carolina State University. Chu, M.T. (1990b). Matrix differential equations: A continuous realization process for linear algebra problems. Technical Report. Department of Mathematics, North Carolina State University. Coulon, O., Alexander, D.C., Arridge, S.R. (2001). A regularization scheme for diffusion tensor magnetic resonance images. In: 17th International Conference on Information Processing in Medical Imaging, Lecture Notes in Computer Science, vol. 2082. Springer, pp. 92–105. Criminisi, A., Perez, P., Toyama, K. (2003). Object removal by exemplarbased inpainting. In: IEEE Conference on Computer Vision and Pattern Recognition (CVPR), vol. 2, pp. 721–728. Deriche, R., Faugeras, O. (1997). Les EDP en traitement des images et vision par ordinateur. Traitement Signal 13 (6). Diewald, U., Preusser, T., Rumpf, M. (2000). Anisotropic diffusion in vector field visualization on Euclidian domains and surfaces. IEEE Trans. Vis. Comput. Graphics 6 (2), 139–149. Di Zenzo, S. (1986). A note on the gradient of a multi-image. Comput. Vision, Graphics Image Process. 33, 116–125.

206

TSCHUMPERLÉ AND DERICHE

Gilboa, G., Sochen, N., Zeevi, Y. (2002). Forward-and-backward diffusion processes for adaptative image enhancement and denoising. IEEE Trans. Image Process. Hadamard, J. (1923). Lectures on the Cauchy Problem in Linear Partial Differential Equations. Yale University Press, New Haven. Jia, J., Tang, C.K. (2003). Image repairing: Robust image synthesis by adaptive ND tensor voting. In: IEEE Conference on Computer Vision and Pattern Recognition (CVPR), vol. 1, pp. 643–650. Kichenassamy, S. (1997). The Perona–Malik paradox. SIAM J. Appl. Math. 57 (5), 1328–1342. Kimmel, R., Malladi, R., Sochen, N. (1998). Image processing via the Beltrami operator. In: Proceedings of the 3rd Asian Conference on Computer Vision, vol. 1. Hong Kong, pp. 574–581. Kimmel, R., Malladi, R., Sochen, N. (2000). Images as embedded maps and minimal surfaces: Movies, color, texture, and volumetric medical images. Int. J. Comput. Vision 39 (2), 111–129. Kimmel, R., Sochen, N. (1999). Geometric-variational approach for color image enhancement and segmentation. In: Scale-Space Theories in Computer Vision, Scale-Space’99, Lecture Notes in Computer Science, vol. 1682. Springer, pp. 295–305. Kimmel, R., Sochen, N. (2002). Orientation diffusion or how to comb a porcupine. J. Vis. Commun. Image Represent. 13, 238–248. Koenderink, J.J. (1984). The structure of images. Biol. Cybernet. 50, 363–370. Kornprobst, P. (1998). Contributions à la Restauration d’Images et à l’Analyse de Séquences: Approches Variationnelles et Solutions de Viscosité. PhD thesis, Université de Nice-Sophia Antipolis. Kornprobst, P., Deriche, R., Aubert, G. (1996). Image restoration via PDEs. In: First Annual Symposium on Enabling Technologies for Law Enforcement and Security—SPIE Conference 2942. Boston, Massachusetts. Kornprobst, P., Deriche, R., Aubert, G. (1997). Nonlinear operators in image restoration. In: Proceedings of the IEEE International Conference on Computer Vision and Pattern Recognition. Puerto Rico, pp. 325–331. Koschan, A. (1995). A comparative study on color edge detection. In: Proceedings of the 2nd Asian Conference on Computer Vision, ACCV’95, pp. 574–578. Krissian, K. (2000). Multiscale analysis: Application to medical imaging and 3D vessel detection. PhD Thesis. INRIA, Sophia Antipolis. Lindeberg, T. (1994). Scale-Space Theory in Computer Vision. Kluwer Academic Publishers. Masnou, S., Morel, J.-M. (1998). Level lines based disocclusion. IEEE Int. Conf. Image Process. 3, 259–263.

ANISOTROPIC DIFFUSION PARTIAL DIFFERENTIAL EQUATIONS

207

Meyer, Y. (2001). Oscillatory Patterns in Image Processing and Nonlinear Evolution Equations. University Lecture Series, vol. 22. American Mathematical Society, Providence, RI. Nielsen, M., Florack, L., Deriche, R. (1997). Regularization, scale-space and edge detection filters. J. Math. Imaging Vision 7 (4), 291–308. Nikolova, M., (2001). Image restoration by minimizing objective functions with nonsmooth data-fidelity terms. In: IEEE Workshop on Variational and Level Set Methods. Vancouver, Canada, pp. 11–19. Nikolova, M., Ng, M. (2001). Fast image reconstruction algorithms combining half-quadratic regularization and preconditioning. In: Proceedings of the International Conference on Image Processing. IEEE Signal Processing Society. Papadopoulo, T., Lourakis, M.I.A. (2000). Estimating the Jacobian of the Singular Value Decomposition: Theory and Applications. Research Report 3961. INRIA, Sophia Antipolis. Pardo, A., Sapiro, G. (2000). Vector probability diffusion. In: Proceedings of the International Conference on Image Processing, IEEE Signal Processing Society. Perona, P. (1998). Orientation diffusions. IEEE Trans. Image Process. 7 (3), 457–467. Perona, P., Malik, J. (1990). Scale-space and edge detection using anisotropic diffusion. IEEE Trans. Pattern Anal. Machine Intell. 12 (7), 629–639. Poynton, C.A. (1995). Poynton’s colour FAQ, www.inforamp.net/poynton [Web page]. Preusser, T., Rumpf, M. (1999). Anisotropic nonlinear diffusion in flow visualization. In: IEEE Visualization Conference. Press, W.H., Flannery, B.P., Teukolsky, S.A., Vetterling, W.T. (1992). Runge– Kutta method. In: Numerical Recipes in FORTRAN: The Art of Scientific Computing. Cambridge University Press, pp. 704–716. ter Haar Romeny, B.M. (1994). Geometry-driven diffusion in computer vision. In: Computational Imaging and Vision. Kluwer Academic Publishers. Rudin, L., Osher, S., Fatemi, E. (1992). Nonlinear total variation based noise removal algorithms. Physica D 60, 259–268. Sapiro, G. (1996). Vector-valued active contours. In: Proceedings of the International Conference on Computer Vision and Pattern Recognition. San Francisco, pp. 680–685. Sapiro, G. (1997). Color snakes. Comput. Vision Image Understanding 68 (2). Sapiro, G. (2001). Geometric Partial Differential Equations and Image Analysis. Cambridge University Press. Sapiro, G., Ringach, D.L. (1996). Anisotropic diffusion of multivalued images with applications to color filtering. IEEE Trans. Image Process. 5 (11), 1582–1585.

208

TSCHUMPERLÉ AND DERICHE

Sochen, N., Kimmel, R., Bruckstein, A.M. (2001). Diffusions and confusions in signal and image processing. J. Math. Imaging Vision 14 (3), 195–209. Sochen, N. (2001). On affine invariance in the Beltrami framework for vision. In: IEEE Workshop on Variational and Level Set Methods. Vancouver, Canada, pp. 51–56. Sochen, N., Kimmel, R., Malladi, R. (1998). A geometrical framework for low level vision. IEEE Trans. Image Process. [Special Issue on PDE based Image Processing] 7 (3), 310–318. Stalling, D., Hege, H.C. (1995). Fast and resolution independent line integral convolution. In: ACM SIGGRAPH, 22nd Annual Conference on Computer Graphics and Interactive Technique, pp. 249–256. Tang, B., Sapiro, G., Caselles, V. (1998). Direction diffusion. In: International Conference on Computer Vision. Tang, B., Sapiro, G., Caselles, V. (2000). Diffusion of general data on non-flat manifolds via harmonic maps theory: The direction diffusion case. Int. J. Comput. Vision 36 (2), 149–161. Teboul, S., Blanc-Féraud, L., Aubert, G., Barlaud, M. (1998). Variational approach for edge-preserving regularization using coupled PDEs. IEEE Trans. Image Process. [Special Issue on PDE based Image Processing] 7 (3), 387–397. Tikhonov, A.N. (1963). Regularization of incorrectly posed problems. Sov. Math. Dokl. 4, 1624–1627. Tomasi, C., Manduchi, R. (1998). Bilateral filtering for gray and color images. In: Proceedings of the IEEE International Conference on Computer Vision, pp. 839–846. Tschumperlé, D. (2002). PDE’s Based Regularization of Multivalued Images and Applications. PhD thesis. Université de Nice Sophia Antipolis. Tschumperlé, D. (2006). Fast anisotropic smoothing of multi-valued images using curvature-preserving PDEs. Int. J. Comput. Vision 68 (1), 65–82. Tschumperlé, D., Deriche, R. (2001a). Constrained and unconstrained PDE’s for vector image restoration. In: Proceedings of the 10th Scandinavian Conference on Image Analysis, Bergen, Norway, pp. 153–160. Tschumperlé, D., Deriche, R. (2001b). Diffusion tensor regularization with constraints preservation. In: IEEE Computer Society Conference on Computer Vision and Pattern Recognition. Kauai, Hawaii. Tschumperlé, D., Deriche, R. (2001c). Regularization of orthonormal vector sets using coupled PDEs. In: IEEE Workshop on Variational and Level Set Methods. Vancouver, Canada, pp. 3–10. Tschumperlé, D., Deriche, R. (2002a). Diffusion PDE’s on vector-valued images: Local approach and geometric viewpoint. IEEE Signal Process. Mag. 19 (5), 16–25. Tschumperlé, D., Deriche, R. (2002b). Orthonormal vector sets regularization with PDE’s and aplications. Int. J. Comput. Vision.

ANISOTROPIC DIFFUSION PARTIAL DIFFERENTIAL EQUATIONS

209

Tschumperlé, D., Deriche, R. (2003). Vector-valued image regularization with PDEs: A common framework for different applications. In: IEEE Conference on Computer Vision and Pattern Recognition, vol. 1, pp. 651– 656. Tschumperlé, D., Deriche, R. (2005). Vector-valued image regularization with PDEs: A common framework for different applications. IEEE Trans. Pattern Anal. Machine Intell. 27 (4). Vese, L.A., Osher, S. (2001). Numerical methods for p-harmonic flows and applications to image processing. CAM Report 01-22. UCLA. Weickert, J. (1994), Anisotropic diffusion filters for image processing based quality control. In: 7th European Conference on Mathematics in Industry, pp. 355–362. Weickert J. (1996a). Anisotropic diffusion in image processing. PhD thesis. University of Kaiserslautern, Laboratory of Technomathematics. Germany. Weickert, J. (1996b). Theoretical foundations of anisotropic diffusion in image processing. Comput. Suppl. 11, 221–236. Weickert, J. (1997a). Coherence-enhancing diffusion of colour images. In: 7th National Symposium on Pattern Recognition and Image Analysis. Weickert, J. (1997b). A Review of Nonlinear Diffusion Filtering. In: ScaleSpace Theory in Computer Vision, Lecture Notes in Computer Science, vol. 1252. Springer, Berlin, pp. 3–28. Weickert, J. (1998). Anisotropic Diffusion in Image Processing. TeubnerVerlag, Stuttgart. Weickert, J. (1999). Coherence-enhancing diffusion of colour images. Image Vision Comput. 17, 199–210. Weickert, J., Benhamouda, B. (1997a). A semidiscrete nonlinear scale-space theory and its relation to the Perona–Malik paradox. In: Adv. Comput. Vision. Springer, Wien, pp. 1–10. Weickert, J., Benhamouda, B. (1997). Why the Perona–Malik Filter Works. Technical Report 97/22. Department of Computer Science, University of Copenhagen. Weickert, J., Brox, T. (2002). Diffusion and regularization of vector and matrix-valued images. Inverse Problems, Image Analysis, and Medical Imaging, Contemp. Math. 313, 251–268. Wei, L.Y., Levoy, M. (2000). Fast texture synthesis using tree-structured vector quantization. In: ACM SIGGRAPH, International Conference on Computer Graphics and Interactive Techniques, pp. 479–488. Wesseling, P. (2000). Principles of Computational Fluid Dynamics. Springer, Berlin. Witkin, A.P. (1983). Scale-space filtering. In: International Joint Conference on Artificial Intelligence, pp. 1019–1021. Yezzi, A. (1998). Modified curvature motion for image smoothing and enhancement. IEEE Trans. Image Process. 7 (3), 345–352.