Khalid Koufany

Genkai Zhang
email: genkai@math.chalmers.se

HUA OPERATORS, POISSON TRANSFORM AND RELATIVE DISCRETE SERIES ON LINE BUNDLE OVER BOUNDED SYMMETRIC DOMAINS

Keywords: Bounded symmetric domains, Shilov boundary, invariant differential operators, eigenfunctions, Poisson transform, Hua systems

Let Ω = G/K be a bounded symmetric domain and S = K/L its Shilov boundary. We consider the action of G on sections of a homogeneous line bundle over Ω and the corresponding eigenspaces of G-invariant differential operators. The Poisson transform maps hyperfunctions on S to the eigenspaces. We characterize the image in terms of twisted Hua operators. For some special parameters the Poisson transform is of Szegö type whose image is in a relative discrete series; we compute the corresponding elements in the discrete series.

April 2011

Introduction

Let X = G/K be a Riemannian symmetric space of non compact type. It is known that the Poisson transform maps certain parabolically induced representation spaces into null spaces of some systems of differential equations. For a minimal parabolic subgroup P min ⊂ G, Kashiwara et al. proved [START_REF] Kashiwara | Eigenfunctions of invariant differential operators on a symmetric space[END_REF], that the Poisson transform gives a Gisomorphism from the set of hyperfunctions on the maximal boundary G/P min onto the joint eigenspace of invariant differential operators on X, thus proving the Helgason conjecture [START_REF] Helgason | A duality for symmetric spaces with applications to group representations[END_REF]. We shall be interested in the case of a Hermitian symmetric space G/K and the Poisson transform corresponding a maximal (instead of minimal) parabolic subgroup P max ⊂ G with G/P being the Shilov boundary of G/K. For a certain special parameter of the induced representation the image of the transform is a subspace of harmonic functions on symmetric space G/K. The precise description of the image for tube domains is given in [START_REF] Johnson | Generalized Hua-operators and parabolic subgroups. The cases of SL(n,C) and SL(n,R)[END_REF][START_REF] Johnson | The Hua operators on bounded symmetric domains of tube type[END_REF][START_REF] Lassalle | Les équations de Hua d'un domaine borné symétrique du type tube[END_REF] in terms of Hua-harmonic functions introduced earlier by Hua [START_REF] Hua | Harmonic analysis of functions of several complex variables in the classical domains[END_REF]. Its generalization to non-tube cases is done by Berline and Vergne [START_REF] Berline | Équations de Hua et noyau de Poisson. Noncommutative harmonic analysis and Lie groups[END_REF]. For tube domains with general parameters Shimeno [START_REF] Shimeno | Boundary value problems for the Shilov boundary of a bounded symmetric domain of tube type[END_REF] proved an analogue of Kashiwara et al. theorem for P max ⊂ G. More precisely, he proved that the Poisson transform is a G-isomorphism from the space of hyperfunctions on the Shilov boundary onto the space of eigenfunctions of the Hua operator of the second order. The generalization to the non-tube bounded symmetric domains has been given in our earlier paper [START_REF] Koufany | Hua operators on bounded symmetric domains]Hua operators and Poisson transform for bounded symmetric domains[END_REF].

A more interesting problem is to consider homogeneous line bundles over Ω with the corresponding weighted action of G. In this setting Shimeno [START_REF] Shimeno | Eigenspaces of invariant differential operators on a homogeneous line bundle on a Riemannian symmetric space[END_REF] generalized the Kashiwara et al. theorem to homogeneous line bundles over Hermitian symmetric spaces of tube type G/K for minimal parabolic P min ⊂ G. That is for a given a line bundle E ν (see below) and for a generic parameter depending on ν of the induced representation from a minimal parabolic P min ⊂ G. subgroup, the Poisson transform maps as isomorphisms from hyperfunction-valued sections of a line bundle over G/P min onto a space of eigenfunctions. We shall find characterizations of the Poisson integrals of hyperfunctionvalued sections of a line bundle over the Shilov boundary G/P min of a bounded symmetric domain for generic parameters. The Poisson transform becomes more interesting for larger parameters of ν as there appear relative discrete series, in particular, the weighted Bergman spaces, in the Plancherel formula [START_REF] Shimeno | The Plancherel formula for the spherical functions with onedimensional K-type on a simply connected simple Lie group of hermitian type[END_REF], the Poisson transform on Shilov boundary for the corresponding parameter is obviously not injective. We shall compute explicitly the image for some of the relative discrete series. We proceed with some more precise description of our result.

Let Ω = G/K be a bounded symmetric domain of tube type of rank r and genus p. For ν ∈ pZ we consider the (unique) character τ ν of K and the corresponding homogeneous line bundle E ν over Ω. We identify

C ∞ -sections of E ν with the space C ∞ (G/K, τ ν ) of C ∞ -functions on G such that f (gk) = τ ν (k) -1 f (g).
We consider the generalized Poisson transform (P s,ν f )(z) = S P s,ν (z, u)f (u)du where P s,ν is the generalized Shilov kernel and S = G/P 1 the Shilov boundary of Ω. The subgroup P 1 is a maximal parabolic subgroup of G. For s ∈ C, let B(S, s, ν) be the space of hyperfunction-valued sections on S = G/P 1 associated with the character of P 1 given by man → e (sρ 0 -ρ 1 )(log a) τ ν (m). Then for λ s = ρ + 2n(s -1)ξ * e -νrξ * e , the space B(S, s, ν) can be considered as a subspace of B(G/P ; L λs,ν ) of hyperfunction-valued sections of a line bundle over G/P . We construct certain Hua operators on G/K and we prove (Theorem 5.2) that for generic values of s the Poisson transform P s,ν is a G-isomorphism between B(S, s, ν) and the space of eigenfunctions of the Hua operator of the second order.

If Ω is a non-tube type domain it is known that Hua-type operators of second order will not be sufficient to characterize the image of Poisson transform on the Shilov boundary. However for type one non-tube domains of r × (r + b)-matrices the Lie algebra k C is a direct sum of gl r and sl r+b and one can construct [START_REF] Berline | Équations de Hua et noyau de Poisson. Noncommutative harmonic analysis and Lie groups[END_REF] a second order Hua operator by taking certain projection on the summand gl r . We prove a corresponding result for line bundles in this case; see 8.1.

For singular value of s we prove (Theorem 9.3) that the Poisson transform is a Szegö type map of principal series representation onto the relative discrete series representation. We compute explicitly the Poisson transform on certain spherical polynomials on the Shilov boundary.

The paper is organized as follows. In §2 we recall very briefly the Jordan algebraic characterization of bounded symmetric domains. In §3 we introduce the line bundle over the bounded symmetric domain Ω. The generalized Poisson transform of hyperfunction-valued sections on the maximal and the Shilov boundaries are studied in §4. The characterization of Poisson integrals of hyperfunction-valued sections on the Shilov boundary is given in §5. In this section we also recall our geometric construction of the Hua operator. The necessary condition is proved in §6. In §7 we compute the radial part of the Hua system and prove sufficiency condition. Finally in §8 we show a relationship between the Poisson transform, Hua operator and the relative discrete series representation.

After a preliminary version of this paper was finished we were informed by Professor T. Oshima that he and N. Shimeno have obtained in [START_REF] Oshima | Boundary value problems on Riemannian Symmetric Spaces of the noncompact Type[END_REF] some similar results about Poisson transforms and Hua operators. Professor A. Koranyi communicated also his recent preprint [START_REF] Korányi | Poisson transform for line bundles from the Shilov boundary to bounded symmetric domains[END_REF] to us where he proved the necessity of Theorem 5.2 using different methods. In particular some of the questions posed in that paper are answered here.

Bounded symmetric domains and Jordan triples

We begin with a brief review of necessary facts on bounded symmetric domains and Jordan triple systems. Let V , dim V = n be a complex vector space, Ω ⊂ V a irreducible bounded symmetric domain. Let Aut(Ω) be the group of all biholomorphic automorphisms of Ω, let G be the connected component of the identity of Aut(Ω), and let K be the isotropy subgroup of G at the point 0 ∈ Ω. As a symmetric space, Ω = G/K. The group K acts as linear transformations on V and we can thus identify K also as a subgroup of GL(V ).

Let g be the Lie algebra of G and g C its complexification. The algebra g has the Cartan decomposition g = k ⊕ p and g C = p + + k C + p -the Harish-Chandra decomposition. Denote Z 0 the element in the center of k which defines the complex structure on p + , i.e., ad(Z 0 )v = iv for v ∈ p + . We can thus identify p + with V , p + = V .

