Frédéric De Gournay

Jonas Kahn

Léo Lebrat

4 -discrete optimal transport

This paper deals with the 3 4 -discrete 2-Wasserstein optimal transport between two measures, where one is supported by a set of segment and the other one is supported by a set of Dirac masses. We select the most suitable optimization procedure that computes the optimal transport and provide numerical examples of approximation of cloud data by segments.

Introduction

The numerical computation of optimal transport in the sense of the 2-Wasserstein distance has known several breakthroughs in the last decade. One can distinguish three kinds of methods : The first method is based on underlying PDEs [START_REF] Benamou | A computational fluid mechanics solution to the monge-kantorovich mass transfer problem[END_REF] and is only available when the measures are absolutely continuous with respect to the Lebesgue measure. The second method deals with discrete measures and is known as the Sinkhorn algorithm [START_REF] Cuturi | Sinkhorn distances: Lightspeed computation of optimal transport[END_REF][START_REF] Benamou | Iterative bregman projections for regularized transportation problems[END_REF]. The main idea of this method is to add an entropic regularization in the Kantorovitch formulation. The third method is known as semidiscrete [START_REF] Mérigot | A multiscale approach to optimal transport[END_REF][START_REF] Lévy | A numerical algorithm for l2 semi-discrete optimal transport in 3d[END_REF][START_REF] De Goes | Blue noise through optimal transport[END_REF] optimal transport and is limited to the case where one measure is atomic and the other is absolutely continuous. This method uses tools of computational geometry, the Laguerre tessellation which is an extension of the Voronoi diagram. The aim of this paper is to develop an efficient method to approximate a discrete measure (a point cloud) by a measure carried by a curve. This requires a novel way to compute an optimal transportation distance between two such measures, none of the methods quoted above comply fully to this framework.

Pre-existing methods. The first method, the PDEs formulation, requires the measures to be absolutely continuous, which is not the case.

The Sinkhorn algorithm can be applied to this problem if the curve is sampled by points. But to the best of our knowledge [START_REF] De Gournay | Approches variationnelles pour le stippling: distances l2 ou transport optimal ?[END_REF], taking Dirac masses along anisotropic "objects", here, a curve may dwindle the efficiency of this algorithm. Indeed, in this particular case the parameter of regularization of the Sinkhorn algorithm has to be chosen smaller than the curve-sampling precision, which causes numerical issues.

The semi-discrete optimal transport is the more favorable setting for such problem. In a previous paper [START_REF] De Gournay | Optimal transport approximation of measures[END_REF] we described method to approximate an image by a curve. The method developed in [START_REF] De Gournay | Optimal transport approximation of measures[END_REF] relies on a sampling of the curve by points and computes the Wasserstein distance between a density and a point cloud. This approach suffers from several drawbacks, the first one is that solving the optimal transport problem is more and more difficult as discretization step decreases. Indeed the closest the Dirac masses are taken along the curve the more stretched the Laguerre tessellation is. The second flaw is the poor quality of approximation of a curve by a sum of Dirac masses. The third problem is the need to transform the discrete measure into a measure with regular density.

Using ideas from the semi-discrete optimal transport, a solution to the above objections is to compute the optimal transport between an absolutely continuous measure and a measure carried by segments. The main tool of such computations is an extension of the Laguerre tessellation with conic boundaries. The difficulty arises when it comes to integrates the continuous density over theses Laguerre cells. In fact, a robust algorithm with exact integration is complex to develop. This is the reason why computer graphic community [START_REF] Hoff | Fast computation of generalized voronoi diagrams using graphics hardware[END_REF][START_REF] Hiller | Beyond stipplingmethods for distributing objects on the plane[END_REF] implements numerical approximations for those Voronoi cells (Laguerre cells with equals weights) using shape primitives and rasterisation 1 with graphic hardware. However numeric precision is intrinsically tied to grid (pixel). Hence, refinement for this method scales badly when the dimension increases.

Based on these observations we prefer to approximate the measure supported by curve by a measure supported by segments and to compute the optimal transport between an atomic measure and a measure supported on a set of segments. We coined this problem the 3 4 -discrete optimal transport.

Contributions. This paper provides and studies an efficient algorithm to compute the 2-Wasserstein distance between a discrete measure and a measure supported by a set of segments. This algorithm is scalable in 2D and 3D and is parallelized. We also study the problem of optimizing the parameter of the measure carried by the set of segments and provide the formula of the gradient with respect to the parameter of the measure. The 3 4 -discrete optimal transport benefits from the strength of Laguerre tessellation. Moreover it avoids the integration of the density over convex polygons (2D) or polyhedrons (3D), a common difficulty of semi-discrete optimal transport. Indeed the integrations needed for the computation of the cost function boils down to computing intersections between polyhedrons and to integrate the moments of ν over segments.

The ideas of this paper owes a much to [START_REF] De Goes | Blue noise through optimal transport[END_REF][START_REF] Mérigot | A multiscale approach to optimal transport[END_REF][START_REF] Lévy | A numerical algorithm for l2 semi-discrete optimal transport in 3d[END_REF] and the semi-discrete approach but it considers a measure supported by a set of segments instead of a measure with regular density. This modification annihilates the convergence theory developed in [START_REF] Kitagawa | Convergence of a newton algorithm for semi-discrete optimal transport[END_REF]. The optimal transport plan is no longer unique and Kantorovitch functional is no longer C 2 . But under slight geometric condition (see (H)) one can show that the dual functional is C 1 with respect to the dual variables.

Outline of the paper. Section 1 is devoted to set up the notations and the known results used in this paper. In Section 2 we compute the different derivatives of the cost function, in Section 3 we discuss the effective numerical implementation of the algorithm. In Section 4 we settle for the optimization procedure that solves the optimal transport, and in Section 5 we showcase numerical approximations of cloud data by segments.

Setting

Reminders on 2-Wasserstein distance

The 2-Wasserstein distance is a special instance of the optimal transport distance between two probabilities measures.

It is defined as follow:

Given Ω ⊂ R d , µ ∈ P(Ω) and ν ∈ P(Ω), the 2-Wasserstein distance between µ and ν, W 2 (ν, µ) is given by:

W 2 2 (ν, µ) = inf γ∈Π(ν,µ) Ω 2 x -y 2 2 dγ(x, y), (MK)
where Π(ν, µ) is the set of coupling between µ and ν, that is, the set of measures whose marginals are µ and ν :

γ ∈ Π(ν, µ) ⇔          Ω ψ(x)dν(x) = Ω 2 ψ(x)dγ(x, y) ∀ψ ∈ L 1 (ν)
and

