Andrea Bressan

Stefan Takacs
email: austriastefan.takacs@ricam.oeaw.ac.at

Sum-factorization techniques in Isogeometric Analysis

come

Introduction

Isogeometric Analysis, [START_REF] Hughes | Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement, Comput[END_REF], was proposed around a decade ago as a new approach for the discretization of partial differential equations (PDEs) and has gained much interest since than. Spline spaces, such as spaces spanned by tensor product B-splines or NURBS, are typically used for geometry representation in standard CAD systems. In Isogeometric Analysis, one uses such spaces for the geometry representation of the computational domain and as space of ansatz functions for solution of the PDE.

The fast assembling of stiffness and mass matrices is a key issue in Isogeometric Analysis, particularly if the spline degree is increased. If the assembling is done in a naive way, and p d quadrature points are used per element, the computational complexity of assembling a standard mass or stiffness matrix has order N p 3d , where N is the number of unknowns, p is the spline order (degree + 1), d is the domain dimension. In recent years, much effort was set on proposing Low rank assembling is based on a completely different idea. It is observed that for practical problems, evaluation of the geometry function on the quadrature points does not have the maximum possible rank. Note that often the geometry is already represented exactly on coarse grid levels and the grid is only refined for a better resolution of the solution of the PDE. In this case, the rank of the geometry function is unchanged during refinement and is small compared to the maximum possible rank. It was observed that low rank geometry functions yield to mass and stiffness matrices which typically have a small tensor rank. Mantzaflaris, Jüttler, Khoromskij, and Langer [START_REF] Mantzaflaris | Low rank tensor methods in Galerkin-based isogeometric analysis[END_REF] have discussed how to set up an assembling algorithm based on this approach. Hofreither [START_REF] Hofreither | A black-box low-rank approximation algorithm for fast matrix assembly in isogeometric analysis[END_REF] has shown that by rewriting the problem accordingly, standard adaptive cross approximation algorithms can be used as black-box methods to determine the mass or stiffness matrices. For these approaches, the overall computational complexity has order N Rp d , where R is the (unknown) tensor-rank of the resulting matrix. If the coefficients of the matrix are not actually needed, it is possible to compute a low-tensor rank representation of the matrix and to compute the corresponding matrix-vector products with a computational complexity of order N Rp.

In this note, we want to have again a look onto the approach of sum factorization. We will show that by performing a global variant of sum factorization, the computational complexity is reduced to order N p d+2 . Still, this algorithm yields exactly the same matrix which would be obtained with straight-forward assembling and the algorithm is as general and as extendable to other differential equation as straight-forward assembling. We will then see that the computational complexity is still preserved for a localized approach. We will observe that also here, like for the low rank assembling, the formation of the matrix is the most expensive part and that the corresponding matrix-vector products has a computational complexity of order N p d+1 . This article is organized as follows. In Section 2, we state the abstract formulation of the problem and give examples of bilinear forms falling into the class. We introduce a fast assembling procedure in Section 3. In Section 4, we discuss the localization of sum factorization and its applications. Then, in Section 5, we give an algebraic description of the proposed algorithm. In Section 6, we explain why the evaluation of the basis functions and the geometry function is non trivial. We comment on matrix-free implementations in Section 7. Finally, in Section 8 we give numerical experiments and draw conclusions.

Preliminaries

In the following L 2 , L ∞ and H 1 denote the standard Lebesgue and Sobolev spaces with standard scalar products. Consider a bilinear form of the type

a : L 2 ([0, 1] d ) × L 2 ([0, 1] d ) → R, a(u, v) := (Fu, v) L 2 ([0,1] d ) = [0,1] d F(x)u(x)v(x) dx, (1)
where d ∈ N is the domain dimension and F ∈ L ∞ ([0, 1] d ) is a given coefficient function. At a first glance, this model problem looks rather restrictive as it is only an L 2 -like scalar product. However, we will see in Section 2.1 that we can write the variational form of any PDE, like Poisson, or system, like Stokes, as a sum of bilinear operators of the form [START_REF] Ainsworth | Bernstein-Bézier finite elements of arbitrary order and optimal assembly procedures[END_REF].

The goal of this article is to formulate fast algorithms that compute the matrix corresponding to a Petrov-Galerkin discretization of a. (This, of course, also covers the special case of standard Galerkin discretizations.) Let Φ = (φ n ) N n=1 , Ψ = (ψ m ) M m=1 be two generating systems of functions that span corresponding subspaces of L 2 ([0, 1] d ). Following the Petrov-Galerkin approach, the restriction of a to span Φ × span Ψ is represented by the matrix

A := [a(φ n , ψ m )] m=1,...,M n=1,...,N . (2)
In practice, the coefficients of A are not computed exactly, but they are approximated using a quadrature rule. This means that A is approximated by

A := [a(φ n , ψ m )] m=1,...,M n=1,...,N , (3)
with a(φ n , ψ m ) := x∈X ω(x) F(x) φ n (x) ψ m (x), ( 4
)
where X is the set of quadrature points and ω(x) are the corresponding quadrature weights. The techniques in this note require a tensor-product structure.

• Tensor-product discretization: For each direction δ = 1, . . . , d, there exist Φ δ , Ψ δ such that

Φ = Φ 1 ⊗ . . . ⊗ Φ d , Φ δ = (φ =δ n ) N δ n=1 , Ψ = Ψ 1 ⊗ . . . ⊗ Ψ d , Ψ δ = (ψ =δ m ) M δ m=1 , (5)
• Tensor-product quadrature: For each direction δ = 1, . . . , d, there exist X δ and ω δ such that

X = X 1 × . . . × X d , ω(x) = ω 1 (x 1 ) . . . ω d (x d ). (6)
The meaning of ( 5) is that the d-variate functions in Φ are the products of the univariate functions in Φ 1 , . . . , Φ d (and analogously the functions in Ψ). By convention, we assume the lexicographic ordering. For a rigorous definition, let π : {1, . . . , N } → {1, . . . , N 1 } × • • • × {1, . . . , N d } be a lexicographic ordering, i.e., the bijection defined by

π := (π 1 , . . . , π d ), π δ (n) := (n -1) mod N ≤δ N ≤δ-1 + 1, N ≤δ := δ i=1 N i . Define σ : {1, . . . , M } → {1, . . . , M 1 } × • • • × {1, . . . , M d } analogously.
Using these orderings, we denote the generating functions for Φ and Ψ as follows:

φ n (x 1 , . . . , x d ) = φ =1 π1(n) (x 1 ) • • • φ =d π d (n) (x d ), n = 1, . . . , N, ψ m (x 1 , . . . , x d ) = ψ =1 σ1(m) (x 1 ) • • • ψ =d σ d (m) (x d ), m = 1, . . . , M.
Analogous to N ≤δ , it is convenient to define

M ≤δ := δ i=1 M i and X ≤δ := δ ą i=1 X i .
Moreover, for the complexity analysis, we need some bounds on the supports of the functions in Φ δ and Ψ δ . As usual in the context of Sobolev spaces, supports are always closed, i.e., we define supp f := {x : f (x) = 0}.

• Bounded number of active functions: for all δ ∈ {1, . . . , d} and x δ ∈ X δ :

#{n : x δ ∈ supp φ =δ n } ≤ p δ , #{m : x δ ∈ supp ψ =δ m } ≤ p δ . (7)
• Bounded overlap: for all δ ∈ {1, . . . , d} and all φ =δ n ∈ Φ δ :

#{m : supp φ =δ n ∩ supp ψ =δ m ∩ X δ = ∅} ≤ k δ p δ . (8)
For the convenience of the reader, we will often express the main results in simplified form using p := max{p 1 , . . . , p d }.

Remark 1. The condition [START_REF] Giannelli | THB-splines: the truncated basis for hierarchical splines[END_REF] bounds the number of test functions that interact with any of the trial functions; the presented complexity analysis could also be done if the roles of the trial functions and the test functions as well as of N and M are interchanged.

To keep the notation tight, here and in what follows, a b means that there is a constant c > 0 (which does not depend on any of the quantities discussed in this paper) such that a ≤ c b and a b means that a b and b a.

Bilinear forms that depend on derivatives

The described algorithms are not limited to the computation of the mass matrix, but can be employed for Isogeometric Analysis of linear and non-linear PDEs or systems of PDEs. In this subsection we show how an assembling procedure for Petrov-Galerkin discretization of bilinear form of (1) extends naturally to bilinear forms that depend on the derivatives of their arguments.

Any bilinear form a :

H r ([0, 1] d ) × H r ([0, 1] d ) → R can be written as: a(φ, ψ) = θ,η∈Θr (F θ,η ∂ θ φ, ∂ η ψ) L 2 =: θ,η∈Θr a θ,η (∂ θ φ, ∂ η ψ), (9)
where Θ r ⊂ N d is the set of multiindices θ = (θ 1 , . . . , θ d ) corresponding to partial derivatives of order up to r. Analogously, the matrix A representing its restriction to Φ × Ψ decomposes as

A = θ,η∈Θr A θ,η ,
where A θ,η is the matrix representing the restriction of a θ,η to

∂ θ Φ × ∂ η Ψ with ∂ θ Φ := ( ∂ θ 1 ∂x θ 1 • • • ∂ θ d ∂x θ d φ n ) N n=1 , ∂ η Ψ := ( ∂ η 1 ∂x η 1 • • • ∂ η d ∂x η d ψ m ) M m=1 .
To extend any assembling procedure from a as in (1) to a as in [START_REF] Hofreither | A black-box low-rank approximation algorithm for fast matrix assembly in isogeometric analysis[END_REF], it is sufficient to prove that the assumptions on Φ and Ψ hold for the associated spaces ∂ θ Φ and ∂ η Ψ so that the procedure can be applied on each a θ,η . This is done in the next lemma.

Lemma 1. If Φ, Ψ ⊂ H r ([0, 1] d
) satisfy any of the conditions (5), ( 7), (8) then also ∂ η Φ, ∂ θ Ψ satisfies that condition for all θ, η ∈ Θ r .

Proof. First of all the tensor product structure carries over because

∂ θ Φ = ( ∂ θ 1 ∂x θ 1 • • • ∂ θ d ∂x θ d φ n ) N n=1 = ∂ θ1 Φ 1 ⊗ • • • ⊗ ∂ θ d Φ d ,
where

∂ θ δ Φ δ := ( ∂ θ δ ∂x θ δ φ =δ n ) N δ n=1 .
The same applies to ∂ η Ψ. For a sufficiently smooth function f : [0, 1] → R and using that the supports are closed, we have supp ∂ ∂x f ⊆ supp f . In particular, we deduce supp

∂ θ δ φ =δ ⊆ supp φ =δ , supp ∂ η δ ψ =δ ⊆ supp ψ =δ ,
so that ( 7) and ( 8) carry over to ∂ η Φ, ∂ θ Ψ with the same values of p δ and k δ .

Remark 2. Note that considering a Petrov-Galerkin approach (different spaces for test and trial functions) is required for the application of the algorithms to higher order problems even if Φ = Ψ because for any θ = η, typically the two spaces ∂ θ Φ and ∂ η Ψ differ.

Examples: B-spline based discretizations

Lemma 2. If Φ δ := (φ =δ n ) N δ n=1 and Ψ δ := (ψ =δ m ) M δ m=1 are univariate B-spline bases of order not more than p δ (degree not more than p δ -1) and sharing the same knot vectors such that no knot coincides with a quadrature point in X δ , then [START_REF] Calabrò | Fast formation of isogeometric Galerkin matrices by weighted quadrature[END_REF] and (8) hold with k δ = 2.

The same holds if the knot vector of Φ δ is a refinement of the knot vector of Ψ δ , i.e., it contains additional knots or some knots have a larger multiplicity.