There exists a quadratic map Q : V → End( V , V ) (where V is the complex conjugate of V ), such that p = {ξ v ; v ∈ V } as holomorphic vector fields, where

ξ v (z) := v -Q(z)v. Define D(z, v)w := (Q z+w - Q z -Q w )v. It satisfies D(z, v)w = D(w, v)z , [D(u, v), D(z, w)] = D({u v z}, w)-D(z, {w ū v}) so V is a Jordan triple system. Furthermore we have [X, ξ z ] = ξ Xz for X ∈ k, z ∈ V and [ξ z , ξ v ] = D(z, v) -D(v, z) for all z, v ∈ V .
In this realization elements in p -are of the form {-Q(z)v} which we write as v. Thus [START_REF] Berline | Équations de Hua et noyau de Poisson. Noncommutative harmonic analysis and Lie groups[END_REF] [v, w] = D(z, w).

We define

(2) z, w = 1 p tr D(z, w),
where tr is the trace functional on End(V ) and p is the genus defined below. It is a K-invariant Hermitian product on V . The group K acts on V by unitary transformations. The domain Ω is realized as the open unit ball of V with respect to the spectral norm,

(3) Ω = {z ∈ V : D(z, z) 2 < 2},
where D(z, z) is the operator norm of D(z, z) on the Hilbert space (V, •, • ).

An element e ∈ V is a tripotent if {e ē e} = e. The subspaces V λ (e) = ker(D(e, ē) -λid) are called Pierce λ-spaces. Then we have

V = V 0 (e) ⊕ V 1 (e) ⊕ V 2 (e).
Two tripotents e and c are orthogonal if D(e, c) = 0. A tripotent e is minimal if it cannot be written as the sum of two non-zero orthogonal tripotents. With the above normalization of inner product we have e, e = 1 for minimal tripotents e. The tripotent e is maximal if V 0 (e) = 0.

A frame is a maximal family of pairwise orthogonal, minimal tripotents. It is known that the group K acts transitively on frames. In particular, the cardinality of all frames is the same, and it is equal to the rank r of Ω.

Let us choose and fix a frame {e j } r j=1 in V . Then, by transitivity of K on the frames, each element z ∈ V admits a polar decomposition z = k r j=1 s j e j , where k ∈ K and s j = s j (z) are the singular numbers of z. Denote e the maximal tripotent e = e 1 + . . . + e r . The Shilov boundary of Ω is S = K/K 1 where K 1 = {k ∈ K : k e = e}. It is know that S coincides with the set of maximal tripotents of V .

The joint Peirce spaces are

(4) V ij = {z ∈ Z : D(e k , ēk )z = (δ ik + δ jk )z, ∀k}, 0 ≤ i ≤ j ≤ r then V = 0≤i≤j≤r V ij , V 00 = 0, and V ii = Ce i (i > 0).
The triple of integers (r, a, b) with ( 5)

a := dim V jk (1 ≤ j < k ≤ r); b := dim V 0j (1 ≤ j ≤ r)
is independent of the choice of the frame and uniquely determines the Jordan triple. Notice that b = 0 exactly if V is a Jordan algebra which is equivalent to say that Ω is of tube type.

The Peirce decomposition associated with e is then

V = V 2 ⊕ V 1 where (6) V 2 = 1≤j≤k≤r V jk V 1 = r j=1 V 0j .
Let

n 1 = dim V 1 and n 2 = dim V 2 , then n 1 = rn, n 2 = r + r(r -1) 2 a and n = n 1 + n 2 .
The genus of Ω is

p = p(Ω) = 1 r tr D(e, ē) = (r -1)a + b + 2.
Let a = Rξ e 1 + . . . + Rξ er . Then a is a maximal abelian subspace of p with basis vectors {ξ e 1 , . . . , ξ er }. Its dual basis in a * will be denoted by {

β j 2 } r j=1 ⊂ a * , i.e., (7)
β j (ξ e j ) = 2δ jk , 1 ≤ j, k ≤ r.
We define an ordering on a * via (8)

β r > β r-1 > • • • > β 1 > 0.
It is known that the restricted roots system Σ(g, a) of g relative to a is of type C r or BC r , it consists of the roots ±β j (1 ≤ j ≤ r) with multiplicity 1, the roots

± 1 2 β j ± 1 2 β k (1 ≤ j = k ≤ r)
with multiplicity a, and possibly the roots

± 1 2 β j (1 ≤ j ≤ k) with multiplicity 2b. The set of positive roots Σ + (g, a) consists of 1 2 (β j ± β k ) (1 ≤ j < k ≤ r), β j and 1 2 β j (1 ≤ j ≤ r).
The half sum of positive roots is given by

(9) ρ = r j=1 ρ j β j = r j=1 b + 1 + a(j -1) 2 β j .
Let A be the analytic subgroup of G corresponding to a. Let n + = α∈Σ + g α and n -= θ(n + ). Let N + and N -be the corresponding analytic subgroups of G. Let M be the centralizer of a in K. The subgroup P = MAN + is a minimal parabolic subgroup of G.

The set Λ = {α 1 , . . . , α r-1 , α r } of simple roots in Σ + is such that

α j = 1 2 (β r-j+1 -β r-j ), 1 ≤ j ≤ r -1 and α r = β 1
for the tube case 1 2 β 1 for the non-tube case. Let Λ 1 = {α 1 , . . . , α r-1 } and P 1 the corresponding standard parabolic subgroup of G with the Langlands decomposition

P 1 = M 1 A 1 N + 1 such that A 1 ⊂ A. Then the Lie algebra of A 1 is a 1 = Rξ e where ξ e = ξ e 1 + . . . + ξ er .
The spaces G/P = K/M is the Furstenberg or maximal boundary, and

G/P 1 = K/K 1 , with K 1 = M 1 ∩ K, is the Shilov boundary of G/K.

Line bundle over Ω

Denote Z = p n Z 0 where Z 0 is the center element defined in §2. The group K is factorized as K = exp(RZ)K s , where K s is the analytic subgroup of K with Lie algebra

k s = [k, k].
For a fixed ν ∈ pZ consider the character τ ν of K defined by

τ ν (k) = e itν if k = exp(tZ) ∈ exp(RZ) 1 if k ∈ K s .
In particular we have J k (z)

ν p = e itν for k = exp(tZ), z ∈ Ω. Thus J k (z) ν p = τ ν (k), k ∈ K.
Here we denote J g , g ∈ G, the Jacobian of the holomorphic mappings g on Ω. See [START_REF] Schlichtkrull | One-dimensional K-types in finite dimensional representations of semisimple Lie groups : A generalization of Helgason's theorem[END_REF] and [START_REF] Dooley | Relative discrete series of line bundles over bounded symmetric domains[END_REF].

Let E ν be the homogeneous line bundle G × K C over G/K = Ω, where K acts on C via the one dimensional representation τ ν . The space C ∞ (Ω, E ν ) of smooth sections of E ν is by definition the space C ∞ (Ω, E ν ) = C ∞ (G/K; τ ν ) of C ∞ -functions F on G such that F (gk) = τ ν (k) -1 F (g).
We will trivialize the bundle via the map [(g, c)] ∈ E ν → (g•0, J g (0) ν p c) ∈ Ω×C and thus identifies C ∞ (G/K; τ ν ) also as the space of C ∞ -functions on Ω with G acting as [START_REF] Johnson | The Hua operators on bounded symmetric domains of tube type[END_REF] g

∈ G : f (z) → J g -1 (z) ν p f (g -1 z).
Let D ν (G/K) denote the space of G-invariant differential operators on G/K acting on C ∞ (G/K; τ ν ). We have the Harish-Chandra isomorphism [START_REF] Shimeno | Eigenspaces of invariant differential operators on a homogeneous line bundle on a Riemannian symmetric space[END_REF] γ

ν : D ν (G/K) ≃ U(a) W
where U(a) W denote the set of W -invariant elements in the enveloping algebra U(a).

The characters of D ν (G/K) are given by

χ λ,ν (D) = γ ν (D)(λ), D ∈ D ν (G/K), λ ∈ a * C . For λ ∈ a * C we define A(G/K, M λ,ν
) to be the space of functions ϕ ∈ C ∞ (G/K; τ ν ) satisfying the system of differential equations [START_REF] Kashiwara | Eigenfunctions of invariant differential operators on a symmetric space[END_REF] M λ,ν :

Dϕ = χ λ,ν (D)ϕ D ∈ D ν (G/K).
) (P λ,ν f )(g) = K f (gk)τ ν (k)dk g ∈ G for any f ∈ B(G/P, L λ,ν ).
As elements f ∈ B(G/P, L λ,ν ) are uniquely determined by its restriction to K it is natural to express the integral above as on K. Indeed, for g ∈ G denote κ(g) ∈ K and H(g) ∈ a to be elements uniquely determined by