Ω φ(y)dµ(y) = Ω 2 φ(y)dγ(x, y) ∀φ ∈ L 1 (µ) (1)
If both (Ω, µ), (Ω, ν) are Polish spaces, an elegant way [START_REF] Villani | Optimal transport: old and new[END_REF] to solve (MK) is via its dual sup

ψ∈L 1 (ν),φ∈L 1 (µ) Ω ψdν + Ω φdµ (2)
s.t ∀(x, y) ∈ Ω 2 : ψ(x) + φ(y) ≤ x -y 2 2 ,
where φ and ψ are the Lagrange multipliers for (1), the marginals constraints of (MK). Introducing the c-transform of φ as :

φ c (x) = inf y∈Ω x -y 2 2 -φ(y).
The problem (2) can be rewritten as :

sup φ∈Φ c (Ω) Ω φ c dν + Ω φdµ, (3)
where Φ c (Ω) is the set of c-concave function on Ω, see [START_REF] Villani | Optimal transport: old and new[END_REF]. Consider now the case where the measure µ is atomic :

µ(x) = n i=1 m i δ xi , x i ∈ R d , m i ∈ R.
In this case φ belongs to R n , and

φ c (x) = min i∈ 1,n x -x i 2
2 -φ i . This naturally leads to the definition of the i-th Laguerre cell [START_REF] Aurenhammer | Power diagrams: properties, algorithms and applications[END_REF] :

L i (x, φ) = {x ∈ Ω such that x -x i 2 2 -φ i ≤ x -x j 2 2 -φ j ∀j ∈ 1, n }. (4)
Provided that ν(L i ∩ L j ) = 0 for every i = j, the final problem states as :

W 2 2 (µ, ν) = sup φ∈R n g(φ, x), with g(φ, x) = i Li(x,φ) x -x i 2 2 -φ i dν + i m i φ i . (5)
Denoting φ as the solution of maximization problem (5) the physical interpretation of the Laguerre cell L i (x, φ ) is that the Dirac mass located at x i is transported to supp(ν) ∩ L i . There is equality between problem (MK) and problem (5) since the quadratic cost is continuous and (Ω, ν),(Ω, µ) are two Polish spaces [START_REF] Villani | Optimal transport: old and new[END_REF].

Setting

In this paper µ will denote a n-atomic probability measure on R d :

µ(x) = n i=1 m i δ xi , s.t i m i = 1 with x i ∈ R d , m i ∈ R + , for all 1 ≤ i ≤ n.
By contrast, ν will be a probability measure supported by a polyline,

ν = p α=1 ρ α l α # λ [0,1] s.t. l α : [0, 1] → R d t → (1 -t)P α + tP α+1 , ρ α ∈ R + , P α ∈ R d , (6)
where the measure l α # λ [0,1] is the push-forward through l α of λ [0,1] , the [0, 1]-Lebesgue measure. The measure l α # λ [0,1] is defined as :

l α # λ [0,1] (A) := λ [0,1] (l α ) -1 (A)
for every Borelian set A. The fact that ν is a probability measure translates into α ρ α = 1. Notice that ρ α = 0 is admissible and in this case the support of the resulting probability ν will be composed of disjoint polylines.

In order to ensure the regularity of function g defined in [START_REF] Beygu | An interacting galaxy system along a filament in a void[END_REF] we enforce hypothesis (H) throughout the paper.

∀α ∈ 1, p , ∀(i, j) ∈ 1, n 2 , i = j : P α+1 -P α , x i -x j = 0. (H)
Under (H), for all φ we have ν(L i (x, φ) ∩ L j (x, φ)) = 0. Indeed if L i and L j are two Laguerre cells with a common boundary L i ∩ L j , the boundary has to be orthogonal to x i -x j . Hypothesis (H) prevents tiny perturbation of φ to harshly shift the affectation of the segment's mass from one Dirac mass to another, see Figure 1. Under this assumption g is a C 1 function of φ, see [START_REF] De Gournay | Differentiation and regularity of semi-discrete optimal transport with respect to the parameters of the discrete measure[END_REF]. Let I be the set of indices (i, α) such that the i th Laguerre cell meets the α th segment, that is:

I = {(i, α) s.t. L i (x, φ) ∩ [P α , P α+1 ] = ∅} . x i x j P α P α+1
Figure 1: When the hypothesis (H) is violated one can cook up a φ such that the functional g in ( 5) is no longer C 1 with respect to φ. Here a small perturbation on φ i or φ j induces a discontinuity on ∇ φ g.

If the space Ω is convex then the Laguerre cells defined in (4) are convex and L i ∩ [P α , P α+1 ] is a segment. For every (i, α) ∈ I denote t iα s , t iα e the starting (respectively ending) time of the α th segment [P α , P α+1 ] in the i th Laguerre cell :

0 ≤ t iα s ≤ t iα e ≤ 1 and L i ∩ [P α , P α+1 ] = [l α (t iα s ), l α (t iα e )]. (7)
Then the cost function g(φ, x) defined in ( 5) can be re-written as :

g(φ, x) = (i,α)∈I t iα e t iα s r iα (t)dt + n i=1 φ i m i (8) with r iα (t) = ( l α (t) -x i 2 2 -φ i )ρ α
2 Derivatives of the cost function

First order derivative with respect to φ

Hypothesis (H) ensures that ν(L i (x, φ)∩L j (x, φ)) = 0, ∀i = j. A direct application of [START_REF] De Gournay | Differentiation and regularity of semi-discrete optimal transport with respect to the parameters of the discrete measure[END_REF] shows that g is differentiable with respect to φ and :

∂g ∂φ i (φ, x) = m i - Li(x,φ) dν. (9)

Computation of the second order derivative with respect to φ

If there exists a segment [P α , P α+1 ] that passes through the intersection of at least three Laguerre cells then the functional g fails to be twice differentiable. On the one hand it is impossible to design an hypothesis in the spirit of hypothesis (H) that can prevent such a pathological case to happen during the optimization in φ. On the other hand such a baneful case almost surely never happens. Hence in this section the Hessian is computed without proof of existence by following a cumbersome calculus. Denote (e i ) i the canonical basis. Following the computation of the first derivative given in ( 9), the second order derivative is given by :

