Vincent Pilaud
email: vincent.pilaud@lix.polytechnique.fr

Pierre-Guy Plamondon
email: pierre-guy.plamondon@math.u-psud.fr

Salvatore Stella
email: stella@math.haifa.ac.il

A τ -Tilting Approach to Dissections of Polygons

Keywords: dissections of polygons, accordion complexes, τ -tilting theory, representations of finite dimensional algebras, g-vectors 2010 Mathematics Subject Classification: 16G10, 16G20, 05E10

We show that any accordion complex associated to a dissection of a convex polygon is isomorphic to the support τ -tilting simplicial complex of an explicit finite dimensional algebra. To this end, we prove a property of some induced subcomplexes of support τ -tilting simplicial complexes of finite dimensional algebras.

Introduction

The theory of cluster algebras gave rise to several interpretations of associahedra [START_REF] Stasheff | Homotopy associativity of H-spaces. I[END_REF][START_REF] Stasheff | Homotopy associativity of H-spaces. II[END_REF][START_REF] Tamari | Monoïdes préordonnés et chaînes de Malcev[END_REF]. Fig. 1 shows two such interpretations for the rank 3 associahedron: as the exchange graph of triangulations of a hexagon and as the exchange graph of support τ -tilting modules over the cluster tilted algebra whose quiver with relations is as depicted. This follows from results in the setting of the "additive categorification of cluster algebras" that was initiated in [START_REF] Buan | Tilting theory and cluster combinatorics[END_REF][START_REF] Caldero | Quivers with relations arising from clusters (An case)[END_REF].

F. Chapoton observed a similar isomorphism between the exchange graph of certain dissections of a heptagon and that of support τ -tilting modules over the path algebra of the quiver 3 , subject to the relation βα = 0. Fig. 2 shows these two exchange graphs, which can be found in [START_REF] Chapoton | Stokes posets and serpent nests[END_REF]Fig. 7] and in [START_REF] Adachi | τ -tilting theory[END_REF]Example 6.4].

The purpose of this note is to show that this isomorphism is an avatar of a more general result in the theory of τ -tilting modules. Any basic finite dimensional algebra Λ gives rise to an exchange graph on support τ -tilting modules. This exchange graph is the dual graph of the support τ -tilting complex [START_REF] Adachi | τ -tilting theory[END_REF] (see Section 2). Let {e 1 , . . . , e n } be a complete set of primitive pairwise orthogonal idempotents of Λ. Let J be a non-empty subset of [n] and e J := j∈J e j . The following result forms the algebraic core of the paper (see Section 2 for definitions and Section 3 for the proof).

Theorem 1.1. The support τ -tilting complex of e J Λe J is isomorphic to the subcomplex of the support τ -tilting complex of Λ induced by the support τ -tilting modules whose g-vectors' coordinates vanish outside of J.

D • = 2 1 ⊕ 3 2 ⊕ 1 3 3 ⊕ 3 2 ⊕ 1 3 2 1 ⊕ 1 ⊕ 1 3 2 1 ⊕ 3 2 ⊕ 2 3 ⊕ 3 2 3 2 ⊕ 3 ⊕ 1 3 1 ⊕ 1 3 2 1 ⊕ 1 2 1 ⊕ 2 3 1 2 ∅ Q = 1 2 3 Figure 1.
The exchange graph on triangulations of a hexagon (left) and the exchange graph on support τ -tilting modules of the quiver with relations Q (right). With this algebraic result, we can explain the isomorphism of Fig. 2, and extend it to dissections of any polygon. Any reference dissection of a polygon gives rise to an exchange graph on certain dissections. This exchange graph is the dual graph of the accordion complex studied in [START_REF] Chapoton | Stokes posets and serpent nests[END_REF][START_REF] Garver | Oriented flip graphs and noncrossing tree partitions[END_REF][START_REF] Manneville | Geometric realizations of the accordion complex of a dissection[END_REF] (see Section 4). Theorem 1.2. Any accordion complex is isomorphic to the support τ -tilting simplicial complex of an explicit finite dimensional algebra. Thus, the corresponding exchange graphs are isomorphic.

D • = 1 ⊕ 2 1 ⊕ 3 2 2 ⊕ 2 1 ⊕ 3 2 1 ⊕ 3 2 ⊕ 3 1 ⊕ 2 1 2 ⊕ 2 1 2 ⊕ 3 2 3 ⊕ 3 2 1 ⊕ 3 2 3 1 ∅ Q(D • ) = 1 2 3
Note that this statement could be seen as a consequence of a more general result [13, Proposition 2.44]. However, the latter combines several isomorphisms (from dissections, via non-kissing paths, through non-attracting chords, to support τ -tilting modules) and the proof is quite intricate. In the present paper, we will easily deduce Theorem 1.2 from Theorem 1.1 and from the Figure 1. The exchange graph on triangulations of a hexagon (left) and the exchange graph on support τ -tilting modules of the quiver with relations Q (right).

D • = 2 1 ⊕ 3 2 ⊕ 1 3 3 ⊕ 3 2 ⊕ 1 3 2 1 ⊕ 1 ⊕ 1 3 2 1 ⊕ 3 2 ⊕ 2 3 ⊕ 3 2 3 2 ⊕ 2 3 ⊕ 1 3 1 ⊕ 1 3 2 1 ⊕ 1 2 1 ⊕ 2 3 1 2 ∅ Q = 1 2 3 Figure 1.
The exchange graph on triangulations of a hexagon (left) and the exchange graph on support τ -tilting modules of the quiver with relations Q (right). With this algebraic result, we can explain the isomorphism of Fig. 2, and extend it to dissections of any polygon. Any reference dissection of a polygon gives rise to an exchange graph on certain dissections. This exchange graph is the dual graph of the accordion complex studied in [START_REF] Chapoton | Stokes posets and serpent nests[END_REF][START_REF] Garver | Oriented flip graphs and noncrossing tree partitions[END_REF][START_REF] Manneville | Geometric realizations of the accordion complex of a dissection[END_REF] (see Section 4). Theorem 1.2. Any accordion complex is isomorphic to the support τ -tilting simplicial complex of an explicit finite dimensional algebra. Thus, the corresponding exchange graphs are isomorphic.

D • = 1 ⊕ 2 1 ⊕ 3 2 2 ⊕ 2 1 ⊕ 3 2 1 ⊕ 3 2 ⊕ 3 1 ⊕ 2 ⊕ 2 1 2 ⊕ 3 2 3 ⊕ 3 2 1 ⊕ 3 2 3 1 ∅ Q(D • ) = 1 2 3
Note that this statement could be seen as a consequence of a more general result [13, Proposition 2.44]. However, the latter combines several isomorphisms (from dissections, via non-kissing paths, through non-attracting chords, to support τ -tilting modules) and the proof is quite intricate. In the present paper, we will easily deduce Theorem 1.2 from Theorem 1.1 and from the With this algebraic result, we can explain the isomorphism of Fig. 2, and extend it to dissections of any polygon. Any reference dissection of a polygon gives rise to an exchange graph on certain dissections. This exchange graph is the dual graph of the accordion complex studied in [START_REF] Chapoton | Stokes posets and serpent nests[END_REF][START_REF] Garver | Oriented flip graphs and noncrossing tree partitions[END_REF][START_REF] Manneville | Geometric realizations of the accordion complex of a dissection[END_REF] (see Section 4). Theorem 1.2. Any accordion complex is isomorphic to the support τ -tilting simplicial complex of an explicit finite dimensional algebra. Thus, the corresponding exchange graphs are isomorphic.

Note that this statement could be seen as a consequence of a more general result [13, Proposition 2.44]. However, the latter combines several isomorphisms (from dissections, via non-kissing paths, through non-attracting chords, to support τ -tilting modules) and the proof is quite intricate. In the present paper, we will easily deduce Theorem 1.2 from Theorem 1.1 and from the known case of triangulations of polygons (see Section 4).

The explicit finite dimensional algebras that appear in Theorem 1.2 had previously appeared in [START_REF] David-Roesler | Derived equivalence of surface algebras in genus 0 via graded equivalence[END_REF][START_REF] David-Roesler | Algebras from surfaces without punctures[END_REF], where gentle algebras are associated to dissections of any surface without puncture. Applying Theorem 1.1 to these algebras, one could obtain an analogue of the accordion complex for these surfaces.

Recollections on τ -tilting theory

