Laurent Condat

Discrete Total Variation: New Definition and Minimization

Keywords: total variation, variational image processing, coarea formula, finite-difference schemes AMS subject classifications. 68Q25, 68R10, 68U05

We propose a new definition for the gradient field of a discrete image, defined on a twice finer grid. The differentiation process from the image to its gradient field is viewed as the inverse operation of linear integration, and the proposed mapping is nonlinear. Then, we define the total variation of an image as the 1 norm of its gradient field amplitude. This new definition of the total variation yields sharp edges and has better isotropy than the classical definition.

1. Introduction. In their seminal paper, Rudin, Osher, and Fatemi [START_REF] Rudin | Nonlinear total variation based noise removal algorithms[END_REF] introduced the total variation (TV) regularization functional for imaging problems. Since then, a variety of papers has demonstrated the effectiveness of TV minimization to recover sharp images, by preserving strong discontinuities, while removing noise and other artifacts [START_REF] Chambolle | An introduction to total variation for image analysis[END_REF][START_REF] Chambolle | A first-order primal-dual algorithm for convex problems with applications to imaging[END_REF][START_REF] Burger | A guide to the TV zoo[END_REF]. TV minimization also appears in clustering and segmentation problems, by virtue of the coarea formula [START_REF] Goldstein | Geometric applications of the split Bregman method: Segmentation and surface reconstruction[END_REF][START_REF] Chambolle | A convex approach to minimal partitions[END_REF]. The TV can be defined in other settings than image processing, for instance on graphs [START_REF] Couprie | Dual constrained TV-based regularization on graphs[END_REF]. Numerical minimization of the TV has long been challenging, but recent advances in large-scale convex nonsmooth optimization, with efficient primal-dual splitting schemes and alternating directions methods, have made the implementation of TV minimization relatively easy and efficient [START_REF] Chambolle | A first-order primal-dual algorithm for convex problems with applications to imaging[END_REF][START_REF] Goldstein | The split Bregman method for L1-regularized problems[END_REF][START_REF] Ng | Solving constrained total-variation image restoration and reconstruction problems via alternating direction methods[END_REF][START_REF] Afonso | Fast image recovery using variable splitting and constrained optimization[END_REF][START_REF] Combettes | Dualization of signal recovery problems[END_REF][START_REF] Zhang | A unified primal-dual algorithm framework based on Bregman iteration[END_REF][START_REF] Briceño-Arias | A monotone+skew splitting model for composite monotone inclusions in duality[END_REF][START_REF] Briceño-Arias | Proximal algorithms for multicomponent image recovery problems[END_REF][START_REF] Combettes | Proximity for sums of composite functions[END_REF][START_REF]A splitting algorithm for dual monotone inclusions involving cocoercive operators[END_REF][START_REF] Condat | A generic proximal algorithm for convex optimization-Application to total variation minimization[END_REF][START_REF] Combettes | A forward-backward view of some primal-dual optimization methods in image recovery[END_REF][START_REF] Komodakis | Playing with duality: An overview of recent primal-dual approaches for solving large-scale optimization problems[END_REF]. Yet, the rigorous definition of the TV for discrete images has received little attention. For continuously defined two-dimensional functions, the TV is simply the L 1 norm of the gradient amplitude. But for discrete images, it is a nontrivial task to properly define the gradient using finite differences, as is well known in the community of computer graphics and visualization [START_REF] Alim | Gradient estimation revitalized[END_REF][START_REF] Hossain | Towards high quality gradient estimation on regular lattices[END_REF]. The classical, so called "isotropic" definition of the discrete TV, is actually far from being isotropic, but it performs reasonably well in practice. In this paper, we propose a new definition of the discrete TV, which corrects some drawbacks of the classical definition and yields sharper edges and structures. The key idea is to associate, in a nonlinear way, an image with a gradient field on a twice finer grid. The TV of the image is then simply the 1 norm of this gradient field amplitude.

In section 2, we review the classical definitions of the discrete TV and their properties. In section 3, we introduce our new definition of the TV in the dual domain and in section 4, we study the equivalent formulation in the primal domain. An algorithm to solve problems regularized with the proposed TV is presented in section 5. The good performances of the proposed TV on some test imaging problems are demonstrated in section 6.

2. Classical definitions of the discrete TV and their properties. A function s(t 1 , t 2 ) defined in the plane R 2 , under some regularity assumptions, has a gradient field ∇s(t 1 , t 2 ) = ∂s ∂t 1 (t 1 , t 2 ), ∂s ∂t 2 (t 1 , t 2 ) , defined in R 2 as well. We can then define the TV of s as the L 1,2 norm of the gradient: TV(s) = R 2 |∇s(t 1 , t 2 )| dt 1 dt 2 , where |(a, b)| is a shorthand notation for the 2-norm √ a 2 + b 2 . The TV has the desirable property of being isotropic, or rotation-invariant: a rotation of s in the plane does not change the value of its TV.

A (grayscale) discrete image x of size N 1 × N 2 has its pixel values x[n 1 , n 2 ] defined at the locations (n 1 , n 2 ) in the domain Ω = {1, . . . , N 1 } × {1, . . . , N 2 }, where n 1 and n 2 are the row and column indices, respectively, and the pixel with index (1, 1) is at the top left image corner. The pixel values are supposed to lie between 0 (black) and 1 (white). The challenge is then to define the discrete TV of x, using only its pixel values, while retaining the mathematical properties of the continuous TV. The so-called anisotropic TV is defined as [START_REF] Rudin | Nonlinear total variation based noise removal algorithms[END_REF] TV a (x) =

N 1 n 1 =1 N 2 n 2 =1 x[n 1 + 1, n 2 ] -x[n 1 , n 2 ] + x[n 1 , n 2 + 1] -x[n 1 , n 2 ] ,
assuming Neumann (symmetric) boundary conditions: a finite difference across a boundary, like

x[N 1 + 1, n 2 ] -x[N 1 , n 2 ]
, is assumed to be zero. The anisotropic TV is well known to be a poor definition of the discrete TV, as it yields metrication artifacts: its minimization favors horizontal and vertical structures, because oblique structures make the TV value larger as it should be. Therefore, one usually uses the so-called isotropic TV defined as

(2)

TV i (x) = N 1 n 1 =1 N 2 n 2 =1 x[n 1 + 1, n 2 ] -x[n 1 , n 2 ] 2 + x[n 1 , n 2 + 1] -x[n 1 , n 2 ] 2 ,
using Neumann boundary conditions as well.

It is hard to quantify the isotropy of a functional like the TV, since the grid Z 2 is not isotropic and there is no unique way of defining the rotation of a discrete image. However, it is natural to require, at least, that after a rotation of ±90 o , or a horizontal or vertical flip, the TV of the image remains unchanged. It turns out that this is not the case with the isotropic TV, with a change factor as large as √ 2 after a horizontal flip, see in Table 1 the TV of an edge at +45 o and at -45 o . In spite of this significant drawback, the isotropic TV is widely used, for its simplicity. We can note that a straightforward way to restore the four-fold symmetry is to define the TV as the average of TV i applied to the image rotated by 0 o , 90 o , -90 o , 180 o . But the drawbacks of TV i , stressed below, would be maintained, like the tendency to blur oblique edges and the too low value for an isolated pixel or a checkerboard.

An attempt to define a more isotropic TV has been made with the upwind TV [START_REF] Chambolle | An upwind finite-difference method for total variationbased image smoothing[END_REF], defined as

(3) TV u (x) = N 1 n 1 =1 N 2 n 2 =1 x[n 1 , n 2 ] -x[n 1 + 1, n 2 ] 2 + + x[n 1 , n 2 ] -x[n 1 -1, n 2 ] 2 + + x[n 1 , n 2 ] -x[n 1 , n 2 + 1] 2 + + x[n 1 , n 2 ] -x[n 1 , n 2 -1] 2 +
, where (a) + means max(x, 0). The upwind TV is indeed more isotropic and produces sharp oblique edges, but as shown below, it is not invariant by taking the image negative, i.e. replacing the image x by 1 -x. Since TV u (x) = TV u (1 -x) = TV u (-x), the upwind TV is not a seminorm, contrary to the other forms considered in this paper. In practice, it penalizes correctly small dark structures over a light background, but not the opposite, see the striking example in Figure 10 (e).

(I) (II) (IIf) (III) (IV) (V) (Vn) (VI) (VII) (VIII) (IX) (X)
Figure 1. Some patterns, for which we report the value of the TV in Table 1. Black and white correspond to 0 and 1, respectively. In (III), the transition goes through the levels 0, 1/8, 7/8, 1. In (IV), the transition goes through the levels 0, 1/2, 1. In (VII), the transition goes through the levels 0, 1/2, 1, 1/2, 0.

Another formulation of the discrete TV, called "Shannon Total Variation" was proposed recently [START_REF] Abergel | The Shannon total variation[END_REF], at the time the present paper was finalized; so, this formulation, which has good isotropy properties, is not included in our comparisons. It aims at estimating the continuous TV of the Shannon interpolate of the image, by using a Riemann sum approximation of the corresponding integral. This way, aliasing is removed from the images, at the price of slightly more blurred edges.

To evaluate the different definitions of the discrete TV, we consider typical patterns of size N ×N , depicted in Figure 1, and we report the corresponding value of the TV in Table 1, when N is large, i.e. ignoring the influence of the image boundaries. For some patterns, we consider its horizontally flipped version, denoted by a 'f', see patterns (II) and (IIf) in Figure 1, and its negative version, denoted by a 'n', see patterns (V) and (Vn). In Table 1, the value is in green if it is an appropriate value for this case, and in red if not. In this respect, some considerations must be reported. An isolated pixel, like in patterns (VIII) or (VIIIn), can be viewed as the discretization by cell-averaging, i.e.

x[n 1 , n 2 ] = n 1 +1/2 n 1 -1/2 n 2 +1/2 n 2 -1/2 s(t 1 , t 2 )dt 1 dt 2 ,
of the indicator function (1 inside, 0 outside) s(t) of a square of size one pixel. According to the coarea formula, the continuous TV of the indicator function of a set is equal to the perimeter of that set. So, it is natural to ask the TV value in pattern (VIII) to be equal to 4. The isotropic TV Asymptotic value, when the image is of size N × N and N → +∞, of the TV, for the examples depicted in Figure 1. A 'f ' means a horizontal flip and a 'n' means taking the image negative. TVa, TVi, TVu, TVp are the anisotropic, isotropic, upwind, proposed TV, defined in (1), ( 2), ( 3), [START_REF] Goldstein | The split Bregman method for L1-regularized problems[END_REF], respectively.