∂ 2 g ∂φ i φ j = -lim εi→0 Li(x,φ+εej ) dν -Li(x,φ) dν ε . ( 10
)
Denote n il the outer normal of L i on the facet L i ∩L l , that is n il = x l -xi x l -xi . Denote by δ il the first order approximation of the evolution of the facet L i ∩ L l in the normal direction n il when we change the Lagrange multiplier of the j-th cell, that is, when we change φ into φ + εe j , see Figure 2 (left). If both i and l are different from j then δ il = 0 otherwise δ ij is given by :

y + εδ ij n ij + o(ε) ∈ (L i ∩ L j ) (x, φ + εe j ) ⇔ y + εδ ij n ij -x j 2 2 -φ j -ε + o(ε) = y + εδ ij n ij -x i 2 2 -φ i ( 11
)
⇒δ ij = 1 2 x j -x i . x j x l x k x m δ jk δ jl δ jm L m L j L k L l P α+1 P α P α-1 n jk δ jk ∂g ∂φj ∂φm ∂g ∂φj ∂φk θ jk α δj m nj m θ jm k-1 L m L j L k
∂g ∂φ i φ j (φ, x) =          - θ∈Θij ρ α 2 x j -x i cos(θ ij α ) if i = j, - k =i ∂g ∂φ i φ k otherwise. ( 12
)
2.3 Computation of the first order derivative with respect to P and ρ

In this section we compute the derivatives of the 2-Wasserstein distance with respect to the parameters of the measure ν. We denote generically as ∂ ν the derivative with respect to one of the parameters that define ν, that is either one of the positions P α or one of the densities ρ α . Let φ be the optimal Lagrange multiplier in ( 5)

G(ν) := g(φ , ν) = W 2 2 (µ, ν)
The derivative of g(φ , ν) with respect to the parameters of ν are given using the chain rule formula.

dG dν = ∂g ∂φ ∂φ ∂ν + ∂g ∂ν ( 13
)
Since φ is a solution of ( 5), the derivative of g with respect to φ is zero at φ . Hence computing the differential of G consists in differentiating the cost function g while keeping φ fixed at φ . Differentiating [START_REF] De Goes | Blue noise through optimal transport[END_REF] with respect to ν we obtain :

∂ ν G = (i,α)∈I (∂ ν t iα e )r iα (t iα e ) -(∂ ν t iα s )r iα (t iα s ) + t iα e t iα s ∂ ν r iα (t)dt
When considering ∂ ν t iα only three cases can occur :

1. If P α belongs to the interior of the ith Laguerre cell then, t iα s = 0 and ∂ ν t iα s = 0. 2. If P α+1 belongs to the interior of the ith Laguerre cell then, t iα e = 1 and ∂ ν t iα e = 0. 3. For all i there exists exactly one j such that t iα s = t jα e . The segment [P α , P α+1 ] intersects L i ∩ L j at point l α (t iα s ). For all x in L i ∩ L j , by definition of the Laguerre cell (4), the following equality holds x -

x i 2 2 -φ i = x -x j 2 2 -φ j then r iα (t iα e ) = r jα (t jα s ). It follows that ∂ ν G = (i,α)∈I t iα e t iα s ∂ ν r iα (t)dt (14)

Derivative with respect to P

We are first interested in differentiating the Wasserstein distance with respect to the position of the endpoints (P α ) α∈ 1,p of the polyline. A direct application of ( 14) yields :

∂G ∂P α = t i(α-1) e t i(α-1) s ∂ Pα r i(α-1) (t)dt + t iα e t iα s ∂ Pα r iα (t)dt
The differential of r iα (t) with respect to P α amounts to differentiate l α (t) -x i 2 , we recall that :

           ∂ l α (t) -x i 2 ∂P α = 2(1 -t) (l α (t) -x i ) ∂ l (α-1) (t) -x i 2 ∂P α = 2t l (α-1) (t) -x i .
Let K(α) the set of indices i such that the i-th Laguerre cell has a non-empty intersection with the α-th segment:

K(α) = {i s.t. L i (x, φ) ∩ [P α , P α+1 ] = ∅} = {i such that (i, α) ∈ I}.
A straightforward computation yields

∂G ∂P α =   P α + P α-1 -P α 3 - i∈K(α-1) t i(α-1) e t i(α-1) s tx i dt   ρ α-1 +   P α + P α+1 -P α 3 - i∈K(α) † t iα e † t iα s tx i dt   ρ α ( 15
)
where † t is the reverse parameterization of the segment [P α , P α+1 ], that is † t iα s = 1 -t iα e and † t iα e = 1 -t iα s . We may interpret [START_REF] Hiller | Beyond stipplingmethods for distributing objects on the plane[END_REF] as the sum of the torques of the surrounding segments [P α-1 , P α ] and [P α , P α+1 ] around the point P α .

Derivative with respect to ρ

The functional g (8) is linear in ρ α , and its derivative is therefore trivial. Throughout this paper, we consider that the density is constant by segment. As a consequence

ρ α = P α -P α+1 α P α -P α+1 , (16)
hence the derivative of ρ with respect to P is easily computed.

Optimization algorithm

The aim of this section is to describe the algorithm which optimizes the Wasserstein distance G with respect to the position of the endpoints P α of the polyline. The proposed algorithm is a gradient descent method with a metric Σ, we recall that the gradient of G with respect to P is given by :

∇ P G = ∂G ∂P + ∂G ∂ρ ∂ρ ∂P , (17)
The differential of G with respect to P and ρ are discussed in Section 2.3.1 and Section 2.3.2 respectively. The main goal of this section is to discuss the choice of the metric Σ.

Let us consider an isolated segment [P α , P α+1 ], that is ρ α-1 = ρ α+1 = 0. The average of the gradient for the segment [P α , P α+1 ] is given by :

1 2 ∂G ∂P α + ∂G ∂P α+1 = 1 2 P α + 1 2 P α+1 - i∈Kα t iα e t iα s x i dt ρ α = ρ α (c α -xα ) (18)
where c α is the center of the α-th segment and xα is the weighted average of points x i seen by the α-th segment :

c α = P α + P α+1 2 and xα = i∈Kα t iα e t iα s x i dt.
Equation ( 18) is reminiscent of the formula of the derivative of G in the semi-discrete setting, see [START_REF] Lévy | A numerical algorithm for l2 semi-discrete optimal transport in 3d[END_REF][START_REF] De Goes | Blue noise through optimal transport[END_REF][START_REF] Mérigot | A multiscale approach to optimal transport[END_REF]. In this setting, µ is a measure with density and the approximating measure ν is a sum of Diracs :

ν = p α=1 m α δ yα
In this case, the derivative of G with respect to the position of the points is given by :