The theory of τ -tilting modules was introduced in [START_REF] Adachi | τ -tilting theory[END_REF], and we mainly follow this source. Let k be an algebraically closed field, let Λ be a basic finite-dimensional k-algebra, and let {e 1 , . . . , e n } be a complete set of pairwise orthogonal idempotents in Λ. Denote by mod Λ the category of finite-dimensional right Λ-modules, and by proj Λ its full subcategory of projective modules. We denote by τ the Auslander-Reiten translation of mod Λ (see, for instance, [2, Chapter IV]). For any Λ-module M , we denote by |M | the number of pairwise non-isomorphic direct summands appearing in any decomposition of M into indecomposable modules. We will say that the pair (M, P ) is basic if both M and P are basic Λ-modules. We define direct sums of pairs componentwise.

Support τ -tilting pairs

One of the main theorems of [START_REF] Adachi | τ -tilting theory[END_REF] is the following. Definition 2.2. The support τ -tilting complex of Λ is the simplicial complex sτ C(Λ) whose vertices are the isomorphism classes of indecomposable τ -rigid pairs and whose faces are sets of τ -rigid pairs whose direct sum is rigid. The exchange graph sτ -tilt(Λ) is the dual graph of sτ C(Λ), i.e., the graph whose vertices are isomorphism classes of basic support τ -tilting pairs, and where two vertices are joined by an edge whenever the corresponding support τ -tilting pairs differ by exactly one direct summand.

2-term silting objects

The study of support τ -tilting pairs turns out to be equivalent to that of another class of objects: 2-term silting objects [1, Section 3]. Let K b (proj Λ) be the homotopy category of bounded complexes of projective Λ-modules. Let 2-cpx(Λ) be the full subcategory of K b (proj Λ) consisting of 2-term objects, that is, complexes

P = • • • → P m+1 → P m → P m-1 → • • •
such that P m is zero unless m ∈ {0, 1}. We will write P 1 → P 0 to denote the complex

• • • → 0 → P 1 → P 0 → 0 → • • • A 2-term object P is rigid if Hom K b (P, P[1]) = 0. It is silting if
it is rigid, and

|P| = |Λ|.
This is a special case of a more general definition of silting objects, see [START_REF] Keller | Aisles in derived categories[END_REF]. Examples of 2-term silting objects include 0 → Λ and Λ → 0.

Definition 2.3. The 2-term silting complex of Λ is the simplicial complex SC(Λ) whose vertices are isomorphism classes of indecomposable rigid 2-term objects in K b (proj Λ) and whose faces are sets of such objects whose direct sum is rigid. The exchange graph 2-silt(Λ) is the dual graph of SC(Λ), i.e., the graph whose vertices are isomorphism classes of basic 2-term silting objects in K b (proj Λ), and where two vertices are joined by an edge whenever the corresponding objects differ by exactly one direct summand.

For any Λ-module M , denote by P M 1 → P M 0 a minimal projective presentation of M .

The g-vector of a 2-term object

The results of this note rely on the following definition. Definition 2.5. Let P be a 2-term object in 2-cpx(Λ). The g-vector of P, denoted by g(P), is the class of P in the Grothendieck group K 0 K b (proj Λ) .

We will usually denote g-vectors as integer vectors by using the basis of the abelian group K 0 K b (proj Λ) given by the classes of the indecomposable projective modules e 1 Λ, . . . , e n Λ concentrated in degree 0. Thus, if P is the 2-term object

i∈[n] (e i Λ) ⊕b i - → i∈[n]
(e i Λ) ⊕a i , then its g-vector is g

(P) = (a i -b i ) i∈[n] .
In contrast to arbitrary 2-term objects, rigid 2-term objects are determined by their g-vector in the following sense. (i) If g(P) = g(Q), then P and Q are isomorphic.

(ii) The object P is isomorphic to an object of the form (P

1 → P 0 ) ⊕ Q id Q
→ Q , where P 1 and P 0 do not have non-zero direct summands in common.

Note that

Q id Q → Q is isomorphic to zero in K b (proj Λ).

Algebraic result

We use the same notations as in the previous section. In particular, Λ is a basic finite-dimensional k-algebra with complete set of pairwise orthogonal idempotents {e 1 , . . . , e n }.

Let J be a subset of [n]. We will study 2-term objects that only involve the indecomposable projective modules e j Λ with j ∈ J. Definition 3.1. Let 2-cpx J (Λ) be the full subcategory of 2-cpx(Λ) whose objects are the 2term objects P 1 → P 0 such that all the indecomposable direct summands of P 1 and P 0 have the form e j Λ with j ∈ J.

Our main interest will lie in the rigid objects in 2-cpx J (Λ). Definition 3.2. Let SC J (Λ) be the subcomplex of SC(Λ) induced by J, that is, the subcomplex whose vertices are rigid objects in 2-cpx J (Λ). Let 2-silt J (Λ) be the dual graph of SC J (Λ). Its vertices are isomorphism classes of basic objects P in 2-cpx J (Λ) satisfying P is rigid; if P ∈ 2-cpx J (Λ) and P ⊕ P is rigid, then P is a direct sum of direct summands of P.

Two vertices are joined by an edge whenever the corresponding objects differ by exactly one indecomposable direct summand.

In other words, the faces of SC J (Λ) correspond to basic rigid objects whose g-vectors have zero coefficients in entries corresponding to elements not in J. In this sense, SC J (Λ) is a representation-theoretic analogue of the accordion complex [START_REF] Chapoton | Stokes posets and serpent nests[END_REF][START_REF] Garver | Oriented flip graphs and noncrossing tree partitions[END_REF][START_REF] Manneville | Geometric realizations of the accordion complex of a dissection[END_REF] (see Theorem 4.3). This is the main motivation for the introduction of this object.

Let e J := j∈J e j , and consider the algebra e J Λe J . Observe that e J Λe J is isomorphic to End Λ (Λe J ). This has the following consequence. Let proj J (Λ) be the full subcategory of proj(Λ) whose objects are isomorphic to direct sums of the indecomposable modules e j Λ, with j ∈ J.

Lemma 3.3. The k-linear categories proj J (Λ) and proj(e J Λe J ) are equivalent. In particular, the categories K b (proj J (Λ)) and K b (proj(e J Λe J )) are equivalent. This lemma immediately implies the following statement.

Theorem 3.4. The simplicial complexes SC J (Λ) and SC(e J Λe J ) are isomorphic. In particular, their dual graphs 2-silt J (Λ) and 2-silt(e J Λe J ) are isomorphic.

Corollary 3.5. The simplicial complex SC J (Λ) is a pseudomanifold of dimension |J| -1. In particular, its dual graph 2-silt J (Λ) is |J|-regular.

Application: accordion complexes of dissections

Let P be a convex polygon. We call diagonals of P the segments connecting two non-consecutive vertices of P. A dissection of P is a set D of non-crossing diagonals. It dissects the polygon into cells. We denote by Q(D) the quiver with relations whose vertices are the diagonals of D, whose arrows connect any two counterclockwise consecutive edges of a cell of D, and whose relations are given by triples of counterclockwise consecutive edges of a cell of D. See Fig. 3 for an example.

We now consider 2m points on the unit circle alternately colored black and white, and let P • (resp. P • ) denote the convex hull of the white (resp. black) points. We fix an arbitrary reference dissection D • of P • . A diagonal δ • of P • is a D • -accordion diagonal if it crosses either none or two consecutive edges of any cell of D • . In other words, the diagonals of D • crossed by δ • together with the two boundary edges of P • crossed by δ • form an accordion. A D • -accordion dissection is a set of non-crossing D • -accordion diagonals. See Fig. 3 for an example. We call D • -accordion complex the simplicial complex AC(D • ) of D • -accordion dissections. This complex was studied in recent works of F. Chapoton [START_REF] Chapoton | Stokes posets and serpent nests[END_REF], A. Garver and T. McConville [START_REF] Garver | Oriented flip graphs and noncrossing tree partitions[END_REF], and T. Manneville and V. Pilaud [START_REF] Manneville | Geometric realizations of the accordion complex of a dissection[END_REF].

(D • | δ • ) ∈ R D• whose δ • -coordinate is ε δ • ∈ D • | δ • .
Example 4.1. When the reference dissection D • is a triangulation of P • , any diagonal of P • is a D • -accordion diagonal. The D • -accordion complex is thus an n-dimensional associahedron (of type A), where n = m -3. As explained in [START_REF] Caldero | Quivers with relations arising from clusters (An case)[END_REF], the D • -accordion complex is isomorphic to the 2-term silting complex of the quiver Q(D • ) of the triangulation D • . The isomorphism sends a diagonal of P • to the 2-term silting object with the same g-vector. See Fig. 1 for an illustration.

With the notations we introduced, we can now restate Theorem 1.2 more precisely. Theorem 4.2. For any reference dissection D • , the D • -accordion complex is isomorphic to the 2-term silting complex of the quiver Q(D • ).