TV a TV i TV u TV p (I) N N N N (II) 2N √ 2N √ 2N 2N (IIf) 2N 2N √ 2N 2N (III) 2N √ 2N √ 2N √ 2N (IIIf) 2N ( √ 37 + 1)N/4 √ 2N √ 2N (IV) 2N √ 2N √ 2N √ 2N (IVf) 2N (1 + 1/ √ 2)N √ 2N √ 2N (V) 2N 2N √ 2N 2N (Vn) 2N 2N 2N 2N (VI) 4N 2 √ 2N 2N 4N (VIf) 4N (2 + √ 2)N 2N 4N (VIn) 4N 2 √ 2N 2 √ 2N 4N (VII) 4N 2 √ 2N ( √ 2 + 1)N 2 √ 2N (VIIf) 4N (3/ √ 2 + 1)N ( √ 2 + 1)N 2 √ 2N (VIIn) 4N 2 √ 2N 2 √ 2N 2 √ 2N (VIII) 4 2 + √ 2 2 4 (VIIIn) 4 2 + √ 2 4 4 (IX) N 2 N 2 N 2 / √ 2 N 2 (X) 2N 2 √ 2N 2 N 2 2N 2
and upwind TV take too small values. This is a serious drawback, since they do not penalize noise as much as they should, and penalizing noise is the most important property of a functional used to regularize ill-posed problems. For the checkerboard (X), it is natural to expect a value of 2N 2 . It is important that this value is not lower, because an inverse problem like demosaicking consists in demultiplexing luminance information and chrominance information modulated at this highest frequency [START_REF] Alleysson | Linear demosaicing inspired by the human visual system[END_REF][START_REF] Condat | Joint demosaicking and denoising by total variation minimization[END_REF]. Interpolation on a quincunx grid also requires penalizing the checkerboard sufficiently. The isotropic TV gives a value of √ 2N 2 , which is too small, and the upwind TV gives an even smaller value of N 2 . Then, an important property of the TV is to be convex and one-homogeneous, so that the TV of a sum of images is less or equal than the sum of their TV. Consequently, viewing the checkerboard as a sum of diagonal lines, like the one in (VI), disposed at every two pixels, the TV of the diagonal line (VI) cannot be lower than 4N . That is, the lower value of 2 √ 2N , achieved by the isotropic TV, is not compatible with the value of 2N 2 for the checkerboard and with convexity of the TV. We can notice that the line in (VI) cannot be explained as the discretization by cell-averaging of a continuously defined diagonal ridge. So, it is coherent that its jaggy nature is penalized. By contrast, the pattern in (VII) can be viewed as the discretization by cell-averaging of a diagonal ridge, depicted in Figure 2 (c). So, a TV value of 2 √ 2N is appropriate for this case. Further on, the line in (VI) can be viewed as the difference of two edges like in (II), one of (a) (b) (c) which shifted by one pixel. So, by convexity, the value of the TV for the edge in (II) cannot be lower than 2N . The value of √ 2N , we could hope for by viewing (II) as a diagonal edge discretized by point sampling, is not accessible. Again, after a small blur, the discrete edges in (III) and (IV) become compatible with a diagonal edge discretized by cell-averaging, see the edges in Figure 2 (a) and (b), respectively. So, the expected value of the TV is √ 2N in these cases. It is true that a TV value of √ 2N would be nice for the binary edge (II), especially for partitioning applications [START_REF] Chambolle | A convex approach to minimal partitions[END_REF], and that the isotropic TV achieves this value, but the price to pay with the isotropic TV is a higher value of 2N for the flipped case (IIf), which does not decrease much by blurring the edge to (IIIf) or (IVf). Therefore, minimizing the isotropic TV yields nice binary edges at the diagonal orientation like in (II), but significantly blurred edges for the opposite orientation, as can be observed in Figure 8 We can mention, mainly in the literature of computational fluid or solid mechanics, the use of staggered grid discretizations of partial differential equations, or marker and cell method [START_REF] Harlow | Numerical calculation of time-dependent viscous incompressible flow of fluid with free surface[END_REF], wherein different variables, like the pressure and velocity, are located at different positions on the grid, i.e. at cell centers or at cell edges. This idea is also applied in so-called mimetic finite difference methods [START_REF] Bochev | Principles of mimetic discretizations of differential operators[END_REF][START_REF] Bazan | Mimetic finite difference methods in image processing[END_REF]. Transposed to the present context, pixel values are located at the pixel centers, whereas a finite difference like x[n 1 + 1, n 2 ] -x[n 1 , n 2 ] is viewed as the vertical component of the gradient at the spatial position (n 1 + 1 2 , n 2 ), i.e. at an edge between two pixels [START_REF] Garamendi | Box relaxation schemes in staggered discretizations for the dual formulation of total variation minimization[END_REF]. This interpretation is insightful, but it does not specify how to define the norm of the gradient. The proposed approach is different from this framework in two respects. First, we define the image gradient field not only at the pixel edges, but also at the pixel centers. Second, a finite difference like x[n 1 + 1, n 2 ] -x[n 1 , n 2 ] is not viewed as an estimate of a partial derivative, but as its local integral; we develop this interpretation in section 4.

3. Proposed discrete TV: dual formulation. It is well known that in the continuous domain, the TV of a function s can be defined by duality as (4)

TV(s) = sup s, -div(u) : u ∈ C 1 c (R 2 , R 2 ), |u(t)| ≤ 1 (∀t ∈ R 2 ) ,
where

C 1 c (R 2 , R 2 )
is the set of continuously differentiable functions from R 2 to R 2 with compact support and div is the divergence operator. So, the dual variable u has its amplitude bounded by one everywhere.

In the discrete domain, the TV can be defined by duality as well. First, let us define the discrete operator D, which maps an image x ∈ R N 1 ×N 2 to the vector field Dx ∈ (R 2 ) N 1 ×N 2 made of forward finite differences of x; that is,

(Dx) 1 [n 1 , n 2 ] = x[n 1 + 1, n 2 ] -x[n 1 , n 2 ], (5)
(Dx) 2 [n 1 , n 2 ] = x[n 1 , n 2 + 1] -x[n 1 , n 2 ], (6)
for every (n 1 , n 2 ) ∈ Ω, with Neumann boundary conditions. Note that for ease of implementation, it is convenient to have all images and vector fields of same size N 1 × N 2 , indexed by (n 1 , n 2 ) ∈ Ω, keeping in mind that for some of them, the last row or column is made of dummy values equal to zero, which are constant and should not be viewed as variables; for instance,

(Dx) 1 [N 1 , n 2 ] = (Dx) 2 [n 1 , N 2 ] = 0, for every (n 1 , n 2 ) ∈ Ω. So, TV i (x) = Dx 1,2 , where the 1,2 norm is the sum over the indices n 1 , n 2 of the 2-norm |(Dx)[n 1 , n 2 ]|.
Then, the isotropic TV of an image x can be defined by duality as

(7) TV i (x) = max u∈(R 2 ) N 1 ×N 2 Dx, u : |u[n 1 , n 2 ]| ≤ 1, ∀(n 1 , n 2 ) ∈ Ω ,
with the usual Euclidean inner product.

The scalar dual variables

u 1 [n 1 , n 2 ] and u 2 [n 1 , n 2 ], like the finite differences (Dx) 1 [n 1 , n 2 ] and (Dx) 2 [n 1 , n 2 ]
, can be viewed as located at the points (n 1 + 1 2 , n 2 ) and (n 1 , n 2 + 1 2 ), respectively. So, the anisotropy of the isotropic TV can be explained by the fact that these variables, which are combined in the constraint |u[n 1 , n 2 ]| ≤ 1, are located at different positions. We propose to correct this half-pixel shift by interpolation: we look for the dual images u 1 and u 2 , whose values

u 1 [n 1 , n 2 ] and u 2 [n 1 , n 2 ] are located at the pixel edges (n 1 + 1 2 , n 2 ) and (n 1 , n 2 + 1 2 )
, respectively, such that, when interpolated, the constraint |u[n 1 , n 2 ]| ≤ 1 is satisfied both at pixel centers and at pixel edges. So, the proposed TV, denoted TV p , is defined in the dual domain as

TV p (x) = max u∈(R 2 ) N 1 ×N 2
Dx, u :

|(L u)[n 1 , n 2 ]| ≤ 1, |(L ↔ u)[n 1 , n 2 ]| ≤ 1, |(L • u)[n 1 , n 2 ]| ≤ 1, ∀(n 1 , n 2 ) ∈ Ω , (8)
where the three operators L , L ↔ , L • interpolate bilinearly the image pair u = (u 1 , u 2 ) on the grids

(n 1 + 1 2 , n 2 ), (n 1 , n 2 + 1 2 ), (n 1 , n 2 ), for (n 1 , n 2 ) ∈ Ω, respectively. That is, (L u) 1 [n 1 , n 2 ] = u 1 [n 1 , n 2 ], (9)
(L u) 2 [n 1 , n 2 ] = (u 2 [n 1 , n 2 ] + u 2 [n 1 , n 2 -1] + u 2 [n 1 + 1, n 2 ] + u 2 [n 1 + 1, n 2 -1])/4, (10)
(L ↔ u) 1 [n 1 , n 2 ] = (u 1 [n 1 , n 2 ] + u 1 [n 1 -1, n 2 ] + u 1 [n 1 , n 2 + 1] + u 1 [n 1 -1, n 2 + 1])/4, (11)
(L ↔ u) 2 [n 1 , n 2 ] = u 2 [n 1 , n 2 ], (12)
(L • u) 1 [n 1 , n 2 ] = (u 1 [n 1 , n 2 ] + u 1 [n 1 -1, n 2 ])/2, ( 13
) (L • u) 2 [n 1 , n 2 ] = (u 2 [n 1 , n 2 ] + u 2 [n 1 , n 2 -1])/2, (14) for every (n 1 , n 2 ) ∈ Ω, replacing the dummy values u 1 [0, n 2 ], u 2 [n 1 , 0], u 1 [N 1 , n 2 ], u 2 [n 1 , N 1 ], (L u) 1 [N 1 , n 2 ], (L u) 2 [N 1 , n 2 ], (L ↔ u) 1 [n 1 , N 2 ], (L ↔ u) 2 [n 1 , N 2 ] by zero.
Thus, we mimic the continuous definition [START_REF] Burger | A guide to the TV zoo[END_REF], where the dual variable is bounded everywhere, by imposing that it is bounded on a grid three times more dense than the pixel grid. Actually, for the dual variable to be bounded everywhere after bilinear interpolation, the fourth lattice of pixel corners (n 1 + 1 2 , n 2 + 1 2 ) must be added; that is, we can define a variant of the proposed approach, in which we add to the constraint set in [START_REF] Goldstein | The split Bregman method for L1-regularized problems[END_REF], the additional constraints

|(L + u)[n 1 , n 2 ]| ≤ 1,
where the operator L + interpolates bilinearly the image pair u on the grid

(n 1 + 1 2 , n 2 + 1 2 ): (L + u) 1 [n 1 , n 2 ] = (u 1 [n 1 , n 2 ] + u 1 [n 1 , n 2 + 1])/2, (15)
(L + u) 2 [n 1 , n 2 ] = (u 2 [n 1 , n 2 ] + u 2 [n 1 + 1, n 2 ])/2. ( 16
)
In the Matlab code accompanying the paper, available on the author's webpage, this variant is implemented, as well. The author observed empirically that this variant, in general, brings very minor changes in the images, which are not worth the extra computational burden. That is why, in the rest of the paper, we focus on the version with the dual variables and the gradient defined on three lattices, and not on this variant with four lattices.

Our definition of the discrete TV, using interpolation in the dual domain, is not new: it was proposed in [START_REF] Hintermüller | Functional-analytic and numerical issues in splitting methods for total variation-based image reconstruction[END_REF] and called staggered grid discretization of the TV. With the isotropic TV, the projection of the image pair u onto the l ∞,2 norm ball, which amounts to simple pixelwise shrinkage, can be used. But using the same algorithms with the proposed TV requires projecting u onto the set {u :

L u ∞,2 ≤ 1, L ↔ u ∞,2 ≤ 1, L • u ∞,2 ≤ 1}.
There is no closed form for this projection. We emphasize that in [START_REF] Hintermüller | Functional-analytic and numerical issues in splitting methods for total variation-based image reconstruction[END_REF], and certainly in other papers using this dual staggered grid discretization, this projection is not implemented, and is replaced by an approximate shrinkage, see [START_REF] Hintermüller | Functional-analytic and numerical issues in splitting methods for total variation-based image reconstruction[END_REF]Eq. (64)]. This operation is not a projection onto the set above, since it is not guaranteed to yield an image pair satisfying the bound constraints, and it is not a firmly nonexpansive operator [START_REF] Bauschke | Convex Analysis and Monotone Operator Theory in Hilbert Spaces[END_REF]; this means that the convergence guarantees of usual iterative fixed-point algorithms are lost, and that if convergence occurs, there is no way to characterize the obtained solution, which depends on the algorithm, the initial conditions, and the parameters. By contrast, we will propose a generic splitting algorithm, with proved convergence to exact solutions of problems involving the proposed TV, in section 5.

