Renaud Vilmart
email: renaud.vilmart@loria.fr

A Near-Optimal Axiomatisation of ZX-Calculus for Pure Qubit Quantum Mechanics

Recent developments in the ZX-Calculus have resulted in complete axiomatisations first for an approximately universal restriction of the language, and then for the whole language. The main drawbacks were that the axioms that were added to achieve completeness were numerous, tedious to manipulate and lacked a physical interpretation. We present in this paper two complete axiomatisations for the general ZX-Calculus, that we believe are optimal, in that all their equations are necessary and moreover have a nice physical interpretation.

Introduction

The ZX-Calculus is a powerful graphical language for quantum computing and reasoning [START_REF] Coecke | Interacting quantum observables: categorical algebra and diagrammatics[END_REF]. The objects manipulated are open graphs, also called diagrams, that represent quantum evolutions through the standard interpretation. One of the most important features of the language is that the graphs can be considered unoriented, that is, any two isomorphic graphs will yield the same result. Isomorphism between diagrams are not the only transformations that preserve the interpretation though, so the ZX-Calculus comes equipped with a set of axioms: transformations between diagrams that, when applied locally, preserve the interpretation.

The language is universal: any 2 n × 2 m matrix can be represented by a ZX-diagram with respect to the standard interpretation. Hence, it has already been used in numerous applications [START_REF] Coecke | Picturing Quantum Processes: A First Course in Quantum Theory and Diagrammatic Reasoning[END_REF], ranging from measurement-based quantum computing [START_REF] Duncan | A Graphical Approach to Measurement-BasedQuantum Computing[END_REF][START_REF] Duncan | Rewriting measurement-based quantum computations with generalised flow[END_REF][START_REF] Horsman | Quantum picturalism for topological cluster-state computing[END_REF] and quantum codes [START_REF] De | The ZX calculus is a language for surface code lattice surgery[END_REF][START_REF] Chancellor | Graphical Structures for Design and Verification of Quantum Error Correction[END_REF][START_REF] Duncan | Verifying the Smallest Interesting Colour Code with Quantomatic[END_REF][START_REF] Duncan | Verifying the Steane code with Quantomatic[END_REF], to protocols [START_REF] Hillebrand | Quantum Protocols involving Multiparticle Entanglement and their Representations[END_REF] and foundations [START_REF] Backens | A complete graphical calculus for Spekkens' toy bit theory[END_REF][START_REF] Duncan | Interacting Frobenius Algebras Are Hopf[END_REF]. The language itself can be manipulated through tools such as Quantomatic [START_REF] Kissinger | Quantomatic[END_REF] or PYZX [START_REF] Kissinger | PyZX[END_REF].

A broader use of the ZX-Calculus was limited though, because of a question that remained open for a while: completeness. The language would be complete if, for any two diagrams that represent the same quantum evolution, they could be transformed into one another by mere application of the axioms. The question has been answered for gradually more expressive restrictions of the language. In 2014, complete axiomatisations were provided for the stabiliser [START_REF] Backens | The ZX-calculus is complete for stabilizer quantum mechanics[END_REF] and the real stabiliser [START_REF] Duncan | Pivoting makes the ZX-calculus complete for real stabilizers[END_REF], then for the one-qubit Clifford+T case [START_REF] Backens | The ZX-calculus is complete for the single-qubit Clifford+T group[END_REF]. However, none of these restrictions are approximately universal. The first complete axiomatisation for an approximately universal restriction -the manyqubit Clifford+T -was recently provided [START_REF] Jeandel | A Complete Axiomatisation of the ZX-Calculus for Clifford+T Quantum Mechanics[END_REF], and soon followed two complete axiomatisations for the general -universal -ZX-Calculus [START_REF] Hadzihasanovic | Two Complete Axiomatisations of Pure-state Qubit Quantum Computing[END_REF][START_REF] Jeandel | Diagrammatic Reasoning Beyond Clif-ford+T Quantum Mechanics[END_REF].

Up to the one-qubit Clifford+T case, all the axioms provided were natural and had a relevant interpretation, however, the axiomatisations for (approximately) universal ZX-Calculus introduced rules that are hard to manipulate, mainly because of their size, and that moreover can not be naturally justified.

We give in this paper a simpler axiomatisation of the general ZX-Calculus, and prove that it is complete for the general ZX-Calculus. It is basically composed of the axioms that make the Clifford -or stabiliser -fragment complete, and of an additional axiom, denoted (EU):

= β2 β1 β3 α2 α1 α3 π γ (EU)
with a side condition that links the angles on the right to those on the left. In ZX-Calculus, the green node with angle α represents a rotation of angle α around the Z axis (denoted R Z (α)), and the red one a rotation around the orthogonal axis, X (denoted R X (α)). This axiom, which is an application of the Euler angles, essentially gives a normal form for one-qubit unitaries, as a sequences of rotations around the axes X, Z and X again. This equality between diagrams has been used in [START_REF] Schröder De Witt | The ZX-calculus is incomplete for quantum mechanics[END_REF] to prove that the then version of ZX-Calculus was not complete, and is part of the axiomatisation of [START_REF] Coecke | ZX-Rules for 2-qubit Clifford+T Quantum Circuits[END_REF].

To prove that the new axiomatisation is complete, we simply derive the rules of the former axiomatisation [START_REF] Jeandel | A Complete Axiomatisation of the ZX-Calculus for Clifford+T Quantum Mechanics[END_REF]. However, since all the power of "beyond-Clifford" is contained in the rule (EU), we will end up using it a lot, which would cause a lot of side computation, for the angles on one side of the rule are not defined from the others in a linear fashion. So to avoid having to go through all this tedious process, we use another kind of normal form for ZX-diagrams, which is the graphical version of the singular-value decomposition of a matrix. Hence, instead of showing that a sound equation is derivable, we will show that we can transform the diagrams on both sides into a particular form, which is essentially unique.

We also provide a second axiomatisation, which is not very far from the other. Indeed, in the first, we may notice a rule (HD) that we call the Euler decomposition of Hadamard, which essentially gives the unitary normal form of the Hadamard gate. The second axiomatisation replaces the rules (HD) and (EU) by another single rule that unifies them.

In Section 2, we formally introduce the language ZX-Calculus, as well as the two aforementioned axiomatisations, and we discuss their minimality. In Section 3, we recover a known complete axiomatisation for the Clifford fragment, hence directly giving us access to already proven lemmas from it. In Section 4, we introduce the singular-value decompositions of cycle-free 0 → 1 and 1 → 1 ZX-diagrams, and show they are essentially unique. Finally, in Section 5, we use these decompositions to show the completeness of the axiomatisations for the Clifford+T and for the unrestricted ZX-Calculus.

ZX-Calculus

In this section, we introduce the ZX-diagrams together with a new simple axiomatisation that we prove complete in the following sections. The definition of the ZX-diagrams and their interpretation is standard.

Diagrams and standard interpretation

A ZX-diagram D : k → l with k inputs and l outputs is generated by: R .

(n,m) Z (α) : n → m α • • • • • • n m R (n,m) X (α) : n → m α • • • • • • n m H : 1 → 1 e : 0 → 0 I : 1 → 1 σ : 2 → 2 : 2 → 0 η : 0 → 2 where n, m ∈ N, α ∈ R,
D 1 ⊗ D 2 := D 1 ⊗ D 2 D 2 • D 1 := D 2 • D 1 := 1 := 1 0 0 1 := 1 √ 2 1 1 1 -1 :=     1 0 0 0 0 0 1 0 0 1 0 0 0 0 0 1     :=     1 0 0 1     := 1 0 0 1 α := 1 + e iα α • • • • • • n m := 2 m              2 n        1 0 • • • 0 0 0 0 • • • 0 0 . . . . . . . . . . . . . . . 0 0 • • • 0 0 0 0 • • • 0 e iα        n + m > 0
For any n, m ≥ 0 and α ∈ R:

α • • • • • • n m = ⊗m • α • • • • • • n m • ⊗n where M ⊗0 = 1 and M ⊗k = M ⊗ M ⊗k-1 for any k ∈ N * .
To simplify, the red and green nodes will be represented empty when holding a 0 angle:

• • • 0 := • • • • • • • • • and 0 := • • • • • • • • • • • • ZX-Diagrams are universal: ∀A ∈ C 2 n × C 2 m , ∃D : n → m, D = A
However, it is customary to restrict the language to a countable or finite set of angles. Some of these restrictions, or fragments, correspond to well-known restrictions of quantum computing: The π 2 -fragment -the restriction where all the angles are multiples of π 2 -corresponds to Clifford; while the π 4 -fragment corresponds to Clifford+T. In the following, we may refer to the π 2 -fragment using the term Clifford, and similarly for the π 4 -fragment.

Calculus

The diagrammatic representation of a matrix is not unique in the ZX-Calculus. As a consequence the language comes with a set of axioms. Additionally to the axioms of the language described in Figure 1, one can:

bend any wire of a ZX-diagram at will, without changing its semantics. This paradigm -the so-called Only Connectivity Matters -can be derived from the following axioms: ) denote one or more wires. In rule (EU), β1, β2, β3 and γ can be determined as follows:

= = = = = = = = • • • = α+β β • • • α • • • (S) • • • • • • • • • • • • = (Ig) = (Ir) -π 4 π 4 = (E) = (CP) = (B) = β2 β1 β3 α2 α1 α3 π γ (EU) π 2 π 2 -π 2 = (HD) α • • • = α • • • • • • • • • (H)
x + := α 1 +α 3 2 , x -:= x + -α3, z := cos α 2 2 cos x + + i sin α 2 2 cos x - and z := cos α 2 2 sin x + -i sin α 2 2 sin x -, then β1 = arg z + arg z , β2 = 2 arg i + z z , β3 = arg z -arg z , γ = x + -arg(z) + α 2 -β 2 2
where by convention arg(0) := 0 and z = 0 =⇒ β2 = 0.

apply the axioms to sub-diagrams. If ZX D 1 = D 2 then, for any diagram D with the appropriate number of inputs and outputs:

• ZX D 1 • D = D 2 • D • ZX D • D 1 = D • D 2 • ZX D 1 ⊗ D = D 2 ⊗ D • ZX D ⊗ D 1 = D ⊗ D 2
where ZX D 1 = D 2 means that D 1 can be transformed into D 2 using the axioms of the ZX-Calculus.

All the axioms of Figure 1, but (EU), are standard in the ZX-calculus. Roughly speaking: (S) and (I) correspond to the axiomatisation of orthonormal basis [START_REF] Coecke | A new description of orthogonal bases[END_REF], each color being associated with an orthonormal basis; (CP) and (B) capture the fact that the two bases are strongly complementary [START_REF] Coecke | Interacting quantum observables: categorical algebra and diagrammatics[END_REF]; (H) means that Hadamard can be used to exchange the colours and (HD) means that Hadamard can be decomposed using π 2 -rotations [START_REF] Duncan | Graphs States and the necessity of Euler Decomposition[END_REF]; (E) states that some particular scalars (ZX-diagram with no input/output) can vanish, which means that their interpretation is one [START_REF] Jeandel | ZX-Calculus: Cyclotomic Supplementarity and Incompleteness for Clifford+T Quantum Mechanics[END_REF]. In the following we investigate the properties of (EU).