One possible approach to Theorem 4.2 would be to provide an explicit bijective map between D • -accordion diagonals and 2-term silting objects for Q(D • ). Such a map is easy to guess using g-vectors, but the proof that it is actually a bijection and that it preserves compatibility is intricate. This approach was developed in the more general context of non-kissing complexes of gentle quivers with relations in [START_REF] Palu | Non-kissing complexes and tau-tilting for gentle algebras[END_REF]Propostion 2.44]. In this note, we use an alternative simpler strategy to obtain Theorem 4.2 by using Theorem 1.1 and understanding accordion complexes as certain subcomplexes of the associahedron.

For that, consider two nested dissections D In order to prove Theorem 4.2 we now turn to associative algebras. Let Q = (Q, I) be any gentle quiver with relations [START_REF] Butler | Auslander-Reiten sequences with few middle terms and applications to string algebras[END_REF] and J be any subset of vertices of Q. We call shortcut quiver the quiver with relations Q J = (Q J , I J ) whose vertices are the elements of J, whose arrows are the paths in Q with endpoints in J but no internal vertex in J, and whose relations are inherited from those of Q. Then the quotient kQ J /I J of the path algebra kQ J is gentle and is isomorphic to the algebra e J (kQ/I)e J .

Figure 2 .

 2 Figure 2. The D • -accordion complex of the dissection D • of Fig. 3 (left) and the 2-term silting complex of the quiver Q(D • ) (right).

Figure 2 .

 2 Figure 2. The D • -accordion complex of the dissection D • of Fig. 3 (left) and the 2-term silting complex of the quiver Q(D • ) (right).

Figure 2 .

 2 Figure 2. The D • -accordion complex of the dissection D • of Fig. 3 (left) and the 2-term silting complex of the quiver Q(D • ) (right).

Following [ 1 ,

 1 Definition 0.1], we say that a Λ-module M is τ -rigid if Hom Λ (M, τ M ) = 0; τ -tilting if it is τ -rigid and |M | = |Λ|; support τ -tilting if there exists an idempotent e of Λ such that e is in the annihilator of M and M is a τ -tilting Λ/(e)-module. Support τ -tilting modules always exist: Λ itself and the zero module are two examples. It is useful to keep track of the idempotents in the annihilator of a support τ -tilting module. For this reason, we will follow [1, Definition 0.3] and call a pair (M, P ), with M ∈ mod Λ and P ∈ proj Λ, a τ -rigid pair if M is τ -rigid and Hom Λ (P, M ) = 0; support τ -tilting pair if it is a τ -rigid pair and |M | + |P | = |Λ|; almost complete support τ -tilting pair if it is a τ -rigid pair and |M | + |P | = |Λ| -1.

Theorem 2 . 1 ([ 1 ,

 211 Theorem 0.4]). A basic almost complete support τ -tilting pair is a direct summand of exactly two basic support τ -tilting pairs.

Theorem 2 . 4 ([ 1 ,

 241 Theorem 3.2]). The map (M, P ) → P M 1 → P M 0 ⊕ (P → 0) induces an isomorphism of simplicial complexes sτ C(Λ) ∼ = SC(Λ), and thus of exchange graphs sτ -tilt(Λ) ∼ = 2-silt(Λ).

Theorem 2 . 6 ([ 9 ,

 269 Sections 2.3 and 2.4]). Let P and Q be two rigid 2-term objects.

Figure 3 .

 3 Figure 3. A dissection D • with its quiver Q(D • ) (left), a D • -accordion diagonal (middle) and a D •accordion dissection (right).

  • ⊂ D • . Observe that any D • -accordion diagonal is a D • -accordion diagonal. Conversely a D • -accordion diagonal δ • is a D • -accordion diagonal if and only if it does not cross any diagonal δ • of D • D • as a Z or a Z, that is if and only if the δ • -coordinate of its g-vector g(D • | δ • ) vanishes for any δ • ∈ D • D • .This observation shows the following statement.

Theorem 4 . 3 (

 43 [START_REF] Manneville | Geometric realizations of the accordion complex of a dissection[END_REF] Section 4.2]). For any two nested dissectionsD • ⊂ D • , the accordion complex AC(D • ) is isomorphic to the subcomplex of AC(D • ) induced by D • -accordion diagonals δ • whose g-vectors g(D • | δ • )lie in the coordinate subspace spanned by elements in D • .

Acknowledgements

We are grateful to F. Chapoton for pointing out to us the isomorphism between the two graphs of Fig. 2 which gave us the motivation for the present note. We also thank R. Schiffler for his comments on a previous version. Finally, we are grateful to an anonymous referee for helpful suggestions on the presentation of this note. The first two authors are partially supported by the French ANR grant SC3A (15 CE40 0004 01). The last author is supported by the ISF grant 1144/16.

In particular, for any dissection D • , the quiver Q(D • ) is a shortcut quiver of the quiver with relations of a cluster tilted algebra of type A.

The following statement is an application of Theorem 1.1 to gentle algebras. Theorem 4.5. For any gentle quiver with relations Q and any subset J of vertices of Q, the 2-term silting complex SC(Q J ) for the shortcut quiver Q J is isomorphic to the subcomplex of the 2-term silting complex SC(Q) induced by 2-term silting objects whose g-vectors lie in the coordinate subspace spanned by vertices in J. [START_REF] Manneville | Geometric realizations of the accordion complex of a dissection[END_REF]. The crucial feature of this fan is that no coordinate hyperplane meets the interior of any of its maximal cones. This is often referred to as the sign-coherence property of g-vectors. It implies that for any two nested dissections D

Combining

). The content of Theorem 4.3 is that this subfan is the g-vector fan F g (D • ). A similar statement holds for Theorem 4.5.

Remark 5.2. In the theory of cluster algebras, a standard operation consists of freezing a subset of the initial cluster. This corresponds to taking a section of the d-vector fan by a coordinate subspace. To the best of our knowledge, the same operation on the g-vector fan studied in this note was not considered before in the literature.

Remark 5.3. The connection between representation theory and accordion complexes was already considered by A. Garver and T. McConville in [START_REF] Garver | Oriented flip graphs and noncrossing tree partitions[END_REF]Section 8]. However, their approach deals with c-vectors and simple-minded collections while our approach deals with g-vectors and silting objects.


Andrea Petreto

Thomas Romera

Ian Masliah

Boris Gaillard
email: bgaillard@lheritier-alcen.com

Manuel Bouyer

Quentin L Meunier

Lionel Lacassagne

Florian Lemaitre

A New Real-Time Embedded Video Denoising Algorithm

Keywords: Embedded system, Video denoising, Real-time, Image processing, night-vision, SIMD

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

I. INTRODUCTION

Image denoising is a major research area for computer vision. It can be based on methods such as pixel-wise filters [START_REF] Zhang | Multiresolution bilateral filtering for image denoising[END_REF]- [START_REF] Buades | A review of image denoising algorithms, with a new one[END_REF], patch-based methods [START_REF]A non-local algorithm for image denoising[END_REF], [START_REF] Lebrun | A nonlocal bayesian image denoising algorithm[END_REF] or Convolutional Neural Networks (CNN) [START_REF] Burger | Image denoising: Can plain neural networks compete with bm3d?[END_REF]- [START_REF] Zhang | Beyond a gaussian denoiser: Residual learning of deep cnn for image denoising[END_REF]. Video denoising can take advantage of the temporal dimension to reduce the noise more efficiently.

Denoising algorithms can be classified according to two criteria: the processing time and the denoising efficiency. In this work, we focus on video with heavy noise and we target real-time processing at 25 frames per second (fps) on embedded systems. Concerning the image quality, we focus on details enhancement rather than aesthetic rendering. Some video denoising algorithms are able to deal with heavy noise [START_REF] Zuo | Video denoising based on a spatiotemporal kalman-bilateral mixture model[END_REF]- [START_REF] Arias | Towards a bayesian video denoising method[END_REF]. However, all these algorithms are slow and cannot be used for real-time denoising. Moreover, they have a high latency as they need future frames to process the current one (typically 5 to 7).

There are algorithms designed for real-time denoising [START_REF] Pfleger | Real-time video denoising on multicores and gpus with kalman-based and bilateral filters fusion[END_REF] on embedded architectures [START_REF] Ehmann | Real-time video denoising on mobile phones[END_REF], but only for light noise like compression artifacts, which give poor results on heavy noise.

As far as we know, there is no real-time algorithm able to handle heavily noisy video as in this work.