Proof. Let q Φ,δ and q Ψ,δ the orders of Φ δ and Ψ δ , respectively. Denote the knot vector of Φ δ by Ξ := (ξ -q Φ,δ +1 , . . . , ξ N δ +q Φ,δ -1 ). Note that supp φ =δ n = [ξ n , ξ n+q Φ,δ ]. Consider some quadrature point x ∈ (ξ j , ξ j+1 ) and observe that the supports only of the functions φ =δ j+1-q Φ,δ , . . . , φ =δ j include the quadrature point x, which are q Φ,δ functions. The same arguments hold for Ψ δ . This shows [START_REF] Calabrò | Fast formation of isogeometric Galerkin matrices by weighted quadrature[END_REF]. Observe that the supports only of the functions ψ =δ n-q Ψ,δ +1 , . . . , ψ =δ n+q Φ,δ contribute to (ξ n , ξ n+q Φ,δ ) ⊇ supp φ =δ n ∩ X δ . This shows condition [START_REF] Giannelli | THB-splines: the truncated basis for hierarchical splines[END_REF] with k δ = 2. If now the functions in Φ δ are based on a refinement of Ξ, we know that the supports of the functions in Φ δ shrink. More precisely, we have supp φ =δ n ∩ X δ ⊆ (ξ n , ξ n +q Φ,δ ) for some accordingly chosen n . So, (8) stays true. Remark 3. If some quadrature points coincide with a knot, we obtain an analogous result with k δ = 3 if only B-splines with orders up to p δ -1 are considered.

Example 1 (Convection diffusion equation). We consider a standard singlepatch isogeometric discretization, so we assume that the computational domain Ω is parameterized by a diffeomorphism

G : Ω := [0, 1] d → Ω := G( Ω) ⊂ R d .
We assume that a source function f ∈ L 2 (Ω) and coefficient-functions

A ∈ L ∞ (Ω, R d×d ), b ∈ L ∞ (Ω, R d ) and c ∈ L ∞ (Ω) are given. The boundary value problem reads as follows. Find u ∈ H 1 (Ω) such that -∇ • (A∇u) + b • ∇u + cu = f in Ω, ∂u ∂n = 0 on ∂Ω.
The variational formulation reads as follows. Find u ∈ H 1 (Ω) such that

a(u, v) := c 0 b A u ∇u , v ∇v
L2(Ω) = (f, v) L2(Ω) for all v ∈ H 1 (Ω).
A standard isogeometric discretization is set up on the parameter domain, i.e., we first define the spline space

V h := span{ B 1 , . . . , B N } := d δ=1 span{ B (δ) 1 , . . . , B (δ)
N δ } ⊂ H 1 ( Ω),
where

B (δ)
n are the standard B-spline basis functions as given by the Cox-de Boor formula. Then, the ansatz functions are transferred to the physical domain using the pull-back principle

V h := V h • G -1 , B n := B n • G -1 .
For the computation of the stiffness matrix A, we transfer the functions of interest to the parameter domain and obtain

A = [a(B n , B m )] m=1,...,N n=1...,N = |J G | I J - G c 0 b A I J -1 G F := 1 ∇ D := B n , 1 ∇ B m L2( Ω) m n = d θ=0 d η=0 (F θ,η D θ B n , D η B m ) L2( Ω) a θ,η (D θ B n , D η B m ) := m=1,...,N n=1,...,N , (10)
where J G is the Jacobi matrix of the geometry function, |J G | the absolute value of its determinant and

D α := 1 for α = 0 ∂ ∂xα for α ∈ {1, . . . , d}.
Note that for each (θ, η), ∂ θ B n and ∂ η B m are tensor-products of B-spline basis functions and/or their derivatives. Using Lemmas 1 and 2, we obtain that (∂ θ B n , ∂ η B m ) satisfies ( 7) and [START_REF] Giannelli | THB-splines: the truncated basis for hierarchical splines[END_REF]. So, we immediately observe that the computation of A θ,η,h belongs to the abstract problem (1) -(2) with F := F θ,η .

Example 2 (Stokes equation). As a second example, consider the Stokes system, which reads in variational form as follows

(∇u, ∇w) L 2 (Ω) + (∇p, w) L 2 (Ω) = (f, w) L 2 (Ω) (∇ • u, q) L 2 (Ω) = 0 . (11)
A Galerkin discretization of such a system yields a block-matrix, which looks for two dimensions as follows

A = A B B 0 =   A 11 A 12 B 1 A 21 A 22 B 2 B 1 B 2 0   .
Each of the blocks A ij and B i is obtained by a Galerkin discretization of a scalar valued differential equation, which can be treated analogously to Example 1.

The discretization has to be chosen such that inf-sup stability holds; in [START_REF] Bressan | Isogeometric discretizations of the Stokes problem: stability analysis by the macroelement technique[END_REF], inf-sup stability has been shown for the isogeometric Taylor-Hood and the isogeometric sub-grid method.

In the isogeometric Taylor-Hood method, the pressure variable is discretized with B-splines of some order p with maximum smoothness. Each of the components of the velocity variable is discretized with B-splines of order p + 1 and reduced smoothness, p -2, obtained by taking each knot in the knot vector twice.

In the isogeometric sub-grid method, the pressure variable is again discretized with B-splines of some order p with maximum smoothness. Each of the components of the velocity variable is discretized with B-splines of order p + 1 with maximum smoothness on a grid obtained from the pressure grid by one uniform refinement step.

Here, the idea of sum factorization can be employed for each of the blocks of the overall system, where due to symmetry it suffices to discretize the Laplace operator A and the divergence operator B. The conditions of Lemmas 1 and 2 hold for A because Φ = Ψ, i.e., the spaces of the trial and the test functions coincide; they hold for B because the knot vector of the trial function space Φ is a refinement of that of the space of test functions Ψ.

Sum factorization

The method of sum factorization is based on the observation that expanding each term in (4) with respect to its components yields

a(φ n , ψ m ) = x1∈X1 • • • x d ∈X d d δ=1 ω δ (x δ ) d δ=1 φ =δ π δ (n) (x δ ) d δ=1 ψ =δ σ δ (m) (x δ ) F(x) = x d ∈X d φ =d π d (n) (x d ) ψ =d σ d (m) (x d )ω d (x d ) x1∈X1 • • • x d-1 ∈X d-1 d-1 δ=1 ω δ (x δ ) d-1 δ=1 φ =δ π δ (n) (x δ ) d-1 δ=1 ψ =δ σ δ (m) (x δ ) F(x) =: a ≤d-1 x d (φ ≤d-1 n , ψ ≤d-1 m ) , that is a(φ n , ψ m ) = x d ∈X d φ =d π d (n) (x d ) ψ =d σ d (m) (x d )ω d (x d ) a ≤d-1 x d (φ ≤d-1 n , ψ ≤d-1 m ), (12)
where

φ ≤δ n := δ i=1 φ =i πi(n)
and

ψ ≤δ m := δ i=1 φ =i σi(m) .
The term a ≤d-1

x d (φ ≤d-1 n , ψ ≤d-1 m
) is independent of the d-th components of Φ, Ψ, ω and it appears for many (φ n , ψ m ) pairs. The advantage of ( 12) over ( 4) is that it shows that a ≤d-1

x d (φ ≤d-1 n , ψ ≤d-1 m
) can be computed once and used many times. By rewriting a ≤d-1

x d (φ ≤d-1 n , ψ ≤d-1 m ) as a ≤d-1 x d (φ ≤d-1 m , ψ ≤d-1 n ) = x∈X ≤d-1 ω ≤d-1 (x)φ ≤d-1 n (x)ψ ≤d-1 m (x)F(x, x d ), (13)
we see that, analogously to (4), it is an approximation with the quadrature

(X ≤d-1 , ω ≤d-1 ) of the bilinear form on L 2 ([0, 1] d-1
) defined by

a ≤d-1 x d (φ, ψ) = [0,1] d-1 φ(x)ψ(x)F(x, x d ) dx.
Let A ≤d-1

x d be the matrix representing the Petrov-Galerkin restriction of a ≤d-1

x d to span Φ ≤d-1 × span Ψ ≤d-1 , i.e., A ≤d-1 x d = [a ≤d-1 x d (φ ≤d-1 n , ψ ≤d-1 m )] n=1,...,N ≤d-1 m=1,...,M ≤d-1 .
According to [START_REF] Mantzaflaris | Low rank tensor methods in Galerkin-based isogeometric analysis[END_REF], the entries of A are linear combinations of the entries in A ≤d-1

x d for different x d ∈ X d , but for the same φ ≤d-1 n and ψ ≤d-1 n
. This suggests the decomposition of A into blocks corresponding to the d-th components of φ n and ψ m , i.e. φ =d π d (n) and ψ =d σ d (m) . They identify a block of A = A ≤d of size M ≤d-1 × N ≤d-1 . Let B =d i,j denote these blocks as in

A ≤d =    B =d 1,1 . . . B =d 1,N d . . . . . . B =d M d ,1 . . . B =d M d ,N d    . (14)
Combining ( 14) and ( 12) we obtain the following recursive assembling procedure.

Algorithm 1 Recursive sum-factorization

1: procedure Assemble([X δ ] d δ=1 , [ω δ ] d δ=1 , [Φ δ ] d δ=1 , [Ψ δ ] d δ=1 , [F(x)] x∈X ) 2:
if d=0 then

3: return [F(x)] x∈X 4:
end if

5:

A ≤d ← 0 start with null matrix

6:
for all x d ∈ X d do sum all A ≤d-1

x d 7:
A ≤d-1

x d ← Assemble([X δ ] d-1 δ=1 , [ω δ ] d-1 δ=1 , [Φ δ ] d-1 δ=1 , [Ψ δ ] d-1 δ=1 , [F(x, x d )] x∈X ≤d-1 ) 8:
for all n d with φ =d

n d (x d ) = 0 do 9:
for all m d with ψ =d

m d (x d ) = 0 do 10: B =d m d ,n d ← B =d m d ,n d + ω d (x d )φ =d n d (x d )ψ =d m d (x d )A ≤d-1 x d cf. ( 14
)
11:

end for 12:

end for 13:

end for 14:

return A = A ≤d 15: end procedure For the complexity analysis, we assume that the values of F, φ n and ψ m are already available; cf. Section 6 on evaluating them and the related computational complexity. The recursive assembling procedure provides a recursive formula for the computational costs for assembling:

C ≤d #X d C ≤d-1 + p 2 d nnz(A ≤d-1 x d ) ,
where C ≤d-1 is the cost of assembling A ≤d-1

x d and nnz(A ≤d-1

x d
) is the number of nonzero entries of A ≤d-1

x d
. Here, the factor p 2 d comes from the number of iterations at line 8 and 9 that are bounded using [START_REF] Calabrò | Fast formation of isogeometric Galerkin matrices by weighted quadrature[END_REF] . Letting A ≤d-2

x d-1 ,x d be the analogous of A ≤d-1

x d and so on till A ≤0

x = F(x) and expanding the recursive cost formula yields

C ≤d d i=1 d δ=i+1 #X δ (p 2 i #X i ) nnz(A ≤i-1 xi,...,x d ). (15)
Note that (15) was obtained using (7) alone. The assumption (8) allow us to bound the number of non-zero entries as follows:

nnz(A ≤i-1 xi,...,x d ) ≤ i-1 δ=1 k δ p δ N δ , which yields C ≤d d i=1 p 2 i d δ=i #X δ i-1 δ=1 k δ p δ N δ ≤ d i=1 p 2 i #X d d-1 δ=1 max{#X δ , k δ p δ N δ }. (16)
Depending on the maximum, we obtain

∀δ, #X δ ≤ k δ p δ N δ ⇒ C ≤d d 2 d p d+2 N, (17) ∀δ, #X δ ≤ k δ N δ ⇒ C ≤d d 2 d p d+1 N, (18) ∀δ, #X δ ≥ k δ p δ N δ ⇒ C ≤d d p 2 #X, (19)
where we assume in all cases that k δ = 2. Remark 4. To achieve the above cost, it is necessary to implement Algorithm 1 properly. From the tensor product structure of the spaces, it follows that the sparsity pattern of A ≤d-1

x d
and that of the non-zero blocks B =δ m δ ,n δ are the same. By pre-computing the sparsity pattern, the sums at line 10 of the algorithm can be performed at a cost proportional to the number of non-zeros entries in A ≤d-1

x d .
Example 3 (Gauss quadrature). In Isogeometric Analysis, typically spline spaces are used for discretization (cf. Examples 1 and 2). The breakpoints of these spline spaces introduce a decomposition of the parameter domain into

K δ ≤ N δ -p δ + 1
elements on which the splines are polynomial. For assembling, typically Gauss quadrature with p δ quadrature nodes per element and direction is performed, so using (17) we obtain

