Marie Trin

THURSTON'S COMPACTIFICATION VIA GEODESIC CURRENTS: THE CASE OF NON-COMPACT FINITE AREA SURFACES

Keywords: Teichmüller space, Thurston's compactification, geodesic currents, sequences of random geodesics

In 1988, Bonahon gave a construction of Thurston's compactification of Teichmüller space using geodesic currents. His argument only applies in the case of closed surfaces, and there are good reasons for that. We present a variant which applies to surfaces of finite area and to do so we prove a control theorem for sequences of random geodesics. Note that this theorem may be of independant interest, especially when the surface is non-compact.

Introduction

The Teichmüller space T(S) of a surface S of finite topological type, with no boundary and of negative Euler characteristic χ(S) is the space of isotopy classes of (complete and finite volume) Riemannian metrics on S of constant curvature -1. Teichmüller space is not compact but Thurston showed in [START_REF] Thurston | On the geometry and dynamics of diffeomorphisms of surfaces[END_REF] how it can be compactified by the space P + ML(S) of projective measured laminations on S. The starting point of Thurston's compactification is the embedding of T(S) into the space P + (R C(S)

+ ) = R + R C(S) + : : T(S) → P + (R C(S) + ) X
→ R + X (•); Here X is the length function associated to the hyperbolic structure X on S and C(S) is the set of free homotopy classes of essential closed curves of S. What Thurston did is to prove that the image of is locally compact and to identify the boundary of T(S) in P + (R C(S) + ) with P + ML(S).

Theorem (Thurston's compactification). If S is a finite analytic type surface with negative Euler characteristic, then the accumulation points of T(S) in P + (R C(S) + ) are the projective classes of functions γ → i(λ, γ) where λ ∈ ML(S) is a measured lamination on S.

Thurston's original proof is explained in [START_REF] Fathi | Travaux de Thurston sur les surfaces[END_REF]. Some versions using realtrees are given by Morgan-Shalen [START_REF] Morgan | Valuations, trees, and degenerations of hyperbolic structures[END_REF], Bestvina [START_REF] Bestvina | R-trees in topology, geometry and group theory[END_REF] or Paulin [START_REF] Paulin | Topologie de Gromov équivariante, structures hyperboliques et arbres réels[END_REF]. An overview of the different compactification methods is availlable in [START_REF] Paulin | Sur la compactification de Turston de l'espace de Teichmüller[END_REF] or [START_REF] Ohshika | Compactifications of Teichmüller spaces iv[END_REF]. A compactification for the set of flat-structures and using geodesic currents is done in [START_REF] Duchin | Length spectra and degeneration of flat metrics[END_REF], note that this article is interested in both compact and non-compact surfaces. Here, we will be mostly interested in a very elegant argument, for closed surfaces, due to Bonahon [START_REF] Bonahon | The geometry of Teichmüller space via geodesic currents[END_REF]. Let's sketch the proof. Recall that geodesic currents are π 1 (S)-invariant Radon measures on the set of bi-infinite geodesics of the universal cover of S. Bonahon embeds T(S) into the space C(S) of geodesic currents of S, sending each element X ∈ T(S) of the Teichmüller space to the associated Liouville current L X ∈ C(S). The Liouville current satisfies two important properties: i(L X , γ) = X (γ) for every essential closed curve γ, and (1.1)

i(L X , L X ) = π 2 |χ(S)|. (1.2)
Here, i : C(S) × C(S) → R + is the intersection form, a continuous bilinear map extending the usual geometric intersection number between curves. Compactness of S implies compactness of the space P + C(S) of projective currents. It follows that each sequence (X n ) n∈N in Teichmüller space admits a subsequence, say the whole sequence, which projectively converges to a non-zero current µ, meaning that there are positive reals ε n such that lim n→∞ ε n L Xn = µ. The continuity of i and property (1.1) ensure that the length functions Xn (•) converge projectively to i(µ, •). Moreover, ε n tends to zero unless X n converges in T(S). Knowing that ε n ---→ n→∞ 0, property (1.2) ensures that i(µ, µ) = 0, meaning that µ is a measured lamination, as we needed to prove.

We stress that Bonahon's argument, with all its simplicity, only applies to closed surfaces. We will come back later to this specificity and to the obstructions to a direct extension of his argument. Recently, Bonahon and Šarić have given another proof of this theorem using geodesic currents. The arguments in [START_REF] Bonahon | A Thurston boundary for infinite-dimensional Teichmüller spaces[END_REF] are geared to infinite type surfaces, it is worth noticing that working in such a general context implies the lost of the simplicity of Bonahon's original proof.

Our goal here is to adapt Bonahon's original argument to be able to deal with non-compact surfaces of finite analytic type .

Let's look at the difficulties that prevent the extension of Bonahon's proof to the non-compact case. The intersection form, especially its continuity, is the linchpin of Bonahon's original proof. However, continuity fails when the surface is not compact, even if it has finite area (see [START_REF] Sasaki | Currents on cusped hyperbolic surfaces and denseness property[END_REF] or Example 2.1 below). We will therefore change our point of view to allow us to benefit from the continuity of i. We will consider currents on Σ instead of S, where Σ is a compact hyperbolic surface with geodesic boundary whose interior is homeomorphic to S, that is S = Σ \ ∂Σ. The second key ingredient of Bonahon's proof is the existence of the Liouville current but, as we will see, when working with currents on Σ we lose the Liouville current. Proposition 2.2. Let Σ be a compact hyperbolic surface with non-empty boundary and X a hyperbolic structure on S = Σ \ ∂Σ. There is no current L X on Σ which satisfies i(L X , γ) = X (γ) for every essential closed curve γ ∈ C(Σ).

In order to recover a version of properties (1.1) and (1.2), we will, for every hyperbolic structure X on S, replace the Liouville current L X by specific sequences of random geodesics (γ

(X) n ) n∈N ,
that is sequences of essential closed geodesics whose associated probability measures in T 1 X converge to the Liouville measure with respect to the weak- * topology. They will be chosen to satisfy (1.1) and (1.2) asymptotically, that is:

lim n→∞ i γ n X (γ n ) , γ = X (γ) π 2 |χ(S)| for all essential closed curve γ, (1.3) lim n→∞ i γ n X (γ n ) , γ n X (γ n ) = 1 π 2 |χ(S)| . (1.4)
As discussed in [START_REF] Erlandsson | Mirzakhani's curve counting and geodesic currents[END_REF], any sequence of random geodesics (γ n ) n∈N satisfies (1.3). Moreover, if the surface is compact then (1.4) is ensured for every sequence of random geodesics. However, for a non-compact surface, arbitrary sequences of random geodesics do not necessarily satisfy (1.4), see Example 3.2 below. Indeed, a large part of this article will be dedicated to building sequences of random geodesics satisfying this property for non-compact surfaces.

Theorem 1.1. For every complete and finite area hyperbolic structure X on a finite analytic type surface S of negative Euler characteristic χ(S), there is a sequence (γ (X) N ) n∈N of random geodesics such that:

lim n→∞ i γ (X) N X (γ (X) N ) , γ (X) N X (γ (X) N ) = 1 π 2 |χ(S)| .
Theorem 1.1 is actually part of a more technical result, Theorem 3.5, that we will prove in section 3. The main additional content of Theorem 3.5 is to ensure that the convergence rates in (1.3) and (1.4) hold with no dependance on the structure X. This uniformity will be important to achieve the proof of Thurston's compactification in section 4. Moreover, the proof of the theorem also ensures that we can control the behavior of sequences of random geodesics into some cusp's neighborhoods.

Acknowledgements. I am grateful to Juan Souto for our discussions and for all his suggestions about this paper. I would like to thank people who gave me their feedbacks on the first version of this document, especially Francis Bonahon, Frédéric Paulin, Beatrice Pozzetti, and Dylan Thurston for their remarks and Didac Martinez-Granado and Arya Vadnere for their questions. I also want to thank the PhD students of the IRMAR in ergodic theory for our debates on hyperbolic geometry, Barbara Schapira for her help with computations and Jing Tao for the time she gave me, her advices and her comments on my work. Finally, I thank the anonymous referee for their careful review of this document.

Preliminaries

In this section, we give some technical results and definitions. We refer the reder to [START_REF] Bonahon | The geometry of Teichmüller space via geodesic currents[END_REF], [START_REF] Bonahon | Bouts des variétés hyperboliques de dimension 3[END_REF] and [START_REF] Erlandsson | Mirzakhani's curve counting and geodesic currents[END_REF] for details. From now on, let S be a non-compact surface of finite analytic type, with negative Euler characteristic χ = χ(S) < 0. We denote by X, X , X n ... points in the Teichmüller space of S, or maybe just the underlying complete and finite area hyperbolic structure. Note that, although not specified, all the hyperbolic structures are complete and finite area. Moreover, we will write Z to refer indifferently to any finite area hyperbolic surface, possibly with punctures or with geodesic boundaries. If S is endowed with a hyperbolic structure X then every free homotopy class of essential closed curves contains a unique geodesic representative, so we identify a class with its geodesic representative when the hyperbolic structure is fixed. We will denote by C(S) the set of free homotopy classes of essential closed curves -by essential we mean non-null-homotopic and nonperipheral-or equivalently the set of essential closed geodesics. Let also Σ be a compact hyperbolic surface with geodesic boundary whose interior is homeomorphic to S. We fix a homeomorphism between S and Σ \ ∂Σ. This homeomorphism immediately induces a correspondance between the essential closed curves of S and the ones of Σ, that is

(2.1) C(S) = C(Σ).
The homeomorphism S = Σ \ ∂Σ also gives an identification between measured laminations of S and the ones of Σ supported by Σ \ ∂Σ:

(2.2) ML(S) = {λ ∈ ML(Σ)|λ supported by Σ \ ∂Σ}.

The identifications (2.1) and (2.2) will allow us to work on Σ rather than on S.

2.1. Currents on surfaces. We recall now a few properties of currents that we will need in the following. A geodesic current on Z is a π 1 (Z)invariant Radon measure on the set of bi-infinite geodesics on the universal cover Z of Z (even if the surface has non-empty boundary). The space C(Z) of geodesic currents on Z was introduced by Bonahon in [START_REF] Bonahon | Bouts des variétés hyperboliques de dimension 3[END_REF] and is endowed with the weak- * topology. For more information on currents we refer to [START_REF] Bonahon | Bouts des variétés hyperboliques de dimension 3[END_REF], [START_REF] Bonahon | The geometry of Teichmüller space via geodesic currents[END_REF], [START_REF] Aramayona | Hyperbolic Structures on Surfaces and Geodesic Currents[END_REF] and, [START_REF] Erlandsson | Mirzakhani's curve counting and geodesic currents[END_REF]Chap. 3].

The currents we will be mainly interested in are weighted multicurves and measured laminations and we will always consider currents on the compact surface Σ. An advantage of doing so is that when Z is compact, the topological space C(Z) is locally compact, and the associated projective space

P + C(Z) = R + C(Z) \ {0} is compact.
Moreover, in the compact case, the geometric intersection number between curves extends to a continuous bilinear map i : C(Z)×C(Z) → R + . It will be important later on to know that this form gives us a characterisation of the measured laminations as being the currents µ ∈ C(Z) satisfying i(µ, µ) = 0. We can also notice that the boundary curves are characterised by a zero intersection form with every current. As mentioned earlier, the reason why we want to work with the currents on the compact surface Σ, rather than with the currents on S, is that the continuity of the intersection number fails in the latter case. See [START_REF] Sasaki | Currents on cusped hyperbolic surfaces and denseness property[END_REF]Prop. 5.1] for a more detailed discussion on that obstruction to a continuous extension of the intersection number on the space of currents for non-compact surfaces.

Example 2.1 shows that there is no continuous extension of the intersection number for currents on S -it is the reason why we chose to work with currents on the compact surface Σ instead of the currents on S = Σ \ ∂Σ. This solves the problem of continuity of i(•, •) but raises a new problem: we won't be able to consider the Liouville current anymore. Proposition 2.2. Let Σ be a compact hyperbolic surface with non-empty boundary and X a hyperbolic structure on S = Σ \ ∂Σ. There is no current

L X on Σ which satisfies i(L X , γ) = X (γ) for every essential closed curve γ ∈ C(Σ).
Proof. If γ is a closed geodesic and µ a weighted multicurve of Σ then

(2.3) i(γ, µ) = min (γ ∩ µ), γ piecewise geodesic homotopic to γ in µ-general position ,
where a piecewise geodesic homotopic to γ is in µ-general position if the set of geodesics passing through the corners has vanishing µ measure.