g ∈ κ(g) exp(H(g))N + ⊂ KAN + = G. Then (13) (P λ,ν f )(g) = K f (k)τ ν (κ(g -1 k))e -(λ+ρ)(H(g -1 k)) dk and maps B(G/P, L λ,ν ) into A(G/K, M λ,ν ), see [23, Theorem 5.2]
The Harish-Chandra c function

c ν (λ) = N - e (λ+ρ)(H(n)) τ ν (κ(n))dn
can be written as

c ν (λ) = f ν (λ) e ν (λ) ,
where the denominator e ν (λ) is given, in our normalization, by

e ν (λ) = j>k Γ( 1 2 a+λ j ±λ k ) j Γ( 1 2 (b+1+2λ j +ν))Γ( 1 2 (b+1+2λ j -ν));
see [START_REF] Schlichtkrull | One-dimensional K-types in finite dimensional representations of semisimple Lie groups : A generalization of Helgason's theorem[END_REF] and [23, page 227].

The following theorem characterizes the range of the Poisson transform. (1 -t 2 j ).

Let h(z, w) be its polarization, i.e. h(z, w) is holomorphic on z and antiholomorphic in w such that h(z, z) = h(z). For any complex number s and for any ν we define the generalized Poisson kernel ( 16)

P s,ν (z, u) = h(z, z) |h(z, u)| 2 s n r h(z, u) -ν , z ∈ Ω, u ∈ S.
We also define the generalized Poisson transform [START_REF] Oshima | A definition of boundary values of solutions of partial differential equations with regular singularities[END_REF] (P s,ν f

)(z) = S P s,ν (z, u)f (u)du, for f ∈ B(S) Let a ⊤ 1 = r j=1 R(ξ j -ξ j+1 )
be the orthogonal complement of a 1 = Rξ e in a with respect to the Killing form. We denote ξ * e the dual vector such that ξ * e (ξ e ) = 1 and we extend ξ * e to a by the orthogonal projection defined above. Then ρ 1 = ρ |a 1 the restriction of ρ to a 1 is given by ρ 1 = nξ * e . Define ρ 0 to be the linear form on a 1 such that ρ 0 = rξ * e .

Consider the following representation of P 1 = M 1 A 1 N 1 given by σ s,ν = τ ν | M 1 ⊗ e sρ 0 -ρ ⊗ 1 and let B(G/P 1 , s, ν) be the space of hyperfunction valued sections of the line bundle on S = G/P 1 corresponding to σ s,ν , i.e., the space of hyperfunctions f on G satisfying

f (gman) = e (sρ 0 -ρ 1 )(log a) τ ν (m) -1 f (g) for all g ∈ G, m ∈ M 1 , a ∈ A 1 , n ∈ N + 1 .
We may fix e ∈ S as a base point and identify also B(G/P 1 , s, ν) with B(S).

For s ∈ C we define λ s ∈ a * by ( 18)

λ s = ρ + 2n(s -1)ξ * e -νrξ * e
Under the identification of C ∞ (G/K; τ ν ) (and thus its subspace A(G/K, M λ,ν )) as smooth functions on Ω we have P λs,ν coincides with P s,ν . We omit the routine computations. Let ν ∈ pZ, then we have

(19) B(G/P 1 ; s, ν) ⊂ B(G/P ; L λs,ν ) therefore ( 20
) P s,ν (B(G/P 1 ; s, ν)) ⊂ A(G/K; M λs,ν ).

Hua operators

The Bergman reproducing kernel of Ω is h(z, z) -p up to a constant. It is also [START_REF] Satake | Algebraic structures of symmetric domains[END_REF] h(z, z)

-p = det B(z, z) -1 where B(z, w) = I -D(z, w) + Q(z)Q( w)
is the Bergman operator. Thus Ω = G/K is a Kähler manifold with the (normalized) Bergman metric

(22) (u, v) z = 1 p ∂ u ∂v log h(z, z) -p = B(z, z) -1 u, v .
Let τ be a finite-dimensional holomorphic representation of K C . Let E be the Hermitian vector bundle over Ω associated with τ . Then there exists a unique connection operator ∇ :

C ∞ (Ω, E) → C ∞ (Ω, E ⊗ T ′
) compatible with the Hermitian structure and the anti-holomorphic differentiation, where T ′ = T ′ Ω is the cotangent bundle over Ω. That is, under the splitting in holomorphic and antiholomorphic parts, T ′ = T ′(1,0) ⊕ T ′(0,1) , we have ∇ = D + ∂ where D is the differentiation operator on E,

(23) D : C ∞ (Ω, E) → C ∞ (Ω, E ⊗ T ′(1,0) ).
Let D be the invariant Cauchy-Riemann operator on E defined in [START_REF] Engliš | Covariant Cauchy-Riemann operators and higher Laplacians on Kähler manifolds[END_REF][START_REF] Zhang | Shimura invariant differential operators and their eigenvalues[END_REF], by [START_REF] Shimeno | The Plancherel formula for the spherical functions with onedimensional K-type on a simply connected simple Lie group of hermitian type[END_REF] Df = B(z, z) ∂f.

Then (25) D : C ∞ (Ω, E) → C ∞ (Ω, E ⊗ T (1,0) ).
We will use the following identifications, T

′(1,0) z = p -= V ′ = V and T (1,0) z = T ′(0,1) = p + = V .
We now specialize the above to the line bundle E ν associated with the one-dimensional representation τ ν . The Hua operator H ν is then defined as the resulting operator of the following diagram

C ∞ (Ω, E ν ⊗ p + ) D / / O O D C ∞ (Ω; E ν ⊗ p + ⊗ p -) Ad p + ⊗p - C ∞ (Ω, E ν ) Hν / / C ∞ (Ω, E ν ⊗ k C ) .
We may call H = H ν the twisted Hua operator to differ it from the trivial case ν = 0. Remarks 5.1.

(1) We can change the order of D and D and define another Hua operator,

H ′ f = Ad p + ⊗p -DDf,
and one can prove, by direct computations, the following relation [START_REF] Shimeno | Boundary value problems for various boundaries of Hermitian symmetric spaces[END_REF] (H -H ′ )F = -2n r νF Id.

(2) On other hand, following Johnson-Korànyi [START_REF] Johnson | The Hua operators on bounded symmetric domains of tube type[END_REF], the Hua operator can also be described using the enveloping algebra : Let {E α } be a basis of p + and {E * α } be the dual basis of p + with respect to the Killing form. Then the Hua operator is defined as element of U(g

) C ⊗ k C by (27) H = α,β E α E * β ⊗ [E β , E * α ],
as operator acting from

C ∞ (E ν ) to C ∞ (E ν ⊗ k C ).
However the previous definition using D and D allows direct computation using coordinates on Ω and hence has some computational (and also conceptual) advantage.

Note that H is by definition G-invariant with respect to the actions of G on two holomorphic bundles, i.e., [START_REF] Upmeier | Toeplitz operators on bounded symmetric domains[END_REF] H(J g (z)

ν p f (gz)) = J g (z) ν p Ad(dg(z) -1 )(Hf )(gz),
where Ad(dg(z) -1 ) stands for the adjoint action on k C of dg(z) -1 : p + = T

(1,0) gz → p + = T

(1,0) z on p + , via the defining action of k C on p + . (Indeed the element dg(z) -1 is in the group K C ⊂ GL(V ) = GL(p + ) via the defining realization K ⊂ GL(V ), thus Ad(dg(z) -1 ) acts on k C ; see [21, Chapt II, Lemma 5.3].)

We state now our main result for case of tube domains. The proof is given in the next two sections. Theorem 5.2. Let Ω be a bounded symmetric domain of tube type. Suppose s ∈ C satisfies the following condition

(29) 4n(1 -s) r ∈ Λ 1 ∪ Λ 2
where

Λ 1 = Z + -2ν + 2, Λ 2 = 2Z ≥ -4ν + 4.
Then the Poisson transform P λs,ν is a G-isomorphism of B(S; s, ν) onto the space of analytic functions F on Ω such that

(30) HF = 2 n r s( n r (s -1) + ν)F Id.
Here Id stands for the element -iZ 0 , which acts on p + as identity.

The necessary condition of the Hua equations

The necessity in the above theorem is a consequence of the following Theorem 6.1. For any u ∈ S the function z → P s,ν (z, ū) satisfies the following system of equations

(31) HP s,ν = 2 n r s[ n r (s -1) + ν]P s,ν Id
Proof. We prove first the claim for z = 0, in which case we have