∂G ∂y α = m α (y α -b α ), ( 19
)
where b α = 1 mα Lα xdµ is the barycenter of the αth Laguerre cell. In the semi-discrete setting, the most commonly used algorithm when minimizing G with respect to y is to update the point position y α to their barycenters b α . This procedure is known as Lloyd's algorithm [START_REF] De Goes | Blue noise through optimal transport[END_REF][START_REF] Mérigot | A multiscale approach to optimal transport[END_REF][START_REF] Lévy | A numerical algorithm for l2 semi-discrete optimal transport in 3d[END_REF]. In view of the formula of the gradient [START_REF] Mérigot | A multiscale approach to optimal transport[END_REF], this method is a gradient descent with metric Σ = diag(m α ). In our case it is natural to consider a metric Σ defined via ρ α as an analogy to the semi-discrete setting. We define :

Σ = diag( ρ α-1 + ρ α 2 ). ( 20
)
The main algorithm that minimizes G with respect to P is given by algorithm 1 below Algorithm 1 Optimization polyline position 1: procedure Optimization in P 2:

∇ P G ← 0 3: do 4: [START_REF] Nesterov | Introductory lectures on convex optimization: A basic course[END_REF], ∇ P G defined in [START_REF] Kitagawa | Convergence of a newton algorithm for semi-discrete optimal transport[END_REF] 5:

P ← P -Σ -1 • ∇ P G with Σ defined in
ρ ← normalizedDensity(P ) as in [START_REF] Hoff | Fast computation of generalized voronoi diagrams using graphics hardware[END_REF] 6:

φ * , ∇ P G ← computeOptimalTransport(x, m, P, ρ) see Algorithm 2 while ∇P ∞ ≥ 10 -3 8: end procedure

Numerical implementation

The Laguerre cells are computed using the computational geometry library CGAL [6]. In 2D the algorithm is fast and scalable, the average complexity for n sites randomly drawn is linear in time and memory. For 3D triangulation the worst case complexity is quadratic, but for random point configurations the complexity is observed to be almost linear [START_REF] Dwyer | Higher-dimensional voronoi diagrams in linear expected time[END_REF]. In this section we discuss the computation of g and its parallelization.

Integration computation

The main issue when computing g and its derivative is to evaluate the intersections [START_REF] Cuturi | Sinkhorn distances: Lightspeed computation of optimal transport[END_REF] between the polylines and the Laguerre cells. The intersections are computed by following each segment l α (t) with t increasing. We first focus our attention on computing the exit time of the Laguerre cell j knowing the starting time t iα s . To clarify things suppose that at the known time t iα e = t jα s the segment exits the Laguerre cell L i and enters the Laguerre cell L j . The objective is to compute the time t jα e = t kα s and the index k such that the segment exits the Laguerre cell L j at time t jα e and enters the Laguerre cell L k . Such a time is computed by solving the following minimization problem

k = argmin m s.t. m =i, t jα s <tm≤1 t m , with t m = 2 P α , x j -x m + x m 2 -φ m + φ j -x j 2 2 2 P α -P α+1 , x j -x m . (21)
The exit time t jα e is now equal to t k . Note that the choice of m can be restricted to the indexes of adjacent Laguerre cells, which is a small set in practice (a dozen of indexes). If the set {m = i, t jα s < t m ≤ 1} is empty, the segment ends in the Laguerre cell L j and we set t jα e = 1 and we stop the procedure for the segment. In the case when t jα s = 0, that is the segment starts in the Laguerre cell L j , we consider the same minimization problem as [START_REF] Polyak | Introduction to optimization[END_REF] where the constraint set is replaced with {m, 0 ≤ t m ≤ 1} .

We still have to compute the index of the Laguerre cell where the segment [P α , P α+1 ] begins. If α > 0, it is obviously the index of the ending Laguerre cell of [P α-1 , P α ]. For the case α = 0, we add a dummy segment [x l , P α ], where x l is the position of the l-th Dirac mass corresponding to the largest multiplier φ l = max k φ k . By definition of the Laguerre diagram (4) the point x l belongs to L l .

Given the starting and ending times, the computation of g(φ, x) and its derivative abridge to integrate polynomials within [t jα s , t jα e ] using Gaussian quadrature. There might exist several solutions k to the minimization problem [START_REF] Polyak | Introduction to optimization[END_REF]. In this case, the segment encounters a corner, the intersection of at least 3 Laguerre cells. In practice this case never occurs thanks to floating point arithmetic but the rounding error can elect a non-suitable candidate. In this case, the algorithm assigns a segment of negligible length to the candidate, see Figure 3.

Parallelism

The evaluation of g(φ, x) and its derivative with respect to φ are discussed in Section 3.1. For each segment of the polyline the intersection times are computed sequentially by [START_REF] Polyak | Introduction to optimization[END_REF] and the integration is performed on the fly. The

L i L l L j P α P α+1 L i L l L j ε L i L l L j
Figure 3: Example of non-uniqueness to the minimization problem [START_REF] Polyak | Introduction to optimization[END_REF]. At the left the polyline [P α , P α+1 ] intersects the junction of 3 Laguerre cells (i, j, l). In the center k = l chosen as expected. At the right k = j chosen and t jα e -t jα s = ε, a small error occurs assigning a bit of the segment to the jth Laguerre cell.

parallelism is unfurled at the highest level of the algorithm for the segments. Remember that the computation of the intersection of the polyline and the Laguerre tessellation requires for each segment the knowledge of the cell of its starting point. This cell can be inferred from the calculation of the previous segment. Hence, we slice the polyline into contiguous chunks of segments with equivalent size, each worker dealing with one chunk. For each of these chunks the sequential integration procedure is launched and then the results are merged. We carry out both OpenMP and c++11 thread class implementation. For a standard chunk size and for a regular processor the c++11 performance overtakes the OpenMP's one. We suppose that this difference of performance (up to a factor two) is due to OpenMP overhead.

In Figure 4, the evolution of the computation time of the cost function g and its gradient with respect to φ is benchmarked. The speedup unit is defined as the ratio of the execution time of the non-parallelized task over the parallelized task. OpenMP and c++11 thread class are both implementations of shared memory parallelization, the memory is simultaneously accessible for every thread. The performance of a process depends intrinsically on how close the data is. If the data is on processor cache the latency is a dozens of CPU cycles, if the data is located on RAM, the latency is 40ns. In our experiment, for large number of threads, 30% of the data was not located on cache. This explains why performance drops as the number of threads increases.

Computation of the optimal transport