4. Proposed discrete TV: primal formulation. We have defined the proposed TV implicitly in [START_REF] Goldstein | The split Bregman method for L1-regularized problems[END_REF], as the optimal value of an optimization problem, expressed in terms of the dual image pair u. In the frame of the Fenchel-Rockafellar duality [START_REF] Bauschke | Convex Analysis and Monotone Operator Theory in Hilbert Spaces[END_REF], we can define the proposed TV as the optimal value of an equivalent optimization problem, expressed in terms of what we will consider as the gradient field of the image.

Proposition 1. Given an image x, the proposed TV has the following primal formulation, equivalent to the dual formulation (8):

(17) TV p (x) = min v ,v↔,v•∈(R 2 ) N 1 ×N 2 v 1,2 + v ↔ 1,2 + v • 1,2 : L * v +L * ↔ v ↔ +L * • v • = Dx ,
where • * denotes the adjoint operator.

Before proving Proposition 1, we give more compact forms of the primal and dual definitions of the proposed TV. For this, let us define the linear operator L as the concatenation of L , L ↔ , L • , and the ∞,∞,2 norm • ∞,∞,2 of a field as the maximum over the three components and the pixels of the 2-norm of its vectors. Then we can rewrite (8) as [START_REF] Combettes | A forward-backward view of some primal-dual optimization methods in image recovery[END_REF] TV p (x) = max

u∈(R 2 ) N 1 ×N 2 Dx, u : Lu ∞,∞,2 ≤ 1 .
Let the vector field v be the concatenation of the three vector fields v , v ↔ , and v • , which appear in [START_REF] Condat | A generic proximal algorithm for convex optimization-Application to total variation minimization[END_REF]. Let the 1,1,2 norm of v be the sum of the 1,2 norm of its three components

v , v ↔ , v • . We have L * v = L * v + L * ↔ v ↔ + L * • v • .
Then we can rewrite [START_REF] Condat | A generic proximal algorithm for convex optimization-Application to total variation minimization[END_REF] as

(19) TV p (x) = min v∈((R 2 ) N 1 ×N 2 ) 3 v 1,1,2 : L * v = Dx .
Proof of Proposition 1. A (primal) convex optimization problem of the form: minimize v {F (L * v) + G(v)}, for two convex, lower semicontinuous functions F and G and a linear operator L * , has a Fenchel-Rockafellar dual problem of the form: maximize u {-F * (u) -G * (-Lu)}, where F * and G * are the Legendre-Fenchel conjugates of F and G, respectively [START_REF] Bauschke | Convex Analysis and Monotone Operator Theory in Hilbert Spaces[END_REF]. Moreover, strong duality holds, and the primal and dual problems have the same optimal value; that is, if a minimizer v of the primal problem and a maximizer û of the dual problem exist, we have

-F * (û) -G * (-Lû) = F (L * v) + G(v).
In our case, F is the convex indicator function of the set {v : L * v = Dx}; that is, the function which maps its variable to 0 if it belongs to this set, to +∞ else. G is the 1,1,2 norm. Then it is well known that F * maps u to u, Dx and that the Legendre-Fenchel conjugate of the 1,2 norm is the convex indicator function of the ∞,2 norm ball [START_REF] Chambolle | An introduction to total variation for image analysis[END_REF][START_REF] Chambolle | A first-order primal-dual algorithm for convex problems with applications to imaging[END_REF]. So, we see that [START_REF] Goldstein | The split Bregman method for L1-regularized problems[END_REF], up to an unimportant change of sign of u, is indeed the dual problem associated to the primal problem [START_REF] Condat | A generic proximal algorithm for convex optimization-Application to total variation minimization[END_REF]; they share the same optimal value, which is TV p (x).

In the following, given an image x, we denote by v , v ↔ , and v • , the vector fields solution to [START_REF] Condat | A generic proximal algorithm for convex optimization-Application to total variation minimization[END_REF] (or any solution if it is not unique). We denote by v the vector field, which is the concatenation of v , v ↔ , and v • . So, for every

(n 1 , n 2 ) ∈ Ω, its elements v [n 1 , n 2 ], v ↔ [n 1 , n 2 ], v • [n 1 , n 2 ] are vectors of R 2 , located at the positions (n 1 + 1 2 , n 2 ), (n 1 , n 2 + 1 2
), (n 1 , n 2 ), respectively. Then we call v the gradient field of x. Thus, the proposed TV is the 1,2 norm of the gradient field v associated to the image x, solution to [START_REF] Condat | A generic proximal algorithm for convex optimization-Application to total variation minimization[END_REF] and defined on a grid three times more dense than the one of x. The mapping from x to its gradient field v is nonlinear and implicit: given x, one has to solve the optimization problem [START_REF] Condat | A generic proximal algorithm for convex optimization-Application to total variation minimization[END_REF] to obtain its gradient field and the value TV p (x). We can notice that the feasible set in ( 17) is nonempty, since the constraint is satisfied by the vector field defined by

v ,1 = (Dx) 1 , v ,2 = 0 (20) v ↔,1 = 0, v ↔,2 = (Dx) 2 (21) v •,1 = 0, v •,2 = 0. ( 22
)
This vector field has a 1,2 norm equal to (Dx) 1 1 + (Dx) 2 1 , which is exactly TV a (x), the value of the anisotropic TV of x. Therefore, we have the property: for every image x, [START_REF] Abergel | The Shannon total variation[END_REF] TV p (x) ≤ TV a (x).

Further on, we have

(L * v + L * ↔ v ↔ + L * • v • ) 1 [n 1 , n 2 ] = v ,1 [n 1 , n 2 ] + (v ↔,1 [n 1 , n 2 ] + v ↔,1 [n 1 , n 2 -1] + v ↔,1 [n 1 + 1, n 2 ] + v ↔,1 [n 1 + 1, n 2 -1])/4 + (24) (v •,1 [n 1 , n 2 ] + v •,1 [n 1 + 1, n 2 ])/2, (L * v + L * ↔ v ↔ + L * • v • ) 2 [n 1 , n 2 ] = v ↔,2 [n 1 , n 2 ] + (v ,2 [n 1 , n 2 ] + v ,2 [n 1 , n 2 + 1] + v ,2 [n 1 -1, n 2 ] + v ,2 [n 1 -1, n 2 + 1])/4 + (25) (v •,2 [n 1 , n 2 ] + v •,2 [n 1 , n 2 + 1])/2,
using, again, zero boundary conditions. So, the quantity

(L * v + L * ↔ v ↔ + L * • v • ) 1 [n 1 , n 2 ]
is the sum of the vertical part of the elements of the vector field v falling into the square

[n 1 , n 1 + 1] × [n 2 -1 2 , n 2 + 1 2 ],
weighted by 1/2 if they are on an edge of the square, and by 1/4 if they are at one of its corners. Similarly,

(L * v + L * ↔ v ↔ + L * • v • ) 2 [n 1 , n 2 ]
is the sum of the horizontal part of the elements of v falling into the square

[n 1 -1 2 , n 1 + 1 2 ] × [n 2 , n 2 + 1]. Equating these two values to (Dx) 1 [n 1 , n 2 ] and (Dx) 2 [n 1 , n 2 ]
, respectively, is nothing but a discrete and 2-D version of the fundamental theorem of calculus, according to which the integral of a function on an interval is equal to the difference of its antiderivative at the interval bounds. So, we have defined the differentiation process from an image x to its gradient field v as the linear inverse problem of integration: integrating the gradient field v allows to recover the image x. Among all vector fields consistent with x in this sense, the gradient field v is selected as the simplest one, i.e. the one of minimal 1,2 norm.

Let us be more precise about this integration property connecting v to x. We first note that it is incorrect to interpret the pixel value x[n 1 , n 2 ] as a point sample of an unknown function s(t 1 , t 2 ), i.e. x[n 1 , n 2 ] = s(n 1 , n 2 ), and the values v

,1 [n 1 , n 2 ], v ↔,1 [n 1 , n 2 ], v •,1 [n 1 , n 2 ]
as point samples of ∂s/∂t 1 at (n 1 + 1 2 , n 2 ), (n 1 , n 2 + 1 2 ), (n 1 , n 2 ), respectively. Indeed, if it were the case, and viewing [START_REF] Alleysson | Linear demosaicing inspired by the human visual system[END_REF] as a kind of extended trapezoidal rule for numerical integration, the right-hand side of (24) would be divided by three. Instead, one can view x as the cell-average discretization of an unknown function s(t 1 , t 2 ), i.e. x[n 1 , n 2 ] =

n 1 +1/2 n 1 -1/2 n 2 +1/2 n 2 -1/2 s(t 1 , t 2 )dt 1 dt 2 ,
and v as the gradient field of s, in a distributional sense. For this, let us define the 1-D box and hat functions ( 26)

Π(t) =    1 if t ∈ (-1 2 , 1 2 ), 1 2 if t = ± 1 2 , 0 else , Λ(t) = Π(t) * Π(t) = max(1 -|t|, 0),
where * denotes the convolution. We also define the 2-D box function Π(t 1 , t 2 ) = Π(t 1 )Π(t 2 ) and the function ψ(t 1 , t 2 ) = Λ(t 1 )Π(t 2 ). The function or distribution ∂s/∂t 1 is such that

(Dx) 1 [n 1 , n 2 ] = x[n 1 + 1, n 2 ] -x[n 1 , n 2 ] = (s * Π)(n 1 + 1, n 2 ) -(s * Π)(n 1 , n 2 ) (27) = n 1 +1 n 1 ∂s ∂t 1 * Π (t 1 , n 2 )dt 1 (28) = ∂s ∂t 1 * ψ (n 1 + 1 2 , n 2 ). (29)
Then, the same equality holds, when replacing ∂s/∂t 1 by the distribution

v 1 (t 1 , t 2 ) = (n 1 ,n 2 )∈Ω v ,1 [n 1 , n 2 ]δ(t 1 -n 1 -1 2 , t 2 -n 2 ) + v ↔,1 δ(t 1 -n 1 , t 2 -n 2 -1 2 ) + v •,1 [n 1 , n 2 ]δ(t 1 -n 1 , t 2 -n 2 ), (30)
where δ(t 1 , t 2 ) is the 2-D Dirac distribution. Indeed,

( v 1 * ψ)(n 1 + 1 2 , n 2 ) = v ,1 [n 1 , n 2 ] + (v ↔,1 [n 1 , n 2 ] + v ↔,1 [n 1 , n 2 -1] + v ↔,1 [n 1 + 1, n 2 ] + v ↔,1 [n 1 + 1, n 2 -1])/4 + (31) (v •,1 [n 1 , n 2 ] + v •,1 [n 1 + 1, n 2 ])/2,
which, according to [START_REF] Alleysson | Linear demosaicing inspired by the human visual system[END_REF], is equal to

(L * v + L * ↔ v ↔ + L * • v • ) 1 [n 1 , n 2 ], which in turn is equal to (Dx) 1 [n 1 , n 2 ]
, by definition of v in [START_REF] Condat | A generic proximal algorithm for convex optimization-Application to total variation minimization[END_REF]. Altogether, the scalar field v 1 , the vertical component of the gradient field v, identified to the distribution v 1 , plays the same role as the partial derivative ∂s/∂t 1 of s, in the sense that they both yield the pixel values of x by integration. The same relationship holds between v 2 and ∂s/∂t 2 . To summarize, v is the discrete counterpart of the gradient of the unknown continuously-defined scene s, whose cell-average discretization yields the image x. So, it is legitimate to call v the gradient field of x. Note that there exists no function s such that ∇s is the Dirac brush ( v 1 , v 2 ), so v is no more than a discrete equivalent of ∇s.

We can notice that, given the image x, the gradient field v solution to [START_REF] Condat | A generic proximal algorithm for convex optimization-Application to total variation minimization[END_REF] is not always unique. For instance, for the 1-D signal x = (0, 0, 1/2, 1, 1), viewed as an image with only one row, one can set v = 0, v ↔ = 0, v • = (0, 0, 1, 0, 0). Another possibility is to take v = 0, v • = 0, v ↔ = (0, 1/2, 1/2, 0, 0). This possible nonuniqueness of v, which is very rare in practice, does not have any impact on the images obtained by TV minimization. We leave the study of theoretical aspects of the proposed TV and gradient field for future work, like showing Gamma-convergence of the proposed TV.

