Guillaume Ch Èze
email: guillaume.cheze@math.univ-toulouse.fr

BOUNDING THE NUMBER OF REMARKABLE VALUES VIA JOUANOLOU'S THEOREM

In this article we bound the number of remarkable values of a polynomial vector field. The proof is short and based on Jouanolou's theorem about rational first integrals of planar polynomial derivations. Our bound is given in term of the size of a Newton polygon associated to the vector field. We prove that this bound is almost reached.

Introduction

In this paper we study polynomial differential systems in C 2 of this form:

dX dt = A(X, Y ), dY dt = B(X, Y ),
where A, B ∈ C[X, Y ] are coprime and deg(A), deg(B) ≤ k.

We associate to these polynomial differential systems the polynomial derivations

D = A(X, Y )∂ X + B(X, Y )∂ Y .
The computation and the study of rational first integrals of such polynomial differential systems is an old and classical problem. We recall that a rational first integral is a rational function f /g ∈ C(X, Y ) such that the curves λfµg = 0, where (λ : µ) ∈ P 1 (C), give orbits of the differential system. Thus it is a function f /g ∈ C(X, Y ) such that D(f /g) = 0.

When we study rational first integrals, we can always consider a rational first integral f /g with a minimum degree. We recall that the degree of f /g is equal to the maximum of deg(f ) and deg(g). Rational first integrals with minimum degree are indecomposable rational functions. We recall that a rational function is decomposable when it can be written u(h) where u ∈ C(T ), h ∈ C(X, Y ) and deg(u) ≥ 2, otherwise f /g is said to be indecomposable. Sometimes indecomposable rational functions are called minimal rational functions. In 1891, Poincaré has shown that if we consider an indecomposable rational first integral f /g, then the number of (λ : µ) ∈ P 1 (C) such that λfµg is reducible in C[X, Y ] or deg(λfµg) < deg(f /g), is finite. Poincaré called these kind of values (λ : µ) : remarkable values. Poincaré was interested by the intersection of different level sets of a given rational first integral. He has also shown how to use remarkable values in order to study the inverse integrating factor, see [START_REF] Poincaré | Sur l'intégration algébrique des équations différentielles du premier ordre et du premier degré[END_REF]. Recently, new results have been given in this direction, see [START_REF] Chavarriga | Darboux integrability and the inverse integrating factor[END_REF][START_REF] Ferragut | On the remarkable values of the rational first integrals of polynomial vector fields[END_REF][START_REF] Ferragut | Polynomial inverse integrating factors for polynomial vector fields[END_REF][START_REF] Coll | Polynomial inverse integrating factors for quadratic differential systems[END_REF][START_REF] Ferragut | Some new results on Darboux integrable differential systems[END_REF]. In [START_REF] García | On the polynomial differential systems having polynomial first integrals[END_REF][START_REF] García | Planar polynomial vector fields having a polynomial first integral can be obtained from linear systems[END_REF], remarkable values are used to study systems with polynomial first integrals. In [START_REF] Jesús | Local behavior of planar analytic vector fields via integrability[END_REF] they are used to study degenerate singular points. They also plays a role in an algorithm computing the decomposition of multivariate rational functions, see [START_REF] Chèze | Nearly optimal algorithms for the decomposition of multivariate rational functions and the extended Lüroth theorem[END_REF], and in another algorithm computing rational first integrals with bounded degree, see [START_REF] Bostan | Efficient algorithms for computing rational first integrals and Darboux polynomials of polynomial vector fields[END_REF].

Poincaré has shown, when all the singular points of the polynomial vector field are distinct, that the number of remarkable values is bounded by the number of saddle points plus 2, see [START_REF] Poincaré | Sur l'intégration algébrique des équations différentielles du premier ordre et du premier degré[END_REF]. Then a direct application of Bezout's theorem shows that the number of remarkable values is bounded by k 2 + 2.

Since Poincaré, a lot of authors, see below, have given bounds about the number of remarkable values and have studied the problem of reducibility in a pencil of algebraic curves. To the best of our knowledge all these bounds are given in term of the degree d of the rational function f /g. Even if there exists a relation between k and d, see below, this means that a geometric point of view about the pencil λfµg is generally used. In this note, we consider the pencil λfµg as level sets of a rational first integrals. With this point of view we can use results about rational first integrals and get easily a new bound about the number of remarkable values. Furthermore, our strategy allows us to give a bound on the total order of reducibility. This number is defined in the following way: We denote by f h and g h the homogeneous polynomial of degree deg(f /g) in C[X, Y, Z] associated to f and g and we set:

σ(f, g) = {(λ : µ) ∈ P 1 (C) | λf h -µg h is reducible in C[X, Y, Z]}.
This set is called the spectrum of f /g. We have introduced the homogeneous polynomial f h and g h in order to have a uniform definition for the remarkable values. The situation where deg(λfµg) is smaller than deg(f /g) corresponds to the situation where Z divides λf hµg h . Thus, by definition of the spectrum, the spectrum is the set of all remarkable values.

If (λ : µ) ∈ σ(f, g) then we have

λf h -µg h = n(λ:µ) i=1 f h (λ:µ),i e (λ:µ),i ,
where f h (λ:µ),i is the homogeneous polynomial in C[X, Y, Z] associated to the irreducible factor f (λ:µ),i of λfµg. The total order of reducibility is

ρ(f, g) = (λ:µ)∈σ(f,g) n(λ : µ) -1 .
For example, if f = Y and g = X 2 then σ(f, g) = {(0 : 1), (1 : 0)}. Indeed, X 2 is reducible and 1 = deg(f ) < deg(f /g) = 2. Moreover, f h = Y.Z, then n(1 : 0) = 2 and g h = x 2 then n(0 : 1) = 1. Thus ρ(f, g) = (2 -1) + (1 -1) = 1 and with this example we have: ρ(f, g) < |σ(f, g)|.

We get this inequality because X 2 is reducible and has not distinct irreducible factors. Thus if we want to bound the number of remarkable values with ρ(f, g) we have to also add the number of (λ : µ) ∈ P 1 (C) such that λf fµg h is a pure power. Then we introduce:

γ(f, g) = {(λ : µ) ∈ P 1 (C) | λf h -µg h = P e where e > 1 and P ∈ C[X, Y, Z]}.
Furthemore the number of elements in γ(f, g) is smaller or equal to 3, see [START_REF] Jouanolou | Équations de Pfaff algébriques[END_REF][START_REF] Busé | On the total order of reducibility of a pencil of algebraic plane curves[END_REF][START_REF] Shreeram | Translates of polynomials[END_REF] for a proof and [START_REF] Shreeram | Translates of polynomials[END_REF]Remark 8] for an example where the bound is reached.

Remark 1. We have:

|σ(f, g)| ≤ ρ(f, g) + |γ(f, g)| ≤ ρ(f, g) + 3.
Thus we are going to bound ρ(f, g) and it will give a bound on the number of remarkable values.

As our theorem uses Newton polygons, we recall that the Newton polygon of a Laurent polynomial f (X, Y ) = α,β c α,β X α Y β , where (α, β) ∈ Z 2 , is the convex hull in R 2 of the exponents (α, β) of all nonzero terms of f .

Theorem 2. Let D = A(X, Y )∂ X + B(X, Y )∂ Y be
points in N D ∩ N 2 , then: ρ(f, g) < B + 2.
Theorem 2 is well suited for "sparse" derivations i.e. when some coefficients of A and B are equal to zero. For example when we consider polynomials A and B of this form: c e,e X e Y e + c e-1,e X e-1 Y e + c e,e-1 X e Y e-1 + c 0,0 , where e is an integer strictly bigger than 1, then B = 3e + 2, and k = deg(A) = deg(B) = 2e. Thus for such examples Theorem 2 gives a linear bound in e. Then the bound is linear in the degree of the derivation. The Newton polygon associated to these polynomials A and B is given in Figure 1. The Newton polygon N D associated to this example is given in Figure 2. In these figures, small points correspond to coefficients equal to zero.

N c 3,3 X 3 Y 3 + c 2,3 X 2 Y 3 + c 3,2 X 3 Y 2 + c 0,0 0 - 6 Y
X q q q q q p q p q p p p p q q p p q q p p p p If we consider dense polynomials A, B with degree k, that is to say each coefficient of A and B is nonzero, then B = k(k + 1)/2. We get then a quadratic bound (in term of the degree k of the derivation) on ρ(f, g).

As the number of remarkable values is |σ(f, g)| then by Remark 1 we get:

Corollary 3. With the previous notations, we have: The number of remarkable values is smaller than

B + 2 + |γ(f, g)|, more precisely |σ(f, g)| < B + 2 + |γ(f, g)| ≤ B + 5.
This bound is almost reached. Indeed, consider:

D = (X k -1)∂ X -(kX k-1 Y + 1)∂ Y An indecomposable polynomial first integral of D is f (X, Y ) = Y (X k -1)+X.
The remarkable values associated to D are (1 : ω) where ω k = 1 and (0 : 1). Indeed, here the denominator g is equal to 1, thus deg(0f -1g) = deg(1) < deg(f /g). We deduce that in this situation k + 1 ≤ |σ(f, 1)|. Furthermore, for this derivation we have B = k, then by Remark 1, this gives

B + 1 ≤ |σ(f, 1)| ≤ ρ(f, 1) + |γ(f, 1)|.
Here |γ(f, 1)| = 1, because 0f h -1g h = Z k+1 . The point (0 : 1) is the only possible point in γ(f, 1). Indeed, if (λ : µ) = (0 : 1) and (λ : µ) ∈ γ(f, 1) then by definition we have λfµ = P e . As this implies f decomposable, we deduce that |γ(f, 1)| = 1. Then Corollary 3 gives

B + 1 ≤ |σ(f, 1)| < B + 3.
Thus the bound is almost reached. We can also remarked that in this example we have ρ(f, 1) = B, and Theorem 2 give ρ(f, 1) ≤ B + 1.