The Euler Angles

The rule (EU) is really all about unitaries. Indeed, we have the following result:

Proposition 1. Any one-qubit unitary can be decomposed as e iγ R Z (α 3 )R X (α 2 )R Z (α 1 ), which can be represented in ZX as:

α2 α1 α3 π γ
If the unitary is not diagonal or anti-diagonal (i.e. if α 2 = 0 mod π), then this decomposition can be made unique if we impose α 1 ∈ [0, π)

In 1775, Euler proved what is now called Euler's rotation theorem [START_REF] Euler | Formulae Generales Pro Translatione Quacunque Corporum Rigidorum[END_REF], stating that there are several ways to decompose a rotation into several rotations around elementary axes. In quantum mechanics, a consequence is that any unitary operator on one qubit can be seen as either a composition of rotations around Z, X, Z; or around X, Z, X. On the one hand, the rule (HD) says -in a distorted, ZX-style way -that the Hadamard gate can be decomposed as a series of rotations, while on the other hand, the rule (EU) gives the equality between two different decompositions of the same unitary:

= β2 β1 β3 α2 α1 α3 π γ
(EU) where

                   x + := α1+α3 2 x -:= x + -α 3 z := cos α2 2 cos (x + ) + i sin α2 2 cos (x -) z := cos α2 2 sin (x + ) -i sin α2 2 sin (x -) β 1 = arg z + arg z β 2 = 2 arg i + z z β 3 = arg z -arg z γ = x + -arg(z) + α2-β2 2
The angles β i and γ seem to not always be defined. Indeed, arg is not defined in 0, and β 2 is not defined when z = 0. By convention, we decide that arg(0) = 0 and that β 2 = 0 when z = 0.

The first proof of incompleteness [START_REF] Schröder De Witt | The ZX-calculus is incomplete for quantum mechanics[END_REF] relied on an euler decomposition, but adding it to the set of ZX axioms has been avoided for a while because of its non-linearity. However, a non-linear axiom is necessary to get the completeness for the general ZX-Calculus [START_REF] Jeandel | Diagrammatic Reasoning Beyond Clif-ford+T Quantum Mechanics[END_REF]. And so, it has been used in [START_REF] Coecke | ZX-Rules for 2-qubit Clifford+T Quantum Circuits[END_REF] to prove the completeness of the 2-qubit π 4 -fragment of the ZX-Calculus. The rule (EU) is actually much more powerful than this, for, as we will prove in the following: Theorem 1. The ZX-Calculus -with axioms in Figure 1 -is complete for pure qubit quantum mechanics. For any two diagrams D 1 and D 2 of the ZX-Calculus:

D 1 = D 2 ⇐⇒ ZX D 1 = D 2

On Minimality

We call an axiomatisation minimal when there is no redundancy in the axioms. Particularly, we want a proof that none of the axioms are derivable from the others. We conjecture that all the axioms in Figure 1 are necessary. Indeed, in [START_REF] Backens | Towards a Minimal Stabilizer ZX-calculus[END_REF], nearly all the rules for Clifford -i.e. all of the axioms in Figure 1 except (E) and (EU)-are proven to be necessary, and all arguments stand here:

-(S): It is the only axiom that can transform a node of degree four or higher into a diagram containing lower-degree nodes -(I g ) or (I r ): These are the only two axioms that can transform a diagram with nodes connected to a boundary to a node-free diagram -(CP): It is the only axiom that can transform a diagram with two connected outputs into one with two disconnected outputs -(HD): The necessity of this axiom requires a non-trivial interpretation given in [START_REF] Duncan | Graphs States and the necessity of Euler Decomposition[END_REF][START_REF] Duncan | Pivoting makes the ZX-calculus complete for real stabilizers[END_REF], and given again in the Appendix at page 15. -(H): It is the only axiom that matches red nodes with 4+ degree to green nodes of the same degree However, (E) and (EU) can also be proven to be necessary:

- In a nutshell, all the axioms are proven to be necessary, except (B) and one of the (I). Another aspect of minimality, is whether a rule can be made "simpler" thanks to the others, according to some measure, be it arbitrary or well-defined. In the previous axiomatisation, we have two rules that are closely related to how unitaries can be decomposed: (HD) and (EU). It so happens that we can fuse them into one, of the same size as (EU), and doing so allows us to simplify the scalar rule:

Theorem 2.      π 2 π 2 -π 2 = (HD) , = β2 β1 β3 α2 α1 α3 π γ (EU) , -π 4 π 4 = (E)      can be replaced by      = β2 β1 β3 α1 α2 π γ (EU') , α = (IV)      where                    x + := α1+α2 2 x -:= x + -α 2 z := -sin (x + ) + i cos (x -) z := cos (x + ) -i sin (x -) β 1 = arg z + arg z β 2 = 2 arg i + z z β 3 = arg z -arg z γ = x + -arg(z) + π-β2 2
As a consequence, the axiomatisation given in Figure 2 is complete for universal quantum mechanics. ) denote one or more wires. In rule (EU'), β1, β2, β3 and γ can be determined as follows:

• • • = α+β β • • • α • • • (S) • • • • • • • • • • • • = (Ig) = (Ir) α = (IV) = (CP) = (B) α • • • = α • • • • • • • • • (H) = β2 β1 β3 α1 α2 π γ (EU')
x + := α 1 +α 2 2 , x -:= x + -α2, z := -sin x + + i cos x -and z := cos x + -i sin x -, then β1 = arg z+arg z , β2 = 2 arg i + z z , β3 = arg z-arg z , γ = x + -arg(z)+ π-β 2 2
where by convention arg(0) := 0 and z = 0 =⇒ β2 = 0.

Proof. The proof, done at the end of the appendix, at page 25, consists in showing that all the rules in Figure 1 are derivable.

On the one hand, this new axiomatisation is one axiom shorter, and (EU') and (IV) can be considered simpler than (EU) and (E). On the other hand, the axiomatisation in Figure 1 has the nice property that it suffices to remove (EU) and (E) to get a complete axiomatisation for the scalar-free Clifford fragment. Moreover, (EU) is arguably more natural, and has already been given for instance in [START_REF] Coecke | ZX-Rules for 2-qubit Clifford+T Quantum Circuits[END_REF].

The following of the paper is dedicated to the proof of Theorem 1. Since [START_REF] Hadzihasanovic | Two Complete Axiomatisations of Pure-state Qubit Quantum Computing[END_REF][START_REF] Jeandel | Diagrammatic Reasoning Beyond Clif-ford+T Quantum Mechanics[END_REF] provided us with two complete axiomatisations for the general ZX-Calculus, all we have to do is prove all the equations used as axioms in either one of these two axiomatisations. As the axiomatisation in [START_REF] Hadzihasanovic | Two Complete Axiomatisations of Pure-state Qubit Quantum Computing[END_REF] requires additional generators and more axioms, we will use the axiomatisation of [START_REF] Jeandel | Diagrammatic Reasoning Beyond Clif-ford+T Quantum Mechanics[END_REF] as a reference which consists in all the axioms of Figure 1 but (EU), together with the following axioms, we call obsolete, as we are proving in the following that they can be derived using the (EU) rule:

Obsolete ZX-rules = π α -α π α π (K) α α+π = 2α+π (SUP) β α π β γ -γ α = α α π β -γ γ β (C) π 4 π 4 π 4 -π 2 π 4 π 4 π 4 = π 4 π π 2 π 4 π 4 π π 4 (BW) θ2 θ1 α -α β -β = γ -γ 2e iθ3 cos (γ) = e iθ1 cos (α) + e iθ2 cos (β) π 2 π 4 π 4 θ3 π 4 π 4 (A) π = (Z0) π = (IV')
Remark 1. The last two equations, (ZO) and (IV'), are actually derivable from (K), (SUP) and the Clifford axiomatisation [START_REF] Jeandel | ZX-Calculus: Cyclotomic Supplementarity and Incompleteness for Clifford+T Quantum Mechanics[END_REF]. However, they are given here, because together with (S), (I), (CP), (B), (HD) and (H), they make the Clifford fragment complete, which will be our first milestone.

Clifford

As we just said, a first and easy step to do is to show that we can recover the rules that are known to make the language complete for Clifford [START_REF] Backens | Towards a Minimal Stabilizer ZX-calculus[END_REF]. This will allow us to freely use in the following all the equations of the π 2 -fragment that are sound. We already have most of these rules that make the ZX-calculus complete for Clifford. We only lack two: the zero (ZO) and the inverse (IV') rules. A first very well known lemma we will use for both proofs is the Hopf law: Lemma 1 (Hopf Law).

ZX =

From there, it is fairly easy to recover the inverse rule:

Proposition 2. The inverse rule is derivable:

ZX =
To prove the zero rule, we will use another well known equation, called π-commutation, which is also one of the now obsolete rules. Proposition 3. The π-commutation is derivable:

ZX π α = -α π α π
and, with some effort, the rule (ZO), which only deals with null diagrams, can be recovered: Proposition 4. The zero rule is derivable:

ZX π = π
As a result: Theorem 3. For any diagrams D 1 , D 2 of the π 2 -fragment:

D 1 = D 2 ⇐⇒ ZX D 1 = D 2
From this first milestone, we get all the sound equations in Clifford, but actually also a bit more. For instance, the following lemmas are known to be derivable from the Clifford axiomatisation (see Appendix): The next step is logically to get the completeness for Clifford+T quantum mechanics, i.e. the completeness of the π 4 -fragment of the ZX-calculus. Now that we are seeking to prove equations that are out of Clifford, we will begin to use (EU) to its full potential. However, we would like, as much as possible, to avoid computing the angles, because, since we work on the problem of completeness, we need to formally prove the equality between two diagrams, and hence to formally write what the angles resulting from (EU) are, which becomes tedious after a few number of application of the rule.

Lemma 2.
To simplify this task, instead of showing directly that two diagrams can be turned into one another, we will define a normal form for them, show that it is unique, and show that there is an algorithm to turn them in this normal form.

First, we show another version of the rule (EU) with dangling branches:

Lemma 4. α2 α1 α3 = γ β2 π β3 β1
where β 1 , β 2 , β 3 , γ can be determined as in rule (EU).

In a particular case, it implies:

Corollary 1. α1 α3 = γ β2 π β3 β1 π 2 -π 2
where β 1 , β 2 , β 3 , γ can be determined as in rule (EU) with α 2 ← π 2 . We can also derive a kind of inverse operation:

Lemma 5. α2 α1 = β3 β1 π 2 β2-π 2 π γ -π 2
where β 1 , β 2 , β 3 , γ can be determined as in rule (EU) applied with the angles α 2 ← α 2 + π 2 and

α 3 ← π 2 .
Then, we show that any diagram in the form of the left hand side of (SUP)-but with arbitrary angles -can be transformed in a state with no branching: Lemma 6. Right now, we have proven all the equations that do not really need a unique normal form. For the rest, we present the singular-value decomposition of a matrix, and introduce it to ZX-diagrams. Definition 1. We call a singular value decomposition (SVD) of a matrix a decomposition of the form

M = U ΣV †
where U and V are unitary, and Σ is diagonal. Notice that M needs not be square (in this case Σ has the same dimensions as M ).

To justify the use of SVDs, we give some of it interesting properties [START_REF] Roger | Positive definite matrices[END_REF]: Proposition 6. The SVD M = U ΣV † of a matrix M has the following properties:

-It exists whatever M -Σ can be made unique if we impose that its diagonal entries are decreasing non-negative real numbers -U and V are not unique in general, though:

-If M is square with distinct and non-zero singular values, then U and V are essentially unique:

U ΣV † = U ΣV † ⇐⇒ (∃d, (U = U d) ∧ (V = V d))
where d is diagonal with diagonal entries some roots of unity.

Even though the singular-value decomposition is relevant for any diagram, we are only going to give its derivation for a particular family of diagrams: Definition 2. We call a cycle-free diagram a diagram composed only of , ,

α • • • n , α • • • n
where n ∈ N and α ∈ R.

Remark 3. Some diagrams that do not strictly follow the conditions of the previous definition will still be considered cycle-free if they are equal to a cycle-free diagram by mere application of the "only topology matters" paradigm, i.e. if they are isomorphic to a cycle-free diagram. E.g.:

α = α
We can now easily give a normal form for one-qubit states, using the SVD. In this case Σ = s |0 , U is a one-qubit unitary, which can be expressed as in Proposition 1, and V is a 0 → 0 unitary, i.e. a global phase. where γ ∈ [0, π 2 ], and α 1 , α 5 , α 1 , α 5 ∈ [0, π). We denote the two forms respectively SVD g and SVD r . Remark 4. We gave two conventions for the SVDs of 0 → 1 and 1 → 1 diagrams. These two depend on the basis in which we consider the decomposition. SVD g corresponds to the computational basis, while SVG r corresponds to the diagonal basis. If M = U ΣV † with Σ diagonal in the computational basis, M = (U H)

• HΣH • (V H) † .
-γ = γ = 0 -γ = γ = π 2 -α i = β i and γ = γ

Clifford+T and Beyond

The point now is to exploit the SVD of ZX-diagrams and their uniqueness. A rule that can directly use these results is (BW):

Proposition 11. ZX π 4 π 4 π 4 -π 2 π 4 π 4 π 4 = π 4 π π 2 π 4 π 4 π π 4
The results on SVDs can not be directly used to prove the equation (C), for its diagrams have 4 inputs/outputs, and have a cycle. However, the SVDs can be used to prove a first intermediary result: Lemma 7. Remark 5. From Lemma 7, (C) can be derived using only the Clifford rules. However, the provided proof requires using half angles. Hence, whenever the considered fragment contains all its half angles, the equation in Lemma 7 should be preferred to (C).

We have derived all the rules necessary for the completeness of the Clifford+T fragment of the ZX-Calculus, which means: Theorem 4. For any diagrams D 1 , D 2 of the π 4 -fragment:

D 1 = D 2 ⇐⇒ ZX D 1 = D 2
Finally, it remains to derive the equation (A). Notice that the diagram on the left hand side contains a cycle, which implies we can not use the results on SVDs. However, the cycle can be easily removed, and we are able to prove:

Proposition 13. ZX θ2 θ1 α -α β -β = γ -γ 2e iθ3 cos (γ) = e iθ1 cos (α) + e iθ2 cos (β) π 2 π 4 π 4 θ3 π 4 π 4
This last proposition ends the proof of Theorem 1.

Discussion and Further Work

We have provided two simple but complete axiomatisations of the ZX-Calculus for universal quantum mechanics. By doing so, we have restored intuitiveness -one of its the first aims -to the language (at least on the structural level, computing the angles in (EU) remains tedious if done formally). This step forward should simplify axiom-related problems such as verification or compilation.

To simplify the task of proving the derivability of equations, we introduced singular-value decomposition of 0 → 1 and 1 → 1 diagrams, and proved that there exists an algorithm to turn any 0 → 1 and 1 → 1 cycle-free diagram into its SVD form. We did not need SVD form for diagrams with cycle, and leave as a further development the extension of the algorithm to arbitrary 0 → 1 and 1 → 1 diagrams, which should be possible by completeness and universality.

We did not need to define the SVD form for larger diagrams either. A problem would arise in ZX, for instance for a diagram with 3 inputs/outputs: do we decompose the diagram as a 0 → 3, or a 1 → 2 diagram and then use the map/state duality? This would result in two completely different decompositions. Still, defining SVDs for diagrams of any arity could prove interesting.

Concerning the result itself, we have proven that, in ZX-Calculus: many-qubit Clifford completeness + completeness for 1-qubit unitaries

   = many-qubit completeness
This formulation is a bit excessive, since we actually have several rules that operate beyond the Clifford fragment, namely (S) and (H), where the angles can take any value in R -and this feature is actually needed for the completeness. Still, since it is not absurd to imagine we can always find similar rules for the considered language, this raises two questions:

-Is it true for fragments of the ZX-Calculus?

The answer in general is no. Indeed, in the case of Clifford+T, the axiomatisation for Clifford is enough to get the 1-qubit completeness [START_REF] Backens | The ZX-calculus is complete for the single-qubit Clifford+T group[END_REF]. However, it has been proven that rules (SUP) and (E) are necessary [START_REF] Jeandel | ZX-Calculus: Cyclotomic Supplementarity and Incompleteness for Clifford+T Quantum Mechanics[END_REF][START_REF] Perdrix | Supplementarity is Necessary for Quantum Diagram Reasoning[END_REF]. Hence, the previous statement does not stand for Clifford+T. -How far from this statement are we in other languages?

For instance, we know a complete presentation for the many-qubit Clifford fragment of quantum circuits [START_REF] Selinger | Generators and relations for n-qubit Clifford operators[END_REF]. Moreover, the rule (EU) has an obvious equivalent in circuits, and is the only needed axiom for 1-qubit completeness. So what do we lack to get the universal completeness?

A Appendix

Proof (Prop. 1).

• Existence: Any element of U (2) can be decomposed as:

e iϕ/2 e iψ0 0 0 e -iψ0 cos (θ) sin (θ) -sin (θ) cos (θ)

e iψ1 0 0 e -iψ1
Hence, the existence is given by:

α2 α1 α3 π γ = e i(γ+ α 2 2 ) 1 0 0 e iα3 cos α2 2 -i sin α2 2 -i sin α2 2 cos α2 2 1 0 0 e iα1 = e i(γ+ α 2 2 ) 1 0 0 e i(α3+ π 2 ) cos α2 2 sin α2 2 -sin α2 2 cos α2 2 1 0 0 e i(α1-π 2 ) = e i(γ+ α 1 +α 2 +α 3 2 ) e -i( α 3 2 + π 4 ) 0 0 e i( α 3 2 + π 4 ) cos α2 2 sin α2 2 -sin α2 2 cos α2 2 e -i( α 1 2 -π 4 ) 0 0 e i( α 1 2 -π 4 )
• Uniqueness:

Suppose α2 α1 α3 π γ = α 2 α 1 α 3 π γ
. The first diagram yields:

e i(γ+ α 2 2 ) cos α2 2 -ie iα1 sin α2 2 -ie iα3 sin α2 2 e i(α1+α3) cos α2 2
and similarly for the second one. If α 2 = 0 mod π, then neither cos α2 2 nor sin α2 2 is null. Hence, dividing element (1,1) by element (0,0) on both sides gives e i(α1+α3) = e i(α 1 +α 3 ) and dividing element (0,1) by element (1,0) on both sides gives e i(α1-α3) = e i(α 1 -α 3 ) . In other words, α 1 +α 3 = α 1 +α 3 mod 2π and α 1 -α 3 = α 1 -α 3 mod 2π, so 2α 1 = 2α 1 mod 2π i.e. α 1 = α 1 mod π. Since we required α 1 , α 1 ∈ [0, π), we get α 1 = α 1 . It then follows easily that α 3 = α 3 , α 2 = α 2 and γ = γ .

Proof (Necessity of Rule (HD)). We define the non-standard interpretation . as follows:

→ → → α • • • • • • n m → α α • • • • • • α • • • • • • n m → α α • • • • • • → D 1 • D 2 → D 1 • D 2 D 1 ⊗ D 2 → D 1 ⊗ D 2
It is then easy to see that all the rules but (HD) hold under this interpretation, hence proving that (HD) could not be derived from the other rules.

A.1 Proofs for Clifford

Proof (Lemma 1).

= = (I) (S) = (B) = (CP) = (S) (I)
Proof (Prop. 2). (2)

= (E)

α π = 1 α π = (S) α 2 π α 2 = 3 α 2 -α 2 α 2 π π = (S) α 2 π π = 2 1 π π α 2
(3)

α π π = (3) π π α 2 π = (2) π α 2 π = 2 α 2 π π (4)
Now, if α ∈ Dπ (where D := Z 1 2 ), then there exists n such that 2 n α = 0 mod 2π. Hence, in this case the scalar on the right hand side of (3) can be removed by applying (4) from right to left n + 1 times then using (2) and 2 to remove it. Hence:

∀α ∈ Dπ, π = α π (5) π = (5) π π 4 -π 4 = (E) π (6) -π 2 = (H) -π 2 = (HD) (S) 2 π 2 -π 2 = (CP) (S) π 2 -π 2 (7) π = (5) π -π 2 = (7) π π 2 -π 2 = (5) π = (1) 2 π = (6)
π Proof (Lem. 2). We are going to use the following equation we get from the Clifford completeness:

π = π π (8) so:
α1 α3 = (S) α1 α3 -π 2 π 2 = π 2 C α3 α1 π 2 π 2 -π 2 = 4 γ β2 π β3 β1 π 2 -π 2 Proof (Lem. 5). α2 α1 = (B) α2 α1 = π 2 C α2+ π 2 α1 -π 2 π 2 π 2 = (EU) (S) β1 π γ β3+ π 2 β2 -π 2 = (H) (S) β3+ π 2 β1 β2 π γ -π 2 = π 2 C β3 β1 β2-π 2 π 2 π γ -π 2 = (B) β3 β1 π 2 β2-π 2 π γ -π 2 Proof (Lem. 6). α1 α3 = (S) α1 α3 = 1 β2 β3 β1 π 2 π γ -π 2 = (CP) (S) (I) π 2 β1+β3 β2 π γ -π 2
Proof (Prop. 5). We first use Lemma 6, where α 3 = α 1 + π. In this case, it can be computed that β 1 + β 3 = 0, so we end up with:

α+π α = β2 π γ -π 2
. From this, we can easily specify the scalar on the right part:

β2 π γ -π 2 = π 2 C β2 π γ -π 2 = α+π α = (S) (I) 2α+π
So finally:

α α+π = β2 π γ -π 2 = 2α+π = 1 2α+π
Proof (Rem. 2).

π 4 π = (E) π π 4 π 4 -π 4 = (H) π -π 4 π 4 π 4 = (HD) 3π 4 π 4 π 4 π -π 2 = 3 2 π 3π 4 π π 2 -π 2 -π 4 = 5 1 2 π π 2 π 2 π -π 2 = 2 3 π 2
Proof (Prop. 7). First, notice that a state in the previous form can easily be transformed into an SVD. Indeed, if β ∈ [π, 2π):

s β α = 3 s β-π -α π α
and similarly for the SVD r . We can show that we can transform an SVD r into an SVD g and vice-versa: : (α) can be obtained as a combination of the last two. Then :

α s β = (H) α s β = π 2 C (S) 2 s α-π 2 β-π 2 -π 2 π 4 π = (EU)
D γ = α s β γ = (S) α s β+γ D 1 D 2 = (S) α2 s 2 α1 s 1 β1+β2 = 6 β1+β2 -π 2 π γ β2 β s 1 s 2 = π 2 C γ -π 4 β2 β1+β2 + π 2 π β + π 2 s 1 s 2 = (H) π γ -π 4 s 2 β2 β1+β2 + π 2 s 1 β + π 2
Finally, the generator R (n,1) Z (α) can be obtained by composition of R Z (α) and R

(2,1) Z (α); and R (n,1) X (α) can be obtained by composition of R (n,1) Z (α) and H.

Proof (Prop. 8). The equality reads s

1 1 + e iα1 e iβ1 (1 -e iα1 ) = s 2 1 + e iα2 e iβ2 (1 -e iα2 )
. If α 1 = π, then it is easy to see that α 2 = π and s 1 e iβ1 = s 2 e iβ2 . If α i = π, then the upper coefficient is non-null, hence we can divide the lower coefficient by the upper one, which yields:

e iβ1 1 -e iα1 1 + e iα1 = e iβ2 1 -e iα2 1 + e iα2 ⇔ e iβ1 tan α 1 2 = e iβ2 tan α 2 2
If α 1 = 0 mod 2π then α 2 = 0 mod 2π. Otherwise, since

β 1 , β 2 ∈ [0, π), β 1 = β 2 and α 1 = α 2 mod 2π.
Proof (Prop. 9). First, if D is in the form SVD g , but where the constraints on the angles are not met, we can transform it into an actual SVD g :

-If α 1 ∈ [π, 2π
) (and similarly for α 5 ): The two 1 → 1 generators R (1,1) Z (α) and H can be put in SVD:

α2 α1 . . . = (S) 3 -α2 α1-π . . . π π α2 -If γ ∈ [-π 2 , 0[: α3 γ = (S) 3 π+α3 -γ π γ -If γ ∈ [-π, -π 2 [: α3 γ = (S) 3 -α3 γ+π π π π α3 -If γ ∈ [ π 2 , π): α3 γ = (
α = (S) α+ π 2 -π 2 = π 2 C 2 α+ π 2 π 2 π -π 4 = π 2 C (I) π 2 π 2 -π 2
The composition of two SVDs can be put in SVD (here, ignoring the scalars): where u, v, u and v have been chosen so that Σ and Σ are real matrices. Notice that U , V † , U , V † are unitaries. We have two SVDs that represent the same matrix: Since β 5 and α 5 are in [0, π), the representation of the unitary is unique by Proposition 1, so β 5 = α 5 , β 4 = α 4 , and x = x + ϕ 1 -ϕ 0 . Similarly, the second equation yields

U • Σ • V † = D 1 = D 2 = U • Σ • V †
α 1 = β 1 , α 2 = β 2 and β 3 -x -π 2 = α 3 -x -π 2 + ϕ 0 -ϕ 1 .
Together, the equations on x and x imply that α 3 = β 3 .

A.3 Proofs for Clifford+T and Beyond

Proof (Prop. 11). Using Proposition 9, we can put both sides of the equation in SVD, and thanks to Proposition 10, the two forms have the same angles. We can even compute:

-π 2 π s 1 π 4 π π 4 π 2 π 4 π π 4 = -π 2 π s 1 π 2 β1 π 2 γ π 2 β1 s 2 = -π 2 π s 2 π 4 π 4 π 4 -π 2 π 4 π 4 π 4 = (CP) (S) -π 2 π s 2 π 4 π 4 π 4 -π 2 π 4 π 4 π 4 = (I) (S) -π 2 π s 2 π 2 π 2 = 2 π 2 C -π 2 π s 2 π 4 π π 4 π = 2 3 2 s 2
Proof (Lem. 7). We prove the equality by simplifying both sides of the equation. The left hand side yields, when ignoring the scalars: where n and m are chosen in {0, 1} so that γ 1 + nπ and β + π 2 + mπ are in [0, π). Similarly, the right hand side yields:

β α π β α = π 2 C α+ π 2 β β+ π 2 -α π 2 = (B) -α β β+ π 2 π 2 α+ π 2 = 1 β+ π 2 -α β π 2 α+ π 2 = (EU)
β α π β α = (-1) m β3 β+ π 2 +mπ β2 γ3+ (n+m)π (-1) n γ2 γ1+nπ
Notice that, due to the symmetry of the two diagrams, the resulting scalars (that we ignored) are equal (and non null). If β 2 = 0 mod π, then we can compute that both α and β are multiples of π, and in this case the equation is trivially derivable. Else, notice that β2 is invertible, its inverse is

1 1-e 2iβ 2
1 -e iβ2 -e iβ2 1

. Hence, we get:

(-1) m β3 β+ π 2 +mπ γ3+ (n+m)π (-1) n γ2 γ1+nπ = (-1) m β3 β+ π 2 +mπ γ3+ (n+m)π (-1) n γ2 γ1+nπ
We can then plug any red dot with angle ∈]0, π 2 [, say π 4 , on the lower branch. We can now use Proposition 10, match the angles γ 1 + nπ = β + π 2 + mπ and (-1) n γ 2 = (-1) m β 3 , so the two initial diagrams are equal.

Proof (Prop. 12).

γ β β -γ α α π = (H) (CP) (B) β α α π β γ -γ -γ 2 γ 2 = (S) 7 -γ 2 α γ 2 π β -γ 2 β γ 2 α = 7 γ 2 β β -γ 2 α π -γ 2 α γ 2 = (H) (B) γ 2 α β π -γ 2 α β -γ 2 γ 2 = (H) 7 β α π α γ 2 -γ 2 γ 2 -γ 2 β = (B) (H) β α γ 2 -γ 2 β π γ 2 α -γ 2 = 7 γ β β -γ α -γ 2 α π γ 2 = (B) (CP) (H) π β α -γ γ β α
Proof (Prop. 13). The idea of the proof is here again to use the SVD, but this time of a state, the equation being between two states. The diagram on the right hand side has an SVD by direct use of Proposition 7. However, the one on the left hand side has a cycle, so we have to work a bit more. Notice that, using Proposition 9: From a direct use of Proposition 7:

γ -γ θ3 π 4 π 4 = s 2 2 δ2
Since the initial equation is sound under the constraint 2e iθ3 cos (γ) = e iθ1 cos (α) + e iθ2 cos (β), we have:

s 1 1 δ1 = s 2 2 δ2
By Proposition 8, either δ 1 = δ 2 = 0, or δ 1 = δ 2 = π, or δ 1 = δ 2 = 0 mod π and 1 = 2 . Notice however that in our case, δ 1 and δ 2 cannot be π because

γ -γ θ3 π 4 π 4 = (1+i) 1 √ 2e iθ3 cos (γ)
and

s π = 0 √ 2se i .
If δ 1 = δ 2 = 0, then 1 and 2 can be discarded: s = s . Hence, in any case, it only remains to prove that the two scalar diagrams s 1 and s 2 can be transformed into one another. Proof (Thm. 2). We are going to prove here that all the equations in Figure 1 can be recovered from the ones in 2. First of all, we try to recover a complete axiomatisation for Clifford. To do so, we simply need to show that the Hadamard gate can be decomposed.

s 1 = 3 2 s 1 π 4 π -π 4 π = 2 s 1 π 2 -π 4 π = (S) (I) s 1 γ π 4 -π 4 π -γ π 4 = (CP)
We have recovered a complete axiomatisation for Clifford. We now have access to all the lemmas deriving from this axiomatisation. The next step is to prove the equation (E) is derivable. where x is considered as a variable, and hence, all the computed angles depend on it, while the angles α i are fixed. We want to find x 0 such that β 3 (x 0 ) + γ 1 (x 0 ) = 0 mod π. Let the functions f and g be defined as: Hence, g -π 2 g π 2 ≤ 0. Since g is continuous, by the intermediate value theorem, there exists x 0 ∈ [ -π 2 , π 2 ] such that g(x 0 ) = 0. Notice now that f (x 0 ) = arctan 0 1+tan(α1) 2 cos(α2-x0) 2 = 0. Also, it can be computed that f = β 3 + γ 1 mod π. Hence, β 3 (x 0 ) + γ 1 (x 0 ) = 0 mod π i.e. β 3 (x 0 ) + γ 1 (x 0 ) = nπ. Hence, denoting β i ← β i (x 0 ) and γ i ← γ i (x 0 ): Since, thank to Proposition 1, the unitary representation is unique if β 1 + mπ ∈ [0, π[ (m has been chosen for this purpose), then the previous diagram is provably equivalent to the one resulting directly from (EU).

2 The

 2 and the generator e is the empty diagram. and the two compositions: -Spacial Composition: for any D 1 : a → b and D 2 : c → d, D 1 ⊗ D 2 : a + c → b + d consists in placing D 1 and D 2 side by side, D 2 on the right of D 1 . -Sequential Composition: for any D 1 : a → b and D 2 : b → c, D 2 • D 1 : a → c consists in placing D 1 on the top of D 2 , connecting the outputs of D 1 to the inputs of D 2 . standard interpretation of the ZX-diagrams associates to any diagram D : n → m a linear map D : C 2 n → C 2 m inductively defined as follows:

Fig. 1 .

 1 Fig. 1. Set of rules ZX for the ZX-Calculus with scalars. The right-hand side of (E) is an empty diagram. (...) denote zero or more wires, while ( • • •) denote one or more wires. In rule (EU), β1, β2, β3 and γ can be determined as follows: x + := α 1 +α 3

  E): It is the only axiom that can transform a non-empty diagram into an empty one -(EU): It is the only non-linear axiom

Fig. 2 .

 2 Fig. 2. Set of rules ZX' for the ZX-Calculus with scalars. The right-hand side of (E) is an empty diagram. (...) denote zero or more wires, while ( • • •) denote one or more wires. In rule (EU'), β1, β2, β3 and γ can be determined as follows: x + := α 1 +α 2

2 where β 1 Proposition 5 .Remark 2 .

 2152 , β 2 , β 3 , γ can be determined as in rule (EU) with α 2 ← π 2 . Now, by specialising the angles to α and α + π, we shall recover (SUP): The supplementarity is derivable: The supplementarity allows us to prove:

Proposition 7 (Proposition 8 (

 78 SVD of a State). Any cycle-free state D : 0 → 1 can be put in the following forms: β ∈ [0, π), and where s and s are 0 → 0 diagrams, i.e. scalars. We call these two forms respectively SVD g and SVD r . SVDs of states are essentially unique). If D 1 = α1 β1 s 1 and D 2 = α2 β2 s 2 are in SVD, and if D 1 = D 2 = 0, then either: α 1 = α 2 mod 2π and α i = 0 mod π α 1 = α 2 and β 1 = β 2 We can have basically the same results for 1 → 1 operators: Proposition 9 (SVD of a 1 → 1 diagram). Any cycle-free diagram D : 1 → 1 can be written in the forms:

Proposition 10 ( 1 → 1

 1011 SVDs are essentially unique). , and that D 1 = D 2 = 0. Then, either:

  From which we can deduce the equation (C) itself:

2 A. 2

 22 Proofs for Singular-Value DecompositionsProof (Lem. 4).

  show that the two decompositions are equivalent:

s

  by composition, the 1 → 1 generator R (1,1) X (α) can be put in SVD. If the 1 → 1 diagram has no cycle, there can still be branching. Hence, there can be a state D : 0 → 1 in tree-like form attached to the "main wire" by a node, say green, as follows:Branching made by a red node can be deduced by composing the green one and Hadamard nodes.Proof (Prop. 10). First we decompose D 1 and D 2 as:

First 2 .

 2 off, let us show that Σ and Σ are essentially the same. One could compute Σ = s (1 + Since γ, γ ∈ [0, π 2 ], tan γ 2 and tan γ 2 are smaller than 1, and since the diagrams are non-null, we get Σ = Σ by Proposition 6, which implies γ = γ . If γ = γ = 0, then Σ and Σ have full rank. Moreover, if γ = γ = π 2 , then Σ and Σ are not colinear to the identity. Hence, if γ = γ ∈]0, π 2 [, then we can apply Proposition 6. By Proposition 6, there exists d = e iϕ0 0 0 e iϕ1 such that U = U • d and V † = d † • V † .

2 with γ = π 2 -2 arctan 1 √ 5 and β 1 = arctan 2 .

 221512 It remains to prove that the two scalar diagrams are equal:

  It now remains to prove the rule (EU) can be derived. We decompose the left hand side diagram as such:

f

  : x → arctan tan (α 1 ) cos (x) + tan (α 3 ) cos (α 2 -x) 1 -tan (α 1 ) cos (x) tan (α 3 ) cos (α 2 -x) g : x → tan (α 1 ) cos (x) + tan (α 3 ) cos (α 2 -x) α 3 ) cos α 2 -π 2

  2 +γ 4 + (n+(-1) n m)γ 2 γ3+ (n+m)π


Charley Gros

Benjamin De Leener

Sara M Dupont

Allan R Martin

Michael G Fehlings

Rohit Bakshi

Subhash Tummala

Vincent Auclair

Donald G Mclaren

Virginie Callot

Julien Cohen-Adad

Michaël Sdika
email: michael.sdika@creatis.insa-lyon.fr

Automatic spinal cord localization, robust to MRI contrasts using global curve optimization

Keywords: Spinal cord, MRI, Detection, Segmentation, Global Optimization, Machine learning Spinal Cord PMJ: Pontomedullary junction FOV: Field of View (S-I) axis: Superior-Inferior axis HOG: Histogram of Oriented Gradient Degenerative Cervical Myelopathy NMO: Neuromyelitis Optica SCI: Spinal Cord Injury ALS: Amyotrophic Lateral Sclerosis SYR: Syringomyelia

come   L'archive ouverte pluridisciplinaire

Introduction

The spinal cord (SC) plays a key role in the central nervous system by ensuring the conduction of both motor and sensory signaling between the brain and the peripheral nervous systems. Although SC magnetic resonance imaging (MRI) has long been technically challenging, MRI has been increasingly used in the last two decades to provide valuable quantitative information through SC morphometry [START_REF] Fradet | Morphometrics of the entire human spinal cord and spinal canal measured from in vivo high-resolution anatomical magnetic resonance imaging[END_REF]Martin et al., 2017b;[START_REF] Papinutto | Age, gender and normalization covariates for spinal cord gray matter and total crosssectional areas at cervical and thoracic levels: A 2D phase sensitive inversion recovery imaging study[END_REF] and to evaluate SC damage a range of neurologic disorders such as multiple sclerosis, traumatic, and neurodegenerative diseases [START_REF] Cawley | Spinal cord atrophy as a primary outcome measure in phase II trials of progressive multiple sclerosis[END_REF][START_REF] Cohen-Adad | Involvement of spinal sensory pathway in ALS and specificity of cord atrophy to lower motor neuron degeneration[END_REF][START_REF] Cohen-Adad | Demyelination and degeneration in the injured human spinal cord detected with diffusion and magnetization transfer MRI[END_REF][START_REF] Ellingson | Morphology and morphometry of human chronic spinal cord injury using diffusion tensor imaging and fuzzy logic[END_REF][START_REF] Grabher | Neurodegeneration in the spinal ventral horn prior to motor impairment in cervical spondylotic myelopathy[END_REF]Martin et al., 2017a;[START_REF] Nakamura | Preferential spinal central gray matter involvement in neuromyelitis optica[END_REF]. With newer acquisition techniques that enhance the quality of SC MRI data [START_REF] Stroman | The current state-of-the-art of spinal cord imaging: Methods[END_REF], there is a need for improved image processing methods to propose robust automated quantitative analysis pipelines.

Localizing the SC on MRI scans is a key step for automating quantitative analysis pipelines such as SC [START_REF] De Leener | Robust, accurate and fast automatic segmentation of the spinal cord[END_REF][START_REF] Horsfield | Rapid semi-automatic segmentation of the spinal cord from magnetic resonance images: application in multiple sclerosis[END_REF] and grey matter [START_REF] Dupont | Fully-integrated framework for the segmentation and registration of the spinal cord white and gray matter[END_REF][START_REF] Prados | Spinal cord grey matter segmentation challenge[END_REF][START_REF] Prados | Fully automated grey and white matter spinal cord segmentation[END_REF] segmentations, template registration (De Leener et al., 2017b;[START_REF] Stroman | Spatial normalization, bulk motion correction and coregistration for functional magnetic resonance imaging of the human cervical spinal cord and brainstem[END_REF] and B0 susceptibility-related distortion correction [START_REF] Vannesjo | Slice-wise rst-order shimming of the human spinal cord at 7T[END_REF][START_REF] Topfer | A 24-channel shim array for the human spinal cord: Design, evaluation, and application[END_REF]. While localizing the SC might appear as a rudimentary computerized task, it is much more challenging to achieve it robustly and accurately across a broad range of SC shapes, craniocaudal vertebral length, pathologies, image field of view (FOV), image resolution and orientation, types of contrast and image artifacts (e.g. susceptibility, motion, chemical shift, ghosting, blurring, Gibbs). For instance, about 20 recently-published SC segmentation methods still require manual intervention for initialization and optimization [START_REF] De Leener | Segmentation of the human spinal cord[END_REF], especially to identify specific anatomic landmarks within the SC. Hence, SC segmentation methods could be made fully automatic if initialized with a robust automated SC centerline localization module.

Most of the existing automatic SC localization methods took advantage of the ellipsoid shape of the SC in cross-section and are based on the Hough transform with vesselness filtering [START_REF] De Leener | Robust, accurate and fast automatic segmentation of the spinal cord[END_REF], active contour [START_REF] Koh | Automatic segmentation of the spinal cord and the dural sac in lumbar MR images using gradient vector flow field[END_REF], or continuous max-flow with cross-sectional similarity prior [START_REF] Pezold | Automatic Segmentation of the Spinal Cord Using Continuous Max Flow with Cross-sectional Similarity Prior and Tubularity Features[END_REF]. Other methods are based on energy-minimization methods through image-based strategies, by pattern registration with partial volume effect modelization [START_REF] Carbonell-Caballero | Accurate quantification methods to evaluate cervical cord atrophy in multiple sclerosis patients[END_REF] or atlas registration with topology constraint [START_REF] Chen | Automatic magnetic resonance spinal cord segmentation with topology constraints for variable fields of view[END_REF]. While these automatic algorithms have shown good performance [START_REF] De Leener | Segmentation of the human spinal cord[END_REF], they often require a specific region of interest or are limited to specific contrast and resolution.

This paper introduces a novel fully automatic, accurate and fast algorithm (OptiC) to detect the center of the SC. OptiC is based on a supervised machine learning detection method combined with a distance-transform-based global optimization. Its performance was evaluated on non-processed data (i.e., not cropped, filtered, reoriented, intensity-normalized, etc.) across a large multi-center dataset from 20 sites (n=20) including a variety of image quality, resolution and pathologies across four image contrasts: T 2 -weighted (T 2 w, n=287), T 1weighted (T 1 w, n=120), T 2 * -weighted (T 2 * w, n=307) and diffusion weighted (DWI, n=90) 3D images. For the sake of transparency and reproducibility, we have made this method publicly available as part of the Spinal Cord Toolbox 1 (v3.0.2 and above) free and open-source software (De Leener et al., 2017a) and on the webpage of the CREATIS lab 2 .

The main contributions of this work are (i) an original framework to robustly find the SC centerline, based on machine learning and global curve optimization, (ii) an efficient algorithm to solve the optimization problem in a non-iterative way with linear complexity, (iii) a robust ponto-medullary junction (PMJ) detector to separate brain and spine regions and (iv) an extensive validation on a broad dataset involving several centers, pathologies, MR contrasts and image resolutions.

In the following sections, the theoretical aspects of OptiC will be detailed, along with the optimization of its hyperparameters and the evaluation across the testing dataset. Second, the results of the training and testing steps will be presented. Finally, the key contributions of OptiC method will be discussed.

1 sourceforge.net/projects/spinalcordtoolbox/ 2 creatis.insa-lyon.fr/site7/en/ctrDetect 2. Material and Methods

Theory

In the following subsections, the OptiC approach will be detailed and a new algorithm will be proposed to efficiently find a global minimum to the optimization problem. The final section presents an automatic post-processing method to differentiate brain and spinal cord regions in order to prevent analysis of the brain region.

OptiC SC centerline detection

The OptiC method relies on two prior anatomical sets of knowledge: (i) the SC has a specific tubular shape and (ii) the SC centerline is continuous along the superior-inferior (S-I) (cranio-caudal) axis. Consequently, the OptiC approach encodes these priors in 2 distinct steps (see Figure 1): (i) step 1: a SC localization map is computed, using standard object detection tools, (ii) step 2: an optimization problem is solved to take into account the global tubular shape of the SC. As a result, the SC centerline is modeled as a curve compromising the tubular similarity and the curve continuity. The goal of the first step is to compute a localization map of the SC, in a supervised manner, based on the detection of its specific ellipsoidal shape. It is achieved by predicting the presence of the SC in a given 2D image-patch with a classification model, as commonly done for object detection. Trained on labeled image-patches, the classification model assigns to each voxel a confidence of SC presence within the corresponding patch, resulting in a SC localization map.

Step2: Spinal cord centerline modeling

From the previously computed SC localization map, the goal of the second step is to model the SC as a curve which constrains the centerline continuity. Thus, the resulting centerline C is a regularized curve compromising the localization map S values and the SC continuity, which could be expressed as the solution of the following optimization problem (Eq. 1):

min c 0 ,... ,c n-1 C n-1 (c 0 , ... , c n-1 )=-∑ z=0 n-1 S z (c z )+ λ ∑ z =0 n-2 || c z+ 1 -c z| | ❑ 2 (Eq. 1)
with λ , a regularization coefficient; c z , the SC center and S z , the localization map on the z th slice; n , the number of slice in the 3D input MRI volume. Thus, the optimization problem has to make a trade-off between the classifier confidence in the SC localization, and the SC continuity between adjacent slices, respectively encoded by the first and second terms of (Eq 1).

Optimization problem resolution

The optimization problem (Eq. 1) could be solved using a gradient descent approach. This would be efficient but would not ensure finding the global solution. A brute force algorithm would find the global solution but would be impractical with a O( s n ) complexity, where s is the number of voxels per slice and n the number of slices in the 3D MRI volume).

In this section, a novel algorithm is detailed to efficiently solve (Eq. 1) and guarantee robustness of the optimization problem since it provides a global solution. Based on distancetransforms, this approach has a complexity linear with the number of voxels, in O( s .n) .

To solve the problem, an auxiliary sequence of 2D images ( M k ) is defined as:

M k (x k )=min x 0 , ..., x k-1 C k (x 0 , ... , x k ).
(Eq. 2) with k the slice number. Slices above the index k are removed, and the problem is seen as a function of the center on slice k only. If the M n-1 image is known, the relation min x n-1 M n-1 (x n-1 )=min x 0 ,..., x n-1 C n-1 ( x 0 , ... , x n-1 ) , (Eq. 3) enables finding the center solution on the last slice:

c n-1 =arg min x n -1 M n-1 (x n-1
) . However, computing the M n-1 image is as difficult as solving the original problem. Fortunately, the sequence ( M k ) satisfies the following recurrence relationship:

M k (x k )=min x 0 , ..., x k-1 -∑ z=0 k S z ( x z )+ λ ∑ z =0 k-1 || x z+1 -x z| | ❑ 2 ¿-S k ( x k )+ λ min x k-1 { M k -1 (x k-1 ) λ + || x k -x k-1 || ❑ 2 } ¿-S k ( x k )+ λ D( M k-1 λ )(x k )
with D defined by the image operator D(f )(x )=min y f ( y )+||x-y||❑ 2 . This recurrence relation not only enables us to practically compute the (M k ), it also provides a way to infer the full centerline, solution of (Eq. 1), from the top to the bottom. Indeed, given c k , c k-1 can be expressed as:

c k-1 =argmin x k-1 { M k-1 ( x k-1 ) λ + || c k -x k-1 || ❑ 2 } ¿ N ( M k-1 λ )(c k ) (Eq. 4)
where N is the proximity operator defined as N ( f )( x )=argmin y {f ( y)+||x-y||❑ 2 } .

D and N are important operators for convex optimization and image processing. They can be used to compute the distance transform and the nearest neighbors transform of an image. As proposed previously [START_REF] Felzenszwalb | Distance transforms of sampled functions[END_REF][START_REF] Meijster | A General Algorithm for Computing Distance Transforms in Linear Time[END_REF], these two transforms can be computed efficiently (with linear complexity) when the input function is uniformly sampled, i.e. is an image.

The resolution of (Eq. 1) is formalized in Algorithm 1. Fast distance and nearestneighbor transforms are first used to compute the sequences ( M k ) and (N k ) from the bottom to the top of the volume. On the top slice, the SC center is then found as the minimum of M n-1 . The complete centerline sequence

c (¿¿ k ) ¿
is then inferred by back-tracking on the nearest-neighbor-transforms ( N k ) . As the computation of M k and N k is linear, in O(s) , the complexity of OptiC is also linear, in O( s .n) .

The proof that the Algorithm 1 finds a global solution of the problem defined by (Eq. 1) is fairly simple. Using (Eq. 3), the minimal value of C n -1 , Ĉn-1 , is given by the global minimization of M n-1 . This global minimization is performed by looping on the pixel of the M n-1 2D image. The backtracking procedure given by (Eq. 4) ensures finding a path, a

c (¿¿ k ) ¿ sequence, satisfying C n -1 ( c 0 , ... , c n -1 )= Ĉn-
M z =-S z + λ D( M z-1 λ ) 5: N z =N ( M z-1 λ ) 6: c n-1 =argmin x n-1 M n -1 (x n-1 ) 7:
for z=n-2:0 do 8:

c z =N z+1 (c z+1 )
2.1.1.4. Implementation details In the optimization problem, the z coordinate was assumed to be the (S-I) axis for all input volumes. Thus, for each MRI volume, preprocessing steps involved imposing a common image orientation (RPI, i.e. Right-to-left, Posterior-to-anterior, Inferior-to-superior) based on the image header (NIFTI input).

In order to be robust to the native image resolution, the input image is re-interpolated to the resolution 1×1×p z mm using trilinear interpolation, p z being the native pixel size in z direction. To account for the image resolution in the (S-I) axis, the regularization coefficient λ is parameterized with a resolution independent parameter λ ' and the interslice distance p z :

λ= λ ' p z 2
The SC localization map was computed using a Support Vector Machine (SVM) trained with Histogram of Oriented Gradient features (HOG) as described previously [START_REF] Dalal | Histograms of oriented gradients for human detection[END_REF]. The HOG features were computed on patches of 32×32 size, from 2D axial slices, using the OpenCV library [START_REF] Bradski | Learning OpenCV: Computer Vision with the OpenCV Library[END_REF] with default parameters and signed gradients. Note that the use of a small patch size as well as the sophisticated normalization scheme of HOG makes the computation of the localization map robust to inhomogeneity artifacts and intensity scale changes.

Pontomedullary junction detection

Typical cervical SC scans include a part (or all) of the brain, which therefore needs to be excluded in the SC centerline output of the proposed algorithm. Here, an automatic method is introduced to detect the presence and position of the PMJ and subsequently crop the SC centerline above it. This method is included in OptiC as a post-processing module. The choice of the PMJ is driven by previous studies showing that this anatomical landmark could be used as a reference for SC template-based analyses [START_REF] Stroman | Spatial normalization, bulk motion correction and coregistration for functional magnetic resonance imaging of the human cervical spinal cord and brainstem[END_REF]. Moreover, this anatomical landmark can be precisely identified from multiple image contrasts and its neighboring voxels exhibit appropriate features for machine learning methods.

Figure 2 details the post-processing pipeline to detect the PMJ. The SC is first straightened along the extracted centerline, with a robust algorithm that preserves the SC topology and its internal and adjacent structure (De Leener et al., 2017b). Straightening of the SC is used to normalize the shape of the PMJ regardless of the patient position during the MRI scan. To detect the presence of the PMJ on the straightened volumes, a HOG+SVM 2D classifier (as described in section 2.1.1.4.) is trained on the midsagittal slice of the straightened training data, resampled to a 1×1mm resolution. When the PMJ is detected on the testing data, the centerline is cropped above it.

Materials

To make sure the algorithm performs well in a large variety of data (i.e., robust to MRI vendor, contrast, sequence parameters, artifacts, types of pathology, etc.), retrospective data from 20 research and clinical centers were collected. The dataset was composed of 501 subjects, including 173 patients with neurologic conditions [multiple sclerosis (MS) n=79, degenerative cervical myelopathy (DCM) n=63, neuromyelitis optica (NMO) n=19, traumatic spinal cord injury (SCI) n=5, amyotrophic lateral sclerosis (ALS) n=5, syringomyelia (SYR) n=2]. Volumes were acquired on 3T systems using standard coils, spanning a broad range of FOVs (e.g. cervical, thoracic, including or not brain sections) and four different contrasts: T 2 w (n=287), T 1 w (n=120), T 2 * w (n=307) and DWI (n=90). Spatial resolutions included isotropic (n=366, from 0.7 to 1.3mm) and anisotropic data with axial (n=430, in plane: from 0.3 to 1.5mm, slice thickness: from 1.0 to 24.5mm), or sagittal orientation (n=8, in plane: from 0.4 to 0.7mm, slice thickness: from 0.8 to 2.8mm). Demography and resolution are summarized in Figure 3. For each of the four contrasts, the dataset (N tot volumes) was randomly split into two independent subsets, the validation (N valid volumes) and the testing (N test volumes) datasets, as illustrated in Figure 4. The training-validation dataset (used to train and validate the model) involved 40 healthy controls (HC), and the testing dataset involved the rest of the HC, and patients. The reason for not including patient data for the training-validation was for providing an unbiased and generic model, which would then be tested on patients (ultimately, this algorithm could be applied for diagnosis purposes, hence without prior knowledge of the possible pathology). The center of the SC was manually localized using FSLView (Jenkinson et al., 2012) on each slice of each volume as ground truth for the SC centerline. Similarly, the PMJ position was manually localized on the midsagittal slice of each straightened T 1 w and T 2 w volume with brainstem section. The training-validation data were also split into 2 independent subsets: the first (n valid (1) ) for model training, the second (n valid (2) ) for evaluation.

Hyperparameters optimization

Three hyperparameters were investigated: the number of training volumes, the effect of rotation augmentation in training, and the regularization coefficient λ ' . To do so, for each hyperparameter, 100 bootstrap iterations were run on the training-validation dataset (N valid ), split into two independent data subsets: one for training (n valid (1) from which the training volumes were selected), the other for evaluation (n valid (2) , common for all the tested values of the current hyperparameter), as detailed in Figure 4. The best model was selected by computing the Mean Square Error (MSE) between the predicted Ĉ and the ground-truth C centerlines: MSE=√❑ with n is the number of 2D axial slices. For each hyperparameter, significant differences between the tested values were investigated with a Kruskal-Wallis H-test. After optimization, each hyperparameter was fixed for the remainder of the study.

Number of training images

In this part, we investigated if increasing the number of 3D images used to train the model yielded a more robust SC detection. Based on metrics measured on a validation dataset, models with the following number of training volumes were compared: 1, 5, 10, 15. The number of training 3D images (n train ) achieving the best results were retained for the rest of the study (i.e. n valid (1) =n train and n valid (2) =N valid -n train ).

Rotation training augmentation

Data augmentation is a common strategy to increase detection robustness in machine learning. This study investigated whether the model can be made rotation invariant or whether adding a rotation to the input volume improved the model generalization and prevented overfitting. Rotations around the (S-I) axis were applied to the training data (n train ) with different angular ranges (±30° and ±180°) and angular increments (each 5° and each 10°).

Regularization coefficient

The regularization coefficient λ ' is a hyperparameter that affects the (S-I) axis regularization of the extracted centerline. The higher λ ' is, the more the centerline points from two adjacent slices tend to converge in the X-Y plane. To optimize λ ' , discrete values were investigated from 0 to 4. Note that, when λ ' is null, the centerline is given by the maximal values of the SC localization map, and when λ ' goes to infinity, the centerline will tend to be a vertical straight line.

Evaluation

Once the OptiC hyperparameters were selected, OptiC performance was evaluated on the testing dataset (N test ), using both HC and patients (see Figure 4).

OptiC SC centerline detection

For each contrast, an OptiC optimized model was evaluated on the testing dataset and compared to a recently-published automatic SC localization method [START_REF] De Leener | Automatic Segmentation of the Spinal Cord and Spinal Canal Coupled With Vertebral Labeling[END_REF], based on the Hough transform detection with vesselness filtering (referred as Hough in the following paragraphs). Validation metrics were MSE between the predicted and ground truth SC centerlines and the localization rate, defined as the percentage of slices for which the predicted centerline was included in the manually-segmented SC.

Pontomedullary junction detection

To achieve fully automatic pipelines involving multi-center datasets with a variety of FOVs, the PMJ detection is a mandatory post-processing step to restrict the centerline curve to SC regions. The proposed method to detect the PMJ was evaluated on T 2 w and T 1 w datasets that included a large number of volumes with brain sections: PMJ appeared on 45% of T 2 w and 90% of T 1 w volumes of the testing dataset. The training and testing datasets were the same as in the previous section.

The detection performance was assessed with recall (measure of completeness) and precision (measure of exactness) metrics: recall= TP TP+ FN precision= TP TP+ FP where TP (true positives) and FN (false negatives) are the number of volumes with brain section where the PMJ is respectively detected and not detected; FP (false positives) are the number of volumes on which PMJ is detected even though the brain section is not present.

For volumes with detected PMJ, the localization accuracy is evaluated by computing the distance between the predicted and true PMJ position along the (S-I) axis.

Proof-of-concept: application to SC segmentation

SC segmentation provides relevant morphometric information about SC atrophy [START_REF] Yiannakas | Fully automated segmentation of the cervical cord from T1-weighted MRI using PropSeg: Application to multiple sclerosis[END_REF] and is an important step in template-based analysis pipelines (De Leener et al., 2017a). Thus, as a proof-of-concept application, a SC segmentation method PropSeg [START_REF] De Leener | Automatic Segmentation of the Spinal Cord and Spinal Canal Coupled With Vertebral Labeling[END_REF], which uses propagation of 3D meshes, was initialized using an OptiC centerline to compare its benefit over the default initialization. Dice coefficient [START_REF] Dice | Measures of the Amount of Ecologic Association Between Species[END_REF], between the predicted and ground-truth segmentations (done by experienced users), evaluated the quality of the segmentation of the testing subjects. The computation time was also measured for the SC segmentation process when initialized with OptiC versus with the default version.

Results

Figure 5 presents the results of the hyperparameter optimization for each contrast. The selected hyperparameters were chosen through a bootstrap analysis, as a trade-off between detection accuracy (measured by the MSE) and simplicity of the model. The OptiC hyperparameter (filled in green) was the same for all contrast, selected as either the best tested value or not significantly different from the best tested value, according to Kruskal-Wallis H-tests.

OptiC SC centerline detection

Figure 6 shows qualitative examples of SC centerline detection on healthy controls (HC) and patient volumes from the testing dataset. On both axial and sagittal views, SC centerlines were predicted closely to the ground-truth centerlines. OptiC performance was evaluated on testing datasets containing four contrasts, including HC and a large number of patients. OptiC was compared to the Hough method in Figure 7 in terms of localization rate and MSE. The localization rates were 98.77 ± 4.49% for OptiC and 63.75 ± 48.79% for Hough (averaged across contrasts n=644, Wilcoxon signedrank test p-value<0.01). The MSE results, for subjects with a detected cord (i.e. localization rate > 0%), were 1.02 ± 1.10mm for OptiC and 9.95 ± 13.76mm for Hough (averaged across contrasts n=485, Wilcoxon signed-rank test p-value<0.01). Note that the Hough method performed adequately on HC but failed on most patients, whereas OptiC performed similarly well on both HC and patients. For each metric, the best value is indicated in green in the y-axis. Results are compared between the Hough method [START_REF] De Leener | Automatic Segmentation of the Spinal Cord and Spinal Canal Coupled With Vertebral Labeling[END_REF] (left distribution plots) and the OptiC method (right distribution plots), where each point represents a healthy control (blue) or patient (red). Each method was evaluated on the 4 contrast datasets: T 2 w, T 1 w, T 2 * w and DWI (from left to right). The table presents the average metrics per contrast for each method, and for the HC and the patient (Pat.) subgroups.

Pontomedullary junction detection

Table 1 shows the results of the PMJ automatic detection on the T 2 w and T 1 w testing datasets. The proposed detection method provides good results in terms of recall (T 2 w: 84%; T 1 w: 100%), with few FN detections, as well as high precision (T 2 w: 100%; T 1 w: 99%), with few FP detection. Moreover, the averaged signed distance error between the predicted and the true PMJ position along the (S-I) axis is in the range of 1 cm below the PMJ (T 2 w: -11.41 ± 5.71mm; T 1 w: -5.75 ± 10.29), resulting in an adequate separation of brain and SC structures around the medulla oblongata (situated just caudal to the PMJ).

Recall [%]

(Best value: 100) Recall and precision evaluated the ability of the proposed method to automatically detect the presence of the PMJ in the MRI volumes, referred to as True Positive volumes (TP). The signed distance error measures the distance between the predicted and the true PMJ position along the (S-I) axis.

Proof of concept: application to SC segmentation

Automatic SC segmentation tasks remain a challenging problem on a variety of MRI data [START_REF] De Leener | Segmentation of the human spinal cord[END_REF], mainly due to the difficulty to robustly detect the SC. This localization is often the first step of automatic SC segmentation pipelines. Hence, as a proof of concept of OptiC's application, we used OptiC to initialize a publicly available SC segmentation method, PropSeg [START_REF] De Leener | Automatic Segmentation of the Spinal Cord and Spinal Canal Coupled With Vertebral Labeling[END_REF]. Figure 8 compared the results of PropSeg when initialized with the default method (called "Hough+PropSeg") and with the SC centerline provided by OptiC (called "OptiC+PropSeg"). Dice coefficients calculated between the predicted and the ground-truth segmentations were 0.87 ± 0.14 when using OptiC initialization versus 0.44 ± 0.43 when using the default initialization (averaged across contrasts n=644, Wilcoxon signed-rank test p-value<0.01). OptiC was particularly valuable in pathological cases with a Dice coefficient of 0.84 ± 0.14 versus 0.36 ± 0.41 when using the default initialization (averaged across contrasts n=259, Wilcoxon signed-rank test p-value<0.01).

The "OptiC+PropSeg" results for T 2 w and T 1 w datasets were obtained by using the centerline automatically cropped at the PMJ position. This additional post-processing slightly improved the Dice coefficient of segmentations on MRI volumes with PMJ sections: Dice coefficients were 0.89 ± 0.13 without cropping versus 0.90 ± 0.11 with PMJ cropping (averaged across contrasts n=183). Moreover, the above results were obtained efficiently when using OptiC, 0.41 ± 0.46 seconds per 2D axial slice (including reading/writing tasks, on an iMac i7 4-cores 3.4 GHz 8Gb RAM) versus 8.05 ± 82.90 seconds without OptiC initialization (averaged across contrasts n=644, Wilcoxon signed-rank test p-value<0.01). Note that, OptiC centerline extraction is not time consuming since it is achieved in 10.00±6.98 milliseconds per 2D axial slice (averaged across contrasts n=644). 2015) segmentations when centerline initialization is performed with the default method (left distribution plots) vs. the OptiC method (right distribution plots), where each point represents a healthy control (blue) or patient (red). Each method was evaluated on the 4 contrast datasets: T 2 w, T 1 w, T 2 * w and DWI (from left to right). The table presents the averaged metrics per contrast for each method.

Discussion

In this article, a new automatic and fast method to detect the center of the SC on MRI data was presented. The proposed method is based on a supervised machine learning detection method combined with a novel algorithm that efficiently provides the global solution of an optimization problem. Additionally, a new post-processing feature to detect and correct for the PMJ landmark allows the pipeline to handle datasets with a variety of image FOVs.

Key benefits of OptiC are (i) robustness to image artifacts (e.g. susceptibility, motion, ghosting, chemical shift, blurring, Gibbs) thanks to the sophisticated normalization HOG strategy [START_REF] Dalal | Histograms of oriented gradients for human detection[END_REF], (ii) respect for anatomical cord continuity as a result of the regularization along the (S-I) axis, (iii) robustness of curve extraction ensured by the finding of the global minimum (iv) fast computation time by reason of a linear complexity with the number of voxels. Further efforts to extend our results might consider the use of more sophisticated object detection classifier, such as Convolutional Neural Networks [START_REF] Long | Fully convolutional networks for semantic segmentation[END_REF][START_REF] Ronneberger | U-Net: Convolutional Networks for Biomedical Image Segmentation[END_REF], and/or other engineering features, such as SURF features (Bay et al., 2008/6).

The hyperparameter study shows that training the model with only one image provides good results. It suggests the possibility to fastly and easily create a new model if needed (e.g. adapted to a new contrast or to a specific pathology) without requiring a large amount of training 3D images. Rotation augmentation on the training data did not prove useful to detect the SC, even if HOG features are not inherently rotation invariant. This can probably be explained by the setting of a common image orientation in the preprocessing steps. The regularization coefficient optimization suggests that the solving of the optimization problem ( λ ' ≠ 0 ) provides significant improvements compared to the SC localization map output by the classifier ( λ ' =0 ).

This study evaluated the performance of the OptiC method to automatically detect the SC on a variety of MRI data and to provide a reliable SC centerline. OptiC was validated on a large multi-center dataset with heterogeneous contrast, resolution, FOVs. Both accuracy and robustness of the method have been shown, even on challenging pathological cases with hyperintense lesions or cord deformations, such as in 3D images from patients with MS (n=79) or DCM (n=63): the averaged MSE from patients was 1.08mm vs. 0.97mm from controls. The robustness of the proposed method provides a proof-of-concept to allow more robust automatic analysis on pathological cases with clinical applications.

The automatic detection of the PMJ allows the extraction of a SC centerline restricted to SC regions even on MRI volumes that include the brain. Other anatomical structures such as the top of C1 vertebra or the bottom of the medulla oblongata might also have been appropriate to separate the brain and SC regions, but these structures are much more challenging to label consistently and to detect automatically. The straightening of the SC is a necessary step to robustly detect the PMJ, however, the computation of the warping fields is time consuming. Validation results demonstrate the precision of the PMJ detection method (low rate of FP detections) and its accuracy (<1cm below the PMJ on average). Note that, the suggested separation position is situated on average within the medulla oblongata, below the other (more superior) brain regions. Notably, its benefit for SC segmentation tasks was shown by slightly improving the Dice coefficient of the segmentation. Therefore, delineating the SC below the medulla oblongata would have benefits when integrated on fully automatic pipelines.

This paper shows the benefits of the OptiC centerline to initialize a state-of-the-art automatic SC segmentation technique [START_REF] De Leener | Automatic Segmentation of the Spinal Cord and Spinal Canal Coupled With Vertebral Labeling[END_REF]. Automatic detection of the SC has also recently been proposed to improve B0 field correction in targeted optimization of shim currents [START_REF] Vannesjo | Slice-wise rst-order shimming of the human spinal cord at 7T[END_REF][START_REF] Topfer | A 24-channel shim array for the human spinal cord: Design, evaluation, and application[END_REF] to adapt the shimming to the exact position of the SC in T 2 *-w images. Other applications, such as template registration or shape analysis, could also benefit from this method.

Conclusion

SC centerline detection is a key step to initialize fully automatic pipelines in various applications such as SC segmentation, shape analysis, template registration or B0 field correction. In this paper, a fully automated method for localizing the SC from MRI data was presented, extensively validated on a large multi-center heterogeneously acquired datasets with a large variety of pathologies and successfully compared with a state-of-the-art SC localization technique. The proposed SC centerline detection method is freely available as part of the Spinal Cord Toolbox: http://sourceforge.net/projects/spinalcordtoolbox/.

Figure 1 .

 1 Figure 1. OptiC pipeline overview. OptiC is based on two steps: (1) The spinal cord (SC) probabilistic localization map is computed from the input image, (2) The SC centerline is modeled as a regularized curve, constrained by a trade-off between the localization map values and the SC continuity in the superior-to-inferior (S-I) axis.

1

 1 and thus it enables finding a global solution of the optimization problem.

Figure 2 .

 2 Figure 2. Automatic post-processing pipeline to restrict the OptiC centerline below the ponto-medullary junction (PMJ). The spinal cord (SC) is first straightened along the regularized SC centerline (blue). A classifier localization map (heat probabilistic map) allows detection of the PMJ on the straightened image (middle). Finally, the SC centerline is restricted below the detected PMJ (right).

Figure 3 .

 3 Figure 3. Demography and resolution. Retrospective data from a large variety of pathologies (1) and contrasts (2) were used. For clarity purposes, the distribution of subjects is only shown for pathologies and contrasts. Volumes were acquired at 3T, spanning a large range of slice orientations, in plane resolutions (3a) and slice thickness (3b). Abbreviations: MS: multiple sclerosis; DCM: degenerative cervical myelopathy; NMO: neuromyelitis optica; SCI: Traumatic spinal cord injury; ALS: amyotrophic lateral sclerosis; SYR: syringomyelia.

Figure 4 .

 4 Figure 4. Training-validation and testing dataset breakdown. The dataset of each contrast (N tot ) was randomly split into two independent subsets: the training-validation (N valid =40) and the testing (N test =N tot -N valid ) datasets. The training-validation dataset only involved healthy controls (HC, in blue) whereas the testing dataset involved both patients (in purple) and HC.The training-validation data were also split into 2 independent subsets: the first (n valid(1) ) for model training, the second (n valid (2) ) for evaluation.

Figure 5 .

 5 Figure 5. Hyperparameter Optimization. Three hyperparameters were investigated (from left to right): the number of training 3D images, the effect of rotation augmentation in training and the regularization coefficient. The optimization was done for the four contrasts (from top to bottom): T 2 w, T 1 w, T 2 * w, DWI. Here, Avg. ± Std. of the Mean Square Error (MSE) between the estimated centerline and the ground truth across 100 bootstrap iterations are plotted in red. For each hyperparameter optimization, (*) indicates significant differences between the selected value (filled in green) and the other tested values, according to Kruskal-Wallis H-tests (*: p≤0.05, **: p≤0.01).

Figure 6 .

 6 Figure 6. Examples of automatic spinal cord (SC) centerline detection on healthy controls (HC) and patients, on T 2 w (top-left), T 2 * w (top-right), T 1 w (bottom-left) and DWI (bottomright) images. This include a comparison between ground-truth (red) and automatic (blue) SC centerlines on axial and sagittal views. MS: multiple sclerosis; DCM: degenerative cervical myelopathy; NMO: neuromyelitis optica; SCI: traumatic spinal cord injury; ALS: amyotrophic lateral sclerosis.

Figure 7 .

 7 Figure 7. Results of the spinal cord (SC) localization in terms of localization rate (top) and Mean Square Error (MSE) between the ground-truth and predicted SC centerline (bottom).For each metric, the best value is indicated in green in the y-axis. Results are compared between the Hough method[START_REF] De Leener | Automatic Segmentation of the Spinal Cord and Spinal Canal Coupled With Vertebral Labeling[END_REF] (left distribution plots) and the OptiC method (right distribution plots), where each point represents a healthy control (blue) or patient (red). Each method was evaluated on the 4 contrast datasets: T 2 w, T 1 w, T 2 * w and DWI (from left to right). The table presents the average metrics per contrast for each method, and for the HC and the patient (Pat.) subgroups.

Figure 8 .

 8 Figure 8. Results of the spinal cord (SC) segmentation in terms of Dice coefficients between the ground-truth and predicted SC segmentation. For each metric, the best value is indicated on green in the y-axis. Results are a comparison of PropSeg (De Leener et al.,2015) segmentations when centerline initialization is performed with the default method (left distribution plots) vs. the OptiC method (right distribution plots), where each point represents a healthy control (blue) or patient (red). Each method was evaluated on the 4 contrast datasets: T 2 w, T 1 w, T 2 * w and DWI (from left to right). The table presents the averaged metrics per contrast for each method.

Table 1 .

 1 Pontomedullary junction (PMJ) detection results on T 2 w and T 1 w testing datasets.

	Precision [%] (Best value: 100)	Distance error [mm] (Best value: 0.00)

Acknowledgements:

-We would like to acknowledge all members from the NeuroPoly laboratory for useful and constructive discussions.

-Acknowledged for sharing data: Manuel Taso (University of Aix-Marseille); Julius Dewald and Haleh Karbasforoushan (Northwestern University); Seth Smith, Bailey Lyttle, Ben Conrad and Bennett Landman (Vanderbilt University); Barry Bedell (Biospective Inc.); Seok Woo Kim and Jisun Song (Hallym University); Emil Ljungberg (University of British Columbia); Eric Klawiter (Massachusetts General Hospital); Maryam Seif and Patrick Freund (University Hospital Zurich); Pierre-François Pradat and Habib Benali (Pitié-Salpêtrière Hospital); Pierre Rainville (Canadian Institutes of Health Research); Rohit Bakshi, Subhash Tummala and Shahamat Tauhid (Brigham and Women's Hospital, Harvard Medical School); Gilles Edan, Anne Kerbrat and Elise Bannier (University Hospital of Rennes); Pierre Labauge (University Hospital of Montpellier); Ferran Prados and Claudia Gandini Wheeler-Kingshott (University College London); Yaou Liu and Zhuoqiong Ren (Xuanwu Hospital); Tom Lillicrap (University of Newcastle). -This work was supported by: Canada Research Chair in Quantitative Magnetic Resonance Imaging (JCA), the Canadian Institute of Health Research [CIHR FDN-143263], the Canada Foundation for Innovation [32454], the Fonds de Recherche du Québec -Santé [28826], the Fonds de Recherche du Québec -Nature et Technologies [2015-PR-182754, 2017-200599], the Natural Sciences and Engineering Research Council of Canada [435897-2013], the Quebec BioImaging Network, the Program "Investissements d'Avenir" operated by the French National Research Agency [ANR-11-INBS-0006], the LABEX PRIMES [ANR-11-LABX-0063] of Université de Lyon, "Investissements d'Avenir" [ANR-11-IDEX-0007], the National Multiple Sclerosis Society NMSS RG-1501-02840 (SAS), the Centre National de la Recherche Scientifique (CNRS), the "Fondation A*midex-Investissements d'Avenir" and the "Fondation Aix-MArseille Université", NIH/NINDS R21 NS087465-01 (SAS), NIH/NEI R01 EY023240 (SAS), DoD W81XWH-13-0073 (SAS), Grant MOP-13034.


Jean-Claude Bermond
email: jean-claude.bermond@inria.fr

Takako Kodate
email: kodate@lab.twcu.ac.jp

Joseph Yu
email: joseph.yu@ufv.ca

Gossiping with interference in radio ring networks

come

Introduction

In this paper, we study the problem of gossiping with interference constraint in radio ring networks. Gossiping (or total exchange information) is a protocol where each node in the network has a message and wants to distribute its own message to every other node in the network. The gossiping problem consists in finding the minimum running time (makespan) of a gossiping protocol and efficient algorithms that attain this makespan.

Transmission model A radio network consists of communication devices equipped with an half duplex interface. The network is assumed to be synchronous and the time is slotted into rounds. The half-duplex hypothesis implies that a node can transmit or receive at most one message during a round. The network is modeled as a digraph, where the vertices represent the nodes and the arcs represent the possible communications. The messages are transmitted through the communication over the arcs and we will denote a call such a transmission.

Interference model

We use a binary asymmetric model of interference based on the distance in the communication digraph like the ones used in [START_REF] Bermond | Optimal gathering protocols on paths under interference constraints[END_REF][START_REF] Bermond | Hardness and approximation of gathering in static radio networks[END_REF][START_REF] Klasing | On the complexity of bandwidth allocation in radio networks[END_REF]. Let d(u, v) denote the distance, that is the length of a shortest directed path, from u to v in G and d I be a non negative integer. We assume that when a node u transmits, all nodes v such that d(u, v) ≤ d I are subject to the interference from u's transmission. So two calls (u, v) and (u , v ) do not interfere if d(u, v ) > d I and d(u , v) > d I .

This problem has been studied in [START_REF] Gasieniec | Gossiping with unit messages in known radio networks[END_REF] where approximation results are given (see also the survey [START_REF] Gasieniec | On efficient gossiping in radio networks[END_REF]).

Here we focus on the case where the transmisson network is a ring network C n on n nodes with the interference distance d I = 1. We presented some partial results at JCDCG^G 2013, and we have now solved completely the gossiping problem by giving the minimum running time (makespan). We show lower bounds and then give gossiping algorithms which meet these lower bounds and so are optimal.

It appears that the determination of the minimum gossiping time needs various tools. The optimal time depends on the congruence of n modulo 12. Our results are summarized in the following theorem.

Theorem 1 The minimum number of rounds R needed to achieve a gossiping in a ring network C n (n ≥ 3), with the interference model d I = 1 is :

      
2n -3 if n ≡ 0(mod 12) 2n -2 if n ≡ 4, 8(mod 12) 2n -1 if n is odd, except when n = 3 for which R = 3, and n = 5 for which R = 10 2n + 1 if n ≡ 2, 6, 10(mod 12) except for n = 6 for which R = 12

* Funded by ANR STINT