We end this section with a remark about the fact that the grid for the gradient field is twice finer than the one of the image. This factor of two appears naturally, according to the following sampling-theoretic argument. Let us consider a 2-D sine function s(t 1 , t 2 ) = sin(at 1 + bt 2 + c), for some a, b, c in (-π, π), which is sampled to give the image x, with

x[n 1 , n 2 ] = s(n 1 , n 2 ). We have |∇s(t 1 , t 2 )| 2 = (a 2 + b 2 ) cos 2 (at 1 + bt 2 + c) = (a 2 + b 2 ) cos(2at 1 + 2bt 2 + 2c)/2 + (a 2 + b 2 )/2.
So, by taking the squared amplitude of the gradient, the frequency of the sine is doubled. According to Shannon's theorem, the function |∇s| 2 must be sampled on a grid twice finer than the one of x, for its information content to be kept. Since, by virtue of the Fourier transform, every sufficiently regular function can be decomposed in terms of sines, this argument applies to an arbitrary 2-D function s, not only to a sine. The picture does not change by applying the square root, passing from |∇s| 2 to |∇s|, the integral of which is the TV of s. Thus, as long as the amplitude of the gradient is the information of interest, it must be represented on a twice finer grid; else aliasing occurs and the value of the TV becomes unreliable.

5. Algorithms for TV minimization. In this section, we focus on the generic convex optimization problem: [START_REF] Condat | A primal-dual splitting method for convex optimization involving Lipschitzian, proximable and linear composite terms[END_REF] Find x ∈ arg min

x∈R N 1 ×N 2 {F (x) + λ TV(x)} ,
where the sought-after image x has size N 1 × N 2 , λ > 0 is the regularization parameter, F is a convex, proper, lower semicontinuous function [START_REF] Bauschke | Convex Analysis and Monotone Operator Theory in Hilbert Spaces[END_REF]. A particular instance of this problem is image denoising or smoothing: given the image y, one solves:

(33) Find x ∈ arg min x∈R N 1 ×N 2 1 2 x -y 2 + λ TV(x) ,
where the norm is the Euclidean norm. This problem is a particular case of (32) with F (x) = 1 2 x -y 2 . More generally, many inverse problems in imaging can be written as: given the data y and the linear operator A, (34) Find x ∈ arg min

x∈R N 1 ×N 2 1 2 Ax -y 2 + λ TV(x) .
Again, this problem is a particular case of (32) with F (x) = 1 2 Ax -y 2 . Another instance is TV minimization subject to a linear constraint, for instance to regularize an ill-posed inverse problem in absence of noise: given the data y and the linear operator A, one solves: [START_REF] Deng | On the global and linear convergence of the generalized alternating direction method of multipliers[END_REF] Find x ∈ arg min

x∈R N 1 ×N 2 {TV : Ax = y} .
This problem is a particular case of (32) with λ = 1 and F (x) = ı {x : Ax=y} (x), where the convex indicator function ı Γ of a set Γ maps its variable x to 0 if x ∈ Γ, to +∞ else. When the TV is the anisotropic, isotropic, or upwind TV, which is a simple function composed with the finite differentiation operator D, there are efficient primal-dual algorithms to solve a large class of problems of the form [START_REF] Condat | A primal-dual splitting method for convex optimization involving Lipschitzian, proximable and linear composite terms[END_REF], see e.g. [START_REF] Chambolle | A first-order primal-dual algorithm for convex problems with applications to imaging[END_REF][START_REF] Condat | A generic proximal algorithm for convex optimization-Application to total variation minimization[END_REF][START_REF] Combettes | A forward-backward view of some primal-dual optimization methods in image recovery[END_REF] and references therein. In section 6, we use the overrelaxed version [START_REF] Condat | A primal-dual splitting method for convex optimization involving Lipschitzian, proximable and linear composite terms[END_REF] of the Chambolle-Pock algorithm [START_REF] Chambolle | A first-order primal-dual algorithm for convex problems with applications to imaging[END_REF]. With the proposed TV, it is not straightforward to apply these algorithms. In fact, (32) can be rewritten as: [START_REF] Combettes | Image restoration subject to a total variation constraint[END_REF] Find (x, v) ∈ arg min

x∈R N 1 ×N 2 ,v∈((R 2 ) N 1 ×N 2 ) 3 {F (x) + λ v 1,1,2 : L * v = Dx} .
So, one has to find not only the image x, but also its gradient field v, minimizing a separable function, under a linear coupling constraint. Let us introduce the function G(v) = λ v 1,1,2 and the linear operator C = -L * , so that we can put [START_REF] Combettes | Image restoration subject to a total variation constraint[END_REF] under the standard form:

(37) Find (x, v) ∈ arg min x∈R N 1 ×N 2 ,v∈((R 2 ) N 1 ×N 2 ) 3 {F (x) + G(v) : Cv + Dx = 0} .
The dual problem is [START_REF] Ring | Structural properties of solutions to total variation regularization problems[END_REF] Find û ∈ arg min

x∈u∈(R 2 ) N 1 ×N 2 {F * (-D * u) + G * (-C * u)} ,
which, in our case, is

(39) Find û ∈ arg min x∈u∈(R 2 ) N 1 ×N 2 {F * (-D * u) : Lu ∞,∞,2 ≤ λ} .
We now assume that the function F is simple, in the sense that it is easy to apply the proximity operator [START_REF] Bauschke | Convex Analysis and Monotone Operator Theory in Hilbert Spaces[END_REF][START_REF] Combettes | Proximal splitting methods in signal processing[END_REF] prox αF of αF , for any parameter α > 0. For the denoising problem [START_REF] Combettes | Proximal splitting methods in signal processing[END_REF], prox αF (x) = (x + αy)/(1 + α). For the regularized least-squares problem [START_REF] Ma | Alternating proximal gradient method for convex minimization[END_REF], prox αF (x) = (Id + αA * A) -1 (x + αA * y). For the constrained problem [START_REF] Deng | On the global and linear convergence of the generalized alternating direction method of multipliers[END_REF], prox αF (x) = x + A † (y -Ax), where A † is the Moore-Penrose pseudo-inverse of A. We also need the proximity operator of αG = αλ • 1,1,2 , which is (40)

prox αG (v) c [n 1 , n 2 ] = v c [n 1 , n 2 ] - v c [n 1 , n 2 ] max(|v c [n 1 , n 2 ]|/(αλ), 1) , ∀(n 1 , n 2 ) ∈ Ω, ∀c ∈ { , ↔, •}.
We can notice that D 2 ≤ 8 [2] and C 2 = L 2 ≤ 3. So, we have all the ingredients to use the Alternating Proximal Gradient Method [START_REF] Ma | Alternating proximal gradient method for convex minimization[END_REF], a particular case of the Generalized Alternating Direction Method of Multipliers [START_REF] Deng | On the global and linear convergence of the generalized alternating direction method of multipliers[END_REF]:

Algorithm 1 to solve [START_REF] Combettes | Image restoration subject to a total variation constraint[END_REF] Choose the parameters 0 < τ < 1/ D 2 , 0 < γ < 1/ C 2 , µ > 0, and the initial estimates x (0) , v (0) , u (0) . Then iterate, for i = 0, 1, . . .

      x (i+1) := prox τ µF x (i) -τ D * (Dx (i) + Cv (i) + µu (i) ) , v (i+1) := prox γµG v (i) -γC * (Dx (i+1) + Cv (i) + µu (i) ) , u (i+1) := u (i) + (Dx (i+1) + Cv (i+1) )/µ.
Assuming that there exists a solution to [START_REF] Combettes | Image restoration subject to a total variation constraint[END_REF], for which a sufficient condition is that there exists a minimizer of F , Algorithm 1 is proved to converge [START_REF] Ma | Alternating proximal gradient method for convex minimization[END_REF][START_REF] Deng | On the global and linear convergence of the generalized alternating direction method of multipliers[END_REF]: the variables x (i) , v (i) , u (i) converge respectively to some x, v, û, solution to [START_REF] Combettes | Image restoration subject to a total variation constraint[END_REF] and [START_REF] Alter | Evolution of characteristic functions of convex sets in the plane by the minimizing total variation flow[END_REF].

It is easy to show that the same algorithm can be used to compute the gradient field v of an image x, solution to (17); we simply replace x (i) by x. This yields Algorithm 2 to find v solution to [START_REF] Condat | A generic proximal algorithm for convex optimization-Application to total variation minimization[END_REF], given x Choose the parameters 0 < γ < 1/ C 2 , µ > 0, and the initial estimates v (0) , u (0) . Then iterate, for i = 0, 1, . . .

v (i+1) := prox γµG v (i) -γC * (Dx + Cv (i) + µu (i) ) , u (i+1) := u (i) + (Dx + Cv (i+1) )/µ.
In practice, we recommend setting τ = 0.99/8 and γ = 0.99/3 in Algorithm 1 and Algorithm 2, so that there only remains to tune the parameter µ.

Further on, let us consider the regularized least-squares problem [START_REF] Ma | Alternating proximal gradient method for convex minimization[END_REF], in the case where the proximity operator of the quadratic term cannot be computed. It is possible to modify Algorithm 1, by changing the metric in the Generalized Alternating Direction Method of Multipliers [START_REF] Deng | On the global and linear convergence of the generalized alternating direction method of multipliers[END_REF], to obtain a fully split algorithm, which only applies A and A * at every iteration, without having to solve any linear system. So, we consider the more general problem [START_REF] Condat | Fast projection onto the simplex and the l1 ball[END_REF] Find x ∈ arg min

x∈R N 1 ×N 2 F (x) + 1 2 Ax -y 2 + λ TV p (x) ,
or, equivalently,

(42) Find (x, v) ∈ arg min x∈R N 1 ×N 2 ,v∈((R 2 ) N 1 ×N 2 ) 3 F (x) + 1 2 Ax -y 2 + G(v) : Cv + Dx = 0 ,
where, again, G(v) = λ v 1,1,2 and C = -L * . The algorithm, with proved convergence to exact solutions of ( 41) and its dual, is:

Algorithm 3 to solve (42) Choose the parameters τ > 0, µ > 0, such that τ < 1/( D 2 + µ A 2 ), 0 < γ < 1/ C 2 ,
and the initial estimates x (0) , v (0) , u (0) . Then iterate, for i = 0, 1, . . .

      x (i+1) := prox τ µF x (i) -τ D * (Dx (i) + Cv (i) + µu (i) ) -τ µA * (Ax (i) -y) , v (i+1) := prox γµG v (i) -γC * (Dx (i+1) + Cv (i) + µu (i) ) , u (i+1) := u (i) + (Dx (i+1) + Cv (i+1) )/µ.
The proposed TV, like the other forms, could be used as a constraint, instead of being used as a functional to minimize [START_REF] Combettes | Image restoration subject to a total variation constraint[END_REF][START_REF] Fadili | Total variation projection with first order schemes[END_REF].

Many other algorithms could be applied to solve problems involving the proposed TV. The most appropriate algorithm for a particular problem must be designed on a case-by-case basis. So, it is beyond the scope of this paper to do any comparison of algorithms in terms of convergence speed.

Experiments.

In this section, we evaluate the proposed TV on several test problems. First, we report in Figure 1 the value of the proposed TV for the patterns shown in Figure 1. For each image, the value was determined by computing the associated gradient field, using Algorithm 2; these gradient fields are depicted in Figure 3. According to the discussion in section 2 and section 4, the proposed TV, which is a seminorm, takes appropriate values in all cases. We observe that, for binary patterns, the gradient field, and thus the value of the TV, is the same as with the anisotropic TV; that is, it is given by Eqs. ( 20)- [START_REF] Chambolle | An upwind finite-difference method for total variationbased image smoothing[END_REF]. Thus, the staircased nature of oblique binary patterns is penalized.

In the remainder of this section, we study the behavior of the proposed TV in several applications, based on TV minimization. Matlab code implementing the corresponding optimization algorithms and generating the images in Figure 3 to Figure 12 is available on the author's webpage.

