David Tschumperlé

Rachid Deriche

Anisotropic Diffusion Partial Differential Equations in Multi-Channel Image Processing : Framework and Applications

Keywords: Multi-Channel Images Regularization, Anisotropic Smoothing, Diffusion PDE's, Tensor-Valued Geometry, Denoising, Inpainting, Nonlinear Interpolation. Preliminary

published or not. The documents may come

D = a b b c = λ 1 u 1 u T 1 + λ 2 u 2 u T 2 λ 1u1 λ 2 u 2 o
When λ 2 >> λ 1 (lenghtened ellipse), the tensor D is said to be anisotropic and has u 2 has its principal orientation. When λ 1 = λ 2 = β, the tensor D is isotropic and thus equal to a weighted version of the 2x2 identity matrix I d

Introduction

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 pre-processing 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 is very interesting as well, when considering the analysis of the images at multiple scales. In a more general manner, image regularization is one of the key stage of most computer vision algorithms since it plays a fundamental role for solving ill-posed computer vision problems [START_REF] Hadamard | Lectures on the Cauchy Problem in Linear Partial Differential Equations[END_REF], including restoration, segmentation, registration, surface reconstruction, etc. This explains why a lot of image regularization formalisms have been already proposed and studied in the literature. Perona & Malik in their pioneering work [START_REF] Perona | Scale-space and edge detection using anisotropic diffusion[END_REF] in the early 90's were the first to imagine image regularization in terms of anisotropic diffusion PDE's (Partial Differential Equations). 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 has been proved later to be unstable [START_REF] Weickert | Why the Perona-Malik filter works[END_REF]. Firstly created to describe physical laws and natural motions of mechanic objects and fluids (strings, water, wind [START_REF] Wesseling | Principles of Computational Fluid Dynamics[END_REF]), diffusion PDE's had been already widely studied and interesting theoretical results coming from the fields of physics and mathematics have found interesting implications for the purpose of data regularization. Actually, PDE's are local formulations and thus, they are well adapted to deal with degraded images where sources of data corruption are local or semi-local too. 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 & Malik, many authors have proposed variants of diffusion PDE's for image regularization since then, mostly 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 [START_REF] Krissian | Multiscale Analysis : Application to Medical Imaging and 3D Vessel Detection[END_REF][START_REF] Perona | Scale-space and edge detection using anisotropic diffusion[END_REF][START_REF] Sapiro | Geometric Partial Differential Equations and Image Analysis[END_REF][START_REF] Weickert | Anisotropic Diffusion in Image Processing[END_REF], how diffusion PDE's may be seen as gradient descents of various energy functionals [START_REF] Aubert | Mathematical Problems in Image Processing: Partial Differential Equations and the Calculus of Variations[END_REF][START_REF] Chambolle | Image recovery via total variation minimization and related problems[END_REF][START_REF] Charbonnier | Deterministic edge-preserving regularization in computed imaging[END_REF][START_REF] Kimmel | Images as embedded maps and minimal surfaces: movies, color, texture, and volumetric medical images[END_REF][START_REF] Rudin | Nonlinear total variation based noise removal algorithms[END_REF], and the link between regularization PDE's and the concept of nonlinear scale spaces [START_REF] Alvarez | Axioms and fundamental equations of image processing[END_REF][START_REF] Lindeberg | Scale-Space Theory in Computer Vision[END_REF][START_REF] Nielsen | Regularization, scale-space and edge detection filters[END_REF]. Extensions of these techniques to deal with color images and more generally multi-channel datasets have been more recently tackled in [START_REF] Blomgren | Color tv: Total variation methods for restoration of vector-valued images[END_REF][START_REF] Kimmel | Images as embedded maps and minimal surfaces: movies, color, texture, and volumetric medical images[END_REF][START_REF] Pardo | Vector probability diffusion[END_REF][START_REF] Sapiro | Geometric Partial Differential Equations and Image Analysis[END_REF][START_REF] Sapiro | Anisotropic diffusion of multivalued images with applications to color filtering[END_REF][START_REF] Tschumperlé | PDE's Based Regularization of Multivalued Images and Applications[END_REF][START_REF] Tschumperlé | Deriche Vector-Valued Image Regularization with PDE's : A Common Framework for Different Applications[END_REF][START_REF] Tschumperlé | Deriche Vector-Valued Image Regularization with PDE's : A Common Framework for Different Applications[END_REF][START_REF] Weickert | Anisotropic Diffusion in Image Processing[END_REF][START_REF] Weickert | Coherence-Enhancing Diffusion of Colour Images[END_REF] (among others), leading to more elaborated expressions : a coupling term between image channels generally appears in the equations. Diffusion equations dealing with constrained multi-dimensional datasets have been also proposed, allowing to regularize images of unit vectors [START_REF] Coulon | A regularization scheme for diffusion tensor magnetic resonance images[END_REF][START_REF] Kimmel | Orientation diffusion or how to comb a porcupine[END_REF][START_REF] Perona | Orientation diffusions[END_REF][START_REF] Tang | Direction diffusion. International Conference on Computer Vision (ICCV)[END_REF], orthonormal matrices [START_REF] Chefd'hotel | Regularizing Flows for Constrained Matrix-Valued Images[END_REF][START_REF] Tschumperlé | Regularization of orthonormal vector sets using coupled PDE's[END_REF], positive-definite matrices [START_REF] Chefd'hotel | Regularizing Flows for Constrained Matrix-Valued Images[END_REF][START_REF] Tschumperlé | Diffusion tensor regularization with constraints preservation[END_REF], or image data defined on implicit surfaces [START_REF] Bertalmio | Variational Problems and Partial Differential Equations on Implicit Surfaces[END_REF][START_REF] Chan | Variational restoration of non-flat image features : Models and algorithms[END_REF][START_REF] Tang | Diffusion of general data on non-flat manifolds via harmonic maps theory : The direction diffusion case[END_REF]. Usually, this kind of constrained PDE's simply get an extra constraint term added to the corresponding unconstrained equation and will not be discussed here.

Despite this wide range of existing constrained and unconstrained PDE formalisms for scalar and multi-channel images, all proposed methods have something in common : a nonlinear regularization PDE such as ∂I ∂t = R locally smoothes 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 scale-space : as the image data are gently regularized step-by-step, a continuous sequence of smoother images I (t) is generated whereas 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 [START_REF] Alvarez | Axioms and fundamental equations of image processing[END_REF][START_REF] Lindeberg | Scale-Space Theory in Computer Vision[END_REF][START_REF] Nielsen | Regularization, scale-space and edge detection filters[END_REF][START_REF] Perona | Scale-space and edge detection using anisotropic diffusion[END_REF][START_REF] Witkin | Scale-space filtering[END_REF]. Roughly speaking, regularization PDE's may be seen as iterative and nonlinear filters that simplify the image little by little and minimize then the image variations (Fig. 1). Note therefore that such equations generally do not converge towards a very interesting solution. Basically, the image obtained at convergence (t → ∞) is constant everywhere, corresponding to an image without any variations. This is indeed the most simplified image we can obtain. To avoid this undesired over-simplification, regularization algorithms are usually based on a modified PDE velocity R ′ = R+α (I noisy -I) including a so-called data fidelity term weighted by a user-defined parameter α ∈ R + . It avoids the expected solution (regularized image) at convergence to be 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 becomes thus a parameter of the method). Here, we are mainly interested in the regularization term R itself rather than the one containing the fidelity term R ′ . For a vast mathematical study about linear or nonlinear fidelity terms, please refer to [START_REF] Meyer | Oscillatory patterns in image processing and nonlinear evolution equations[END_REF][START_REF] Nikolova | Image restoration by minimizing objective functions with nonsmooth data-fidelity terms[END_REF][START_REF] Nikolova | Fast image reconstruction algorithms combining half-quadratic regularization and preconditioning[END_REF].

As it is clear that local and oriented image smoothing is one of the key idea used by most PDE-based regularization methods, it naturally leads to the problem of defining a coherent geometry from a multi-channel image. It must be the first aim of a good regularization algorithm. Following this simple and general principle, recent contributions [START_REF] Tschumperlé | Deriche Vector-Valued Image Regularization with PDE's : A Common Framework for Different Applications[END_REF][START_REF] Tschumperlé | Deriche Vector-Valued Image Regularization with PDE's : A Common Framework for Different Applications[END_REF][START_REF] Weickert | Anisotropic Diffusion in Image Processing[END_REF] 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 one hand, they unify a lot of previously proposed equations into generic diffusion PDE's 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, one retrieves the geometry of the structures inside the image (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 regularized enough.

In this chapter, we will first discuss the definition of a local geometry for multi-channel images, by reviewing and comparing proposed solutions in the literature [START_REF] Blomgren | Total Variation Methods for Restoration of Vector Valued Images[END_REF][START_REF] Blomgren | Color tv: Total variation methods for restoration of vector-valued images[END_REF][START_REF] Zenzo | A note on the gradient of a multi-image[END_REF][START_REF] Sapiro | Vector-Valued Active Contours[END_REF][START_REF] Weickert | Anisotropic Diffusion in Image Processing[END_REF] (section 1). Then, we will review important works already proposed for scalar and multi-channel image regularization within a diffusion PDE framework. These methods may be classified into three different approaches which are (1) variational formulations, (2) divergence expressions and (3) oriented Laplacians. We will mainly try to focus on the interpretation of the algorithms in terms of local smoothing (section 2). We particularly point out 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 multi-channel images while taking specific curvature constraints into account (section 3). This formulation is mathematically positioned between previous existing equations, in a way that it solves most issues encountered with classical regularization methods. Moreover, we show that a theoretical interpretation of the curvatureconstrained formalism exists in terms of Line Integral Convolutions which is a simple filtering technique originally proposed by Cabral and Leedom in [START_REF] Cabral | Imaging vector fields using line integral convolution[END_REF]. 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 4). This iterative scheme has two main advantages compared to classical PDE implementations : on one hand, it preserves thin image structures remarkably well, since it naturally works at a sub-pixel 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 PDE's a generic and very efficient approach for solving image processing problems needing multi-channel image regularization. We finally illustrate this effectiveness, in terms of computational speed and visual quality, with results on color image restoration, color image inpainting and non-linear resizing, among all possible applications in the area of image regularization (section 5).

1 Defining a Local Geometry for Multi-Channel Images

Local geometric features