In this paper, we introduce a new Real-Time Embedded Video Denoising (RTE-VD) algorithm targeting very noisy video in section II, which has been implemented for general purpose processors. We then compare our algorithm to the state-of-the-art in terms of denoising efficiency in section III. We then present the main optimizations applied to achieve real-time performance in section IV. Finally in section V, we compare the execution time of RTE-VD to other state-of-theart algorithms and study its speed and power consumption depending on various embedded CPUs and frequencies.

II. DENOISING ALGORITHM

Most state-of-the-art video denoising algorithms rely on patch based filtering methods [START_REF] Arias | A comparison of patch-based models in video denoising[END_REF], [START_REF] Maggioni | Video denoising using separable 4d nonlocal spatiotemporal transforms[END_REF], [START_REF] Arias | Towards a bayesian video denoising method[END_REF] which tend to be effective to reduce noise but are very time consuming. In addition, their memory access patterns generate many cache misses, making it difficult to speed the computation up. Therefore, we did not consider this technique for real-time denoising.

A crucial aspect of video denoising is to keep temporal coherence between frames. While we cannot ensure this coherence using patch search, it is possible to do so with an optical flow estimation [START_REF] Buades | Video denoising with optical flow estimation[END_REF], [START_REF] Liu | A high-quality video denoising algorithm based on reliable motion estimation[END_REF].

Optical flow algorithms are also time consuming but are sensitive to code transformations and can be highly accelerated [START_REF] Petreto | Energy and execution time comparison of optical flow algorithms on SIMD and GPU architectures[END_REF]- [START_REF] Bako | Embedded implementation of a real-time motion estimation method in video sequences[END_REF]. We thus decided to use an optical flow algorithm in our real-time denoising chain.

While there exists many optical flow algorithms [START_REF] Middlebury | Optical flow database[END_REF], optical flow estimation for video denoising should be dense, robust to noise and handle flow discontinuities [START_REF] Liu | A high-quality video denoising algorithm based on reliable motion estimation[END_REF]. Therefore, we decided to use the TV-L1 algorithm to compute the optical flow, since it satisfies all these constraints [START_REF] Zach | A duality based approach for realtime tv-l 1 optical flow[END_REF]. Our complete denoising application is composed of three steps (Figure 1). First, the video stream is stabilized using a global one-pass Lucas-Kanade approach [START_REF] Lucas | An iterative image registration technique with an application to stereo vision[END_REF] in order to compensate camera motion. Then, the dense optical flow is computed using the TV-L1 algorithm. The flow enables to extrapolate denoised past frames at the current frame, and match the two together. Finally, a 3D spatio-temporal bilateral filter is applied [START_REF] Tomasi | Bilateral filtering for gray and color images[END_REF].

Fig. 2. Visual comparison on pedestrians sequence with a standard deviation noise of 40 (PSNR in Table I). In order to evaluate the denoising efficiency of RTE-VD, we reproduced the conditions of the comparison made in [START_REF] Davy | Non-local video denoising by cnn[END_REF]. We ran our algorithm on seven sequences from the Derf's Test Media Collection [START_REF] Derf | Derf's test media collection[END_REF]. The sequences are originally colored 1080p video. We down-scaled them to 960×540 pixels in grayscale by keeping only brightness. A Gaussian noise is added to each pixel with a standard deviation of 20 and 40. Then, the algorithms are ran on the first hundred frames of each sequence.

However, we were not able to reproduce the same results as in [START_REF] Davy | Non-local video denoising by cnn[END_REF]: all the PSNR values that we have obtained are significantly lower for a given method.

We compare RTE-VD to VBM3D [START_REF] Dabov | Video denoising by sparse 3d transform-domain collaborative filtering[END_REF], VBM4D [START_REF] Maggioni | Video denoising using separable 4d nonlocal spatiotemporal transforms[END_REF] and STMFK [START_REF] Pfleger | Real-time video denoising on multicores and gpus with kalman-based and bilateral filters fusion[END_REF]. VBM3D is one of the most popular video denoising method. It is also one of the fastest high performance denoising algorithm even if its computation time is still over a second for a 960×540 pixel frame on a Xeon server. VBM4D has a more powerful filtering than VBM3D but at a higher computational cost. STMKF is a state-of-the-art realtime video denoising method. For each algorithm, the source code provided on the authors website is used with the default parameters depending on the noise intensity.

The quantitative results for the test sequences are presented in Table I. Our real-time constraint limits the denoising efficiency of RTE-VD compared to VBM3D and VBM4D; however, RTE-VD outperforms STMKF. RTE-VD and STMKF are able to perform real-time denoising while VBM3D and VMB4D require offline computation. Processing times are studied in more details further in this article.

Some visual results are presented in Figure 2: we can see that RTE-VD is a clear improvement compared to the noisy input. In most cases, VBM3D and VBM4D are visually better. Also STMKF is much less effective than the other considered methods.

For the two noise deviation tested, the RTE-VD PSNR is overall less than 3 dB below the ones from VBM3D and VBM4D. RTE-VD PSNR is overall 2 dB above STMKF for a noise deviation of 20 and more than 4 dB above for a noise deviation of 40. It confirms that RTE-VD handles heavy noisy situations better than STMKF, as can be seen in Figure 2.

RTE-VD tends to trade noise reduction off for information keeping. Figure 2 shows that RTE-VD is able to produce more detailed results than STMKF on moving objects but encounters difficulties on still objects. The loss of sharpness in static scenes can be explained by a bad optical flow estimation on non-moving objects. Error estimations are indeed proportionally more important on small movements and bad compensations may generate a little blur effect on non moving objects. Filtering the estimated optical flow may reduce this effect and will be considered in our futur works.

IV. ALGORITHM OPTIMIZATIONS

In order to run RTE-VD in real-time on an embedded system, several optimizations and tradeoffs have to be done.

We apply the following high-level transformations (described in [START_REF] Lacassagne | High level transforms for SIMD and low-level computer vision algorithms[END_REF]) to each part of the chain:

• SIMDization (handcrafted Neon instructions as the compiler vectorization is inefficient), • Multi-thread parallelization (with OpenMP),

• Operator fusion to reduce memory accesses, • Operator pipeline to increase memory access locality, • Modular memory allocation (for cache blocking) to reduce the memory footprint and enforce memory access locality.

A. Global Lucas-Kanade stabilization

To stabilize the video stream, we use a modified version of the Lucas-Kanade method [START_REF] Lucas | An iterative image registration technique with an application to stereo vision[END_REF] with a global approach. the input image is convolved with a gate function, (a function that is 0 outside a specified interval and 1 inside it) in order to compute its average over a neighborhood. As the gate size depends on the largest movement to compensate, the convolution kernels can be very large. To speed the convolution up, integral images (also known as summed area table [START_REF] Crow | Summed-area tables for texture mapping[END_REF]) are used. The convolution is computed by multiple threads, each one processing a strip. Consequently, each thread only computes the partial integral image it needs. We then apply the Lucas-Kanade flow estimation to a neighborhood of the same size as the input images.

B. TV-L1 Dense Optical Flow Estimation

The TV-L1 optical flow algorithm is the main step to optimize as it takes more than 80% of the total time. The key transformations of TV-L1 are presented in [START_REF] Petreto | Energy and execution time comparison of optical flow algorithms on SIMD and GPU architectures[END_REF] and result in a processing 5× faster and 6× less power consuming on an embedded ARM Cortex A57 architecture.

With n scales and a zoom factor of 1: 2, the largest movement between two images that the algorithm can handle is equal to 2 n -1. Therefore we fixed the scale number to 3 to be able to estimate displacements as large as 7 pixels.

Instead of using a convergence criteria to stop the iterations, a fixed number of iterations per scale has been chosen. In most situations, three iterations per scale are enough. However, more complex situations, like large displacements or discontinuities, require more iterations. Therefore, we set 10 iterations per scale as the standard configuration. This provides a more robust estimation than with only three iterations per scale, while keeping a rather low computation time.

However, as the most important movements to compensate are the largest ones, it is interesting to do more iterations on the most zoomed-out scales (the smallest ones). Therefore, we compare the computation of 10 iterations per scale on 3 scales (10-10-10) with the computation of 3 iterations at the largest scale, 20 iterations at the medium scale and 80 iterations at the smallest one (3-20-80). While the two have roughly the same number of operations, the (3-20-80) version is actually faster as it spends more time on the smallest scale which is cache friendlier. The denoising quality is similar.