6.1. Smoothing of a binary edge. We consider the smoothing problem [START_REF] Combettes | Proximal splitting methods in signal processing[END_REF] with the proposed TV, where the initial image y (N 1 = N 2 = 256), is an oblique binary edge, obtained by point sampling a continuously defined straight edge with slope 5/16. The central part of y is depicted in Figure 7 (a). So, we solve

(I) (II) (III) (IV) (V) (VI) (VII) (VIII) (IX) (X)
(43) Find (x, v) ∈ arg min x∈R N 1 ×N 2 ,v∈((R 2 ) N 1 ×N 2 ) 3 1 2 x -y 2 + λ v 1,1,2 : L * v = Dx ,
using Algorithm 1 (µ = 0.05, 2000 iterations). The central part of the smoothed image x, as well as the corresponding gradient field v, are depicted in Figure 7 (b), for λ = 2; see the caption of Figure 3 for the representation of the gradient field by colored arrows. The result for stronger smoothing with λ = 20 is depicted in Figure 7 (c).

We observe that the edge undergoes a slight blur, which remains concentrated over one or two pixels vertically, even for a strong smoothing parameter λ. This is expected, since such a slightly blurred edge has a lower TV value than the binary edge in y. Importantly, the minimization of the proposed TV tends to make all the gradient vectors of the field v aligned with the same orientation, which is exactly perpendicular to the underlying edge with slope 5/16. This shows that not only the amplitude, but also the orientation of the gradient vectors obtained with the proposed approach, is meaningful. 6.2. Smoothing of a disk. We consider the smoothing problem [START_REF] Combettes | Proximal splitting methods in signal processing[END_REF], with λ = 6, where y is the image of a white disk, of radius 32, over a black background (N 1 = N 2 = 99), depicted in Figure 8 (a). To simulate cell-average discretization, a larger (16N 1 ) × (16N 2 ) binary image was constructed by point sampling a 16 times larger disk, and then y was obtained by averaging over the 16 × 16 blocks of this image. In the continuous domain, it is known [START_REF] Ring | Structural properties of solutions to total variation regularization problems[END_REF] that TV smoothing of a disk of radius R and amplitude one over a zero background, with zero/Dirichlet boundary conditions, gives the same disk, with lower amplitude 1 -2λ/R, assuming λ < R/2. Here, we consider a square domain of size N 1 × N 2 with symmetric/Neumann boundary conditions, so the background is expected to become lighter after smoothing, with amplitude 2πλR/(N 1 N 2 -πR 2 ). We can notice that the total intensity remains unchanged and equal to πR 2 after smoothing. Moreover, according to the coarea formula, the TV of the image of a disk is 2πR, the perimeter of the disk, multiplied by the difference of amplitude between the disk and the background. Thus, in the discrete domain, we expect the smoothed image x to be similar to y, after an affine transform on the pixel values, so that the pixel values in the interior of the disk and in the background are 1 -2λ/R = 0.625 and 2πλR/(N 1 N 2 -πR 2 ) ≈ 0.183, respectively; this reference image is depicted in Figure 8 (b).

The images x obtained by solving [START_REF] Combettes | Proximal splitting methods in signal processing[END_REF] with the anisotropic, isotropic, upwind, and proposed TV (using 2000 iterations of Algorithm 1 with µ = 0.1), are shown in Figure 8.

-With the anisotropic TV the perimeter of the disk is evaluated in the sense of the Manhattan distance, and not the Euclidean distance. So, the TV of the disk is over-estimated. Since blurring an edge does not decrease the TV, TV minimization lets the TV value decrease by shrinking the shape of the disk and attenuating the amplitude of the edge more than it should.

-With the isotropic TV, the bottom, right, and top-left parts of the edge are sharp, but the other parts are significantly blurred. Contrary to the three other forms, the isotropic TV does not yield a symmetric image; the image is only symmetric with respect to the diagonal at -45 o .

-The upwind TV performs relatively well.

-The proposed TV outperforms the three other forms. Except at the top, bottom, left, right ends, the edge is sharper than with the upwind TV. The edge has the same spread everywhere, independently of the local orientation, which is a clear sign of the superior isotropy of the proposed approach. Since the proposed TV does not blur a horizontal or vertical edge after smoothing, the fact that the top, bottom, left, right ends of the disk edge are blurred here shows the truly nonlocal nature of the proposed TV; this is due to the higher number of degrees of freedom optimized during TV minimization, with not only the image but also its three gradient subfields. The other forms of the TV have less flexibility, with the gradient fully determined by local finite differences on the image.

The gradient field v, solution to (43), is depicted in Figure 4. We can observe its quality, with all the arrows pointing towards the disk center, showing that the gradient orientation is perpendicular to the underlying circular edge everywhere. 6.3. Smoothing of a square. We consider the smoothing problem [START_REF] Combettes | Proximal splitting methods in signal processing[END_REF], with λ = 6, where y is the image of a white square, of size 64 × 64, over a black background (N 1 = N 2 = 100), depicted in Figure 9 (a). In the continuous domain, the solution of the smoothing problem, when the function y is equal to 1 inside the square [-1, 1] 2 and 0 outside, λ < 1/(1 + √ π/2), and with zero boundary conditions, contains a square of same size, but with rounded and blurred corners, and lower amplitude [START_REF] Alter | Evolution of characteristic functions of convex sets in the plane by the minimizing total variation flow[END_REF][START_REF] Allard | Total variation regularization for image denoising. I. Geometric theory[END_REF]. The following closed-form expression can be derived: where

(44) x(t 1 , t 2 ) =        0 if |t 1 | > 1 or |t 2 | > 1, 0 else, if r ≤ λ, 1 -λ(1 + √ π/2) else, if r ≥ 1/(1 + √ π/2), 1 -λ/r else,
r = 2 -|t 1 | -|t 2 | + 2(1 -|t 1 |)(1 -|t 2 |).
Since symmetric, instead of zero, boundary conditions are considered here, x(t 1 , t 2 ) is actually the maximum of this expression and a constant, which can be calculated. So, the reference result in the discrete case was simulated by point sampling this function x(t 1 , t 2 ) on a fine grid, with λ = 6/32, in a large 1600 × 1600 image, which was then reduced by averaging over its 16 × 16 blocks. This reference image is depicted in Figure 9 (b).

The image x, solution to [START_REF] Combettes | Proximal splitting methods in signal processing[END_REF] with the anisotropic, isotropic, upwind, and proposed TV (using 2000 iterations of Algorithm 1 with µ = 0.3), is shown in Figure 9. The anisotropic TV yields a square, without any rounding of the corners. This shows again that the metric underlying anisotropic TV minimization is not the Euclidean one. With the isotropic TV, the asymmetric blur of the corners contaminates the top and left sides of the square. Only the top-left corner has the correct aspect. With the upwind TV, the level lines at the corners are more straight than circular. The proposed TV yields the closest image to the reference image.

6.4. Denoising of the Bike. We consider the application of the smoothing/denoising problem [START_REF] Combettes | Proximal splitting methods in signal processing[END_REF], or (43) with the proposed TV, to remove noise in a natural image. The initial image y, depicted in Figure 10 (a), is a part of the classical Bike image, depicted in Figure 10 (b), corrupted by additive white Gaussian noise of standard deviation 0.18. λ is set to 0.16. With the anisotropic TV, the noise is removed, but the contrast of the spokes is more attenuated than with the other forms of the TV. With the isotropic TV, the noise is less attenuated and some small clusters of noise remain. This is also the case, to a much larger extent, with the upwind TV: the dark part of the noise is removed, but not the light part, and a lot of small light clusters of noise remain. This drawback of the isotropic and upwind TV can be explained by the too low penalization of a single isolated pixel, as reported in Table 1 and in section 2. The proposed TV (using 1000 iterations of Algorithm 1 with µ = 1) yields the best result: the noise is removed, the spokes have an elongated shape with less artifacts and a good contrast. With the anisotropic TV, the solution is not unique, and every image with nondecreasing pixels values horizontally and vertically is a solution of the TV minimization problem. One solution, equal to y, is shown in Figure 5 (b). The result with the isotropic, upwind, and proposed TV (using 1000 iterations of Algorithm 1, with µ = 1) is the same, and corresponds to what is expected; it is shown in Figure 5 (c). The gradient field v associated to the solution with the proposed TV is not shown, but it is the same as in Figure 3 (III).

We also consider the flipped case, where y is flipped horizontally. The solution with the isotropic TV is shown in Figure 5 (d). It suffers from a strong blur. Indeed, as reported in Table 1, the value of the isotropic TV for slightly blurred edges at this orientation, like in the cases (IIIf) and (IVf), is too high. So, when minimizing the TV, the TV value is decreased by the introduction of an important blur. By contrast, the anisotropic, upwind, and proposed TV are symmetric, so they yield flipped versions of the images shown in Figure 5 (b) and (c).

6.6. Upscaling of a disk. We consider the upscaling problem, which consists in increasing the resolution of the image y of a disk, shown in Figure 11 (a), by a factor of 4 in both directions. Upscaling is viewed as the inverse problem of downscaling: the downscaling operator A maps an image to the image of its averages over 4 × 4 blocks, and we suppose that y = Ax , for some reference image x , that we want to estimate. Here, y is of size 23 × 23 and the reference image x , shown in Figure 11 (b), of size 92 × 92, was constructed like in subsection 6.2: to approximate cell-average discretization, a larger 1472×1472 image x 0 was constructed by point sampling a 16 times larger disk, and x was obtained by averaging over the 16 × 16 blocks of this image; that is, x = AAx 0 . Then y was obtained as y = Ax . Hence, the upscaled image is defined as the solution to the constrained TV minimization problem [START_REF] Deng | On the global and linear convergence of the generalized alternating direction method of multipliers[END_REF]. We have

A † = 16A * .
The results with the anisotropic, isotropic, upwind, and proposed TV (using 2000 iterations of Algorithm 1, with µ = 1) are shown in Figure 11 (c)-(f). With the anisotropic TV, the result is very blocky. With the isotropic TV, the disk edge is jagged, except at the top-left and bottom-right ends. The result is much better with the upwind TV, and even better with the proposed TV, which has the most regular disk edge. The distance x -x between the upscaled image and the reference image is 2.91, 1.59, 1.23, with the isotropic, upwind, proposed TV, respectively. So, this error is 23% lower with the proposed TV than with the upwind TV. 6.7. Deconvolution of a disk. We consider the deconvolution, a.k.a. deblurring, problem, which consists in estimating an image x , given its blurred and noisy version y. x is the image of a disk, constructed like in subsection 6.2 and shown in Figure 6 (b). The initial image y, depicted in Figure 6 (a), was obtained by applying a Gaussian filter, of standard deviation (spread) 3.54 pixels, to x and adding white Gaussian noise of standard deviation 0.05. The image was restored by solving, for the proposed TV, (42) with F = 0 and, for the other TV forms, [START_REF] Ma | Alternating proximal gradient method for convex minimization[END_REF]. In all cases, A = A * is convolution with the Gaussian filter, with symmetric boundary conditions, and λ is set to 0.1. Algorithm 3 was used in all cases, for simplicity, with C = -Id for all but the proposed TV. The distance x -x between the restored image and the reference image is 4.10, 3.56, 2.78, 1.46, with the anisotropic, isotropic, upwind, proposed TV, respectively. We observe in Figure 6 (c)-(f) that the noise is well removed in all cases.

Again, the proposed TV provides the roundest and least blurred edge of the disk. 6.8. Segmentation of the Parrot. We consider a convex approach to color image segmentation. Given the set Σ = {c k ∈ [0, 1] 3 : k = 1, . . . , K} of K ≥ 2 colors c k , expressed as triplets of R,G,B values, and the color image y ∈ (R 3 ) N 1 ×N 2 , we would like to find the segmented image (45)

x = arg min

x∈Σ N 1 ×N 2 1 2 x -y 2 + λ 2 K k=1 per(Ω k ) ,
for some λ > 0, where