As stated in the introduction, image regularization may be seen as a filter that reduces local pixel variations. More precisely, one wants to smooth a multi-channel image I : Ω → R n while preserving its edges (discontinuities in the image intensities), i.e. performs a local smoothing mostly along directions of the edges, avoiding a smoothing orthogonal to these edges. At a first glance, a naive idea would be to apply a scalar-valued regularization filter on each channel I i of the multi-channel image I, doing this independently for each i = 1...n. But, the correlation between image channels would be ignored in this case, and it might cause important disparities in the smoothing behavior, since local smoothing directions and amplitudes could be very different from each channel to another. Such decoupled regularization methods generally lead to undesirable over-smoothing effects destroying significant edge structures in the image. Multi-channel image regularization is rather based on a coherent image smoothing which locally uses the same smoothing directions and amplitudes for all image channels I i . Naturally, this means that one has first to measure the local geometry of a multi-channel image I. Such a geometry consists actually in the definition of these important features at each image point

X = (x, y) ∈ Ω of I : -Two orthogonal directions θ + (X) , θ - (X) ∈ S 1 (unit vectors of R 2 )
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 multi-channel images.

-A corresponding variation norm N (X) measuring the local strength of an edge. This is the extension of the vector gradient norm for multi-channel images.

In order to construct such a vector geometry, different approaches have been considered so far and are detailed below.

Geometry from a scalar feature

One simple method consists in computing first a scalar image f (I), using a vector to scalar function f : R n → 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 [START_REF] Poynton | Poynton's colour FAQ[END_REF] :

f = L * = 116 g(Y ) -16 with Y = 0.2125R + 0.7154G + 0.0721B
where g : R → R is defined by 16 116 else Thus, we may define a vector-valued local vector geometry {N , θ + , θ -} of I by choosing

g(s) = 3 √ s if s > 0.008856 g(s) = 7.787s +
       θ + = ∇f (I) ∇f (I) θ -⊥ θ + and N = ∇f (I)
However, this method has two major drawbacks. First, this is not always possible to easily define a significant function f for multi-channel 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 will not be able to detect iso-lightness vector contours in a color image. It is the case for the image illustrated on Fig. 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 here will be either isotropic or oriented in a wrong direction : the existing color edges inside the yin-yang symbol will be probably blurred.

Di Zenzo multi-valued geometry

In order to overcome this limitation, a very elegant solution has been proposed by Di Zenzo in [START_REF] Zenzo | A note on the gradient of a multi-image[END_REF]. He considers a multichannel image I : Ω → R n as a vector field, and looks for the local variations of the vector norm dI 2 , mainly given by a variation matrix G = (g i,j ). We get :

dI = I x dx + I y dy
where

I x = ∂I ∂x and I y = ∂I ∂y (∈ R n ) then dI 2 = dI T dI = I x 2 dx 2 + 2 I T x I y dxdy + I y 2 dy 2
i.e.

dI 2 = dX T G dX where G = n i=1 ∇I i ∇I T i and dX = dx dy
G is denoted as the structure tensor. It sums variation contributions from each image channel I i . It is easy to see that G is a 2 × 2 symmetric and semi positive-definite matrix. Its coefficients (g i,j ) are simply :

           g 11 = n i=1 I 2 ix g 12 = g 21 = n i=1 I ix I iy g 22 = n i=1 I 2 iy
In the common case of color images I = (R, G, B), G is defined as :

G = R 2 x + G 2 x + B 2 x R x R y + G x G y + B x B y R x R y + G x G y + B x B y R 2 y + G 2 y + B 2 y ( 1
)
The interesting point about G is that its positive eigenvalues λ +/-give the maximum and the minimum values of dI 2 while the orthogonal eigenvectors θ + and θ -are the corresponding orientations of these extrema, and are formally given by :

λ +/-= g 11 + g 22 ± √ ∆ 2 and θ +/-/ / 2 g 12 g 22 -g 11 ± √ ∆ (2)
where ∆ = (g 11g 22 ) 2 + 4 g 2 12 . The vectors θ ± are normalized to the unit 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 in [START_REF] Weickert | Anisotropic Diffusion in Image Processing[END_REF]. He rather proposed to study the eigenvalues and eigenvectors of a Gaussian-smoothed version G σ of the structure tensor G :

G σ = n i=1 ∇I iα ∇I T iα * G σ where ∇I iα = ∇(I i * G α ) (3)
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 worth to notice 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 (it is the case for the inside of the strips in Fig. 3a). For this configuration, the variation norm N we have to define should be low.

-When λ + ≫ λ -, there are a lot of vector variations. The current point may be located on a vector edge (it is the case for the edges of the strips in Fig. 3a). For this configuration, the variation norm N should be high.

-When λ + ≃ λ -≫ 0, we are located on a saddle point of the vector surface, which can be possibly a corner structure in the image (for instance, the intersections of the strips in Fig. 3a). In this case N should be even higher than for the previous configuration. Regularization algorithms have indeed a tendency to smooth corners fastly. A very high variation measure estimated on corner points would attenuate the smoothing there, which is often a desired effect.

Actually, a lot of proposed regularization algorithms acting on multi-channel images have implicitly or explicitly based their smoothing behavior from these Di Zenzo's 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 [START_REF] Blomgren | Total Variation Methods for Restoration of Vector Valued Images[END_REF][START_REF] Sapiro | Vector-Valued Active Contours[END_REF][START_REF] Sapiro | Color snakes[END_REF] (Fig. 3b and Fig. 4b). This norm will not particularly give importance to corners compared to straight edges.

-N -= λ +λ -, also called coherence norm, have been chosen in [START_REF] Sapiro | Anisotropic diffusion of multivalued images with applications to color filtering[END_REF][START_REF] Weickert | Anisotropic Diffusion in Image Processing[END_REF][START_REF] Weickert | Coherence-enhancing diffusion of colour images[END_REF]. 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 (Fig. 3c), as well as in the center and left-right parts of the child's eye (Fig. 4c). This will mainly perturb any regularization process that uses this norm since some colored sharp corners, considered as homogeneous regions, will be probably over-smoothed.

-N + = λ + + λ -, also denoted by ∇I is often chosen [START_REF] Bertalmio | Variational problems and partial differential equations on implicit surfaces: The framework and examples in image processing and pattern formation[END_REF][START_REF] Blomgren | Color tv: Total variation methods for restoration of vector-valued images[END_REF][START_REF] Pardo | Vector probability diffusion[END_REF][START_REF] Tang | Direction diffusion. International Conference on Computer Vision (ICCV)[END_REF][START_REF] Tschumperlé | Diffusion tensor regularization with constraints preservation[END_REF][START_REF] Tschumperlé | Regularization of orthonormal vector sets using coupled PDE's[END_REF] since it detects edges and corners in a good way, and is easy to compute. Indeed, it does not require an eigenvalue decomposition of G as the other norms did, because

N + = ∇I = trace(G) = n i=1 ∇I i 2 (4)
Moreover, the norm N + has the interesting property of giving preferences to certain corners (Fig. 3d). This is very valuable for image restoration purposes, since the smoothing can be attenuated on high-curvature structures which are classically hard to preserve.

Note that for the scalar case (n = 1), the structure tensor calculus reduces to : In this case, the eigenvectors θ 1 +/-and the eigenvalues λ 1 +/-of G 1 are :

when n = 1 , dI 2 = dX G 1 dX where G 1 = ∇I∇I T = I 2 x I x I y I x I y I 2 y (a) Color checkerboard (size 40 × 40) (b) N = √ λ+ (c) N-= √ λ+ -λ- (d) N+ = √ λ+ + λ-
           θ 1 -= ξ = ∇I ⊥ ∇I θ 1 + = η = ∇I ∇I associated to    λ 1 -= 0 λ 1 + = ∇I 2
Basically, it means that 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, we can use it as a measure in many image analysis processes involving multi-channel images (not only for regularization algorithms). For instance, color edge detection may be performed by finding thresholded local maxima of the N + norm (Fig. 5 and [START_REF] Koschan | A comparative study on color edge detection[END_REF][START_REF] Tschumperlé | Constrained and unconstrained PDE's for vector image restoration[END_REF][START_REF] Tschumperlé | Diffusion PDE's on Vector-Valued images[END_REF]). This vector geometry computation has also been integrated as a measure of contours in some multi-channel image segmentation methods [START_REF] Sapiro | Vector-Valued Active Contours[END_REF][START_REF] Sapiro | Color snakes[END_REF]. For all reasons given above, the norm N + = λ + + λ -associated to the Di Zenzo geometry is probably one of the best measure for detecting local variations in multi-channel images and will be considered as it in the next parts of this chapter.

PDE-based Smoothing of Multi-Valued Images : A Review

We review and propose a classification of classical smoothing methods based on diffusion PDE's into three different approaches, related to different interpretation levels of the regularization processes, from the most global to the most local ones. For each section, we will start describing the original idea for scalar-valued images, then extending it for multi-channel datasets.

Variational methods

Contrary to the formulation of the original Perona-Malik equation, several methods have been proposed to apprehend the problem of image regularization as a global minimization procedure, within a variational framework. Formalisms described in [START_REF] Aubert | Mathematical Problems in Image Processing: Partial Differential Equations and the Calculus of Variations[END_REF][START_REF] Chambolle | Image recovery via total variation minimization and related problems[END_REF][START_REF] Charbonnier | Deterministic edge-preserving regularization in computed imaging[END_REF][START_REF] Sapiro | Geometric Partial Differential Equations and Image Analysis[END_REF][START_REF] Weickert | Anisotropic Diffusion in Image Processing[END_REF] among numerous references, contributed to define generic energy functionals measuring global image variations. The idea is that minimizing adapted variation functionals will 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 gives a very unifying way to proceed : A noisy scalar image I noisy can be regularized by minimizing the following φ-functional :

min I:Ω→R E(I) = Ω φ( ∇I ) dΩ (5)
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) = I noisy ∂I ∂t = div φ ′ ( ∇I ) ∇I ∇I (6)
Different choices of functions φ lead to different proposed regularization methods. One especially finds the simple isotropic smoothing (equivalent to a Gaussian convolution), as introduced by Tikhonov [START_REF] Tikhonov | Regularization of incorrectly posed problems[END_REF], as well as the well-known Perona-Malik [START_REF] Perona | Scale-space and edge detection using anisotropic diffusion[END_REF] and Total Variation (TV) anisotropic flows [START_REF] Rudin | Nonlinear total variation based noise removal algorithms[END_REF]. Lot of regularization methods acting on scalar-valued images have been unified by the φ-function formalism (Fig. 6).

Function name

φ(s) Reference Tikhonov s 2 [99] Perona-Malik 1 -exp(-s 2 /K 2 ) [80] Minimal surfaces 2 √ 1 + s 2 -2 [31] Geman-McClure s 2 /(1 + s 2 ) [46] Total Variation s [85] Green 2log(cosh(s)) [52]
Figure 6: List of different φ-functions and corresponding references.

A natural extension of the φ-functionals for the regularization of multi-channel images I could consist in minimizing the following cost functional E(I) measuring a global multi-channel image variation :