Another critical step of TV-L1 to ensure stability is the warping, which is in our case a motion compensation using bicubic interpolation. Several warps can be applied at each scale. They take a long time to compute, but gives better results. This leads to important design choices, in particular regarding the number of warps at each scale. Thus, we compared the (10-10-10) single warp per scale configuration with a configuration comprising 3 iterations with 1 warp at the largest scale, 20 iterations with 2 warps at the medium scale and 80 iterations with 4 warps at the smallest scale. This latter configuration is the one used in the remaining experiments, as it significantly achieves better results with a similar computation time for 960×540 pixel images.

C. Spatio-temporal filter

To actually filter out the noise, we compose a spatial bilateral filter [START_REF] Tomasi | Bilateral filtering for gray and color images[END_REF] with a unilateral filter to handle the temporal dimension and speed the processing up. Thus, we are able to decorrelate the strength of the filter in the spatial domain from the temporal domain. The filter is defined by the equations 1, 2 and 3, where I f is the filtered image, I p the previous compensated and filtered image, I the current image and x the coordinates of the current pixel. Ω is the filter kernel domain. σ i , σ d and σ t are the smoothing parameters of the filter respectively of the spatial intensity difference, the distance and the temporal intensity difference between pixels.

I f (x) = 1 W p xi∈Ω I t (x i )e -[I t (x i )-I t (x)] 2 2 σ 2 i × e -[x i -x] 2 2 σ 2 d (1)
Where:

W p = xi∈Ω e -[I t (x i )-I t (x)] 2 2 σ 2 i × e -[x i -x] 2 2 σ 2 d (2)
and:

I t (x) = I p (x)e -[Ip(x)-I(x)] 2 2 σ 2 t + I(x) 1 -e -[Ip(x)-I(x)] 2 2 σ 2 t
(3) While the bilateral filter is not separable, we approximate it by a separable filter [START_REF] Pham | Separable bilateral filtering for fast video preprocessing[END_REF] which is faster than the exact one. We also used a fast approximation of the exponential function which manipulates the bit representation of floats as defined in the IEEE-754 standard. Such an approximation is described in [START_REF] Schraudolph | A fast, compact approximation of the exponential function[END_REF] and is accurate enough.

V. PROCESSING TIME AND ENERGY CONSUMPTION

In this section we first compare the execution time of RTE-VD to the algorithms previously introduced. We then briefly study the impact of our optimizations on computation speed. Finally, we evaluate the performances of RTE-VD on various embedded CPUs and frequencies.

We consider 4 different platforms. The first one is an Intel Xeon Silver 4114 2×10C/20T@2.20GHz. The 3 others are the latest Nvidia Jetson embedded platforms; for those platforms, we only consider the ARM CPUs and not the GPUs. Their names and specifications are given in Table II. 2.0 Xeon VBM4D [START_REF] Maggioni | Video denoising using separable 4d nonlocal spatiotemporal transforms[END_REF] 45 Xeon STMKF [START_REF] Pfleger | Real-time video denoising on multicores and gpus with kalman-based and bilateral filters fusion[END_REF] 0.015 AGX RTE-VD (this work) 0.037 AGX

A. Processing time analysis

We measured the computation time for the different methods considered in section III. We ran RTE-VD, STMKF, VBM3D and VBM4D on the Xeon platform. We also tested RTE-VD and STMKF on the AGX platform. The results for 960×540 pixel images are presented in Table III. On the same platform, RTE-VD is more than 200× faster than VBM3D and more than 4600× faster than VBM4D. On the AGX platform, RTE-VD is only 2.5× slower than STMKF but still achieves realtime processing with 26.7 frames per second.

To exhibit the efficiency of our optimizations we also compare our Fast optimized implementation to a Slow straightforward implementation. The integral image is hard to parallelize without using major transformations. Those transformations have only been applied to the Fast version. Thus, in order to be fair, we compare the Fast and Slow versions in mono-thread and then analyze only the Fast version in multi-thread. The images used for the experiments are square images, with a width varying from 200 to 1500 pixels. The size of the image has negligible impact on the speed of the processing chain. Thus, 1000×1000 images are considered for the following experiments.

The results for all algorithms involved in RTE-VD on the AGX are presented in Table IV. It shows that in mono-thread, we have an overall speedup of ×18, and in multi-thread, a speedup of ×70 using the 8 cores of the AGX. We can also observe that the major part of the speedup is obtained on the filtering, mainly due to its approximation with a separable filter. As a consequence, the optical flow estimation now takes almost 98% of the total time of the Fast version, while filtering was the most time consuming part of the processing chain without the optimizations. Since the optical flow estimation is, in the end, the most consuming step, we produced a more detailed study in [START_REF] Petreto | Energy and execution time comparison of optical flow algorithms on SIMD and GPU architectures[END_REF]. In this previous work the impact of each optimisation is discussed. Both time and power consumption are considered.

B. Time and energy efficiency of RTE-VD

Since we target embedded systems, we have to consider not only the computation speed but also the energy consumption. Therefore, we have run RTE-VD at various frequencies on the three Nvidia embedded systems. The frequencies have been taken among the available frequencies of each board and the external memory frequency has been set to its maximum. For the AGX and NANO, we have used a multi-threaded version on all the physical cores, respectively 8 and 4. For the TX2, we have first used only the two Denver cores, then only the four A57 cores and finally all 6 cores. The cooling system of each target board has been set to the maximum and the energy saving policies of the Operating System have been deactivated. We have simultaneously measured time and power consumption for various frequencies and image sizes.

In order to perform simple and reproducible power measurements, the electrical consumption of the entire system has been measured. A board was developed to this effect and has been inserted between the power source and the target system. The board samples both voltage and current at 5 kHz. Measurements and code executions are synchronized using GPIOs.

We can define four different metrics related to power consumption:

• Static energy: the energy associated to the static power when the system is idle (here, when the power consumption comes from the leakage and running the operating system) • Dynamic energy: the energy associated to the dynamic power: that is the extra energy consumed by the computation: the total energy minus the static energy. • Compute energy: the energy consumed by the whole system, that is the sum of static and dynamic energy. • Period energy: the energy consumed by the whole system between two execution starts, including the waiting time. The period energy is more relevant if we consider a complete system for which we know all of its applications and behaviours. Since this is not yet the case here, we only consider compute and dynamic energy in the following.

Our measurements show that considering the compute energy, the maximum frequency is always the most efficient for all configurations, being both the fastest and the least Figure 3 represents the computation time in nanosecond per pixel and the dynamic energy consumption in nanojoule per pixel for all the boards and frequencies. The NANO is the least energy consuming; however, it is possible to be 2.7× faster while consuming almost as little energy using the AGX. The AGX is the fastest since it is 2.3× faster than the fastest configuration of the TX2 and 4.4× faster than the fastest NANO configuration. Due to its etching process (16nm), the TX2 is often less energy efficient and slower than the two others (12nm). This is especially visible if we compare the NANO to the TX2 using only its A57 quad core at equivalent frequency. Since the NANO also possesses a quad core A57, it is able to compute as fast as the TX2 while consuming 1.6× less. However the maximum frequency of the TX2 being higher, it is able to be faster than the NANO even using only the A57 cores.

The Table V synthesizes the best configurations minimizing either energy consumption or computation time for each platform. It gives for each configuration the largest image size possibly processed at 25 frames per second.

VI. CONCLUSION

In this article, we proposed a novel real-time embedded video denoising chain called RTE-VD, which is able to restore details on very noisy video while achieving high performance on embedded systems.

In order to achieve real-time processing, we applied several code transformations (like SIMDization, multi-threading, operator fusion and pipelining) and were able to get an implementation 70× faster than a naive one. The PSNR and visual results from our experiments validates our approach.

We thus compared RTE-VD to other state-of-the-art algorithms: VBM3D, VBM4D and STMKF. RTE-VD always denoises better than STMKF while being less effective than the more costly algorithms VBM3D and VBM4D. On heavy noise situations (σ = 40), RTE-VD achieves an overall PSNR more than 4 dB above STMKF. RTE-VD is also 200× faster than VBM3D and 4600× faster than VBM4D.

While RTE-VD is 2.5× slower than STMKF, it is still able to process 960×540 pixel video at 25 fps on a Nvidia Tegra AGX. Given these results, we believe that RTE-VD is a denoising algorithm particularly well positioned for speed/accuracy tradeoff.

Since we are targeting embedded systems, we also studied the link between time and energy consumption on various embedded CPUs. Within the tested platforms, it appears that the Nvidia Jetson Nano is the least power consuming platform while the Nvidia Tegra AGX is the fastest. We were also able to determine for each platform multiple efficient frequencies minimizing either computation time or energy consumption.

As future work, we plan to increase the whole performance by balancing the load on the CPU and the GPU and by studying what parts of the algorithm can be hybridized with 32-and 16-bit computation.