Ω k = {(n 1 , n 2 ) ∈ Ω : x[n 1 , n 2 ] = c k } and per denotes the perimeter.
That is, we want a color image, whose color at every pixel is one of the c k , close to y, but at the same time having homogeneous regions. However, this nonconvex "Potts" problem is very difficult, and even NP-hard [START_REF] Chambolle | A convex approach to minimal partitions[END_REF]. And a rigorous definition of the perimeter of a discrete region is a difficulty in itself. So, we consider a convex relaxation of this problem [START_REF] Chambolle | A convex approach to minimal partitions[END_REF]: we look for the object ẑ ∈ ∆ N 1 ×N 2 , such that, at every pixel, ẑ

[n 1 , n 2 ] = (ẑ k [n 1 , n 2 ]) K k=1 is an assignment vector in the simplex ∆ = {(a k ) K k=1 : K k=1 a k = 1 and a k ≥ 0, ∀k}. The elements ẑk [n 1 , n 2 ] ∈ [0, 1] are the proportions of the colors c k at pixel (n 1 , n 2 ); that is, the segmented image x is obtained from ẑ as (46) x[n 1 , n 2 ] = K k=1 ẑk [n 1 , n 2 ]c k , ∀(n 1 , n 2 ) ∈ Ω.
Now, by virtue of the coarea formula, the segmentation problem can be reformulated as [START_REF] Chambolle | A convex approach to minimal partitions[END_REF] (47) Find ẑ = arg min

z∈∆ N 1 ×N 2 z, p + λ K k=1 TV(z k ) ,
where the Euclidean inner product is

(48) z, p = (n 1 ,n 2 )∈Ω K k=1 z k [n 1 , n 2 ]p k [n 1 , n 2 ], (49)
with p k [n 1 , n 2 ] = y[n 1 , n 2 ] -c k 2 .
The problem (47) can be put under a form similar to (32):

(50) Find ẑ ∈ arg min

z∈(R K ) N 1 ×N 2 F(z) + λ TV(z) ,
with the TV of z having a separable form with respect to k, i.e. TV(z) = K k=1 TV(z k ), and F(z) having a separable form with respect to the pixels, i.e.

F(z) = (n 1 ,n 2 )∈Ω F n 1 ,n 2 (z[n 1 , n 2 ]), where (51) F n 1 ,n 2 (a) = ı ∆ (a) + a, p[n 1 , n 2 ] .
For any α > 0, we have prox αFn 1 ,n 2 (a) = P ∆ (a -αp[n 1 , n 2 ]), where P ∆ is the projection onto the simplex, which can be computed efficiently [START_REF] Condat | Fast projection onto the simplex and the l1 ball[END_REF]. So, the primal-dual algorithms described in section 5 can be used for the segmentation problem, as well. With the proposed TV, we must introduce K gradient fields v k , associated to the images z k . We used 1000 iterations of Algorithm 1, with µ = 50.

We compare the performances of the anisotropic, isotropic, upwind, proposed TV on this problem, with y a part, of size 399 × 400, of the classical Parrot image, shown in Figure 12 (a). We set λ = 0.09 and we set the K = 6 colors as some kind of black, white, yellow, blue, green, brown, visible in Figure 12 (b)-(e). In this respect, we would like the edges, which are the interfaces between the regions Ω k , to be sharp, and their perimeter to be correctly measured by the TV of the assignment images ẑk . But these two goals are antagonist: the coarea formula is not well satisfied for discrete binary shapes, as we have seen in section 2: the length of oblique binary edges is overestimated by the anisotropic, isotropic, and proposed TV, and the length of small structures, like in the extreme case of a single isolated pixel, is underestimated by the upwind TV. This seems like an intrinsic limitation and the price to pay for convexity, in a spatially discrete setting. As visible in Figure 12 (b), the anisotropic TV yields sharp edges, but their length is measured with the Manhattan distance, not the Euclidean one. So, the edges tend to be vertical and horizontal. With the isotropic TV, for half of the orientations, the edges are significantly blurred, as is visible on the dark region over a green background, in the bottom-left part of the image in Figure 12 (c). The upwind TV tends to introduce more regions made of a few pixels, because their perimeter is underestimated, see the eye of the parrot in Figure 12 (d). The best tradeoff is obtained with the proposed TV: there is a slight, one or two pixel wide blur at the edges, but this blur cannot be avoided, for the perimeter of the regions to be correctly evaluated.

Conclusion.

We proposed a new formulation for the discrete total variation (TV) seminorm of an image. Indeed, the classical, so-called isotropic, TV suffers from a poor behavior on oblique structures, for half of the possible orientations. It is important to have a sound definition of the TV, not least to be able to compare different convex regularizers for imaging problems, based on their intrinsic variational and geometrical properties, and not on the quality of their implementation.

Our new definition of the gradient field of an image has potential applications going far beyond TV minimization; for instance, one can consider edge detection based on the gradient amplitude, nonlinear diffusion and PDE flows based on the gradient orientation, one can define higher order differential quantities. . . We will explore some of these problematics in future work. The extension of the proposed TV to color or multichannel images will be investigated, as well.

Figure 2 .

 2 Figure 2. In (a), (b), (c), continuously-defined images whose cell-average discretization yields Figure 1 (III), (IV), (VII), respectively.

  (c), Figure 5 (d), Figure 12 (b).

Figure 3 .

 3 Figure 3. Same patterns as in Figure 1, with the associated gradient fields, solutions to (17). The vectors v [n1, n2], v↔[n1, n2], v•[n1, n2], are represented by red, blue, green arrows, starting at (n1 + 1 2 , n2), (n1, n2 + 1 2 ), (n1, n2), respectively.

Figure 4 .

 4 Figure 4. Zoom on the top-left part of the disk edge in Figure 8 (f ), with the associated gradient field.

Figure 5 .

 5 Figure 5. Inpainting experiment, see subsection 6.5. The region to reconstruct is in blue in (a). In (b), one solution of anisotropic TV minimization. In (c), solution of isotropic, upwind, and proposed TV minimization. In (d), solution of isotropic TV minimization, for the flipped case.

6. 5 .

 5 Inpainting of an edge. We consider an inpainting problem, which consists in reconstructing missing pixels by TV minimization. The image is shown in Figure5(a), with the missing pixels in blue. We solve the constrained TV minimization problem[START_REF] Deng | On the global and linear convergence of the generalized alternating direction method of multipliers[END_REF], where A is a masking operator, which sets to zero the pixel values in the missing region and keeps the other pixels values unchanged. We have A † = A * = A. The image y, shown in Figure5(b), has its pixel values in the missing region equal to zero.

  (a) Initial image (b) Reference image (c) Restored, anisotropic TV (d) Restored, isotropic TV (e) Restored, upwind TV (f) Restored, proposed TV

Figure 6 .

 6 Figure 6. Deconvolution experiment, see subsection 6.7. The initial image in (a) is obtained from the ground-truth image in (b), by convolution with a Gaussian filter and addition of white Gaussian noise.

  (a) Initial image (central part) (b) Smoothed image (central part), proposed TV, λ = 2 (c) Smoothed image (central part), proposed TV, λ = 20

Figure 7 .

 7 Figure 7. Smoothing experiment, see subsection 6.1. In (b) and (c), central part of the images and their gradient fields obtained by smoothing the binary edge in (a), with λ = 2 and λ = 20, respectively.

  Smoothed image, anisotropic TV

Figure 8 .

 8 Figure 8. Smoothing experiment, see subsection 6.2. In (c), the image obtained by smoothing the image in (a), using the anisotropic TV. In (b), the ideal result one would like to obtain. Every image is represented in grayscale on the left and in false colors on the right, to better show the spread of the edges. The fact that in (b) the disk interior and the background are rendered with the same blue false color is a coincidence, due to the limited number of colors in the colormap 'prism' of Matlab.

  (d) Smoothed image, isotropic TV (e) Smoothed image, upwind TV (f) Smoothed image, proposed TV

Figure 8 ,

 8 Figure 8, continued. In (d), (e), (f ), the images obtained by smoothing the image in (a), using the isotropic TV, upwind TV, proposed TV, respectively.

  Smoothed image, anisotropic TV

Figure 9 .

 9 Figure 9. Smoothing experiment, see subsection 6.3. In (c), the image obtained by smoothing the image in (a), using the anisotropic TV. In (b), the ideal result one would like to obtain. Every image is represented in grayscale on the left and in false colors on the right, to better show the spread of the corners.

  (d) Smoothed image, isotropic TV (e) Smoothed image, upwind TV (f) Smoothed image, proposed TV

Figure 9 ,

 9 Figure 9, continued. In (d), (e), (f ), the images obtained by smoothing the image in (a), using the isotropic TV, upwind TV, proposed TV, respectively. The fact that in (e) the square interior and the background are rendered with the same blue false color is a coincidence, due to the limited number of colors in the colormap 'prism' of Matlab.

  (c) Denoised image, anisotropic TV (d) Denoised image, isotropic TV

Figure 10 ,

 10 Figure 10, continued. In (c), (d), the images obtained by denoising the image in (a), using the anisotropic TV and isotropic TV, respectively.

  (e) Denoised image, upwind TV (f) Denoised image, proposed TV

Figure 10 ,

 10 Figure 10, continued. In (e), (f ), the images obtained by denoising the image in (a), using the upwind TV and proposed TV, respectively.

  Upscaled image, anisotropic TV (d) Upscaled image, isotropic TV (e) Upscaled image, upwind TV (f) Upscaled image, proposed TV

Figure 11 .

 11 Figure 11. Upscaling experiment, see subsection 6.6. The images in (b)-(f ), when reduced by averaging over 4 × 4 blocks, yield the image in (a), exactly.

  (a) Initial image (b) Segmented image, anisotropic TV (c) Segmented image, isotropic TV (d) Segmented image, upwind TV (e) Segmented image, proposed TV

Figure 12 .

 12 Figure 12. Segmentation experiment, see subsection 6.8.

Table 1

 1


Fairouz Zobiri
email: fairouz.zobiri@imag.fr

Nacim Meslem
email: nacim.meslem@gipsa-lab.fr

Brigitte Bidegaray
email: brigitte.bidegaray@imag.fr

Event-Based Sampling Algorithm for Setpoint Tracking Using a State-Feedback Controller

Event-based control techniques are investigated for output reference tracking in the case of linear time-invariant systems. In event-based control, the controller remains at rest if the system is behaving according to some predefined conditions, the feedback loop being closed only when the system states violate these conditions. In this work a reference system, which consists in the continuously-controlled version of the system under study, is employed. Based on the difference between the state of the event-triggered system and that of the reference system, we define a Lyapunov-like function, and show that if we can keep this function confined to a certain region, the tracking error would also be bounded. The trespassing of this function outside of the desired region is used as an event-triggering condition.

I. INTRODUCTION

Event-based control is a control strategy in which the control law is updated following the occurence of an event related to the state of the system under study. It has appeared in the last decades as an appealing alternative to periodic control. The ever-growing number of publications on the topic ( [START_REF] Årzén | A simple event-based PID controller[END_REF], [START_REF] Lehmann | Event-based control: A state-feedback approach[END_REF], [START_REF] Lunze | A state-feedback approach to event-based control[END_REF], [START_REF] Heemels | Periodic eventtriggered control for linear systems[END_REF], [START_REF] Meslem | Event-based controller synthesis by bounding methods[END_REF] and references therein) testifies of this shift of interest away from periodic control.

By periodic control, we refer to the classical approach that consists in uniformly sampling a control law in order to transfer it onto a digital platform. The discretization is thus carried out by taking a sample value of the continuous-time signal at fixed time intervals. Of course, the popularity of this approach is due, apart from its orderly nature, to the fact that it has been exhaustively investigated for several decades.