At last, we remark that we can always consider a rational function f /g as a rational first integral. Indeed, f /g is a rational first integral of the jacobian derivative

D f /g = A f /g ∂ X -B f /g ∂ Y ,
where

A f /g = ∂ Y (f )g -f ∂ Y (g) gcd ∂ Y (f )g -f ∂ Y (g), ∂ X (f )g -f ∂ X (g) and B f /g = ∂ X (f )g -f ∂ X (g) gcd ∂ Y (f )g -f ∂ Y (g), ∂ X (f )g -f ∂ X (g)
.

Thus if we want to study the reducibility in the pencil λfµg then we can always suppose that f /g is an indecomposable rational first integral of a derivation.

Related results. The study of remarkable values corresponds to the study of the irreducibility in a pencil of algebraic plane curves. It is an old problem and it has been widely studied since Bertini, see [START_REF] Steven | Bertini and his two fundamental theorems[END_REF]. It seems that Bertini and Poincaré has proved independently the following result: Theorem 4. The spectrum associated to an indecomposable rational function is finite.

To the best of our knowledge, the first author, after Poincaré, who has given a bound on the spectrum was Ruppert, in [START_REF] Ruppert | Reduzibilität Ebener Kurven[END_REF]. He has shown that the number of remarkable values associated to an indecomposable rational function of degree d is smaller than d 2 -1. Ruppert's strategy is based on the computation of the first de Rham's cohomology group of the complementary of a plane curve. This approach gives an effective method to compute the spectrum. Lorenzini has also studied the spectrum of a rational function in [START_REF] Lorenzini | Reducibility of polynomials in two variables[END_REF] and he has shown that ρ(f, g) ≤ d 2 -1. In [START_REF] Busé | On the total order of reducibility of a pencil of algebraic plane curves[END_REF], the authors have used Ruppert's approach and have shown that the bound d 2 -1 still works if we take into account the multiplicities of the factors. Theorem 2 can be seen as a counterpart of this bound with multiplicities. Indeed, in [START_REF] Poincaré | Sur l'intégration algébrique des équations différentielles du premier ordre et du premier degré[END_REF], Poincaré has introduced the remarkable factor

R(X, Y ) = (λ:µ)∈σ(f,g) n(λ:µ) i=1 f e (λ:µ),i -1 (λ:µ),i
and has also given the following relation:

deg(f ) + deg(g) -1 = k + deg(R).
For a proof of this equality, we can read [START_REF] Ferragut | On the remarkable values of the rational first integrals of polynomial vector fields[END_REF]. Thus a bound in term of the degree of the rational function which takes into account the multiplicities of the irreducible factors is in the same vein than a bound in term of the degree k of the vector field. The contribution of this present paper is to give a short proof which gives a nearly optimal bound in the sparse case. Stein, in [START_REF] Stein | The total reducibility order of a polynomial in two variables[END_REF], considers the polynomial case. He has shown that ρ(f, 1) ≤ d. It seems that the word "spectrum" has been introduced by Stein. This expression is also used by other authors. For example this expression is used in [START_REF] Busé | On the total order of reducibility of a pencil of algebraic plane curves[END_REF] because in this paper the set of remarkable values corresponds to the spectrum of a pencil of matrices. The strategy used by Stein was the following: first construct a rational first integral with some factors f (1:µ),i and then get a contradiction. The construction of the first integral was obtained with geometrical arguments specific to the polynomial case (i.e. g = 1). This approach has been extended to the rational case by Bodin in [START_REF] Bodin | Reducibility of rational functions in several variables[END_REF]. The bound obtained is ρ(f, g) ≤ d 2 + d. In this note, we are going to use the same approach. However, the construction of the rational first integral will be a direct consequence of a theorem due to Jouanolou. This allows to get a bound in term of k, the degree of the derivation, and not in term of d, the degree of the rational function. Furthermore, our proof is direct and gives a nearly optimal result for sparse derivations. Some authors have already used the Darboux theory of integrability in order to show that the spectrum is finite. Moulin-Ollagnier, in [START_REF] Moulin | Algebraic closure of a rational function[END_REF], has proved the finiteness of the spectrum by studying the number of distinct cofactors of a derivation. Moulin-Ollagnier calls the factor f (λ:µ),i "small Darboux polynomials". We can also mention the paper [START_REF] Chavarriga | Darboux integrability and the inverse integrating factor[END_REF], where the authors show with a simple proof using Darboux theory of integrability that the spectrum is finite. Unfortunately, these approachs do not give sharp bounds.

There exist many other papers about the spectrum, for more details read e.g. [START_REF] Bodin | Indecomposable polynomials and their spectrum[END_REF][START_REF] Bodin | Irreducibility of hypersurfaces[END_REF][START_REF] Najib | Sur le spectre d'un polynôme à plusieurs variables[END_REF][START_REF] Najib | Une généralisation de l'inégalité de Stein-Lorenzini[END_REF][START_REF] Busé | Noether forms for the study of noncomposite rational functions and their spectrum[END_REF][START_REF] Shreeram | Translates of polynomials[END_REF][START_REF] Kaliman | Two remarks on polynomials in two variables[END_REF][START_REF] Vistoli | The number of reducible hypersurfaces in a pencil[END_REF].

Structure of this paper. In Section 1 we recall some results about invariant algebraic curves and Jouanolou's theorem. In Section 2 we recall some classical results used in the proof of Theorem 2. In Section 3, we prove Theorem 2.

Notations. N = {0; 1; 2 . . .} is the set of integer numbers. |S| is the number of elements in the set S. In the following k will denote the degree of the derivation and d the degree of a rational first integral. In this paper, when we consider a rational function f /g, we always suppose that f and g are coprime. Furthermore, we recall that deg(f /g) means max{deg(f ), deg(g)}.

Invariant algebraic curves and Jouanolou's theorem

In 1878, G. Darboux [START_REF] Darboux | Mémoire sur les équations différentielles du premier ordre et du premier degré[END_REF] has given a strategy to find first integrals. One of the tools developed by G. Darboux is now called invariant algebraic curves and it will be the main ingredient in our proof. Definition 5. A polynomial f is said to be an invariant algebraic curve associated to D, if D(f ) = Λ.f , where Λ is a polynomial. The polynomial Λ is called the cofactor of f . There exist a lot of different names in the literature for invariant algebraic curves, for example we can find: Darboux polynomials, special integrals, eigenpolynomials, special polynomials, or second integrals. A lot of properties of a polynomial differential system are related to invariant algebraic curves of the corresponding derivation D, see e.g. [START_REF] Goriely | Integrability and nonintegrability of dynamical systems[END_REF][START_REF] Dumortier | Qualitative theory of planar differential systems[END_REF]. G. Darboux shows in [START_REF] Darboux | Mémoire sur les équations différentielles du premier ordre et du premier degré[END_REF] that if the derivation D has at least k(k + 1)/2 + 1 irreducible invariant algebraic curves then D has a first integral which can be expressed by means of these polynomials. More precisely the first integral has the following form: i f ci i where f i are invariant algebraic curves and c i are complex numbers. This kind of integral is called nowadays a Darboux first integral. If the all the c i belong to Z then we have a rational first integral. That is to say a first integral which belongs to C(X, Y ). The relation between rational first integral and invariant algebraic curves is given in the following proposition: Proposition 6. If f /g is a rational first integral then for all (λ : µ) ∈ P 1 (C), λfµg is an invariant algebraic curve and this curve corresponds to the level set f /g = µ/λ. Furtermore, the irreducible factors f (λ:µ),i of λfµg are also invariant algebraic curves.

Proof. The first part is a straightforward computation. The second part come from the fact that irreducible factors of an invariant algebraic curve are also invariant algebraic curves, see e.g. [START_REF] Dumortier | Qualitative theory of planar differential systems[END_REF]Proposition 8.4] or [START_REF] Goriely | Integrability and nonintegrability of dynamical systems[END_REF]Proposition 2.5].

As mentioned before we are going to use a theorem due to Jouanolou. This theorem is the following, see [START_REF] Jouanolou | Équations de Pfaff algébriques[END_REF]: Theorem 7. If a derivation has at least k(k+1)/2+2 irreducible invariant algebraic curves then there exist integers

n i ∈ Z such that i f ni i ∈ C(X, Y ) is a rational first integral.
Several authors have given simplified proof of Jouanolou's theorem. M. Singer proves this result in [START_REF] Michael | Liouvillian first integrals of differential equations[END_REF]. A direct proof of Jouanolou's result is also given in [START_REF] Cerveau | Formes intégrables holomorphes singulières[END_REF] and in [START_REF] Dumortier | Qualitative theory of planar differential systems[END_REF]. Darboux and Jouanolou theorems are improved in [START_REF] Llibre | Darboux theory of integrability for polynomial vector fields in R n taking into account the multiplicity at infinity[END_REF][START_REF] Llibre | Darboux theory of integrability in C n taking into account the multiplicity[END_REF]. The authors show that we get the same kind of result if we take into account the multiplicity of invariant algebraic curves. The multiplicity of an invariant algebraic curve is defined and studied in [START_REF] Christopher | Multiplicity of invariant algebraic curves in polynomial vector fields[END_REF]. This notion of multiplicity does not correspond to the multiplicity e (λ:µ),i of the irreducible factors of λfµg.

Darboux and Jouanolou's theorems have been also studied in the sparse case, i.e. bounds are given in terms of a Newton polytope in [START_REF] Chèze | Darboux theory of integrability in the sparse case[END_REF]. More precisely, the result is the following: Theorem 8. If a derivation has at least B+2 irreducible invariant algebraic curves, where B is defined as in Theorem 2, then there exist integers n i ∈ Z such that

i f ni i ∈ C(X, Y
) is a rational first integral. The key point in the proof of this theorem is the study of the Newton polygon of the cofactors for a given derivation, for other results about the structure of the cofactors see [START_REF] Ferragut | Seeking Darboux Polynomials[END_REF].