min I:Ω→R n E(I) = Ω φ(N (I)) dΩ (7)
where N (I) is one of the three local variation norms defined in section [START_REF] Abraham | Manifolds, Tensor Analysis, and Applications[END_REF]. 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 rather quite natural to minimize a functional defined by a function ψ : R 2 → R of two variables instead of a single one. The ψ-functional below is thus a more complete extension of the φ-function formulation for multi-channel images.

min I:Ω→R n E(I) = Ω ψ(λ + , λ -) dΩ (8)
The Euler-Lagrange equations of ( 8) can be derived, and reduce to a simple form of divergence-based PDE (see Appendix A for details about this Euler-Lagrange derivation) :

∂I i ∂t = div ∂ψ ∂λ + θ + θ T + + ∂ψ ∂λ - θ -θ T -∇I i (i = 1..n) (9)
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 [START_REF] Blomgren | Color tv: Total variation methods for restoration of vector-valued images[END_REF][START_REF] Pardo | Vector probability diffusion[END_REF][START_REF] Tang | Diffusion of general data on non-flat manifolds via harmonic maps theory : The direction diffusion case[END_REF] :

ψ(λ + , λ -) = φ( λ + + λ -)
or the Beltrami flow framework [START_REF] Kimmel | Image processing via the Beltrami operator[END_REF][START_REF] Kimmel | Images as embedded maps and minimal surfaces: movies, color, texture, and volumetric medical images[END_REF][START_REF] Kimmel | Geometric-variational approach for color image enhancement and segmentation[END_REF][START_REF] Sochen | Diffusions and confusions in signal and image processing[END_REF][START_REF] Sochen | A geometrical framework for low level vision[END_REF][START_REF] Sochen | On affine invariance in the Beltrami framework for vision[END_REF] :

ψ(λ + , λ -) = (1 + λ + )(1 + λ -)
Note that this last approach is also equivalent to define the minimizing functional E(I) as a Polyakov action which is actually 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, here by finding a minimal surface (Fig. 7).

Despite the interesting global geometric interpretation of variational formulations, such methods clearly lacks in flexibility. Indeed, they are formulated as global minimizations processes, despite the local geometric smoothing properties that are intrinsically desired for regularization purposes. Such PDE's are obtained by the Euler-Lagrange derivation of a functional and cannot thus be finely tuned to adapt themselves to local geometric cases (contours, corners, etc.). Unfortunately, this adaptability is primordial in many situations especially when the level of noise is high.

Divergence-based diffusion PDE's

One level of flexibility for designing regularization PDE's has been reached with the introduction of more generic divergence expressions [START_REF] Aubert | Mathematical Problems in Image Processing: Partial Differential Equations and the Calculus of Variations[END_REF][START_REF] Alvarez | Image selective smoothing and edge detection by nonlinear diffusion (II)[END_REF][START_REF] Kornprobst | Nonlinear operators in image restoration[END_REF][START_REF] Sapiro | Geometric Partial Differential Equations and Image Analysis[END_REF][START_REF] Weickert | Anisotropic Diffusion in Image Processing[END_REF]. 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. In one hand, this gives more freedom to design regularization PDE's that better fit local constraints. On the other hand, one often loses the global interpretation of the regularization process : generally, such designed equations do not correspond to a functional minimization anymore. Historically, authors of [START_REF] Alvarez | Image selective smoothing and edge detection by nonlinear diffusion (II)[END_REF] 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 regularization of scalarvalued images :

∂I ∂t = div (g( ∇I * G σ ) ∇I)
where G σ is a 2D normalized Gaussian function. This has initially been done to ensure the well-posedness of the regularization formulations. But, it also appeared that it allowed 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 multi-channel images has been more recently proposed by Weickert in [START_REF] Weickert | Theoretical foundations of anisotropic diffusion in image processing[END_REF][START_REF] Weickert | Coherence-enhancing diffusion of colour images[END_REF][START_REF] Weickert | A review of nonlinear diffusion filtering[END_REF][START_REF] Weickert | Anisotropic Diffusion in Image Processing[END_REF]. Basically, the idea consists in considering image pixels as chemical concentrations or temperatures which diffuse with respect to some physical laws (Fick Law and continuity equations). He proposed this very generic divergence-based equation, parametrized by a field D : Ω → P(2) of 2 × 2 diffusion tensors :

∂I i ∂t = div (D∇I i ) (i = 1..n) (10)
The tensor field D defines a gradient flux and controls then the local diffusion behavior of the smoothing process [START_REF] Ashikhmin | Synthesizing Natural Textures[END_REF]. Note that the ψ-functional formalism described in previous section (2.1) is just a particular case of the PDE [START_REF] Ashikhmin | Synthesizing Natural Textures[END_REF] with D defined as :

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 σ (3) such that :

   u = θ + v = θ - and          λ 1 = β λ 2 = β (if λ + = λ -) β + (1 -β)exp -C (λ+-λ-) 2 else (11) (C > 0 and β ∈ [0, 1] are user-fixed parameters of the method).
The tensor D is computed at each image point as :

D = λ 1 uu T + λ 2 vv T .
It is worth to notice that the tensor field D is the same for all image channels I i , ensuring that all I i are smoothed by a common multi-channel geometry, which takes the correlation between image channels into account (since D depends on G σ ), contrary to a uncorrelated channel-by-channel approach.

Weickert assumed that the tensor shape at each point X = (x, y) of the field D give the preferred smoothing geometry at this point. The idea behind the choice [START_REF] Barash | Bilateral filtering and anisotropic diffusion : Towards a unified viewpoint[END_REF] was then :

-On almost constant regions, we have λ + ≃ λ -≃ 0 and then we will get λ 1 ≃ λ 2 ≃ β, i.e D ≃ α I d . The tensor D is then defined to be isotropic in flat regions.

-Along image contours, we have

λ + ≫ λ -≫ 0 and as a result, λ 2 > λ 1 > 0
Here, the diffusion tensor D will be then anisotropic, mainly directed by the smoothed direction θ -of the image contours.

However, it is important to notice that the amplitudes and directions of the local smoothing performed by the divergencebased PDE [START_REF] Ashikhmin | Synthesizing Natural Textures[END_REF] are actually 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 simple following example. Suppose we want to anisotropically smooth a scalar image I : Ω → R everywhere along the gradient direction ∇I ∇I with a constant strength of 1. This is of course 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 :

∀X ∈ Ω, D (X) = ∇I ∇I ∇I ∇I
T leading to the simplification of [START_REF] Ashikhmin | Synthesizing Natural Textures[END_REF] as

∂I ∂t = div 1 ∇I 2 ∇I∇I T ∇I = div (∇I) = ∆I
where ∆I = ∂ 2 I ∂x 2 + ∂ 2 I ∂y 2 stands for the Laplacian of I. As noticed in [START_REF] Koenderink | The structure of images[END_REF], 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 σ = √ 2 t. So, this particular choice of anisotropic tensors D leads to an isotropic smoothing behavior, without preferred smoothing orientations. Note that choosing D = I d (identity matrix) would give exactly 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 the equation [START_REF] Ashikhmin | Synthesizing Natural Textures[END_REF] implicitly depends on the spatial variations of the tensor field D. Clearly, the divergence equation [START_REF] Ashikhmin | Synthesizing Natural Textures[END_REF] hampers the design of a significant pointwise smoothing behavior.

Oriented heat flows

Oriented heat flows, also named oriented Laplacians formulations consider that a local smoothing process can be decomposed into two orthogonal and uni-dimensional heat flows respectively oriented along two directions u 1 and u 2 (these vectors forming an orthonormal basis) associated with two smoothing amplitudes c 1 and c 2 . The smoothing amplitudes and orientations are naturally different for each image point, since they adapt themselves to the local configuration of the image (Fig. 8). The resulting equation is written as the sum of these two heat flows :

∂I ∂t = c 1 I u1u1 + c 2 I u2u2 (12)
where u 1 and u 2 are unit orthogonal vectors and c 1 , c 2 ≥ 0. I u1u1 and I u2u2 denote the second derivatives of I in the directions u 1 and u 2 and their vector components are formally defined as : ∀i = 1..n,

I iu 1 u 1 = u T 1 H i u 1 and I iu 2 u 2 = u T 2 H i u 2
where H i is the Hessian of I i , defined on each point X ∈ Ω by Here, the diffusion behavior is entirely defined by the knowledge of the smoothing directions u 1 , u 2 and the corresponding weights c 1 and c 2 .

H i = I ixx I ixy I ixy I iyy = ∂ 2 Ii ∂x 2 ∂ 2 Ii ∂x∂y ∂ 2 Ii ∂x∂y ∂ 2 Ii ∂y 2 (13)
This formulation has been first proposed in [START_REF] Kornprobst | Contributions à la Restauration d'Images et à l'Analyse de Séquences: Approches Variationnelles et Solutions de Viscosité[END_REF][START_REF] Kornprobst | Image restoration via PDE's[END_REF] for the regularization of scalar-valued images I, with the following choices for c 1 , c 2 and u 1 , u 2 :

     u 1 = ξ = ∇I ⊥ ∇I u 2 = η = ∇I ∇I and    c 1 = 1 c 2 = g( ∇I * G σ )
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 the edges ξ, even on very high gradients since c 1 = 1 everywhere in the image. The general formalism of oriented Laplacians [START_REF] Barbaresco | Spatial denoising of statistical parameters estimation by Beltrami diffusion on embedding siegel space[END_REF] allows also to find other well-known equations, such as the mean curvature flow ∂I ∂t = I ξξ , obtained with c 1 = 1, c 2 = 0, u 1 = ξ and v 2 = η [START_REF] Deriche | Les EDP en traitement des images et vision par ordinateur[END_REF]. Note that this was not possible to get with divergencebased expressions [START_REF] Ashikhmin | Synthesizing Natural Textures[END_REF]. Other works for scalar image regularization has used similar variants of this technique [START_REF] Carmona | Adaptive Smoothing Respecting Feature Directions[END_REF][START_REF] Krissian | Multiscale Analysis : Application to Medical Imaging and 3D Vessel Detection[END_REF].

Ringach and Sapiro [START_REF] Sapiro | Anisotropic diffusion of multivalued images with applications to color filtering[END_REF] proposed an extension of the mean curvature flow I t = I ξξ for multi-channel images, using an oriented Laplacian formulation [START_REF] Barbaresco | Spatial denoising of statistical parameters estimation by Beltrami diffusion on embedding siegel space[END_REF]. They naturally used the Di Zenzo attributes to incorporate informations on the multichannel geometry in their proposed equation :

∂I ∂t = g(λ + -λ -) I θ-θ- (14)
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 I i are smoothed along a common vector edge direction with a common intensity. Despite this great idea, some drawbacks subsist :

-The coherence norm N -= λ +λ -is used here as a measure of vector-valued variations, in order to reduce diffusion on image contours. This may not be a good choice since some corner structures will not respond highly to the norm N -, as explained in section (1), and will be then over-smoothed.