Fig. 1 .

 1 Fig. 1. Principle components of the denoising application.

Fig. 3 .

 3 Fig. 3. Speed and energy efficiency of RTE-VD for embedded CPUs at different frequencies

TABLE I PSNR

 I RESULTS FOR SEVEN SEQUENCES OF THE Derf's Test Media Collection ON STATE-OF-THE-ART ALGORITHMS AND OUR METHOD (RTE-VD).

	Noise	Method	crowd	park joy	pedestrians	station	sunflower	touchdown	tractor	overall
		STMKF [15]	26.25	25.59	28.34	26.66	26.97	28.87	25.37	26.70
	σ = 20	RTE-VD VBM3D [12]	26.38 28.75	25.65 27.89	30.58 35.49	30.98 34.19	32.51 35.48	30.17 32.85	29.38 31.44	28.73 31.34
		VBM4D [13]	28.43	27.11	35.91	35.00	35.97	32.73	31.65	31.11
		STMKF [15]	20.80	20.75	20.70	20.41	20.70	20.86	19.80	20.56
	σ = 40	RTE-VD VBM3D [12]	22.55 24.81	21.64 23.78	25.72 30.65	27.76 30.62	27.87 30.88	27.05 30.21	25.99 27.82	24.85 27.43
		VBM4D [13]	24.65	23.22	31.32	31.53	31.39	30.09	28.09	27.35
	III. DENOISING EFFICIENCY COMPARISON

TABLE II TECHNICAL

 II SPECIFICATIONS OF THE TARGET EMBEDDED BOARDS.

	Board Process	CPU	Fmax (GHz)	Idle Power (W)
	TX2	16 nm 4×A57 + 2×Denver2 2.00	2.0
	AGX	12 nm	8×Carmel	2.27	6.3
	NANO 12 nm	4×A57	1.43	1.2

TABLE IV EXECUTION

 IV TIME (MS) AND SPEEDUP OF RTE-VD ON AGX CPU

	Algorithm	Slow 1C	Fast 1C	Fast 8C	speedup
	LK	6.66	1.59	0.37	×18
	Flow	260.73	107.93	27.59	×10
	Filter	1 717.85	1.39	0.25	×6 871
	Total	1 985.24	110.90	28.21	×70

ACKNOWLEDGMENT

The authors would like to thank the French Direction Générale de l'Armement for supporting this work.


Arthur Hennequin
email: arthur.hennequin@lip6.fr

Ben Couturier

Vladimir V Gligorov

Lionel Lacassagne
email: lionel.lacassagne@lip6.fr

SparseCCL: Connected Components Labeling and Analysis for sparse images

Connected components labeling and analysis for dense images have been extensively studied on a wide range of architectures. Some applications, like particles detectors in High Energy Physics, need to analyse many small and sparse images at high throughput. Because they process all pixels of the image, classic algorithms for dense images are inefficient on sparse data. We address this inefficiency by introducing a new algorithm specifically designed for sparse images. We show that we can further improve this sparse algorithm by specializing it for the data input format, avoiding a decoding step and processing multiple pixels at once. A benchmark on Intel and AMD CPUs shows that the algorithm is from ×1.6 to ×2.5 faster on sparse images.

I. INTRODUCTION

In computer vision, Connected Component labeling (CCL) is a common and wide spread algorithm. It consists in assigning a unique label to each group of connected pixels. These groups of pixels, called Connected Components (CC), are then used for higher level tasks, like tracking, motion detection or optical character recognition. First instances of this algorithm were proposed by pioneers like Rosenfeld [START_REF] Rosenfeld | Sequential operator in digital pictures processing[END_REF] or Haralick [START_REF] Haralick | Some neighborhood operations[END_REF]. In High Energy Physics (HEP), CCL is used in the tracking of particles by labeling hits on the detectors' sensors to extract the real impact positions.

A CCL algorithm by itself, only provides the association of pixels, this is why it is followed by an analysis algorithm. The purpose of the analysis is to compute features of each CC, like the bounding box or the first statistical moments i order to to compute the center of gravity. If naive algorithm perform the labeling first and then the analysis, the optimized algorithms do the analysis during the labeling. These algorithms are called Connected Component Analysis (CCA).

Most of CCL algorithms used to be sequential ones developed on single-core processors [START_REF] He | The connectedcomponent labeling problem: a review of state-of-the-art algorithms[END_REF] [START_REF] Bolelli | Toward reliable experiments on the performance of connected components labeling algorithms[END_REF]. Recently, new parallel algorithms were developed for multi-core processors [START_REF] Niknam | A parallel algorithm for connected component labeling of gray-scale images on homogeneous multicore architectures[END_REF] [START_REF] Gupta | A new parallel algorithm for two-pass connected component labeling[END_REF], SIMD processors [START_REF] Wende | Swendsen-wang multi-cluster algorithm for the 2d/3d Ising Model on Xeon Phi and GPU[END_REF] [11] [START_REF] Hennequin | Designing efficient simd algorithms for direct connected component labeling[END_REF] and GPUs [START_REF] Playne | A new algorithm for parallel connectedcomponent labelling on GPUs[END_REF] [START_REF] Hennequin | A new direct connected component labeling and analysis algorithm for GPUs[END_REF]. These algorithms are very efficient for natural images but not for very low density images (very few pixels set to one) like those generated in HEP experiment.

The case considered here is similar to matrix algebra. When a matrix has very few non-zero value, the classical dense structure and the classical dense algorithms turn out to be inefficient. The dense structure is replaced by various flavours of lists that hold the non-zero value and specialized or dedicated algorithms are designed to process these data efficiently. In the case of tracking hits on detectors' sensors, the same phenomenon happens: even if CCA algorithms are very fast, they are inefficient to cluster and label matrices of hits with around 0.5% of hits.

The section II presents some classic connected components labeling algorithms. Section III introduces a new algorithm for connected components labeling of sparse images and its specialization for the pattern recognition of CERN's LHCb experiment. Finally, in section IV, we evaluate this new algorithm and compare it to state-of-the-art.

II. CLASSIC ALGORITHMS FOR DENSE IMAGES

In this section, we present three classes of connected components labeling algorithms.

A. One component at a time

In this first class of algorithm, we process one connected component at a time. The image is scanned one time and, for every foreground pixel encountered, a traversal of the connected component is done to label all the pixels. This algorithm and its variants are often called flood fill or sometimes seed fill. The traversal can be done using a stack in depth-first order or a queue, in breadth-first order. Implementations of algorithms of this class are found in [START_REF] Vincent | Watersheds in digital spaces: An efficient algorithm based on immersion simulations[END_REF] and [START_REF] Abubaker | One scan connected component labeling technique[END_REF]. This algorithm can be optimized by only adding, on top of the stack, the branching pixels -ie. the pixels that have more than one non-visited neighbour -and directly processing the others. Doing so, we avoid a store and a load for these pixels. If the image is sparse and if we have a list of pixel coordinates, we can directly start at known pixel positions avoiding the read of many background pixels. However, this does not prevent the test of every pixel on the contour of the connected component and the visited pixels have to be removed from the list.

B. Iterative algorithms

The second class was introduced by Haralick [START_REF] Haralick | Some neighborhood operations[END_REF]. Each pixel is initialized with a unique temporary label, then this label is propagated to the pixel's neighbors using local minimum or maximum propagation. The propagation step is repeated until the image of labels reaches stabilization, ie. there is no more change within the image. This algorithm was particularly fitted for implementations on parallel architectures, due to its high regularity, before the apparition of fast scatter and gather operations needed for union-find based algorithms. The number of iterations, and thus the processing time, of iterative algorithms depends on the longest path in the image. For an n × n spiral, the number of iterations is equal to n The two-pass CCL algorithms are split into three steps and perform two image scans (like the pioneer algorithm of Rosenfeld [START_REF] Rosenfeld | Sequential operator in digital pictures processing[END_REF]). The first scan (or first labeling) assigns a temporary label to each CC and some label equivalences are built if needed. The second step solves the equivalence table by computing the transitive closure of the graphs associated to the label equivalences. The third step performs a second scan (or second labeling) that replaces temporary labels of each CC by its final label, by doing a simple look up:

I(i, j) ← L[I(i, j)].
Figure 1 defines some notations and gives an example of a classic Rosenfeld algorithm execution. Let p x , e x , the current pixel and its label. Let p a , p b , p c , p d , the neighbor pixels, and a, b, c, d, the associated labels. L is the equivalence table, e a label and r its root. The first scan of Rosenfeld is described in algorithm 1, the transitive closure in algorithm 2, while the classical union-find algorithms are provided in algorithms 3 & 4. In the example (Fig. 1), we can see that the rightmost CC requires three labels (1, 2 and 4). When the mask is in the position seen in figure 1 (in bold type), the equivalence between 2 and 4 is detected and stored in the equivalence table L. At the end of first scan, the equivalence table is complete and applied to the image.