Indecomposable rational functions and remarkable values

Definition 9. We say that a rational function f /g ∈ C(X, Y ) is decomposable if we have f /g = u(h) where u(T ) ∈ C(T ) with deg(u) ≥ 2 and h(X, Y ) ∈ C(X, Y ), otherwise f /g is said to be indecomposable.

Sometimes indecomposable rational functions are called minimal rational functions.

Remark 10. When f /g = u(h) with h indecomposable, we can suppose that h = h 1 /h 2 with h 1 and h 2 irreducible in C[X, Y ] and deg(h 1 ) = deg(h 2 ). Indeed, let (λ 1 : µ 1 ), (λ 2 : µ 2 ) ∈ σ(h 1 , h 2 ), and set

H i = λ i h 1 -µ i h 2 . Then, we have deg(H 1 ) = deg(H 2 ) = deg(h), H 1 , H 2 are irreducible in C[X, Y ], and f = u v(H 1 /H 2 )
, where v ∈ C(T ) is an homography.

The following proposition explains why indecomposable functions are important in our situation. A consequence of this proposition is that any two indecomposable rational first integrals are equal, up to an homography.

In the situation of a polynomial first integral, we can find a proof of this proposition in [START_REF] Ferragut | On the remarkable values of the rational first integrals of polynomial vector fields[END_REF]Corollary 18]. However we have not found a reference for this result in the situation of a rational first integral. Thus we give a proof based on the following lemma, see [START_REF] Michael | Liouvillian first integrals of differential equations[END_REF]Lemma A.1]. This lemma means that at a nonsingular point of a vector field there is at most one algebraic solution. Another proof of Proposition 11 using Lüroth's theorem is posssible.

Lemma 12. Let D = A(X, Y )∂ X +B(X, Y )∂ Y be a derivation and let f 1 , f 2 be two invariant algebraic curves. Suppose that (x 0 , y 0 ) is a nonsingular point of D, i.e. A(x 0 , y 0 ) = 0 or B(x 0 , y

0 ) = 0. If f 1 (x 0 , y 0 ) = f 2 (x 0 , y 0 ) = 0 and f 1 is irreducible then f 1 divides f 2 .
Proof of Proposition 11. If ker D = C then there exists f /g ∈ C(X, Y )\C such that D(f /g) = 0. We can suppose that f /g is indecomposable. Indeed, if f /g = u(h), with deg(u) ≥ 2 and h indecomposable, then we have 0 = D(f /g) = u ′ (h)D(h). As deg(u) ≥ 2 we deduce u ′ (h) = 0 and D(h) = 0. Thus h is indecomposable and h ∈ ker D. Now, as f /g is indecomposable, we can suppose, thanks to Remark 10 and the previous computation, that f and g are irreducible in C[X, Y ] with deg(f ) = deg(g). Then, we can assume that f /g is a rational function with minimal degree in the set of all indecomposable rational functions F/G ∈ ker D such that F and G are

irreducible in C[X, Y ] with deg(F ) = deg(G). Now, let F 1 /F 2 ∈
i , β i ∈ C such that H i = α i f -β i g.
We prove the claim: As the number of singular points of D is finite we can consider (x i , y i ) ∈ C 2 , where i = 1, 2, are nonsingular points of D such that H i (x i , y i ) = 0.

we apply Lemma 12, to H i and g(x i , y i )f (X, Y )f (x i , y i )g(X, Y ) and we deduce that H i divides g(x i , y i )f (X, Y )f (x i , y i )g(X, Y ). As we have supposed the degree of f /g minimal, this proves the claim. Now, thanks to our claim, we get

H 1 /H 2 = V (f /g) where V ∈ C(T ) is an homography. Thus F 1 /F 2 = U V (f /g) . It follows ker D = C(f /g).
It follows that indecomposable first integrals correspond to first integrals with minimal degree.

Proof of Theorem 2

As mentioned before, the strategy will be the following: first construct a rational first integral with factors some of the f (λ,µ),i and then get a contradiction.

Proof of Theorem 2. Suppose that ρ(f, g) ≥ B + 2, we are going to prove that this situation is absurd. We denote by s the number of remarkable values in σ(f, g)\γ(f, g). Furthermore, we use the following convention: if Z divides f h (λi:µi) then we set f h (λi:µi),n(λi,µi) = Z. Then, f (λ1,µ1),1 , . . . , f (λ1,µ1),n(λ1:µ1)-1 , . . . , f (λs,µs),1 , . . . , f (λs,µs),n(λs:µs)-1 are non-trivial invariant algebraic curves by Proposition 6. As s i=1 n(λ i : µ i ) -1 = ρ(f, g) ≥ B + 2, Jouanolou's theorem in the sparse case, see Theorem 8, implies that these curves give a rational first integral. That is to say, there exist integers c (λ1:µ1),1 , . . . , c (λs:µs),n(λs:µs)-1 such that .g e .

As the factorization into irreducible factors in C[X, Y ] is unique, we deduce that (α l : β l ) and (γ m : δ m ) belong to the spectrum σ(f, g). Witout loss of generality we can suppose that (λ 1 , µ 1 ) = (α 1 , β 1 ). Then we deduce that equality ( * ) is impossible since the factor f (λ1:µ1),n(λ1:µ1) does not appear in the left hand side of the equality but must appear in the right hand side. This concludes the proof.

Figure 1 .

 1 Figure 1. Newton polygon N c 3,3 X 3 Y 3 + c 2,3 X 2 Y 3 + c 3,2 X 3 Y 2 + c 0,0

Figure 2 .

 2 Figure 2. The associated Newton polygon N D .

Proposition 11 .

 11 If we denote by ker D the set of rational first integrals of D we have ker D = C(f /g) where f /g ∈ C(X, Y ) is indecomposable, or ker D = C.

  :µj )-1 i=1 f c (λ j :µ j ),i (λj :µj ),i ∈ ker D.By Proposition 11, we gets j=1 n(λj :µj )-1 i=1 f c (λ j :µ j ),i (λj :µj ),i = u f g ,where u ∈ C(T ). Furthermore, u(f /g) is a rational function of this form:u f g = l (α l fβ l g) m (γ m fδ m g) .g e ,where α l , β l , γ m , δ m ∈ C and e ∈ Z. Thus we get the equality j :µ j ),i(λj :µj ),i = l (α l fβ l g) m (γ m fδ m g)

  a derivation, such that deg(A), deg(B) ≤ k, and let f /g ∈ C(X, Y ) be an indecomposable rational function which is a first integral of D. Consider a generic point (x, y) in C 2 and the Newton polygon N D associated to the

	Laurent polynomial x	A(X, Y ) X	+ y	B(X, Y )

Y , and denote by B the number of integer

  ker D, and set F 1 /F 2 = U (H 1 /H 2 ) where H 1 /H 2 is indecomposable. As before, we can show that D(H 1 /H 2 ) = 0 and suppose, thanks to Remark 10, that H 1 and H 2 are irreducible in C[X, Y ], and deg(H 1 ) = deg(H 2 ).

	Now we claim:
	(⋆) there exist α

Acknowledgment

The author thanks Thomas Cluzeau, Jacques-Arthur Weil and Antoni Ferragut for their precious comments during the preparation of this text. Furthermore, the author has greatly appreciate the comments of the anonymous referee. These comments have allow to improve the submitted paper.


Weizhi Lu

Weiyu Li

Kidiyo Kpalma

Joseph Ronsin

Compressed Sensing Performance of Random Bernoulli Matrices with High Compression Ratio

Keywords: compressed sensing, compression ratio, high dimension, random matrix, Bernoulli distribution

This paper studies the sensing performance of random Bernoulli matrices with column size n much larger than row size m. It is observed that as the compression ratio n/m increases, this kind of matrices tends to present a performance floor regarding the guaranteed signal sparsity. The performance floor is effectively estimated with the formula 1 2 ( πm/2 + 1). To the best of our knowledge, it is the first time in compressed sensing, a theoretical estimation is successfully proposed to reflect the real performance.

I. INTRODUCTION

Compressed sensing states that a sparse signal can be well acquired through an underdetermined system of linear equations [START_REF] Donoho | Compressed sensing[END_REF]. This technique successfully transfers the cost of acquisition to the later recovery process, which can be briefly described as follows. Let x ∈ R n be a k-sparse signal with at most k nonzero elements, and A ∈ R m×n be a sensing matrix with m < n. Then x can be recovered from a finite number of linear observations y = Ax. The recovery process can be formulated as an 0 -based minimization problem min ||x|| 0 subject to y = Ax [START_REF] Donoho | Compressed sensing[END_REF] which holds a unique solution, if x is sparse enough such that the sensing matrix A satisfies some well-known conditions, e.g. null space property (NSP) [START_REF] Cohen | Compressed sensing and best k-term approximation[END_REF] and restricted isometry property (RIP) [START_REF] Candes | Decoding by linear programming[END_REF]. Note that, for the convenience of analysis, the columns of A are typically normalized in the following study.

For compressed sensing, it is natural to seek a sensing matrix with high compression ratio n/m. Empirically, the sensing performance will inevitably degrade as the compression ratio increases. A question of practical interests then arise: how fast will the performance degrade with the increasing of compression ratio? This paper is developed to address this question for the random Bernoulli matrix, which is popular in compressed sensing and performs as well as Gaussian ones. Surprisingly, as will be shown in the final simulation, the random Bernoulli matrix approximately presents a 'performance floor' with the fast increasing of compression ratio. To be specific, the decreasing speed of guaranteed sparsity k is very slow, and can even be ignored in the setting where m is fixed while n tends to infinity. This property enables W. Lu, K. Kpalma and J. Ronsin are with INSA de Rennes, IETR, UMR CNRS 6164, F-35708, France (e-mail: weizhi.lu, kidiyo.kpalma and joseph.ronsin@insa-rennes.fr).