, i(γ k , µ) ≤ k (b 1 ∩ µ) + k (b 2 ∩ µ) + 2 (γ ∩ µ) = 2 (γ ∩ µ).
We want to extend the previous inequality for µ a current, to do so we need a well-defined notion of intersection with γ. For this purpose we can embed Σ into the closed doubled surface DΣ, for more details about how to pass from Σ to DΣ the reader can refer to [START_REF] Trin | Counting arcs of the same type[END_REF]. Hence, C(Σ) is a subset of C(DΣ), the double γ of γ is a curve and in C(DΣ) we have i(γ k , µ) ≤ 2i(γ, µ), (2.4) for any µ weighted multicurve of Σ. Moreover, the weighted multicurves are dense in C(Σ) and the intersection number is continuous in C(DΣ) so Eq. (2.4) induces that

(2.5) ∀ν ∈ C(Σ), i(γ k , ν) ≤ 2i(γ, ν) < ∞.
However, lim k→∞ X (γ k ) = ∞ for any hyperbolic structure X on S, so Eq. (2.5) prevents any intersection with a fixed current to produce the length.

Cusps neighborhoods and intersection number.

Everything in the next section relies on a good understanding of the behaviour of geodesics in cusps. More precisely, if X is a hyperbolic structure on S then we denote by H i k the embedded horosphere of length 1/k around the i-th cusp. The horosphere H i k bounds the horoball B i k of area 1/k. We will refer to H i k and B i k as the horosphere and horoball of depth k. We also set X k the compact core of X bounded by the horospheres H i k and B k its complement:

X k = X \ i B i k , B k = i B i k . (2.6)
There is a direct link between the number of times a curve turns around a cusp and the depth it reaches [START_REF] Basmajian | Prime orthogeodesics, concave cores and families of identities on hyperbolic surfaces[END_REF]Prop. 3.4]. It follows that every curve that goes deep into a cusp has a large self-intersection number. To make this link more clear we recall a notion introduced in [8, Def. 2.6]: the peripheral self-intersection number. Definition 2.3. Let Z be a hyperbolic surface (compact or not) and recall that a peripheral subgroup of π 1 (Z) is nothing other than a cyclic subgroup generated by a non-essential closed curve. The peripheral selfintersection number i per (γ, γ) of γ ∈ C(Z) is the supremum over all maximal peripheral subgroups G ⊂ π 1 (Z) of the maximal number of times that the image of a lift γ of γ under Z → G Z meets itself transversely.

The peripheral self-intersection number is a topological invariant. It is thus independent of the metric on S, or more specifically, whether one considers the curves on S or on Σ. Moreover, for every compact subset K of Z \ ∂Z there is a upper bound for the peripheral self-intersection number of the closed geodesics contained in K. Conversely, for every N > 0 there is a compact subset K N of Z \ ∂Z that contains all the geodesics γ with i per (γ, γ) ≤ N [8, Lem. 2.7]. In the absence of boundary, one can easily quantify this property.

Lemma 2.4. Let X be a non-compact finite topological type surface with no boundary, and γ be an essential closed curve on X, this curve has support on X k if and only if i per (γ, γ) ≤ 4k.

Proof. If we think of the curves of π 1 (X) as deck transformations then a peripheral subgroup of π 1 (X) is a subgroup generated by a parabolic element. Let's study a given cusp C i , we can assume that the correspondence between X and H 2 is such that an associated maximal parabolic element is z → z + 1. In that case, H i k lifts to the horizontal line { (z) = k} and if γ is a closed geodesic of X then the number of times that the image of a lift γ under X → < z → z + 1 > X meets itself transversely is {n ∈ Z \ {0}|γ ∩ (γ + n) = ∅}. However, γ stays in X k around C i , if and only if its lifts stay below the line { (z) = k}, if and only if its lifts are half circles of radius at most k. Such a geodesic of H 2 meets at most 4k translations of itself (n = ±1, ±2... ± 2k). The same process applies for every cusps and then to every maximal parabolic subgroup and we obtain the lemma.

Construction of controled sequences of random geodesics

In this section we prove that for all non-compact hyperbolic surfaces of finite volume with no boundary there are sequences of random geodesics satisfying (1.4). However, we will first see with Example 3.2 that in the noncompact case not all the sequences of random geodesics have this property.

3.1. Sequences of random geodesics. As we saw in Proposition 2.2, the Liouville current does not exist anymore in our setting. However, for every (complete and finite area) hyperbolic structure X on S the Liouville measure on T 1 X still exists. Recall that the Liouville measure L X is the measure on the unit tangent bundle T 1 X, obtained by pushing forward the Haar measure on PSL 2 (R) and normalized so that L X (T 1 X) = 2π vol X (S) = 4π 2 |χ(S)|. We are going to consider geodesics approximating the Liouville measure in the following sense. Definition 3.1. A sequence (γ n ) n∈N of essential closed geodesics on X is a sequence of random geodesics if the associated probability measures converge to L X with respect to the weak- * topology, meaning that:

T 1 X f dγ n X (γ n ) -→ n→+∞ T 1 X f dL X 4π 2 |χ(S)| ,
for every f ∈ C 0 c (T 1 X) continuous and compactly supported function on T 1 X.

Remark. We will generally use the notation γ for the renormalisation γ X (γ)

.

The Birkhoff ergodic theorem, together with the ergodicity of the geodesic flow, implies the existence of such sequences of geodesics. We refer to [START_REF] Erlandsson | Mirzakhani's curve counting and geodesic currents[END_REF]Chap. 2] for some facts about sequences of random geodesics that we will use here.

The construction of the Liouville measure ensures that for a compact subsurface K of X we have L X (T 1 K) = 2π vol X (K). Then, if the boundary of K is smooth, the Portmanteau Theorem implies that for every sequence of random geodesics (γ n ) n∈N we have

X (γ n ∩ K) X (γ n ) -→ n→+∞ vol X (K) 2π|χ(S)| .
Applying this property to our compact core X k we have

(3.1) X (γ n ∩ X k ) X (γ n ) -→ n→+∞ vol X (X k ) 2π|χ(S)| , and hence, (3.2) X (γ n ∩ B k ) X (γ n ) -→ n→+∞ vol X (B k ) 2π|χ(S)| .
What is much more surprising is that sequences of random geodesics can also be used to compute lengths. More concretely, we have

(3.3) i(γ n , I) X (γ n ) -→ n→+∞ X (I) π 2 |χ(S)| ,
for every compact geodesic segment I in X. This property is basically due to Bonahon [6, Prop. 14], we also refer the reader to [8, Prop. 2.4] for details. A direct consequence of (3.3) is that we can use random geodesics (γ n ) n∈N to compute the length of any essential geodesic γ ∈ C(S):

(3.4) i(γ n , γ) X (γ n ) -→ n→+∞ X (γ) π 2 |χ(S)| .
Note that in this equation the curve γ is fixed. Meaning that a priori, equation (3.4) does not say anything about i(γ n , γ n ). However, for compact sets (3.3) holds uniformly. As a consequence, cutting the geodesics γ n into geodesic segments we have

(3.5) i γ n X (γ n ) , γ n|K X (γ n|K ) -→ n→+∞ 1 π 2 |χ| .
for K any fixed compact subsurface of X.

All those considerations about sequences of random geodesics apply to compact surfaces, hence, if S were compact, applying (3.5) to K = S, then we would immediatly have that every sequence of random geodesics satisfies (1.4). However, that is not necessarily true in general.

Example 3.2. First, note that an excursion of length into some B k i has between ke /2 -2 and 4ke /2 self-intersections. Consider now a sequence of random geodesics (γ n ) n∈N . Add to γ n an excursion of length 6 log( X (γ n )) at depth k n ---→ n→∞ ∞ and pull it tight into a new geodesic γ n . If we add the excursions in a well-chosen way (for example, gluing it at the deepest point of an excursion) then the (γ n ) n∈N are still random geodesics and

i(γ n , γ n ) X (γ n ) 2 ≈ i(γ n , γ n ) + k n X (γ n ) 3 ( X (γ n ) + 6 log( X (γ n ))) 2 ∼ +∞ i(γ n , γ n ) X (γ n ) 2 + k n X (γ n ) ---→ n→∞ ∞.
One can can also refer to the arguments in Lemma 3.4 below to prove that such sequences of random geodesics exist.

In [START_REF] Steven | Renewal theorems in symbolic dynamics, with applications to geodesic flows, non-Euclidean tessellations and their fractal limits[END_REF]Cor. 11.2] or [START_REF] Lalley | Distribution of periodic orbits of symbolic and Axiom A flows[END_REF], Lalley gives a construction of random geodesics that justifies the use of the term "random": if for all n the geodesic γ n is randomly chosen among the geodesics of length at most n then (γ n ) n∈N is a sequence of random geodesics with probability 1. Hence, we wonder which proportion of sequences of random geodesics satisfies (1.4). This problem might be linked to the study of the length of cusp excursions for random geodesics, see for example [START_REF] Haas | Geodesic cusp excursions and metric diophantine approximation[END_REF], [START_REF] Pollicott | Limiting distributions for geodesics excursions on the modular surface[END_REF] or [START_REF] Sullivan | Disjoint spheres, approximation by imaginary quadratic numbers, and the logarithm law for geodesics[END_REF] and the references therein. Anyway, the above example makes clear that to obtain (1.4) in the noncompact case we have to control the excursions of the sequences of random geodesics into cusps neighborhoods. We will do it through the cutting process described below.

3.2. Cutting process. Suppose that X is a fixed complete and finite area hyperbolic structure for S. Recall that X t denotes the compact core of X bounded by the horospheres of length 1/t around the cusps of S and that B t = X \ X t is its complement. Given two parameters k ∈ N and 0 < θ < π/4, and a curve γ we want to cut the excursions of γ in B k in order to prevent γ from leaving X k/ sin(θ) . To do so, we will study γ through its lifts in the universal cover X of X. We focus here on a given cusp but we apply the same construction around each cusps of X. For t ≥ 1 we denote by H t the horosphere of depth t around this cusp and B t the horoball it bounds. Since X is a hyperbolic surface endowed with a complete hyperbolic metric, its universal cover identifies with H 2 , and we can suppose that the parabolic element associated to the cusp we are interested in is z → z + 1. With this normalization H t lifts to the horizontal line { (z) = t} and we have that if a curve enters H t with some angle α ∈ [0, π/2) then it reaches the horosphere H k/ sin(α) (we measure the non-oriented angle with the normal to the horosphere). We want to cut γ in order to replace its long excursions into B k (ie. the ones which cross H k/ sin(θ) ) by short ones (excursions staying between H k/ sin(2θ) and H k/ sin(θ) ). To make it explicit we make a description of the process on the universal cover.

Figure 3. Cutting process

If γ makes excursions in B k we are going to modify γ explaining the process on a fixed lift γ which makes an excursion in the horoball { (z) > k} bounded by { (z) = k} but the same process applies to all lifts of B k . First, if γ enters with an angle greater than θ then we don't change it. On the other hand, if it enters with an angle smaller than θ then we replace this arc by a geodesic arc I which enters with angle between θ and 2θ and whose exit point coincides with the exit point of a different lift γ of γ (see Fig. 3). This is always possible as long as 2k cotan(θ) -2k cotan(2θ) ≥ 1. If we apply the same process to all the excursions of γ around every cusp then γ is replaced by a closed piecewise geodesic γ . Now, pulling γ tight we obtain a closed geodesic γ * : we refer to γ * as the geodesic obtained by cutting process of parameters k and θ from γ. Note that if θ is small then γ and γ * have basically the same length, more precisely, they can be mapped one to each other through a homotopy with small displacement and without disturbing to much the lengths. For the lengths, it is easy to see that there is some e θ ---→ θ→0 0, idependent from X, such that for every k ≥ 1 and θ small (3.6) X (γ ) ≤ (1 + e θ ) X (γ * ). Here X (γ ) refer to the arc length of γ , we will use again this abuse of notation but its meaning is clear from the context.

Construction of controled sequences of random geodesics. Lemma 3.3.

There is some θ 0 > 0 such that if (γ n ) n∈N is a sequence of random geodesics on X and (γ * n ) n∈N is obtained from the γ n applying the cutting process of parameters k > 1 and θ > θ 0 then there is

µ n ---→ n→∞ 0 such that 1 ≤ X (γ n ) X (γ * n ) ≤ (1 + µ n ) vol X (S) vol X (X k ) (1 + e θ ),
for every n.

Here, e θ is as in (3.6).

Proof. We use the same notation as in the description of the cutting process, and, as above, we denote by X (γ n ) the arc length of the piecewise geodesics. We take θ 0 small enougth such that (3.6) occurs. The γ n being random geodesics, (3.1) ensures that we can find a sequence µ n ---→ n→∞ 0 such that