Hf (0) = α,β ∂ α ∂β f (0)[e α , ēβ ].
Recall the following formulas in [14, Lemma 5.2]: for any fixed w ∈ V and any complex number s, [START_REF] Zhang | Nearly holomorphic functions and relative discrete series of weighted L 2 -spaces on bounded symmetric domains[END_REF] ∂h(w, z) s = -sh(w, z) s w z, ∂h(z, w) s = -sh(z, w) s wz where w z = B(w, z) -1 (w -Q(w)z), wz = w z are called quasi-inverses of w with respect to z and w with respect to z respectively, viewed as (1, 0)-form and (0, 1)-form in terms of the Hermitian inner product (2). Then we have

∂P s,ν (z, ū) = ∂ h(z, u) -ν h(z, z) |h(z, u)| 2 n r s = -( n r s)P s,ν (z, ū)[z z -u z],
and [START_REF] Schlichtkrull | One-dimensional K-types in finite dimensional representations of semisimple Lie groups : A generalization of Helgason's theorem[END_REF] the last derivative can be written as

-( n r s)∂ (P s,ν (z, ū)[z z -u z]) = -( n r s)(ν + n r s)P s,ν (z, ū) [z z -u z] ⊗ ūz +( n r s) 2 P s,ν (z, ū) [z z -u z] ⊗ zz -( n r s)P s,ν (z, ū) ∂[z z -u z ]. Using
∂[z z -u z] = ∂(z z ) = ∂ ∂ log h(z, z) -1 = B(z, z) -1 Id
where Id is the identity form in p + ⊗ p -= p + ⊗ (p + ) ′ . Evaluating at z = 0 we get then, by the commutation relation ( 1) Thus the claim is true for general z by the invariant property (28) of H.

HP s,ν (0, ū) = ( n r s)(ν + n r s)D(u, ū) +0 -( n r s) n α=1 D(e α , ēα

The sufficiency condition of the Hua equations

The aim of this section is to prove that, using the radial part of the Hua operator, each solution F of the system (30) satisfies the system of equations [START_REF] Kashiwara | Eigenfunctions of invariant differential operators on a symmetric space[END_REF]. Then under the condition [START_REF] Yan | A class of generalized hypergeometric functions in several variables[END_REF] it can be proved that the boundary value of F is contained in B(S; s, ν).

We start to show that eigenfunctions of the Hua operator ( 30) are all τ ν -spherical functions. We need only to prove the claim for K-invariant functions F on Ω, i.e., F (kz) = F (z), k ∈ K.

For that purpose we compute the radial part of the Hua operator. The functions F (z) will be identified as permutation invariant and even function F (t) on the diagonal z = r j=1 t j e j . The operator H has the form

HF (z) = r j=1 H j F (z)D(e j , e j ),
for some operators where the scalar-valued operators H j are given by

H j in t = (t 1 , • • • , t r ).
H j =(1 -t 2 j ) 2 ∂ 2 ∂t 2 j + 1 t j (1 + 2ν) ∂ ∂t j + + a 2 k =j (1 -t 2 j )(1 -t 2 k ) 1 t j -t k ∂ ∂t j - ∂ ∂t k + 1 t j + t k ∂ ∂t j + ∂ ∂t k + +(-2ν)(1 -t 2 j ) 1 t j ∂ ∂t j .
Proof. Recall [START_REF] Zhang | Shimura invariant differential operators and their eigenvalues[END_REF] that the operator D acting on T (1,0) -valued function takes the form

DF (z) = h -ν (z, z)B(z, z) k ∂ k (h ν (z, z)B(z, z) -1 F (z)) ⊗ vk which is a p + ⊗ p --valued function. Thus HF (z) = Ad p + ⊗p - k,l h -ν (z, z)B(z, z) k ∂ k (h ν (z, z)( ∂l F (z))v l ) ⊗ vk ,
where we have performed the cancellation of the B -1 -term in D with B in D. Performing the Leibniz rule for ∂ k we get the above as sum of two terms, say I and II where

II = Ad p + ⊗p - k,l B(z, z) k ∂ k (( ∂l F (z))v l ) ⊗ vk and I = Ad p⊗p k,l h -ν (z, z)B(z, z) k (∂ k h ν (z, z))(( ∂l F (z))v l ) ⊗ vk
The second order term II is computed as in [START_REF] Koufany | Hua operators on bounded symmetric domains]Hua operators and Poisson transform for bounded symmetric domains[END_REF] or [START_REF] Faraut | Analysis on Symmetric Cones[END_REF] and we need only to treat I. Recall that

∂ v h ν (z, z) = ∂ v e ν log h(z,z) = -νe ν log h(z,z) ∂ v (-log h(z, z)) = -νh ν (z, z)(v, z z )
where z z = ∂(log h(z, z)) is quasi-inverse of z with respect to z; see [START_REF] Loos | Bounded symmetric domains and Jordan pairs[END_REF], [START_REF] Zhang | Nearly holomorphic functions and relative discrete series of weighted L 2 -spaces on bounded symmetric domains[END_REF]. We have then

I = -ν k,l (v k , z z )( ∂l F (z))[B(z, z)v l , vk ].
For z = r j=1 t j e j we have z z = r j=1 t j 1-t 2 j e j , and

B(z, z) = r j=1 (1 -t 2 j ) 2 D(e j , ēj )
which is diagonalized under the Peirce decomposition of V . In particular the above sum reduces to

I = -ν r j=1 t j 1 -t 2 j ( ∂e j F (z))(1 -t 2 j ) 2 D(e j , ēj ).
Being a K-invariant function F (z) is rotation invariant on the plan Ce j and we have then ∂e j F (z) = 1 2 ∂ t j F (t). Finally we have

I = - ν 2 r j=1 t j (1 -t 2 j )∂ t j F (t)D(e j , ēj ).
To put I +II in a better form we write t j (1-t

2 j ) = (1-t 2 j )( 1-t 2 j t j - (1-t 2 j ) 2 t j
) and we get then the form H j as claimed.

Note that this formula is consistent with the formula for Laplace-Beltrami operator on line bundle with parameter ν where the root multiplicity 1 for the root γ j is replaced by 1 + 2ν and the multiplicity 2b of γ j 2 by 2b -2ν; see e.g. [START_REF] Shimeno | The Plancherel formula for the spherical functions with onedimensional K-type on a simply connected simple Lie group of hermitian type[END_REF].

We prove now the sufficiency of the Hua equation (30) in Theorem 5.2. Let ϕ λs,ν the unique elementary spherical function of type τ ν in E ν ,

ϕ λs,ν (g) = K e (λs+ρ)H τ ν (k -1 κ(g -1 k))dk, g ∈ G
With some slightly abuse of notation we denote ϕ λs,ν (z) the corresponding K-invariant function on Ω via our trivialization.

Suppose F be a K-invariant analytic eigenfunction of the Hua equation [START_REF] Zhang | Shimura invariant differential operators and their eigenvalues[END_REF]. Let g ∈ G. Then the function

Φ g (z) = K F (gkz)dk is K-invariant solution of the differential system H j Φ = 2 n r s( n r (s -1) + ν)Φ, j = 1, . . . , r.
Therefore by a result of Yan [START_REF] Yan | A class of generalized hypergeometric functions in several variables[END_REF], Φ g is proportional to ϕ λs,ν and thus, K F (gkz)dk = ϕ λs,ν (z)F (g • 0).

By Shimeno [24, Theorem 3.2] we see that this integral formula characterizes the joint eigenfunctions in E λ , that is F ∈ A(G/K, M λs,ν ). Now applying Theorem 4.1 we can find f ∈ B(G/P, L λs,ν ) such that F = P λs,ν (f ). The rest of the proof is the same as in [START_REF] Berline | Équations de Hua et noyau de Poisson. Noncommutative harmonic analysis and Lie groups[END_REF][START_REF] Koufany | Hua operators on bounded symmetric domains]Hua operators and Poisson transform for bounded symmetric domains[END_REF], and we get, under the condition [START_REF] Yan | A class of generalized hypergeometric functions in several variables[END_REF], that f is actually a function on S, i.e., f ∈ B(S, s, ν). We shall not repeat the computations.

Type one domains

Let V = M r,r+b (C) be the vector space of complex r×(r+b)-matrices. V is a Jordan triple system for the following triple product {xȳz} = xy * z + zy * x, and the endomorphisms D(z, v) are given by D(z, v)w = {zvw} = zv * w + wv * z. As a homogeneous space, the bounded domain I r,r+b can be identified with SU(r, r + b)/S(U(r) × U(r + b)).