W. Li is with CREST-ENSAI, Bruz, F-35170, France (e-mail: weiyu.li@ensai.fr). the significant compression of high-dimensional signal with sparsity lower than the performance floor. Then it becomes interesting to theoretically evaluate the performance floor, namely the guaranteed sparsity k. Unfortunately, the exact performance estimation is still an open problem in compressed sensing. Currently, the guaranteed sparsity k is often simply estimated with

k < 1 2 (µ m (A) -1 + 1), (2)
where the parameter µ m (A) represents the maximum absolute correlation between distinct columns of A [START_REF] Donoho | Uncertainty principles and ideal atomic decomposition[END_REF]. However, as will be shown later, the formula ( 2) is only a sufficient condition for perfect recovery, such that the estimated k is usually much smaller than the real value we can achieve in practice. In this paper, we will prove that the formula (2) can be modified to be a sufficient and necessary condition in the limit, if the maximum correlation µ m (A) is replaced with the average absolute correlation between distinct columns of A, denoted by µ a (A). This improved estimation allows us to propose a simple formula 1 2 ( πm/2 + 1) to effectively approximate the performance floor of random Bernoulli matrices with fixed row size m. To the best of our knowledge, it is the first time that a theoretical estimation is reported being able to reflect the real sensing performance. Thus the contribution of this paper is of both practical and theoretical interests.

The rest of the paper is organized as follows. In the next section, by analyzing the proof process of formula (2), we demonstrate how the sufficient and necessary condition is approached by the average correlation. In section III, we first calculate the average correlation of random Bernoulli matrix, then estimate its performance floor. The numerical evidence is illustrated and discussed in section IV. Finally, this paper is concluded in section V.

II. ESTIMATION METHODS BASED ON AVERAGE CORRELATION VS. MAXIMUM CORRELATION

This section demonstrates how the formula (2) is modified to be a sufficient and necessary condition in the limit by analyzing its proof process shown in the following Lemma 2. To simplify the proof of Lemma 2, we first review an important conclusion in Lemma 1.

Lemma 1. [5]

A k-sparse signal x can be uniquely recovered from y = Ax with the 0 -minimization problem as shown in formula (1), if and only if the sparsity k < Spark(A)/2, where Spark(A) denotes the smallest number of linearly dependant columns of A.

Lemma 2. [6] If k < 1 2 (µ m (A) -1 + 1),
a k-sparse signal x is the unique recovery of y = Ax with the 0 -minimization problem as shown in formula (1).

Proof. We first need to prove a critical property, that is

Spark(A) ≥ 1 + µ m (A) -1 . Let subset ψ ⊆ {1, 2, ..., n}
with cardinality denoted as |ψ|, A ψ be a submatrix of A with columns indexed by ψ, and A T ψ be the transpose of A ψ . Consider the Gram matrix G = A T ψ A ψ ∈ R |ψ|×|ψ| . Clearly its diagonal elements g ii = 1, and off-diagonal elements

g ij ≤ µ m (A) ≤ 1, i = j.
Recall that the columns of A ψ are linearly independent if and only if Gram matrix G has positive determinant, equivalently each eigenvalue is positive. With Gershgorin circle theorem [START_REF] Horn | Matrix Analysis[END_REF], the i-th eigenvalue of G is bounded in the interval [g ii -r i , g ii + r i ], where

r i = |ψ| j=1;j =i |g ij |.
To render all eigenvalues positive, we only require 1 -(|ψ| -1)µ m (A) > 0. In other words, for any |ψ| < µ m (A) -1 + 1, the columns of A ψ are linear independent. From the definition of Spark, it can be easily deduced that Spark(A) ≥ µ m (A) -1 + 1. By merging the result with Lemma 1, we immediately derive the condition k < 1 2 (µ m (A) -1 + 1) for the unique solution based on 0minimization. Now we focus our attention on the analysis of the proof of Lemma 2. First, we need to show how the maximum correlation µ m (A) is involved. From the proof, it can be observed that the following inequality

1 - |ψ| j=1;j =i |g ij | > 0 (3)
must hold to ensure the Gram matrix A T ψ A ψ being positive definite, where g ij with i = j denotes the correlation between the i-th and j-th columns of submatrix A ψ . Considering

|g ij | ≤ µ m (A), the condition in formula (3) is then simply relaxed to 1 -(|ψ| -1)µ m (A) > 0. (4)
This relaxation process leads to the sufficient but unnecessary property of Lemma 2, which makes the estimation of Lemma 2 far away from the real performance.

To reflect the real performance, it is necessary to reduce the relaxation error between |ψ| j=1;j =i |g ij | and (|ψ| -1)µ m (A), such that the sufficient and necessary condition can be approached for Lemma 2. To this end, we propose to replace the maximum correlation µ m (A) with the average correlation µ a (A). In this case, as proved in Lemma 3, the relaxation error will be close to zero with high probability as the submatrix size |ψ| increases, if the average correlation of submatrix A ψ with high probability can be approximated by the average correlation of A. According to the law of large numbers, the condition above should be well satisfied by the random Bernoulli matrix, if the matrix column size is sufficiently large. Thus in the following part, we are motivated to verify the feasibility of the proposed estimation through evaluating the performance of random Bernoulli matrix with high compression ratio. Lemma 3. Suppose a i ∈ [0, t] is arbitrarily distributed with mean of e < t and variance of σ 2 . Then for k elements i.i.d drawn from the distribution, we have P r( k i=1 a i = ke) → 1 and P r(

k i=1 a i = kt) → 0, if k → ∞.
Proof. Suppose a binary distribution as below a i = 0 with probability 1 -p t with probability p where p = e/t. Then it follows that P r(a i = t) ≤ P r(a i = t), and P r(

k i=1 a i = kt) ≤ P r( k i=1 a i = kt) = p k = ( e t ) k
According to the law of large numbers, it is known that

lim k→∞ P r(| k i=1 a i -ke| < kε) ≥ 1 - σ 2 kε 2
where ε is an arbitrarily small positive constant. Then the conclusion of the lemma can be easily derived.

III. AVERAGE COLUMN CORRELATION OF RANDOM BERNOULLI MATRIX

To evaluate the performance floor, this section calculates the average column correlation of random Bernoulli matrices in Lemma 4. According to the law of large numbers, the average column correlation of random Bernoulli matrix with n m should be equivalent to the expected value of the absolute correlation between two arbitrary Bernoulli vectors. Therefore in Lemma 4 we only calculate the expected value mentioned above.

Lemma 4. Suppose v and w are two distinct normalized column vectors of random Bernoulli matrix A ∈ R m×n with i.i.d elements being ± 1 √ m equiprobably, and f (v, w) = |v T w| denotes the correlation between them, then the expected value of f is derived with the following two forms:

1)

E(f ) = 2 m 1 2 m m 2 m m 2 (5)
2)

lim m→∞ E( √ mf ) = 2 π ( 6
)
where * denotes the minimum integer not smaller than * .

Proof. First, following from

f = |v T w| = | m i=1 (v i w i )|,
f is equivalently written as

f = 1 m | m i=1 z i |,
where z i being ±1 equiprobably. Then with combination probability, the expected value of f can be formulated as

E(f ) = 1 m 1 2 m m i=0 ( m i |m -2i|)
where m i := m! (m-i)!i! . With m i |m -2i| =                  m m-1 0 if i = 0 m m-1 m-i-1 -m m-1 i-1 if 1 ≤ i ≤ m 2 m m-1 i-1 -m m-1 m-i-1 if m 2 < i < m m m-1 m-1 if i = m one can further derive that m i=0 ( m i |m -2i|) =    2m m-1 m 2 -1 if m is even 2m m-1 m-1 2 if m is odd Finally, with m-1 i-1 = i ) Then we have lim m→∞ E( √ mf ) = 2
π , whenever m is even or odd. The proof is completed.

Note that the formula (6) in fact converges very fast, and can be satisfied with a relatively small m (on the order of tens). This implies that the expected value E(f ) in formula [START_REF] Gribonval | Sparse representations in unions of bases[END_REF] allows to be approximately written as E(f ) = 2/(πm). For confirmation, the coincidence between E(f ) and 2/(πm) is illustrated in Figure 3. Then the average correlation µ a (A) of random Bernoulli matrices with n m can also be approximated as 2/(πm). In this case, the guaranteed sparsity on the performance floor of random Bernoulli matrices can be estimated with which is derived by replacing µ m (A) with µ a (A) = 2/(πm) in formula [START_REF] Cohen | Compressed sensing and best k-term approximation[END_REF]. As it is expected in the former section, k avr should be close to the real performance floor. Note that the maximum correlation of random Bernoulli matrices without same columns is equal to 1. For comparison, in the following simulation we also consider the performance floor estimated with µ m (A) = 1, which is simply derived as

k avr = 1 2 ( πm/2 + 1) (7)
k max = 1 2 (1 + 1) = 1. (8)
by incorporating µ m (A) = 1 into formula (2). As will be shown later, k avr performs much better than k max .

IV. NUMERICAL SIMULATIONS

In this section, we will show that the random Bernoulli matrix indeed approximately presents a performance floor, which can be effectively estimated with k avr . To illustrate the performance floor, the guaranteed sparsity k of random Bernoulli matrices with fixed m and increasing n is derived in Table I. Before analyzing the data, we first briefly introduce the simulation setting. The sparse signal x with sparsity k is randomly generated in each simulation, and recovered with subspace pursuit algorithm [START_REF] Dai | Subspace pursuit for compressive sensing signal reconstruction[END_REF]. The recovery rate is measured with 1-||x -x|| 2 /||x|| 2 . Note that here we only consider the largest k that can be recovered with rate larger than 0.99, because the perfect recovery is hard to be validated with simulation, unless we can enumerate all possible distributions of k nonzero elements. Each result in Table I is derived after 10000 simulation runs.