-In flat regions (N -→ 0), the diffusion is made along a single direction θ -, which is mainly directed by the noise since no coherent structures exist in these regions. Undesired texture effects result from this mono-directional 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 I i follow a common direction). Isotropic smoothing would be more adapted in order to remove noise in such flat regions.

Trace-based diffusion PDE's

A simple generalization of oriented Laplacians have been proposed in [START_REF] Tschumperlé | PDE's Based Regularization of Multivalued Images and Applications[END_REF][START_REF] Tschumperlé | Deriche Vector-Valued Image Regularization with PDE's : A Common Framework for Different Applications[END_REF][START_REF] Tschumperlé | Deriche Vector-Valued Image Regularization with PDE's : A Common Framework for Different Applications[END_REF]. The idea relies on the use of a generic diffusion tensor field T : Ω → P(2) to describe the diffusion geometry of the equation ( 12), instead of separately describing local directions θ + , θ -and amplitudes c 1 , c 2 of smoothing. Actually, the proposed equation was just a rewrite of the previous PDE [START_REF] Barbaresco | Spatial denoising of statistical parameters estimation by Beltrami diffusion on embedding siegel space[END_REF], using a trace operator :

∀i = 1, .., n, ∂I i ∂t = trace (TH i ) (15)
where H i stands for the Hessian of I i (13) and the tensor field T is computed as :

∀X ∈ Ω, T (X) = c 1 u 1 u T 1 + c 2 v 2 v T 2
Note that in this case, each channel I i of I is also smoothed with a common tensor field T.

Actually, equations ( 12) and ( 15) are strictly equivalent, but this last one makes clearly appear the separation of the smoothing geometry (defined by the tensor field T) from the smoothing itself. This is actually close to the idea of the Weickert's method that led to the divergence PDE [START_REF] Ashikhmin | Synthesizing Natural Textures[END_REF] : the regularization problem simplifies now to the design of a tensor field T adapted to the considered application. But in the case of trace-based PDE's, 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, equation ( 15) has a simple geometric interpretation in terms of local filtering with oriented Gaussian kernels. Indeed, let consider first that T is a constant tensor field. Then, it can be demonstrated that the formal solution of the PDE ( 15) is (see Appendix B for details) :

I i (t) = I i (t=0) * G (T,t) (i = 1..n) (16)
where * stands for the convolution operator and G (T,t) is an oriented Gaussian kernel, defined by :

G (T,t) (X) = 1 4πt exp - X T T -1 X 4t with X = (x y) T (17)
This is in fact a generalization of the Koenderink's idea [START_REF] Koenderink | The structure of images[END_REF], who proved this property in the field of computer vision for the simpler case of the isotropic diffusion tensor T = I d , resulting in the well-known heat flow equation : ∂Ii ∂t = ∆I i . Fig. 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. It is worth to notice that the Gaussian kernels G (T,t) give 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), i.e. 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 G T,t (X) over the image I. Fig. 10 illustrates three examples of spatially varying tensor fields T, represented with fields of ellipsoids, and the corresponding evolutions of (15) on a color image. As before, the shape of each tensor T gives the exact geometry of the local smoothing process performed by the trace-based PDE (15) point by point. As 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 will code similarly the preferred local smoothing directions as well as the desired smoothing amplitudes along these directions, for each image pixel X. Of course, T (X) must depend 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 in [START_REF] Tschumperlé | PDE's Based Regularization of Multivalued Images and Applications[END_REF][START_REF] Tschumperlé | Deriche Vector-Valued Image Regularization with PDE's : A Common Framework for Different Applications[END_REF][START_REF] Tschumperlé | Deriche Vector-Valued Image Regularization with PDE's : A Common Framework for Different Applications[END_REF] is : where p 1 , p 2 ∈ R are parameters of the proposed method.