The goal of this section is to compare the available methods to optimize g with respect to φ (5). The functional g is concave with respect to φ since it is a dual formulation of the problem (MK) see [START_REF] Villani | Optimal transport: old and new[END_REF]. Under the hypothesis (H) the functional g is C 1 with respect to φ. Note that a higher level of regularity cannot be established as argued in Figure 5. This lack of regularity precludes the use of the convergence framework for second order methods established by Mérigot et al. [START_REF] Kitagawa | Convergence of a newton algorithm for semi-discrete optimal transport[END_REF].

However, we may apply classical convergence results of concave analysis to functional the g. To trigger those arguments, the Lipschitz constant of the gradient L has to be bounded. Under Hypothesis (H), using the notation of Section 2.2 summarized in Figure 2 and the expression of the Hessian matrix [START_REF] Dwyer | Higher-dimensional voronoi diagrams in linear expected time[END_REF], we establish a finite upper bound for L using the Gershgorin circle theorem : 4: Speedup for the computation of the cost function g with c++11 thread class. The parallelism is placed over the segment integration for n = 200K points and p = 80K segments for an increasing number of threads with a super-calculator equipped with Intel Xeon R E5-2680

|L| ≤ 2 max i j =i ∂ 2 g ∂φ i ∂φ j ≤ p 2 max α (ρ α ) min l =m ( x l -x m 2 ) min i,j,α,i =j | cos((θ ij α ))| , (22)
x i P α-2 P α-1 P α+2 P α+1 P α φ i ∇g(φ i )
Figure 5: Counter-example to the smoothness of ∇g. Consider the a Laguerre tessellation and a set of segment as displayed in the left and increase φ i . The Laguerre cell L i increases and the gradient of g with respect to φ exhibits kinks. Those kinks happen each time a Laguerre cell meets a segment for the first time. Note that under Hypothesis (H), the boundary of the Laguerre cell can only cross transversally the segment. If Hypothesis (H) is not met ∇g may fail to exist.

where p is the number of segments composing the polyline, ρ α is the density associated to each segment and θ ij α is the angle between the vector P α+1 -P α and (x j -x i ). For fixed measures ν and µ the evaluation of the bound ( 22) is costly as it involves combinatorial quantities and is not evaluated in practice.

The gradient Lipschitz condition ( 22) is sufficient to ensure convergence of ascent methods. Note that g is bounded from above as µ and ν are compactly supported measures. The gradient method with step s < 2 L converges to a stationary point φ see [START_REF] Polyak | Introduction to optimization[END_REF][START_REF] Bertsekas | Nonlinear programming[END_REF][START_REF] Nesterov | Introductory lectures on convex optimization: A basic course[END_REF]. The same holds for variant step-size method with line search such as Global Brazalai Borwein algorithm [START_REF] Raydan | The barzilai and borwein gradient method for the large scale unconstrained minimization problem[END_REF][START_REF] Fletcher | On the barzilai-borwein method[END_REF].

Because of the counter-example in Figure 5, and the lack of regularity of the gradient of g, quadratic convergence cannot be guaranteed for second order methods. Note also that in a generic setting, the Newton method is impractical since the Hessian matrix defined in equation ( 12) fails to be invertible, see Section 4.1.3.

Choice of optimization method

The different methods are tested against the same benchmark. It consists in drawing uniformly 10K points and 500 segments in 2D. Several realization of the optimization procedure are performed and plotted in the corresponding figures. The methods are benchmarked in Table 1.

First order method

We first implement a gradient ascent method. The generic convergence history is displayed in Figure 6. As explained in Section 4, the Lipschitz constant of the gradient cannot be satisfactorily computed, so the step-size is chosen according to strong Wolfe conditions. In practice the algorithm settles for a constant step of 0.05. As it can be observed in Figure 6 the rate of convergence of the gradient towards zero is too slow for this method to be used in practice. We have tried other first order methods including Polak-Ribire, Fletcher-Reeves, Barzalai Borwein and Nesterov acceleration. Only Nesterov acceleration has behaved differently from the other methods and it is displayed in Figure 7. Note however, that Nesterov acceleration requires an estimation of L, the Lipschitz constant of the gradient. Figure 7 was obtained with an estimation of L that promotes quick decreases of gradient norm at the beginning of the algorithm. Note however that the estimation of L seems too optimistic since the Nesterov algorithm does not converge. A more pessimistic choice of L leads to a convergence rate similar to the gradient algorithm, at least during the first 7K iterations.

Quasi-Newton Method

The main idea behind quasi-Newton methods is to build at the iteration k an approximation of the Hessian matrix of g. We choose the limited memory BFGS method, which only stores a limited amount of vectors determined by the user. The result is displayed in Figure 8. This method converges faster than first order method, the L-BFGS

Newton and quasi-Newton hybridization

In the light of the previous section, the flaw of the quasi-Newton method is its lack of briskness when it reaches the vicinity of a maximizer. Second-order methods are known to converge quadratically in the basin of attraction of a maximizer. The idea of the hybrid algorithm is to start with a quasi-Newton method, and then switch to a second order Newton algorithm when the basin of attraction is reached. The quandary is to determine when to start the Newton method in preference to the L-BFGS method. Indeed far from the basin of attraction the direction given by the Newton method is poor and the natural step 1 is truncated by the line search which causes extra computation time.

We propose to switch from the L-BFGS algorithm to the Newton algorithm as soon as there is no empty Laguerre cell. A Laguerre cell L i is called empty as long as there is no mass assigned to its centroid x i that is ν(L i ) = 0. This criterion is informally motivated by the fact that the second order information relates the competition and the connectivity between Laguerre cells, hence if one of them is empty its information is inconsequential. From a practical point of view, this condition ensures the Hessian on being a full rank matrix and hence invertible. This procedure is described in Algorithm 2 and the result are displayed in Figure 9.

Algorithm 2 Computation of optimal transport Input: Dirac positions x and masses m. Input: Nodes P and density ρ of the polyline. Require: φ init a starting φ for the computation of g Require: gradTol L2 tolerance on the gradient norm Require: outerMax maximum number of iterations Require: wolfeProcedure standard line search with strong Wolfe condition and initial step s = 1. [START_REF] Bertsekas | Nonlinear programming[END_REF] 1: function computeOptimalTransport(x,m,P ,ρ)

2: bf ← LBFGS(memSize) Initialization of L-BFGS 3: φ ← φ init 4:
∇φ,cost,hiddenNumber ← computeIntegration(φ) H ← computeHessian(φ)