Note that in Table I the compression ratio n/m exponentially increases, while the decreasing speed of k is very slow and can even be ignored compared with the fast increasing n/m. This implies that there indeed approximately exists a performance floor regarding the fixed row size m, as graphically illustrated in Figure 2. Note that the gradual performance degradation is inevitable with the exponentially increasing of n. In practice, it seems hard to present a definite range for the performance floor. Here we empirically suppose the performance floor occurs, as the sparsity k decreases in a I, with m=10 2 , 12 2 and 14 2 .

step not greater than 1. The relevant results are underlined in bold in Table I. Note that, due to the limitation of computer memory, as shown in Table I, we cannot test enough samples n to describe the performance floor, especially as m increases. For instance, we obtain only five samples for the largest case of m = 16 2 . Thus to present a uniform metric, the performance floor of each m is quantified only with the mean of the first five k among the results underlined in Table I. The quantified result is denoted with p f . Note that the whole performance floor is hard to be accurately reflected with the metric p f , because it only considers the first five largest values on the performance floor. To be specific, as m increases, p f tends to behave better than the real performance floor, since with five samples it can only consider few relatively large k and ignores most other smaller k on the performance floor. In contrast, for relatively small m, p f probably performs worse than the real performance floor, as the five samples are likely to include some relatively small k. In this case, it is reasonable to infer that the distance between p f and the estimation k avr will increase with the increasing of m, if k avr is close to the real performance. This conjecture is validated by the following simulation results.

In Figure 3, we compare the two estimations k avr and k max against the performance floor measured with p f . As it is expected, the estimation k avr based on average correlation performs much better than the estimation k max based on maximum correlation. Precisely, the estimation k avr is very close to p f while the estimation k max is of no practical use. The infeasibility of k max results from the fact that it is derived by simply assuming all column correlations take the value of the maximum correlation. Obviously, this assumption is far from the truth and suffers from significant errors. Despite obtaining much better estimate than k max , as illustrated in Figure 3, k avr is still not perfect, which becomes smaller than p f as m increases. As discussed before, this is partially due to the inaccuracy of p f on measuring the real performance floor. Precisely, p f tends to behave worse than the real performance floor when m is relatively small, but perform better as m increases. From this fact, it can be derived that the error between k avr and the real performance floor should be smaller than the result shown in Figure 3.

V. CONCLUSION

This paper has shown that the random Bernoulli matrix approximately presents a performance floor with the increasing of compression ratio, which enables the significant compression of high-dimensional sparse signals. More importantly, we successfully estimated the performance floor by exploring the average correlation between distinct columns of random Bernoulli matrix, instead of the traditional maximum column correlation. By this means, the theoretical estimation is built on a sufficient and necessary condition derived in the limit, and its value is thus much closer to the real sensing performance compared with the traditional estimation only based on a sufficient condition. Empirically, the result of this paper is also valid for Gaussian random matrices.

Figure 1 :

 1 Figure 1: The values of E(f ) in formula (5) and 2/(πm) over varying m.

Figure 2 :

 2 Figure 2: The performance curves of three random Bernoulli matrices from TableI, with m=10 2 , 12 2 and 14 2 .

Figure 3 :

 3 Figure3: The performance floor measured with p f , and the theoretical estimations kavr in formula[START_REF] Horn | Matrix Analysis[END_REF] and kmax in formula[START_REF] Dai | Subspace pursuit for compressive sensing signal reconstruction[END_REF].

Table I :

 I The largest k guaranteed by (m, n)-sized random Bernoulli matrices under the condition of recovery rate larger than 0.99. For each m, all k with decreasing step smaller than 2 are underlined in bold.

		n/m	2 1	2 2	2 3	2 4	2 5	2 6	2 7	2 8	2 9	2 10	2 11	2 12	2 13
		6 2	5	3	3	2	2	2	1	1	1	1	1	1	1
		8 2	12	9	7	5	5	4	3	3	3	2	2	2	2
	m	10 2 12 2	22 34	16 24	12 19	10 16	8 13	7 11	6 10	5 9	5 8	4 7	4 6	4 6	3 5
		14 2	50	35	27	22	19	16	14	12	11	10	9	9	8
		16 2	68	50	38	30	26	22		17	15	14	13	12	11


Tiziano D'albis
email: tiziano.dalbis@univ-rennes1.fr

Claire Haegelen

Caroline Essert

Sara Fernández-Vidal

Florent Lalys

Pierre Jannin

Tiziano D ' Albis

PyDBS: An automated image-processing workflow for deep-brain stimulation surgery

Keywords: Surgical planning, Deep Brain Stimulation, 3D visualization

Purpose Deep-brain stimulation (DBS) is a surgical procedure for treating motor-related neurological disorders. DBS clinical efficacy hinges on precise surgical planning and accurate electrode placement, which in turn call upon several image-processing and visualization tasks, such as image registration, image segmentation, image fusion, and 3D visualization. These tasks are often performed by a heterogeneous set of software tools, which adopt differing formats and geometrical conventions, and require patient-specific parameterization or interactive tuning. To overcome these issues, we developed PyDBS, a fully-integrated and automated image-processing workflow for DBS surgery.

Methods PyDBS consists of three image-processing pipelines and three visualization modules assisting clinicians through the entire DBS surgical workflow, from the preoperative planning of electrode trajectories to the the postoperative assessment of electrode placement. The system's robustness, speed and accuracy were assessed by means of a retrospective validation, based on 92 clinical cases.

Results

The complete PyDBS workflow achieved satisfactory results in 92% of tested cases, with a median processing time of 28 minutes per patient.

Conclusion

The results obtained are compatible with the adoption of PyDBS in clinical practice.

Introduction

High-frequency deep brain stimulation (DBS) is a minimally-invasive surgical therapy for treating motor-related disorders (e.g., Parkinson's disease, essential tremor, and dystonia), which has recently been extended to the treatment of severe psychiatric disorders (e.g., obsessive-compulsive disorder and severe depression). DBS consists in implanting electrodes in the basal ganglia, located deep in the patient's brain, and delivering highfrequency electric pulses to the surrounding neural tissues. The reversibility of this therapy, along with its spectacular outcomes, have largely contributed to its success. It is important to note that both the quality of clinical improvement and the occurrence of adverse effects have been found to correlate with electrode placement accuracy [START_REF] Tisch | Effect of electrode contact location on clinical efficacy of pallidal deep brain stimulation in primary generalised dystonia[END_REF][START_REF] York | Relationship between neuropsychological outcome and dbs surgical trajectory and electrode location[END_REF][START_REF] Lalys | Anatomo-clinical atlases correlate clinical data and electrode contact coordinates: Application to subthalamic deep brain stimulation[END_REF]. Nevertheless, accurate positioning of DBS electrodes is still a challenging and laborious procedure, requiring extensive surgical, anatomical, and functional knowledge of the brain, and typically relying on several manual adjustments and tuning [START_REF] Dormont | Neuroimaging and deep brain stimulation[END_REF].

Electrode trajectories are planned prior to surgery based on preoperative magnetic resonance (MR) imaging. This is a complex task, given that the brain structures targeted for stimulation are often not directly visible on conventional MR sequences due to limitations in image contrast and spatial resolution. Surgeons therefore often rely on multiple co-registered MR sequences or digital anatomical atlases for targeting. Once a trajectory has been planned, a frame-based stereotactic procedure is conducted. To express the planned trajectory in stereotactic coordinates, the frame is detected on a MR or computer tomography (CT) scan of the patient's skull (after frame fixation), and this scan is then registered to the MR image used for targeting. During surgery, the planned trajectory is refined by means of neurophysiological microelectrode recordings (MERs) and direct observation of the stimulation's clinical outcomes. Following surgery, the actual position of the electrode contacts is assessed by segmenting the corresponding artifacts on a postoperative CT scan. It is essential to know the exact position of the electrode contacts after surgery for the purposes of both lead programming and clinical research studies. The DBS surgical workflow described above requires several image-processing and visualization tasks, such as intra-modal and mono-modal image registration, intensity-based and atlas-based image segmentation, and image fusion and 3D visualization. The accuracy, robustness, and usability of the software tools employed for these tasks are all critical elements for the success of the entire surgical workflow. At present, most clinical centers rely on commercial software systems for DBS surgery. These systems, however, often restrict their functionalities to the planning phase of the workflow, demonstrate a lack of state-of-the-art registration and segmentation algorithms, and typically require several manual adjustments before the final result can be obtained. In medical-imaging research, great effort is put into the design, implementation, and validation of novel image-processing algorithms, which address very specific computational problems and are recognized for their ever-increasing performance levels. Nevertheless, only in a few instances have research studies taken a broad perspective at the surgical-workflow level. This kind of wide view is necessary to bridge the gap between research and application, with the aim of developing a software system that is effective in clinical practice. In this respect, effort should be spent in the integration and automation of the different algorithms and software tools available, in the interest of developing a coherent and usable software environment. This coherent software environment would then also be beneficial for the research itself by providing a shared platform where different algorithms can be easily compared, or even combined for solving problems at higher levels of abstract complexity. The study presented in this paper contributes to this goal.

Related work

Several solutions have been proposed to assist physicians through the different stages of DBS surgical workflow.

The study by Guo et al. [START_REF] Guo | Visualization and navigation system development and application for stereotactic deep-brain neurosurgeries[END_REF] has presented a visualization and navigation system for stereotactic DBS procedures.

This system aids physicians in the planning phase by registering and fusing anatomical and functional atlases to preoperative patient images and automatically detecting the stereotactic frame on MR or CT scans. Intraoperatively, the system simulates the microelectrode descent along the stereotactic trajectory and collects MER data for the construction of functional atlases. Similar functionalities have also been implemented into the Cicerone system [START_REF] Miocinovic | Cicerone: stereotactic neurophysiological recording and deep brain stimulation electrode placement software system[END_REF]. Despite this system not providing an automatic frame-detection algorithm, it does offer the additional functionality of visualizing a predicted volume of tissue activated by the stimulation. Neither of these systems, however, addresses the postoperative phase of the workflow, where the implanted electrodes and stimulation contacts are to be detected.

Another software environment developed in order to provide assistance for DBS surgery is CranialVault [START_REF] Dhaese | Cranialvault and its crave tools: A clinical computer assistance system for deep brain stimulation (DBS) therapy[END_REF].

This software combines a centralized data repository with a series of software tools for automatic data processing and visualization (CRAVE tools). The planning module includes functionalities for rigid registration of patient preoperative MR and CT scans and for non-linear deformation of multiple anatomical atlases in patient space. These atlases are used for both atlas-based segmentation of brain structures and localization of two anatomical landmarks used in stereotactic neurosurgery: the anterior commissure (AC) and the posterior commissure (PC). Furthermore, a proposed stimulation target is estimated by projecting target centroids from previously implanted patients onto patient space, with electrophysiological and somatotopic maps also provided.

By means of a standard-angle approach, the system provides pre-computed implantation plans that can be modified or validated by the surgeon. CranialVault and the CRAVE tools also include an intraoperative module with similar features to those implemented in the first two studies mentioned [START_REF] Guo | Visualization and navigation system development and application for stereotactic deep-brain neurosurgeries[END_REF][START_REF] Miocinovic | Cicerone: stereotactic neurophysiological recording and deep brain stimulation electrode placement software system[END_REF], along with a postoperative module for electrode contact detection.

Even though it includes a large set of functionalities, CranialVault does not provide segmentations of the brain vessels or cortical sulci. Given that brain vessel damage, particularly that which affects the large arteries lying within the cortical sulci, is a major risk factor in DBS surgery, it is highly desirable to segment these structures and visualize them during preoperative surgical planning. Furthermore, as CranialVault relies on a patient-specific stereotactic frame, it does not provide functionalities for the segmentation and detection of standard stereotactic frames, such as the Leksell G-frame (Elekta AB, Stockolm, Sweden) or CRW frame (Radionics, Burlington, MA, USA). It should also be noted that, to date, only the planning module has been preliminary validated [START_REF] Dhaese | Cranialvault and its crave tools: A clinical computer assistance system for deep brain stimulation (DBS) therapy[END_REF]. This consisted in measuring the accuracy of a) co-registration of patient's preoperative images, b) detection of points AC and PC, and c) prediction of the stimulation targets, by ensuring planning consistency across surgeons [START_REF] Pallavaram | Neurologist consistency in interpreting information provided by an interactive visualization software for deep brain stimulation postoperative programming assistance[END_REF]. Nevertheless, a validation of the segmented brain structures has not been reported.