c 1 = f - (λ+,λ-) = 1 (1 + λ + + λ -) p1 and c 2 = f + (λ+,λ-) = 1 (1 + λ + + λ -) p2 with p 1 < p 2 (18
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 mostly 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 ≃ I d (identity matrix). This is one possible choice for f -, f + in order to satisfy basic image denoising requirements. Actually, this is quite natural to design a smoothing behavior from the image structure before applying the regularization process itself.

The trace-based equation ( 15) has been a great attempt to separate the smoothing geometry from the smoothing process itself, while providing a geometrical interpretation on how the smoothing is performed. It proved some natural links between PDE's and other local filtering techniques, as the Bilateral Filtering [START_REF] Barash | Bilateral filtering and anisotropic diffusion : Towards a unified viewpoint[END_REF][START_REF] Tomasi | Bilateral filtering for gray and color images[END_REF]. Another similar approach based on non-Gaussian convolution kernels has been also proposed for the specific case of the Beltrami Flow in [START_REF] Sochen | Diffusions and confusions in signal and image processing[END_REF].

But the fact that the trace equation ( 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 also best known as the "mean curvature flow" effect, characterized by the PDE ∂I ∂t = ∂ 2 I ∂θ -2 . This problem is illustrated on Fig. 11b and Fig. 12b where [START_REF] Bertalmio | Image inpainting[END_REF] has been applied on synthetic and real color image and T has been defined as [START_REF] Black | Robust anisotropic diffusion[END_REF] (then f -= 0). One can easily see how image structures are subject to the mean curvature flow effect, resulting in rounding the corners of the square in Fig. 11b, or in blending parallel thin curved structures in Fig. 12b. To avoid this over-smoothing effect, one may try to stop the action of the diffusion PDE on corners (by vanishing tensors T (X) there, i.e 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 under-smoothing on edges may occur when limiting the diffusion too much on regions with high intensity variations (Fig. 11c). There is a difficult trade-off between complete noise removal and preservation of curved structures, when using trace-based PDE's [START_REF] Bertalmio | Image inpainting[END_REF]. Actually, this kind of regularization processes does not care about the curvature of the smoothing directions, and by extension, of the curvature of the image contours. Taking this curvature into account is a very desirable goal and has motivated the work presented in the sequels : in section 3, we propose a class of trace-based regularization PDE's that smooth an image I along a tensor field T, while implicitly taking curvatures of specific integral curves of T into account. Roughly speaking, the method will locally filter the image with curved Gaussian kernels when necessary, in order to better preserve image structures. For comparison purposes, results of the curvature-preserving equation is shown on Fig. 11d and Fig. 12c.

Links between existing regularization methods

The link between these three formulations is generally not trivial, especially for vector-valued images. Actually, it is well known for the classical case of φ-functional regularization of scalar images (n = 1) : One can start from a regularizing (A) : min

I:Ω→R Ω φ( ∇I ) dΩ (19)
⇒ (B) :

∂I ∂t = div φ ′ ( ∇I ) ∇I ∇I ⇒ (C) : ∂I ∂t = φ ′ ( ∇I ) ∇I I ξξ + φ ′′ ( ∇I ) I ηη
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 underlying tensor

D = φ ′ ( ∇I ) ∇I I d in the divergence expression.
It is also worth to mention that this global-to-local path (from variational to trace-based equations) can be rarely be followed in the inverse order.

For multi-channel images, this link can be also found :

(A) min

I:Ω→R n Ω ψ(λ + , λ -) dΩ (20)
⇒ (B) ∂I i ∂t = div (D ∇I i ) where D = ∂ψ ∂λ + θ + θ T + + ∂ψ ∂λ - θ -θ T - ⇒ (C) ∂I i ∂t = n j=1 trace (δ ij D + Q ij )H j
where the δ ij is the Kronecker symbol (δ ij = 1 when i = j, and 0 elsewhere) , Q ij designates a family of n 2 tensors (i, j = 1..n), defined as the symmetric parts of the following matrices P ij (i.e, Q ij = (P ij + P ij T )/2 ) :

P ij = α ∇I T i ∇I j I d + 2 ∂α ∂λ + θ + θ T + + ∂α ∂λ - θ -θ T -∇I j ∇I T i G + 2 (α + ∂β ∂λ + )θ + θ T + + (α + ∂β ∂λ - )θ -θ T -∇I j ∇I T i with α = f1(λ+,λ-)-f2(λ+,λ-)
λ+-λ-and β = λ+f2(λ+,λ-)-λ-f1(λ+,λ-)

Curvature-Preserving PDE's

The framework of curvature-preserving PDE's, first introduced in [START_REF] Tschumperlé | Fast Anisotropic Smoothing of Multi-Valued Images using Curvature-Preserving PDE's[END_REF] defines a specific variant of multi-channel diffusion PDE's. Its goal is to provide a generic tensor-driven regularization method as the divergence-based PDE [START_REF] Ashikhmin | Synthesizing Natural Textures[END_REF] and trace-based PDE [START_REF] Bertalmio | Image inpainting[END_REF], 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.

The single direction case

To illustrate the general idea of curvature-preserving PDE's, we first focus on image regularization along a vector field w : Ω → R 2 instead of a tensor field T. We consider then a local smoothing everywhere along a single direction w w , with a smoothing strength w . We denote the two spatial components of w by w (X) = (u (X) v (X) ) T . The curvature-preserving regularization PDE that smoothes I along w is defined by :

∀i = 1, . . . , n, ∂I i ∂t = trace ww T H i + ∇I T i J w w (21)
where J w stands for the Jacobian of w , and H i is the Hessian of I i .

J w =   ∂u ∂x ∂u ∂y ∂v ∂x ∂v ∂y   and H i =    ∂ 2 Ii ∂x 2 ∂ 2 Ii ∂x∂y ∂ 2 Ii ∂x∂y ∂ 2 Ii ∂y 2   
The PDE (21) adds a term ∇I T i J w w to the trace-based equation ( 15) that smoothes I along w with locally oriented Gaussian kernels (see section 2.4). This extra term naturally depends on the variation of the vector field w. Let us explain how ( 21) is related to w.

Let C X (a) be the curve defining the integral curve of w, starting from X and parametrized by a ∈ R :

     C X (0) = X ∂C X (a) ∂a = w(C X (a) ) (22)
When a → +∞ the integral curve C X (a) is tracked forward, and backward when a → -∞ (Fig. 13). We denote by F the family of integral curves of w. A second-order Taylor development of C X (a) around a = 0 is :

C X (h) = C X (0) + h ∂C X (a) ∂a |a=0 + h 2 2 ∂ 2 C X (a) ∂a 2 |a=0 + O(h 3 ) = X + hw (X) + h 2 2 J w (X) w (X) + O(h 3 )
with h → 0, and O(h n ) = h n ǫ n . Then, we can compute a second-order Taylor development of I i (C X (a) ) around a = 0, which corresponds to the variations of the image intensity near X when following the integral curve C X :

I i (C X (h) ) = I i X + hw (X) + h 2 2 J w (X) w (X) + O(h 3 ) = I i (X) + h∇I i T (X) (w (X) + h 2 J w (X) w (X) ) + h 2 2 trace w (X) w T (X) H i(X) + O(h 3 )
The term trace w (X) w T (X) H i(X) = ∂ 2 Ii ∂w 2 corresponds to the second directional derivative of I i along w.

The second derivative of the function a → I i (C X (a) ) at a = 0 is then :

∂ 2 I i (C X (a) ) ∂a 2 |a=0 = lim h→0 1 h 2 I i (C X (h) ) + I i (C X (-h) ) -2I i (C X (0) ) = lim h→0 1 h 2 h 2 ∇I T i J w (X) w (X) + h 2 trace w (X) w T (X) H i(X) + O(h 3 ) = trace w (X) w T (X) H i(X) + ∇I T i J w (X) w (X) (23)
Note that this is exactly the right term in the proposed curvature-preserving PDE [START_REF] Blomgren | Color tv: Total variation methods for restoration of vector-valued images[END_REF].

Actually, ( 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 that Y = C X (ǫ) . Indeed, C X and C Y describe the same curve [START_REF] Buerkle | Transport and diffusion in timedependent flow visualization[END_REF] with different parametrization :

∀a ∈ R, C Y (a) = C X (ǫ+a) . As (21) is verified on Y, then ∂Ii(C X (a) ) ∂t |a=ǫ = ∂ 2 Ii(C X (a) ) ∂a 2
|a=ǫ . This is obviously true for ǫ ∈ R since ( 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, ∂I i (C (a) ) ∂t = ∂ 2 I i (C (a) ) ∂a 2 (24)
We may recognize in (24) a one-dimensional heat flow constrained on C. This is actually very different from a heat-flow oriented by w, as in the formulation ∂Ii ∂t = ∂ 2 Ii ∂w 2 since the curvatures of integral curves of w are now implicitly taken into account. In particular, the constrained equation ( 21) has the interesting property to vanish when image intensities are perfectly constant on the integral curve C, whatever the curvature of C is. 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 PDE's [START_REF] Ashikhmin | Synthesizing Natural Textures[END_REF], [START_REF] Bertalmio | Image inpainting[END_REF]. This curvature-preserving property of ( 21) is illustrated on Fig. 11d and Fig. 12b. The constrained equation ( 21) is an elliptic PDE since the matrix ww T is positive definite. The existence and unicity of the solutions of ( 21) are not directly approached here. Anyway, in section 3.2, we show that its solution can be approximated by the technique of Line Integral Convolutions, which is a well-posed analytical approach.

Curvature-Preserving PDE's and Line Integral Convolutions

Line Integral Convolutions (LIC) have been first introduced in [START_REF] Cabral | Imaging vector fields using line integral convolution[END_REF] as a technique to render a textured image I LIC that represents a vector field w : Ω → R 2 . The idea, originally expressed under a discrete formula, consists in smoothing an image I noise -containing only noise -by averaging its pixel values along the integral curves of w. Actually, a continuous formulation of a LIC is then :

∀X ∈ Ω, I LIC (X) = 1 N +∞ -∞ f (p) I noise (C X (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 [START_REF] Buerkle | Transport and diffusion in timedependent flow visualization[END_REF] 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 noticed in section 3.1, the curvature-preserving PDE (21) can be seen as the one-dimensional heat flow (24) constrained on the integral curve C X ∈ F. Using the variable substitution L (a) = I(C X (a) ), [START_REF] Carmona | Adaptive Smoothing Respecting Feature Directions[END_REF] can be also written as ∂L ∂t (a) = L ′′ (a) . The solution L [t] at time t is known to be the convolution of L [t=0] by a normalized Gaussian kernel G t (see [START_REF] Deriche | Les EDP en traitement des images et vision par ordinateur[END_REF][START_REF] Koenderink | The structure of images[END_REF]) :

L [t] (a) = +∞ -∞ L [t=0] (p) G t(a-p) dp with G t(p) = 1 √ 4πt exp - p 2 4t (26)
Substituting L in (5.4) with a = 0, and remembering that C X (0) = X and G t(-p) = G t(p) :

∀X ∈ Ω, I [t] (X) = +∞ -∞ I [t=0] (C X (p) ) G t(p) dp (27)
The equation ( 27) is a particular form of the continuous LIC-based formulation (25) with a Gaussian weighting function f = G t . Here, the normalization factor is

N = +∞ -∞ G t(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. This kind of anisotropic image smoothing considers then a curved filtering, instead of just an oriented one. Applying this setting on a multi-channel image I, with w being the lowest eigenvector of the structure tensor field G (i.e. the contour direction) allows the anisotropic smoothing of I with edge preservation, even if these edges are curved. This is illustrated on Fig. 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 sub-pixel precision. The streamlines have been computed with a 4 th -order Runge-Kutta scheme.

Note that ( 27) is an analytical solution of (21) when w does not evolve over time. This property is generally not verified when dealing with general nonlinear regularization PDE's, where the smoothing geometry is re-evaluated at each time step (thus defining a temporal non-linearity). In order to get rid of this kind of non-linearity, we will then to perform several successive iterations of the LIC scheme [START_REF] Chan | Total variation denoising and enhancement color images based on the cb and hsv color models[END_REF], where the vector field w is updated at each iteration. This is actually a good way of approximating [START_REF] Blomgren | Color tv: Total variation methods for restoration of vector-valued images[END_REF]. 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 equation ( 21) will be efficiently discretized by several iterations of the LIC formulation ( 27). This will be detailed in section 4.

Between Traces and Divergences

We illustrate here how the curvature-preserving PDE ( 21) may be regarded compared to trace and divergence expressions [START_REF] Ashikhmin | Synthesizing Natural Textures[END_REF], [START_REF] Bertalmio | Image inpainting[END_REF], for the case of single direction smoothing T = ww T . In this case, the divergence PDE (10) may be developed as :

div ww T ∇I i = div   u 2 ∂Ii ∂x + uv ∂Ii ∂y uv ∂Ii ∂x + v 2 ∂Ii ∂y   = u 2 ∂ 2 I i ∂x 2 + 2uv ∂ 2 I i ∂x∂y + v 2 ∂ 2 I i ∂y 2 + ∇I T i     2u ∂u ∂x + u ∂v ∂y + v ∂u ∂y 2v ∂v ∂y + u ∂v ∂x + v ∂u ∂x     = trace ww T H i + ∇I T i         u ∂u ∂x + v ∂u ∂y u ∂v ∂x + v ∂v ∂y     +     u ∂u ∂x + u ∂v ∂y v ∂u ∂x + v ∂v ∂y         = trace ww T H i + ∇I T i J w w + div(w)∇I T i w
Thus, we recognize in these three different terms :

-The first term corresponds to the trace PDE [START_REF] Bertalmio | Image inpainting[END_REF], that smoothes locally I along w, using oriented Gaussian kernels.

-The two first terms correspond to the curvature-constrained regularization PDE [START_REF] Blomgren | Color tv: Total variation methods for restoration of vector-valued images[END_REF], that smoothes 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) that performs local diffusions of I along w. This last term div(w)∇I T i w is mainly responsible for the perturbations of the effective smoothing direction, as described in section 2.2. 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 [START_REF] Bertalmio | Image inpainting[END_REF] and divergence formulations [START_REF] Ashikhmin | Synthesizing Natural Textures[END_REF], and allows at the same time the full respect of the pre-defined smoothing directions w, while preserving curved images structures. Note that we can also write the curvature-preserving PDE ( 21) as a divergence-based PDE minus a constraint term :

trace ww T H i + ∇I T i J w w = div ww T ∇I i -div(w)∇I T i w
Two particular cases of directions w are worth studying, in the case of scalar-valued images (n = 1) :

-When w = ∇I ⊥ ∇I (isophote direction), then ∇I T J w w = -I ww , vanishing then the velocity of the curvaturepreserving evolution equation [START_REF] Blomgren | Color tv: Total variation methods for restoration of vector-valued images[END_REF], by counterbalancing the trace-based term (which is nothing more than the mean curvature motion in this case). No smoothing will be then performed. This is quite natural since pixel 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 ww T ∇I i also vanishes here.

-When w = ∇I ∇I (gradient direction), then ∇I T J w w = 0, and the velocity of the curvature-preserving PDE (21) becomes simply I ww , which really corresponds to a smoothing of the image along the gradient direction (the same as the unconstrained trace-based PDE [START_REF] Bertalmio | Image inpainting[END_REF]). 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 allows to better understand the difference of regularization behaviors between the trace, divergence and curvature-preserving formulations. Note also that in case where w is a divergence free field (i.e div(w) = 0), the divergence-based PDE [START_REF] Ashikhmin | Synthesizing Natural Textures[END_REF] and the curvaturepreserving formulation ( 21) are strictly equivalent. This is very rarely the case anyway.

Extension to multi-directional smoothing

In [START_REF] Tschumperlé | Fast Anisotropic Smoothing of Multi-Valued Images using Curvature-Preserving PDE's[END_REF], the single-direction smoothing PDE [START_REF] Blomgren | Color tv: Total variation methods for restoration of vector-valued images[END_REF] has been extended so that it can deal with 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 represents 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 tensors-valued fields T are not as natural as with direction fields w.

To tackle this problem, we proposed to locally decompose a tensor-driven smoothing process into several vector-driven smoothing processes along different orientations. We first notice that

π α=0 a α a T α dα = π 2 I d where a α =   cos α sin α  
Then, any 2 × 2 tensor T may be written as :

T = 2 π √ T π α=0 a α a T α dα √ T where √ T = f + uu T + f -vv T stands for the square root of T = f + uu T + f -vv T . One can easily verify that ( √ T) 2 = T and ( √ T) T = √ T.
Thus, the tensor T may be decomposed as :

T = 2 π π α=0 √ Ta α a T α √ T T dα = 2 π π α=0 ( √ Ta α )( √ Ta α ) T dα ( 28
)
We have split the tensor T into a sum of atomic tensors ( √ Ta α )( √ Ta α ) T , each being purely anisotropic and directed only along the direction of the vector √ Ta α ∈ R 2 . The equation ( 28) naturally suggests to decompose any tensor-driven regularization PDE into a sum of single direction smoothing processes, each of them respecting the overall geometry T. For instance :

-If T = I d (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 = uu T (where u ∈ S 1 ), the tensor is purely anisotropic and : ∀α ∈ [0, π],

√ Ta α = (u T a α )u. The resulting smoothing will be then performed only along the direction u of the tensor T.

Then, using [START_REF] Chan | Non-texture inpainting by curvature-driven diffusions (CDD)[END_REF] and considering that each single direction smoothing must be done with a curvature-preserving approach (21), we end up with the following curvature-constrained regularization PDE, acting on a multi-channel image I : Ω → R n and driven by a tensor-valued smoothing geometry T :

∀i = 1, . . . , n, ∂I i ∂t = 2 π π α=0 trace ( √ Ta α )( √ Ta α ) T H i + ∇I T i J √ Taα √ Ta α dα
which can be simplified as :

∀i = 1, . . . , n, ∂I i ∂t = trace(TH i ) + 2 π ∇I T i π α=0 J √ Taα √ Ta α dα (29)
where a α = (cos α sin α) T , and J √ Taα stands for the Jacobian of the vector field Ω → √ Ta α . Note that this kind of smoothing decomposition along all orientations of the plane can be also found in [START_REF] Weickert | Anisotropic Diffusion Filters for Image Processing Based Quality Control[END_REF]. As in the single direction smoothing case, (29) may be seen as a trace-based equation [START_REF] Bertalmio | Image inpainting[END_REF], where an extra term has been added in order to respect the curvature of all integral lines passing through the tensor-valued geometry T.

Implementation considerations

In order to implement the regularization method [START_REF] Chan | Variational restoration of non-flat image features : Models and algorithms[END_REF], one can benefit from the LIC-based interpretation of curvaturepreserving PDE's presented in section 3.2. Indeed, we can explicitly discretize (29) by the following Euler scheme :

I [t+dt] = I [t] + 2dt N N -1 k=0 R( √ Ta α )
where α = kπ/N (in the interval [0, π]), dt is the usual temporal discretization step and R(w) represents a discretization of the mono-directional smoothing PDE velocity ( 21) that preserve curvatures along a vector field w. If we write this expression as :

I [t+dt] = 1 N N -1 k=0 I [t] + 2dt R( √ Ta α )
, we may express it as the averaging of different Gaussian-pondered LIC's along vector fields √ Ta α :

I [t+dt] = 1 N N -1 k=0 I [t] LIC( √ Taα)
,

where each Gaussian variance has a standard deviation dt.

Basically, the difficulty here is the LIC computation, which needs the tracking of integral curves of a vector field. Here, we used a very simple method based on the classical Runge-Kutta [START_REF] Press | Runge-Kutta Method in Numerical Recipes in FOR-TRAN: The Art of Scientific Computing[END_REF] integration scheme. Faster LIC implementations have been proposed in [START_REF] Stalling | Fast and Resolution Independent Line Integral Convolution[END_REF] 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 curvaturepreserving PDE (29) :

-Compute the smoothed structure tensor field G σ from I [t] :

G σ = G σ * n i=1     ∂I [t] i ∂x 2 ∂I [t] i ∂x ∂I [t] i ∂y ∂I [t] i ∂x ∂I [t] i ∂y ∂I [t] i ∂y 2    
σ will depend on the noise scale. We used relatively low values (between 0 and 1.5) for our experiments in section 5.

-Compute the eigenvalues λ + , λ -and eigenvectors θ + , θ -of G σ .

-Compute the smoothing geometry tensor field T from

G σ : T = 1 (1+λ + +λ -) p 1 θ -θ -T + 1 (1+λ + +λ -) p 2 θ + θ + T
-For all α in [0, π] (discretized with a user-fixed step d α ) :

-Compute the vector field w = √ T a α .

-Perform a Line Integral Convolution of I [t] along C X in the forward and backward directions.

-Average all LIC's computed in step 4.

The main parameters of the algorithm are p 1 , p 2 , σ, dt and the number of PDE iterations nb that are applied. The characteristics of this scheme, compared to the classical finite-difference one is :

-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 sub-pixel accuracy. Precise 4 th Runge-Kutta interpolation [START_REF] Press | Runge-Kutta Method in Numerical Recipes in FOR-TRAN: The Art of Scientific Computing[END_REF] is used to track the integral curves C in the image.

-It allows to choose 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 get the result, even if each iteration is more time-consuming. For our applications, presented in section 5, we were even able to choose only nb = 1 iteration with very large time steps dt. In fact, this leads to a rough approximation of ( 29), since we lost the temporal non-linearity property of the PDE. But for images with few noise, this gave surprisingly good results. Actually, the spatial non-linearity seems to play a more important role than the temporal non-linearity in the PDE evolution.

The smoothing is done as an averaging of multiple LIC's in different directions α. The choice of the discretization step d α is important in this context. Actually, in regions where the smoothing needs to be mostly anisotropic, only few values of α are necessary since in all cases, the smoothing will be done along the same single direction. But in homogeneous regions needing isotropic smoothing, a smaller d α will give much better results. Practically speaking, we chose d α = 45 o which is enough to get a good precision for isotropic smoothing. On Fig. 14, we illustrate the efficiency of the scheme, compared to the classical finite-difference one. A synthetic noisy image is anisotropically smoothed with the PDE (29), with p 1 = 0.01 and p 2 = 100 (smoothing mostly along isophotes θ -, with a strength of 1). The LIC-based scheme (Fig. 14c) clearly better preserves the structure along time t. This is due to the important role played by the sub-pixel accuracy property of the underlying LIC computation.

Applications

We illustrate the wide variety of problems that can be tackled by the regularization PDE's proposed throughout this chapter. Particularly, we use the curvature-preserving regularization approach described in section 3 in the following examples, dealing with color images. We show results of color image denoising, inpainting and resizing by nonlinear interpolation. Given processing time have been obtained on a 2.8Ghz i686 Intel processor.

Color image denoising and artifacts removal

Image denoising is of course 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. It generally leads to small random variations that affect pixels of the image and that must be removed.

In Fig. 15-19, we illustrate how the curvature-preserving PDE framework (29) can be successfully applied to remove such artifacts while preserving the essential structures of the processed images.

-Fig. 15 shows an application of the regularization method on the famous 512 × 512 "Baboon" color image, artificially degraded by adding uncorrelated Gaussian noise on (R, G, B). This color image has been then denoised with equation [START_REF] Chan | Variational restoration of non-flat image features : Models and algorithms[END_REF]. Thanks to the proposed LIC-based numerical implementation, only one PDE iteration has been necessary to denoise the image, with parameters p 1 = 0.5, p 2 = 0.7, σ = 1.5 and dt = 50. Processing time is 19.3 seconds for the entire image.

-Fig. 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 (29) allows to clearly remove the grains while preserving quite fine structures (palm tree leafs). Shown image is a 152 × 133 portion of the original one. Only one PDE iteration has been necessary, with p 1 = 0.5, p 2 = 0.7, σ = 1 and dt = 10. Processing time is 11 seconds for the entire image (size 586 × 367).

-Fig. 17 shows the restoration of a digital photograph shot under low luminosity conditions by a cellular phone. Such devices have usually poor quality cameras, leading to the apparition of important digital noise (more precisely, Poisson noise) on the acquired color images. Processed color image has a size of 262 × 280 and has been restored in 4 seconds (one PDE iteration), with parameters p 1 = 0.2, p 2 = 0.5, σ = 2, dt = 120. Note how the curvature-preserving PDE ( 29) is able to adapt itself locally to the multi-channel 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 image is such an example of interest. Fig. 18 illustrate the suppression of compression artifacts in color images. A JPEG compressed color image of size 283 × 249, where the JPEG quality ratio has been set to 10% is processed by the multi-channel image regularization algorithm. Usual block effects inherent to the DCT compression are visible on the compressed image. One PDE iteration is applied then, with p 1 = 0.5, p 2 = 0.9, σ = 2, dt = 200, in order to get the regularized result (right). Processing time is 5.4 seconds for the entire image.

-Fig. 19 illustrates another regularization with a different type of noise. This time, 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 p 1 = 0.5, p 2 = 0.8, σ = 1, 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 a true color images, even by storing them with a quantified palette.

Color image inpainting

Image inpainting is a very new and challenging application, which consists in 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 : Ω → R 3 , and one mask image M : Ω → {0, 1}. The inpainting algorithm must fill-in the regions where M (X) = 1, by the mean of some intelligent interpolations. Inpainting algorithms can be used for instance to remove various structures in images (scratches, logos or real objects) that are usually bigger than other image artifacts. Image inpainting has been first proposed as a method based on a variational formulation by Masnou and Morel [START_REF] Masnou | Level Lines Based Disocclusion[END_REF], followed by many solutions based on diffusion or transport PDE's [START_REF] Bertalmio | Image inpainting[END_REF][START_REF] Bertalmio | Simultaneous Structure and Texture Image Inpainting[END_REF][START_REF] Chan | Mathematical models for local deterministic inpaintings[END_REF][START_REF] Tschumperlé | Deriche Vector-Valued Image Regularization with PDE's : A Common Framework for Different Applications[END_REF][START_REF] Tschumperlé | Deriche Vector-Valued Image Regularization with PDE's : A Common Framework for Different Applications[END_REF]. It is also worth to cite some papers related to inpainting without use of PDE's [START_REF] Criminisi | Object Removal by Exemplar-based Inpainting IEEE Conference on Computer Vision and Pattern Recognition[END_REF][START_REF] Jia | Tang Image Repairing : Robust Image Synthesis by Adaptive N D Tensor Voting[END_REF], among others.

In this article, we see the inpainting process as a direct application of the proposed curvature-preserving PDE [START_REF] Chan | Variational restoration of non-flat image features : Models and algorithms[END_REF]. 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 isophotes directions is thus naturally done, reconstructing the missing parts of the image, since the performed smoothing tries to follow a coherent multi-valued image geometry computed from the known parts of the image. We illustrate the concept of inpainting with PDE's in Fig. 2021222324.

-Fig. 20 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 when it comes to reconstruct textured regions. For instance, the algorithm has not been able to reconstruct one eye of the woman on the left. This is anyway a very hard task, and it is not surprising that diffusion equations that perform only local smoothing fail in such complex cases.

-Fig. 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 [START_REF] Chan | Variational restoration of non-flat image features : Models and algorithms[END_REF] with parameters p 1 = 0.001, p 2 = 100, σ = 4, dt = 150. Note that p 1 << p 2 encourages smoothing only along the isophote directions with a strength of 1 everywhere. Another similar example of real object removal can be found in Fig. 24, where a cage automatically disappear in a color photograph of a parrot.

-Fig. 22-23 illustrate the reconstruction capabilities of our inpainting technique. Here, half of the pixels of two color images have been suppressed by masking them with checkerboard-shaped mask with squares 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 p 1 = 0.001, p 2 = 100, σ = 4 and dt = 50 This kind of 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 in this article, the initialization of the pixel values inside the inpainting masks at t = 0 has been done by white noise. Actually, the inpainting algorithm is not much dependent of the initialization step : the equation ( 29) diffuses neighborhood values inside the inpainting mask until convergence, and there is then a strong border condition. We didn't see much difference between different types of initialization (noise, zero-filling or linear interpolation).

Color image interpolation

With the same kind of techniques, one can easily perform image magnification by edge-preserving interpolation. Starting from a linear or bicubic interpolation of a small image, and applying the PDE [START_REF] Chan | Variational restoration of non-flat image features : Models and algorithms[END_REF] on the image (excepted on the original known pixels that form a sparse inpainting mask), we can compute magnified images that have been regularized while taking their local geometry into account. It allows to remove usual bloc or jagging effects inherent to classical linear interpolation techniques. This is a technique which is actually very similar to image inpainting with a very sparse grid for the inpainting mask.

-Fig. 25 illustrates one example of image resizing. An original 195 × 173 color image is resized by a factor ×4 with classical nearest-neighbor, linear and bicubic interpolations, then by the PDE-based technique [START_REF] Chan | Variational restoration of non-flat image features : Models and algorithms[END_REF]. It is quite clear that the non-linear regularization filter, driven by the image geometry, allows to remove the aliasing effects usually encountered with simple interpolation methods, while correctly preserving the small structures of the image.

Notice that the original known points of the images are always preserved during the regularizing PDE flow, ensuring that resizing back the image to its original dimension (sub-sampling) always results in the original input data.

Flow visualization

We present here a last application of regularization PDE's for visualization purposes. Considering a 2D vector field F : Ω → R 2 , we have several ways to visualize it. We can first use vector graphics (Fig. 26) (left), but we need to subsample the field since this kind 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 (29) but where T is directed by the directions of F, instead of the local geometry of I. It means that we have to define the tensor field D in (29) as :

∀X ∈ Ω, D (X) = 1 F F F T (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 [START_REF] Cabral | Imaging vector fields using line integral convolution[END_REF]. But it is interesting to see that when the PDE evolution time t goes by, a visualization scale-space of F is explicitly constructed (Fig. 27). Here, the used regularization equation [START_REF] Chan | Mathematical models for local deterministic inpaintings[END_REF] ensures that the smoothing of the pixels is done exactly in the direction of the flow F. This is not the case in [START_REF] Becker | PDE methods in flow simulation post processing[END_REF][START_REF] Buerkle | Transport and diffusion in timedependent flow visualization[END_REF][START_REF] Diewald | Anisotropic diffusion in vector field visualization on euclidian domains and surfaces[END_REF], where the authors proposed a similar idea using a divergence-based expression. Using similar divergence-based techniques would raise a risk of smoothing the image in false directions, as this has been pointed out in section 2.2.

Conclusion

Multi-channel image regularization is a fundamental process in many image processing and computer vision algorithms and it is then primordial to get a full control of this process, as well as to understand what the used equations are exactly doing from a geometric point of view. In this chapter, we have described the most classical PDE-based methods proposed for the regularization of multi-channel images and introduced a very efficient curvature-preserving framework that generally outperforms its competitors. This is not only due 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 multi-channel image regularization technique can play a role in a lot of image processing applications. The processing time, which was one of the famous drawback of PDE-based methods, is not a problem anymore. All these reasons makes the framework of multi-valued diffusion PDE's a very good choice for image regularization purposes. This has been illustrated in this chapter with results on color image denoising, inpainting and resizing. But many other applications may benefit from the proposed curvature-preserving framework.

It is worth to notice that 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 [START_REF] Tschumperlé | The C++ Template Image Processing Library[END_REF].

Appendix A

In this appendix, we demonstrate how the minimization of the functional

min I:Ω→R n E(I) = Ω ψ(λ + , λ -) dΩ (31)
can be performed by a gradient descent and its corresponding PDE flow. The Euler-Lagrange equations corresponding to the functional (31) are :

∂I i ∂t = div       ∂ψ ∂I ix ∂ψ ∂I iy       (i = 1..n) (32)
Actually, the vector ( ∂ψ ∂Ii x , ∂ψ ∂Ii y

) T can be written in a more comprehensive form. From the chain-rule property of the derivation, we have :

      ∂ψ ∂I ix ∂ψ ∂I iy       =       ∂λ + ∂I ix ∂λ - ∂I ix ∂λ + ∂I iy ∂λ - ∂I iy            ∂ψ ∂λ + ∂ψ ∂λ -      (33)
We know formally the expressions ∂ψ ∂λ± since the function ψ is directly defined from the λ ± . Finding the ∂λ± ∂Ii x and ∂λ± ∂Ii y is more tricky. Here is a simple way to proceed :

As the λ ± are the eigenvalues of the structure tensor G = (g kl ), we may decompose its derivatives (with respect to I ix and I iy ), in terms of derivatives with respect to the g kl : i.e (34) can be written as :

∂λ ± ∂I ix = k,
      ∂λ ± ∂I ix ∂λ ± ∂I iy       =      2 ∂λ ± ∂g 11 ∂λ ± ∂g 12 ∂λ ± ∂g 12 2 ∂λ ± ∂g 22      ∇I i (35)
Thus, one last obstacle remains to be crossed, that is finding the formal expressions of ∂λ± ∂g kl . Remind that the λ ± and θ ± are the eigenvalues and eigenvectors of the structure tensor G :

G = λ + θ + θ T + + λ -θ -θ T -
The derivation of this tensor, with respect to one of its coefficient g kl is :

∂G ∂g kl = ∂λ + ∂g kl θ + θ T + + ∂λ - ∂g kl θ -θ T - (36)
+ λ + ∂θ + ∂g kl θ T + + λ - ∂θ - ∂g kl θ T - + λ + θ + ∂θ T + ∂g kl + λ -θ - ∂θ T - ∂g kl
Moreover, as the θ ± are unitary and orthogonal eigenvectors, we have :

   θ T + θ + = θ T -θ -= 1 θ T + θ -= θ T -θ + = 0 and            ∂θ T + ∂g kl θ + = θ T + ∂θ + ∂g kl = 0 ∂θ T - ∂g kl θ -= θ T - ∂θ - ∂g kl = 0 (37)
We first multiply the equation ( 36) by θ T ± at the left, by θ ± at the right, then use the properties [START_REF] Coulon | A regularization scheme for diffusion tensor magnetic resonance images[END_REF]. It allows high simplifications, and leads to these two relations :

∂λ + ∂g kl = θ T + ∂G ∂g kl θ + and ∂λ - ∂g kl = θ T - ∂G ∂g kl θ - (38)
Equations ( 38) formally tell us how eigenvalues of a diffusion tensor G vary with respect to a particular coefficient g kl of G.

Actually, this interesting property can be proved for any symmetric matrix. For instance, authors of [START_REF] Papadopoulo | Estimating the jacobian of the singular value decomposition: Theory and applications[END_REF] proposed a similar demonstration in a purely matrix form, leading to the same result. They used it to deal with general covariance matrices.

Moreover in our case, the matrices ∂G ∂g kl are very simple to write :

∂G ∂g 11 = 1 0 0 0 , ∂G ∂g 12 = 0 1 1 0 and ∂G ∂g 22 = 0 0 0 1
With all these elements, we can express [START_REF] Chu | Matrix differential equations : A continuous realization process for linear algebra problems[END_REF] as :

      ∂λ + ∂I ix ∂λ + ∂I iy       = 2 θ + θ T + ∇I i and       ∂λ - ∂I ix ∂λ - ∂I iy       = 2 θ -θ T -∇I i (39)
Finally, replacing [START_REF] Deriche | Les EDP en traitement des images et vision par ordinateur[END_REF] in the Euler-Lagrange equations ( 33) and (32), gives the vector-valued gradient descent of the functional (8) :

min I:Ω→R n Ω ψ(λ + , λ -) dΩ =⇒ ∂I i ∂t = 2 div ∂ψ ∂λ + θ + θ T + + ∂ψ ∂λ - θ -θ T -∇I i (40)
(for i = 1..n) Note that ( 40) is a divergence-based equation such that : 2) is then a 2 × 2 diffusion tensor, whose eigenvalues are :

∂I i ∂t = div (D∇I i ) where D = 2 ∂ψ ∂λ + θ + θ T + + 2 ∂ψ ∂λ - θ -θ T - D ∈ P(
λ 1 = 2 ∂ψ ∂λ + and λ 2 = 2 ∂ψ ∂λ -
associated to these corresponding orthonormal eigenvectors :

u 1 = θ + and u 2 = θ -
It is also worth to mention that computing this gradient descent is done exactly in the same way, when dealing with image domains Ω defined in higher dimensional spaces (Ω ⊂ R p where p > 2) More particularly, the case of 3D volume regularization (p = 3) can be written as :

min I:Ω→R n Ω ψ(λ 1 , λ 2 , λ 3 ) dΩ =⇒ ∂I i ∂t = 2 div ∂ψ ∂λ 1 θ 1 θ T 1 + ∂ψ ∂λ 2 θ 2 θ T 2 + ∂ψ ∂λ 3 θ 3 θ T 3 ∇I i
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.

where ∇G (T,t) and H G (T,t) are respectively the gradient and the Hessian of G (T,t) . It means that

trace (T H G (T,t) ) = - 1 8πt 2 exp - X T T -1 X 4t trace I d - XX T T -1 2t = - 1 8πt 2 exp - X T T -1 X 4t 2 - X T T -1 X 2t = ∂G (T,t) ∂t
And as the convolution is a linear operation, we have (I i (t) * G (T,t) ) = I i0 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.

Appendix C

In this appendix, we develop tensor-driven divergence PDE's into their trace-based counterpart. Most divergence-based regularization PDE's acting on multivalued images have the following form :

∂I i ∂t = div (D∇I i ) (i = 1..n) (41)
where D is a diffusion tensor based only on first order operators. The fact is that D is often computed from the structure tensor G = n j=1 ∇I j ∇I T j and depends mainly on the spatial derivatives I ix and I iy . Intuitively, as the divergence div () = ∂ ∂x + ∂ ∂y is itself a first order derivative operator, we should be able to write [START_REF] Zenzo | A note on the gradient of a multi-image[END_REF] only with first and second spatial derivatives I ix , I iy , I ixx , I ixy and I iyy . Thus, it could be expressed with oriented Laplacians in each image channel I i as well, i.e an expression based on the trace operator ∂Ii ∂t = trace (DH i ).

We want to make the link between the two different diffusion tensors D and T in the divergence-based and trace-based regularization PDE's, in the case when D is not constant :

∂I i ∂t = div (D∇I i ) and ∂I i ∂t = trace (TH i )
As we noticed 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∇I i ) = trace (DH Ii ) + ∇I T i div (D) (42)
where div () is the matrix divergence.

A natural idea is then to decompose the additional term ∇I T i div (D) into oriented Laplacians, expressed with additional diffusion tensors Q in the trace operator.

For this purpose, we will consider that the divergence tensor D is defined at each point X ∈ Ω by

D = f 1 (λ + , λ -) θ + θ T + + f 2 (λ + , λ -) θ -θ T -with f 1/2 : R 2 → R (43)
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 :

∂I i ∂t = div (D∇I i ) with (
i.e we separate the isotropic and anisotropic parts of D, with

α = f 1 (λ + , λ -) -f 2 (λ + , λ -) λ + -λ - and β = λ + f 2 (λ + , λ -) -λ -f 1 (λ + , λ -) λ + -λ - (45)
Indeed, we have

αG + βI d = f 1 -f 2 λ + -λ - (λ + θ + θ T + + λ -θ -θ T -) + λ + f 2 -λ -f 1 λ + -λ - (θ + θ T + + θ -θ T -) = 1 λ + -λ - θ + θ T + (λ + f 1 -λ -f 1 ) + θ -θ T -(λ + f 2 -λ -f 2 ) = f 1 θ + θ T + + f 2 θ -θ T - = D
Here we assumed that λ + = λ -(i.e the structure tensor G is anisotropic). Anyway, if G is isotropic, one generally chooses an isotropic diffusion tensor D too, in the divergence operator of (42), i.e f 1 (λ + , λ -) = f 2 (λ + , λ -). In this case, we choose α = 0 and β = f 1 (λ + , λ -).

This decomposition is useful to rewrite the matrix divergence div (D) into :

div (D) = α div (G) + G∇α + ∇β (46)
and the additional term of the equation ( 42) would be computed as :

∇I T div (D) = trace div (D) ∇I T i = αtrace div (G) ∇I T i (47)
+ trace G∇α∇I T i (48)

+ trace ∇β∇I T i (49)
In the following, we propose to find formal expressions of ( 47), ( 48) and ( 49).

• First, remember that the structure tensor G is defined as :

• We finally have to compute ∇α and ∇β, in the expression ( 48) and ( 49). This can be done by the decomposition :

∇α
and as the λ ± , eigenvalues of the structure tensor G, depends on the I jx and I jy : Then,

∇λ ± = n j=1 2H j θ ± θ T ± ∇I j (52)
We can replace (52) into the expressions of [START_REF] Granlund | Signal Processing for Computer Vision[END_REF], in order to find the spatial gradients of α and β :

               ∇α = n j=1
• The final step consists in putting together the equations ( 50) and [START_REF] Hunt | TTe application of constrained least squares estimation to image restoration by digital computer[END_REF], in order to express the additional term ∇I T i div (D) in the PDE [START_REF] Florack | Image Structure[END_REF].

∇I T i div (D) = n j=1 trace H j P ij (55)
where the P ij are the following 2 × 2 matrices :

P ij = α
Note that the indices i, j in the notation P ij do not designate the coefficients of a matrix P, but the parameters of the family consisting of n 2 matrices P ij (each of them is a 2 × 2 matrix).

The matrices P ii are symmetric, but generally not the P ij (where i = j), since the gradients ∇I i and ∇I j are not aligned in the general case. Yet, we want to express the equation ( 55) 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 :

trace (AH) = trace A + A T 2 H
where (A + A T )/2 is a 2 × 2 symmetric matrix (the symmetric part of A).

Thus, we define the symmetric matrices Q ij , corresponding to the symmetric parts of the P ij :

Q ij = P ij + P ij T 2 (57)
and we have :

∇I T i div (D) = n j=1 trace H j Q ij
Finally, the divergence-based PDE (42) can be written as :

div (D∇I i ) = n j=1 trace (δ ij D + Q ij )H j ( 58
)
where δ ij is the Kronecker's symbol :

δ ij = 0 if i = j 1 if i = j
This makes the link between divergence PDE's and sums of atomic trace-based PDE's. A direct geometric interpretation of (58) is not direct anyway.

  (a) Initial image I (t=0) (b) t = 50 (c) t = 250 (d) t = 1000 (e) t = 3000

Figure 1 :

 1 Figure 1: Nonlinear regularization PDE's and the notion of anisotropic scale-space.

  (a) Red channel R (b) Green channel G (c) Blue channel B (d) Color image (R, G, B) (e) Lightness (scalar) image L *

Figure 2 :

 2 Figure 2: Using lightness L * to detect geometry of a color image fail for iso-lightness contours.

Figure 3 :

 3 Figure 3: Comparing possible vector variation norms N , N -and N + for a synthetic color image.

  (a) Color photograph (small portion 60 × 40) (b) N = √ λ+ (c) N-= √ λ+λ-(d) N+ = √ λ+ + λ-

Figure 4 :

 4 Figure 4: Comparing possible vector variation norms N , N -and N + for a real color image.

  Detecting color edges with the norm N+ = √ λ+ + λ-

Figure 5 :

 5 Figure 5: Using a vector variation norm for color edge detection.

Figure 7 :

 7 Figure 7: Example of image denoising by surface area minimization.

Figure 8 :

 8 Figure 8: Principle of oriented Laplacians : Two 1D smoothing are done along adapted directions.

Figure 9 :

 9 Figure 9: Trace-based PDE's (15) viewed as convolutions by oriented 2D Gaussian kernels.

Figure 10 :

 10 Figure 10: Trace-based PDE's (15) with non-constant diffusion tensor fields T.

  (a) Noisy synthetic color image (b) Applying trace-based PDE (15), with p1 = 0.5, p2 = 1.2.(c) Applying trace-based PDE[START_REF] Bertalmio | Image inpainting[END_REF], with p1 = 0.9, p2 = 1.2.(d) Applying our constrained PDE[START_REF] Chan | Variational restoration of non-flat image features : Models and algorithms[END_REF], with p1 = 0.5, p2 = 1.2.

Figure 11 :

 11 Figure 11: Problems encountered when using trace-based PDE's (15) on curved image structures.

  (a) Image of a fingerprint (b) Applying trace-based PDE (15), with p1 = 0.5, p2 = 1.2.(c) Applying our constrained PDE[START_REF] Chan | Variational restoration of non-flat image features : Models and algorithms[END_REF], with p1 = 0.5, p2 = 1.2.

Figure 12 :

 12 Figure 12: Comparisons between trace-based PDE's (15) and our new curvature-preserving PDE's (29) on a real image.

  (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).

Figure 13 :

 13 Figure 13: Integral curve C X of a vector field w : Ω → R 2 .

  Regularization using the LIC-based scheme (stopped at t = 100).

Figure 14 :

 14 Figure 14: Comparisons between classical explicit PDE schemes, and LIC-based implementation of the PDE (29).

  illustrates a simple case where one wants to remove 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. Result of Noisy color image (left), denoised image (middle), zoom on the eye (right).

Figure 15 :

 15 Figure 15: Denoising of the color image "Baboon" corrupted with artificial Gaussian noise.

Figure 16 :

 16 Figure 16: Denoising of the color image "Tunisia Desert", containing watered effects.

Figure 17 :

 17 Figure 17: Denoising of the color image "Lolotte", real color photograph.

Figure 18 :

 18 Figure 18: Removing blocs artifacts from the JPEG compressed color image "Baby".

Figure 19 :

 19 Figure 19: Removing quantization noise in the quantized color image "Penguin".

  Corrupted color image (left), inpainted image (right). Details are shown on the bottom row.

Figure 20 :

 20 Figure 20: Removing undesirable text in color image by an inpainting method.

Figure 21 :

 21 Figure 21: Removing a real object in a color image by an inpainting method.

Figure 22 :

 22 Figure 22: Reconstruction of 50% of a color image by an inpainting method.

Figure 23 :

 23 Figure 23: Reconstruction of 50% of a color image by an inpainting method.

Figure 24 :

 24 Figure 24: Removing a real object in a color image by an inpainting method.

  (a) Thumbnail image. (b) Details from the image resized by nearest-neighbor interpolation. (c) Details from the image resized by linear interpolation. (c) Details from the image resized by bicubic interpolation. (d) Details from the image resized by a non-linear regularization PDE.

Figure 25 :

 25 Figure 25: Comparisons of image resizing methods, Nearest-neighbor (first row), Linear (second row), Bicubic (third row) and PDE-based (last row) interpolations.

Figure 26 :

 26 Figure 26: Visualization of a 2D vector field : using arrows (left), after 5 PDE iter. (middle), after 15 PDE iter. (right).

∂

  (I i0 * G (T,t) ) ∂t = I i0 * ∂G (T,t) ∂t = I i0 * trace (T H G (T,t) ) = trace T H Ii 0 * G (T,t)as well as lim

  t→0

f 1

 1 (λ + , λ -) = 2 ∂ψ ∂λ + f 2 (λ + , λ -) = 2 ∂ψ ∂λ -In order to develop the additional diffusion term ∇I T i div (D) in the equation (42), we propose to write D as a linear combination of G and I d :D = α(λ + , λ -)G + β(λ + , λ -)I d

2 I

 2 jx I jxx + I jx I jyy + I jy I jxy I jx I jxy + I jy I jxx + 2 I jy I jyy = n j=1 I jx (I jxx + I jyy ) I jy (I jxx + I jyy ) + I jx I jxx + I jy I jxy I jx I jxy + I jy I jyy = n j=1 ∆I j ∇I j + H j ∇I jwhere ∆I j and H j are respectively the Laplacian and the Hessian of the image component I j .Then, we can write the expression 47 as :αtrace div (G) ∇I i T = n j=1αtrace H j ∇I i T ∇I j I d + ∇I j ∇I i T

  , we derived eigenvalues of a structure tensor G, with respect to the spatial image derivatives. We ended up with the following relation : ± θ T ± ∇I j

  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 in[START_REF] Tschumperlé | Deriche Vector-Valued Image Regularization with PDE's : A Common Framework for Different Applications[END_REF][START_REF] Tschumperlé | Deriche Vector-Valued Image Regularization with PDE's : A Common Framework for Different Applications[END_REF] unifies a whole range of previously proposed vector-valued regularization algorithms (variational and divergence based PDE's) 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. Though the geometric interpretation of the overall sum of trace equations is not direct, it is interesting to see that additional diffusion tensors Q ij clearly appear in the trace expressions, and contribute to modify the inner tensor D, which is finally not representative of the smoothing behavior. More generally, tensors appearing in traces and divergences generally lead to different smoothing behaviors.

					λ+-λ-
	and	f 1 =	∂ψ ∂λ +	and f 2 =	∂ψ ∂λ -

  ∇I T i ∇I j I d

	+ 2 + 2 (α + ∂α ∂λ +	θ + θ T + + ∂β ∂λ + )θ + θ T ∂α ∂λ -+ + (α + θ -θ T -∇I j ∇I T i G ∂β ∂λ -)θ -θ T -∇I j ∇I T i

Appendix B

In this appendix, we demonstrate that the solution of the generic trace-based PDE :

is the convolution of the image

by an oriented Gaussian kernel G (T,t) defined as :

To demonstrate this, we simply derive the kernel expression G (T,t) in time and in space :