11: d ← -H -1 ∇φ
The Hessian is definite so Newton direction is taken end while 19: end function In the numerical tests, the switch between the two Newton method happens approximatively 30 iterations before termination of the algorithm. During the Newton phase of the hybrid algorithm, the algorithm undergoes two stages. During the twenty first iteration of the Newton method the algorithm stabilizes around the maximum and the Wolfe line-search prevents picking s = 1. During the second stage, quadratic convergence is achieved in few iterations.

Levenberg-Marquartdt's method was also implemented with a regularization parameter that tends to zero when approaching the basin of convergence, however in large scale optimization problems the lightweight hybrid method is preferable since it requires the inversion of the Hessian only in the last iterations. Table 1: Comparison of algorithms for solving the optimal transport problem for 10 4 points and 500 lines. Time required for an iteration in millisecond and the gradient norm after 1000 iterations. The BFGS/Newton method converges up to numerical error.

5 Numerical examples

Representation of picture

The approximation of a measure by a curve has multiple applications of which some are described in [START_REF] De Gournay | Optimal transport approximation of measures[END_REF]. In this section we discuss the representation of a picture by a polyline. The first step is the discretization of the picture, here a landscape, by a sum of Dirac masses of weights m. Several approaches are possible but the most intuitive is to take the Dirac positions x on a Cartesian grid and their weights equal to the pixels intensities. The polyline is then initialized randomly, and Algorithm 1 is launched. With no speed or curvature constraints, the gradient method in P described in Algorithm 1 empirically gives a stationary point. This solution is highly non-smooth, the length of the segments and the angle between them are disparate see Figure 10 top right. This scribbled solution is, in practice, difficult to carry out with laser engravers. Indeed for the same number of segments p, the engraving time of uncontrolled trajectories can be up to five times longer. In order to get around this problem, we follow the method described in [START_REF] De Gournay | Optimal transport approximation of measures[END_REF], that is projecting the polyline after the gradient step of Algorithm 1 on a set of kinematic constraints. The constraints imposed on the speed and the acceleration polyline read as :

K(K 1 , K 2 ) = (P α ) α∈ 1,p s.t P α+1 -P α ≤ K 1 , ∀α ∈ 1, p -1 2P α -P α-1 -P α+1 ≤ K 2 , ∀α ∈ 2, p -1 , (23)
with K 1 the constant controlling the speed of the polyline and K 2 its acceleration. The projection on the set K(K 1 , K 2 ) is performed using the Alternating Direction Method of Multipliers (ADMM) [START_REF] Glowinski | On alternating direction methods of multipliers: a historical perspective[END_REF].

The authors wish to thank Andrew Gibson for letting them use the original picture and Alban Gossard for his priceless help in realizing the wood engraving.

Galaxy filaments

Galaxies are known to cluster along filaments and other low dimensional structures [START_REF] Tempel | Galaxy filaments as pearl necklaces[END_REF][START_REF] Beygu | An interacting galaxy system along a filament in a void[END_REF]. However mathematical extraction of these filaments is a challenging task. We might try to apply refinements of our algorithm to find those filaments. As a naive proof of concept, we have applied our method while setting ρ 2k+1 = 0 for all k, so that the lines are disjoint. We use the data of [START_REF] Tempel | Bisous modeldetecting filamentary patterns in point processes[END_REF] available at : https://github.com/etempel/bisous. The galaxies are represented by Dirac masses. Note that in our tests their masses m are arbitrarily set to 1 n but the code enables other values. In Figure 11, computations are performed for n = 180K galaxies and a decreasing number of filaments 8K,1K,500.

Figure 2 :

 2 Figure 2: (Left) Normal displacements (δ il , δ ik , δ im ) of the boundary of the Laguerre cell when φ i increases and (φ l , φ k , φ m ) remain constant. (Right) geometrical interpretation of the Hessian Denote Θ ij the indices of the segments that intersect the facet L i ∩ L j , it is given by Θ ij = {α such that, (i, α) ∈ I and (j, α) ∈ I} . For every α in Θ ij denote θ ij α the angle between the segment [P α , P α+1 ] and the outer normal n ij . See Figure 2 (right), simple geometrical consideration shows that the Hessian formula boils down to :

Figure 6 :

 6 Figure 6: Gradient ascent method : cost function history (left) and norm of the gradient (right). Method stopped by the max iteration criterion.

Figure 7 :

 7 Figure 7: Nesterov acceleration of gradient ascent : cost function history (left) and norm of the gradient (right). Method stopped by the max iteration criterion.

Figure 8 :

 8 Figure 8: L-BFGS method : cost function history (left) and norm of the gradient (right). Method converges to the targeted gradient norm 5.10 -5 in a average of 1500 iterations.

Figure 9 :

 9 Figure9: Hybrid Newton method : cost function history (left) and norm of the gradient (right). The method converges to the targeted gradient norm 10 -6 in 600 iterations on average.

Figure 10 :

 10 Figure 10: Curvling by 3 4 -discrete optimal transport. Original image (top left), approximation by a polyline with 400K segments without any constraints (top right), approximation by a polyline with kinematic constraints and 250K segments (bottom left), final rendering after wood engraving, the polyline is composed of 80K segments (bottom right). In all tests measure µ is represented by 320K Dirac masses.

Figure 11 :

 11 Figure 11: A thin slice of the global solution for locating Bisous filaments (white), and galaxies (yellow).


Adnan Al Alwani
email: adnan.alalwani@unicaen.fr

Youssef Chahir
email: youssef.chahir@unicaen.fr

Francois Jouen
email: francois.jouen@ephe.sorbonne.fr

Keywords: Events recognition, Non-Redundant Temporal Local Binary Pattern, Thermal signature, Support Vector Machine, Wavelet

In this paper, we propose a method for events recognition from thermal signature based 1-D signal. We use the non-redundant temporal Local Binary Pattern NRTLBP as a descriptor of the Pattern Of Interest (POI) signal. The original signal is extracted directly from local patch in region of interest. We introduce the wavelet decomposition as a pre-processing stage in order to extract the approximation wave-components of the raw signal. Then, NRTLBP is applied on the wave-components which provide wavelet domain descriptor of the raw thermal signature. Finally, we provide an evaluation of our method on the real dataset (Preterm Pain in Infants "PPI" ) composed of thermal videos developed in the context of Infant pain project, a french project

Introduction