Contribution

To the best of our knowledge, a complete, fully-integrated and automatic image-processing workflow for planning and postoperative assessment of DBS interventions has not yet been reported, despite all the previous related works. Furthermore, software's robustness, speed and accuracy have never been assessed in a systematic validation study. In this paper we present PyDBS, an automated image-processing workflow aiding surgeons and clinicians in the planning and assessment of DBS procedures. Leveraging on state-of-the-art algorithms and software tools, PyDBS provides a fully-integrated and automatic software environment for DBS surgery. The system's robustness, speed and accuracy were evaluated through a retrospective study involving 92 clinical cases. The results obtained are compatible with the adoption of PyDBS in clinical practice.

Materials and methods

PyDBS comprises three image-processing pipelines and three visualization/interaction modules. Each pipeline is executed at a different phase of the DBS surgical workflow, with an associated 3D scene collating all results generated at completion. The 3D scenes, specified in the Medical Reality Markup Language (MRML), can be readily visualized by means of a 3D Slicer [START_REF] Fedorov | 3D Slicer as an image computing platform for the quantitative imaging network[END_REF]10]. Using the MRML format, the full set of parameters required for 3D visualization, such as screen layout, colors, transparencies, camera viewpoints, and geometrical transformations, can, in fact, be pre-computed and stored within the scene file itself. Furthermore, the user can easily interact with a generated 3D scene by means of the three PyDBS visualization/interaction modules implemented as 3D Slicer plugins. serves as definition of the space of possible entry points for the implantation, with the cortical sulci and brain ventricles representing structures to be avoided by the electrode trajectory [START_REF] Zrinzo | Avoiding the ventricle: a simple step to improve accuracy of anatomical targeting during deep brain stimulation: Clinical article[END_REF], and the basal ganglia the stimulation targets. The T2-weighted image is then also registered onto the T1-weighted image for an improved visualization of the deep basal ganglia nuclei (e.g., the red nucleus). The pipeline output consists of an MRML scene containing the registered MR images, in addition to the masks and meshes of the segmented anatomical structures. Based on this scene, and by using the dedicated PyDBS plugin, the surgeon validates the results and interactively defines the electrode trajectories for the implantation (Fig. 2).

PyDBS in the surgical workflow

On the day of surgery, a surgical stereotactic frame is fixed on the patient's skull and head CT scan is performed. This scan, along with a geometrical 3D model of the surgical frame, provides the input to the PyDBS preoperative pipeline. At this point, the head CT scan is registered onto the T1-weighted MR volume, the skull segmented, and the surgical frame detected. With the benefit of frame detection, the planned trajectories originally defined in the patient's AC-PC space can be expressed in the coordinates of the stereotactic frame. The pipeline output consists of an MRML scene containing the results of both the inclusion and preoperative pipelines. Based on this scene, the surgeon validates the MR-CT registration and frame detection, then retrieves the stereotactic coordinates of the planned trajectories. The surgeon is also able to visualize the estimated positions of the five micro-electrodes used for intraoperative trajectory refinement. Furthermore, if a trajectory angle has been modified during surgery, the new trajectory can be readily updated and visualized in the space of the preoperative MR image.

Following surgery, a postoperative head CT scan is performed. This image, along with a geometrical 3D model of the implanted electrodes, represents the input to the PyDBS postoperative pipeline. At this point, the postoperative CT is registered onto the preoperative CT and projected onto the patient's AC-PC space.

Furthermore, the implanted electrodes are segmented based on the postoperative CT image, and the stimulation contacts are localized. Lastly, an MRML scene containing the results of the three pipelines is generated. Based on this scene, the surgeon assesses the electrode contact positions and compares the real trajectories with the planned ones.

System perspective

From a computing system perspective, PyDBS operates from input data consisting of patient-specific information, such as patient images and clinical data, and generic models, such as an anatomical atlas, stereotactic frame model, and electrode model, then outputting a patient-specific model for the planning and postoperative assessment of DBS surgery (Fig. 3). The generated patient-specific model consists of images, segmented anatomical structures (i.e., binary masks and meshes), and geometrical transformations (i.e., linear registration matrices and non-linear deformation fields). All images, masks, and meshes are mapped to a shared reference space (the patient's AC-PC space).

Data import

All patient-specific data is managed by an internal relational database. An automatic import procedure populates the database with relevant image metadata extracted from DICOM media (e.g., modality, volume size, and sequence type), converting the DICOM series into NIfTI format [12]. NIfTI images are stored in predefined locations on the file-system and then referenced by the internal database. Clinical data, such as the brain structure targeted for stimulation, the lateralization, and electrode type, can be either automatically imported from a clinical database or interactively provided at the beginning of a pipeline execution. Additional patient clinical data (e.g., clinical scores) can also be imported into the database. All data in the database is easily accessible from a web-based application or a graphical user interface, both of which are provided by PyDBS.

Pipeline execution

PyDBS pipelines are executed from the command-line or by means of a simple graphical user interface. At the beginning of each pipeline, the user is asked to select the patient image(s) to be processed from a list of matching entries in the internal database. Following this initial selection, a pipeline is automatically executed with no need for further user interaction (with the exception of manually localizing AC and PC points in the inclusion pipeline). On completion, the results, namely the images, masks, meshes, and geometrical transformations, are stored in a predefined folder on the local file-system, with an associated MRML scene description file generated for each pipeline.

Inclusion pipeline

The inclusion pipeline is composed of six processing steps (Fig. 4, left). Firstly, the two MR images of the patient brain (T1-weighted and T2-weighted) are imported from the internal database. Secondly, the user localizes the anatomical landmarks AC and PC on the T1-weighted image. This is the only step out of all three pipeline procedures that requires user interaction. We opted for the manual detection of AC and PC due to the high inter-surgeon variability that characterizes this localization [START_REF] Pallavaram | Intersurgeon variability in the selection of anterior and posterior commissures and its potential effects on target localization[END_REF]. Furthermore, there are several localization strategies currently in use in clinical practice, and we decided not to constrain the users in their choice. Once the AC-PC reference frame has been defined, a batch execution of the BrainVISA Morphologist pipeline [START_REF] Rivière | A freely available anatomist/brainvisa package for structural morphometry of the cortical sulci[END_REF][START_REF]BrainVISA Morphologist[END_REF] is triggered. This pipeline performs an intensity-based segmentation of the scalp, brain (gray and white matter), and cortical sulci. We selected the BrainVISA Morphologist pipeline due to its proven robustness for segmenting the cortical sulci [START_REF] Rivière | A freely available anatomist/brainvisa package for structural morphometry of the cortical sulci[END_REF], a particularly important segmentation for the purposes of our work. The T2-weighted and T1weighted MR images are then subject to rigid co-registration. All linear registrations are computed by means of the FMRIB's Linear Image Registration Tool (FLIRT), which has been proven to outperform other commonly-used linear-registration algorithms in terms of accuracy and robustness [START_REF] Jenkinson | A global optimisation method for robust affine registration of brain images[END_REF][START_REF] Jenkinson | Improved optimization for the robust and accurate linear registration and motion correction of brain images[END_REF]. For both mono-and multimodal linear-registrations, we have adopted a cost function based on mutual information.

Atlas-based segmentation

The T1-weighted image of the patient's brain is registered onto an anatomical atlas for performing an atlas-based segmentation of the brain ventricles and basal ganglia. The atlas we have adopted was constructed by averaging 57 T1-weighted MR images of Parkinson's disease patients, defining 24 anatomical structures, all manually segmented by an anatomist [START_REF] Haegelen | Automated segmentation of basal ganglia and deep brain structures in MRI of Parkinson's disease[END_REF]. Atlas registration is performed over five processing steps (Fig. 5), including both linear and non-linear registrations and taking a similar global-to-local approach to the one previously described by Lalys et al. [START_REF] Lalys | Post-operative assessment in deep brain stimulation based on multimodal images: registration workflow and validation[END_REF]. Firstly, all non-brain voxels are masked out of both images, followed by rigid coregistration of both (six degrees of freedom [DOF]). Secondly, the two images are affinely co-registered (12 DOF). Thirdly, a box volume corresponding to a region-of-interest (ROI) around the basal ganglia and defined in the anatomical atlas space is cropped from both images. Fourthly, the two cropped ROIs are linearly coregistered (12 DOF). Finally, the ROIs are non-linearly co-registered using the symmetric diffeomorphic registration algorithm (SyN) with cross-correlation as cost function [START_REF] Avants | Symmetric diffeomorphic image registration with cross-correlation: Evaluating automated labeling of elderly and neurodegenerative brain[END_REF][START_REF]Advanced Normalization ToolS[END_REF]. We chose this deformation algorithm based on a recent validation study that reported SyN to be one of the best performing algorithms for iconic non-linear brain image registration [START_REF] Klein | Evaluation of nonlinear deformation algorithms applied to human brain MRI registration[END_REF]. At the end of the registration process, brain-ventricle and basalganglia segmentations are projected from the atlas space onto the patient space (Fig. 6). In the final step, an MRML scene is generated.

Preoperative pipeline

The preoperative pipeline is composed of six processing steps (Fig. 4, middle). Firstly, the preoperative CT scan of the patient's head (fixed to the stereotactic frame) is imported from the internal database. Secondly, this CT image is rigidly registered onto the preoperative T1-weighted MR image of the patient's brain. Thirdly, an intensity-based segmentation of the patient's skull and the four frame screws is performed by means of a chain of thresholding and connected-component-labeling operations (Fig. 7a,top). In the fourth step, the frame artifact is segmented using an analogous strategy (Fig. 7a,bottom). In the fifth step, the frame model is rigidly registered onto the patient's CT scan. This is done by means of a two-step procedure, the first step consisting of ten landmark points automatically detected on the frame artifact, namely the two most superior points for each frame side and the three points at half of the frame height (Fig. 7b). These landmarks are used to estimate an initial rigid registration matrix between frame model and the patient's CT scan (point-based optimization algorithm [START_REF] Arun | Least-squares fitting of two 3-d point sets[END_REF]). Then, the two images are rigidly-co registered using the estimated matrix as initialization. At the final step, a second MRML scene is generated with these results. unsatisfactory, with respect to their potential use for preoperative planning of electrode trajectories and postoperative assessment of electrode placement. If any element within the 3D scene was considered misleading for the purposes of these two tasks, the corresponding validation step was marked as unsatisfactory. Pipeline results were considered satisfactory only if all their validation steps were satisfactory.

The inclusion pipeline was validated in three steps: 1) validation of the intensity-based brain segmentation, 2) validation of the T2-weighted image registration, and 3) validation of the atlas-based segmentation of the brain ventricles and basal ganglia (Fig. 9 a-c). The preoperative pipeline was validated in two steps: 1) validation of preoperative CT registration, and 2) validation of frame detection (Fig. 9 d-e). Lastly, the postoperative pipeline was validated in two steps: 1) validation of postoperative CT registration, and 2) validation of electrode segmentation and electrode contact detection (Fig. 9 f-g).