#X δ = p δ K δ ≤ p δ N δ and C ≤d d 2 d p d+2 N.
Example 4 (Weighted quadrature). A way to reduce the computational costs is to reduce the number of quadrature points. In [START_REF] Calabrò | Fast formation of isogeometric Galerkin matrices by weighted quadrature[END_REF] it was shown, that this can be done without loosing accuracy using test-function dependent weights:

ω(x, m) = d δ=1 ω δ (x δ , σ δ (m)).
The adaptation of Algorithm 1 to test-function-dependent quadrature formulas is straightforward. Indeed, it can equivalently be thought as a quadrature with constant weights ω(x) = 1 and a different test space Ψ whose functions are pre-multiplied by the quadrature weights. It is interesting to note that the computational cost reported in [START_REF] Calabrò | Fast formation of isogeometric Galerkin matrices by weighted quadrature[END_REF] is explained by (18). Indeed with their construction #X δ N δ . Note however that this strategy breaks the symmetry between test and trial functions and consequently the symmetry of the assembled matrix.

Localized sum factorization

We have described sum factorization as a global assembling procedure that requires the pre-computation of the weights F on the whole domain. This is in contrast to the usual FEM assembling procedure which is performed element by element. The idea of this section is to consider a domain that is a union of disjoint boxes and to apply sum factorization on each box. Such a localization makes sense even if a global tensor-product discretization is chosen because the localization reduces the memory footprint required to store precomputed values of F and improves locality of the data which might lead to positive caching effects. The case of localizing a global tensor-product discretization is discussed in Section 4.1. It is worth noting that the possibilities opened by localized sum factorization go beyond that case and allow interesting applications, like multipatch domains and adaptive methods, see Section 4.2.

We start our discussion with a given collection P of pairwise disjoined, axis aligned boxes

D = Ś d δ=1 D δ ⊂ R d . For each box D, let Φ D = (φ D,n ) N D n=1
and Ψ D = (ψ D,m ) M D m=1 be generating systems of functions D → R and X D a quadrature on D fulfilling ( 5), ( 6), ( 7) and [START_REF] Giannelli | THB-splines: the truncated basis for hierarchical splines[END_REF].

Let Φ = (φ n ) N
n=1 and Ψ = (ψ m ) M m=1 be two generating systems of functions defined on the union of the boxes in P, which satisfy the following assumptions.

• Representability by box-local generators: For each box D ∈ P, we have

∀ n = 1, . . . , N, φ n | D ∈ Φ D ∪ {0}, ∀ m = 1, . . . , M, ψ m | D ∈ Ψ D ∪ {0}. (20)
• Bounded number of box-local generators: There is some R > 0 such that

D∈P N D ≤ RN where N D = # Φ D , (21)
D∈P #{n : supp φ n ∩ X D = ∅} ≤ RN. ( 22
)
Analogously to Lemma 1, we obtain the following statement.

Lemma 3. If Φ, Ψ ⊂ H r ([0, 1] d ) and Φ D , Ψ D ⊂ H r (D d ) satisfy any of the con- ditions (20), (21) or (22) then also ∂ η Φ, ∂ θ Ψ, ∂ η Φ D , ∂ θ Ψ D satisfy that condition for all θ, η ∈ Θ r . Remark 5. Note that (22) is equivalent to N n=1 #{D ∈ P : supp φ n ∩ X D = ∅} ≤ RN,
which follows from the stronger condition

∀ φ n ∈ Φ, #{D ∈ P : supp φ n ∩ D = ∅} ≤ R. ( 23
)
Remark 6. Provided that (20) holds and that for all D ∈ P we have

m = n ⇒ φ n | D = φ m | D ∨ φ n | D = φ m | D = 0,
then the condition (21) implies the condition (22) with the same value of R.

Let a : L 2 ([0, 1] d ) × L 2 ([0, 1] d )
→ R be a bilinear form of the form (1) and A be the matrix associated to the restriction of a to span Φ × span Ψ. Then, it follows from (20) that A can be decomposed as in

A = D∈P M Ψ,D A D N Φ,D . (24)
where the matrices N Φ,D and M Ψ,D are selection matrices whose columns contain at most a non zero coefficient. The cost decomposes in two parts:

C = C ass + C acc , (25)
where C ass is the cost of assembling the matrices A D for all D ∈ P and C acc is the cost of accumulating the A D into A.

The matrix products in (24) realize only a selection and relabeling of the non zero entries to sum to A, thus the cost C acc of accumulating the result in A is bounded by the number of selected non zero coefficients. As the number of non zero columns in

A D N Φ,D is ≤ #{n : supp φ n ∩ X D = ∅}
Concluding, C is dominated by C ass and thus we have the following analogues to (17), ( 18) and ( 19):

∀D, δ, #X D,δ ≤ N D,δ k δ p δ ⇒ C d 2 d R p d+2 N, ( 28
) ∀D, δ, #X D,δ ≤ N D,δ k δ ⇒ C d 2 d R p d+1 N, ( 29
) ∀D, δ, #X D,δ ≥ N D,δ k δ p δ ⇒ C d R p 2 #X. (30)

Localized sum factorization for tensor-product B-spline bases

In this subsection, we again consider the discretization from Example 3: Φ is a B-spline basis defined on [0, 1] d .

Again, the breakpoints of Φ define a partitioning of the domain [0, 1] d into elements on which the basis functions are polynomials. On each element, tensor product Gauss quadrature of order p δ is employed. Now, we localize sum-factorization to a partition P in which each box D is aligned with the element boundaries and we study the dependence of the assembling cost on the size of the boxes measured by the number of contained elements. The main advantages of partitioning the assembling process are reducing the memory requirement and allow for parallelization. By assembling a box at at time only F(X D ) needs to be pre-computed and stored instead of F(X). Assembling on different boxes can be done in parallel provided that data-races are avoided when the result is gathered to the final matrix.

First, we observe that the restriction to axes-aligned boxes does not destroy the tensor-product structure. Proof. We have Proof. The support of φ =δ n contains at most p δ elements in direction δ. The the number of boxes containing s δ elements intersecting supp φ =δ n is maximized if the end of the leftmost box coincide with the end of the leftmost element. In this case the total number of boxes intersecting supp φ =δ n is

φ n (x) = d δ=1 φ =δ π δ (n) (x δ ), and D = D 1 × • • • × D d . So, we obtain (φ n | D ) (x) = d δ=1 φ =δ | D δ (x δ ) and, therefore, φ n | D ∈ Φ 1 | D1 ⊗ • • • ⊗ Φ d | D d , and if φ n | D = 0 φ n | D ∈ (Φ 1 | D1 \{0}) ⊗ • • • ⊗ (Φ d | D d \{0}),
1 + p δ -1 s δ = s δ + p δ -1 s δ .
Using Remark 5, we obtain (22) and using Lemma 4 further (21).

Finally, we give an estimate for the number of quadrature points per box.

Lemma 6. If Φ is a tensor-product B-spline basis of order p 1 , . . . , p d , and we perform Gauss-quadrature with p δ quadrature points per element and direction, then we have for each box D of size

s 1 × . . . × s d ≥ 1 elements that #X D,δ ≤ s δ p δ s δ + p δ -1 ∈ [1, p δ [ N D,δ .
Proof. We have on the one hand for Gauss quadrature #X D,δ = s δ p δ and on the other hand s δ + p δ -1 ≤ N D,δ , where equality is obtained for B-splines of maximum smoothness.

Lemmas 5 and 6 show two competing effects on the cost of localized sum factorization. On one hand, the upper bound for R from Lemma 5 deteriorates for small boxes. On the other hand, Lemma 6 states that the ratio between quadrature points and basis function decreases for small boxes allowing one to apply the bound in (29) instead of that in (28). This observation allows to conclude as follows.

Example 5 (Per-element sum-factorization). Let Φ be a B-spline space of order p 1 , . . . , p d . In this example, P is the collection of the elements associated to Φ. The local spaces Φ D , Ψ D and the local quadrature X D are defined by Note that this example employs the sum factorization approach which was proposed in [START_REF] Antolin | Efficient matrix computation for tensor-product isogeometric analysis: The use of sum factorization[END_REF].

     Φ D,δ = {φ =δ n | D δ : φ =δ n ∈ Φ δ }\{0}, Ψ D,δ = {ψ =δ m | D δ : ψ =δ m ∈ Ψ δ }\{0}, X D is
Example 6 (Per macro-element sum-factorization). Let Φ be again a B-spline space of order p 1 , . . . , p d . In this example, we now consider a partition P of [0, 1] d such that all boxes have size of at least p 1 × • • • × p δ elements. Again, we assume (31). For this case, we obtain from Lemmas 5 and 6 and (28)

C d 4 d p d+2 N, (32)
i.e., asymptotically the same costs as for the global approach.

This description can be refined leading to a compromise that allows for more freedom in the choice of the partitioning P. The following Lemma is a refinement of Lemma 5 that allows for macro-elements that have size 1 in one direction, possibly different for each macro-element. It will be used in Example 7.

Lemma 7. For all 2 ≤ p 1 , . . . , p d ∈ N and covering P of