Feature extraction is an important challenge in pattern recognition, specifically if the pattern is represented by a temporal or spatial signal such as speech, acoustic signal or 2D signal. The temporal signals have finite durations and are sampled into digital format. In such case, it is more appropriate to represent a pattern as a finite time sequence {x[0], x [START_REF] Chekmenev | Non-contact, wavelet-based measurement of vital signs using thermal imaging[END_REF], ..., x[N -1]}. Providing this sequence directly to a classifier is not typical due to the huge number of inputs and due to the randomness of the signal. Therefore, the sequence x[n] must be encoded into a low-dimension feature vector Y = (y 1 , y 2 , ..., y D ), D << N , which better characterizes the pattern and provides a robust classifier training and classification. Feature extraction from One-Dimension (1-D) signal is generally categorized into two approaches : temporal and spectral approach. Both approaches are deemed as nonparametric approaches. As well as parametric approaches, both approaches have been used in classification of many patterns examples. Spectral approaches have gained a big popularity in the last two decades, which are out of scope of this paper. In contrast, in the temporal approach, the features are extracted directly from the temporal sequence. Moreover, use of multifunction allows to characterize the time sequence by the following quantities : mean value, slope absolute value, number of zero crossings, slope sign changes, and waveform length. Although there are techniques that claimed a successful discrimination in various applications, the above temporal characterizations didn't help much in the other applications. Techniques from the fields of machine learning and pattern recognition might be valuable for thermal pattern examination. Thermal pattern may have sufficient knowledge to understand the underlying characteristics of the patterns that we observe. However, reliable classification of thermal signature would help us to learn about an unknown physical aspects, discriminate the normal and abnormal states in the condition monitoring. The realization of such a tool can be difficult task in some cases. This paper presents an approach for event recognition from temporal thermal signal. The recognition process is based on the Non-Redundant Temporal Local Binary Pattern-based features (NRTLBP) over the time series thermal sequences which are extracted from the thermal video sequences. A traditional way to characterize the properties of a time series sequence x[n] is through its raw temperature value extracted from local patches that are defined over the objects face in image sequences. It has been shown that if a signal has band-limited characteristics, significant improvement in the performance of pattern recognition can be readily made by a relatively simple preprocessing of the signal in the time-frequency domain. The redundant representation of a 1D signal in a 2D time-frequency domain can provide an additional degree of freedom for signal analysis. Such pre-processing effectively separates the intertwined time domain features of the signal, allowing the important characteristics to be exposed in the time-frequency domain, resulting in more effective pattern matching. Hence, the signal being analyzed needs to be conducted only around in selected regions of interest in the time-frequency domain. For this, we use the wavelet transform to decompose the signal into pyramid structure in the form of approximation component and detail component. The approximation wave components are also known to be a robust and distortion tolerant feature space for many pattern recognition applications. Moreover, we consider the approximation component as the feature vector of raw signal, which is a necessary requirement for extracting NRTLBP-based features. As the temporal LBP produces efficient feature descriptor, we employ a Support vector Machine SVM to evaluate our method on the real Preterm dataset. Extracting features from one-dimensional signal based on thermal video monitoring methods, have been considered in a number of approaches. The body work of these approaches is based on measurement, quantifying and tracking of thermal signature. Many other approaches in the field of the thermal signature exam-ination have been proposed in the literature. In [START_REF] Chekmenev | Non-contact, wavelet-based measurement of vital signs using thermal imaging[END_REF] the authors examined the temperature changes around the neck and the nasal regions. These regions were identified manually in the recorded images, then the respiration rate signals were extracted from the thermal image by performing wavelet transform. Respiration signal rate monitoring is also explored in [START_REF] Ding | Derivation of respiratory signal from singlechannel ECGs based on source statistics[END_REF]. A tracking algorithm is developed to follow three facial features in order to produce the respiration signal rate. New non-contact respiration signal modality for neonatal was developed and investigated by [START_REF] Abbas | Neonatal non-contact respiratory monitoring based on real time infrared thermography[END_REF], using infrared thermography imaging. This development includes subsequent image processing (region of interest (ROI) detection) and optimization based wavelet technique. In [START_REF] Farah | An Evaluation of Thermal Imaging Based Respiration Rate Monitoring in Children[END_REF] F. AL-Khalidi et al. developed a technique to model the respiration signal from nasal region. A feature extracted from this region enabled a respiration signal to be produced. The ROI under the tip of the nose was segmented into eight parts. The pixel values within each segment were averaged to obtain a single value representing that segment. Additionally to the direction of vital signal examination. The approach in [START_REF] Sun | Imaging the cardiovascular pulse[END_REF] used mid-wave IR sensors for the measurement of cardiac and breathing rates at a distance. They used the Fourier Transform to isolate the relevant frequency components from skin temperature in the vicinity of the major arteries. Furthermore, they showed that the frequency components from skin temperature is directly related to the pulse waveform related to the heart signal. Our method differs from other approaches by taking the root of pattern recognition to recognize the event from thermal video. We investigate the discrimination power of NRTLBP in the field of the thermal pattern descriptor in Section 3, whereas Section 4 illustrates its performance using wavelet transform to decompose the raw statistical features and build new feature representation. In Section 5, the experimental protocol is recalled and we conduct comprehensive experiments to demonstrate the effectiveness of our approach in recognizing if an infant is in pain. Finally, conclusions are drawn in Section 6.

A Framework for Event Recognition

We propose a framework of using thermal attributes extracted from the video sequences. We assume that the layout of subjects is considered for all the cases in the front view, and ROI is defined over the subjects face. In addition, we build raw thermal signatures for all subjects samples. This is achieved first by defining local patch as shown in the Fig1. Then the maximum M x and the minimum M n values are computed from the local patches along the video sequences. For more illustration, we compute the M x and M n temperature values and defined it as the raw thermal signals (we call it raw signals for abbreviation). Furthermore, we establish the raw signals for the three condition events, which characterize the monitoring state during daily care. For instance, normal event, disease event and post-disease event are three events adopted in our study. The time series of each event is shown in the Fig2. As illustrated in this figure, the response nature of each event specify the discriminative features of that event. Then, we apply the NRTLBP descriptor on the raw thermal signal in order to extract efficient feature vector. Another extension consists to approximate components, before using NRTLB, via wavelet decomposition of the temporal raw signal. A video sequence is then described by two descriptors which are directly applied into SVM for event recognition.

Non-Redundant Temporal Local Binary Pattern-based features

Once the thermal features are extracted, the video sequences is considered as a collection of raw signals. We denote by