The complex Lie algebra k C is given by the set of all matrices

a 0 0 d , a ∈ M r,r (C), d ∈ M r+b,r+b (C), tr(a) + tr(d) = 0.
Hence, k C can be written as the sum

k C = k (1) C ⊕ k (2)
C , where k

C = { a 0 0 -tr(a) r+b I r+b , a ∈ M r,r (C)}, k (2)
C = { 0 0 0 d , d ∈ M r+b,r+b (C), tr(d) = 0}.
Let H (1) be the projection of the Hua operator H on k

(1)

C . It maps C ∞ (G/K, τ ν ) into C ∞ (G/K, τ ν ⊗ k (1) C ).
The main result of this section is the following theorem which is the analogue of Theorem 5.2 for non tube domains of type one. Theorem 8.1. Let Ω = I r,r+b be a bounded symmetric domain of type one. Suppose s ∈ C satisfies the condition [START_REF] Yan | A class of generalized hypergeometric functions in several variables[END_REF]. Then the Poisson transform P λs,ν is a SU(r, r + b)-isomorphism of B(S r,r+b ; s, ν) onto the space of analytic functions F on I r,r+b such that (34)

H (1) F = (r + b)s((r + b)(s -1) + ν)F ⊗ I r .
The bounded domain I r,r+b is of non tube type, however the characterization of Poisson integrals involves a Hua operator of the second order. As a consequence, the proof of Theorem 8.1 is similar to Theorem 5. Indeed, using the computations of the proof of Theorem 6.1 and some arguments in the proof of [START_REF] Koufany | Hua operators on bounded symmetric domains]Hua operators and Poisson transform for bounded symmetric domains[END_REF]Theorem 5.3] we can extend the formula [START_REF] Zhang | Invariant differential operators on hermitian symmetric spaces and their eigenvalues[END_REF] to any (not necessarily tube) bounded symmetric domain:

H 1 P s,ν (z, ū) = ( n r s) 2 D(B(z, z)(z z -ūz ), z z -u z) -( n r s)pZ 0 -( n r s
)νD(B(z, z)(ū z -ūz ), z z -u z ) P s,ν (z, ū)I r . Specifying this formula to the domain I r,r+b we get proposition.

On the other hand, to prove the sufficiency condition, the key point is to compute the radial part of the Hua operator H (1) . This follows immediately from the proof of Theorem 7.1.

Theorem 8.3. Let f be C 2 and K-invariant function, then for a = r j=1 t j e j , (35) 4H (1) f (a) = r j=1 H j F (t 1 , . . . , t r )D(e j , ēj ) (1) ,

where the scalar-valued operators H j are given by

H j =(1 -t 2 j ) 2 ∂ 2 ∂t 2 j + 1 t j (1 + 2ν) ∂ ∂t j + + k =j (1 -t 2 j )(1 -t 2 k ) 1 t j -t k ∂ ∂t j - ∂ ∂t k + 1 t j + t k ∂ ∂t j + ∂ ∂t k + +(2b -2ν)(1 -t 2 j ) 1 t j ∂ ∂t j
and D(e j , ēj ) (1) is the k

C component of D(e j , ēj ) ∈ k C .

Relative discrete series

The Poisson transform is not injective for singular s being in the set Λ 1 ∪Λ 2 in [START_REF] Yan | A class of generalized hypergeometric functions in several variables[END_REF]. It arises thus a question of understanding the image. A finite subset of such s corresponds to the relative discrete series, i.e. the images constitute discrete components in the decomposition of the L 2 space for the bundle. In this final section we find the precise parameters and compute explicitly the corresponding Poisson transforms on some distinguished functions, thus producing elements in the relative discrete series.

Fix ν > p -1 and ν ∈ pZ. Let α = ν -p > -1 and consider the weighted probability measure on Ω (36)

dµ α (z) = c α h(z, z) α dm(z)
The group G acts unitarily on L 2 (Ω, µ α ) via [START_REF] Johnson | The Hua operators on bounded symmetric domains of tube type[END_REF].

The irreducible decomposition of L 2 (Ω, µ α ) under the G-action has been given by Shimeno in [START_REF] Shimeno | The Plancherel formula for the spherical functions with onedimensional K-type on a simply connected simple Lie group of hermitian type[END_REF]Theorem 5.10] where he proved abstractly that all discrete parts called relative discrete series appearing in the decomposition are holomorphic discrete series. In this section we need their explicit realization given by the second author in [START_REF] Zhang | Nearly holomorphic functions and relative discrete series of weighted L 2 -spaces on bounded symmetric domains[END_REF].

Let us introduce the conical functions, see [START_REF] Faraut | Analysis on Symmetric Cones[END_REF]. Let c 1 , . . . , c r be the fixed Jordan frame. Put e j = c 1 + . . . + c j , for j = 1, . . . , r. Let

U j = {z ∈ V : D(c j , c j )z = z}. Then U j is a Jordan subalgebra of V 1 = V 1 (e)
with a determinant polynomial ∆ j . We extend the principal minors ∆ j to all V via ∆ j (z) := ∆ j (P U j (z)), where P U j is the orthogonal projection onto U j . Notice that ∆ r (z) = ∆(z) = det(z).

For any

m = (m 1 , . . . , m r ) ∈ C r , consider the associated conical function ∆ m (z) := ∆ m 1 -m 2 1 (z)∆ m 2 -m 3 2 (z) • • • ∆ mr r (z). If z = r j=1 z j c j then ∆ m (z) = r j=1 z m j j . Denote ℓ = α+1 2 -1 = ν-p-1 2 if α is an odd integer [ α+1 2 ] = [ ν-p+1 2
] otherwise here [t] stands for the integer part of t ∈ R. Define (37)

D ν = {m = r j=1 m j γ j , 0 ≤ m 1 ≤ • • • ≤ m r ≤ ℓ, m j ∈ Z}.
We let A 2,α m to be the subspace of L 2 (Ω, µ α ) generated by the function ∆m (q(z)), for m ∈ D ν , where q(z) = zz is quasi-inverse of z with respect to z.

We reformulate [24, Theorem 5.10] and [32, Theorem 4.7, remark 4.8] in the following.

Theorem 9.1. The relative discrete series representations appearing in L 2 (Ω, µ α ) are all holomorphic discrete series of the form A 2,α m with m ∈ D ν . The highest weight vector of A 2,α m is given by ∆m (q(z)).

We can now state the main theorem of this section.

Theorem 9.2. Let δ = (δ, δ, . . . , δ) such that s n r = n r + δ -ν and 0 ≤ δ < [ ν-p 2 ]. Then the Poisson transform P λs,ν is a G-equivariant Szgeö type map from the space B(S, s, ν) K of K-finite elements onto the K-finite elements in the relative discrete series A 2,α δ . This theorem is consequence of the following proposition. Finally since ∆(z) h(z,z) = ∆(q(z)), see [START_REF] Zhang | Nearly holomorphic functions and relative discrete series of weighted L 2 -spaces on bounded symmetric domains[END_REF]Corollary 4.4], the proof is completed.

Finally

  , let B(G/P ; L λ,ν ) be the space of hyperfunctions f on G satisfying f (gman) = e (λ-ρ)(log a) τ ν (m) -1 f (g) for all g ∈ G, m ∈ M, a ∈ A, n ∈ N + . The space B(G/P ; L λ,ν ) is a G-submodule of B(G) and can be identified with the space of hyperfunction valued sections of the line bundle L λ,ν on G/P associated with the character of P given by man → e (ρ-λ)(log a) τ ν (m), m ∈ M, a ∈ A, n ∈ N + . 4. Poisson transform 4.1. Poisson integrals on G/P for a minimal P . Let λ ∈ a * C and ν ∈ C. We define the Poisson transform P λ,ν by (12

Theorem 4 . 1 ( 2 , 3 , 4 . 2 .

 412342 [START_REF] Shimeno | Eigenspaces of invariant differential operators on a homogeneous line bundle on a Riemannian symmetric space[END_REF], Theorem 8.1). If λ ∈ a * C and ν ∈ C satisfy the conditions . . .} for any α ∈ Σ +[START_REF] Lassalle | Les équations de Hua d'un domaine borné symétrique du type tube[END_REF] e ν (λ) = 0,then the Poisson transform P λ,ν is a G-isomorphism of B(G/P, L λ,ν ) onto A(G/K, M λ,ν ).Poisson integrals on the Shilov boundary. Let h be the unique K-invariant polynomial on p + = V whose restriction on Re 1 ⊕ • • • ⊕ Re r is given by h(

  ). Since D(u, ū) = 2I and n α=1 D(e α , ēα ) = 2 n r I by [14, Lemma 5.1], the claim is proved for z = 0. Note furthermore that the Poisson kernel satisfies P s,ν (gz, gu) = J g (z) ν p P s,ν (gz, gu)J g (u) ν p .