However, lately, the control community has been pondering over the shortcomings of the periodic implementation in modern control facilities. First, periodic sampling requires the use of a very small time-step in order to satisfy the Nyquist-Shannon sampling theorem and keep the properties of the original, continuous-time signal. But, this is often met with hardware limitations, as the sampling period can only get as small as the hardware resolution allows it to. Then, in a modern industrial environment, applications are often set to communicate through a bandlimited (sometimes wireless) network, where transferring data at a very high rate may lead to the congestion of the network [START_REF] Zhang | Survey on recent advances in networked control systems[END_REF], possibly resulting in disastrous consequences on the behavior of the interconnected systems. It also represents a considerable waste of energy, along with all its environmental and economical implications.

From such conclusions stemmed the idea that control should only be applied when needed. In event-triggered control, the value of the control signal is kept constant as long as the system does not undergo any significant change. A new sample of the control signal is taken when the state of the system reaches unacceptable values or when some stability condition is violated. During the time the control is held unchanged, the controller is at rest, and no state measurement is fed back to it, reducing thus energy consumption and reducing the load on the communication channels. It has also the merit of being a more intuitive approach, as humans tend to turn their attention on tasks that need intervention rather than on those which are running smoothly. However, it does present the disadvantage of having to monitor the state of the system continually.

Up to now, the research on event-based strategies has been focused on stability and stabilization issues ( [START_REF] Anta | To sample or not to sample: Self-triggered control for nonlinear systems[END_REF], [START_REF] Girard | Dynamic triggering mechanisms for event-triggered control[END_REF], [START_REF] Postoyan | A unifying Lyapunovbased framework for the event triggered control of nonlinear systems[END_REF], to name a few). Inspired by results developed in the literature, we introduce an event-based approach for output tracking of a reference trajectory by a linear time-invariant (LTI) system.

A few works on event-based control for reference tracking do exist throughout the literature. For instance, in [START_REF] Durand | Event-based LQR with integral action[END_REF], an event-based LQR controller is implemented, while the reference tracking is achieved through integral action. In [START_REF] Tallapragada | On event triggered trajectory tracking for control affine nonlinear systems[END_REF], [START_REF]On event triggered tracking for nonlinear systems[END_REF], [START_REF] Shao | Event-triggered and self-triggered asymptotically tracking control of linear systems[END_REF], an additional fictitious system that generates the desired trajectory is defined. The event-triggering conditions are then established according to the error between the state of the reference system and that of the actual system.

In all these approaches, the reference system is chosen to be an external system, that would suit the purpose of analysis. In this work though, we propose to select the reference system as the continuously-controlled version of the event-triggered system being considered. It is indeed reasonable to assume that the behavior of the continuouslycontrolled system should serve as a role model for the performance of the event-based algorithm. Our main objective then would be to maintain the tracking error between the states of the two systems under a certain bound. This is done by noticing that this error can be viewed as the state of a new virtual system, and will ultimately decay to zero if the system is stable. The practical stability of the error-related system is then guaranteed by associating a Lyapunov-like function to it, which would be kept from increasing above a certain level, defining thus our triggering conditions. This work is organized as follows: in the next section we expose the mathematical framework of the problem at hand by explaining the event-based scheme. In section III, we define the event-triggering conditions that ensure the boundedness of the error. Some simulation results are then presented in section VI.

II. PROBLEM STATEMENT

Let us consider the LTI system described as follows ẋ(t ) = Ax(t ) + Bu(t ),

y(t ) = C x(t ), (1)
x(t ) ∈ R n is the state vector, y(t ) ∈ R p is the system output and u(t ) ∈ R m is the control input. We wish to drive the output of system (1) to asymptotically follow a trajectory determined by an exogenous signal r (t ), denoted as the reference input. We first define the reference system which consists in system (1) for which the control law is classically defined as

u r (t ) = -K x(t ) +Gr (t ), (2)
where K ∈ R m×n is the feedback gain and G ∈ R m×p is a calibration matrix, where

G = (-C (A -B K ) -1 B ) -1 .
As a result, we obtain the following closed-loop system

ẋr (t ) = (A -B K )x r (t ) + BGr (t ), y r (t ) = C x r (t ), (3)
such that 1) (A -B K ) is Hurwitz with desired eigenvalues, 2) lim t →∞ y r (t )r (t ) = 0.

The control law u r (t ) is applied continuously to system (1) and classical control theory affirms that achieving 1) and 2) is always possible when the pair (A, B ) is controllable. We assume in what follows that it is the case for system (1).

In event-based control, however, the control is updated only at some time instants t k , k ∈ N, when certain conditions on the behavior of the plant are no longer satisfied, resulting in time instants t k which are not necessarily evenly spaced. The control input, now denoted by ū(t ), is defined

at time t k by ū(t k ) = -K x(t k ) +Gr (t k ). (4)
It is maintained constant on time intervals [t k , t k+1 ):

ū(t ) = ū(t k ) ∀t ∈ [t k , t k+1 ). ( 5
)
Then system (1) becomes

ẋ(t ) = Ax(t ) + B ū(t k ), y(t ) = C x(t ), t ∈ [t k , t k+1 ), k ∈ N. (6)
As mentioned previously, we consider the continuouslycontrolled system (3) as a reference system, meaning that, at each time instant, the behavior of the event-based system will be compared to that of the reference system. From this comparison, we will determine whether the behavior of system ( 6) is acceptable or whether the control should be updated.

Remark 1. The principle of a reference system in the eventbased scheme has also been used in references [START_REF] Tallapragada | On event triggered trajectory tracking for control affine nonlinear systems[END_REF], [START_REF]On event triggered tracking for nonlinear systems[END_REF], [START_REF] Shao | Event-triggered and self-triggered asymptotically tracking control of linear systems[END_REF]. Our approach differs in the fact that the reference system is provided within the method, whereas in [START_REF] Tallapragada | On event triggered trajectory tracking for control affine nonlinear systems[END_REF] and [START_REF]On event triggered tracking for nonlinear systems[END_REF] the choice of the reference system is left for the user to make. The work presented in [START_REF] Shao | Event-triggered and self-triggered asymptotically tracking control of linear systems[END_REF] deals with a stability problem since the states are driven to zero.

Since we have assumed that all the conditions are united to achieve lim t →∞ y r (t )r (t ) = 0, it is sufficient to guarantee that lim t →∞ y(t )y r (t ) = 0, to ensure that the output y(t ) of the event-triggered system follows the trajectory of r (t ).

In this work, however, we consider the more general problem of driving the state error x(t )x r (t ) to zero as time t tends to ∞. This principle is depicted in Fig. 1. From the relationship between x(t ) and y(t ) shown in system (1), we deduce that lim t →∞ x(t

) -x r (t ) = 0 implies lim t →∞ y(t ) -y r (t ) = 0.
Notice that for the plants where the states are not available for measurement, techniques for obtaining an estimation of these states (such as observers) do exist and can be used alongside our approach with only a small modification.

Defining the event-triggering conditions is the topic of the next section.

III. EVENT-TRIGGERING CONDITIONS FOR TRACKING A REFERENCE INPUT

In this section, we define a new system based on the reference and the event-based systems. The state of this new system, the error between the real state x(t ) and the reference state x r (t ),

e(t ) = x(t ) -x r (t ), (7)
will serve as a basis for updating the control ū(t ). The behavior of the error signal e can be described as follows

1) when t ∈ (t k , t k+1 ) ė(t ) = (A -B K )e(t ) -B K ∆ k x(t ) + BG∆ k r (t ), (8a)
where

∆ k x(t ) = x(t k ) -x(t ), ∆ k r (t ) = r (t k ) -r (t ), 2) when t = t k , k = 0, 1, 2, . . . ė(t ) = (A -B K )e(t ). ( 8b
)
Remark 2. Throughout this work, we allow the exogenous signal r (t ) to be only piecewise Lipschitz, and we authorize it to contain jumps, provided that there exists a minimum interval of time between two successive jumps where r (t ) is continuous. Moreover, we require the signals used as exogenous inputs r (t ) to be uniformly bounded.

Remark 3. Note also that despite the jumps in r (t ), the states x r (t ), x(t ) and e(t ) remain Lipschitz. However, their first derivatives do exhibit jumps. Indeed, jumps in ẋ(t ) occur when an event happens and jumps in ẋr (t ) are due to jumps in r (t ). Therefore ė(t ) exhibits jumps both when either ẋ(t ) or ẋr (t ) exhibit jumps. Based on equation (8b), we can associate to the system (8) a Lyapunov-like function of the following form V (e(t )) = e(t ) T P e(t ), [START_REF] Postoyan | A unifying Lyapunovbased framework for the event triggered control of nonlinear systems[END_REF] where P , a symmetric positive-definite matrix, is a solution to the Lyapunov equation

r (t ) Mux ẋr = (A -B K )x r + BGr y r = C x r u(t ) = 0 ū(t k ) = -K x +Gr ẋ = Ax + B ū y = C x y x r Event Generator V (e(t )) = δ x x(t )
(A -B K ) T P + P (A -B K ) = -Q, ( 10
)
where Q is a positive-definite matrix. Equations ( 9) and ( 10) stem from the fact that V (e) should be a positive-definite function, whose first timederivative is non-positive. These conditions, according to Lyapunov's theory, are sufficient to guarantee the stability of system [START_REF] Girard | Dynamic triggering mechanisms for event-triggered control[END_REF].

In our approach though, we do not try to prevent V (e) from increasing, but we show that we are capable of maintaining V (e) confined within a certain region. In other words, between two time instants t k and t k+1 , V (e) may increase; but, whenever V (e) reaches a certain upper bound, the control is updated, forcing V (e) back inside the chosen region. This concept is illustrated in Fig. 2.

Let δ > 0 be the upper limit that we want to achieve on V (e), 0 ≤ V (e) ≤ δ. [START_REF] Tallapragada | On event triggered trajectory tracking for control affine nonlinear systems[END_REF] From equation ( 9) we can deduce that the Lyapunov-like function V (e) is naturally bounded as follows

λ min P e(t ) 2 ≤ V (e) ≤ λ max P e(t ) 2 , ( 12
)
where λ min P and λ max P are, respectively, the smallest and largest eigenvalues of matrix P . Since P is positive-definite, its eigenvalues are real and positive.

By requiring that V (e) be bounded from above, and from inequality [START_REF]On event triggered tracking for nonlinear systems[END_REF], we can guarantee an upper bound on the error e(t ), expressed as

e(t ) ≤ δ λ min P = . ( 13
)
Consequently we can define our triggering conditions based on this premise. Definition 1. We define the time-instant t k+1 (k ∈ N) at which the control-law ū(t ) is updated as the minimum time instant t > t k for which V (e(t )) = δ:

t k+1 = inf{t > t k ,V (e(t )) = δ}. ( 14
)
Theorem 1. If the event-based control ū(t ) is updated according to the event-triggering condition defined in [START_REF]Yaw damper design for a 747® jet aircraft[END_REF], then the tracking error e(t ) remains confined in the ball of radius , i.e. e(t ) ≤ ,

for defined by [START_REF] Shao | Event-triggered and self-triggered asymptotically tracking control of linear systems[END_REF].

Proof: To prove Theorem 1, we study the time derivative of the Lyapunov-like function V defined by equation [START_REF] Postoyan | A unifying Lyapunovbased framework for the event triggered control of nonlinear systems[END_REF]. Since the matrix P is positive-definite, for all e ∈ R n , V (e) ≥ 0. At time t 0 , e(t ) is chosen to be zero (we assume x(t 0 ) = x r (t 0 )). Then, at time t = t k , when V (e(t )) reaches the value δ the control is updated such that V (e(t )

) |t =t k = ė(t ) T Pe(t ) + e(t ) T P ė(t ) |t =t k = e(t k ) T (A -B K ) T Pe(t k ) +e(t k ) T P (A -B K )e(t k ) = -e(t k ) T Qe(t k ).
From equation [START_REF] Durand | Event-based LQR with integral action[END_REF], the matrix Q is symmetric positivedefinite, meaning that V (e) is negative-definite, V (e(t )) decreases at time t = t k and is pushed back to the region comprised between the lines V (e(t )) = 0 and V (e(t )) = δ. Since Q is a positive-definite matrix, its largest and smallest eigenvalues, λ max Q and λ min Q , respectively, are real and positive. Then, we can set an upper and lower limit on the descent of V (e(t )),