III. CONTEXT AND SPARSECCL

In this section we present SparseCCL, a parameterizable connected components labeling and analysis algorithm for sparse images. Then, we present a specialization of the algorithm in the context of the LHCb experiment, where the Algorithm 1: Rosenfeld algorithm -first labeling (step 1)

Input: a, b, c, d, four labels, px, the current pixel in (i, j)

1 if px = 0 then 2 a ← E[i -1][j -1], b ← E[i -1][j] ; 3 c ← E[i -1][j + 1], d ← E[i][j -1] ; 4 if (a = b = c = d = 0) then 5 ne ← ne + 1, ex ← ne ; 6 else 7 ra ← F ind(L, a), r b ← F ind(L, b) ; 8 rc ← F ind(L, c), r d ← F ind(L, d) ;
1 for e ∈ [1 : n] do 2 L[e] ← L[L[e]]
very few hits of high-energy particles are scattered across the detector's sensors.

A. General parameterizable ordered SparseCCL

In this first version of the algorithm, we assume that the image is represented by a list of active pixels ordered by their coordinates. This kind of representation allows the algorithm to take advantage of the sparse nature of the data. This is due to the size of the list scaling directly with the number of pixels to label and not the total number of pixels. Another case where this representation is useful is when the coordinate ranges are too large or if the number of dimensions makes the storage not practical. Figure 2 gives an example of such an image and its list representation.

The algorithm parameterization is done through the functions is_adjacent / is_far_enough for the labeling and init_features / accumulate_features for the analysis. The algorithm is generic enough to be adapted to n-dimension and every type of connectivity and pixel format. Algorithm 5 gives the complete algorithm. (1,0), [START_REF] Hennequin | Designing efficient simd algorithms for direct connected component labeling[END_REF][START_REF] Cabaret | What is the world's fastest connected component labeling algorithm[END_REF], [START_REF] Hennessy | LHCb VELO upgrade[END_REF][START_REF] Cabaret | What is the world's fastest connected component labeling algorithm[END_REF], [START_REF] Hennequin | A new direct connected component labeling and analysis algorithm for GPUs[END_REF][START_REF] Gupta | A new parallel algorithm for two-pass connected component labeling[END_REF], [START_REF] Hennessy | LHCb VELO upgrade[END_REF][START_REF] Gupta | A new parallel algorithm for two-pass connected component labeling[END_REF], [START_REF]LHCb VELO Upgrade[END_REF][START_REF]LHCb VELO Upgrade[END_REF], [START_REF] Bolelli | Toward reliable experiments on the performance of connected components labeling algorithms[END_REF][START_REF] Niknam | A parallel algorithm for connected component labeling of gray-scale images on homogeneous multicore architectures[END_REF], [START_REF] Lacassagne | A new SIMD iterative connected component labeling algorithm[END_REF][START_REF] Niknam | A parallel algorithm for connected component labeling of gray-scale images on homogeneous multicore architectures[END_REF], [START_REF] Cabaret | What is the world's fastest connected component labeling algorithm[END_REF][START_REF] Playne | A new algorithm for parallel connectedcomponent labelling on GPUs[END_REF] { } Algorithm 5: SparseCCL SparseCCL is designed to minimize the memory footprint in order to have the best data locality and to fit in the L1 cache. The only internal memory it needs is an integer table of size n to store the equivalences. The input table containing the pixels and the output tables containing the connected components features are allocated outside of the algorithm. The algorithm is divided into two parts: the first scan where pixels are associated using an equivalence table and the second scan where we resolve equivalences by performing a transitive closure of the graph embedded in the equivalence table. We can also do an on-the-fly analysis of connected components in the second scan.

1 // First scan: pixel association 2 start j ← 0 3 for i ← 0 to n -1 do 4 L[i] ← i 5 ai ← i 6 for j ← start j to i -1 do 7 if is_adjacent(pixel[i], pixel[j]) then
The first scan iterates over the pixels in the list and adds them one by one to the equivalence table. The equivalence table is an index table implementing a forest of equivalence trees. Each cell of the table corresponds to one pixel, the content of the cell is the index of the parent pixel. A pixel is a root if its entry in the equivalence table is its own index. For each pixel, the algorithm checks on previously added pixels for adjacency and merge their two equivalence trees if they are adjacent. The merging is done by calling the union and find functions described in algorithms 4 and 3. Because the list is ordered, we can keep track of a start index to avoid testing pixels that are too far from each other. With this optimization, the first scan complexity becomes O (kn) instead of O (n(n -1)/2), with k << n. The is_adjacent and is_far_enough parameterization for 2-dimensions 8connectivity labeling is described in algorithms 6 and 7.

Algorithm 6: is adjacent(p 1 , p 2 ) 1 return |p 1 .x -p 2 .x| ≤ 1 and |p 1 .y -p 2 .y| ≤ 1 Algorithm 7: is far enough(p 1 , p 2 ) 1 return p 1 .y -p 2 .y > 1
The second scan iterates over each temporary label in the equivalence table. If the label is a root, it creates a new connected component label by incrementing the labels counter and initialises the features for this connected component. If the label has a parent, it takes the parent's label and accumulates the features. Because the temporary label of a parent is always smaller than the one of the child, we know that the parent is already processed. Before continuing to the next label, we update the equivalence table with the new label. Algorithms 8 and 9 give an example of the features needed to compute the connected components center of gravity (G x , G y ) = (S x /S, S y /S). (S x , S y ) are the sums of x and y coordinates and S the number of pixels.

Algorithm 8: init features(label, pixel)

1 sumx[label] ← pixel.x 2 sumy[label] ← pixel.y 3 sumn[label] ← 1
Algorithm 9: accumulate features(label, pixel)

1 sumx[label] ← sumx[label] + pixel.x 2 sumy[label] ← sumy[label] + pixel.y 3 sumn[label] ← sumn[label] + 1

B. Acceleration structure for unordered pixels

In some scenarios, we might not receive an ordered list of pixels as input and sorting them would already take too much time. We also cannot afford accessing a full pixel image buffer because of data locality and the time it will take to reset such buffer between two images labeling. We compromise by doing dimension reduction: we use a table for each row and add the pixels to its corresponding row when we encounter it in the first scan. Now, when checking for adjacency, we only have to check the previous, the current and the next row table. The pixels within each row are not sorted so they have to be all checked. Each row table has a size property (N row ) that keeps track of how many pixels were added to the row. When resetting the tables, we only have to set the size of used rows to zero. Table I

C. Case study: specialization for LHCb VELO Upgrade

LHCb, one of the four major experiments at the Large Hadron Collider (LHC) [START_REF] Alves | The LHCb Detector at the LHC[END_REF], is a general purpose spectrometer optimized for the study of particles containing b and anti-b quarks (B mesons). The experiment's detector is specifically designed to filter out these particles and the products of their decay. Each LHCb sub-detector is specialized in measuring a different characteristic of the particles produced by colliding protons. Collectively, the detector's components gather information about the identity, trajectory, momentum and energy of each generated particle, and can single out individual particles from the billions that spray out from the collision point each second. After the upcoming high luminosity upgrade, the LHC will collide protons 30 million times per second at LHCb. These collision "events" must be analysed in real time by a system called the High Level Trigger 1 (HLT1), in order to find and keep the small fraction of these events which contain B mesons.

The VELO (VErtex LOcator) sub-detector, shown in figure 3, is a high precision pixel detector surrounding the beamline where the collision occurs. It is divided into 52 Lshaped modules. Each module is itself composed by 4 sensors of 3 chips each. The chips have 256×256 pixels, so the sensors have 256 rows and 768 columns. Each pixel is a square with a length of 55 microns. The sensor pixels are packed into Super-Pixels (SP) of size 2×4 pixels, so the sensors have 64 SP rows and 384 SP columns [START_REF] Hennessy | LHCb VELO upgrade[END_REF] [START_REF] Poikela | VeloPix: the pixel ASIC for the LHCb upgrade[END_REF].