= 1 4 (∂ 2 r + 1 rH

 11 It's convenient to find the radial part of 4H (due to the usual convention that ∂∂ = 1 4 (∂ 2x + ∂ 2 y ) ∂ r ) for radial functions in z ∈ C). Theorem 7.1. Let Ω be the tube type domain. Let F be C 2 and K-invariant function, then for a = r j=1 t j e j , j F (t 1 , . . . , t r )D(e j , ēj ),

Let

  I r,r+b = {z ∈ M r,r+b (C) : I r -z * z ≫ 0} where I r denote the unit matrix of rank r. Then I r,r+b is a bounded symmetric domain of dimension r(r + b), rank r and genus 2r + b. The multiplicities are 2b and a = 2 if 2 ≤ r, a = 0 if r = 1. The domain I r,r+b is of tube type if and only if b = 0. Its Shilov boundary is S r,r+b = {z ∈ M r,r+b (C) : z * z = I r }.

  2 and [14, Theorem 6.1]. The necessarily condition is guaranteed by the following Proposition 8.2. For any u ∈ S r,r+b the function z → P s,ν (z, ū) satisfies HP s,ν (z, ū) = (r + b)s[(r + b)(s -1) + ν]P s,ν (z, ū)I r .

Proposition 9 . 3 .

 93 Under the same conditions as in Theorem 9.[START_REF] Engliš | Covariant Cauchy-Riemann operators and higher Laplacians on Kähler manifolds[END_REF] we have(38) (P s,ν ∆δ )(z) = S P s,ν (z, u) ∆δ (u)du = (s n r ) δ ( n r ) δ ∆δ (q(z)).Let us first prove the formula below Lemma 9.4. For any highest weight m the following formula holds(39) ∆(z) ∆(w) δ K m (z, w) = c(m, δ) K m+δ (z, w) where c(m, δ) = ( n r + m) δ := r j=1 ( a 2 (r -j) + 1 + m j ) δ Proof. Since f (z) → ∆(z) δ f (z) is a K-intertwining map from P m onto P m+δ we have that ∆(z) δ ∆(w) δ K m (z, w) is equal to K m+δ (z,w) up to a constant c(m, δ). Now taking z = w = e and using [5, Lemma 3.1, Theorem 3.4] we find that the constant is c(m, δ) = K m (e, e) K m+δ (e, e) = d m (n/r) m (n/r) m+δ d m+δ = ( n r + m) δ .Proof of Proposition 9.3. Put σ = s n r . We compute the image P s,ν ∆δ ,(P s,ν ∆δ )(z) = S h(z, u) -ν h(z, z) σ h(z, u) -σ h(u, z) -σ ∆δ (u)du =h(z, z) σ S h(z, u) -ν-σ h(u, z) -σ ∆δ (u)duWe use the Faraut-Koranyi expansion [4, Theorem 3.8] of the reproducing kernels h(z, u) -ν-σ and h(u, z) -σ so that,(P s,ν ∆δ )(z) =h(z, z) σ S m≥0 (ν + σ) m K m (z, u)× × m ′ ≥0 (σ) m ′ K m ′ (u, z) ∆δ (u)duwhich by Lemma 9.4 is(P s,ν ∆δ )(z) =h(z, z) σ m≥0 (ν + σ) m K m+δ (z, u)∆ δ (z) -1 c(m, δ)× × m ′ ≥0 (σ) m ′ K m ′ (u, z) du =h(z, z) σ m≥0 (ν + σ) m (σ) m+δ c(m, δ)∆ δ (z) -1 × × K m+δ (z, u)K m+δ (u, z)duThe last equality follows from the Schur orthogonality relation. Furthermore, since the ratio of the Fischer inner product and the standard K-invariant inner product of L 2 (S) is constant, see[START_REF] Faraut | Analysis on Symmetric Cones[END_REF] Corollary 3.5] or[START_REF] Upmeier | Toeplitz operators on bounded symmetric domains[END_REF], we obtain(P s,ν ∆δ )(z) =h(z, z) σ m≥0 (ν + σ) m (σ) m+δ c(m, δ)∆ δ (z) m+δ K m+δ (z,z) which again by Lemma 9.4 is (P s,ν ∆δ )(z) =h(z, z) σ ∆δ (z) m≥0 (ν + σ) m (σ) m+δ ( n r ) -1 m+δ × × ∆δ (z)K m (z, z) =h(z, z) σ (σ) δ ( n r ) δ ∆δ (z) m≥0 (σ + δ) m K m (z, z) =h(z, z) σ (σ) δ ( n r ) δ ∆δ (z)h(z, z) -(σ+δ)

Research by G. Zhang supported by Swedish Science Council (VR). .


Rémi Agier

Sébastien Valette

Laurent Fanton

Pierre Croisille

Rémy Prost

Hubless 3D medical image bundle registration

Keywords: Big Data, Medical Imaging, Points of Interest, Registration

We propose a hubless medical image registration scheme, able to conjointly register massive amounts of images. Exploiting 3D points of interest combined with global optimization, our algorithm allows partial matches, does not need any prior information (full body image as a central patient model) and exhibits very good robustness by exploiting inter-volume relationships. We show the efficiency of our approach with the rigid registration of 400 CT volumes, and we provide an eye-detection application as a first step to patient image anonymization.

Introduction

The increasing availability of digital medical imaging techniques such as Magnetic Resonance Imaging (MRI), Computed Tomography (CT) and Ultrasound (US) and the amount of data to process in healthcare networks has grown exponentially, thus illustrating the Big Data challenge in medecine. To perform early screening, monitoring or to follow-up thousands of patients throughout their healing, stress has to be put on the ability to perform fast and robust image analysis. This sums up to one question: how to process several thousand volumes automatically and robustly? In this context, we propose a versatile feature-based co-registration framework which can serve as a first step towards volume collections processing. The aim of this paper is to register the volumes altogether. this paper currently focus on robustness rather than accuracy, with the belief that the challenge is in inter-patient variability.

The Picture Archiving and Communication System (PACS) is at the center of image management in healthcare networks. Medical images are stored in the Digital imaging and communications in medicine (DICOM) format which contains, in addition to images, informations such as acquisition parameters and patient data. With cohort studies [START_REF] Bild | Multi-ethnic study of atherosclerosis: objectives and design[END_REF] and mass computation, spatial consistency between images is crucial. Despite the fact that DICOM contains data about spatial positions, we can see in figure 1 that they are not consistent. So, when one wants to deal with large medical image datasets, an initial Figure 1: Mass 3D image registration. (Top) a bundle of heterogeneous images arranged according to DICOM embeded data. This shows that DICOM metadata are generally not consistent and may hinder subsequent computations. (Bottom) After mass hubless registration, images are well aligned and the bundle may be used in further computation. Registration is robust over heterogeneous images, with, for example, only small body parts like heads. For visualization purposes, images are distributed over an horizontal grid.

global registration is needed. We propose to solve this global challenge using a hubless approach. To our knowledge, this is the first time that an approach which jointly process medical images is proposed. This paper is organized as follows : Section 2 gives a brief overview of related approaches both in Medical Imaging and Computer Vision fields. Section 3 presents our hubless registration proposal. Section 4 shows experimental results and a conclusion follows.

Related works

Points of interest

During the last decades points of interest [START_REF] Harris | A combined corner and edge detector[END_REF][START_REF] Lowe | Distinctive image features from scaleinvariant keypoints[END_REF] have successfully been exploited for tasks such as object recognition [START_REF] Lowe | Object recognition from local scaleinvariant features[END_REF][START_REF] Lowe | Distinctive image features from scaleinvariant keypoints[END_REF], action recognition, robotic navigation, panorama creation, etc... They aim at being fast while reducing the amount of data to process, mainly to deal with realtime processing or tasks involving large amounts of data. Their suitability to medical imaging has been evaluated in [START_REF] López | Evaluation of methods for ridge and valley detection[END_REF], and various applications have been proposed in this context, like image annotation [START_REF] Datta | Content-based image retrieval: approaches and trends of the new age[END_REF], image retrieval [START_REF] Zheng | Interest point based medical image retrieval[END_REF]. Most of the applications relate to image registration and matching. There are some 3D points of interest development, like [START_REF] Knopp | Hough transform and 3d surf for robust three dimensional classification[END_REF], which use rasterized meshes to describe 3D shapes.

Image registration