Results

Validation results are reported in Fig. 10a. The inclusion pipeline achieved satisfactory results for all patients except two (#75 and #127), corresponding to a 98% success rate. For patient #75, the atlas-based segmentation produced unsatisfactory results. This was due to the patient exhibiting very large ventricles compared to the anatomical atlas, resulting in a mismatch between the region-of-interest cropped from the patient image and the corresponding region on the atlas. For patient #127, the intensity-based brain segmentation terminated with errors. This can probably be attributed to the presence of movement artifacts in the original T1-weighted MR volume. Given this pipeline's unsuccessful termination, the preoperative and postoperative pipelines could not be executed for this patient.

The preoperative pipeline produced satisfactory results for 86 patients out of the 91 tested, corresponding to a 94% success rate. The pipeline failed on one occasion due to an unsatisfactory registration of the preoperative CT with the preoperative MR (patient #79), and on four occasions due to an unsatisfactory frame detection (patients #39, #54, #89, and #116). Where frame detection failed, the problem was owing to an inaccurate segmentation of the frame artifact caused by non-frame voxels that were not masked out from the image. In all cases, the preoperative pipeline terminated with no reported errors. Satisfactory results were obtained with the postoperative pipeline for 90 patients out of the 91 tested, corresponding to a 99% success rate. For one patient (#89), the pipeline failed due to an unsatisfactory registration of the postoperative CT with the preoperative CT. For this patient, despite the electrode contacts being correctly detected, their coordinates could not be correctly expressed in the AC-PC reference frame. In all cases, the postoperative pipeline terminated with no reported errors.

Combining the results of the three pipelines, satisfactory results were obtained for 85 patients out of the 92 tested, corresponding to a 92% success rate. Pipeline processing times are reported in Fig. 10b. Median processing times were approximately 15 minutes for the inclusion pipeline, 7 minutes for the preoperative to serve as a comparison. This was the case for our study, in which we aimed to validate a complete system at the surgical-workflow level. The algorithms and software tools that are used by the three PyDBS pipelines had, in fact, been individually validated for the most part in separate studies, though their level of integration and automation, as well as their effectiveness in the specific context of DBS surgery, had yet to be assessed. This could have been conducted prospectively by evaluating the potential effects of adopting the proposed system on the DBS surgical workflow, in terms of its total duration or the associated clinical outcomes, for example. A prospective validation, however, also implies several clinical risks, and it would be preferably postponed to the late stages of a system evaluation when its performances have been already proved otherwise. We therefore opted for a retrospective validation, and measured our system's performances in terms of robustness, speed, and accuracy. In terms of accuracy, however, only a qualitative evaluation was possible owing to the lack of groundtruth for comparison. An alternative could have been to quantitatively assess system accuracy in comparison with a gold-standard. This would have been very difficult in our study, however, due to the heterogeneity of surgical procedures adopted at different clinical centers, and the substantial engineering effort that would have been required to reproduce and adapt alternative solutions to our data.

For qualitative validation, one would ideally involve a pool of experts and evaluate the results at the root of the intra-observer and inter-observer variability in the assessment. This was also difficult in our study, however, owing to the cost of involving several experts in such a demanding validation study, and the fact that only one expert was available at our site. We thus believe that a multi-observer, and possibly even multi-site validation study, is a necessary milestone we shall aim for in the future, even though the validation proposed in this paper proved very informative for future developments of the project.

Overall, the results we obtained proved PyDBS to be a valuable tool for assisting clinicians through the DBS surgical workflow. Furthermore, we plan to employ the PyDBS inclusion pipeline output as input data for an automatic-trajectory-planning software, such as the one proposed by Essert et al. [START_REF] Essert | Automatic computation of electrode trajectories for deep brain stimulation: a hybrid symbolic and numerical approach[END_REF]. With the advantage of its modularity and flexibility, PyDBS is also well-suited for supporting clinical research studies. For example, the PyDBS postoperative pipeline, supplemented by additional functionalities not described here, has already been employed in a clinical research study on electrode deformations [25], as well as in a study focusing on the correlation between stimulation targets and clinical scores [START_REF] Lalys | Anatomo-clinical atlases correlate clinical data and electrode contact coordinates: Application to subthalamic deep brain stimulation[END_REF].

Future projects also include the integration of anatomo-clinical atlases [START_REF] Lalys | Anatomo-clinical atlases correlate clinical data and electrode contact coordinates: Application to subthalamic deep brain stimulation[END_REF] into the PyDBS inclusion pipeline.

Once registration with preoperative patient's images has been conducted, anatomo-clinical atlases could, in fact, be used in order to identify the best stimulation targets according to the clinical outcomes of previously implanted patients. In this setting, PyDBS could be used both to visualize the atlases during planning and to build and update them whenever a new patient is treated. Finally, it is our objective to validate PyDBS prospectively by using it in parallel with the commercial DBS-planning system currently in use at our institute.

Fig. 1

 1 Fig. 1 displays typical usage of PyDBS within the DBS surgical workflow. Prior to surgery, T1-weighted and T2weighted MR images of the patient's brain are acquired. The two images, along with a digital anatomical atlas, are the input to the PyDBS inclusion pipeline. The inclusion pipeline, following interactive localization of the anatomical landmarks AC and PC, then performs an intensity-based segmentation of the scalp, brain, and cortical sulci, as well as an atlas-based segmentation of the brain ventricles and basal ganglia. The scalp surface

  25.Lalys, F., Haegelen, C., Dalbis, T., Jannin, P.: Analysis of electrode deformations in deep brain stimulation surgery. Int J Comput Assist Radiol Surg 9(1),107-117 (2014) .

Fig. 1

 1 Fig. 1 Integration of PyDBS within the DBS surgical workflow.

Fig. 2

 2 Fig. 2 MRML scene generated by the PyDBS inclusion pipeline. The three slice views are presented in "probeeye" mode, that is, projected perpendicularly to the electrode trajectory (red cylinder). The left panel shows the PyDBS plugin for interactive trajectory planning. The user defines a trajectory by selecting an entry and target point on the MR image (red spheres, E: entry point, T: target point).A trajectory can be refined by dragging the two points on the 3D view, such as for maximizing the distance from the cortical sulci (gray meshes). Multiple trajectory plans can be defined, and plans can be saved to an XML file for later inspection.

Fig. 3

 3 Fig. 3 PyDBS from a system perspective. PyDBS inputs are patient-specific data and generic models. The PyDBS output is a patient-specific model for DBS surgery.

Fig. 4

 4 Fig. 4 Detailed diagram of the PyDBS pipelines: processing modules and dependencies. Dashed arrows depict inter-pipeline dependencies and indicate the specific result on which the target processing module depends. Within-pipeline dependencies are not shown.

Fig. 5

 5 Fig. 5 Registration between the T1-weighted patient image and anatomical atlas. The contour of the atlas brain (red line) is superimposed onto an axial slice of the patient image. The darker square in d-f indicates the region of interest (ROI). Initially, the two volumes are not registered (a). The registration process is composed of five phases (b-f): rigid registration (b), affine registration (c), crop of the ROI (d), affine registration of the two ROIs (e), and non-linear registration of the two ROIs (f).

Fig. 6

 6 Fig. 6 Atlas-based segmentation of brain ventricles (blue) and basal ganglia: caudate (cyan), putamen (magenta), globus pallidus externum (green), globus pallidus internum (dark orange), thalamus (yellow), sub-thalamic nucleus (light orange), red nucleus (red), and substantia nigra (gray). Two common targets for DBS are the globus pallidus internum (arrows in a-c) and sub-thalamic nucleus (arrows in d-f ). The first three columns show axial, sagittal, and coronal slices, respectively. In a-c only the T1-weighted MR volume is shown. In d-f the T2weighted volume is overlaid onto the T1-weighted volume for improved visualization of deep structures. g) 3D view of brain ventricles and basal ganglia on the T1-weighted volume.

Fig. 8

 8 Fig. 8 Electrode segmentation and contact detection. Panels a-c depict three slices of one patient's postoperative CT image magnified onto the electrode artifact and overlaid with the contours of detected objects, namely the electrode artifact (blue), electrode contacts (red), and three surrounding anatomical structures: putamen (magenta), globus pallidus externum (green), and globus pallidus internum (orange). The slice in a) is perpendicular to the electrode axis and centered at the lowest contact. The slices in b) and c) are aligned to the electrode axis and perpendicular to each other. In d) a 3D view is displayed.

Fig. 9

 9 Fig. 9 Validation steps for the three pipelines. Each step is illustrated with a screen-shot of the 3D scene visualized by the observer. a) Skin (magenta), gray matter (white), white matter (red), and cortical sulci (yellow) are superimposed onto the T1-weighted volume. b) The T2-weighted volume is superimposed onto the T1weighted volume. c) Brain ventricles (blue) and basal ganglia are superimposed onto the T1-weighted volume (basal-ganglia colors have been presented in Fig. 5). d) The skull segmented on the preoperative CT (yellow) is superimposed onto the T1-weighted volume. e) The frame artifact (yellow) and the registered frame model (red) are superimposed onto the preoperative CT. f) The skull segmented on the preoperative CT (red) along with the gray matter segmented on the T1-weighted volume (yellow) are superimposed onto the postoperative CT. g) The electrode artifact (blue) and electrode contacts (red) are superimposed onto the postoperative CT. Registration and segmentation are abbreviated as "reg." and "seg.", respectively.