The modules are positioned along the z axis. More information about the geometry can be found in the LHCb VELO Upgrade Technical Design Report [START_REF]LHCb VELO Upgrade[END_REF]. Figure 4 shows the format of a Super-Pixel (SP) encoded in a 32 bit integer. The less significant byte is a bitmask representing the pixels. Then we find, from bit 8 to 13, the row of the SP and, from bit 14 To take advantage of the data format, we specialized the SparseCCL algorithm to label SP instead of pixels. This allows to further reduce the amount of memory needed and to skip a decoding step. We first start by preparing the data: we remove the SP that are known to be isolated and resolve them using lookup tables, for the remaining SP, we test if there is more than one CC inside and split them if necessary. Figure 5 shows the two possible configurations for a SP: one CC or two CCs. Once the SP list is prepared, we run the algorithm using a combination of bitwise operations and a lookup table to test the adjacency. Another lookup table is used for a fast computation of the first statistical moment and the number of pixels within a SP. In 8-connectivity CCL there are eight directions of adjacency, but by using symmetries we can reduce their number to four. Figure 6 shows the configurations of SP and the pixels we have to test. Configurations a, b and c are quickly tested using only bitwise operations. While configuration d could be tested the same way, it was found faster to use a 256-entry lookup table using the dark pixels bit pattern as the address. Configuration e shows the pixels required to take the decision to split the SP in two.

IV. EXPERIMENTAL EVALUATION

Evaluating CCL algorithms has always been a challenge as the speed of such algorithms is data-dependent. To model natural images, we follow the approach which is proposed in [START_REF] Cabaret | What is the world's fastest connected component labeling algorithm[END_REF] and generate pseudo-random noise images of varying densities and granularity. The density parameter controls, at low density, the number of pixels in the image. The granularity is the size of a macro-pixel side, it controls how clustered the pixels are and thus the minimum size of a connected component.

In the case of LHCb's VELO detector, simulation data have shown that densities of hits in the sensors are very low and granularity is around 2. This is due to charge sharing in the silicon, when a particle hits the border between 2 or 4 pixels.

Table II shows the time, in microseconds to process one image of 768 × 256 pixels, for state-of-the-art dense CCL algorithms [START_REF] Cabaret | What is the world's fastest connected component labeling algorithm[END_REF] and sparse CCL algorithms. We observe, that while these algorithms are well optimized, a simple flood fill looking only at active pixels is 10 times faster at a density of 1%. In this benchmark, we measure the time in cycles per pixel (cpp). Normalizing by the number of active pixels allows to see the real impact of the increasing density: the more the connection of pixels impacts the speed, the bigger the slope of the plot will be. Normalizing by the frequency of the machine allows to abstract the frequency of the CPU for a better comparison of architectures.

We evaluate four algorithms. The first one is a flood fill algorithm as described in section II-A. While the flood fill is generally inefficient on dense images, we found that it outperforms fast implementations of iterative and two-pass algorithms on sparse images. This is due to its ability to use the pixel list information as a starting point for its connected component mapping. The other three algorithms we evaluated are variant of our algorithm: SparseCCL. The ordered by row variant is the simple case where the input is an ordered list of single pixels. The row table variant is the algorithm described in section III-B that can take an unordered list of single pixels as input. The last variant is the specialization of the algorithm for super-pixel encoding described in section III-C where we assume the list of super-pixel ordered.

Figure 7 shows the measured time in cpp for the four algorithms, for 364 × 768 pixels images of varying densities from 0% to 2.5% and a granularity g=1. The number of pixel n is given by the formula n = d 100 × w × h, where d is the density, w the number of columns and h the number of rows. In our test configuration, the number of pixels ranges from 0 to 6988. We observe that for a granularity of 1, the ordered SparseCCL working on pixels has the best behavior at low density. The Super-Pixels variant is slower at low densities because each SP is more likely to contain only 1 pixel. It presents no advantage over the pixel versions. The row table variant starts higher than the ordered one because of the cost of table reset. The flood fill algorithm is significantly slower than other version at low density, but scales better with the number of pixels and eventually becomes faster for densities > 1.8% on Intel and > 2.5% on AMD. Floodfill (sparse) SparseCCL (ordered) SparseCCL (row table) SparseCCL (Super-Pixels) Fig. 7. Cycles Per Pixel (cpp) for the four algorithms depending on the density of the image, with a granularity g=1 on an Intel Xeon Gold 6126 @2.6GHz

Figure 8 shows the same algorithms, but with a granularity g=2. The flood fill algorithm benefits from the data locality induced by the increased granularity. On the contrary, the pixel based SparseCCL variants are slowed down by it as the number of tests they have to perform is slightly higher. Thanks to the Super-Pixel encoding, the last variant of SparseCCL specialized for the LHCb experiment speeds up. Floodfill (sparse) SparseCCL (ordered) SparseCCL (row table) SparseCCL (Super-Pixels) Fig. 8. Cycles Per Pixel (cpp) for the four algorithms depending on the density of the image, with a granularity g=2 on an Intel Xeon Gold 6126 @2.6GHz Following the method described in [START_REF] Hennequin | Designing efficient simd algorithms for direct connected component labeling[END_REF], we wrote an AVX-512 version of sparseCCL. The SIMD version of the algorithm uses AVX-512's scatter and gather instructions to implement the union-find operations. For low densities, the cost of control flow is greater than the gain due to SIMD parallelism, making this SIMD version slower than the scalar one.

V. CONCLUSION

This article introduced a new parameterizable connected component labeling and analysis algorithm for sparse images, and a specialization of this algorithm for a practical application. We studied the impact of different versions on images of varying densities and granularities and showed that at low densities our algorithm performed better than state-of-the-art dense CCL algorithms and flood fill optimized for sparse images.

Fig. 1 .

 1 Fig. 1. Example of 8-connected CCL with Rosenfeld algorithm: binary image (top), image of temporary labels (bottom left), image of final labels (bottom right) after the transitive closure of the equivalence table.

9 exAlgorithm 2 :

 92 ← min + (ra, r b , rc, r d ) ; 10 if (ra = 0 and ra = ex) then U nion(L, ex, ra); 11 if (r b = 0 and r b = ex) then U nion(L, ex, r b ); 12 if (rc = 0 and rc = ex) then U nion(L, ex, rc); 13 if (r d = 0 and r d = ex) then U nion(L, ex, r d ); Sequential solve of equivalences (step 2)

Algorithm 3 :

 3 find(L, e)Input: e a label, L an equivalence tableResult: r, the root of e 1 r ← e 2 while L[r] = r do 3 r ← L[r]4 return r Algorithm 4: union(L, e 1 , e 2 ) Input: e 1 , e 2 two labels, L an equivalence table Result: e, the least common ancestor of the e's 1 if e 1 < e 2 then 2 e ← e 1 , L[e 2 ] ← e 3 else 4 e ← e 2 , L[e 1 ] ← e

Fig. 2 .

 2 Fig.2. Sparse binary image and its list representation. Each entry in the list is a tuple of pixel coordinates (x, y). The list is sorted in column major order and contains n = 9 pixels.

8 ai

 8 ←union(L, ai, find(L, j)) 9 else if is_far_enough(pixel[i], pixel[j]) then 10 start j ← start j + 1 11 // Second scan: transitive closure and analysis 12 labels ← 0 13 for i ← 0 to n -1 do 14 if L[i] = i then 15 labels ← labels + 1 16 l ← labels 17 init_features(l, pixel[i]) 18 else 19 l ← L[L[i]] 20 accumulate_features(l, pixel[i]) 21 L[i] ← l

Fig. 3 .

 3 Fig. 3. CERN LHCb VELO geometry to 22, the column of the SP. The 31 th bit is a flag indicating if the SP is isolated, ie. if it doesn't have any neighbor. The SP are delivered in raw banks. There is one raw bank per sensor and each one contains the number of SP in the bank followed by the encoded SP.

Fig. 4 .

 4 Fig. 4. CERN LHCb VELO Super-pixel format

Fig. 5 .

 5 Fig. 5. A Super-Pixel containing one CC (left) and a Super-Pixel containing two CCs, split in two (right)

Fig. 6 .

 6 Fig. 6. a. Diagonal "forward" link, b. Diagonal "backward" link, c. Vertical link, d. Horizontal link, e. Two clusters condition

  shows an example of such structure.

	Row	Nrow	Pixels (index, column)
	0	1	(5, 1)
	1..3	0
	4	2	(9, 8), (7, 10)
	5	2	(4, 9), (1, 10)
	6..11 0
	12	1	(8, 12)
	13	2	(3, 3), (6, 11)
	14	1	(2, 4)
	15	0
			TABLE I
	STRUCTURE REPRESENTING THE RECEIVED PIXELS

ACKNOWLEDGEMENTS

The authors would like to thank the LHCb computing and simulation teams for their support and for producing the simulated LHCb samples used in the paper. VVG is supported by ERC-CoG-724777 "RECEPT".