X (γn) X (γ n|X k ) = (1 + µ n ) vol X (S) vol X (X k ) . The construction of γ n ensures that γ n|X k = γ n|X k , thus X (γ n|X k ) X (γ n ) ≤ 1 and if θ > θ 0 then X (γ n ) X (γ * n ) ≤ (1 + e θ ).
The upper bound follows from those three inequalities. Now, γ n and γ n coincide on X k but γ n has shorter excursions than γ n in B k , hence, X (γn) X (γ n ) ≥ 1. The geodesic γ * n is the unique geodesic representative of the free homotopy class of γ n which proves that X (γ n ) X (γ * n ) ≥ 1 and the lower bound follows.

Lemma 3.4. Let (γ n ) n∈N be a sequence of random geodesics. If (γ * n ) n∈N is obtained from (γ n ) n∈N applying the cutting processes of parameters k n ---→ n→∞ ∞ and θ n ---→ n→∞ 0, then (γ * n ) n∈N is a sequence of random geodesics.
Proof. In this proof, we denote by γ the canonical lift of a geodesic γ to the unit tangent bundle of X.

Let f ∈ C 0 c (T 1 X) be a continuous and compactly supported function on

T 1 X, there is K a compact core of X such that Supp(f ) ⊂ T 1 K. Since k n ---→ n→∞ ∞ then there is n 0 ∈ N such that for all n ≥ n 0 , γ n|K = γ n|K .
The homotopy between γ n and γ * n induces that the arcs of γ n|K are freely homotopic to geodesic arcs of γ * n . Such a homotopy induces a projection from γ n|K to γ * n and lifts to Ψ n : γn|K → γ * n , which is a homeomorphism on its image. The homotopy can be chosen to have low displacement, that is d(p, Ψ n (p)) ≤ ε n ---→ n→∞ 0 for every p ∈ γn|K , and not to distort too much the lengths. Moreover, we can find

ϕ n : [0, X (γ n|K )] → R + a piecewise smooth reparametrization of [0, X (γ n|K )] such that for all t ∈ [0, X (γ n|K )], Ψ n (γ n|K (t)) = γ * n (ϕ n (t)
). The homotopy between γ n and γ * n does not distort too much the lengths, hence, we have some δ n ---→ n→∞ 0 such that 1 -δ n ≤ ϕ n ≤ 1 + δ n where it is defined.

Fix some µ > 0. A compactly supported continuous function is uniformly continuous, thus, there is ε µ > 0 such that if d(p, q) ≤ ε µ then |f (p)-f (q)| ≤ µ. We can suppose that for every n ≥ n 0 , ε n ≤ ε µ . We have

T 1 X f dγ * n = X (Ψn(γ n|K )) 0 f • γ * n (t)dt = X (γ n|K ) 0 f • γ * n (ϕ n (s))ϕ n (s)ds, it follows that (1 -δn) X (γ n|K ) 0 f (Ψn(γ n|K (s))ds ≤ T 1 X f dγ * n ≤ (1 + δn) X (γ n|K ) 0 f (Ψn(γ n|K (s))ds ⇒ (1 -δn)    T 1 X f dγn -µ X (γ n|K )    ≤ T 1 X f dγ * n ≤ (1 + δn)    T 1 X f dγn + µ X (γ n|K )    ⇒ (1 -δn) X (γn) X (γ * n )    T 1 X f dγn -µ    ≤ T 1 X f dγ * n ≤ (1 + δn) X (γn) X (γ * n )    T 1 X f dγn + µ   
Adapting the proof of Lemma 3.3 we have

X (γ n ) X (γ * n ) ---→ n→∞
1, and passing to the limit in n we obtain

T 1 X f dL X 4π 2 |χ(S)| -µ ≤ lim n T 1 X f dγ * n ≤ lim n T 1 X f dγ * n ≤ T 1 X f dL X 4π 2 |χ(S)| + µ.
This is true for all µ, hence, lim

n→∞ T 1 X f dγ * n = T 1 X f dL X 4π 2 |χ(S)|
and we have proved that (γ * n ) n∈N is a sequence of random geodesics. Now, for every hyperbolic structure X on S, we will be able to build sequences (γ (X) n ) n∈N of random geodesics satisfying (1.4). Moreover, we will build them in such a way that neither the converging rates in (1.4) and (3.4), nor the peripheral self-intersection numbers i per (γ

(X) n , γ (X) n ) depend on X.
Theorem 3.5. For every complete and finite area hyperbolic structure X on a finite analytic type surface of negative Euler characteristic S, there is a sequence (γ (X) n ) n∈N of random geodesics such that :

lim n→∞ i γ (X) n X (γ (X) n ) , γ (X) n X (γ (X) n ) = 1 π 2 |χ(S)| .
More precisely, they can be chosen such that

(1) i(γ

(X) n , γ(X) n ) ≤ 1 π 2 |χ(S)| 1 + 1 n , ∀n ∈ N, (2) ∀α ∈ C(S), ∃n α ∈ N : i(γ (X) n , α) X (α) π 2 |χ| -1 -1 ≤ 3 n , ∀n ≥ n α , ( 3
) i per (γ (X) n , γ (X)
n ) ≤ C n , ∀n ∈ N, where C n and n α do not depend on X.

Proof. To obtain the desired sequence (γ (X) n ) n∈N we start with an arbitrary sequence of random geodesics (γ n ) n∈N . For every p we set k p = e p/2 and θ p = e -p/2 , if we apply the cutting process with parameters k p and θ p to the sequence (γ n ) n∈N then we obtain a sequence (γ p n ) n∈N of piecewise geodesics and by pulling it tight a sequence (γ p n ) n∈N of geodesics. We will chose the (γ (X) N ) N ∈N among the γ p n . First, study the self-intersection number of those γ p n . As γ p n is the geodesic representative of γp n , its self-intersection number is lower than the number of self-intersections of γp n . To count it, we divide X into two parts, the compact core X k and its complement B k . On X k , the geodesic arcs γp n|X k and γ n|X k are identical so γp n has i(γ n|X k , γ n ) self-intersections. On the complement, we count the self-intersections of γp n considering its different excursions in

B k : i(γ p n , γ p n ) ≤ i(γ n ∩ X k , γ n ) + I,J excursions in B k i(I, J).
We can distinguish two types of pairs (I, J): the ones where at least one of the excursions stays in B k ∩ X k/ sin(2θ) , and the ones where both I and J reach B k/ sin(2θ) . In the first case, I and J meet at most as many times as the corresponding excursions of γ n and then:

i(γ p n , γ p n ) ≤ i(γ n ∩ X k/ sin(2θ) , γ n ) + I,J excursions in B k which reach B k/ sin(2θ)
i(I, J).

Moreover, an excursion of γp n in B k which reaches B k/ sin(2θ) has a length of at least ln(1/θ), a lower bound for the length of the geodesic arc which enters with angle 2θ. It follows that there is at most X (γn∩B k ) ln(1/θ) such excursions. Also, the intersection number of two excursions reaching B k/ sin(2θ) is at most 4k/θ, the self-intersection number of the excursion which enters with angle θ. All in all,

i(γ p n , γ p n ) ≤ i(γ n ∩ X k/2θ , γ n ) + X (γ n ∩ B k ) ln(1/θ) 2 4k
θ .

Applying equations (3.5) and (3.2) we have

i(γ n ∩ X k/2θ , γ n ) = (1 + ε p n ) X (γ n ) X (γ n ∩ X k/2θ ) π 2 |χ| and X (γ n ∩ B k ) = (1 + δ p n ) X (γ n ) 2π|χ| C k
where C is the number of cusps of S, where ε p n ---→ n→∞ 0 and δ p n ---→ n→∞ 0 depend on X. As a consequence,

i(γ p n , γ p n ) ≤ (1 + ε p n ) X (γ n ) 2 π 2 |χ| + (1 + δ p n ) C X (γ n ) 2π|χ| • k • ln(1/θ) 2 4k
θ and we obtain a upper bound for the self-intersection number of the normalized curves:

i γ p n X (γ p n ) , γ p n X (γ p n ) ≤ 1 π 2 |χ| (1 + ε p n ) + (1 + δ p n ) 2 C 2 |χ| 4 p 2 X (γ n ) X (γ p n ) 2 . (3.7)
We next study the intersection number of the γ p n with closed curves. The set C(S) is infinite and can be enumerated with C(S) = {α q |q ∈ N} in such a way that i per (α q , α q ) ≤ 4q for every q. This enumeration is fixed whatever the structure X. Recall that for every p we have k = k p = e p/2 and θ = θ p = e -p/2 . Hence, when p is big enough, for q ≤ p the curve α q is included in X q ⊂ X k . However in

X k we have i(γ n , •) = i(γ p n , •) thus (3.8) i(γ p n , α q ) = X (γ n ) X (γ p n ) i(γ n , α q ).
Now, applying Lemma 3.3, for every p there is µ p n ---→ n→∞ 0, depending on X, such that

(3.9) 1 ≤ X (γ n ) X (γ p n ) ≤ (1 + µ p n ) vol X (S) vol X (X k ) (1 + e p ).
with e p = e e -p/2 with the notation of (3.6). Therefore, there are m p large enough such that ε p mp , δ p mp , µ p mp ≤ 1 p , and i(γ mp , α q ) X (α q )/π 2 |χ| -1 ≤ 1 p for every q ≤ p. Thus (3.9) and (3.7) give us

1 ≤ X (γ mp ) X (γ p mp ) ≤ vol X (S) vol X (X k ) (1 + 1 p )(1 + e p ) ---→ p→∞ 1, (3.10) (3.11) i(γ p mp , γp mp ) ≤ 1 π 2 |χ| 1 + (1 + 1 p ) 4C 2 p 2 |χ| vol X (S) vol X (X k ) (1 + 1 p ) 2 (1 + e p ) ---→ p→∞ 1 π 2 |χ| .
The terms on the right in inequalities (3.10) and (3.11) do not depend on X anymore so, for N an integer there is p N , independent from X and with p

N > p N -1 , such that 1 ≤ X (γ p N ) X (γ p N mp N ) ≤ 1 + 1 N and i(γ p N mp N , γp N mp N ) ≤ 1 π 2 |χ| (1 + 1 N ).
As a consequence, we can take γ (X)

N = γ p N mp N . The previous constructions ensure that i(γ (X) N , γ(X) N ) ≤ 1 π 2 |χ| (1 + 1 N
), and we have proved [START_REF] Aramayona | Hyperbolic Structures on Surfaces and Geodesic Currents[END_REF] in the statement of the theorem. and a sequence (ε k ) k∈N of positive real numbers such that lim k→∞ ε k k (•) = F pointwise (we have written k for X k ). We will prove that F is given by taking the intersection number with a suitable measured lamination.

Fix a filling curve β on S, that is a closed curve such that the connected components of S \ β are balls and annular neighborhoods of the cusps. Such a curve gives us a bound on the length of every curve γ ∈ C(S), namely, (4.1)

X (γ) ≤ X (β)i(γ, β)(1 + i(γ, γ)) for every hyperbolic structure X [19, Lem. 2.1]. Since F = lim k→∞ ε k k (•)
is non-zero, there is γ ∈ C(S) with F (γ) = 0. We obtain from (4.1) that 0 < F (γ) ≤ F (β)(1 + i(γ, γ))i(γ, β) and hence that F (β) = 0. Since we are only interested in convergence in P + (R

C(S)
+ ), we can assume that F (β) = 1, meaning that

lim k→∞ δ k k (•) π 2 |χ| = F,
where δ k = π 2 |χ| k (β) . We will now prove that F is of the form i(µ, •) where µ is a measured lamination on S.

Applying Theorem 3.5 to each X k , we obtain some sequences of essential closed geodesics (γ We want F to be the pointwise limit of (i(µ n , •)) n∈N . To prove it, it is sufficient to show that the convergence δ k n i(γ

Figure 1 .

 1 Figure 1. Obstruction to the continuity of i Example 2.1 (Discontinuity of the intersection form in the non-compact case). Take a hyperbolic surface with at least two cusps, fix an embedded horocycle around each of them, and a simple geodesic arc between those curves which meet them orthogonally. Note that this arc is part of a cuspsto-cusps geodesic arc γ. Consider a sequence of closed curves (γ n ) n∈N , where γ n is the geodesic homotopic to the closed curve which runs the geodesic arc mentioned above, turns n times around the first cusp following the fixed horocycle, goes back along the geodesic arc and turns n times around the second cusp as in Fig.1. The self-intersection number of such a sequence is going to grow without bound. On the other hand, it approaches the weight 2 current associated to γ which has 0 self-intersection number.

Figure 2 .

 2 Figure 2. Obstruction to the existence of the Liouville current Now, consider b 1 and b 2 two boundary components of Σ, maybe the same, and γ a non-trivial geodesic arc joining them. For every k, we define γ k as the unique closed geodesic homotopic to the piecewise geodesic which follows γ, turns k times around b 1 , follows back γ and turns k times around b 2 . We obtain from Eq. (2.3) that for any weighted multicurve µ,

3 n

 3 n ), •) = k (•)/π 2 |χ|. As all along, let Σ be a compact complete hyperbolic surface with boundary whose interior is homeomorphic to S and let's identify C(S) with C(Σ). In particular, we can consider the weighted curves γ(k) n ) as currents of Σ. The space P + C(Σ) being compact each (γ(k) n ) k∈N projectively converges to a non-zero current µ n ∈ C(Σ).We first want to show that the µ n are measured laminations. Consider the sequence (γ (k) n ) k∈N for n fixed, there are some ε k n > 0 such that ε k n γ(k) n tends to µ n up to a subsequence in k. So, by diagonal extraction we can suppose that ε k n γ(k) n ---→ k→∞ µ n for every n. What we have to show is that lim k→∞ ε k n = 0 for every n. The sequence (X k ) k∈N leaves every compact set of T(S) so there is a simple closed curve α such that lim k→∞ k (α) = ∞. Recall that to prove Theorem 3.5 we have enumerated C(S) = {α n |n ∈ N} such that i per (α n , α n ) ≤ 4n, since α is a simple curve we can suppose that α = α 1 . The γ (k) n come from Theorem 3.5 thus i(γ (k) n , α) k (α) whatever k and n. By hypothesis k (α) ---→ k→∞ ∞ and we can suppose, up to a shift in n, that for every n, i(γ

  n , α) ---→ k→∞ ∞. However, ∞ > i(µ n , α) = lim k→∞ ε k n i(γ (k) n , α)

n 1 - 1 ≤

 11 , γ) ---→ n→∞ δ k k (γ) π 2 |χ| is uniform in k when γ ∈ C(S) is fixed. If γ ∈ C(S)is fixed then Theorem 3.5 ensures that δ k n δ k -1 ≤ n and n for every k and for n large enough (n γ and n β do not depend on k) with n ---→ n→∞ 0. Moreover, fixing γ we know thatδ k k (γ) π 2 |χ| ---→ k→∞ F (γ) hence the sequence (δ k k (γ) π 2 |χ|) k∈N is bounded by some d γ and we obtainδ k n i(γ (k) n , γ) -δ k k (γ) π 2 |χ| ≤ v n d γ ---→ n→∞ 0.Hence, the convergence holds uniformly in k, and limn→∞ lim k→∞ δ k n i(γ (k) n , γ) = lim k→∞ lim n→∞ δ k n i(γ (k) n , γ), which implies that F (γ) = lim n→∞ i(µ n , γ). Moreover, ML(S) is a closed subset of R C(S) + hence F (•) = lim n→∞ i(µ n , •) is of the form F (•) = i(µ,•) where µ ∈ ML(S), which was what we needed to prove.