Fig. 10

 10 Fig. 10 Validation results. a) For each validation step and for each patient an empty (filled) circle indicates a satisfactory (unsatisfactory) result. For the three pipelines we reported the fraction and percentage of patients with entirely satisfactory validation steps. When one of the pipelines terminated with a runtime error, the corresponding patient code was marked in red. b) Pipeline processing times. Colored bars represent median values; error bars span from the lower to the upper quartile of the distribution. Patient #127, for whom runtime errors occurred, was excluded from time performance estimation.

Acknowledgements

The authors thank the French National Research Agency (ANR) who founded this work through the ACouStiC project grant (ANR 2010 BLAN 0209 01).

Postoperative pipeline

The postoperative pipeline is composed of six processing steps (Fig. 4, right). Firstly, the postoperative CT is imported in from the internal database. Secondly, this image is thresholded in order to obtain a coarse skull segmentation. In the third module, the postoperative CT is rigidly registered onto the preoperative CT, a process which takes two steps, namely a global-search registration between the two segmented skulls, followed by a local-search registration between the two original volumes. In the fourth step, the electrode artifacts are segmented on the postoperative CT by applying the brain masks computed on the T1-weighted MR image (one per hemisphere) and by thresholding the resulting volume. In the fifth step, the electrode contacts are detected by registering the geometrical model of the implanted electrodes to the segmented artifacts. The electrode tip position is calculated using the coordinates of the most inferior voxel in the electrode artifact, and the direction of the electrode axis is estimated by applying principal-component analysis to the artifact voxels within a 15mm radius from the electrode tip (Fig. 8).

Validation

The three pipelines were retrospectively validated against a dataset comprising 92 patients available at our site.

We included all patients having been implanted from September 2006 to September 2012 at the Neurosurgical Department of the University Hospital of Rennes (France). The inclusion criteria consisted of availability of: 1) 3T T1-weighted MR image with gadolinium injection (1mm x 1mm x 1mm, Philips Medical system), 2) 3T T2weighted MR image (1mm x 1mm x 1.1mm, Philips Medical system), 3) preoperative CT scan with stereotactic frame (0.55mm x 0.55mm x 0.6mm, GE Healthcare VCT 6), and 4) postoperative CT scan with implanted electrodes (0.44mm x 0.44mm x 0.6mm, GE Healthcare VCT 64). Informed consent was obtained from all patients. Despite there being only T1-weighted MR images acquired after frame fixation available for a subset of patients (n=37), in all cases, the stereotactic frame was detected from a preoperative CT scan. The dataset comprised patients having undergone unilateral (n=14) or bilateral (n=78) DBS surgery. The targeted anatomical structures were the sub-thalamic nucleus (n=37), globus pallidus internum (n=32), or the caudal part of the ventro-lateral thalamic nucleus (n=23). Patient pathologies included Parkinson's disease (n=63), essential tremor (n=12), dystonia (n=9), Tourette syndrome (n=4), and obsessive-compulsive disorder (n=4). For all patients, the Leksell G stereotactic frame (Elekta Instruments AB, Stockholm, Sweden) has been used, and the implanted electrodes were Medtronic 3389 or 3387 electrodes (Medtronic Sofamor Danek).

The pipelines were executed on a 64-bit laptop with an Intel® i7-2860QM CPU (8 cores, 2.50GHz clock frequency) equipped with 16GB of RAM. The results were qualitatively validated by an expert, with the processing times measured for each pipeline. This qualitative validation consisted of seven steps (three for the inclusion pipeline, two for the preoperative pipeline, and two for postoperative pipeline), each supplemented by a 3D scene visualized by the expert (Fig. 9). The corresponding results were evaluated as either satisfactory or pipeline, and 4 minutes for the postoperative pipeline, with a total processing time of 28 minutes per patient. The processing times for the most time-consuming modules of each pipeline have also been reported. For the inclusion pipeline, the most time-consuming module was the intensity-based brain segmentation (Morphologist pipeline), which took approximately 7 minutes to complete, followed by the atlas registration (approximately 6 minutes), and the interactive detection of the AC and PC points (approximately 2 minutes). For both the preoperative and postoperative pipelines, the most time-consuming module was the CT volume registration, which took a median time of 3.5 minutes to complete.

Discussion and future work

In this article, we have presented PyDBS, an integrated image-processing workflow for DBS surgery. PyDBS includes three image-processing pipelines and three visualization/interaction modules. The three pipelines are fully automatic (exception made for the manual localization of the points AC and PC) and completely general, meaning that no patient-specific parameter tuning is needed. Each pipeline addresses a specific phase of the surgical workflow, generating as output images, meshes, and geometrical transformations, all collected in a 3D scene readily available for visualization. Furthermore, for each 3D scene, a PyDBS visualization/interaction module provides an easy means to access and manipulate the displayed data.

We retrospectively tested all three pipelines on a dataset of 92 patients available at our site, with the results subjected to qualitative validation by an expert. Satisfactory results were obtained across the entire workflow for 85 patients out of the 92 tested, corresponding to a 92% success rate. The cases where pipelines failed were analyzed in order to further improve our software. Particular mention should be made that the frame-detection module was found to be the least robust, as it resulted in the highest error rate (4%). In this respect, we believe that a new semi-automatic procedure for frame detection should be implemented for the cases where the automatic algorithm fails. In more general terms, we plan to supplement the three PyDBS visualization/interaction modules with additional error-correction functionalities, ensuring that the user is able to interactively adjust the automatically-generated results, whenever this may be required.

PyDBS is implemented in the Python programming language, hence its name, language which was chosen based on the clarity of its syntax and its rapidity of development, making it well-suited for research prototyping.

The processing speed was not, however, compromised, as all low-level image-processing tasks, such as registration and segmentation, were delegated to pre-compiled routines. We were, in fact, able to obtain a median processing time of about 28 minutes per patient, and the preoperative pipeline, the most critical in terms of time performance, completed in about 7 minutes. This is a time performance that is compatible with clinical adoption. Next we will discuss some of the limitations facing our work. Firstly, the system performances were evaluated by means of a qualitative validation study. Qualitative validation is commonly used in medical-imaging research when quantitative performance measures are difficult to define, or when a ground truth is not directly available

Conflict of interest

The authors declare that they have no conflict of interest. All images are binary, except for the input CT volume. Thresholding and connected-component-labeling operations are abbreviated as "th" and "cc", respectively. For segmenting skull and screws we retained only the largest connected component, while for the frame artifact we retained the six largest ones. This is on account of the frame artifact being composed of three separate parts (i.e., one anterior part and two lateral parts), each of which may be split into another two parts if the frame is not entirely contained in the CT volume. When segmenting the frame artifact, the dilated brain mask is used to clean up small skull components. b) Segmentation results for one patient. The patient's skull (gray), frame screws (blue), and frame artifact (yellow) are all displayed. The anterior part of the frame artifact is not shown, since it was not included in the frame model. The red spheres (five per side) indicate the automatically-detected landmark points used for the pointbased registration between the frame model and the CT scan.