Image registration in medical application is a wide research field [START_REF] Hill | Medical image registration[END_REF][START_REF] Sotiras | Deformable medical image registration: A survey[END_REF] where a lot of different approaches have emerged. We split the approaches in two categories:

• dense -or voxelwise-registration, with minimization or maximization of an energy function (for example mutual information [START_REF] Pluim | Mutual-information-based registration of medical images: a survey[END_REF]).

The main advantage is to provide a dense registration, with an information about deformation at each point of the space. However, voxelwise computation is time-consuming.

• sparse registration, using points of interest [START_REF] Allaire | Full orientation invariance and improved feature selectivity of 3d sift with application to medical image analysis[END_REF][START_REF] Cheung | N-sift: Ndimensional scale invariant feature transform for matching medical images[END_REF][START_REF] Khaissidi | A fast medical image registration using feature points[END_REF]. This approach provides a less accurate method because points of interest do not always span the whole space. But if one has to perform multiple registrations, this approach is several orders of magnitude faster than dense registration, because points of interest have to be extracted only once per image.

Registration can be used directly for medical applications, as instance for therapy planning [START_REF] Pelizzari | Accurate threedimensional registration of ct, pet, and/or mr images of the brain[END_REF][START_REF] Rosenman | Image registration: an essential part of radiation therapy treatment planning[END_REF]. It can also be used as an essential step for other algorithms such as atlas-based approaches [START_REF] Gass | Multi-atlas segmentation and landmark localization in images with large field of view[END_REF].

Bundle optimization

As first applications in computer vision, points of interest are generally used in the medical field to match two images. But novel approaches in computer vision have appeared, using multiple images in order to tackle problems such as real time 3D reconstruction [START_REF] Triggs | Bundle adjustment -a modern synthesis[END_REF], efficient tracking using low-end cameras [START_REF] Karlsson | The vslam algorithm for robust localization and mapping[END_REF]. More recent works, like [START_REF] Frahm | Building rome on a cloudless day[END_REF], deal with large amount of data to reconstruct a town.

Bundle optimization is a promising paradigm, that paved the way to augmented reality and virtual reality [START_REF] Klein | Parallel tracking and mapping for small ar workspaces[END_REF], and can contribute to emerging challenges in medical image processing:

• Large amounts of images : nowadays, medical imaging is a very spread technology and more and more images are produced each day.

• Multiple modalities : in order to be more accurate, multiple modalities (CT, MRI, US, etc.) may be used in order to establish a diagnostic. Algorithms have to follow this trend and manage multiple modalities.

Note that some papers already have proposed the use of multiple medical images such as Multi-Atlas approaches [START_REF] Gass | Multi-atlas segmentation and landmark localization in images with large field of view[END_REF]. But these approaches most often carry out several applications of one-toone matching with a hub model [START_REF] Bartoli | Stratified generalized procrustes analysis[END_REF] and are limited to completely overlapping input images [START_REF] Marsland | A minimum description length objective function for groupwise non-rigid image registration[END_REF] contrary to image bundle optimization. We call these approach hub-based, in contrast with our hubless model (figure 2).

3 Hubless registration

Multiple views vs Multiple patients

There is a fundamental difference between scene reconstruction and body registration. More precisely, 3D reconstruction assumes a unique scene acquired with different points of view in contrast with multiple patients, variability breaks the hypothesis of a unique scene. The challenge shifts from 3D estimation with 2D data to inter-patient variability handling among 3D images.

The "hub-based model" issue

If one wants to register two images I 1 and I 2 , one containing an upper body part and one a lower body part, one faces the overlap limit (August and Kanade, The hubless approach : no image is picked as a model. This allows to compute 10 links between the 5 images.

2005

). One solution is to use a third image I re f , a full body image (the hub) and register it with I 1 and I 2 . Afterwards, the two registrations are composed in order to obtain the upper body to lower body registration.

In the case of image bundles, an image used as a hub to register two images may not be the best hub for other images. As a result, picking the right hub is a complex task. Moreover, this "one size fits all" model may be difficult to register with some images due to patient anatomy variability.

Challenges and Contributions

Our approach of registration is to deal with the whole image bundle at once, which brings major improvements compared to one-to-one registration. This paper only deals with rigid registration, as a proof of concept, with the intend to generalize it to non-rigid registration in future works. A first challenge is to deal with a big amount of data. In this paper, each image is a scanner acquisition which resolution ranges between 200 3 and 600 3 , each voxel being encoded with 16 bits. In the last experimentation, we have a bundle of 400 images representing more than 100GB of image data. We overcome this challenge by converting the whole image bundle into a compact representation such as Speeded Up Robust Features (SURF [START_REF] Bay | Surf: Speeded up robust features[END_REF])) that we extend to 3D in section 3.4.

Moreover, we need to be able to process images containing not only full-body scans but also body parts which is usually a problem. Indeed, incomplete data (body parts) generally hinder the registration process, due to the overlap limitation. In our case, we overcome this difficulty, and use partial matches to improve our results.

Also, as seen on figure 2, registering a group of n images with a hub-based approach, is performed with n registrations. On the other hand, with a hubless approach, one can benefit from a much higher number of registrations. This number can be as high as n(n -1)/2, depending on the overlap between images. The advantage of this difference is the problem which becomes more and more overdetermined with the number of images. We exploit this fact to increase robustness and accuracy of our approach. A second challenge is to face the overdetemined nature of our problem, for which we propose a novel solution in section 3.6.

As a consequence, partially overlapping image sets can easily be processed with our approach, and incorrect registrations do not significantly impact accuracy as long as the problem remains overdetermined.

3D surf

The SURF approach comes as a fast and efficient points of interest extractor, but was originally created for 2D images. We then developed a generalization of the SURF descriptor to deal with 3D medical images.

Blob-like structure detection

The first step is to extend the 3D scale-space of a 2D image into a 4D function of a 3D image. In spirit with 2D SURF, we compute a box-filter approximate Hessian matrix H(x) at each point x = (x, y, z) of the image. This results in a 3x3 Matrix H. Blob detection is carried out by analysing the sign of the eigenvalues of H. But in contrast with 2D SURF, computing the determinant of H is not sufficient to check that its eigenvalues are all negative or positive. In spirit with what Allaire et al. proposed for 3D SIFT [START_REF] Allaire | Full orientation invariance and improved feature selectivity of 3d sift with application to medical image analysis[END_REF], we use the trace of the Hessian tr(H) and the sum of principal second-order minors ∑ det P 2 (H) in addition to the determinant to have sufficient knowledge on the eigenvalues.

Description

For a 2D image, SURF splits the neighborhood of each point of interest into smaller 4x4 square subregions. For each sub-region, a set of features is computed using Haar wavelet responses (2 responses per direction). All responses are concatenated into a 64element vector which is normalized to be contrastinvariant. We apply the same method but take benefit of the third dimension, by splitting neighborhoods into 2x2x2 cubic-sub-regions, and extracting 3x2 responses. We build a 48-element vector, which is smaller compared to 2D SURF but contains more local information thanks to the third dimension.

Upright SURF

We chose to not extend the rotation-invariance of 2D SURF to our case, as with most problems subject to the curse of dimensionality, giving the descriptor invariance to rotations would decrease robustness. Various invariant solutions exist, as in [START_REF] Allaire | Full orientation invariance and improved feature selectivity of 3d sift with application to medical image analysis[END_REF][START_REF] Cheung | N-sift: Ndimensional scale invariant feature transform for matching medical images[END_REF]. But in our current application, we have the prior information that during image acquisition, the patient is standing, which make the principal orientation computing unnecessary as long as inter-patient deformation remains reasonable.

Pairwise registration

Once SURF descriptors are extracted from the whole image bundle, we register all possible image pairs using RANdom SAmple Consensus [START_REF] Fischler | Random sample consensus: a paradigm for model fitting with applications to image analysis and automated cartography[END_REF]. We currently deal with a rigid transform model with 4 degrees of freedom : 3 translations and 1 isotropic scale.

f : x → s.x + t (1)
The output of this step is a set of n * (n -1)/2 transforms which link the images altogether. However this set contains incorrect transforms due to several issues:

• non-overlapping images pairs.

• pairs with too few matches due to variability between patients.

Y =   1.5 -1 -1 -0.8 -1 1.8   t 3 2 = [-1, -0.8] t 1 3 = [-1, 1.8] t 2 1 = [1.5, -1] Y =   2 -1 -1 -1 -1 2   n 1 n 2 n 3

Hubless bundle registration

Problem Statement