Ś d δ=1 [0, p δ [ such that for all D ∈ P D ∩ d ą δ=1 [0, p δ [ = ∅, (33)
D = [a 1 , b 1 [× • • • × [a d , b d [, a i , b i ∈ Z, ( 34
) b δ -a δ = p δ δ = j D 1 δ = j D ( 35
)
for some j D ∈ {1, . . . , d} that depends on D, the following sharp inequality holds

#P ≤ 2 d-1 d δ=1 (p δ -2) + 2 d .
Proof. We prove the bound by induction on the dimension d. For d = 1 any partition of [0, p] in segments of length 1 has exactly p elements. Let P be a partition realizing the maximum cardinality in dimension d and define P B to be the subset of P containing the boxes that intersect some box B. Then for

L = [0, 1[× Ś d δ=2 [0, p δ [, R = [p 1 -1, p 1 [× Ś d δ=2 [0, p δ [ we have #P = #P L + #P R + #(P \ P L∪R ).
A box D ∈ P \ P L∪R can not have length p 1 in the first direction. Thus, D has size

1×p 2 ו • •×p d and is contained in [i, i+1[×R d-1
for some i ∈ {1, . . . , p 1 -2}.

Applying Lemma 5 yields that for each i = 1, . . . , p 1 -2, there are at most 2 d-1 such boxes so that we have

#(P S \ P L∪R ) ≤ 2 d-1 (p 1 -2).
The projection (x 1 , . . . , x d ) → (x 2 , . . . , x d ) maps the boxes in

P L to disjoint boxes in R d-1 that cover Ś d δ=2 [0, p δ [. By induction hypothesis they are at most 2 d-2 d δ=2 (p δ -2) + 2 d-1
. The same reasoning applies to P R so that we obtain Narrow macro-elements are particularly of interest for the construction of hierarchical LR splines, [START_REF] Bressan | A hierarchical construction of LR meshes in 2D[END_REF], where local refinement is done in alternating directions. Note that additional costs are caused by the reordering of F and by the permutation of the coefficients in A D .

#P S ≤ 2[2 d-2 d δ=2 (p δ -2) + 2 d-1 ] + 2 d-1 (p 1 -2) = 2 d-1 d δ=1 (p δ -2) + 2 d . For d =
P L = {[1 -p 1 , 1[× D : D ∈ P} , P R = {[p 1 -1, 2p 1 -1[× D : D ∈ P} ,

Localized sum factorization for some non-tensor-product bases

Example 8 (Multipatch domains). Another example of localized sum factorization is Isogeometric Analysis on multipatch geometries. The domain Ω of the PDE is the union of the subdomains Ω , = 1, . . . , L each parameterized by a map G defined on Ω = [0, 1] d :

Ω = L =1 G ( Ω ), where Ω = [0, 1] d .
This means that the assembling is done on L copies of the parametric domain [0, 1] d . This can be expressed on the framework form the beginning of this section by writing

Ω := {1, . . . , L} × [0, 1] d .
The partition P contains the parametric patches

D := { } × [0, 1] d ,
where we have a tensor product spline space Φ D on each of those patches. The global space Φ is obtained by identifying matching functions that are non-zero on the glued boundaries of the patches.

The condition ( 21) and ( 22) are fulfilled according to Remark 5 with R that is the maximum number of patches on which a function φ ∈ Φ is active. In cases of practical interest the patches are big compared to the supports of the basis functions and their angles at the vertices (and edges in 3D) are close to π/2. Consequently the value of R is small ≈ 2 d . However, in contrived examples, it can be made arbitrarily big by constructing a multipatch in which many patches contain a common vertex.

Example 9 (HB-splines). The Hierarchical B-spline (HB) is a basis that breaks the global tensor product structure and allows for adaptive methods in which only a part of the domain is refined, see [START_REF] Buffa | Adaptive isogeometric methods with hierarchical splines: error estimator and convergence[END_REF]. The basis is obtained by selecting functions from different tensor product B-spline bases on different regions of the domain.

Let B 1 , . . . , B L be the tensor product B-spline bases of the same degree defined on the domain Ω that generate nested spaces

i < j ⇒ span B i ⊂ span B j . ( 36
) Let Ω ⊇ Ω 1 ⊇ • • • ⊇ Ω L = ∅ be corresponding closed domains.
For simplicity we further assume that the B i are a sequence of dyadically refined spaces and that each of the Ω i is a union of elements associated to B i . The hierarchical basis (HB-splines) is defined by Kraft's selection criteria [START_REF] Kraft | Adaptive and linearly independent multilevel B-splines, Surface Fitting and Multiresolution Methods[END_REF]:

H = L i=1 {γ ∈ B i : supp γ ⊆ Ω i and supp γ ∩ (Ω i \ Ω i+1 ) = ∅} . ( 37
)
A HB basis is said to be β-admissible if for all x ∈ Ω, max{i :

∃γ ∈ B i ∩ H : x ∈ supp γ} -min{i : ∃γ ∈ B i ∩ H : x ∈ supp γ} ≤ β + 1. ( 38
)
The localized sum factorization applies to β-admissible (usually 2-admissible) HB bases Φ and Ψ if the collection P is chosen such that each box D is contained in one of the rings

∆ i := Ω i \ Ω i+1 .
In this case, the restriction of the functions in H to any of the boxes D is contained in the union of β tensor product spaces. If each ring ∆ i allows a partition into macro-elements of size p δ elements per direction, then ( 21) and ( 22) hold with R = 3 d . For Gauss quadrature, we obtain using (28) the following bound for the costs C, bounded as follows:

C d6 d p d+2 N.
Note that narrow macro-elements are also possible as a consequence of Lemma 7.

Remark 7. Note that localized sum factorization cannot be applied directly to Truncated Hierarchical B-spline bases, cf. [START_REF] Giannelli | THB-splines: the truncated basis for hierarchical splines[END_REF], even if it spans the same space.

The problem is that truncation breaks the tensor structure, i.e. that the Φ D corresponding to the truncated basis do not fulfill ( 5). Nevertheless, it is possible to express the matrix A trunc corresponding to the truncated basis as

A trunc = M Ψ AN Φ
where A is the system matrix for the hierarchical basis and M Ψ , N Φ are the change of basis. Since the number of non zero entries in each column of M Ψ , N Φ is bounded by 2 βd+1 p d for β-admissible hierarchical splines we obtain that

C trunc = C + C post , and C post 2 βd+1 p 2d N Φ .
Interestingly it has the same asymptotic complexity as the cost for hierarchical splines for two dimensional domains (2d = d + 2). Remark 8. The extension to trimmed domains would be a challenge. So far, there is no clear solution for such discretizations. A straight-forward approach would be to compute the discretization ignoring the trimming first and recompute the matrix elements corresponding to the basis functions whose support intersects with the trimming curve. In reasonable cases, one could assume that the number of those matrix elements is small compared to the total number. Remark 9. The presented examples are not exhaustive. Surely, different generating systems fulfill the assumptions, for instance the hierarchical LR splines [START_REF] Bressan | A hierarchical construction of LR meshes in 2D[END_REF] do. Moreover, multipatch domains (Example 8) and adaptivity (Example 9) can be combined effortless.

Two algebraic descriptions

We called the described recursive algorithm sum factorization, but it differ from the original description in [START_REF] Antolin | Efficient matrix computation for tensor-product isogeometric analysis: The use of sum factorization[END_REF]. The unifying idea is that the map from the value of the coefficients on the quadrature points to the system matrix is linear and can be represented as a tensor. Moreover it has a Kronecker-like structure that is inherited from the tensor structure of the discrete spaces and of the quadrature.

Let X be lexicographically ordered and its points be x i1,...,i d with i δ = 1, . . . , #X δ and F be the d tensor defined as

F i1,...,i d = F(x i1,...,xi d ).
Then the 2d tensor T that corresponds to A through π and σ

T σ1(m),...,σ d (m),π1(n),...,π d (n) = a(φ n , ψ m ) = A m,n
can be written as

T m1,...,m d ,n1,...,n d = F i1,...,i d Q m1,...,m d ,n1,...,n d ,i1,...,i d
where, using Einstein notation, repeated indices denote sums and Q is the 3d tensor corresponding to the linear maps from F to T h , i.e., Q σ1(m),...,σ d (m),π1(n),...,π d (n),i1,...,i d = φ n (x i1,...,i d )ψ m (x i1,...,i d )ω(x i1,...,i d ).

The Cartesian structure of the quadrature points, together with the tensor structure of the generators allow for a factorisation of the above as

T m1,...,m d ,n1,...,n d = F i1,...,i d Q =1 m1,n1,i1 . . . Q =d m d ,n d ,i d .
The original description of sum factorisation corresponds to the partial evaluations T ≤δ m1,...,m δ ,n1,...,n δ ,i δ+1 ,...,

i d = F i1,...,i d Q =1 m1,n1,i1 . . . Q =δ m δ ,n δ ,i δ
where the temporaries for all i δ+1 , . . . , i d are computed at each step. This achieves a cost reduction because nnz(Q =δ ) ≈ k δ p δ N δ while nnz(Q) ≈ d δ=1 k δ p δ N δ and this impact the cost of computing the tensor contraction. The recursive description corresponds to the computation of a slice of T ≤δ m1,...,m δ ,n1,...,n δ ,i δ+1 ,...,i d corresponding to fixed indices i δ+1 , . . . , i d at a time.

They do exactly the same operation, although in a different order, but the recursive implementation requires less temporary memory. The complexity gain with respect to a standard per-element quadrature is given by the fact that in standard per-element assembling, the coefficients of the composition of the Q =δ are recomputed for each element.

The same factorisation can be explained, without recurring tensors, using a Kronecker construction. This second description correspond to the first by the flattening of the indices i 1 , . . . , i d . Let X be lexicographically ordered and define

Q Φ = [φ n (x)] x∈X n=1,...,N , Q Ψ = [ω(x)ψ m (x)] x∈X m=1,...,M
and F as the diagonal matrix containing F(x) for x ∈ X. Then

A = (Q Ψ ) F Q Φ . (39)
The factorization can be now expressed by decomposing Q Φ and Q Ψ as Kronecker products:

Q Φ = d K δ=1 Q =δ Φ , Q =δ Φ = [φ =δ n (x δ )] x δ ∈X δ n=1,...,N δ Q Ψ = d K δ=1 Q =δ Ψ , Q =δ Ψ = [ψ =δ m (x δ )ω δ (x δ )] x δ ∈X δ m=1,...,M δ .
Note that we have the following correspondence with the previous representation:

∀m δ , n δ , i δ Q =δ+1 m δ ,n δ ,i δ = Q =δ Ψ m δ ,i δ Q =δ Φ m δ ,i δ .
Since for all matrices C, D

C K D = (C K I)(I K D)
we have

A = (I K . . . K Q =d Ψ ) . . . (Q =1 Ψ K . . . K I) F (Q =1 Φ K . . . K I) . . . (I K . . . K Q =d Φ ). (40) Note that (Q =1 Ψ K . . . K I) F (Q =1 Φ K . . . K I) =     . . . A =1 x2,...,x d . . .     xi∈Xi
where the right hand side is the block diagonal matrix with blocks A =1 x2,...,x d for x i ∈ X i , i = 2, . . . , d. The blocks are lexicographically ordered along the diagonal. More generally after δ products we obtain

(I K . . . Q =δ Ψ . . . K I) . . . F . . . (I K . . . Q =δ Φ . . . K I) =     . . . A =δ x δ+1 ,...,x d . . .     xi∈Xi
The improved performance can then easily explained by the cost of matrixmatrix multiplication:

nnz(I K . . . Q =δ Φ . . . K I) ≤ #Xp δ while nnz(Q Φ ) ≤ #X d δ=1 p δ .
By computing the product step by step the cost δ p δ , while if the product is computed in one step the cost is δ p δ .

6 Pre-assembling costs

The overall costs presented in the last sections only covers the costs of assembling the matrix, assuming that the values of the functions in the generating sets Φ δ and Ψ δ and the values of the coefficient function F are already known or can be computed in constant time. This follows the presentation in the literature, however it is worth to have a closer look onto the corresponding computational costs. Our aim is to show algorithms for which the computation of Φ δ , Ψ δ and F can be performed at a cost that is inferior to the assembling cost. For simplicity, we focus on the global sum-factorization approach, but the results hold also for the localized version as they apply to each box separately.

Evaluation of the functions in the generating sets Φ and Ψ

The evaluation of a B-spline function φ n or ψ m cannot be done in constant time.

As we have a tensor-product structure, an efficient approach is to pre-compute the function values for the corresponding univariate functions φ =δ n or ψ =δ m . For a standard B-spline basis, this yields costs C, which are bounded as follows

C d δ=1 p 2 δ #X δ ≤ 2p 3 d δ=1 N δ , (41)
where the second estimate holds for #X δ ≤ k δ p δ N δ with k δ = 2 as in ( 17) and (18). The resulting costs are always much smaller than p d+2 N , the costs of assembling.

Evaluation of the coefficient function F

As mentioned in Section 2, we have to evaluate the function F for every quadrature point. The computational costs related to this task are often not addressed in standard papers on assembling in Isogeometric Analysis. If the computational costs for evaluating F at one quadrature point, grows faster than p 2 , then a direct approach for evaluating F would lead to computational costs which exceed the costs of the assembling strategy.

In classical Isogeometric Analysis, we assume that the computational domain Ω is parameterized by a diffeomorphism

G : Ω := [0, 1] d → Ω = G( Ω) ⊂ R s ,
which is an element of span Φ and has the form

G = N n=1 c n φ n .
In this case, also the evaluation of F involves the evaluation of G and/or its derivatives, cf. [START_REF] Hughes | Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement, Comput[END_REF]. For ease of the presentation, in the following we only present the evaluation of G, the evaluation of its derivatives is completely analogous.

For the evaluation of G at the quadrature points X, we can again exploit the tensor product structure of Φ. This yields the following algorithm:

Algorithm 2 Recursive function computation

1: procedure Eval([X δ ] d δ=1 , [Φ δ ] d δ=1 , [c n ] N n=1 ) 2:
if d=0 then return G ≤d 13: end procedure Again, we derive the number of floating point operations, assuming that the functions in the generating set Φ δ have already been evaluated. By counting the number of invocations of lines 7 and 9 and using [START_REF] Calabrò | Fast formation of isogeometric Galerkin matrices by weighted quadrature[END_REF], we obtain that the costs C ≤d satisfy

G ≤d-1 n d ←Eval([X δ ] d-1 δ=1 , [Φ δ ] d-1 δ=1 , [c n ] n:π d (n)=n d ) 8: for all x d ∈ X d : φ =d n d (x d ) = 0 do 9: G ≤d (x, x d ) ← G ≤d (x, x d ) + φ =d n d (x d )G ≤d-1 n d (x)
C ≤d N d C ≤d-1 +sp d #X ≤d ,
where s comes from the fact that the c i are vectors in R s . Recursively plugging in yields

C ≤d s d δ=1 d i=δ+1 N i p δ #X ≤δ .
Depending on the ratio between #X δ and N δ and assuming k δ = 1, as with Gauss quadrature, we obtain the following costs

∀δ, #X δ ≤ k δ p δ N δ ⇒ C ≤d sd p d+1 N, (42)
∀δ, #X δ ≤ k δ N δ ⇒ C ≤d sd p N, (43)
∀δ, #X δ ≥ k δ p δ N δ ⇒ C ≤d sd p #X. (44)
Here, d is the dimension of parametric domain and s is the dimension of the target space. So, for a surface embedded in R 3 , we have d = 2 and s = 3.

The computation of the map is not sufficient to IGA applications: as shown in Example 1, it is often the case that the Jacobian of G, its pseudoinverse and the determinant of the Jacobian are required. The i-th column of J G is computed using the presented algorithm by replacing Φ i with ∂Φ i . Remembering that J G has d columns, we have the following costs:

∀δ, #X δ ≤ p δ N δ ⇒ C jac sd 2 p d+1 N, (45) ∀δ, #X δ ≤ N δ ⇒ C jac sd 2 p N, (46) ∀δ
, #X δ ≥ p δ N δ ⇒ C jac sd 2 p #X. (47)
For second order PDEs, the coefficients F are computed as in [START_REF] Hughes | Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement, Comput[END_REF]. If s = d, i.e. G : R d → R d then J -1 G can be computed in ≈ d 3 operations per quadrature point using LU factorization. If G : R d → R s , s > d then the formula in [START_REF] Hughes | Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement, Comput[END_REF] is still applicable by denoting with J -1 G the left pseudoinverse of J G and with |J G | the area density, i.e.

J -1 G := (J G J G ) -1 J G , (48)
|J G | := det(J G J G ) 1 2 . ( 49
)
Both can be computed with ≈ sd 2 operations per quadrature point using LU factorization to invert J G J G giving

C inv sd 2 #X.
The costs C mult of evaluating the product in ( 10) is of the same order as the costs for evaluating the block J -T G AJ -1 G . This has costs of order ≈ s 2 d per quadrature point. Thus the costs of computing F according to [START_REF] Hughes | Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement, Comput[END_REF] decomposes as

C F = C jac + C inv + C mult C jac + s 2 d #X and, since s ≥ d, we conclude ∀δ, #X δ ≤ p δ N δ ⇒ C F s 2 d p d+1 N, (50) ∀δ, #X δ ≤ N δ ⇒ C F s 2 d p N, (51)
∀δ, #X δ ≥ p δ N δ ⇒ C F s 2 d p #X. ( 52
)
7 A note about matrix-free implementations

In this section, we remark on matrix-free implementation. Matrix-free methods are often combined with iterative solvers. Iterative solvers typically require only the ability to compute matrix-vector products v = Au, while a direct access to the entries of A is not required. By definition

Au = [0,1] d ψ m (x)u(x) M m=0
, where

u(x) = N n=1 u n φ n (x).
For this case, unfortunately, the cost of matrix-free setting must be counted based on M , i.e., on the dimension of the test space. So, we assume N = M for simplicity.

Algorithm 1 can be changed to a matrix-free setting by removing the loop at line 8 and pre-multiplying F(X) by u(X). This gives the following algorithm, where the blocks w ≤d m d are the blocks of

M ≤d-1 coefficients in v = v ≤d : v = (v 1 , . . . , v M ≤d-1 w1 , . . . , v M ≤d -M ≤d-1 +1 , . . . , v M ≤d w M d
).

Algorithm 3 Matrix-free sum-factorization

1: procedure Apply([X δ ] d δ=1 , [ω δ ] d δ=1 , [Ψ δ ] d δ=1 , [F(x)u(x)] x∈X ,) 2:
if d=0 then v ≤d ← 0 start with null vector 6:

for all x d ∈ X d do sum all v ≤d-1

x d 7: v ≤d-1 x d ←Apply([X δ ] d-1 δ=1 , [ω δ ] d-1 δ=1 , [Ψ δ ] d-1 δ=1 , [F(x, x d )u(x, x d )] x∈X ≤d-1 ) 8:
for all m d with ψ =d m d (x d ) = 0 do 9:

w ≤d m d ← w ≤d m d + ω d (x d )ψ =d m d (x d )v ≤d-1 x d 10:
end for 11:

end for 12:

return v = v ≤d 13: end procedure

The above yields a cost that, excluding the computation of u(X), is

C ≤d app #X d C ≤d-1 app + p d M ≤d-1 .
This, remembering the simplification M δ = N δ , leads analogously to (15) to the following bound for the costs:

C ≤d app d i=1 d δ=i+1 #X δ (p i #X i ) N ≤i-1 . (53)
Depending on the relative ratio between N δ and #X δ , we obtain

∀δ, #X δ ≤ p δ N δ ⇒ C ≤d app d p d+1 N, (54) ∀δ
, #X δ ≤ N δ ⇒ C ≤d app d p N, (55) ∀δ
, #X δ ≥ p δ N δ ⇒ C ≤d app d p #X. (56)
By comparing the estimate (55) with ( 18) and (51), we observe that a matrix free approach can achieve significant saving if weighted quadrature is employed as proposed in [START_REF] Sangalli | Matrix-free weighted quadrature for a computationally efficient isogeometric k-method[END_REF]. Indeed in this case, both the computation of the coefficients and the application of the matrix are computed with a cost ≈ pN . Note that this is a factor of p d less than the assembling of the matrix and the same factor less than the cost of the multiplication of the assembled matrix with some given vector.

However, by comparing the estimates estimate ( 54) and ( 56) with ( 17), ( 19), ( 50) and (52), we observe that for standard Gauss quadrature, the potential saving is "only" for a factor of p.

Numerical experiments and conclusions

We have implemented the global strategy and the per-(macro)-element strategies and have tested their behavior for a few sample problems. Our implementation is a C++ code which has been carefully optimized. It is available online1 as an stand-alone assembling library.

In the numerical experiments, we assemble a standard stiffness matrix for the 2D and 3D domains depicted in Fig. 1. The 2D domain is decomposed into 200 × 200 elements, the 3D domain into 22 × 22 × 22 elements. On each domain, we assemble the stiffness matrix for splines of several orders p using Gauss quadrature of the same order and we compare the time used by different algorithms. This was done on a single socket machine with an Intel(R) Core(TM) i3-8100 CPU running at 3.60GHz.

In Fig. 2 assembling uses the same Gauss quadrature and was performed using the freely available IGA library G+Smo [START_REF] Bressan | G+Smo[END_REF].

In Fig. 4 and 5, we focus on the effect of the macro-element size. The considered approaches are: standard macro-elements (macroS) corresponding to a size of p × p in 2D and p × p × p in 3D, narrow macro-elements (macroN) having size p × 1 in 2D and p × p × 1 in 3D and rotated narrow macro-elements (macroR) having size 1 × p in 2D and 1 × p × p in 3D. According to Example 7, the narrow dimension should come last in sum factorization. As our code does not reorder the dimensions we expect the same behavior for macroS and macroN, but a cost higher by a factor p for macroR. This means that macroR behaves as element in 2D and as p 6 N in 3D. We observe that in any case, the sum factorization approaches are faster than the standard approach. Conforming with the theory, we see that global is significantly faster than element. Moreover, we obtain that the macro-element approaches macro and macroN are indeed almost as fast as the global sumfactorization approach global.

The results are used for fitting the parameters c and e in the formula t = cp e N (p), where t is the measured time, N (p) = d δ=1 (p+K δ -1) is the number of degrees of freedom and K δ is the number of elements in the corresponding direction, i.e., K 1 = K 2 = 200 for d = 2 and K 1 = K 2 = K 3 = 22 for d = 3. The fitted curves are dashed in Fig. 2 Table 1: Speed-up factors for assembling the stiffness matrix using a parallel macro-element implementation.

We observe that the exponents e obtained in our experiments are significantly lower than those predicted by the theory. We believe that the huge difference between the speed in performing computation and the speed in accessing memory in modern processors, is masking one order in p. For the 3D example, there is also another reason: the fitting is distorted by the small number of elements, i.e., the approximation p d N ≈ #X does not apply: L δ = 22 so that pN δ = 22p + p 2 -p is twice bigger than #X δ = 22p for p = 20.

In Fig. 6 and 7, we report the assembling times using a proof-of-concept multithreaded approach based on macro-elements of size p × p in 2D and p × p × p in 3D (macroS). To have a sufficient number of macro-elements the 2D domain is split in 2000×2000 elements and the 3D domain is split in 64×64×64 elements.

The lines correspond to the polynomial orders p = 4, 6 and 8, and the abscissa is the number of parallel threads. The tests have been executed on a dual socket machine with Intel(R) Xeon(R) CPU E5-2695 v4 running at 2.10GHz. The figures show that the macro-element approach is viable to parallelization. Table 1 shows the corresponding speed-up factors, i.e. the ratio between the time for the single-threaded execution and the time for the multi-threaded execution.

In this paper, we have developed a unified complexity analysis for sum factorization approaches. The theory shows that in Isogeometric Analysis, the computational costs can be reduced significantly by using the technique of sum factorization. This might be of interest because sum factorization can use the standard Gauss quadrature yielding exactly the same matrix as the standard approach.

We show that sum factorization does not yield its optimal complexity if it is applied on each element separately. However, we do not need to consider the global approach to obtain its optimal complexity: it is sufficient to apply it to blocks of at least p elements in each direction.

Moreover, we have shown that parallel implementations of localized sum-factorization are a viable strategy for fast assembling of the system matrix in IgA applications.

  and each column contains at most d δ=1 k δ p δ non zeros because of (8), we obtain C acc D∈P d δ=1 k δ p δ #{n : supp φ n ∩ X D = ∅}. Using (22) and assuming k δ = 2, we further obtain C acc 2 d p d RN. If #X D,δ ≤ k δ p δ N D,δ for all D ∈ P and δ = 1, . . . , d, then (17) and (21) yield C ass D∈P d2 d p d+2 N D ≤ d2 d Rp d+2 N. (26) Similarly, under the stronger hypothesis #X D,δ ≤ k δ N δ , (18) and (21), we obtain C ass D∈P ≤ d2 d Rp d+1 N.

Lemma 4 .

 4 Assume that Φ and Ψ satisfy[START_REF] Bressan | G+Smo[END_REF] andP is a partition of [0, 1[ d into boxes of the type [a 1 , b 1 [× • • • × [a d , b d [, thenΦ D := {φ n | D : n = 1, . . . , N }\{0}, Ψ D := {ψ m | D : m = 1, . . . , M }\{0} satisfies (5) for all D. For this choice, (22) implies (21).

  which shows (20). The same is possible for Ψ.If supp φ n ∩ X D = ∅ then 0 = φ n | D and φ n | D ∈ Φ D . Then, we have Φ D = {φ n | D : φ n | D = 0}, so (22) implies (21).As next step, we show (22) for B-spline bases. Lemma 5. If Φ is a tensor-product B-spline basis of order p 1 , . . . , p d , then taking P such that each D ∈ P contains at least s 1 × . . . × s d ≥ 1 elements, yields (21) and (22) with R = d δ=1 s δ + p δ -1 s δ .

  the standard Gauss quadrature on the box D.

  (22) with R = p d and, by Lemma 4, we have (21). In this case #X D,δ = N D,δ = p δ so that (29) holds and, as proved in[START_REF] Antolin | Efficient matrix computation for tensor-product isogeometric analysis: The use of sum factorization[END_REF],C d 2 d p 2d+1 N.

  and completing P L ∪ P R to a covering P of Ś d δ=1 [0, p δ [ using the boxesD = [i, i + 1[× d ą δ=2 [t i p δ + 1, (t i + 1)p δ + 1[ for i = 1, . . . , p 1 -2 and t i ∈ {-1, 0}.Example 7 (Narrow macro-elements). Using the last lemma, we obtain that macro-elements D ∈ P of size s 1 × • • • × s d with s δ ≥ p δ for δ = j D and s j D ≥ 1 achieve the same asymptotic costs in p as we have obtained in Example 6. From Lemma 7 and Remark 5, we have (21) and (22) with R = d2 d-1 p. On each macro-element D the cost C D d2 d p d+1 N D can be achieved by permuting the coordinates so that j D , the direction of size one element, comes last. Indeed, if j D = d, we have #X D,d = #N D,d = p δ and the cost in (18) is achieved. Since the estimate holds for all boxes, we have (29) and C d 2 4 d p d+2 N.

6 :

 6 for all n d ∈ {1, . . . , N d } do 7:

and 3 ,

 3 we report the assembling times for standard per-element assembling (standard), per-element sum factorization (element), per-macro-element sum factorization (macroS) and global sum-factorization (global). The macroelements have size p × p in 2D and p × p × p in 3D. Standard per-element 3D domain: bended and twisted box

Figure 1 :

 1 Figure 1: Computational domains

Figure 2 :Figure 3 :

 23 Figure 2: Time in seconds for assembling the stiffness matrix on the 2D domain depending on the polynomial order p.

Figure 4 :

 4 Figure 4: Time in seconds for assembling the stiffness matrix on the 2D domain with macro-elements depending on the polynomial order p.

Figure 5 :

 5 Figure5: Time in seconds for assembling the stiffness matrix on the 3D domain using macro-elements depending on the polynomial order p.

Figure 6 :

 6 Figure 6: Times in seconds for assembling the stiffness matrix on the 2D domain using a parallel macro-element implementation.

Figure 7 :

 7 Figure 7: Times in seconds for assembling the stiffness matrix on the 3D domain using a parallel macro-element implementation.

  1 the result is sharp as we can cover [0, p[ with p segments of length 1.

	Given a covering P of we construct a maximal cardinality covering of Ś d δ=2 [0, p δ [ with cardinality 2 d-2 d δ=2 (p δ -2) + 2 d-1 Ś d δ=1 [0, p

δ [ by setting

  -5. The fitting procedure yields the following values for the exponent e: .0 3.7 5.2 6.6 9.0 10.9 12.[START_REF] Giannelli | THB-splines: the truncated basis for hierarchical splines[END_REF] 14.3 15.2 15.6 17.6 3 7 1 2.0 3.7 5.0 6.6 8.3 10.9 11.4 13.6 13.6 17.4 17.5

							threads
	dim p 1	2	4	6	8	12	16	20	24	28	32	36
	2	3 1 1.9 3.5 4.7 5.9 7.6	9.3 10.6 11.3 12.4 13.5 12.8
	2	5 1 2.0 3.6 4.9 6.2 8.5 10.3 11.9 14.0 15.2 15.8 16.1
	2	7 1 2.0 3.8 5.1 6.6 8.9 11.0 13.2 15.2 16.9 18.4 19.8
	3	3 1 1.9 3.7 5.2 6.6 8.9 11.0 12.9 14.2 16.7 15.0 17.4
	3	5 1 2
				2D domain			3D domain
				theory experiments theory experiments
		standard		6		4.98	9		8.05
		element		5		3.92	7		5.18
		macroN		4		3.26	5		3.75
		macroR		5		3.99	6		4.51
		macroS		4		3.30	5		3.76
		global		4		3.47	5		3.70
							29

https://github.com/IgASF/IgASF

Acknowledgments

The first author has received funding from the European Research Council under the European Union's Seventh Framework Programme (FP7/2007-2013) / ERC grant agreement 339643. The second author was partially funded by the Austrian Science Fund (FWF) under grant NFN S117-03


Zulma Sandoval

Miguel Castro

Javad Alirezaie

Francis Bessière

Cyril Lafon

Jean-Louis Dillenseger
email: jean-louis.dillenseger@univ-rennes1.fr

Transesophageal 2D ultrasound to 3D computed tomography registration for the guidance of a cardiac arrhythmia therapy

Introduction

Catheter ablation is regularly used to treat cardiac arrhythmia. This intervention is required when antiarrhythmic medications are ineffective. Depending on the arrhythmia, the success rate varies from 30% to more than 95% [START_REF] Huang | Catheter Ablation of Cardiac Arrhythmias[END_REF]. One of the critical points is the ability to diffuse via the catheter enough energy at the arrhythmogenic area that is sometimes deeply located into the muscle. The more challenging the arrhythmia is to ablate, the more complications can occur. The complication rate is about 1% to 3% [START_REF] Huang | Catheter Ablation of Cardiac Arrhythmias[END_REF]. Strokes represent the most serious one, and vascular injuries are the most frequent. The more challenging the arrhythmia is to ablate As high intensity focused ultrasound (HIFU) offers the possibility to create a lesion at a distance from a transducer using dynamic focusing, methods are developed to treat cardiac arrhythmia non-invasively and more efficiently by introducing a transducer into the esophagus [START_REF] Yin | Noninvasive transesophageal cardiac thermal ablation using a 2-D focused, ultrasound phased array: a simulation study[END_REF] [START_REF] Constanciel | Design and evaluation of a transesophageal HIFU probe for ultrasound-guided cardiac ablation: simulation of a hifu mini-maze procedure and preliminary ex vivo trials[END_REF]. Since the esophagus is located just behind the heart, it offers a perfect acoustic window so that HIFU can be directed toward the heart without many acoustic fences [START_REF] Ho | Anatomy for Cardiac Electrophysiologists: A Practical Handbook[END_REF]) so this technology could allow treating foci without intravascular catheter. One of these approach is actually under feasibility study on animal testing [START_REF] Bessiere | Ultrasound-guided transesophageal high-intensity focused ultrasound cardiac ablation in a beating heart: A pilot feasibility study in pigs[END_REF]. The probe is combining both a multielement HIFU and ultrasound (US) transducers so that HIFU therapy can be delivered in real time under US 2D anatomy imaging. During the therapy, the probe would be introduced into the esophagus and placed behind the heart. From this location, HIFU waves delivered from the esophagus ablate the inner myocardial tissue without damaging intervening and surrounding tissues. The US imaging transducer aims to enable optimal positioning of the transducer as well as to verify the absence of obstacles in the acoustic beam path of the HIFU transducer. However, the limited field of view of the 2D US image, the image quality and both cardiac and respiratory motions make the guidance a very complex task. In order to enhance the guidance of the transesophageal HIFU cardiac ablation, we propose an integration of the high-resolution CT anatomical preoperative 3D imaging in which the ablation path has been defined with the intraoperative 2D US images. This would help the therapist to adjust the HIFU focal point to the the planed ablation path.

Registration of 2D US to 3D CT, using electromagnetic tracking system is reported in the literature (Linte, Moore, Wedlake, Bainbridge, Guiraudon, Jones & Peters 2009) [START_REF] Huang | Dynamic 2D ultrasound and 3D CT image registration of the beating beart[END_REF] [START_REF] Lang | Surface-based CT-TEE registration of the aortic root[END_REF]. In these approaches, preliminary registration is obtained using information of the tracking system (i.e. position and orientation of the probe tip in an external reference coordinate system). Then, final registration is performed between the CT volume and US volume reconstructed using the tracking information of the acquired 2D US images. The slice to volume registration problem is converted into a volume to volume registration. However, disadvantages such as accuracy, sensitivity to metallic interferences, requirement of additional hardware on the transesophageal probe, sensor calibration before intervention and accuracy losses in reused sensors make the integration of electromagnetic tracking during the clinical routine workflow very complex [START_REF] Franz | Electromagnetic tracking in medicine: a review of technology, validation, and applications[END_REF] [START_REF] Gao | Registration of 3D trans-esophageal echocardiography to X-ray fluoroscopy using image-based probe tracking[END_REF]. In addition, no accuracy improvement is reported by Lang et. al. when comparing registration methods in phantom studies using electromagnetic tracking system with clinical studies without tracking information [START_REF] Lang | Surface-based CT-TEE registration of the aortic root[END_REF].

In this paper, we propose to perform a 2D-3D (slice-volume) registration of the intraoperative 2D US and the preoperative CT without any external tracking system. More precisely the 2D-3D registration consists of finding the 3D pose (location and orientation) of the US image slice inside the preoperative 3D volume using only imagebased information.

Methods

To estimate the pose of the 2D US image slice inside the preoperative 3D volume, we have to find a way to reduce the amount of degrees of freedom. Our proposal is to find the pose of the 2D US image slice inside the preoperative 3D volume, exploiting the specific geometrical restrictions involved in this HIFU therapy with a transesophageal approach.

The first consideration is that the HIFU probe's trajectory is constrained by the esophagus. In fact, the position of the esophagus in the human body is highly constrained by its surrounding organs and tissues such as the vertebral column, trachea, carotid-jugular vascularization, aortic arch, right pulmonary artery, left main bronchus, left atrium and diaphragm [START_REF] Ho | Anatomy for Cardiac Electrophysiologists: A Practical Handbook[END_REF]. Thus, we made the assumption that the probe insertion does not modify the esophageal position in the body. This constraint has a direct impact on the potential geometrical position of ultrasound probe in the human body. Consequently, for a specific US image, its origin (the apex of the image sector) is given by the position of the probe along the esophagus. A second geometrical constraint is imposed by our specific transesophageal HIFU probe. Counter to classical transesophageal echocardiography (TEE) imaging probes in which the US transducer that can rotate from 0 degree to 180 degrees to obtain different viewing planes, our prototype of the HIFU therapy probe used a US transducer which is not able to rotate. Thus, the plane acquired with this prototype corresponds with the 0 degrees view. Because the tip of the HIFU therapy probe is inside of a water-filled balloon that is used as acoustic coupling medium, we considered that the US image is perpendicular to the esophagus central axis, matching the 0 degrees orientation of TEE.

Assuming that, we were able to reduce the number of parameters to i) the depth d of the US probe along the esophagus and ii) its orientation θ around the centerline of the esophagus. In order to compensate for deviation of the esophagus centerline estimation, we also allowed a light shift (a translation) of the transesophageal probe tip position from the centerline. This choice will drive the global framework into two stages (Fig. 1): 1) the reformatting of the preoperative CT dataset according to the esophagus topology. The main idea is to provide CT 2D slices which potentially have the same spatial location (and so the same information) as the future TEE US 2D images; 2) an intraoperative intensity-based registration between the US image and the reformatted CT 2D slices obtained in the previous stage.

Esophagus-based volume reformatting

The goal of this resampling or reconstruction is to observe the anatomical structures in CT in the same order and orientation as the future 2D US TEE images that would be acquired during the procedure.

The preoperative volume is denoted V pre with an associated world coordinate system (O pre ; x pre , y pre , z pre ) (Fig. 2-left). The TEE probe is inserted in the esophagus at a depth d from a specific insertion point O end . We assume, for tracking purpose, that 2D US TEE images are acquired with an orientation of zero degrees, therefore the images are located on an axial plane perpendicular to the probe axis. Based on our assumption that the probe's trajectory is constrained by the esophagus, during the insertion the US images will be a set of planes perpendicular to the esophagus centerline. From this ascertainment we proposed to create a therapy-adapted preoperative volume V ther which is made by all the sections perpendicular to the axis of the esophagus. The reference of this new volume is set as (O end ; x ther , y ther , d) with O end the transesophageal probe insertion point; axis d corresponds to the esophagus' centerline; and x ther y ther define the image plane (Fig. 2-right).

The procedure to reconstruct the therapy-adapted preoperative CT image was divided into three steps: 1) the segmentation of the esophagus from the preoperative volume; 2) the definition of the centerline as a smooth spline function; and 3) the resampling of the CT using the centerline information. 2.1.1. Segmentation of the esophagus An automatic segmentation of the esophagus was out of our scope. So in a first step we chose to extract the inner surface of the esophagus from the preoperative volume using the software ITK-SNAP [START_REF] Yushkevich | User-guided 3D active contour segmentation of anatomical structures: significantly improved efficiency and reliability[END_REF] which proposes an implementation of the active contours segmentation in a semi-automatic mode. The user defines first a region of interest and intensity thresholds to extract the potential esophagus voxels. The active contour segmentation is initialized by spherical surfaces covering the esophagus along its major axis. The snake evolution can be manually supervised in order to avoid leakages toward other structures. A closing operator is finally used to fill the small holes.

Esophagus centerline function

On slice i of V pre , let E i , the barycenter of the segmented esophagus, be a point of the centerline. This gave us a set E = {E 1 , E 2 , . . . , E N } of N centerline points. In order to obtain a smooth centerline, we fit a cubic smoothing spline (De Boor 1978) parametric function f (t) to E. On this spline, a smoothing parameter controls the trade-off between fidelity to the data and smoothness of the function estimate. In our case, we suggest to tune this parameter to get a very smooth curve to prevent the impact of the segmentation inaccuracy. The spline function was then regularly sampled using the natural parameter d(t) = t t 0 f (u) du. This defined the d axis of the therapy-adapted preoperative volume (Fig. 2-right). We had now a set of M regularly sampled coordinates {d 1 , d 2 , . . . , d M } along the d axis in V ther for which we have a set of M corresponding 3D points {E(d 1 ), E(d 2 ), . . . , E(d M )} and tangents to f at these points {f (d 1 ), f (d 2 ), . . . , f (d M )} along the esophagus curve in V pre .

2.1.3. Resampling of the CT For a specific coordinate d i , the plane (x ther , y ther , d i ) in V ther corresponds to the 3D plane P (d i ) in V pre which includes the 3D points E(d i ) and is perpendicular to the esophagus curve at E(d i ). More precisely, on this plane, we were able to define a coordinate system R i = (E(d i ); u i , v i , w i ) where w i is the tangent f of the curve f at P (d i ) ; u i is arbitrary set in a plane perpendicular to y pre : u i = w i × y pre ; and v i is perpendicular simultaneously to u i and w i : v i = u i × w i . We can now define a rigid transformation matrix T(d i ) between the reference R i and the reference (O pre ; x pre , y pre , z pre ) of V pre .

The value at any point (x t , y t , d i ) of V ther can now be given by the corresponding value at point (x, y, z) in V pre with:

[x, y, z, 1] T = T(d i )[x t , y t , 0, 1] T (1)
It should be noted that we have a specific transformation matrix T(d i ) for each slice d = d i of the therapy-adapted preoperative volume.

Intensity-based image registration

In this section, we will perform the registration between a specific US image (denoted I U S ) and V ther . Because the TEE is ECG-gated, we were able to select from the US sequence the frames corresponding to the same cardiac phase as that of the CT volume (in our specific case 75% of the cardiac cycle (R-R interval)). The US image I U S to register was concurrently manually selected by an expert from the acquired sequences. The automatic selection was out of scope of this work, but we expect that such a selection becomes available in the future.

Because of the specificity of V ther , the registration consists of performing a 2D-2D registration between I U S and a (x t , y t ) slice of V ther at a specific d i . We separate the overall process into three steps: 1) finding a candidate zone of slices in V ther ; 2) performing the US 2D-CT 2D registration for all the slices of this zone and; 3) estimating the most probable depth position of the probe in the zone.

Initialization and candidate zone definition

We used various information about the transesophageal probe (inserted length, visual analysis of the image sequence during navigation, etc.) to define a candidate zone along the esophagus centerline where the image transducer center can potentially be. The size of this zone corresponded broadly to 10 mm along the esophagus. According to the depth sampling step along the esophagus, this zone covered a set of 30 to 80 adjacent slices of V ther .

2D-2D registration selection

In this section we will consider the 2D-2D registration between I U S and a specific slice of coordinate d i in the candidate zone of V ther .

Similar to the initialization step, we first aligned the US image transducer center (the vertex of the field of view in I U S ) and the esophagus center point (0, 0, d i ) in the slice d i of V ther .

The aim of 2D-2D registration is to estimate the transformation matrix which best aligns the information in I U S to the information in the slice d i of V ther .

To perform this alignment it was necessary to define three important issues: i) the kind of alignment (rigid or elastic transform), ii) the constraints of the alignment (the parameters of the transformation) and iii) the measure of the alignment (similarity metric).

The choice for the "kind of alignment" was driven by the study of the movements and deformations introduced by the cardiac and respiration cycles [START_REF] Sra | Cardiac image registration of the left atrium and pulmonary veins[END_REF]. In our specific case, the registered images correspond to the same cardiac phase, the deformations caused by respiration motion are small in the region of interest for therapy, and the movement of the heart during respiration is in the same direction of the probe insertion. For these reasons we chose to perform a rigid registration because the deformations caused by the cardiac cycle are mostly avoided in our application.

To define the "constraints of the alignment" in this 2D registration, three parameters have to be estimated including the rotation angle and the two 2D coordinate positions of the rotation center.

The selection of the similarity measure depends on the nature of the image information of the involved modalities. For CT the value of the voxel intensity is related to the Hounsfield scale which gives a relatively global and homogeneous information about each tissue. In contrast, beside some specular reflections on the boundaries, in US imaging the information is mainly presented as speckle that shows some boundaries between tissues but also produces an inhomogeneous value distribution inside each tissue. There are two ways to compare image information in the process of alignment: feature-based or intensity-based. Feature-based comparison needs an additional segmentation step of the objects of interest in both image modalities. This process is usually time or computational expensive and the segmentation errors could also propagate to the posterior registration and which could have a direct impact on the final accuracy. For these reasons we chose rather to use an intensity-based approach which extracts information from gray levels without any additional process. In cardiac imaging such as US-CT, intensity-based multimodal registration has been reported by [START_REF] Huang | Dynamic 2D ultrasound and 3D CT image registration of the beating beart[END_REF]). The authors used Mutual Information (MI) to directly (or after a simple thresholding) compare CT and US cardiac images. In a previous work, we evaluated the use of different intensity-based measures reported in the literature in our application [START_REF] Sandoval | Intensity-based similarity measures evaluation for CT to ultrasound 2D registration[END_REF]. In this study, the best similarity measures were the Woods Criterion (WC) and the Mutual Information (MI). Between these 2 methods, MI was selected because it is implemented in most of the medical image libraries and packages.

Therefore, for a slice d i we performed a 2D-2D rigid intensity-based registration with MI as similarity measure. The resulting transformation matrix is denominated

T 2DU S-2DCT (d i ).
2.2.3. Best registration result After the 2D-2D registration process on all the slices of the candidate zone, the slice of depth d which presents the highest similarity measure value after registration was considered as the final result of the registration scheme. The estimated 2D transform T 2DU S2DCT (d) associated with the geometrical transform T(d) (eq. 1) gave the final 2D-3D transformation matrix :

T 2DU S-3DCT = T(d) • T 2DU S-2DCT (d) (2)
Any 2D point (x U S , y U S ) of I U S can be associated to a a 3D point (x pre , y pre , z pre ) of V pre by :

[x pre , y pre , z pre , 1] T = T 2DU S-3DCT [x U S , y U S , 1] T (3)
The algorithm was implemented in C++ using the Insight Segmentation and Registration Toolkit (ITK) library [START_REF] Ibanez | The ITK Software Guide, second edn[END_REF].

Results

Experimental protocol

Two experiments were conducted to assess the accuracy of the proposed registration method, one on a physical phantom data, the second on real clinical data.

On either case, to perform validation for a specific US image I U S we compared the following:

• the 3D pose of I U S estimated using our proposed method (section 2). The result of this estimation is a 2D-3D transformation matrix (eq. 2) we called T RecEval .

• a ground truth 3D pose of I U S . In both experiments, because no localization device could be placed on the intracavitary probe, we had to estimate directly the pose in the original CT volume V pre by a human expert using visual referential points. The result of this estimation is a 2D-3D transformation matrix we called T Ref .

The validation can be done by estimating the errors between these two 3D poses of I U S . To quantify the error, we defined specific landmarks in I U S and used the two transformations matrices T RecEval and T Ref to project these points in V pre . The estimated error, called the Target Registration Error (T RE), will be the distance between these two projected points [START_REF] Fitzpatrick | The distribution of target registration error in rigid-body point-based registration[END_REF]:

T RE(P j ) = T Ref P j -T RecEval P j , (4)
where P j denotes a landmark on the US image.

The global accuracy of the CT reconstruction can then be evaluated using the mean Target Registration Error (mT RE) defined on L landmarks in I U S :

mT RE = 1 L L j=1
T RE(P j ) (5)

Phantom study -Validation of the esophagus-based volume reformatting

This study aimed to prove that it is possible to reconstruct CT volumes based on the esophagus position with slices which correspond to TEE images (section 2.1). Hence, we compared the reconstructed CT planes and TEE images of a phantom with a known geometrical setup.

Data

We used a commercial heart model (PVAH-01 Medical Imaging Technologies, SHELLEY Automation Inc. London, Canada) limited to the ventricles. The size of this phantom (left ventricle volume: 76 mL and right ventricle volume: 66 mL) is in the same order as a woman heart. We modified this phantom in order to include a tubular structure that simulates the esophagus (Fig. 3). We opened four windows in the tubular structure to allow the US probe to acquire some images at certain defined positions. Unfortunately, only 3 windows were usable, the fourth one was directed on the "apex" part of the phantom with only little usable information for registration. For technical reasons the esophagus mimicking structure could not be set as close to the phantom as in the real anatomy. It was 4 -5 cm distant from the heart instead to be adjacent. However, this allowed to include all the phantom (and so more details) in the US image field of view. The heart and esophagus models were immersed into water filling the internal region of the heart cavities and the inside of the esophagus. The water ensured the US wave propagation to image the heart phantom appropriately. We also fixed three rubber markers on the outer surface of the heart cavities to provide landmarks enabling to compute a spatial reference in the images.

The TEE acquisitions were acquired on a Vivid Q 2D real-time cardiovascular echo system with a 6VT-D 4D TEE probe (General Electric, Healthcare, Wauwatosa, USA). On this probe we fixed the orientation angle of the acquisition plane to zero degrees. The preoperative volumes were acquired by a Cone Beam CT (CBCT) fixed on an Artis zeego C-Arm (Siemens Healthcare, Erlangen, Germany) (Fig. 4). Acquired ultrasound images had a dimension of 636 × 422 pixels with pixel size of 0.4316 × 0.4316 mm 2 . The CT volumes were composed by 512 × 512 × 372 voxels with a voxel size of 0.4864 × 0.4864 × 0.4864 mm 3 .

The TEE probe was inserted in the esophagus of the phantom and moved forward to reach the first window. On this probe position, we simultaneously acquired a 2D US image using the TEE and a 3D CT volume with the CBCT arm. The TEE probe can therefore be seen inside the esophagus at the US acquisition time in the CBCT volume. The cavities anatomy, the three rubber makers and the origin point of the TEE image can be seen for both 2D US and 3D CT. This acquisition procedure was repeated for the other two windows.

Esophagus-based volume V ther

The physical phantom CT image was processed as described in 2.1. A smoothed centerline was fitted from a set of 180 extracted points of segmented esophagus. From this esophagus centerline, the physical phantom CT data (Fig. 5-a andb) was processed to obtain the esophagus-based volume reconstruction (Fig. 5-c andd). It can be seen that on this new volume, the esophagus centerline is a straight line. All the axial slices are perpendicular to the esophagus centerline. The location of the first and last slice of the reconstructed data can be seen as blue lines in Fig. 5-a.

Ground truth

The direct estimation of the I U S 3D pose in V pre as well as the transformation matrix T Ref was performed by a human expert using the software Slicer [START_REF] Pieper | 3D slicer[END_REF]. Slicer allows to easily extract any oblique 3D slice from a volume and to superimpose the information from two different sources. The procedure consists of the following: i) The vertex of the I U S field of view is manually set at the center C of the TEE probe seen in V pre ; ii) The human expert performs a visual registration between the 2D US image and the CT volume. The 6 parameters of a 3D rigid transformation (3 translations and 3 rotations) that align the moving CT volume with the US image were manually tuned using the Slicer software. At each step, the superimposition of the US image and the transformed CT slice served as a visual similarity measure. This is used to estimate the correct alignment, especially the myocardium form and the position of the fixed markers; iii) The transformation matrix provided an estimate of T Ref after convergence. This process was repeated three times resulting in three transformation matrices. The final reference matrix T Ref was constructed using the mean parameters of these three matrices.

Registration

In order to prove the validity and accuracy of our concept, we chose to reduce the number of geometrical parameters to be estimated. The transformation T RecEval is estimated in three steps. i) For a specific I U S (Fig. 6-a), the slice of coordinate d i in V ther , which contains the segmented TEE tip center C, was considered the corresponding slice in the therapy-adapted preoperative volume (Fig. 6-b). On this slice, C was a common point between I U S and V ther . ii) The orientation of I U S in the corresponding slice was estimated by a rotation of I U S around C. This rotation was estimated using a 2D-2D fiducial registration procedure based on the rubber markers (Fig. 6-c). iii). The estimated 2D rotation matrix combined with T(d i ) gave the transformation matrix T RecEval . In order to integrate the segmentation errors, we compute T RecEval for the five slices around the presumed transducer center.

Results of the phantom study

As shown in Fig. 6, the myocardium form in the CT reconstruction is similar to the imaged myocardium in the US image. However, due to the speckle size, the acquired US has high boundary dispersion. The wall of the myocardium is perceived larger in the US image than in the CT volume. Moreover, the brightness of the US image seems to be too high. Although, it is difficult to see the alignment due to the quality of the US image, an error around the fiducial points is visible.

On each I U S the centers of the three rubber markers are chosen as the fiducial points P i used to measure the T RE (eq. 4). Since the rubber markers form a large spot on the US image, each fiducial point was manually selected ten times. This allows us to include intra-observers variability in our measures.

To summarize, each fiducial point P i , T RE(P i ) was estimated 50 times (5 slices and 10 manual selections). Fig. 7 shows the box plots of the T RE for each fiducial points for every TEE probe position. The mT RE are always lower than 5 mm. For all the three probe positions in the esophagus, we have a minimum error lower than 3.5 mm. The global mean error of all the fiducial markers in all TEE probe positions was measured as 4.03 mm.

Clinical study -Validation of 2D US to 3D CT registration

This study aimed to measure the accuracy of the proposed 2D US to 3D CT registration on real human images. To do this, we have a CT volume and a recorded TEE sequence of the same patient. This study demonstrated that the performance of the proposed registration scheme achieved clinical success with current technology and imaging protocols.

3.3.1. Data Three intraoperative US sequences and one preoperative 3D CT of a patient suffering from atrial fibrillation pathology were obtained. The US sequences were acquired from three different depths inside the esophagus using a conventional TEE probe. The rotation and inclination angles of the TEE probe was fixed on zero degrees in order to reproduce the same imaging condition. The datasets were provided by Louis Pradel University Hospital, Lyon, France.

The ECG-gated US acquisition was carried out on a Vivid-5 2D real-time cardiovascular echo system (General Electric, Healthcare, Wauwatosa, USA). One cardiac cycle was recorded as a sequence of 46 frames. The dimension of each frame was 636 × 420 voxels with an image spacing of 0.35278 × 0.35278 mm 2 .

An ECG-gated cardiac multislice CT image was acquired after injection of contrast agent with a Philips 64-slice scanner (Brilliance CT, Philips Healtcare) at 75% of cardiac cycle (R-R interval). The dimension of the reconstructed image was 512 × 512 × 395 voxels with an image spacing of 0.546875 × 0.546875 × 0.44998169 mm 3 .

Esophagus-based volume V ther

The clinical CT image was processed as described in 2.2. The result of the esophagus segmentation can be seen on Fig. 9-a. The Esophagus-based reformatted volume V ther is shown in Fig. 9-b.

Ground truth

The procedure to obtain the ground truth and T Ref is similar to the one followed for the phantom study. For a specific I U S , the medical expert (a senior radiologist) first marked some fiducial points in both 2D US and CT volume. This allowed an initial registration. The expert then aligned interactively and visually the 2D US and on the CT volume using the software Slicer. The six parameters of a rigid transformation (3 translations and 3 rotations) were also set manually until visual convergence. This process was repeated three times resulting in three transformation matrices. The final reference matrix T Ref was constructed using the mean parameters of these three matrices.

3.3.4. Registration For a specific I U S of the clinical study, we chose interactively in V ther an interval of 20 slices along the esophagus centerline as candidate zone. The registration was then performed as described in 2.2: 2D-2D rigid registration for each [START_REF] Tobon-Gomez | Benchmark for algorithms segmenting the left atrium from 3D CT and MRI datasets[END_REF].

On a classical PC without any optimization, the mean computation time for the registration process is 5.11 s with a standard deviation of 0.98 s for each reconstructed image.

3.3.5. Results of the clinical study For each of the three US images, some corresponding fiducial points were marked by our medical expert in both the US image (Fig. 11-a) and the 3D CT volume V pre (Fig. 11-b). The registration is then performed as described previously. The transformation matrix of the best 2D-2D registration result of the US image and the reformatted CT slice of V ther (Fig. 11-c) was used to project the fiducial points of I U S to V pre . This allowed to determine the T RE for each fiducial point (Fig. 11-d). The framework has been repeated 10 times for the three TEE positions.

The registration process was able to estimate an insertion length for each TEE position, 15.9, 62.1 and 78.9 mm respectively, from a specific reference entrance point O end . The box plots of the T RE at the three different depths inside the esophagus are shown in Fig. 12. For the two US images close to the insertion point of the TEE (Image 1 and 2), the mT RE was 3.5 mm, and for the third image which was farther from the insertion point, the mT RE increased to 9.5 mm.

The global mT RE including all the fiducial points is 5.6 mm.

Discussions

Phantom study -Accuracy of Esophagus-based volume reformatting

The error for the fiducial points in each window shows that the position of the reformatted slice is around the position of the corresponding slice of the ground truth. However, with the orientation sometimes towards the left or right, or similar to window 2, the slice is oriented upwards. Also, we can observe that the error of the points closest to the vertex (fiducial 1 and 2) are not the smallest for all windows as expected given the mechanical constrictions of the probe. The transesophageal probe is flexible and follows the anatomical curves, but the tip of the probe is rigid. So, we assume that the differences between the reconstructed slice and the TEE image should be due to light displacement of the tip extremity, rather than an inclination of the tip. As the error shows the inclination of the reconstructed slice but without a global tendency, we believe that the reconstruction is sensitive to the small changes in the centerline. One of the principal sources of error in the centerline extraction, that could produce these small changes, is the segmentation process of the esophagus in the phantom images.

The segmentation is challenging because the acquired images of the phantom present some artifacts introduced by the frame that supports the phantom. Also, the local morphological operations used to fill the holes of the windows can modify the real position of the centerline. However these are not sources of error in human data.

As displayed in Fig. 8, we can see that the error increases in relation to the depth. The bigger error corresponds with the probe that is farthest from the insertion place. We can conclude from the measures that if the depth increases, the probe moves away from the centerline. We believe that the probe tip gains some freedom, thus it can move away from the estimated esophagus centerline.

Clinical study -Accuracy of 2D US to 3D CT registration

From Fig. 12, we can observe that the error increases while the TEE tip is progressing in the esophagus. This phenomenon has already be reported in the phantom study. Some more accuracy could probably be gained by considering the estimated pose as a starting point for performing a 6 degrees of freedom registration around it. The consequence Concerning the registration accuracy, The global mT RE of 5.6 mm is on the same order of magnitude as those reported in the literature: less than 5 mm for [START_REF] Lang | Surface-based CT-TEE registration of the aortic root[END_REF]) and 1.5 -4.2 mm for [START_REF] Gao | Registration of 3D trans-esophageal echocardiography to X-ray fluoroscopy using image-based probe tracking[END_REF]. However, in those approaches, more information was available than in our case. For example Lang. et al. used a combination of two cross sections and/or 3D-TEE images to perform the registration and the target was segmented manually in the US images. In addition, a manual selection of the cross section was needed for a good initialization [START_REF] Lang | Surface-based CT-TEE registration of the aortic root[END_REF]. In their approach, Gao et al. used X-Ray fluoroscopy as an intraoperative intermediate imaging system combined with a 3D model of the US imaging tip. They also needed a 3D US probe and a manual pre-positioning of the 3D pose of the 3D model before registration [START_REF] Gao | Registration of 3D trans-esophageal echocardiography to X-ray fluoroscopy using image-based probe tracking[END_REF].

Beside the fact that the registration will give a 3D intraoperative view of the cardiac and extra cardiac structures to improve US navigation and to avoid sites where acoustic fences might avoid a proper focus of the US beam, its accuracy seems also to be in the range of the acceptable errors for the therapy. The goal of the therapy is to juxtapose thermal lesions to create a line of electrical isolation of a pathological zone. RF catheter ablation usually requires a 6 mm point by point inter-lesion distance to get a proper isolating line [START_REF] Wolf | Evaluation of left atrial linear ablation using contiguous and optimized radiofrequency lesions: the aline study[END_REF]. We can thus suppose that the correlation in between US and CT should need to be under 1 cm, but some dosimetric study on transesophageal HIFU and cardiac ablation must be conducted to confirm this assumption.

The primary limitation of our method is the esophagus segmentation step. It is actually performed using the semi-automatic active-contour segmentation method. The viability of our registration step is directly related to a precise V ther which depends on an accurate esophagus centerline estimation. We believe that the smoothing splines compensates for the segmentation uncertainties provided by the human expert. More tests should be performed to verify this hypothesis. On the other hand, a more automatic and accurate segmentation method could overcome this limitation.

On a standard PC (HP Z800, 8xIntel R Xeon R CPU E5620@2.40GHz, 24 GB RAM), without optimization, the time spend for the registration process is about 5.11 ± 0.98 s which could be considered as prohibitive for a real time intraoperative usage. This time is mostly spent by the 2D/2D registration procedure which could be speed up by the implementation of the latest optimized methods and the use of GPU (MI is a good candidate for parallel computing).

Conclusion

We have proposed a novel 2D-US to 3D-CT registration approach adapted to the guidance of the atrial fibrillation transesophageal HIFU therapy. Our approach takes advantage of the anatomical constraints imposed at the transesophageal HIFU probe to simplify the registration process. The method was successfully evaluated on a human dataset acquired from a patient with cardiac arrhythmia pathology. Our method exhibited level of performance and accuracy equivalent to previous approaches. However, compared to the other methods, our procedure requires only a light manual initialization (choice of the cardiac phase, very rough estimation of the probe along the esophagus) and no time consuming manual pre-positioning or manual selection of cross sections for proper initialization. Moreover, there is no need for any extra intermediate intraoperative external ionizing imaging device for the registration in our procedure.

Figure 1 .

 1 Figure 1. The general framework of our approach. During the preoperative stage, CT datasets are reformatted following the esophagus topology. During the intraoperative stage an intensity-based registration centered in the esophagus center is performed between the US image and reformatted CT images obtained in the previous stage.

Figure 2 .

 2 Figure 2. Coordinate systems of the preoperative and the therapeutic-adapted volumes

Figure 3 .

 3 Figure 3. Physical phantom of the heart and esophagus. Three windows are opened in the simulated esophagus for TEE acquisition. Two of three markers can be seen on the outer surface of the heart.

Figure 4 .

 4 Figure 4. CT volume of the phantom acquired with the Artis zeego C-Arm. Axial (a), coronal (b) and sagittal (c) views.

Figure 5 .

 5 Figure 5. Original data with segmented heart: a) mid-sagittal slice and b) a 3D view. Esophagus-based reconstructed data: c) mid-sagittal slice and d) a 3D view.

Figure 6 .

 6 Figure 6. a) US image acquired at the TEE probe position 1 from the esophagus phantom. b) Reformatted CT image from the volume obtained after the segmentation of the heart phantom, the markers and the tip of the probe at the TEE probe position 1. c) Result of the fiducial registration.

Figure 7 .

 7 Figure 7. Box plots of the Target Registration Error (T RE) at each fiducial points for each TEE probe position.

Figure 8 .

 8 Figure 8. Box plots of the Target Registration Error (T RE) for each TEE probe position.

Figure 9 .

 9 Figure 9. Esophagus-based volume V ther construction steps: a) Esophagus segmentation in V pre ; b) Esophagus-based reformatted volume V ther

Figure 10 .

 10 Figure 10. Registration results. a) superimposition of the 2D US to the corresponding 2D reformatted CT slice. b) final 2D US to 3D CT registration.

Figure 11 .

 11 Figure 11. Fiducial points marked in the a) US image and in the b) CT volume by the expert are presented. c) The best 2D-2D registration result of the US image and a reformatted CT slice. d) T RE for each fiducial point.

Figure 12 .

 12 Figure 12. Box plots of the Target Registration Error (T RE) at three different depths inside the esophagus

Acknowledgments

This work was part of the CardioUSgHIFU (ANR-11-TECS-0004) and CHORUS (ANR-17-CE19-0017) projects which have been supported by the French National Research Agency (ANR).