The author is supported by Centre Henri Lebesgue, program ANR-11-LABX-0020-

Applying Proposition 2.4 we have i per (γ (X) N , γ (X) N ) ≤ 4e p N where p N does not depend on X, which gives us the third point. At last, (3.8) and the choice of p N and m p induces that

hence, we obtain the second point with n α = q when α = α q . Moreover, up to passing to a subsequence, the (γ

0. As a consequence, Lemma 3.4 ensures that we have built a sequence of random geodesics. At last, for K a compact subsurface of X we have

and if we pass to the limit, using (3.5), we obtain that lim

Proof of Thurston's compactification

Armed with Theorem 3.5, we are now able to prove Thurston's compactification. As we already mentioned in the introduction, the starting point of this compactification is the embedding of T(S) and P + ML(S) into the space P + (R

The image of T(S) in P + (R

+ ) is included into a compact set (use Eq. (4.1) for instance), thus, the closure T(S) of T(S) is compact. The boundary of this set is given by the following theorem.

Theorem (Thurston's compactification). If S is a finite analytic type surface with negative Euler characteristic then the accumulation points of T(S) in P + (R C(S) + ) are the projective classes of functions γ → i(λ, γ) where λ ∈ ML(S) is a measured lamination on S.

Our arguments apply to the compact case, but for the sake of concreteness we will focus on non-compact surfaces.

Let X k ∈ T(S) be a sequence which converges in P + (R

C(S)

+ ) and leaves all compact sets of T(S), meaning that there are a non-zero element F of R C(S) + thus ε k n ---→ k→∞ 0 for every n, and i(γ

n , γ(k) n ) = 0 and µ n is a measured lamination on Σ. By construction, i per (γ

n ) ≤ C n for every k and n, as mentioned earlier (or in [START_REF] Erlandsson | Mirzakhani's curve counting and geodesic currents[END_REF]Lem. 2.7]) it ensures that for n fixed the γ (k) n are all included in the same compact subsurface of Σ \ ∂Σ. It follows that µ n is supported on a compact set of Σ \ ∂Σ and by (2.2) it is a measured lamination of S.

Recall that β is a filling curve of S, as a consequence, i(µ n , β) = 0 and hence, we can suppose that i(µ n , β) = 1 for every n and we obtain

where

is well-defined.

To sum up, we have the following convergence diagram, where all the convergences are pointwise.

n , •)

1 , •) δ 2 1 i(γ


Marion Savanier
email: marion.savanier@gmail.com

Emilie Chouzenoux

Jean-Christophe Pesquet

Cyril Riddell
email: cyril.riddell@ge.com

Jean-Christophe Pesquet

J.-C Pesquet

Deep Unfolding of the DBFB Algorithm with Application to ROI CT Imaging with Limited Angular Density

Keywords: region-of-interest, computed tomography, angular sub-sampling, deep unfolding, forward-backward, iterative reweighted scheme

come   L'archive ouverte pluridisciplinaire

I. INTRODUCTION

CT imaging is commonly used for diagnostic purposes and image guidance in interventional radiology and surgery. In interventions such as follow-up examinations of deployed stents and needle biopsies, only a small region of the patient is of interest. Irradiating only the ROI by X-rays involves focusing the X-ray beam with collimation techniques using radio-opaque blades before passing through the patient [START_REF] Chityala | Region of interest (ROI) computed tomography[END_REF]. A focused irradiation results in a substantial reduction in patient dose [START_REF] Yu | Radiation dose reduction in computed tomography: techniques and future perspective[END_REF] and truncated measurements (or projections). The inverse problem of reconstructing an ROI from a set of truncated projections is ill-posed. In practice, it is often combined with angular sub-sampling for fast or further dosesaving acquisitions.

Previous research demonstrated the ability of model-based iterative reconstruction (MBIR) for CT reconstruction from few-view measurements [START_REF] Sidky | Image reconstruction in circular cone-beam computed tomography by constrained, total-variation minimization[END_REF], [START_REF] Bian | Evaluation of sparse-view reconstruction from flat-panel-detector cone-beam CT[END_REF]. MBIR implements nonlinear iterative algorithms aiming at minimizing a penalized cost function. In CT, most works rely on total variation (TV) regularization [START_REF] Rudin | Nonlinear total variation based noise removal algorithms[END_REF]. Although these works achieve a reduction in angular sub-sampling artifacts compared to analytical reconstruction methods, truncated projections still challenge MBIR [START_REF] Langet | Compressed-sensing-based content-driven hierarchical reconstruction: Theory and application to C-arm cone-beam tomography[END_REF], especially in terms of computation time.

To implement MBIR, one must choose the reconstruction grid, i.e., the support of the reconstructed area. When the grid matches the support of the ROI, the data will not agree with the reprojection of the ROI. When the reconstruction grid includes the support of the entire object [START_REF] Yu | Compressed sensing based interior tomography[END_REF], the reconstruction becomes computationally expensive and less stable due to the increase of unknowns for the same amount of data. Truncated data only allows a rough estimation of the exterior anatomical background, which holds no clinical value. Thus, one more practical solution is to consider an intermediate smaller grid size with a "margin" outside of the ROI [START_REF] Paleo | Efficient implementation of a local tomography reconstruction algorithm[END_REF]. This achieves a faster and more stable ROI reconstruction in general. Yet, when dense objects such as metallic cables or needles are outside the reconstruction grid, the reprojection of the extended ROI contains high-frequency errors. Moreover, when too few projections are used for reconstruction, such objects suffer from aliasing, and additional streak artifacts can degrade the reconstructed ROI [START_REF] Langet | Compressed-sensing-based content-driven hierarchical reconstruction: Theory and application to C-arm cone-beam tomography[END_REF]. Another approach is to reconstruct the entire object with large voxels and then subtract the reprojection of the exterior from the data before reconstructing only the ROI from the subtracted data [START_REF] Hamelin | Design of iterative ROI transmission tomography reconstruction procedures and image quality analysis[END_REF], [START_REF] Ziegler | Iterative reconstruction of a region of interest for transmission tomography[END_REF]. This approach produces a low-frequency approximation of the exterior of the ROI. However, such an approximation is poor in the presence of dense objects in the exterior of the ROI, and unwanted high-frequency content remains after subtraction that must again be dealt with. Apart from the choice for the reconstruction grid, a way of reducing reconstruction time in MBIR is through preconditioning techniques. However, these techniques are often too restrictive for CT reconstruction when proximal algorithms are used [START_REF] Savanier | Unmatched preconditioning of the proximal gradient algorithm[END_REF]; hence the emergence of heuristics to accelerate TVbased methods. In [START_REF] Langet | Compressed-sensing-based content-driven hierarchical reconstruction: Theory and application to C-arm cone-beam tomography[END_REF], for example, substantial acceleration of the forward-backward algorithm is reached thanks to two modifications. The first modification involves weighing a leastsquares data fidelity term by the ramp filter. Each gradient step then performs an approximate inversion of the forward projection, and each proximity step acts as a post-processing filter. The use of an approximate inverse of the measurement operator in the data fidelity term is also advocated in [START_REF] Tirer | On the convergence rate of projected gradient descent for a back-projection based objective[END_REF]- [START_REF] Nilchian | Iterative FBP for improved reconstruction of X-ray differential phasecontrast tomograms[END_REF]. The second modification is the linear decrease of the regularization parameter strength along with the iterations. Despite these modifications, the proximity operator of TV regularization does not have a closed form, so the method still requires many sub-iterations. Thus, designing MBIR methods with a low computational cost for CT imaging is still challenging. An increasingly popular alternative to MBIR is convolutional neural networks (CNN) due to their increased expressivity and fast inference. CNNs, and in particular U-net [START_REF] Ronneberger | U-Net: Convolutional Networks for Biomedical Image Segmentation[END_REF], have already been used for removing sub-sampling streaks in CT reconstructions obtained using analytical methods [START_REF] Han | Deep residual learning for compressed sensing CT reconstruction via persistent homology analysis[END_REF], [START_REF] Jin | Deep convolutional neural network for inverse problems in imaging[END_REF]. However, there are concerns about the lack of guarantees and capacity for generalization of post-processing CNNs, because these networks do not ensure data consistency [START_REF] Sidky | Do CNNs solve the CT inverse problem[END_REF]. The deep unfolding paradigm [START_REF] Monga | Algorithm unrolling: Interpretable, efficient deep learning for signal and image processing[END_REF], [START_REF] Wang | Deep learning for tomographic image reconstruction[END_REF] circumvents this issue by offering a way to include a priori information in a neural network. Unfolding networks have been applied to many inverse problems, such as denoising [START_REF] Wang | Proximal Deep Structured Models[END_REF], [START_REF] Le | The faster proximal algorithm, the better unfolded deep learning architecture ? The study case of image denoising[END_REF], deblurring [START_REF] Bertocchi | Deep unfolding of a proximal interior point method for image restoration[END_REF], MRI reconstruction [START_REF] Yang | Deep ADMM-Net for Compressive Sensing MRI[END_REF], and CT reconstruction from fewview data [START_REF] Adler | Solving ill-posed inverse problems using iterative deep neural networks[END_REF]. Unfolding consists of untying each iteration of an optimization algorithm for MBIR, defining a set of learnable parameters, and training each iteration (or layer) in an end-to-end manner. Some authors allow the learning of the optimization algorithm hyperparameters [START_REF] Bertocchi | Deep unfolding of a proximal interior point method for image restoration[END_REF] as well as linear operators in the regularization, such as convolution kernels in ISTA-net [START_REF] Zhang | ISTA-Net: Interpretable optimization-inspired deep network for image compressive sensing[END_REF] and in [START_REF] Pustelnik | Wavelet-Based Image Deconvolution and Reconstruction[END_REF]. Others use CNNs to replace proximity operators, as in PD-net [START_REF] Adler | Solving ill-posed inverse problems using iterative deep neural networks[END_REF] and ADMMnet [START_REF] Yang | Deep ADMM-Net for Compressive Sensing MRI[END_REF]. Deep unfolding networks automatically inherit from the feedback mechanism of MBIR for data consistency. In this paper, we design a deep unfolding network for addressing the problem of ROI image reconstruction from few-view truncated measurements. Specifically:

• We first introduce an original cost function involving a Cauchy fidelity term and a semi-local total-variation regularization to limit sub-sampling streaks from objects inside and outside the reconstruction grid. • We then advocate for an iterative optimization algorithm combining for the first time an instance of the dual block forward-backward algorithm (DBFB) [START_REF] Abboud | Dual Block Coordinate Forward-Backward Algorithm with Application to Deconvolution and Deinterlacing of Video Sequences[END_REF] with an iterative reweighted scheme. • We propose a neural network architecture inspired by this algorithm, allowing supervised parameter learning and fast reconstruction on the GPU. • We evaluate the performance and generalization of our network on three datasets against state-of-the-art methods such as MBIR, other deep unfolding networks, and classical deep learning post-processing.

This paper is organized as follows: Section II introduces our notation, and Section III provides a mathematical formulation of our reconstruction problem. Section IV presents our cost function and introduces a convergent iterative algorithm to minimize it. Section V then explains how this algorithm is unfolded into a deep learning architecture and discusses our training strategy. This is followed by experiments (Section VI), results (Section VII), and discussions (Section VIII).

II. NOTATION

Throughout the paper, the underlying image space is the Euclidean space R L equipped with the standard scalar product ⟨•, •⟩, and the norm ∥ • ∥. Moreover, |||L||| denotes the spectral norm of squared matrix L ∈ R L×L . Let S + L be the set of symmetric positive definite matrices in R L×L . For

Q ∈ S + L , ∥ • ∥ Q denotes the Q-weighted norm, i.e., for every x ∈ R L , ∥x∥ Q = ⟨x | Qx⟩. The diagonal matrix with diagonal entries equal to vector z = (z l ) L ℓ=1 ∈ R L is denoted diag((z ℓ ) L ℓ=1
). The class of functions which are proper, convex, lower-semicontinuous, defined on R L and taking values in R ∪ {+∞} is denoted by Γ 0 (R L ). The proximity operator of g ∈ Γ 0 (R L ) at x ∈ R L is uniquely defined as [START_REF] Bauschke | Convex Analysis and Monotone Operator Theory in Hilbert Spaces[END_REF]

prox g (x) = argmin z∈R L g(z) + 1 2 ∥x -z∥ 2 .
The indicator function of a nonempty subset C of R L is the function ι C equal to 0 on C and +∞ outside of C. If C is closed and convex, the projection onto C is denoted by proj C .

III. PROBLEM FORMULATION

We consider a 1D detector array of B ∈ N bins rotating around an object. The detector is too short to measure the projections of the entire object; its size defines a circular ROI we aim to reconstruct. Let S ∈ N be the number of projection angles. The vector of truncated sub-sampled tomographic data is y ∈ R T with T = B S. A reconstruction of the object attenuation map in the ROI can be obtained by considering a model of the form:

Hx G = y + n, (1)
where n ∈ R T accounts for some acquisition noise, x G ∈ R L is the scanned image restricted to a grid G with support larger than the ROI but smaller than that of the entire object, and H ∈ R T ×L is the projector that models projection over this intermediary grid G. Operator H contains a subset of the columns of the projector on the entire space, or equivalently it corresponds to setting a subset of the columns of the full projector, corresponding to the pixels outside of G, to zero.

When the grid G corresponds to the ROI, (1) assumes that the image values are 0 outside of the ROI. This assumption is not necessarily valid for truncated data, so (1) does not hold.

Hereafter, we suppose that the grid is extended beyond the ROI so that no assumption is made about the values outside the ROI. MBIR finds an estimate of x G by computing a minimizer of a penalized cost function that can, for instance, be written as the sum of a data fidelity term f involving H and y, and a regularization term r, as

argmin x∈R L f (x) + r(x). (2)
Our proposed neural network will unfold a baseline algorithm for solving the MBIR problem (2) for a specific choice of f and r. The next section presents the two terms of the cost function and the baseline nested-loop minimization algorithm.

IV. MODEL-BASED ITERATIVE RECONSTRUCTION A. Cost function 1) Cauchy data fidelity: When objects with high gradients (metallic wires or needles) do not belong to the reconstruction grid G, the error between the data and the reprojection of the estimate over G contains outliers at the projections of those high gradients. Angular sub-sampling of these outliers leads to streaks originating from these objects, i.e., from outside the grid. This means the data should not be trusted equally but through a statistical analysis different from measurement noise. To avoid the streaks, we propose to decrease the influence of the largest errors between y and the reprojection of Hx G using an M-estimator. Here, we focus on the Cauchy estimator ϕ, which is a redescending M-estimator [START_REF] Andrews | Robust estimates of location[END_REF] i.e., its derivative decreases to zero on ] -∞, κ] ∩ [κ, +∞[. It is defined as

(∀ζ ∈ R) ϕ(ζ) = βκ 2 2 ln 1 + ζ κ 2 , (3)
where β > 0 is a weighting term and κ > 0 monitors the sensitivity to outliers: the lower κ, the lower the influence of the outliers.

A graphical comparison of the Cauchy function (3) and the quadratic function ϕ(•) = β 2 (•) 2 is displayed in Figure 1 for β = 1 and various values of κ. Our data fidelity term f reads

(∀x ∈ R L ) f (x) = g(Hx -y), (4)
with

(∀z = (z t ) 1⩽t⩽T ∈ R T ) g(z) = T t=1 ϕ(z t ). (5)
2) Semi-local total variation: Sub-sampling streaks are commonly handled by total variation (TV) regularization [START_REF] Sidky | Image reconstruction in circular cone-beam computed tomography by constrained, total-variation minimization[END_REF], [START_REF] Rudin | Nonlinear total variation based noise removal algorithms[END_REF], [START_REF] Langet | Compressed-sensing-based content-driven hierarchical reconstruction: Theory and application to C-arm cone-beam tomography[END_REF]. Semi-local variants (STV) [START_REF] Kindermann | Deblurring and denoising of images by nonlocal functionals[END_REF] extend TV in a neighborhood of pixels indexed in Λ J = {-J, . . . , J} \ {0}:

(∀x ∈ R L ) r STV (x) = J j=1 L ℓ=1 α j,ℓ (x -V j x) 2 ℓ + (x -V -j x) 2 ℓ = J j=1 r j (∇ j x). (6)
Hereinabove ℓ ∈ {1, . . . , L} is the spatial index and V j , V -j ∈ R L×L are shift operators as shown in Figure 2 for j ∈ {1, . . . , J} and J = 6. Moreover, for every j ∈ {1, • • • , J}, we define

∇ j = V ⊤ j V ⊤ -j ⊤ ∈ R 2L×L and, for every z = (z 1 , z 2 ) ∈ R 2L , r j (z) = L ℓ=1 α j,ℓ (z 1 ) 2 ℓ + (z 2 ) 2
ℓ . Parameters (α j,ℓ ) 1⩽j⩽J,1⩽ℓ⩽L are nonnegative weights that can be chosen to vary spatially, so making STV adaptive to the spatial contents [START_REF] Gilboa | Nonlocal operators with applications to image processing[END_REF]. We recover the standard TV regularization for constant values of these parameters and J = 1. We add a nonnegativity constraint on the pixel values and a quadratic term 1 2 ∥x∥ 2 M = x ⊤ Mx to the STV regularization. Here, matrix M = diag((m ℓ ) L ℓ=1 ) ∈ S + L is such that, for every ℓ ∈ {1, . . . , L}, m ℓ = 1 if the ℓ-th entry x ℓ of vector x belongs to the ROI, and m ℓ = ξ > 1 otherwise. Thus M ∈ S + L acts as a mask, limiting high values outside of the ROI.

Altogether, our regularization function in (2) reads

(∀x ∈ R L ) r(x) = J j=1 r j (∇ j x) + 1 2 ∥x∥ 2 M + ι [0,+∞[ L (x).
(7) To our knowledge, STV regularization has not been used in CT for handling sub-sampling streaks. In this work, STV provides extra capacity compared to TV for learning. The Cauchy fidelity term has been used in ultrasound imaging [START_REF] Ouzir | Robust optical flow estimation in cardiac ultrasound images using a sparse representation[END_REF] and in CT imaging [START_REF] Kazantsev | A novel tomographic reconstruction method based on the robust student's t function for suppressing data outliers[END_REF] for mitigating the ring artifacts that appear due to defective detector bins only.

B. Minimization algorithm 1)

Reweighting for non-convex data fidelity: Given our choices for r and f , Problem (2) becomes argmin

x∈R L g(Hx -y) + J j=1 r j (∇ j x) + 1 2 ∥x∥ 2 M + ι [0,+∞[ L (x).
(8) Because of the non-convexity of g, we adopt an iterative reweighting strategy where Problem ( 8) is replaced by a sequence of surrogate convex problems built following a majoration principle. As shown in Appendix A, the associated Majorization-Minimization strategy is given by Algorithm 1,

Input: Number of iterations K ∈ N * , x 0 ∈ R L for k = 0 to K -1 do
Define majorant function Q(x; x) using ( 10)

x k+1 = argmin x∈R L Q(x, x k ) (9)
end for Output: x K approximating the solution to [START_REF] Paleo | Efficient implementation of a local tomography reconstruction algorithm[END_REF] where, for every x ∈ R L , the surrogate function

Q(•, x) is Q(x, x) = ι [0,+∞[ L (x)+h 0 (B 0 x; B 0 x)+h 1 (B 1 x)+ 1 2 ∥x∥ 2 M , (10) where ∀(x, x) ∈ (R L ) 2 , with B 0 = H ∈ R T ×L , B 1 = ∇ ⊤ 1 • • • ∇ ⊤ j ⊤ ∈ R 2JL×L h 0 (•; B 0 x) = g(• -y; B 0 x -y) h 1 (B 1 x) = J j=1 r j (∇ j x). (11)
Let (x k ) k∈N be generated by Algorithm 1. The cost sequence value (defined from ( 8)) monotonically converges.

2) Dual block coordinate forward-backward algorithm: The k-th iteration of Algorithm 1 requires to solve [START_REF] Hamelin | Design of iterative ROI transmission tomography reconstruction procedures and image quality analysis[END_REF], which amounts to minimizing the function Q(•, x), with x equals to the current iterate x k . Since M ∈ S + L , Q(•, x) is strongly convex for every x ∈ R L . The minimization ( 9) is hence well-defined, with a unique solution that can be conveniently obtained using the dual forward-backward algorithm [START_REF] Le | The faster proximal algorithm, the better unfolded deep learning architecture ? The study case of image denoising[END_REF], [START_REF] Combettes | Proximal splitting methods in signal processing[END_REF]. An accelerated and light version of this algorithm is its block coordinate version (DBFB) [START_REF] Abboud | Dual block-coordinate forward-backward algorithm with application to deconvolution and deinterlacing of video sequences[END_REF], which allows for accessing the proximity operators of h 0 and h 1 separately. Algorithm 2 describes N ∈ N * iterations of DBFB. The output x N generated by DBFB with input x k then defines x K+1 in Algorithm 1. For every n ∈ {1, . . . , N }, DBFB updates the main primal variable x n as well as two dual variables z 0 n ∈ R T and s 1 n ∈ R 2JL , associated with the data fidelity (data step (D)) or the regularization (regularization step (R)) terms, respectively. Each dual variable is activated (or not) at iteration n according to a binary variable ε n . When N → ∞, the DBFB sequence (x n ) n∈N converge to the solution to (9) under the following assumptions on the algorithm parameters [START_REF] Abboud | Dual block-coordinate forward-backward algorithm with application to deconvolution and deinterlacing of video sequences[END_REF]: [START_REF] Langet | Compressed-sensing-based content-driven hierarchical reconstruction: Theory and application to C-arm cone-beam tomography[END_REF], operators B 0 and B 1 defined in [START_REF] Savanier | Unmatched preconditioning of the proximal gradient algorithm[END_REF], stepsizes (σ, τ 1 , . . . , τ J ) ∈]0, +∞[ J+1 .

         σ ⩾ |||B 0 M -1 B ⊤ 0 |||, (∀j ∈ {1, . . . , J}) τ j ⩾ |||∇ j M -1 ∇ ⊤ j |||, γ n ∈ [ϵ, 2 -ϵ] with ϵ ∈]0, 1] (∃M ∈ N \ {0, 1})(∀n ∈ N) 0 < n+M -1 n ′ =n ε n ′ < M. Algorithm 2 DBFB algorithm to solve (9) with x k = x ∈ R L Input: Number of iterations N ∈ N * , tangent point x ∈ R L , initial dual variables z 0 0 ∈ R T , (∀j ∈ {1, . . . , J}) z j 0 ∈ R 2L with constant J defined in
s 1 0 = (z 1 0 , . . . , z J 0 ) Σ = diag(τ 1 , . . . , τ J ) w 0 = -M -1 (B ⊤ 0 z 0 0 + B ⊤ 1 s 1 0 ) For n = 0, 1 . . . , N                         (z 1 n , . . . , z J n ) ≡ s 1 n x n = proj [0,+∞[ L (w n ) Select ε n ∈ {0, 1} and γ n ∈]0, +∞[ If ε n = 0 (D)       z0 n = z 0 n + γ n σ -1 B 0 x n z 0 n+1 = z0 n -γ n σ -1 prox γ -1 n σh0(•;B0x) (γ -1 n σz 0 n ) w n+1 = w n -M -1 B ⊤ 0 (z 0 n+1 -z 0 n ) s 1 n+1 = s 1 n If ε n = 1 (R)       s1 n = s 1 n + γ n Σ -1 B 1 x n s 1 n+1 = s1 n -γ n Σ -1 prox γ -1 n Σh1 (γ -1 n Σs 1 n ) w n+1 = w n -M -1 B ⊤ 1 (s 1 n+1 -s 1 n ) z 0 n+1 = z 0 n Output: x N approximating the minimizer of Q(•, x).
The first three assumptions are stepsize range conditions. The last one means that each step (D) and (R) is performed at least once every M iterations. The practical implementation of the proximal operators involved in DBFB for our choices for h 0 and h 1 is discussed in Appendix B.

The overall iterative strategy for approximating the solution to (8) consists of applying Algorithm 1, where, for every k ∈ N * , N k ∈ N * iterations of Algorithm 2 are used as an inner solver with x = x k to compute x k+1 in [START_REF] Hamelin | Design of iterative ROI transmission tomography reconstruction procedures and image quality analysis[END_REF]. We call the resulting iterations reweighted DBFB (RDBFB) algorithm (see Figure 3). Such a combination of an iterative reweighted algorithm with dual ascent steps is original to our knowledge.

In the context of CT, iterative reweighted algorithms usually involve surrogates to the regularization term [START_REF] Yu | Edge-preserving reconstruction from sparse projections of limited-angle computed tomography using ℓ 0 -regularized gradient prior[END_REF], [START_REF] Wang | Reweighted Anisotropic Total Variation Minimization for Limited-Angle CT Reconstruction[END_REF] rather than to the data fidelity, as done here.

V. UNFOLDED RECONSTRUCTION Hereinafter, we present a deep neural network, designated as U-RDBFB (Unfolded Reweighted DBFB), by unfolding all the steps of RDBFB. Specifically, the network mimics the application of K iterations of Algorithm 1, as K main layers, each of them grouping N k ∈ N iterations of Algorithm 2. This yields an architecture with K-1 k=0 N k layers in total.

A. From RDBFB iterations to U-RDBFB layers

The deep unfolding paradigm recasts every step of Algorithm 2 as one neural network layer: step (D) becomes L D (ε n = 0) and step (R) becomes L R (ε n = 1). It requires truncating the number of layers drastically. To optimize the depth of our network, we propose two modifications of the steps (D) and (R) to construct the corresponding layers. These modifications stem from recent works [START_REF] Chouzenoux | Convergence of proximal gradient algorithm in the presence of adjoint mismatch[END_REF]- [START_REF] Chouzenoux | Convergence Results for Primal-Dual Algorithms in the Presence of Adjoint Mismatch[END_REF] that have shown the potential to replace adjoint operators with surrogates to accelerate the convergence of MBIR methods [START_REF] Zeng | Counter examples for unmatched projector/backprojector in an iterative algorithm[END_REF]- [START_REF] Sidky | Iterative image reconstruction for CT with unmatched projection matrices using the generalized minimal residual algorithm[END_REF]. In CT, a frequently encountered operator is the ramp filter F, which satisfies FHH ⊤ ≈ Id . Thus, to improve conditioning and allow larger values for the step sizes associated with layer L D , hence a lower number of such a layer, we replace B 0 in Algorithm 2 with FH. By setting ν n,0 = γ n σ -1 and by using the relation

prox νn,0h0(•;FHx) = prox νn,0 g(•;Hx-y) (• -Fy) + Fy, we define layer L D as Data layer (L D ):         x n = proj [0,+∞[ L (w n ) u n = z 0 n + ν n,0 F(Hx n -y) z 0 n+1 = u n -ν -1 n,0 prox νn,0 g(•;FHx-y) (ν n,0 u n ) w n+1 = w n -M -1 H ⊤ (z 0 n+1 -z 0 n ) z j n+1 = z j n (∀j ∈ {1, . . . , J}). (12)
Similarly, for step (R), we unfold by replacing the adjoint of the regularization operator

B 1 with B1 = ∇⊤ 1 • • • ∇⊤ j ⊤ .
Setting, for every j ∈ {1, . . . , K}, ν n,j = γ n τ -1 j yield the following regularization layer L R :

Regularization layer (L R ) :               x n = proj [0,+∞[ L (w n ) For j ∈ {1, . . . , J}      (∀ℓ ∈ {1, . . . , L}) (z j n+1 ) ℓ = z j n + ν n,j ∇ j x n ℓ max 1, ∥ z j n + ν n,j ∇ j x n ℓ ∥ 2 /α j,ℓ w n+1 = w n -M -1 J j=1 ∇j (z j n+1 -z j n ) z 0 n+1 = z 0 n . (13)
Note that L R does not involve x.

B. Total architecture

The total architecture of U-RDBFB, denoted A, can be summarized as

A = L K-1 • • • • • L 0 . ( 14
)
For 0 ⩽ k ⩽ K -1, L k corresponds to a sequence of N k layers L D or L R and implements the following update:

(z 0,k+1 , x k+1 ) = L k (z 0,k , x k ; Θ k ), (15)
where

• x k is the current reweighted estimate x in L R -L D (x 0 = H ⊤ Fy). • x k+1 is the next reweighted estimate; it is equal to x N k -1
given by the For layer L R , the regularization parameters (α k,n,j,ℓ ) 1⩽j⩽J,1⩽ℓ⩽L , κ k,n , ξ k,n , and step sizes (ν k,n,j ) J j=1 are learned as well as the surrogates ( ∇k,n j ) J j=1 to the adjoints of operators (∇ ⊤ j ) J j=1 . To infer all these parameters, we introduce learning modules in each layer (L

N k -th layer L R -L D . • z 0,k ∈ R T × (R 2L ) J is the initial value of the variables (z j 0 ) 0⩽j⩽J for layers L R -L D (for k = 0, (z j 0 ) J j=1 are initialized to zero while z 0 0 is set to -Fy). • z 0,k+1 ∈ R T × (R 2L ) J is equal to (z j N k-1 ) 0⩽j⩽J given by the N k -th layer L R -L D • Θ k is
(θ) k,n ) N k -1 n=0 for θ ∈ Θ k (implementation discussed in Appendix C).
Schematic views of layers L D and L R can be found in Figure 9a and Figure 9b, and a composition A of such layers is displayed in Figure 9c.

Here we propose using K = 7 in [START_REF] Nilchian | Iterative FBP for improved reconstruction of X-ray differential phasecontrast tomograms[END_REF] with N k = 4 for each k ∈ {0, • • • , K -1}, resulting in a total of 28 layers:

L 0 = • • • = L 6 = (L R • L D ) 2 . ( 16
)

C. Incremental training strategy

For each learning module, we circumvent the issue of optimizing the initial values of their parameters using an incremental training strategy, as sometimes advocated for when initializing the weights of recurrent neural networks [START_REF] Safran | On the quality of the initial basin in overspecified neural networks[END_REF], [START_REF] Martens | Learning recurrent neural networks with hessian-free optimization[END_REF]. The learning in each layer (k, n) (n ⩾ 1) starts by considering all the previous layers from (0, 0) to (k, n -1) with their past trained parameters as an initialization. This means that an increasing number of layers are trained simultaneously. In the last step, all layers are trained end-to-end.

VI. EXPERIMENTS

We now describe our experimental setup. First, we illustrate the benefits of using a Cauchy-based data fidelity function. We compare the results of our RDBFB iterative approach to a simpler DBFB scheme minimizing the same cost function but with g in ( 8) replaced with the ℓ 2 norm. Second, we comment on the improved performance brought by our unfolding strategy (learning of adjoints and parameters as well as the use of the ramp filter), and so U-RDBFB (Section V) is compared to the original RDBFB algorithm.

Third, we compare U-RDBFB to several state-of-the-art reconstruction methods and comment on the transfer of performance over different synthetic datasets.

A. Datasets

We used three datasets for evaluation. 1) Abdomen dataset: Our first dataset consists of 2D images obtained from 60 CT volumes of size 512 × 512 × 512 from the lower lungs to the lower abdomen of 60 patients, which were extracted from the public dataset CT Lymph Nodes from https://www.cancerimagingarchive.net/. These volumes correspond to fully sampled CT reconstructions. They were made isotropic by interpolating the axial slices. A total of 50 out of 512 slices were kept per volume. We randomly added intense metallic wires between 3000 and 5000 Hounsfield units (HU) of varying sizes on the axial slices. We shifted the HU values of the images by 1000 so that air is 0 HU and water is 1000 HU, using a tissue → (a tissue -µ water ) × (1000/µ water ), where µ water is the value of the attenuation coefficient of water equal to 0.017 mm -1 and a tissue ∈ R Q is the initial vector of attenuation values. Finally the 512 × 512 slices x P were normalized between [0, 1] (a value of 1 corresponding to an object of HU intensity equal to 5000).

To eliminate bias with respect to model discretization, projections were simulated for each slice of each volume in a 2D parallel geometry with a short detector of 600 bins (bin size equal to half a pixel size, i.e., 0.5 mm) and an angular density of 110 projections over 180 °through to the forward projector H P ∈ R 512 2 ×(110×600) . The projections were then rebinned by a factor 2 (operator R ∈ R 300×600 ). Noisy projections y = (y t ) T t=1 are computed as (∀t ∈ {1, . . . , T }) y t = µ log I 0 P(I 0 exp (-µ(RH P x P ) t ))

,

where we set µ = µ water /1000, I 0 = 10 4 , and, for some δ > 0, P(δ) denotes a realization of a Poisson law with mean δ. In this context, the ROI was a centered disk of diameter 300.

The resulting pairs of axial slice/projections (x P , y) were split into a training of 2500 pairs from a pool of 50 patients and a testing set of 500 pairs from 10 other patients.

2) Head dataset: We used a second dataset containing 2D images extracted from 10 CT high-dose brain reconstructions. These volumes are from the public repository 2016 Low Dose CT Grand Challenge from https://www.cancerimagingarchive.net/. After extracting 50 slices of size 512 × 512 per volume (pixel size of 0.5 mm), we performed the same processing as for the Abdomen dataset (addition of intense wires, normalization, projection, rebinning) for generating a testing set of 500 pairs of axial slices/projections (x P , y).

3) Geometrical dataset: Our third dataset was created using the toolbox TomoPhantom [START_REF] Kazantsev | Tomo-Phantom, a software package to generate 2D-4D analytical phantoms for CT image reconstruction algorithm benchmarks[END_REF]. 500 geometrical 2D piecewise-constant phantoms were randomly generated on a 512 × 512 grid and normalized between 0 and 1. Again, we performed the same processing as for the Abdomen dataset for generating a testing set of 500 pairs of axial slices/projections (x P , y).

B. Training details for U-RDBFB

Let i ∈ {1, . . . , I} be the index covering all I = 2500 instances of the training set. The reconstruction grid (G) is a disk of diameter 400. Let x * G,i ∈ R L be the output of U-RDBFB for a given projection input y i . Our network is thus designed to minimize

I i=1 ℓ(C G x * G,i , C P x P,i )
, where C G is a cropping operator which extracts the ROI from the grid G, C P is a cropping operator which extracts the ROI from the entire 512 × 512 grid, and ℓ is the loss retained for training the network. For all instances of the training set,

ℓ(C G •, C P x P,i ) = 1 I ∥C G • -C P x P,i ∥ 2 , (17)
corresponding to the MSE loss. We implemented U-RDBFB following [START_REF] Han | Deep residual learning for compressed sensing CT reconstruction via persistent homology analysis[END_REF] in Pytorch, using a Tesla V100 32 Gb GPU. We used six epochs for training each layer L R and ten epochs for training each layer L D ; the only exception was the last layer, for which we used 20 epochs. The learning rate is decreased with a step decay by a factor of 0.99 from 10 -2 every four epochs. The batch size for each epoch varied from 20 to 8 as the number of trained layers increased. We employed the toolbox TorchRadon [START_REF] Ronchetti | TorchRadon: Fast differentiable routines for computed tomography[END_REF] to include Pytorch-compatible parallel-beam tomographic operators in all architectures. Standard auto-differentiation tools can compute all necessary derivatives for backpropagation. The training procedure takes about one day and a half.

C. Competing methods

The quantitative metric used to assess the reconstruction quality of C G x * G,i is the PSNR. We also evaluate the reconstruction performance using the structural similarity index (SSIM), the PieApp value [START_REF] Prashnani | Pieapp: Perceptual imageerror assessment through pairwise preference[END_REF], and the Mean Absolute Error (MAE) of the difference between C G x * G,i and C P x P,i . We compare U-RDBFB with FBP, an iterative method, and four deep-learning methods that we describe hereinafter.

1) FBP: This analytical method consists of computing H ⊤ ROI Fy, where H ROI ∈ R 300 2 ×(110×300) . As is commonly the case when applying FBP on truncated data, we extrapolated the projections prior to ramp filtering (using anti-symmetric padding).

2) RDBFB algorithm: For completeness, we perform comparisons with the iterative method proposed in Section IV. For each reweighted iteration k, we used N k = 10 DBFB iterations alternating between data and regularization steps (1:1 correspondence). For an easier manual tuning of the hyperparameters, instead of using J = 6 as in U-RDBFB, we set J = 1 so that STV reduces to TV and α 1,ℓ ≡ α 1 , for all ℓ ∈ {1, • • • , L}. The remaining cost function parameters (ξ, κ, β) are selected by optimizing PSNR on the training set via a grid search.

3) Post-processing U-net: The third competing method is the CNN proposed in [START_REF] Han | Deep residual learning for compressed sensing CT reconstruction via persistent homology analysis[END_REF], [START_REF] Jin | Deep convolutional neural network for inverse problems in imaging[END_REF], which is a post-processing of FBP. It relies on a trained residual U-net, with a depth of 4 levels, filters of size 32, and batch normalization to improve the stability of training.

4) Preconditioned Neumann Network (PNN):

Our fourth competing method is a preconditioned Neumann network (PNN) initially introduced in [START_REF] Gilton | Neumann networks for linear inverse problems in imaging[END_REF] for MRI reconstruction. It builds on a method for solving Problem (2) with f (x) = 1 2 ∥Hx -y∥ 2 . For a differentiable function r, the resulting minimizer reads

(H ⊤ H + ∇r)x = H ⊤ y, (18)
which can be rewritten as

(H ⊤ H + λ Id )x + (∇r -λ Id )x = H ⊤ y. ( 19
)
Setting

T λ = (H ⊤ H + λ Id ) -1 yields (Id -λT λ + T λ ∇r)x = T λ H ⊤ y. ( 20
)
Using the Neumann identity

B -1 = ∞ n=0 (Id -B) n
, the authors derive the architecture of PNN with N ∈ N * layers (see Figure 4)

(λT λ -T λ ∇r) N • T λ (H ⊤ y). ( 21
)
All instances of T λ are applied approximately using an unrolling of 10 iterations of the conjugate gradient algorithm. The operator T λ ∇r is replaced by a U-net, denoted by Ψ, which has the same architecture as the aforementioned U-net without the residual connection. The weights of the U-net are shared for all layers. Following [START_REF] Gilton | Neumann networks for linear inverse problems in imaging[END_REF], no batch normalization is used. The inner U-net has a depth of 4, the learning rate is set to 10 -4 , and the initial value for λ is 0.01. We choose N = 3. One feature of PNN compared to other deep unfolding networks is that it contains skip connections. 5) ISTA-net: Our fifth competing method is ISTA-net, derived from the work of [START_REF] Zhang | ISTA-Net: Interpretable optimization-inspired deep network for image compressive sensing[END_REF]. ISTA-net is designed to solve Problem (2) for f (x) = 1 2 ∥Hx -y∥ 2 , and r(x) = λ∥Wx∥ 1 (λ > 0), where operator W is not known a priori but learned. W is an orthogonal linear operator in the initial ISTA algorithm, whose iteration reads

x n+1 = W ⊤ soft W(x n -τ H ⊤ (Hx n -y), λτ , ( 22
)
where soft is the soft-thresholding operation and τ > 0 is the gradient step size. In ISTA-net, the authors replace W and W ⊤ by two decoupled nonlinear operators namely 5a). The property of orthogonality of W is not imposed but favored during training by adding a term, weighted by χ ∈ ]0, +∞[, penalizing the difference between

A n • ReLU •B n and C n • ReLU •D n (see Figure
(C n • ReLU •D n ) • (A n • ReLU •B n )x

6) PD-net:

The last competing method is the learned Primal-Dual (PD-net) introduced in [START_REF] Adler | Learned primal-dual reconstruction[END_REF] by unrolling the Primal-Dual Hybrid Gradient (PDHG) optimization algorithm [START_REF] Cheng | Model learning: Primal dual networks for fast MR imaging[END_REF]. The authors consider Problem (2) with a more generic data fidelity term f (x) = G(Hx; y). They replace both the proximity operators of G and r in PDHG by residual CNN so that one layer n of their network reads

z n+1 = CNN(z n + σHx n ; y) ( 23
)
x n+1 = CNN(x n -τ H ⊤ z n+1 ) (24) xn+1 = x n+1 + γ(x n+1 -x n ). ( 25
)
The CNNs act both in the image and projection domains. Furthermore, buffers of previous iterates of size N p ∈ N in the primal domain (image) and of size N d ∈ N in the dual domain (projection) are kept to enable the network to learn an acceleration. We used 9 layers, N d = N p = 3, and 32 filters in the convolutional layers. This network is illustrated in Figure 5b. The competing networks were also trained with the MSE loss (using [START_REF] Jin | Deep convolutional neural network for inverse problems in imaging[END_REF] for unfolding networks and a regularization term for ISTA-net weighted by χ) in a standard end-to-end manner. The number of epochs was chosen such that the training and testing losses have stabilized. Note that codes are publicly available for these networks. We re-implemented them in Pytorch and kept the setting of the parameters advocated by the authors, except for PNN, for which we reduced the number of layers to 3 to obtain a stable behavior for training. The total number of parameters of each network is reported in Table I. Ψ is a trained network and the scale parameter λ is also trained.

VII. RESULTS

A. Assessing the benefits of the Cauchy fidelity term

Figure 6c shows the full reconstructed images on grid G of size 400 × 400 obtained using the DBFB algorithm with a quadratic fidelity term and the reweighed DBFB algorithm with a Cauchy fidelity on a test instance of the Abdomen dataset (shown in Figure 6a). Since the two data fidelity terms can be put into our optimization framework Alg. 1-Alg. 2, the comparison is straightforward. Figure 6d shows The full image contains two intense objects out of the ROI and at the border of the reconstruction grid. In the solution obtained using the quadratic data fidelity term, the reduction of sub-sampling streaks is selective; only the streaks originating from objects within G have been eliminated in the ROI. When trading the quadratic term with a Cauchy term, as we proposed, the intensity of these streaks is reduced. This artifact reduction translates into an improvement of the PSNR as shown in Figure 6b. Figure 6e shows the ROIs obtained using the same reconstruction methods and grid size when increasing the number of projections from 110 to 600. The images now look identical and close to the ground truth. This observation highlights that, for relatively 'clean' data (no modeling of beam hardening and scattering), the benefits of using a Cauchy fidelity over a quadratic fidelity emerge when data is subsampled. We have shown that by using a regularized cost function with a Cauchy fidelity term and, thus, a more complex optimization framework, we can successfully reconstruct truncated data on a short grid and that the reconstruction is at least as good as the one obtained with quadratic fidelity and even better when the data are sub-sampled.

B. Comparing iterative RDBFB algorithm with U-RDBFB

Figure 7a shows the evolution of the PSNR along the iterations when inserting the ramp filter in the reweighted DBFB algorithm, re-tuning the regularization strength, and still using the same data. The PSNR stagnates around 300 iterations with the ramp filter while it stagnates around 1250 iterations without it (see Figure 6b). Thus applying the ramp filter on the reprojection error before backprojection can empirically accelerate convergence without degrading the solution (reconstructed ROI displayed in Figure 7c) in an early stopping scenario. This motivates our translation of a data iteration of DBFB to a data layer of U-RDBFB which embeds the ramp filter. It also provides empirical evidence that performance can be optimized by introducing mismatched adjoints without learning. U-RDBFB also includes learned parameters, especially adjoints to the STV operators. It performs a total of 28 RDBFB iterations. Figure 7c compares the reconstruction ROI obtained with U-RDBFB, 500 and 28 iterations of the reweighed DBFB algorithm with the ramp filter. We see that after 28 iterations of the reweighed DBFB algorithm, there is a local offset near the intense object, and some streaks remain. On the contrary, the image obtained with U-RDBFB does not contain these artifacts. It is similar to the ground truth and slightly smoother than the reconstruction after 500 iterations of the reweighed DBFB algorithm but with approximately 18 times fewer iterations.

C. Comparing U-RDBFB with deep learning methods on the Abdomen dataset

Table II reports the performance of U-RDBFB compared to U-net and other deep unfolding networks on the testing set of the Abdomen dataset and Figures 8a-8b display the training and testing losses as a function of the number of epochs for all these networks. U-RDBFB performs, on average, better than the other unfolding networks (PNN, PD-net, and ISTAnet) and U-net for all considered metrics. We note that the peaks in the training and testing losses associated with U-RDBFB correspond to the addition of a new data layer during incremental training. Figure 10c illustrates the reconstructed ROIs for four examples from the test set of the Abdomen dataset. The FBP reconstruction is also displayed as it is also the input of U-net. The figure confirms that U-RDBFB reduces streaks more effectively than the other unfolding networks. At first sight, the images produced by U-net have fewer artifacts than most deep unfolding networks. However, in the secondrow and fourth-row images, U-net introduces an artificial dark background. This observation highlights that U-net can hallucinate structures under the sub-sampling streaks of the FBP input. Unfolding networks avoid these hallucinations; by simply alternating between U-net and several consistency layers, PNN already minimizes this effect. Figure 10b shows the complete reconstruction on grid G for all unfolding networks. In all cases, since the training loss acts on the ROI only, the exterior is always poorly reconstructed (with ISTA-net, it is very sparse).

D. Changing the testing set

We now evaluate the generalization ability of the methods trained on the Abdomen dataset, using examples from the Head and Geometrical datasets. In Table III and Table IV, we report the performance of the trained networks when tested on the Head and the Geometrical datasets. U-RDBFB outperforms the competing unfolding networks for both datasets. ISTA-net is second-best on the Head dataset, and the iterative algorithm ranks second-best on the Geometrical dataset. The four reconstructed images displayed in Figures 11a-11b confirm this trend. For the Head dataset, U-net performs poorly relative to U-RDBFB on all metrics except for SSIM, where the two methods are rather close. One explanation is that it introduces an offset in some images while limiting the streaks (cf second head image). However, when applied to the geometrical images, it yields unwanted background patterns that strongly degrade our metrics. Offsets and background artifacts are also visible with most unfolding networks, especially PNN and PD-net, except for U-RDBFB. We still note that the head images and the first geometrical image obtained with U-RDBFB have a slightly patchy look, often characteristic of TV regularization. This suggests that U-RDBFB retains the characteristics of the original optimization problem, avoiding the generation of unexpected content as is possible with U-net.

VIII. DISCUSSION

The results show that U-RDBFB outperforms its iterative counterpart as measured in PSNR, SSIM, and PieApp. It not only presents results similar to its iterative counterpart RDBFB in terms of streaks reduction, but it also recovers image details in a much lower number of iterations. We note that U-RDBFB leads to smoother images: the noise is reduced, but the resolution of the image is somehow decreased. This could be explained by the fact that U-RDBFB is designed to minimize the ℓ 2 norm of the error with respect to the noiseless ground truth, using very few regularization layers; it tends to selectively smooth some parts of the image to remove remaining artifacts. All metrics agree that U-RDBFB improves upon learned postprocessing U-net and other unfolding networks for our Abdomen dataset. U-net was often associated with a high PSNR compared to the other reconstructions, but this was not always reflected by the PieApp metric. This may be explained by the learned post-processing being limited by the information content of the FBP input, while the unfolding networks act directly with the information content of the data, which is greater than that of the FBP. The computation time for U-RDBFB was about 200 ms in GPU for a 400 × 400 reconstruction grid. This is much faster than the iterative reconstruction, which, in our case, requires around 180 s after the regularization parameters have been selected, but slower than other deep learning methods (38 ms for U-net, 68 ms for PNN, 74 ms for ISTA-net, 85 ms for PD-net). U-RDBFB contains fewer learnable parameters than all the other networks. Thanks to our incremental strategy, training U-RDBFB was also found to be easier than other unfolding networks, such as PNN, whose stability highly depends on the initialization for parameter λ and the learning rate. Optimizing the architecture and, more precisely, the number of parameters of a neural network is key to transferring its performance to out-of-distribution examples, as shown on the Head and Geometrical datasets. Generally, deep unfolding networks are introduced mainly to ensure data consistency through H and embed a fast optimization scheme for fast inference. Our results suggest that including additional a priori knowledge can further boost the performance of deep-learning-based techniques. Note that the structure of U-RDBFB (K, N k and distribution of L D and L R ) was not extensively fine-tuned. Our results illustrate that the most straightforward choices work well in our context of ROI imaging from angularly subsampled data. Our results also hinted that even without learning, trading a quadratic fidelity for a Cauchy fidelity, discarding the data incompatible with the a priori support of the object, and including some preconditioning through the ramp filter are still of interest to improve reconstruction and reduce the number of iterations.

IX. CONCLUSION

In this paper, we introduced an iterative reweighted algorithm where each inner optimization problem is solved using dual block coordinate forward-backward iterations, and we proposed an unfolded version of it, yielding a neural network for ROI reconstruction from a few measurements. These methods include a convex surrogate to a Cauchy data fidelity and a TV-based regularization to limit sub-sampling streaks originating from inside and outside the reconstruction grid. Our experiments demonstrated the benefits of the unfolded strategy over the original iterative algorithm. By balancing the capacity of the network and the use of prior knowledge, our architecture displayed high generalization ability compared to various neural networks, including U-net and other deep unfolding networks. Future work will investigate loss functions that better preserve image resolution and further optimization of the structure of our network for application to other geometries (e.g., cone-beam) using real CT data.

APPENDIX A REWEIGHTING STRATEGY

Let ϕ be given by (3). It was shown in [START_REF] Chouzenoux | A majorizeminimize subspace approach for ℓ 2 -ℓ 0 image regularization[END_REF] that, for every ζ ∈ R, the following convex quadratic function φ(•, ζ), defined for every ζ ∈ R as Finally, for every x ∈ R L , we set (∀x ∈ R L ) f (x, x) = g(Hx -y; Hxy), [START_REF] Abboud | Dual Block Coordinate Forward-Backward Algorithm with Application to Deconvolution and Deinterlacing of Video Sequences[END_REF] that satisfies f (x, x) ⩾ g(Hx -y) = f (x). Given this majoration, the iterative reweighting strategy approximates the solution to (8) by the estimate produced by Algorithm 1, where

φ(ζ, ζ) = ϕ(ζ) + β (ζ -ζ)ζ 1 + (ζ/κ) 2 + β 2 (ζ -ζ) 2 (1 + (ζ/κ) 2 ) , (26)
(∀(x, x) ∈ (R L ) 2 ) Q(x, x) = f (x, x) + r(x). (29)
(a) Schematic view of a N k -th layer L D [START_REF] Tirer | On the convergence rate of projected gradient descent for a back-projection based objective[END_REF]. The layer relies on x = x k , the k-th reweighted iterate to generate the next reweighted iterate x k+1 . The layer takes as inputs wn, z 0 n from the previous layers. The projections y are also used as input.

(b) Schematic view of a N k -th layer L R [START_REF] Tirer | Back-projection based fidelity term for ill-posed linear inverse problems[END_REF]. The layer takes as inputs wn, (z j n ) J 1 from the previous layer and generates the next reweighted iterate x k+1 . The update of parameters for j ∈ {2, . . . , 5} is hidden in the orange block for the sake of readability. Only parameters α j depend on the input.

Fig. 1 :

 1 Fig. 1: Comparison between the Cauchy and the quadratic functions.

Fig. 2 :

 2 Fig. 2: Shift operators (Vj)j∈Λ 6 applied to a given pixel position ℓ

Fig. 3 :Algorithm 1

 31 Fig. 3: Baseline RDBFB algorithm for solving (8): K steps of Alg. 1, where each step k consists of applying N k iterations of Alg. 2.

  the vector of trainable parameters. The parameters are layer-dependent, so we index them by k and n. For layer L D , the parameters are those of the Cauchy function (β k,n , κ k,n ), the one of the quadratic regularization ξ k,n , and a single step size ν k,n,0 .

  n and x n in the loss function. Each A n , B n , C n and D n is a 2D convolutional operator. B n and C n are associated with a kernel of size 3 × 3 and 32 input and output channels; A n has 1 input channel and 32 output channels and vice-versa for D n . As suggested by the authors, we learn these convolutional operators as well as λ and τ , which are allowed to vary at each iteration. Experiments are carried out with 10 layers, χ = 0.1, x 0 is the FBP reconstruction, λ and τ are initialized to 100 and 0.1 respectively.

Fig. 4 :

 4 Fig.4: Architecture of PNN[START_REF] Gilton | Neumann networks for linear inverse problems in imaging[END_REF]: The network maps a linear function of the measurements T λ H ⊤ y to a reconstruction xn by successive applications of an operator of the form λT λ -Ψ, while summing the intermediate outputs of each block. All instances of T λ are replaced by an unrolling of 10 iterations of the conjugate gradient algorithm. Ψ is a trained network and the scale parameter λ is also trained.

  (a) Architecture of ISTA-net[START_REF] Zhang | ISTA-Net: Interpretable optimization-inspired deep network for image compressive sensing[END_REF]: Each layer is composed of a gradient step followed by the application of a nonlinear operator, which is the combination of two learnable linear convolutional operators (An, Bn) separated by a ReLU, a soft-thresholding operation and then two other learnable linear convolutional operators (Cn, Dn) separated by a ReLU. The property (Cn • ReLU •Dn) • (An • ReLU •Bn) = Id is favored during training. (b) Architecture of PD-net [53]: The red and blue boxes represent the primal and dual networks, respectively. Buffers of 3 primal (x 1 n , x 2 n , x 3 n ) and dual (z 1 n , z 2 n , z 3 n ) estimates are used at each iteration. The initial primal estimates are set to the FBP reconstruction given by H ⊤ Fy, and the initial dual estimates are set to zero.

Fig. 5 :

 5 Fig. 5: Two competing unfolded proximal algorithms

Fig. 6 :

 6 Fig. 6: (a) Ground truth xP. (b) Evolution of the PSNR along iterations using a Cauchy or quadratic data term for 110 projections. (c) Reconstructed extended ROIs using 110 projections, a Cauchy data fidelity, and a quadratic data fidelity. (d) Reconstructed ROIs using 110 projections. (e) Reconstructed ROIs using 600 projections. From left to right: Ground truth, reweighted DBFB with Cauchy fidelity, DBFB with quadratic fidelity, FBP.

Fig. 7 :

 7 Fig. 7: (a) Evolution of the PSNR along iterations using a Cauchy fidelity term with and without the ramp filter for the example of Fig 6a. (b) Evolution of the PSNR along layers in U-RDBFB for the example of Fig 6a. (c) Reconstructed ROIs. From left to right: Ground truth, Cauchy with ramp filter (500 iterations of the modified reweighted DBFB), Cauchy with ramp filter (28 iterations of the modified reweighted DBFB), U-RDBFB.

Fig. 8 :

 8 Fig. 8: (a) MSE on the training set as a function of the epoch number. (b) MSE on the testing set as a function of the epoch number.

  is a tangent majorant approximation to ϕ at ζ, that is(∀ζ ∈ R) φ(ζ, ζ) ⩾ ϕ(ζ) and φ(ζ, ζ) = ϕ(ζ). (27)This allows us to deduce a tangent majorant function g of function g at any point z ∈ R T :(∀z ∈ R T ), g(z, z) = T t=1 ϕ(z t , z t ) = g(z) + β diag z t 1 + (z t /κ) 2 T t=1 (z -z) + β 2 (z -z) ⊤ diag 1 1 + (z t /κ) 2 T t=1(z -z) ⩾ g(z).

  (c) U-RDBFB in the case where A = L 1 • L 0 2 where L 0 = (L R ) 2 • L D and L 1 2 = L D (i.e., K = 2, N 0 = 3, N 1 = 1). Red blocks represent the hidden structures to infer all the parameters θ ∈ Θ. When k = 1, the dual variables of DBFB are initialized with the values of the dual variables at the end of the previous N 0 iterations of DBFB (k = 0).

Fig. 9 :Fig. 10 :

 910 Fig. 9: Architecture of U-RDBFB

Fig. 11 :

 11 Fig. 11: Evaluation on the Head and Geometrical datasets.

TABLE I :

 I Number of learnable parameters (Θ)

TABLE II :

 II Quantitative assessment of the reconstructed ROIs. Mean values computed over the test set of the Abdomen dataset.

	Metrics	U-RDBFB	U-net	PPN	PD-net	ISTA-net	RDBFB
	PSNR	31.7	16.8	18.5	17.6	23.1	21.1
	SSIM	0.979	0.902	0.838	0.765	0.957	0.908
	MAE (×10 -3 )	2.17	7.86	4.87	6.47	3.71	6.25
	PieApp	0.276	0.885	0.962	0.984	0.650	0.476

TABLE III :

 III Quantitative assessment of the reconstructed ROIs. Mean values computed over the test set of the Head dataset.

	Metrics	U-RDBFB	U-net	PPN	PD-net	ISTA-net	RDBFB
	PSNR	27.3	25.1	27.0	25.2	25.8	26.1
	SSIM	0.856	0.653	0.733	0.736	0.816	0.848
	MAE (×10 -3 )	16.6	44.1	52.5	25.4	18.4	16.5
	PieApp	0.267	1.170	1.265	1.324	0.897	0.158

TABLE IV :

 IV Quantitative assessment of the reconstructed ROIs. Mean values computed over the testing set of the Geometrical dataset.

This work was supported by the European Research Council Starting Grant MAJORIS ERC-2019-STG-850925, the ANRT CIFRE Convention 2018/1587, and the ANR Research and Teaching Chair in Artificial Intelligence BRIDGEABLE.

APPENDIX B IMPLEMENTATION OF DBFB

Step (D) involves the calculation of the proximity operator prox γ -1 n σh0(•,B0x) , which has a closed-form [29, Example 24.2], for (∀(z, z) ∈ (R T ) 2 ),

Step (R) requires calculating the proximity operator of h 1 scaled by parameter γ ∈ ]0, +∞[. It also has a closed form: for

, where, for every z

where, for every ℓ ∈ {1, . . . , L},

APPENDIX C PARAMETER LEARNING

We discuss our choices for (L

k,n = softplus(a k,n ) where a k,n is a learnable real-valued parameter.

• Parameters of the Cauchy function for

) where c k,n is inferred from a fully connected layer, whose weights are shared across the U-RDBFB network, applied on a histogram of the absolute value of the filtered reprojection error, i.e., H ⊤ (Fx n -y). We implemented the learnable histogram layer proposed in [START_REF] Wang | Learnable histogram: Statistical context features for deep neural networks[END_REF], which is piecewise differentiable. More precisely, we built a cumulated histogram using 100 bins from 0 to the maximum value of the filtered reprojection error.

• Diagonal elements of M involved in [START_REF] Yu | Compressed sensing based interior tomography[END_REF] corresponding to the locations of pixels outside of the ROI for both L D and L R :

) where e k,n is learned.

• step size for L R : For every j ∈ {1, . . . , J}, ν k,n,j = L

). • Parameters of the STV regularization for L R : For every j ∈ {1, . . . , J},

where A k,n is a grouped convolution of 7 groups with size 3 × 3 kernels and J = 7 channels and B k,n is a grouped convolution of 14 groups with size 5 × 5 kernels and 14 channels.

In L D , initial values of a k,n , d k,n are set to 1. In L R , initial values for b k,n,j , e k,n are 1. Normalization scalars W κ , W ν , W β and W α are set to 10 -5 , 10, 10 and 0.05 respectively. Finally, in each layer L R , operations ( Lk,n j ) J j=1 are learned; they are convolutions with a kernel of the same size as for (L k,n j ) J j=1 .