Once we have the set of transforms, we want to be able to consistently register the images. Note that we restrict ourselves to using transforms with 3 degrees of freedom (translations). Extension to 4 degrees of freedom is explained in a subsequent section. In a slightly more abstract formulation, our problem is equivalent to compute point positions given only relative positions between them. Figure 4 shows an example where 3 nodes n 1 , n 2 and n 3 , have to be located given their 3 * (3 -1)/2 = 3 relative positions t 2 1 , t 3 2 and t 1 3 .

Solving Laplacian equations

We propose to solve this problem by writing it as a Laplace equation [START_REF] Cohen-Or | Encoding meshes in differential coordinates[END_REF], using the image bundle complete graph. The graph carries several kinds of information:

• Each node n i corresponds to one patient image I i , with a local reference frame x i which will be adjusted consistently with the bundle.

• Each edge e k j carries the translation t k j between the two images I j and I k as computed using pairwise registration (section 3.3). Note that due to the presence of incorrect registrations (section 3.3), the set of translations is not always consistent. The problem now simplifies to finding the best set of local frames x i given the set of translations t k j . In graph theory, there is an analogy which can be done in relation with relative / absolute information, based on a matrix representation of the graph:

x 1 = [0, 2] x 2 = [1.67, 1] x 3 = [0.83, 0.2] P =   0 2 1.67 1 0.83 0.2   P =   0 2 2 1 1 0   n 1 n 2 n 3
E t .P = Y (2)
where E is the incidence matrix that shows the relationship between vertices. P is the local frame matrix which stores local frames x i , each row representing one vertex, each column referring to one space dimension. Y refers to the observation matrix, where each row contains the translation carried by one edge e k j . In our case, we aim at finding the best estimate of P given an inconsistent observation matrix Y, which we propose to solve using a Laplacian equation. Given a graph, the matrix L = E t .E is known as the graph laplacian operator matrix. This matrix is symmetric, singular and positive semi-definite.

By left-multiplying equation ( 2), by E the Laplacian matrix appear as follow :

L t .P = E.Y (3)
which cannot be solved as it is because the Laplacian matrix is always singular, of rank n -1. In a geometric point-of-view, this result in the fact that absolute reconstruction from relative positions is defined up to an offset. Applying the same displacement to all local frames doesn't change relative positions between them. To make the Laplacian matrix non singu-Figure 6: Example of problem with 3 volumes (legs, head, body) and 1 impossible registration. To successfully register this bundle of 3 volumes, the registration between the legs and head has to be discarded. lar, we fix the absolute coordinates of one vertex, referred to as anchor in (Cohen-Or and Sorkine, 2006). We solve equation 3 using least square minimization. Figure 5 shows the solution to the problem depicted in figure 4 with an anchor fixed on vertex n 1 .

Graph decimation

All possible image pairs are used to perform the registration and construct the graph. As a consequence, non-overlapping body parts create outliers in the graph, as shown in figure 6. In other word, some edges represent a wrong observation. As a reminder, the breakdown of least square approaches is 0 [START_REF] Meer | Robust regression methods for computer vision: A review[END_REF]. Then in our case, the error caused by one outlier edge will not be filtered out but it will be distributed over all the graph nodes. One solution is to discard edges representing the worst registrations, while keeping the graph connected. This requires the definition of a quality criterion for each individual transform. For robustness purpose, given two images and their registration computed using RANSAC, we simply use the number of inliers found by RANSAC. We then sequentially remove edges according to this criterion, avoiding the removal of any edge that would break the graph connectedness. We chose to remove edges until the graph contains k * n edges. We experimentally set k to 3.

Extension to scale computation

Our graph approach can only deal with summable values, allowing us to manage only translations between images.

But scale information between two images depicts a multiplicative relationship, as in Equation ( 1). We address this problem in a separable way, with, on one hand, the translation parameters and on the other hand the scale parameters. Taking the logarithm of equation (1) changes the scale into an additive value.

We exploit this fact by processing log-scale as a fourth dimension, supplementary to the three first coordinates.

Results

We applied this method to register a bundle of heterogeneous Computed Tomography (CT) 3D images. Images exhibit different dimensions, different resolutions and most importantly may contain different body parts. All computation have been carried out on a 24-core workstation with 128 GB of RAM, using the DESK framework [START_REF] Jacinto | A web interface for 3d visualization and interactive segmentation of medical images[END_REF].

In practice, SURF points of interest provide a very compact representation of the data. As an example, one image, which weighs from 200 to 600 MB, can be turned into a 0.1 to 1MB SURF description. Complexity of this part is linear with the number of images and can be easily parallelized. For a typical full body acquisition (1.8*0.6*0.6 meters, resized into an isotropic image with 1.5mm spacing), about 5000 points of interest are extracted in about 2 minutes.

Hubless bundle registration exhibits quadratic complexity, as we have to register each possible image pair. During RANSAC registration, we perform 6000 iterations with a distance threshold of 40mm. As our metric reflects natural units, we add a criterion about difference in scale during match computations: if the scale ratio between two points of interest is larger than 1.5, we forbid the match between these points. Thanks to a very compact representation, computing a transformation between two typical full bodies (with full overlap) is done in a few seconds.

Points of interest approaches were originally used to register two images with different point of views of the same scene. A main criterion to evaluate performance of points extraction is repeatability. In our case, repeatability is less significant because we register images from different patients. Moreover we currently only deal with global translations between images, as we currently focus on robustness rather than accuracy. However we experimentally checked that images containing overlapping body parts are well registered together.

Figure 1 shows the registration of a bundle of 35 CT images, and figure 8 an example with 400 CT images. Our approach is able to register both bundles in a very robust way. Figure 7 shows processing time for a bundle depending on the number of images it contains. This graph exhibits a slight quadratic behavior, as computation time is still dominated by the points of interest extraction stage, which has linear complexity. Quadratic complexity will be more visible with bundles containing more than 400 images.

A first application of our algorithm is automatic eye covering for anonymization. The only needed user input is the head bounding box b head in one of the images, and the eyes bounding box b eyes . Afterwards, we perform these steps:

1. compute the hubless bundle registration 2. extract images containing heads, by propagating a human-made cut out of a head in one image, in spirit with atlas-based approaches 3. apply a second pass of the whole hubless process, using only points of interest located in the head.

4. transport the location of the pre-positioned eyes on every images containing eyes using the previous registration.

Re-computation using only partial set of points of interest allows the algorithm to be more accurate. Conceptually, this can be considered as a first step towards deformable registration, with a locally-rigid registration. First results can be seen in figure 9. One can note that the eyes are correctly covered except for cases where the patient head is tilted. Solving these cases will imply the use of transformation models with more degrees of freedom than simple translations.

Conclusion and Perspectives

We have shown the benefits of computing joint registrations compared to one-to-one registrations. This brings robustness in the case of heterogeneous medical datasets which may contain disjoint body parts. Compact approaches, using point of interest extraction, allow us to deal with large datasets, in reasonable time. Our approach provides an efficient way to background screening of large medical databases, using a simple translation-and-scale transform model. Future work may study the extension of this approach to non-rigid transforms. We also plan to tackle the complexity issues that will arise when dealing with even bigger image bundles. Finally, our image graph is currently connected, but sub-graph extraction is a promising way to classify the patients, analyzing both the graph inner properties and patient data.

Figure 2 :

 2 Figure2: Fundamental differences between hub-based and hubless approaches. Here, 5 images are represented by 5 red dots. (Left) Hub-based approach : one needs an sixth image as a central model. This approach allows only 5 links between images, and its efficiency depends on the choice of the model image. (Right) The hubless approach : no image is picked as a model. This allows to compute 10 links between the 5 images.

Figure 3 :

 3 Figure 3: Points of interest extracted from a CT volume. Clearly, important structures inside the patient such as lungs are well detected. Blue and red circles represent dark and light blobs, respectively.

Figure 4 :

 4 Figure 4: A simple problem with 3 volumes (green dots) and 3 translations: Blue arrows depict observed registrations, orange dotted lines represent registration error. Matrix Y depicts ground truth transforms between all images. Matrix Y represents estimated transforms.

Figure 5 :

 5 Figure 5: Reconstructed positions. Using our algorithm on the problem shown in figure 4, we can reconstruct the positions. Blue arrows represent the corrected registrations between all images. Matrix P depicts ground the truth positions. Matrix P represents the estimated positions.

Figure 7 :

 7 Figure 7: Processing time vs number of volumes

Figure 8 :

 8 Figure 8: Mass registration of 400 volumes. (Top Image) A bundle of volumes displayed with raw information about position contained in DICOM. (Bottom Image) Registered bundle.

Figure 9 :

 9 Figure 9: Covering eyes on meshes extracted from CT-images. Note that here the surface meshes are solely used for visualization purposes, only the image is used for processing.