X i = T 1 , T 2 , • • • , T n . i = M
x or M n, and T is the temperature value at each patch. In order to represent the raw signal in efficient manner we propose the Temporal Local Binary Pattern descriptor, which is an extension of original 2-D LBP operator into the temporal domain. Normal LBP has originally been proposed for texture analysis and classification [START_REF] Ojala | Multiresolution gray-scale and rotation invariant texture classification with local binary patterns[END_REF]. Recently, it has been applied on face recognition [START_REF] Ahonen | Face recognition with local binary patterns[END_REF] and facial expression recognition [START_REF] Ahonen | Face Description with Local Binary Patterns: Application to Face Recognition[END_REF]. The TLBP operator labels the samples of a signal by thresholding a center sample against neighborhood set within defined window. We denote by x[n] the sampled signal, and w is the size of samples window. The TLBP operator on a sample kw in is given by:

T LBP w (f [k]) = w 2 -1 q=0 (f [k] -f [k -j -q]2 q + f [ w 2 + q + 1] -f [k]2 q+ w 2 ( 1
)
Where q is the number of sampled points (neighbor samples of k) whose the distances to k do not exceed the window w. For a M block, a TLBP histogram of 2 q bins is computed for feature representation. After the computation of the LBP for the whole signal, a histogram of the TLBP code is used as feature descriptor. It contains information about the distribution nature over the whole signal and characterizing a statistical description of signal. Another aspect of TLBP is called no-redundant TLBP. Intuitively, the NRTLBP considers a TLBP code and its complement as same. NRLBP is defined by:

N RT LBP (f [k]) = min(T LBP (f [k]), 2 q -1 -T LBP (f [k])) (2)
This is motivated by the advantage of using dynamic shape of the signal at various instances and levels. Finally, the no-redundant local patterns formed from signal f [k] are formed as a histogram distribution vector.

NRTLBP of Approximated Coefficients

We further extend the idea of LBP into a wavelet domain in order to perform the comparison results. The method is based on using multi-resolution wavelet transform to decompose the raw statistical features and build new feature representation [START_REF] Mallat | A wavelet tour of signal processing, the sparse way[END_REF]. The wavelet transform C(j, k) of a finite-energy signal f (t) is defined as its scalar product with the wavelet Ψ jk (t). In other words, the wavelet transform represents the correlation of the signal f (t) and the wavelet Ψ jk (t) as:

a jk = t f (t)Ψ jk (t). (3)
Where Ψ (t) = 2 j Ψ (2 j t -k) is the mother wavelet corresponds to scale j = 1, k is the translation factor, and j is the scale parameter. A wavelet representation of a function consists of, a coarse overall approximation a j , and detail coefficients that influence the function at various scales. Therefore, the approximation coefficients are considered as raw feature vector in a wavelet-domain. Then using NRTLBP descriptor for feature description. The smooth change property of the approximation features is illustrated in the Fig 4 . Where each subfigure shows the feature in wavelet-domain for the corresponding event in the time domain explained in the Fig. 2.

Experiments

In this section, we evaluate the proposed method by experimenting with the real dataset from ANR pretherm project. Due to non available of the benchmark related to our method, we did not attempt a comparison with other methods. Rather, we provide results only to show that our algorithm works well on a real data. The preterm data involves 20 neonates videos. Each includes the subject behaviors during daily condition monitoring. All events are performed for each subject and notated by the clinical during condition monitoring. From each video we select the set of events which includes normal health event (pre-disease), disease event, and post-disease event. A total of 60 clips were extracted ranging from N occurrences of normal event, to M occurrences of disease event and so on. We divide the clips into three groups of 20 clips. For instance, each group divided into 2 subset of 10 clips each. In the training phase, we use half of the group and the rest was used for the test task. We consider a recognition task for three class problems including the normal event, disease event and post-disease event respectively. Furthermore, two features descriptors are used in the experiments. The first one is the descriptor of raw thermal signal, we denote it by NRTLBP and the second one is the descriptor based on the approximation components, and we denote it by WNRTLBP. The number of neighboring samples was set to 4 samples for NRTLBP codes, and the Daubechies wavelet with one level decomposition was used for approximation coefficients extraction. We performed the evaluation of our method using linear SVM. In all experiments, first subset of each group was used for training and the second subset was used for testing. We started the experiment directly by providing the raw features (Mx, and Min raw features) into classifier. The recognition results are in the form of recognition rate, the results of the first experiment are reported in the table 1. We also run the test by using NRTLBP descriptor based on two types of raw features . The results are shown in table 2. In order to be able to do a full comparison of methods, WNRTLBP based experiments are performed in the same manner as above. Results from this experiment are presented in table 3. As we can see, the approach based on the wavelet achieves a good accuracy rate. From a comparison viewpoint, the results that are achieved based on feature descriptors reported that, NRTLBP and WNRTLBP perform best recognition rates when the maximum value of raw feature is used.

Conclusion

In this paper, we have applied temporal NRLBP as a descriptor of thermal signature from thermal video based on local patch for events recognition. And shown NRTLBP to be suitable for the description of thermal signal containing information about dynamic characteristics of thermal signature. We have extended NRTLBP considering the wavelet transform of feature representation.

We have also shown that the power of descriptors for one-dimensional thermal signature. The experimental results reveal that non redondant temporal LBP and its modifications tend to be good descriptor based approach for thermal features description.

Fig. 1 .

 1 Fig. 1. Local temperature extracted from ROI defined over the object's face.

Fig. 2 .

 2 Fig. 2. Raw thermal signals of Mx and Min temperature values. From top to bottom, normal event, disease event and post-disease event signatures

Fig. 3 .

 3 Fig. 3. Approximation component of raw features for three kinds of events. from top to bottom, normal event, disease event, and post-disease event signatures.

Table 1 .

 1 Recognition rate results, from the Mx and Mn raw features that are directly applied into SVM

	Events	Recognition rate% Recognition rate%
		(Mx Temp.value) (Mn Temp. value)
	Normal Response	50	60
	Disease Response	90	60
	Post-disease response	50	70
	Overall	63.33%	63.33%

Table 2 .

 2 Recognition rate results, using NRTLBP based raw feature descriptor

	Events	Recognition rate% Recognition rate%
		(Mx Temp.value) (Mn Temp. value)
	Normal Response	80	50
	Disease Response	90	90
	Post-disease response	80	90
	Overall	83.33%	76,67%

Table 3 .

 3 Recognition rate results, using WNRTLBP based raw feature descriptor

	Events	Recognition rate% Recognition rate%
		(Mx Temp.value) (Mn Temp. value)
	Normal Response	90	100
	Disease Response	100	95
	Post-disease response	85	93
	Overall	91,667%	96%

supported by the French National Research Agency Projects for science (ANR).