-λ max Q e(t k ) 2 ≤ -e(t k ) T Qe(t k ) ≤ -λ min Q e(t k ) 2 .
(16)

We have demonstrated earlier that in this region e(t ) ≤ . Thus V (e(t )

) |t =t k ≤ -λ min Q 2 for all k.
However, the boundedness of the tracking error e(t ) is not sufficient to deem the triggering condition in Definition 1 appropriate. We need also to guarantee that under these triggering conditions, there exists a minimum time lapse between two update instants. This extra requirement avoids, to borrow the terminology of hybrid systems, Zeno phenomena, where an infinite number of updates is needed within a finite interval of time. This also allows to give an estimation of the number of necessary updates of the control, for the purpose of minimizing the communication cost of the event-triggered control. We can prove the existence of this minimum delay in the case when r (t ) Lipschitz.

Theorem 2. Let r (t ) be a Lipschitz input signal. Then, there exists a minimum time τ min > 0, independent of k, such that ∀k ∈ N, t k+1t k > τ min , where the t k , k ∈ N are defined in Definition 1.

Remark 4.

Assuming that r is Lipschitz is not a very strong assumption, since in practice, a low-pass filter is used on the reference input to handle the abrupt changes that it may contain and avoid actuator saturation, eventually rendering it Lipschitz.

Proof: We show that there exists τ independent of k such that V (e(t )) remains negative in the interval t ∈ [t k , t k + τ). We will necessarily end up having τ < τ min . For t > t k , V (e(t )) = -e(t ) T Qe(t )

-2∆ k x(t ) T K T B T Pe(t ) + 2∆ k r (t ) T G T B T Pe(t ) ≡ -e(t ) T Qe(t ) + R k (t ).
We prove that there exist β > 0 and τ 1 > 0, independent of k, such that for all t ∈ (t k , t k + τ 1 ), -e(t ) T Qe(t ) < -β, and 0 < τ 2 < τ 1 such that for all t ∈ (t k ,

t k + τ 2 ), |R k (t )| < β/2. Hence for t ∈ (t k , t k + τ 2 ), V (e(t )
) is uniformly bounded from above by -β/2. To determine τ 1 , we have to prove that e(t ) 2 remains relatively large on some time interval.

From equations ( 11), [START_REF]On event triggered tracking for nonlinear systems[END_REF], and ( 16)

-e(t k ) T Qe(t k ) ≤ - λ min Q δ λ max P .
We therefore consider β = λ min Q δ/2λ max P . To have -e(t ) T Qe(t ) < -β it suffices to show that e(t ) > β/λ min Q = δ/2λ max P . We can rewrite ė(t ) as

ė(t ) = (A -B K )e(t ) + F k (t ), ( 17
)
where

F k (t ) = -B K ∆ k x(t ) + BG∆ k r (t ).
In order to find a bound on F k (t ), we need the following lemma.

Lemma 1. If r (t ) is an L r -Lipschitz function, then the following assumptions are verified.

(i) There exists L r > 0 such that

∆ k r (t ) ≤ L r (t -t k ). ( 18
) (ii)
As noticed in Remark 3, x is Lipschitz. Let L x be its Lipschitz constant. Therefore,

∆ k x(t ) ≤ L x (t -t k ). ( 19
)
From Lemma 1, we can state that

F k (t ) ≤ B K L x (t -t k ) + BG L r (t -t k ) ≤ C τ 1 ,
• being the subordinated matrix-norm in R n×n and R n×p , and C constant. For t ∈ (t k , t k + τ 1 ), we can also write an integral formula for (17), namely

e(t ) = e (A-B K )(t -t k ) e(t k ) + t t k e (A-B K )(t -s) F k (s)d s. Since A -B K is Hurwitz, e(t ) ≥ e(t k ) -C τ 1 2 .
From [START_REF]On event triggered tracking for nonlinear systems[END_REF] we also know that e(t k ) ≥ δ/λ max P . Therefore for a sufficiently small τ 1 , which does not depend on k, we can ensure that e(t ) ≥ δ/2λ max P and hence -e(t ) T Qe(t ) < -β. We then move on to prove that |R k (t )| is bounded from the above by β/2.

|R k (t )| ≤ 2 P B K e(t ) ∆ k x(t ) +2 P BG e(t ) ∆ k r (t ) ≤ 2 δ λ min P τ 2 ( |P B K |L x + P BG L r ) ,
and it is easy to choose τ 2 small enough such that

|R k (t )| < β/2 for t ∈ (t k , t k + τ 2 ).
Remark 5. The proof of Theorem 2 uses in a crucial way the hypothesis that r is Lipschitz. If r experiences jumps as will be the case in our numerical experiments, it is no more possible to give a lower bound of the delay between two samples. Indeed, although we can estimate | V (e(t ))| by 2 P B K e(t ) ∆ k x(t ) + 2 P BG e(t ) ∆ k r (t ) which can also be bounded by a constant on the time interval [t k , t k + τ 1 ], we do not master the times at which jumps in r take place. Typically if it happens very shortly after t k , when V (e(t )) is very close to δ, the time for V (e(t )) to reach δ again can be very short and cannot be estimated in terms of the parameters of the system.

IV. SIMULATION RESULTS

A. Yaw damper example

As an application to our approach, we consider the example introduced in [START_REF]Yaw damper design for a 747® jet aircraft[END_REF]. It is a simplified model of a jet aircraft during cruise flight, with two inputs: rudder and aileron deflections; and two outputs: yaw rate and bank angle. The state-space representation is given in the form of system [START_REF] Årzén | A simple event-based PID controller[END_REF], where

A =    
-0.0558 -0.9968 0.0802 0.0415 0.5980 -1150 -0.0318 0 -3.0500 0.3880 -0.4650 0 0 0.0805 1.0000 0

    , B =    
0.0729 0.0000 -4.7500 0.00775 0.15300 0.1430 0 0

    , C = 0 1 0 0 0 0 0 1 .
We want the outputs to track a certain reference signal r (t ) ∈ R 2 , which consists in two step signals that start at some value and remain constant for half the simulation time, then both jump at half-time to a different value. So, through pole-placement, we design a state feedback control of the form (4) with For a maximum tolerance on the error of = 0.5, and the selected Lyapunov function, we find δ = 0.1233. Fig. 3 shows the evolution of the two outputs with respect to time, with each output represented with the reference signal it is supposed to track. We can see that for a relatively large value of δ, the system is able to track the reference, even though the first output y(1) exhibits some ripples. This is due to the large value of δ which results in the fact that the individual updates of the control law are rather spaced in time as shown in Fig. 5. However, from Fig. 4, it is obvious that the overall behavior of the system is similar to that of the continuously-controlled reference system. This also shows that the large overshoot experienced by the first output is due to the dynamics of the system and not the event-based scheme.

K = 0.
From Fig. 5 we can see that the control updates are unevenly spaced in time. Indeed, during the transient period or when there is an abrupt change in the value of the reference, the control is updated quite often, whereas when the reference settles at a constant value, the updates become less frequent and rather regular in time. Fig. 6 shows the evolution of the Lyapunov-like function V (e(t )) with respect to time. We can notice that it remains enclosed within the region bounded from above by δ. Fig. 7 is a zoom on the graph of V (e) around t = 15, when r jumps to a new value. It can be noticed that V (e) trespasses the boundary set by δ in some instants. This is due to the fact that computer simulations run in discrete-time rather than continuous-time, and therefore, detecting the time at which V (e(t )) reaches the value δ is not possible.

B. Comment on the Discrete Implementation of the System

For obvious practical reasons our simulations are performed using an underlying discrete time. For the above simulation this sampling time is t s = 10 -4 . Besides a numerical scheme is used to solve the differential equation, namely a second order Adams-Bashforth scheme (two-step explicit scheme).

For the sake of simplicity let us assume that we use a simple explicit Euler scheme instead. Let t n , n ∈ N be the uniform time sampling at which the computation is

(A -B K )e n -B K ∆ k x n + BG∆ k r (t n ), t n ∈ (t k , t k+1 ). (20)
Hence, using the Hurwitz property on A -B K , and that in the worst case e n = δ/λ max P e n+1 ≤ (1

+ t s ) δ/λ max P + t s B K L x (t n -t k ) + t s BG B r .
This yields an upper bound on the worst value of V (e n+1 ) ≤ λ max P e n+1 2 , which can be bigger than δ. This explains the fact that the Lyapunov function exceeds δ on Fig. 7. We cannot avoid such overflows with explicit numerical schemes. When this happens, t n+1 is chosen as t k+1 , and for the computation of e n+2 equation (8b) is used. Then V (e n+2 ) = V (e n+1 )t s e n+1 T Qe n+1

+ t 2 s V ((A -B K )e n+1 ), (21)
Which can still be above δ. A few iterations can be necessary to recover V (e n ) < δ.

C. Comparing the number of control updates for different values of δ

The number of updates of the control law in a given time interval depends directly on the value of δ. Selecting a rather large δ implies that we are allowing a large tolerance on the value of the tracking error e(t ), and consequently, the control law will not need to be updated very often. On the other hand, if we require a smaller tolerance on the error e(t ), the control law will have to be updated more frequently, and this is associated with a small δ. The latter case converges to the continuous-control scheme if we require zero tolerance on the tracking error. TABLE I illustrates the variations of the number of updates in the control signal for the previous example as we change δ. We run the simulation for a fixed time-span in each case (T final = 30). We have run the simulation at a time-step of 10 -4 , thus ending up with 300, 000 simulation instants.

The table shows that even for very small values of and δ, the number of necessary updates in the value of the control remains relatively small compared to the overall number of simulation instants. For example, for = 0.05, the number of updates is only 304 which is still very small for 300, 000 simulation instants, or a ratio of 1/1000. However, the quality of the tracking is much improved for such value of the tolerance and its corresponding δ, as depicted in Fig. 8.

V. CONCLUSION

In this work, we have proposed an event-based implementation for an output tracking problem. We have shown that if we update the control only when the value of a defined Lyapunov-like function reaches an upper limit, we can guarantee a good tracking of the reference signal, with a relatively small number of updates. We can thus reduce the communications between the controller and process as well as the number of operations performed by the controller.

As further work, it would be interesting to consider a variable value of δ. Indeed, if we assume that we have a prior knowledge of the reference input, we can adapt the value of δ accordingly, allowing a larger tolerance when r experiences abrupt changes, and decreasing the tolerance when r is smoother, thus allowing more time for the controller to adapt to the changes in r .

Another idea to explore would be the self-triggered approach. In event-triggered control, even though the control is not updated continuously, the triggering condition is monitored continuously. On the other hand, in a selftriggered scheme, at each update instant we compute not only the control, but also the next instant at which the control would have to be updated, removing thus the obligation to monitor a condition continuously.

Fig. 1 .

 1 Fig. 1. Schematic of the proposed event-based tracking controller.

Fig. 2 .

 2 Fig. 2. Example of the evolution of the Lyapunov equation under our triggering conditions.

Fig. 3 .Fig. 4 .Fig. 5 .

 345 Fig. 3. The evolution of the two output signals along with their respective reference signals

Fig. 6 .

 6 Fig. 6. The evolution of the Lyapunov function

Fig. 7 .

 7 Fig. 7. Behavior of the Lyapunov function when a jump in r occurs

Fig. 8 .

 8 Fig. 8. Output signals for = 0.05 and δ = 0.0012

TABLE I .

 I NUMBER OF UPDATES OF THE CONTROL SIGNAL FOR DIFFERENT VALUES OF AND δ

		δ	Number of Updates
	0.8	0.316	66
	0.5	0.1233	71
	0.05	0.0012	304
	0.001	4.93 × 10 -7	19088

ACKNOWLEDGMENT

This work has been partially supported by the LabEx PERSYVAL-Lab (ANR-11-LABX-0025-01).